add plastic materials

tests/user_examples/extra_src/shared/general_continuum/all_continuum.h

@@ -34,3 +34,4 @@
 #include "all_continuum_dynamics.h"
 #include "general_continuum.h"
 #include "continuum_particles.h"
+#include "riemann_solver_extra.h"

tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h

@@ -17,6 +17,11 @@ typedef DataDelegateContact<ContinuumParticles, BaseParticles, DataDelegateEmpty
     ContinuumContactData;
 typedef DataDelegateContact<ContinuumParticles, SolidParticles> FSIContactData;
 
+typedef DataDelegateContact<PlasticContinuumParticles, SolidParticles, DataDelegateEmptyBase>
+PlasticContinuumWallData;
+typedef DataDelegateContact<PlasticContinuumParticles, BaseParticles, DataDelegateEmptyBase>
+PlasticContinuumContactData;
+typedef DataDelegateContact<PlasticContinuumParticles, SolidParticles> PlasticFSIContactData;
 /**
  * @class BaseShearStressRelaxation1stHalfType
  */
@@ -51,6 +56,52 @@ class BaseShearStressRelaxation2ndHalfWithWall : public fluid_dynamics::Interact
 };
 
 using ShearStressRelaxation2ndHalfWithWall = BaseShearStressRelaxation2ndHalfWithWall<ShearStressRelaxation2ndHalf>;
+
+template <class BaseStressRelaxation1stHalfType>
+class BaseStressRelaxation1stHalfWithWall : public fluid_dynamics::InteractionWithWall<BaseStressRelaxation1stHalfType>
+{
+public:
+    template <typename... Args>
+    BaseStressRelaxation1stHalfWithWall(Args &&...args)
+        : fluid_dynamics::InteractionWithWall<BaseStressRelaxation1stHalfType>(std::forward<Args>(args)...) {};
+    virtual ~BaseStressRelaxation1stHalfWithWall() {};
+    void interaction(size_t index_i, Real dt = 0.0);
+
+protected:
+    virtual Vecd computeNonConservativeAcceleration(size_t index_i) override;
+};
+using StressRelaxation1stHalfWithWall = BaseStressRelaxation1stHalfWithWall<StressRelaxation1stHalf>;
+using StressRelaxation1stHalfRiemannWithWall = BaseStressRelaxation1stHalfWithWall<StressRelaxation1stHalfRiemann>;
+using StressRelaxation1stHalfDissipativeRiemannfWithWall = BaseStressRelaxation1stHalfWithWall<StressRelaxation1stHalfDissipativeRiemann>;
+
+template <class BaseStressRelaxation2ndHalfType>
+class BaseStressRelaxation2ndHalfWithWall : public fluid_dynamics::InteractionWithWall<BaseStressRelaxation2ndHalfType>
+{
+public:
+    template <typename... Args>
+    BaseStressRelaxation2ndHalfWithWall(Args &&...args)
+        : fluid_dynamics::InteractionWithWall<BaseStressRelaxation2ndHalfType>(std::forward<Args>(args)...) {};
+    virtual ~BaseStressRelaxation2ndHalfWithWall() {};
+    void interaction(size_t index_i, Real dt = 0.0);
+};
+using StressRelaxation2ndHalfWithWall = BaseStressRelaxation2ndHalfWithWall<StressRelaxation2ndHalf>;
+using StressRelaxation2ndHalfRiemannWithWall = BaseStressRelaxation2ndHalfWithWall<StressRelaxation2ndHalfRiemann>;
+using StressRelaxation2ndHalfDissipativeRiemannWithWall = BaseStressRelaxation2ndHalfWithWall<StressRelaxation2ndHalfDissipativeRiemann>;
+
+/**
+* @class StressDiffusionWithWall
+*/
+template <class BaseStressDiffusionType>
+class BaseStressDiffusionWithWall : public fluid_dynamics::InteractionWithWall<BaseStressDiffusionType>
+{
+public:
+    template <typename... Args>
+    BaseStressDiffusionWithWall(Args &&...args)
+        : fluid_dynamics::InteractionWithWall<BaseStressDiffusionType>(std::forward<Args>(args)...) {};
+    virtual ~BaseStressDiffusionWithWall() {};
+    void interaction(size_t index_i, Real dt = 0.0);
+};
+using StressDiffusionWithWall = BaseStressDiffusionWithWall<StressDiffusion>;
 } // namespace continuum_dynamics
 } // namespace SPH
 

tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp

@@ -70,5 +70,114 @@ void BaseShearStressRelaxation2ndHalfWithWall<BaseShearStressRelaxation2ndHalfTy
     }
     this->velocity_gradient_[index_i] += velocity_gradient;
 }
+
+template <class BaseStressRelaxation1stHalfType>
+Vecd BaseStressRelaxation1stHalfWithWall<BaseStressRelaxation1stHalfType>::computeNonConservativeAcceleration(size_t index_i)
+{
+    return this->acc_prior_[index_i];
+}
+
+template <class BaseStressRelaxation1stHalfType>
+void BaseStressRelaxation1stHalfWithWall<BaseStressRelaxation1stHalfType>::interaction(size_t index_i, Real dt)
+{
+    BaseStressRelaxation1stHalfType::interaction(index_i, dt);
+
+    Vecd acc_prior_i = computeNonConservativeAcceleration(index_i);
+    Vecd acceleration = acc_prior_i;
+    Real rho_dissipation(0);
+
+    Matd stress_tensor_i = this->reduceTensor(this->stress_tensor_3D_[index_i]);
+
+    for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k)
+    {
+        StdLargeVec<Vecd>& acc_ave_k = *(this->wall_acc_ave_[k]);
+        Neighborhood& wall_neighborhood = (*fluid_dynamics::FluidWallData::contact_configuration_[k])[index_i];
+        for (size_t n = 0; n != wall_neighborhood.current_size_; ++n)
+        {
+            size_t index_j = wall_neighborhood.j_[n];
+            Vecd& e_ij = wall_neighborhood.e_ij_[n];
+            Real dW_ijV_j = wall_neighborhood.dW_ijV_j_[n];
+            Real r_ij = wall_neighborhood.r_ij_[n];
+
+            Real face_wall_external_acceleration = (acc_prior_i - acc_ave_k[index_j]).dot(-e_ij);
+            Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(Real(0), face_wall_external_acceleration);
+            acceleration += 2 * stress_tensor_i * wall_neighborhood.dW_ijV_j_[n] * wall_neighborhood.e_ij_[n];
+            rho_dissipation += this->riemann_solver_.DissipativeUJump(this->p_[index_i] - p_in_wall) * dW_ijV_j;
+        }
+    }
+    this->acc_[index_i] += acceleration / this->rho_[index_i];
+    this->drho_dt_[index_i] += rho_dissipation * this->rho_[index_i];
+}
+//=================================================================================================//
+template <class BaseStressRelaxation2ndHalfType>
+void BaseStressRelaxation2ndHalfWithWall<BaseStressRelaxation2ndHalfType>::interaction(size_t index_i, Real dt)
+{
+    BaseStressRelaxation2ndHalfType::interaction(index_i, dt);
+
+    Real density_change_rate = 0.0;
+    Vecd p_dissipation = Vecd::Zero();
+
+    Vecd vel_i = this->vel_[index_i];
+    Matd velocity_gradient = Matd::Zero();
+
+    for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k)
+    {
+        StdLargeVec<Vecd>& vel_ave_k = *(this->wall_vel_ave_[k]);
+        StdLargeVec<Vecd>& n_k = *(this->wall_n_[k]);
+
+        Neighborhood& wall_neighborhood = (*fluid_dynamics::FluidWallData::contact_configuration_[k])[index_i];
+        for (size_t n = 0; n != wall_neighborhood.current_size_; ++n)
+        {
+            size_t index_j = wall_neighborhood.j_[n];
+            Vecd& e_ij = wall_neighborhood.e_ij_[n];
+            Real dW_ijV_j = wall_neighborhood.dW_ijV_j_[n];
+            Vecd nablaW_ijV_j = wall_neighborhood.dW_ijV_j_[n] * wall_neighborhood.e_ij_[n];
+            Vecd vel_in_wall = 1.9 * vel_ave_k[index_j] - 0.9 *vel_i;
+            Matd velocity_gradient_ij = -(vel_i - vel_in_wall) * nablaW_ijV_j.transpose();
+            velocity_gradient += velocity_gradient_ij;
+
+            density_change_rate += (this->vel_[index_i] - vel_in_wall).dot(e_ij) * dW_ijV_j;
+            Vecd u_jump = 2.0 * (this->vel_[index_i] - vel_ave_k[index_j]);
+            p_dissipation += this->riemann_solver_.DissipativePJumpExtra(u_jump, n_k[index_j]) * dW_ijV_j;
+        }
+    }
+    this->drho_dt_[index_i] += density_change_rate * this->rho_[index_i];
+    this->velocity_gradient_[index_i] += velocity_gradient;
+    this->acc_[index_i] += p_dissipation / this->rho_[index_i];
+}
+
+template <class BaseStressDiffusionType>
+void BaseStressDiffusionWithWall<BaseStressDiffusionType>::interaction(size_t index_i, Real dt)
+{
+    BaseStressDiffusionType::interaction(index_i, dt);
+
+    Real rho_in_wall = this->plastic_continuum_.getDensity();
+    Real gravity = abs(this->acc_prior_[index_i](1, 0));
+    Real density = rho_in_wall;
+    Mat3d diffusion_stress_rate_ = Mat3d::Zero();
+    Mat3d diffusion_stress_ = Mat3d::Zero();
+    Mat3d stress_tensor_i = this->stress_tensor_3D_[index_i];
+    for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k)
+    {
+        Neighborhood& wall_neighborhood = (*fluid_dynamics::FluidWallData::contact_configuration_[k])[index_i];
+        for (size_t n = 0; n != wall_neighborhood.current_size_; ++n)
+        {
+            size_t index_j = wall_neighborhood.j_[n];
+            Real r_ij = wall_neighborhood.r_ij_[n];
+            Real dW_ijV_j = wall_neighborhood.dW_ijV_j_[n];
+
+            Real y_ij = this->pos_[index_i](1, 0) - this->pos_[index_j](1, 0);
+            //stress boundary condition
+            Mat3d stress_tensor_j = stress_tensor_i;
+            diffusion_stress_ = stress_tensor_i - stress_tensor_j;
+            diffusion_stress_(0, 0) = diffusion_stress_(0, 0) - (1 - sin(this->fai_)) * density * gravity * y_ij;
+            diffusion_stress_(1, 1) = diffusion_stress_(1, 1) - density * gravity * y_ij;
+            diffusion_stress_(2, 2) = diffusion_stress_(2, 2) - (1 - sin(this->fai_)) * density * gravity * y_ij;
+            diffusion_stress_rate_ += 2 * this->zeta_ * this->smoothing_length_ * this->sound_speed_ * diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * this->smoothing_length_);
+        }
+    }
+    this->stress_rate_3D_[index_i] += diffusion_stress_rate_;
+}
+
 } // namespace continuum_dynamics
 } // namespace SPH

tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp

@@ -5,8 +5,8 @@ namespace SPH
     namespace continuum_dynamics
     {
         ContinuumInitialCondition::ContinuumInitialCondition(SPHBody &sph_body)
-            : LocalDynamics(sph_body), ContinuumDataSimple(sph_body),
-              pos_(particles_->pos_), vel_(particles_->vel_) {}
+            : LocalDynamics(sph_body), PlasticContinuumDataSimple(sph_body),
+              pos_(particles_->pos_), vel_(particles_->vel_), stress_tensor_3D_(particles_->stress_tensor_3D_) {}
         //=================================================================================================//
         ContinuumAcousticTimeStepSize::ContinuumAcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL)
             : fluid_dynamics::AcousticTimeStepSize(sph_body, acousticCFL) {}
@@ -245,7 +245,6 @@ namespace SPH
             Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity();
             von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i);
         }
-
         //=================================================================================================//
         FixedInAxisDirection::FixedInAxisDirection(BodyPartByParticle& body_part, Vecd constrained_axises)
             : BaseMotionConstraint<BodyPartByParticle>(body_part), constrain_matrix_(Matd::Identity())
@@ -258,7 +257,6 @@ namespace SPH
         {
             vel_[index_i] = constrain_matrix_ * vel_[index_i];
         };
-
         //=================================================================================================//
         ConstrainSolidBodyMassCenter::
             ConstrainSolidBodyMassCenter(SPHBody& sph_body, Vecd constrain_direction)
@@ -282,5 +280,67 @@ namespace SPH
         {
             vel_[index_i] -= velocity_correction_;
         }
+        //=================================================================================================//
+        //================================================Plastic==========================================//
+        //=================================================================================================//
+        BaseRelaxationPlastic::BaseRelaxationPlastic(BaseInnerRelation& inner_relation)
+            : LocalDynamics(inner_relation.getSPHBody()), PlasticContinuumDataInner(inner_relation),
+            plastic_continuum_(particles_->plastic_continuum_), rho_(particles_->rho_),
+            p_(*particles_->getVariableByName<Real>("Pressure")), drho_dt_(*particles_->registerSharedVariable<Real>("DensityChangeRate")), pos_(particles_->pos_), vel_(particles_->vel_), acc_(particles_->acc_), acc_prior_(particles_->acc_prior_),
+            stress_tensor_3D_(particles_->stress_tensor_3D_), strain_tensor_3D_(particles_->strain_tensor_3D_),
+            stress_rate_3D_(particles_->stress_rate_3D_), strain_rate_3D_(particles_->strain_rate_3D_),
+            elastic_strain_tensor_3D_(particles_->elastic_strain_tensor_3D_), elastic_strain_rate_3D_(particles_->elastic_strain_rate_3D_) {}
+        Matd BaseRelaxationPlastic::reduceTensor(Mat3d tensor_3d)
+        {
+            Matd tensor_2d;
+            for (int i = 0; i < (Real)Dimensions; i++)
+            {
+                for (int j = 0; j < (Real)Dimensions; j++)
+                {
+                    tensor_2d(i, j) = tensor_3d(i, j);
+                }
+            }
+            return tensor_2d;
+        }
+        Mat3d BaseRelaxationPlastic::increaseTensor(Matd tensor_2d)
+        {
+            Mat3d tensor_3d = Mat3d::Zero();
+            for (int i = 0; i < (Real)Dimensions; i++)
+            {
+                for (int j = 0; j < (Real)Dimensions; j++)
+                {
+                    tensor_3d(i, j) = tensor_2d(i, j);
+                }
+            }
+            return tensor_3d;
+        }
+        //====================================================================================//
+        //===============================StressDiffusion======================================//
+        //====================================================================================//
+        StressDiffusion::StressDiffusion(BaseInnerRelation& inner_relation)
+            : BaseRelaxationPlastic(inner_relation), fai_(DynamicCast<PlasticContinuum>(this, plastic_continuum_).getFrictionAngle()), smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()),
+            sound_speed_(plastic_continuum_.ReferenceSoundSpeed()) {}
+        void StressDiffusion::interaction(size_t index_i, Real dt)
+        {
+            Real gravity = abs(acc_prior_[index_i](1, 0));
+            Real density = plastic_continuum_.getDensity();
+            Mat3d diffusion_stress_rate_ = Mat3d::Zero();
+            Mat3d diffusion_stress_ = Mat3d::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 r_ij = inner_neighborhood.r_ij_[n];
+                Real dW_ijV_j = inner_neighborhood.dW_ijV_j_[n];
+                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) = diffusion_stress_(0, 0) - (1 - sin(fai_)) * density * gravity * y_ij;
+                diffusion_stress_(1, 1) = diffusion_stress_(1, 1) - density * gravity * y_ij;
+                diffusion_stress_(2, 2) = diffusion_stress_(2, 2) - (1 - sin(fai_)) * 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_;
+        }
     } // namespace continuum_dynamics
 } // namespace SPH