contact density non-dimensionlized

src/shared/materials/base_material.h

@@ -103,27 +103,30 @@ class Solid : public BaseMaterial
   public:
     Solid(Real rho0, Real contact_stiffness, Real contact_friction = 0.0)
         : BaseMaterial(rho0), contact_stiffness_(contact_stiffness),
+          contact_impendence_(rho0 * sqrt(contact_stiffness / rho0)),
           contact_friction_(contact_friction)
     {
         material_type_name_ = "Solid";
     };
-    explicit Solid(Real rho0) : Solid(rho0, 1.0){};
+    explicit Solid(Real rho0) : Solid(rho0, rho0){};
     Solid() : Solid(1.0){};
     virtual ~Solid(){};
 
     Real ContactFriction() { return contact_friction_; };
     Real ContactStiffness() { return contact_stiffness_; };
+    Real ContactImpedence() { return contact_impendence_; };
     virtual Solid *ThisObjectPtr() override { return this; };
     /** Get average velocity when interacting with fluid. */
     virtual StdLargeVec<Vecd> *AverageVelocity(BaseParticles *base_particles);
     /** Get average acceleration when interacting with fluid. */
     virtual StdLargeVec<Vecd> *AverageAcceleration(BaseParticles *base_particles);
 
   protected:
-    Real contact_stiffness_; /**< contact-force stiffness related to bulk modulus*/
-    Real contact_friction_;  /**< friction property mimic fluid viscosity*/
-
-    void setContactStiffness(Real c0) { contact_stiffness_ = c0 * c0; };
+    Real contact_stiffness_;  /**< Stiffness has the unit of bulk modules*/
+    Real contact_impendence_; /**< contact-impendence */
+    Real contact_friction_;   /**< friction property mimic fluid viscosity */
+    void setContactStiffness(Real c0) { contact_stiffness_ = rho0_ * c0 * c0; };
+    void setContactImpedence(Real c0) { contact_impendence_ = rho0_ * c0; };
 };
 } // namespace SPH
 #endif // BASE_MATERIAL_H

src/shared/materials/complex_solid.h

@@ -87,7 +87,9 @@ class CompositeSolid : public ElasticSolid
             composite_ptrs_keeper_.createPtr<ElasticSolidType>(std::forward<Args>(args)...);
         composite_materials_.push_back(added_material);
         c0_ = SMAX(c0_, added_material->ReferenceSoundSpeed());
+        rho0_ = SMAX(rho0_, added_material->ReferenceDensity());
         setContactStiffness(c0_);
+        setContactImpedence(c0_);
     };
 };
 

src/shared/materials/elastic_solid.cpp

@@ -46,6 +46,7 @@ LinearElasticSolid::
     lambda0_ = getLambda(youngs_modulus, poisson_ratio);
     setSoundSpeeds();
     setContactStiffness(c0_);
+    setContactImpedence(c0_);
 }
 //=================================================================================================//
 Real LinearElasticSolid::getBulkModulus(Real youngs_modulus, Real poisson_ratio)

src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.cpp

@@ -9,10 +9,9 @@ RepulsionForce<Contact<Inner<>>>::
     RepulsionForce(SelfSurfaceContactRelation &self_contact_relation)
     : RepulsionForce<Base, DataDelegateInner>(self_contact_relation, "SelfRepulsionForce"),
       ForcePrior(particles_, "SelfRepulsionForce"),
-      solid_(DynamicCast<Solid>(this, sph_body_.getBaseMaterial())),
       self_repulsion_density_(*particles_->getVariableDataByName<Real>("SelfRepulsionDensity")),
       vel_(*particles_->getVariableDataByName<Vecd>("Velocity")),
-      contact_impedance_(solid_.ReferenceDensity() * sqrt(solid_.ContactStiffness())) {}
+      contact_impedance_(solid_.ContactImpedence()) {}
 //=================================================================================================//
 void RepulsionForce<Contact<Inner<>>>::interaction(size_t index_i, Real dt)
 {
@@ -34,7 +33,6 @@ void RepulsionForce<Contact<Inner<>>>::interaction(size_t index_i, Real dt)
 RepulsionForce<Contact<>>::RepulsionForce(SurfaceContactRelation &solid_body_contact_relation)
     : RepulsionForce<Base, DataDelegateContact>(solid_body_contact_relation, "RepulsionForce"),
       ForcePrior(particles_, "RepulsionForce"),
-      solid_(DynamicCast<Solid>(this, sph_body_.getBaseMaterial())),
       repulsion_density_(*particles_->getVariableDataByName<Real>("RepulsionDensity"))
 {
     for (size_t k = 0; k != contact_particles_.size(); ++k)
@@ -72,7 +70,6 @@ void RepulsionForce<Contact<>>::interaction(size_t index_i, Real dt)
 RepulsionForce<Contact<Wall>>::RepulsionForce(SurfaceContactRelation &solid_body_contact_relation)
     : RepulsionForce<Base, DataDelegateContact>(solid_body_contact_relation, "RepulsionForce"),
       ForcePrior(particles_, "RepulsionForce"),
-      solid_(DynamicCast<Solid>(this, sph_body_.getBaseMaterial())),
       repulsion_density_(*particles_->getVariableDataByName<Real>("RepulsionDensity"))
 {
     for (size_t k = 0; k < this->contact_configuration_.size(); ++k)

src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.h

@@ -49,11 +49,13 @@ class RepulsionForce<Base, DataDelegationType>
     template <class BaseRelationType>
     RepulsionForce(BaseRelationType &base_relation, const std::string &variable_name)
         : LocalDynamics(base_relation.getSPHBody()), DataDelegationType(base_relation),
+          solid_(DynamicCast<Solid>(this, this->sph_body_.getBaseMaterial())),
           repulsion_force_(*this->particles_->template registerSharedVariable<Vecd>(variable_name)),
           Vol_(*this->particles_->template getVariableDataByName<Real>("VolumetricMeasure")){};
     virtual ~RepulsionForce(){};
 
   protected:
+    Solid &solid_;
     StdLargeVec<Vecd> &repulsion_force_;
     StdLargeVec<Real> &Vol_;
 };
@@ -67,7 +69,6 @@ class RepulsionForce<Contact<Inner<>>> : public RepulsionForce<Base, DataDelegat
     void interaction(size_t index_i, Real dt = 0.0);
 
   protected:
-    Solid &solid_;
     StdLargeVec<Real> &self_repulsion_density_;
     StdLargeVec<Vecd> &vel_;
     Real contact_impedance_;
@@ -83,7 +84,6 @@ class RepulsionForce<Contact<>> : public RepulsionForce<Base, DataDelegateContac
     void interaction(size_t index_i, Real dt = 0.0);
 
   protected:
-    Solid &solid_;
     StdLargeVec<Real> &repulsion_density_;
     StdVec<Solid *> contact_solids_;
     StdVec<StdLargeVec<Real> *> contact_contact_density_, contact_Vol_;
@@ -99,7 +99,6 @@ class RepulsionForce<Contact<Wall>> : public RepulsionForce<Base, DataDelegateCo
     void interaction(size_t index_i, Real dt = 0.0);
 
   protected:
-    Solid &solid_;
     StdLargeVec<Real> &repulsion_density_;
     StdVec<StdLargeVec<Real> *> contact_Vol_;
 };

src/shared/particle_dynamics/solid_dynamics/contact_dynamics/repulsion_density_summation.cpp

@@ -8,7 +8,7 @@ namespace solid_dynamics
 RepulsionDensitySummation<Inner<>>::
     RepulsionDensitySummation(SelfSurfaceContactRelation &self_contact_relation)
     : RepulsionDensitySummation<Base, DataDelegateInner>(self_contact_relation, "SelfRepulsionDensity"),
-      mass_(*particles_->getVariableDataByName<Real>("Mass"))
+      Vol_(*particles_->getVariableDataByName<Real>("VolumetricMeasure"))
 {
     Real dp_1 = self_contact_relation.getSPHBody().sph_adaptation_->ReferenceSpacing();
     offset_W_ij_ = self_contact_relation.getSPHBody().sph_adaptation_->getKernel()->W(dp_1, ZeroVecd);
@@ -20,15 +20,14 @@ void RepulsionDensitySummation<Inner<>>::interaction(size_t index_i, Real dt)
     const Neighborhood &inner_neighborhood = inner_configuration_[index_i];
     for (size_t n = 0; n != inner_neighborhood.current_size_; ++n)
     {
-        sigma += std::max(inner_neighborhood.W_ij_[n] - offset_W_ij_, Real(0));
+        sigma += std::max(inner_neighborhood.W_ij_[n] - offset_W_ij_, Real(0)) * Vol_[inner_neighborhood.j_[n]];
     }
-    repulsion_density_[index_i] = mass_[index_i] * sigma;
+    repulsion_density_[index_i] = sigma;
 }
 //=================================================================================================//
 RepulsionDensitySummation<Contact<>>::
     RepulsionDensitySummation(SurfaceContactRelation &solid_body_contact_relation)
     : RepulsionDensitySummation<Base, DataDelegateContact>(solid_body_contact_relation, "RepulsionDensity"),
-      rho0_(sph_body_.getBaseMaterial().ReferenceDensity()),
       offset_W_ij_(StdVec<Real>(contact_configuration_.size(), 0.0))
 {
     for (size_t k = 0; k != contact_particles_.size(); ++k)
@@ -66,7 +65,7 @@ void RepulsionDensitySummation<Contact<>>::interaction(size_t index_i, Real dt)
             sigma += corrected_W_ij * contact_Vol0_k[contact_neighborhood.j_[n]];
         }
     }
-    repulsion_density_[index_i] = rho0_ * sigma;
+    repulsion_density_[index_i] = sigma;
 };
 //=================================================================================================//
 ShellContactDensity::ShellContactDensity(SurfaceContactRelation &solid_body_contact_relation)
@@ -91,7 +90,7 @@ ShellContactDensity::ShellContactDensity(SurfaceContactRelation &solid_body_cont
             contact_max += Dimensions == 2 ? contact_temp : contact_temp * Pi * temp;
         }
         /** a calibration factor to avoid particle penetration into shell structure */
-        calibration_factor_.push_back(solid_.ReferenceDensity() / (contact_max + Eps));
+        calibration_factor_.push_back(1.0 / (contact_max + Eps));
 
         contact_Vol_.push_back(contact_particles_[k]->getVariableDataByName<Real>("VolumetricMeasure"));
     }

src/shared/particle_dynamics/solid_dynamics/contact_dynamics/repulsion_density_summation.h

@@ -64,7 +64,7 @@ class RepulsionDensitySummation<Inner<>> : public RepulsionDensitySummation<Base
     void interaction(size_t index_i, Real dt = 0.0);
 
   protected:
-    StdLargeVec<Real> &mass_;
+    StdLargeVec<Real> &Vol_;
     Real offset_W_ij_;
 };
 using SelfContactDensitySummation = RepulsionDensitySummation<Inner<>>;
@@ -78,7 +78,6 @@ class RepulsionDensitySummation<Contact<>> : public RepulsionDensitySummation<Ba
     void interaction(size_t index_i, Real dt = 0.0);
 
   protected:
-    Real rho0_;
     StdVec<StdLargeVec<Real> *> contact_Vol0_;
     StdVec<Real> offset_W_ij_;
 };