Generalized non-hourglass formulation for ULSPH solid dynamics

src/shared/materials/general_continuum.cpp

@@ -43,16 +43,15 @@ Mat3d PlasticContinuum::ConstitutiveRelation(Mat3d &velocity_gradient, Mat3d &st
     Mat3d spin_rate = 0.5 * (velocity_gradient - velocity_gradient.transpose());
     Mat3d deviatoric_strain_rate = strain_rate - (1.0 / stress_dimension_) * strain_rate.trace() * Mat3d::Identity();
     Mat3d stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate + K_ * strain_rate.trace() * Mat3d::Identity() + stress_tensor * (spin_rate.transpose()) + spin_rate * stress_tensor;
-    Real 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();
-    Real f = sqrt(stress_tensor_J2) + alpha_phi_ * stress_tensor_I1 - k_c_;
+    Real f = sqrt(stress_tensor_J2) + alpha_phi_ * stress_tensor.trace() - k_c_;
     Real lambda_dot_ = 0;
     Mat3d g = Mat3d::Zero();
     if (f >= TinyReal)
     {
         Real deviatoric_stress_times_strain_rate = (deviatoric_stress_tensor.cwiseProduct(strain_rate)).sum();
-        //non-associate flow rule
+        // non-associate flow rule
         lambda_dot_ = (3.0 * alpha_phi_ * K_ * strain_rate.trace() + (G_ / sqrt(stress_tensor_J2)) * deviatoric_stress_times_strain_rate) / (9.0 * alpha_phi_ * K_ * getDPConstantsA(psi_) + G_);
         g = lambda_dot_ * (3.0 * K_ * getDPConstantsA(psi_) * Mat3d::Identity() + G_ * deviatoric_stress_tensor / (sqrt(stress_tensor_J2)));
     }
@@ -64,9 +63,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 +74,46 @@ 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 f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_);
+    return (f > TinyReal) ? (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal) : 1.0;
+}
+//=================================================================================================//
+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/all_continuum_dynamics.h

@@ -30,3 +30,4 @@
 
 #include "base_continuum_dynamics.h"
 #include "continuum_integration.hpp"
+#include "continuum_dynamics_variable.h"

src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.cpp

@@ -0,0 +1,57 @@
+#include "continuum_dynamics_variable.h"
+
+namespace SPH
+{
+namespace continuum_dynamics
+{
+//=============================================================================================//
+VonMisesStress::VonMisesStress(SPHBody &sph_body)
+    : BaseDerivedVariable<Real>(sph_body, "VonMisesStress"),
+      p_(particles_->getVariableDataByName<Real>("Pressure")),
+      shear_stress_(particles_->getVariableDataByName<Matd>("ShearStress")) {}
+//=============================================================================================//
+void VonMisesStress::update(size_t index_i, Real dt)
+{
+    Matd stress_tensor = shear_stress_[index_i] - p_[index_i] * Matd::Identity();
+    derived_variable_[index_i] = getVonMisesStressFromMatrix(stress_tensor);
+}
+//=============================================================================================//
+VonMisesStrain::VonMisesStrain(SPHBody &sph_body)
+    : BaseDerivedVariable<Real>(sph_body, "VonMisesStrain"),
+      strain_tensor_(particles_->getVariableDataByName<Matd>("StrainTensor")) {}
+//=============================================================================================//
+void VonMisesStrain::update(size_t index_i, Real dt)
+{
+    derived_variable_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]);
+}
+//=============================================================================================//
+VerticalStress::VerticalStress(SPHBody &sph_body)
+    : BaseDerivedVariable<Real>(sph_body, "VerticalStress"),
+      stress_tensor_3D_(particles_->getVariableDataByName<Mat3d>("StressTensor3D")) {}
+//=============================================================================================//
+void VerticalStress::update(size_t index_i, Real dt)
+{
+    derived_variable_[index_i] = stress_tensor_3D_[index_i](1, 1);
+}
+//=============================================================================================//
+AccDeviatoricPlasticStrain::AccDeviatoricPlasticStrain(SPHBody &sph_body)
+    : BaseDerivedVariable<Real>(sph_body, "AccDeviatoricPlasticStrain"),
+      plastic_continuum_(DynamicCast<PlasticContinuum>(this, sph_body_.getBaseMaterial())),
+      stress_tensor_3D_(particles_->getVariableDataByName<Mat3d>("StressTensor3D")),
+      strain_tensor_3D_(particles_->getVariableDataByName<Mat3d>("StrainTensor3D")),
+      E_(plastic_continuum_.getYoungsModulus()), nu_(plastic_continuum_.getPoissonRatio()) {}
+//=============================================================================================//
+void AccDeviatoricPlasticStrain::update(size_t index_i, Real dt)
+{
+    Mat3d deviatoric_stress = stress_tensor_3D_[index_i] - (1.0 / 3.0) * stress_tensor_3D_[index_i].trace() * Mat3d::Identity();
+    Real hydrostatic_pressure = (1.0 / 3.0) * stress_tensor_3D_[index_i].trace();
+    Mat3d elastic_strain_tensor_3D = deviatoric_stress / (2.0 * plastic_continuum_.getShearModulus(E_, nu_)) +
+                                     hydrostatic_pressure * Mat3d::Identity() / (9.0 * plastic_continuum_.getBulkModulus(E_, nu_));
+    Mat3d plastic_strain_tensor_3D = strain_tensor_3D_[index_i] - elastic_strain_tensor_3D;
+    Mat3d deviatoric_strain_tensor = plastic_strain_tensor_3D - (1.0 / (Real)Dimensions) * plastic_strain_tensor_3D.trace() * Mat3d::Identity();
+    Real sum = (deviatoric_strain_tensor.cwiseProduct(deviatoric_strain_tensor)).sum();
+    derived_variable_[index_i] = sqrt(sum * 2.0 / 3.0);
+}
+//=================================================================================================//
+} // namespace continuum_dynamics
+} // namespace SPH

src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.h

@@ -0,0 +1,99 @@
+/* ------------------------------------------------------------------------- *
+ *                                SPHinXsys                                  *
+ * ------------------------------------------------------------------------- *
+ * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle *
+ * Hydrodynamics for industrial compleX systems. It provides C++ APIs for    *
+ * physical accurate simulation and aims to model coupled industrial dynamic *
+ * systems including fluid, solid, multi-body dynamics and beyond with SPH   *
+ * (smoothed particle hydrodynamics), a meshless computational method using  *
+ * particle discretization.                                                  *
+ *                                                                           *
+ * SPHinXsys is partially funded by German Research Foundation               *
+ * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1,            *
+ *  HU1527/12-1 and HU1527/12-4.                                             *
+ *                                                                           *
+ * Portions copyright (c) 2017-2023 Technical University of Munich and       *
+ * the authors' affiliations.                                                *
+ *                                                                           *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may   *
+ * not use this file except in compliance with the License. You may obtain a *
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.        *
+ *                                                                           *
+ * ------------------------------------------------------------------------- */
+/**
+ * @file continuum_dynamics_variable.h
+ * @brief Here, we define the algorithm classes for computing derived solid dynamics variables.
+ * @details These variable can be added into variable list for state output.
+ * @author Shuaihao Zhang and Xiangyu Hu
+ */
+
+#ifndef CONTINUUM_DYNAMICS_VARIABLE_H
+#define CONTINUUM_DYNAMICS_VARIABLE_H
+
+#include "base_general_dynamics.h"
+#include "general_continuum.h"
+
+namespace SPH
+{
+namespace continuum_dynamics
+{
+/**
+ * @class VonMisesStress
+ * @brief computing von_Mises_stress
+ */
+class VonMisesStress : public BaseDerivedVariable<Real>
+{
+  public:
+    explicit VonMisesStress(SPHBody &sph_body);
+    virtual ~VonMisesStress(){};
+    void update(size_t index_i, Real dt = 0.0);
+
+  protected:
+    Real *p_;
+    Matd *shear_stress_;
+};
+/**
+ * @class VonMisesStrain
+ * @brief computing von_Mises_strain
+ */
+class VonMisesStrain : public BaseDerivedVariable<Real>
+{
+  public:
+    explicit VonMisesStrain(SPHBody &sph_body);
+    virtual ~VonMisesStrain(){};
+    void update(size_t index_i, Real dt = 0.0);
+
+  protected:
+    Matd *strain_tensor_;
+};
+/**
+ * @class VerticalStress
+ */
+class VerticalStress : public BaseDerivedVariable<Real>
+{
+  public:
+    explicit VerticalStress(SPHBody &sph_body);
+    virtual ~VerticalStress(){};
+    void update(size_t index_i, Real dt = 0.0);
+
+  protected:
+    Mat3d *stress_tensor_3D_;
+};
+/**
+ * @class AccumulatedDeviatoricPlasticStrain
+ */
+class AccDeviatoricPlasticStrain : public BaseDerivedVariable<Real>
+{
+  public:
+    explicit AccDeviatoricPlasticStrain(SPHBody &sph_body);
+    virtual ~AccDeviatoricPlasticStrain(){};
+    void update(size_t index_i, Real dt = 0.0);
+
+  protected:
+    PlasticContinuum &plastic_continuum_;
+    Mat3d *stress_tensor_3D_, *strain_tensor_3D_;
+    Real E_, nu_;
+};
+} // namespace continuum_dynamics
+} // namespace SPH
+#endif // CONTINUUM_DYNAMICS_VARIABLE_H