non-hourglass formulation for updated Lagrangian SPH

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

@@ -27,7 +27,7 @@
  * @details The fluid dynamics algorithms begin for fluid bulk without boundary condition,
  * then algorithm interacting with wall is defined, further algorithms for multiphase flow interaction
  * built upon these basic algorithms.
- * @author	Chi Zhang and Xiangyu Hu
+ * @author	Chi ZHang and Xiangyu Hu
  */
 #pragma once
 

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

@@ -27,7 +27,7 @@
  * @details The fluid dynamics algorithms begin for fluid bulk without boundary condition,
  * then algorithm interacting with wall is defined, further algorithms for multiphase flow interaction
  * built upon these basic algorithms.
- * @author	Chi Zhang and Xiangyu Hu
+ * @author	Chi ZHang and Xiangyu Hu
  */
 #pragma once
 

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

@@ -3,7 +3,7 @@
 namespace SPH
 {
 //=================================================================================================//
-	namespace continuum_dyannmics
+	namespace continuum_dynamics
 	{
 
 	}

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

@@ -8,7 +8,7 @@
 
 namespace SPH
 {
-namespace continuum_dyannmics
+namespace continuum_dynamics
 {
 
 typedef DataDelegateContact<ContinuumParticles, SolidParticles, DataDelegateEmptyBase>
@@ -51,7 +51,7 @@ class BaseShearStressRelaxation2ndHalfWithWall : public fluid_dynamics::Interact
 };
 
 using ShearStressRelaxation2ndHalfWithWall = BaseShearStressRelaxation2ndHalfWithWall<ShearStressRelaxation2ndHalf>;
-} // namespace continuum_dyannmics
+} // namespace continuum_dynamics
 } // namespace SPH
 
 #endif // CONTINUUM_DYNAMICS_COMPLEX_H

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

@@ -5,7 +5,7 @@
 namespace SPH
 {
 //=================================================================================================//
-namespace continuum_dyannmics
+namespace continuum_dynamics
 {
 //=================================================================================================//
 //==========================BaseShearStressRelaxation1stHalfWithWall================================//
@@ -16,7 +16,7 @@ void BaseShearStressRelaxation1stHalfWithWall<BaseShearStressRelaxation1stHalfTy
     BaseShearStressRelaxation1stHalfType::interaction(index_i, dt);
     Matd shear_stress_i = this->shear_stress_[index_i];
     Real rho_i = this->rho_[index_i];
-    Real rho_in_wall = this->granular_material_.getDensity();
+    Real rho_in_wall = this->continuum_.getDensity();
     Vecd acceleration = Vecd::Zero();
     for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k) // There may be several wall bodies.
     {
@@ -48,7 +48,7 @@ void BaseShearStressRelaxation2ndHalfWithWall<BaseShearStressRelaxation2ndHalfTy
 
     Vecd vel_i = this->vel_[index_i];
     Matd velocity_gradient = Matd::Zero();
-    Real rho_in_wall = this->granular_material_.getDensity();
+    Real rho_in_wall = this->continuum_.getDensity();
     for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k)
     {
 
@@ -70,5 +70,5 @@ void BaseShearStressRelaxation2ndHalfWithWall<BaseShearStressRelaxation2ndHalfTy
     }
     this->velocity_gradient_[index_i] += velocity_gradient;
 }
-} // namespace continuum_dyannmics
+} // namespace continuum_dynamics
 } // namespace SPH

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

@@ -2,7 +2,7 @@
 namespace SPH
 {
 //=================================================================================================//
-    namespace continuum_dyannmics
+    namespace continuum_dynamics
     {
         ContinuumInitialCondition::ContinuumInitialCondition(SPHBody &sph_body)
             : LocalDynamics(sph_body), ContinuumDataSimple(sph_body),
@@ -23,7 +23,7 @@ namespace SPH
         //=================================================================================================//
         BaseRelaxation::BaseRelaxation(BaseInnerRelation &inner_relation)
             : LocalDynamics(inner_relation.getSPHBody()), ContinuumDataInner(inner_relation),
-              granular_material_(particles_->granular_material_), rho_(particles_->rho_),
+              continuum_(particles_->continuum_), rho_(particles_->rho_),
               p_(*particles_->getVariableByName<Real>("Pressure")), drho_dt_(*particles_->getVariableByName<Real>("DensityChangeRate")),
               pos_(particles_->pos_), vel_(particles_->vel_),
               acc_(particles_->acc_), acc_prior_(particles_->acc_prior_) {}
@@ -118,6 +118,7 @@ namespace SPH
                 acceleration += rho_[index_j] * (shear_stress_i / (rho_i * rho_i) + shear_stress_[index_j] / (rho_[index_j] * rho_[index_j])) * nablaW_ijV_j;
             }
             acc_shear_[index_i] += acceleration; // for with artificial stress
+            //acc_shear_[index_i] = acceleration; // for without artificial stress
         }
 
         void ShearStressRelaxation1stHalf::update(size_t index_i, Real dt)
@@ -151,10 +152,138 @@ namespace SPH
 
         void ShearStressRelaxation2ndHalf::update(size_t index_i, Real dt)
         {
-            shear_stress_rate_[index_i] = granular_material_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]);
+            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 * 0.5;
             Matd stress_tensor_i = shear_stress_[index_i] + p_[index_i] * Matd::Identity();
             von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i);
         }
-    } // namespace continuum_dyannmics
+
+        //=================================================================================================//
+        //===================================Non-hourglass formulation=====================================//
+        //=================================================================================================//
+        ShearAccelerationRelaxation::ShearAccelerationRelaxation(BaseInnerRelation& inner_relation)
+            : BaseRelaxation(inner_relation),
+            G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())),
+            smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), shear_stress_(particles_->shear_stress_),
+            B_(*this->particles_->template registerSharedVariable<Matd>("CorrectionMatrix", Matd::Identity())), acc_shear_(particles_->acc_shear_) {}
+        //original
+        void ShearAccelerationRelaxation::interaction(size_t index_i, Real dt)
+        {
+            Real rho_i = rho_[index_i];
+            Vecd acceleration = Vecd::Zero();
+            Vecd vel_i = vel_[index_i];
+            Neighborhood& inner_neighborhood = inner_configuration_[index_i];
+            //Matd B_i = B_[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];
+
+                Vecd v_ij = (vel_i - vel_[index_j]);
+                acceleration += 2 * v_ij * dW_ijV_j / (r_ij + 0.01 * smoothing_length_);
+            }
+            acc_shear_[index_i] += G_ * acceleration * dt / rho_i;
+        }
+        //=================================================================================================//
+        void AngularConservativeShearAccelerationRelaxation::interaction(size_t index_i, Real dt)
+        {
+            Real rho_i = rho_[index_i];
+            Vecd acceleration = 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 r_ij = inner_neighborhood.r_ij_[n];
+                Real dW_ijV_j = inner_neighborhood.dW_ijV_j_[n];
+                Vecd& e_ij = inner_neighborhood.e_ij_[n];
+                Real eta_ij = 2 * (0.7*(Real)Dimensions+2.1) * (vel_[index_i] - vel_[index_j]).dot(e_ij) / (r_ij + TinyReal);
+                acceleration += eta_ij * dW_ijV_j * e_ij;
+            }
+            acc_shear_[index_i] += G_ * acceleration * dt / rho_i;
+        }
+        //=================================================================================================//
+        ShearStressRelaxation ::
+            ShearStressRelaxation(BaseInnerRelation& inner_relation)
+            : BaseRelaxation(inner_relation),
+            shear_stress_(particles_->shear_stress_), shear_stress_rate_(particles_->shear_stress_rate_),
+            velocity_gradient_(particles_->velocity_gradient_), strain_tensor_(particles_->strain_tensor_),
+            strain_tensor_rate_(particles_->strain_tensor_rate_), von_mises_stress_(particles_->von_mises_stress_),
+            von_mises_strain_(particles_->von_mises_strain_), Vol_(particles_->Vol_), 
+            B_(*particles_->getVariableByName<Matd>("CorrectionMatrix")) {}
+        void ShearStressRelaxation::initialization(size_t index_i, Real dt)
+        {
+            strain_tensor_[index_i] += strain_tensor_rate_[index_i] * 0.5 * dt;
+            shear_stress_[index_i] += shear_stress_rate_[index_i] * 0.5 * dt;
+        }
+        void ShearStressRelaxation::interaction(size_t index_i, Real dt)
+        {
+            Matd velocity_gradient = Matd::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_ijV_j_[n];
+                Vecd& e_ij = inner_neighborhood.e_ij_[n];
+
+                //Matd velocity_gradient_ij = - (vel_[index_i] - vel_[index_j]) * nablaW_ijV_j.transpose();
+                Vecd v_ij = vel_[index_i] - vel_[index_j];
+                Matd velocity_gradient_ij = -v_ij * (B_[index_i] * e_ij * dW_ijV_j_).transpose();
+                velocity_gradient += velocity_gradient_ij;
+            }
+            velocity_gradient_[index_i] = velocity_gradient;
+            //calculate strain
+            Matd strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose());
+            strain_tensor_rate_[index_i] = strain_rate;
+            strain_tensor_[index_i] += strain_tensor_rate_[index_i] * 0.5 * dt;
+            Matd strain_i = strain_tensor_[index_i];
+            von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_i);
+        }
+        void ShearStressRelaxation::update(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 * 0.5;
+
+            //VonMises Stress
+            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())
+        {
+            for (int k = 0; k != Dimensions; ++k)
+                constrain_matrix_(k, k) = constrained_axises[k];
+        };
+        //=================================================================================================//
+        void FixedInAxisDirection::update(size_t index_i, Real dt)
+        {
+            vel_[index_i] = constrain_matrix_ * vel_[index_i];
+        };
+
+        //=================================================================================================//
+        ConstrainSolidBodyMassCenter::
+            ConstrainSolidBodyMassCenter(SPHBody& sph_body, Vecd constrain_direction)
+            : LocalDynamics(sph_body), ContinuumDataSimple(sph_body),
+            correction_matrix_(Matd::Identity()), vel_(particles_->vel_),
+            compute_total_momentum_(sph_body, "Velocity")
+        {
+            for (int i = 0; i != Dimensions; ++i)
+                correction_matrix_(i, i) = constrain_direction[i];
+            ReduceDynamics<QuantitySummation<Real>> compute_total_mass_(sph_body, "MassiveMeasure");
+            total_mass_ = compute_total_mass_.exec();
+        }
+        //=================================================================================================//
+        void ConstrainSolidBodyMassCenter::setupDynamics(Real dt)
+        {
+            velocity_correction_ =
+                correction_matrix_ * compute_total_momentum_.exec(dt) / total_mass_;
+        }
+        //=================================================================================================//
+        void ConstrainSolidBodyMassCenter::update(size_t index_i, Real dt)
+        {
+            vel_[index_i] -= velocity_correction_;
+        }
+    } // namespace continuum_dynamics
 } // namespace SPH

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

@@ -7,7 +7,7 @@
 
 namespace SPH
 {
-namespace continuum_dyannmics
+namespace continuum_dynamics
 {
 typedef DataDelegateSimple<ContinuumParticles> ContinuumDataSimple;
 typedef DataDelegateInner<ContinuumParticles> ContinuumDataInner;
@@ -42,7 +42,7 @@ class BaseRelaxation : public LocalDynamics, public ContinuumDataInner
     virtual ~BaseRelaxation(){};
 
   protected:
-    GeneralContinuum &granular_material_;
+    GeneralContinuum &continuum_;
     StdLargeVec<Real> &rho_, &p_, &drho_dt_;
     StdLargeVec<Vecd> &pos_, &vel_, &acc_, &acc_prior_;
 };
@@ -110,6 +110,77 @@ class ShearStressRelaxation2ndHalf : public BaseRelaxation
     StdLargeVec<Real> &von_mises_stress_;
 };
 
+//=================================================================================================//
+//===================================Non-hourglass formulation=====================================//
+//=================================================================================================//
+/**
+    * @class Integration1stHalf
+    * @brief Pressure relaxation scheme with the mostly used Riemann solver.
+*/
+template <class FluidDynamicsType>
+class BaseIntegration1stHalf : public FluidDynamicsType
+{
+public:
+    explicit BaseIntegration1stHalf(BaseInnerRelation& inner_relation);
+    virtual ~BaseIntegration1stHalf() {};
+    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:
+    StdLargeVec<Vecd>& acc_shear_;
+};
+using Integration1stHalf = BaseIntegration1stHalf<fluid_dynamics::Integration1stHalf>;
+using Integration1stHalfRiemann = BaseIntegration1stHalf<fluid_dynamics::Integration1stHalfRiemann>;
+
+/**
+* @class ShearAccelerationRelaxation
+*/
+class ShearAccelerationRelaxation : public BaseRelaxation
+{
+public:
+    explicit ShearAccelerationRelaxation(BaseInnerRelation& inner_relation);
+    virtual ~ShearAccelerationRelaxation() {};
+    void interaction(size_t index_i, Real dt = 0.0);
+protected:
+    Real G_, smoothing_length_;
+    StdLargeVec<Matd>& shear_stress_, & B_;
+    StdLargeVec<Vecd>& acc_shear_;
+};
+
+/**
+ * @class AngularConservativeShearAccelerationRelaxation
+ */
+class AngularConservativeShearAccelerationRelaxation : public ShearAccelerationRelaxation
+{
+public:
+    explicit AngularConservativeShearAccelerationRelaxation(BaseInnerRelation& inner_relation)
+        : ShearAccelerationRelaxation(inner_relation) {};
+    virtual ~AngularConservativeShearAccelerationRelaxation() {};
+
+    void interaction(size_t index_i, Real dt = 0.0);
+    //void update(size_t index_i, Real dt = 0.0);
+
+};
+
+/**
+* @class ShearStressRelaxation
+*/
+class ShearStressRelaxation : public BaseRelaxation
+{
+public:
+
+    explicit ShearStressRelaxation(BaseInnerRelation& inner_relation);
+    virtual ~ShearStressRelaxation() {};
+    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:
+    StdLargeVec<Matd>& shear_stress_, & shear_stress_rate_, & velocity_gradient_, & strain_tensor_, & strain_tensor_rate_;
+    StdLargeVec<Real>& von_mises_stress_, & von_mises_strain_, & Vol_;
+    StdLargeVec<Matd>& B_;
+};
+
 /**
  * @class BaseMotionConstraint
  */
@@ -145,6 +216,44 @@ class FixConstraint : public BaseMotionConstraint<DynamicsIdentifier>
 };
 using FixBodyConstraint = FixConstraint<SPHBody>;
 using FixBodyPartConstraint = FixConstraint<BodyPartByParticle>;
-} // namespace continuum_dyannmics
+
+/**
+ * @class FixedInAxisDirection
+ * @brief Constrain the velocity of a solid body part.
+ */
+class FixedInAxisDirection : public BaseMotionConstraint<BodyPartByParticle>
+{
+public:
+    FixedInAxisDirection(BodyPartByParticle& body_part, Vecd constrained_axises = Vecd::Zero());
+    virtual ~FixedInAxisDirection() {};
+    void update(size_t index_i, Real dt = 0.0);
+
+protected:
+    Matd constrain_matrix_;
+};
+
+/**
+ * @class ConstrainSolidBodyMassCenter
+ * @brief Constrain the mass center of a solid body.
+ */
+class ConstrainSolidBodyMassCenter : public LocalDynamics, public ContinuumDataSimple
+{
+private:
+    Real total_mass_;
+    Matd correction_matrix_;
+    Vecd velocity_correction_;
+    StdLargeVec<Vecd>& vel_;
+    ReduceDynamics<QuantityMoment<Vecd>> compute_total_momentum_;
+
+protected:
+    virtual void setupDynamics(Real dt = 0.0) override;
+
+public:
+    explicit ConstrainSolidBodyMassCenter(SPHBody& sph_body, Vecd constrain_direction = Vecd::Ones());
+    virtual ~ConstrainSolidBodyMassCenter() {};
+
+    void update(size_t index_i, Real dt = 0.0);
+};
+} // namespace continuum_dynamics
 } // namespace SPH
 #endif // CONTINUUM_DYNAMICS_INNER_H