Generalized non-hourglass formulation for ULSPH solid dynamics

src/shared/materials/general_continuum.cpp

@@ -64,9 +64,7 @@ Mat3d PlasticContinuum::ReturnMapping(Mat3d &stress_tensor)
 {
     Real stress_tensor_I1 = stress_tensor.trace();
     if (-alpha_phi_ * stress_tensor_I1 + k_c_ < 0)
-    {
         stress_tensor -= (1.0 / stress_dimension_) * (stress_tensor_I1 - k_c_ / alpha_phi_) * Mat3d::Identity();
-    }
     stress_tensor_I1 = stress_tensor.trace();
     Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / stress_dimension_) * stress_tensor.trace() * Mat3d::Identity();
     Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum();
@@ -77,4 +75,49 @@ Mat3d PlasticContinuum::ReturnMapping(Mat3d &stress_tensor)
     }
     return stress_tensor;
 }
+//=================================================================================================//
+Matd J2Plasticity::ConstitutiveRelationShearStress(Matd &velocity_gradient, Matd &shear_stress, Real &hardening_factor)
+{
+    Matd strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose());
+    Matd deviatoric_strain_rate = strain_rate - (1.0 / (Real)Dimensions) * strain_rate.trace() * Matd::Identity();
+    Matd shear_stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate;
+    Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum();
+    Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_);
+    Real lambda_dot_ = 0;
+    Matd g = Matd::Zero();
+    if (f > TinyReal)
+    {
+        Real deviatoric_stress_times_strain_rate = (shear_stress.cwiseProduct(strain_rate)).sum();
+        lambda_dot_ = deviatoric_stress_times_strain_rate / (sqrt(2.0 * stress_tensor_J2) * (1.0 + hardening_modulus_ / (3.0 * G_)));
+        g = lambda_dot_ * (sqrt(2.0) * G_ * shear_stress / (sqrt(stress_tensor_J2)));
+    }
+    return shear_stress_rate_elastic - g;
+}
+//=================================================================================================//
+Matd J2Plasticity::ReturnMappingShearStress(Matd &shear_stress, Real &hardening_factor)
+{
+    Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum();
+    Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_);
+    Real r = 1.0;
+    if (f > TinyReal)
+        r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal);
+    return r * shear_stress;
+}
+//=================================================================================================//
+Real J2Plasticity::ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor)
+{
+    Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum();
+    Real r = 1.0;
+    Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_);
+    if (f > TinyReal)
+        r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal);
+    return r;
+}
+//=================================================================================================//
+Real J2Plasticity::HardeningFactorRate(const Matd &shear_stress, Real &hardening_factor)
+{
+    Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum();
+    Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_);
+    return (f > TinyReal) ? 0.5 * f / (G_ + hardening_modulus_ / 3.0) : 0.0;
+}
 } // namespace SPH

src/shared/materials/general_continuum.h

@@ -22,7 +22,7 @@
  * ------------------------------------------------------------------------- */
 /**
  * @file 	general_continuum.h
- * @brief 	Describe the linear elastic and Drucker-Prager's plastic model
+ * @brief 	Describe the linear elastic, J2 plasticity, and Drucker-Prager's plastic model
  * @author	Shuaihao Zhang and Xiangyu Hu
  */
 
@@ -43,7 +43,7 @@ class GeneralContinuum : public WeaklyCompressibleFluid
     Real contact_stiffness_; /* contact-force stiffness related to bulk modulus*/
   public:
     explicit GeneralContinuum(Real rho0, Real c0, Real youngs_modulus, Real poisson_ratio)
-        : WeaklyCompressibleFluid(rho0, c0), E_(0.0), G_(0.0), K_(0.0), nu_(0.0), contact_stiffness_(c0 * c0)
+        : WeaklyCompressibleFluid(rho0, c0), E_(0.0), G_(0.0), K_(0.0), nu_(0.0), contact_stiffness_(rho0_ * c0 * c0)
     {
         material_type_name_ = "GeneralContinuum";
         E_ = youngs_modulus;
@@ -98,5 +98,31 @@ class PlasticContinuum : public GeneralContinuum
 
     virtual GeneralContinuum *ThisObjectPtr() override { return this; };
 };
+
+class J2Plasticity : public GeneralContinuum
+{
+  protected:
+    Real yield_stress_;
+    Real hardening_modulus_;
+    const Real sqrt_2_over_3_ = sqrt(2.0 / 3.0);
+
+  public:
+    explicit J2Plasticity(Real rho0, Real c0, Real youngs_modulus, Real poisson_ratio, Real yield_stress, Real hardening_modulus = 0.0)
+        : GeneralContinuum(rho0, c0, youngs_modulus, poisson_ratio),
+          yield_stress_(yield_stress), hardening_modulus_(hardening_modulus)
+    {
+        material_type_name_ = "J2Plasticity";
+    };
+    virtual ~J2Plasticity(){};
+
+    Real YieldStress() { return yield_stress_; };
+    Real HardeningModulus() { return hardening_modulus_; };
+
+    virtual Matd ConstitutiveRelationShearStress(Matd &velocity_gradient, Matd &shear_stress, Real &hardening_factor);
+    virtual Matd ReturnMappingShearStress(Matd &shear_stress,  Real &hardening_factor);
+    virtual Real ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor);
+    virtual Real HardeningFactorRate(const Matd &shear_stress, Real &hardening_factor);
+    virtual J2Plasticity *ThisObjectPtr() override { return this; };
+};
 } // namespace SPH
 #endif // GENERAL_CONTINUUM_H

src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp

@@ -11,6 +11,25 @@ ContinuumInitialCondition::ContinuumInitialCondition(SPHBody &sph_body)
       vel_(particles_->registerStateVariable<Vecd>("Velocity")),
       stress_tensor_3D_(particles_->registerStateVariable<Mat3d>("StressTensor3D")) {}
 //=================================================================================================//
+AcousticTimeStep::AcousticTimeStep(SPHBody &sph_body, Real acousticCFL)
+    : LocalDynamicsReduce<ReduceMax>(sph_body),
+      fluid_(DynamicCast<Fluid>(this, particles_->getBaseMaterial())),
+      rho_(particles_->getVariableDataByName<Real>("Density")),
+      p_(particles_->getVariableDataByName<Real>("Pressure")),
+      vel_(particles_->getVariableDataByName<Vecd>("Velocity")),
+      smoothing_length_min_(sph_body.sph_adaptation_->MinimumSmoothingLength()),
+      acousticCFL_(acousticCFL) {}
+//=================================================================================================//
+Real AcousticTimeStep::reduce(size_t index_i, Real dt)
+{
+    return fluid_.getSoundSpeed(p_[index_i], rho_[index_i]) + vel_[index_i].norm();
+}
+//=================================================================================================//
+Real AcousticTimeStep::outputResult(Real reduced_value)
+{
+    return acousticCFL_ * smoothing_length_min_ / (reduced_value + TinyReal);
+}
+//=================================================================================================//
 ShearAccelerationRelaxation::ShearAccelerationRelaxation(BaseInnerRelation &inner_relation)
     : fluid_dynamics::BaseIntegration<DataDelegateInner>(inner_relation),
       continuum_(DynamicCast<GeneralContinuum>(this, particles_->getBaseMaterial())),
@@ -97,7 +116,7 @@ void ShearStressRelaxation::update(size_t index_i, Real dt)
 //====================================================================================//
 StressDiffusion::StressDiffusion(BaseInnerRelation &inner_relation)
     : BasePlasticIntegration<DataDelegateInner>(inner_relation),
-      fai_(DynamicCast<PlasticContinuum>(this, plastic_continuum_).getFrictionAngle()),
+      phi_(DynamicCast<PlasticContinuum>(this, plastic_continuum_).getFrictionAngle()),
       smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()),
       sound_speed_(plastic_continuum_.ReferenceSoundSpeed()) {}
 //====================================================================================//
@@ -116,14 +135,116 @@ void StressDiffusion::interaction(size_t index_i, Real dt)
         Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j];
         Real y_ij = pos_[index_i](1, 0) - pos_[index_j](1, 0);
         diffusion_stress_ = stress_tensor_3D_[index_i] - stress_tensor_3D_[index_j];
-        diffusion_stress_(0, 0) -= (1 - sin(fai_)) * density * gravity * y_ij;
+        diffusion_stress_(0, 0) -= (1 - sin(phi_)) * density * gravity * y_ij;
         diffusion_stress_(1, 1) -= density * gravity * y_ij;
-        diffusion_stress_(2, 2) -= (1 - sin(fai_)) * density * gravity * y_ij;
+        diffusion_stress_(2, 2) -= (1 - sin(phi_)) * density * gravity * y_ij;
         diffusion_stress_rate_ += 2 * zeta_ * smoothing_length_ * sound_speed_ *
                                   diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_);
     }
     stress_rate_3D_[index_i] = diffusion_stress_rate_;
 }
 //====================================================================================//
+ShearStressRelaxationHourglassControl ::
+    ShearStressRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi)
+    : fluid_dynamics::BaseIntegration<DataDelegateInner>(inner_relation),
+      continuum_(DynamicCast<GeneralContinuum>(this, particles_->getBaseMaterial())),
+      shear_stress_(particles_->registerStateVariable<Matd>("ShearStress")),
+      shear_stress_rate_(particles_->registerStateVariable<Matd>("ShearStressRate")),
+      velocity_gradient_(particles_->registerStateVariable<Matd>("VelocityGradient")),
+      strain_tensor_(particles_->registerStateVariable<Matd>("StrainTensor")),
+      strain_tensor_rate_(particles_->registerStateVariable<Matd>("StrainTensorRate")),
+      B_(particles_->getVariableDataByName<Matd>("LinearGradientCorrectionMatrix")),
+      von_mises_stress_(particles_->registerStateVariable<Real>("VonMisesStress")),
+      von_mises_strain_(particles_->registerStateVariable<Real>("VonMisesStrain")),
+      scale_penalty_force_(particles_->registerStateVariable<Real>("ScalePenaltyForce", Real(1.0))),
+      acc_shear_(particles_->registerStateVariable<Vecd>("AccelerationByShear")),
+      acc_hourglass_(particles_->registerStateVariable<Vecd>("AccelerationHourglass")),
+      G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())), xi_(xi)
+{
+    particles_->addVariableToSort<Matd>("ShearStress");
+    particles_->addVariableToSort<Matd>("ShearStressRate");
+    particles_->addVariableToSort<Matd>("VelocityGradient");
+    particles_->addVariableToSort<Matd>("StrainTensor");
+    particles_->addVariableToSort<Matd>("StrainTensorRate");
+    particles_->addVariableToSort<Real>("VonMisesStress");
+    particles_->addVariableToSort<Real>("VonMisesStrain");
+    particles_->addVariableToSort<Real>("ScalePenaltyForce");
+    particles_->addVariableToSort<Vecd>("AccelerationByShear");
+    particles_->addVariableToSort<Vecd>("AccelerationHourglass");
+}
+//====================================================================================//
+void ShearStressRelaxationHourglassControl::initialization(size_t index_i, Real dt)
+{
+    Matd velocity_gradient = Matd::Zero();
+    Neighborhood &inner_neighborhood = inner_configuration_[index_i];
+    Vecd vel_i = vel_[index_i];
+    for (size_t n = 0; n != inner_neighborhood.current_size_; ++n)
+    {
+        size_t index_j = inner_neighborhood.j_[n];
+        Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j];
+        Vecd &e_ij = inner_neighborhood.e_ij_[n];
+        Matd velocity_gradient_ij;
+        velocity_gradient_ij = -(vel_i - vel_[index_j]) * (B_[index_i] * e_ij * dW_ijV_j).transpose();
+        velocity_gradient += velocity_gradient_ij;
+    }
+    velocity_gradient_[index_i] = velocity_gradient;
+}
+//====================================================================================//
+void ShearStressRelaxationHourglassControl::interaction(size_t index_i, Real dt)
+{
+    shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]);
+    shear_stress_[index_i] += shear_stress_rate_[index_i] * dt;
+    Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity();
+    von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i);
+    strain_tensor_rate_[index_i] = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose());
+    strain_tensor_[index_i] += strain_tensor_rate_[index_i] * dt;
+    von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]);
+}
+//====================================================================================//
+void ShearStressRelaxationHourglassControl::update(size_t index_i, Real dt)
+{
+    Real rho_i = rho_[index_i];
+    Matd shear_stress_i = shear_stress_[index_i];
+    Vecd acceleration = Vecd::Zero();
+    Vecd acceleration_hourglass = Vecd::Zero();
+    Neighborhood &inner_neighborhood = inner_configuration_[index_i];
+    for (size_t n = 0; n != inner_neighborhood.current_size_; ++n)
+    {
+        size_t index_j = inner_neighborhood.j_[n];
+        Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j];
+        Vecd &e_ij = inner_neighborhood.e_ij_[n];
+        acceleration += ((shear_stress_i + shear_stress_[index_j]) / rho_i) * dW_ijV_j * e_ij;
+        Vecd v_ij = vel_[index_i] - vel_[index_j];
+        Real r_ij = inner_neighborhood.r_ij_[n];
+        Vecd v_ij_correction = v_ij - 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_j]) * r_ij * e_ij;
+        acceleration_hourglass += 0.5 * xi_ * (scale_penalty_force_[index_i] + scale_penalty_force_[index_j]) * G_ * v_ij_correction * dW_ijV_j * dt / (rho_i * r_ij);
+    }
+    acc_hourglass_[index_i] += acceleration_hourglass;
+    acc_shear_[index_i] = acceleration + acc_hourglass_[index_i];
+}
+//====================================================================================//
+ShearStressRelaxationHourglassControlJ2Plasticity ::
+    ShearStressRelaxationHourglassControlJ2Plasticity(BaseInnerRelation &inner_relation, Real xi)
+    : ShearStressRelaxationHourglassControl(inner_relation, xi),
+      J2_plasticity_(DynamicCast<J2Plasticity>(this, particles_->getBaseMaterial())),
+      hardening_factor_(particles_->registerStateVariable<Real>("HardeningFactor"))
+{
+    particles_->addVariableToSort<Real>("HardeningFactor");
+}
+//====================================================================================//
+void ShearStressRelaxationHourglassControlJ2Plasticity::interaction(size_t index_i, Real dt)
+{
+    shear_stress_rate_[index_i] = J2_plasticity_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i], hardening_factor_[index_i]);
+    Matd shear_stress_try = shear_stress_[index_i] + shear_stress_rate_[index_i] * dt;
+    Real hardening_factor_increment = J2_plasticity_.HardeningFactorRate(shear_stress_try, hardening_factor_[index_i]);
+    hardening_factor_[index_i] += sqrt(2.0 / 3.0) * hardening_factor_increment;
+    scale_penalty_force_[index_i] = J2_plasticity_.ScalePenaltyForce(shear_stress_try, hardening_factor_[index_i]);
+    shear_stress_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_try, hardening_factor_[index_i]);
+    Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity();
+    von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i);
+    Matd strain_rate = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose());
+    strain_tensor_[index_i] += strain_rate * dt;
+    von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]);
+}
 } // namespace continuum_dynamics
-} // namespace SPH
+} // namespace SPH

src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h

@@ -49,6 +49,22 @@ class ContinuumInitialCondition : public LocalDynamics
     Mat3d *stress_tensor_3D_;
 };
 
+class AcousticTimeStep : public LocalDynamicsReduce<ReduceMax>
+{
+  public:
+    explicit AcousticTimeStep(SPHBody &sph_body, Real acousticCFL = 0.6);
+    virtual ~AcousticTimeStep(){};
+    Real reduce(size_t index_i, Real dt = 0.0);
+    virtual Real outputResult(Real reduced_value) override;
+
+  protected:
+    Fluid &fluid_;
+    Real *rho_, *p_;
+    Vecd *vel_;
+    Real smoothing_length_min_;
+    Real acousticCFL_;
+};
+
 template <class FluidDynamicsType>
 class BaseIntegration1stHalf : public FluidDynamicsType
 {
@@ -200,9 +216,38 @@ class StressDiffusion : public BasePlasticIntegration<DataDelegateInner>
     void interaction(size_t index_i, Real dt = 0.0);
 
   protected:
-    Real zeta_ = 0.1, fai_; /*diffusion coefficient*/
+    Real zeta_ = 0.1, phi_; /*diffusion coefficient*/
     Real smoothing_length_, sound_speed_;
 };
+
+class ShearStressRelaxationHourglassControl : public fluid_dynamics::BaseIntegration<DataDelegateInner>
+{
+  public:
+    explicit ShearStressRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi_ = 4.0);
+    virtual ~ShearStressRelaxationHourglassControl(){};
+    void initialization(size_t index_i, Real dt = 0.0);
+    void interaction(size_t index_i, Real dt = 0.0);
+    void update(size_t index_i, Real dt = 0.0);
+
+  protected:
+    GeneralContinuum &continuum_;
+    Matd *shear_stress_, *shear_stress_rate_, *velocity_gradient_, *strain_tensor_, *strain_tensor_rate_, *B_;
+    Real *von_mises_stress_, *von_mises_strain_, *scale_penalty_force_;
+    Vecd *acc_shear_, *acc_hourglass_;
+    Real G_, xi_;
+};
+
+class ShearStressRelaxationHourglassControlJ2Plasticity : public ShearStressRelaxationHourglassControl
+{
+  public:
+    explicit ShearStressRelaxationHourglassControlJ2Plasticity(BaseInnerRelation &inner_relation, Real xi_ = 0.2);
+    virtual ~ShearStressRelaxationHourglassControlJ2Plasticity(){};
+    void interaction(size_t index_i, Real dt = 0.0);
+
+  protected:
+    J2Plasticity &J2_plasticity_;
+    Real *hardening_factor_;
+};
 } // namespace continuum_dynamics
 } // namespace SPH
 #endif // CONTINUUM_INTEGRATION_H

tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp

@@ -26,7 +26,6 @@ BoundingBox system_domain_bounds(Vec2d(-SL - BW, -PL / 2.0),
 Real rho0_s = 1.0e3;         // reference density
 Real Youngs_modulus = 2.0e6; // reference Youngs modulus
 Real poisson = 0.3975;       // Poisson ratio
-// Real poisson = 0.4;		 //Poisson ratio
 Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s));
 Real gravity_g = 0.0;
 //----------------------------------------------------------------------

tests/2d_examples/test_2d_oscillating_beam_ULGNOG/CMakeLists.txt

@@ -0,0 +1,21 @@
+STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR})
+PROJECT("${CURRENT_FOLDER}")
+
+SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
+SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/")
+SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input")
+SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload")
+
+file(MAKE_DIRECTORY ${BUILD_INPUT_PATH})
+file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/
+	DESTINATION ${BUILD_INPUT_PATH})
+
+add_executable(${PROJECT_NAME})
+aux_source_directory(. DIR_SRCS)
+target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS})
+target_link_libraries(${PROJECT_NAME} sphinxsys_2d)
+set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}")
+
+add_test(NAME ${PROJECT_NAME}
+	COMMAND ${PROJECT_NAME} --state_recording=${TEST_STATE_RECORDING}
+	WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})