From 8dc7f2e657d43c6bed227116b98255e0bedec799 Mon Sep 17 00:00:00 2001 From: WeiyiVirtonomy Date: Fri, 9 Aug 2024 17:04:43 +0200 Subject: [PATCH 1/2] add ContactSiffness(index_i) --- src/shared/materials/base_material.h | 1 + src/shared/materials/complex_solid.cpp | 5 +++++ src/shared/materials/complex_solid.h | 2 ++ 3 files changed, 8 insertions(+) diff --git a/src/shared/materials/base_material.h b/src/shared/materials/base_material.h index f170b0c93b..7ba8784439 100644 --- a/src/shared/materials/base_material.h +++ b/src/shared/materials/base_material.h @@ -113,6 +113,7 @@ class Solid : public BaseMaterial Real ContactFriction() { return contact_friction_; }; Real ContactStiffness() { return contact_stiffness_; }; + virtual Real ContactStiffness(size_t index_i) { return contact_stiffness_; }; virtual Solid *ThisObjectPtr() override { return this; }; /** Get average velocity when interacting with fluid. */ virtual StdLargeVec *AverageVelocity(BaseParticles *base_particles); diff --git a/src/shared/materials/complex_solid.cpp b/src/shared/materials/complex_solid.cpp index f284912598..22e6acc7ec 100644 --- a/src/shared/materials/complex_solid.cpp +++ b/src/shared/materials/complex_solid.cpp @@ -36,6 +36,11 @@ Real CompositeSolid::CompositeDensity(size_t index_i) return composite_materials_[(*material_id_)[index_i]]->ReferenceDensity(); } //=================================================================================================// +Real CompositeSolid::ContactStiffness(size_t index_i) +{ + return composite_materials_[(*material_id_)[index_i]]->ContactStiffness(); +} +//=================================================================================================// void CompositeSolid::initializeLocalParameters(BaseParticles *base_particles) { ElasticSolid::initializeLocalParameters(base_particles); diff --git a/src/shared/materials/complex_solid.h b/src/shared/materials/complex_solid.h index 691a8678a5..029ba4037c 100644 --- a/src/shared/materials/complex_solid.h +++ b/src/shared/materials/complex_solid.h @@ -78,6 +78,8 @@ class CompositeSolid : public ElasticSolid virtual Real VolumetricKirchhoff(Real J) override { return 0.0; }; virtual std::string getRelevantStressMeasureName() override { return "PK2"; }; + Real ContactStiffness(size_t index_i) override; + Real CompositeDensity(size_t index_i); template From d30d16f13062095acf96f94a1987e7517021c9fd Mon Sep 17 00:00:00 2001 From: WeiyiVirtonomy Date: Fri, 9 Aug 2024 17:05:42 +0200 Subject: [PATCH 2/2] use ContactStiffness(index_i) in contact repulsion --- .../contact_dynamics/contact_repulsion.cpp | 22 +++++++++---------- .../contact_dynamics/contact_repulsion.h | 1 - 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.cpp b/src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.cpp index 20d9c27e5e..deb4fc330f 100644 --- a/src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.cpp +++ b/src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.cpp @@ -37,16 +37,11 @@ RepulsionForce>::RepulsionForce(BaseContactRelation &solid_body_contac solid_(DynamicCast(this, sph_body_.getBaseMaterial())), repulsion_factor_(*particles_->getVariableDataByName("RepulsionFactor")) { - const Real contact_stiffness_1 = solid_.ContactStiffness(); - for (size_t k = 0; k != contact_particles_.size(); ++k) { contact_solids_.push_back(&DynamicCast(this, contact_bodies_[k]->getBaseMaterial())); contact_Vol_.push_back(contact_particles_[k]->getVariableDataByName("VolumetricMeasure")); contact_repulsion_factor_.push_back(contact_particles_[k]->getVariableDataByName("RepulsionFactor")); - - const Real contact_stiffness_k = contact_solids_[k]->ContactStiffness(); - contact_stiffness_ave_.push_back(2 * contact_stiffness_1 * contact_stiffness_k / (contact_stiffness_1 + contact_stiffness_k)); } } //=================================================================================================// @@ -54,10 +49,10 @@ void RepulsionForce>::interaction(size_t index_i, Real dt) { Real sigma_i = repulsion_factor_[index_i]; Vecd force = Vecd::Zero(); + Real contact_stiffness_i = solid_.ContactStiffness(index_i); + for (size_t k = 0; k < contact_configuration_.size(); ++k) { - Vecd force_k = Vecd::Zero(); - StdLargeVec &contact_density_k = *(contact_repulsion_factor_[k]); StdLargeVec &Vol_k = *(contact_Vol_[k]); @@ -65,13 +60,16 @@ void RepulsionForce>::interaction(size_t index_i, Real dt) for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) { size_t index_j = contact_neighborhood.j_[n]; + + Real contact_stiffness_j = contact_solids_[k]->ContactStiffness(index_j); + Real contact_stiffness_ave = 2 * contact_stiffness_i * contact_stiffness_j / (contact_stiffness_i + contact_stiffness_j); + Vecd e_ij = contact_neighborhood.e_ij_[n]; - Real sigma_star = 0.5 * (sigma_i + contact_density_k[index_j]); + Real sigma_star = (sigma_i + contact_density_k[index_j]) * contact_stiffness_ave; // force due to pressure - force_k -= 2.0 * sigma_star * e_ij * contact_neighborhood.dW_ij_[n] * Vol_k[index_j]; + force -= sigma_star * e_ij * contact_neighborhood.dW_ij_[n] * Vol_k[index_j]; } - force += force_k * contact_stiffness_ave_[k]; } repulsion_force_[index_i] = force * particles_->ParticleVolume(index_i); } @@ -90,7 +88,7 @@ RepulsionForce>::RepulsionForce(BaseContactRelation &solid_body_co //=================================================================================================// void RepulsionForce>::interaction(size_t index_i, Real dt) { - Real p_i = repulsion_factor_[index_i] * solid_.ContactStiffness(); + Real p_i = repulsion_factor_[index_i] * solid_.ContactStiffness(index_i); /** Contact interaction. */ Vecd force = Vecd::Zero(); for (size_t k = 0; k < contact_configuration_.size(); ++k) @@ -136,7 +134,7 @@ void RepulsionForce>::interaction(size_t index_i, Real dt) size_t index_j = contact_neighborhood.j_[n]; Vecd e_ij = contact_neighborhood.e_ij_[n]; - Real p_star = contact_density_k[index_j] * solid_k->ContactStiffness(); + Real p_star = contact_density_k[index_j] * solid_k->ContactStiffness(index_j); // force due to pressure force -= 2.0 * p_star * e_ij * contact_neighborhood.dW_ij_[n] * Vol_k[index_j]; } diff --git a/src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.h b/src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.h index 762e48da58..3aa010e99e 100644 --- a/src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.h +++ b/src/shared/particle_dynamics/solid_dynamics/contact_dynamics/contact_repulsion.h @@ -86,7 +86,6 @@ class RepulsionForce> : public RepulsionForce &repulsion_factor_; StdVec contact_solids_; - StdVec contact_stiffness_ave_; StdVec *> contact_repulsion_factor_, contact_Vol_; }; using ContactForce = RepulsionForce>;