diff --git a/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp b/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
index fb1e42957ef2de7b349259947d59897fa9d7efee..18adc77137d0591fedcec9e297af97e6b41a18d2 100644
--- a/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
@@ -111,20 +111,7 @@ int main( int argc, char ** argv )
const real_t particleMass = real_t(1) / ss->shapes[sphereShape]->getInvMass();
const real_t Mij = particleMass; // * particleMass / ( real_t(2) * particleMass ); // Mij = M for sphere-wall collision
- const real_t lnDryResCoeff = std::log(restitutionCoeff);
-
- // normal material parameters
- const real_t stiffnessN = Mij * ( math::pi * math::pi + lnDryResCoeff * lnDryResCoeff ) / (collisionTime * collisionTime); // Costa et al., Eq. 18
- const real_t dampingN = - real_t(2) * Mij * lnDryResCoeff / collisionTime; // Costa et al., Eq. 18
-
- WALBERLA_LOG_INFO_ON_ROOT("normal: stiffness = " << stiffnessN << ", damping = " << dampingN);
-
- // from Thornton et al
- const real_t kappa = real_t(2) * ( real_t(1) - nu ) / ( real_t(2) - nu ) ;
- const real_t stiffnessT = kappa * stiffnessN;
- const real_t dampingT = std::sqrt(kappa) * dampingN;
-
- WALBERLA_LOG_INFO_ON_ROOT("tangential: kappa = " << kappa << ", stiffness T = " << stiffnessT << ", damping T = " << dampingT);
+ const real_t kappa = real_t(2) * ( real_t(1) - nu ) / ( real_t(2) - nu ) ; // from Thornton et al
real_t uTin = uNin * impactAngle;
@@ -154,10 +141,7 @@ int main( int argc, char ** argv )
kernel::DoubleCast double_cast;
kernel::LinearSpringDashpot dem(1);
mpi::ReduceContactHistory rch;
- dem.setStiffnessN(0,0,stiffnessN);
- dem.setStiffnessT(0,0,stiffnessT);
- dem.setDampingN(0,0,dampingN);
- dem.setDampingT(0,0,dampingT);
+ dem.setStiffnessAndDamping(0,0,restitutionCoeff,collisionTime,kappa,Mij);
dem.setFrictionCoefficientStatic(0,0,frictionCoeff_s);
dem.setFrictionCoefficientDynamic(0,0,frictionCoeff_d);
diff --git a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
index 966a55925f85fea9a051687ae6836afbfc27ca7b..12ef798b27eb670f4b32135f656be9217999a197 100644
--- a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
@@ -694,15 +694,7 @@ int main( int argc, char **argv )
const real_t particleMass = densitySphere * sphereVolume;
const real_t Mij = particleMass; // * particleMass / ( real_t(2) * particleMass ); // Mij = M for sphere-wall collision
- const real_t lnDryResCoeff = std::log(restitutionCoeff);
-
- const real_t stiffnessN = Mij * ( math::pi * math::pi + lnDryResCoeff * lnDryResCoeff ) / (collisionTime * collisionTime); // Costa et al., Eq. 18
- const real_t dampingN = - real_t(2) * Mij * lnDryResCoeff / collisionTime; // Costa et al., Eq. 18
-
- // from Biegert et al
const real_t kappa = real_t(2) * ( real_t(1) - poissonsRatio ) / ( real_t(2) - poissonsRatio ) ;
- const real_t stiffnessT = kappa * stiffnessN;
- const real_t dampingT = std::sqrt(kappa) * dampingN;
Vector3<uint_t> domainSize( uint_c(domainSizeNonDim[0] * diameter ), uint_c(domainSizeNonDim[1] * diameter ), uint_c(domainSizeNonDim[2] * diameter ) );
@@ -735,10 +727,6 @@ int main( int argc, char **argv )
WALBERLA_LOG_INFO_ON_ROOT("Collision Response properties:" );
WALBERLA_LOG_INFO_ON_ROOT(" - collision time = " << collisionTime );
- WALBERLA_LOG_INFO_ON_ROOT(" - damping coeff n = " << dampingN );
- WALBERLA_LOG_INFO_ON_ROOT(" - stiffness coeff n = " << stiffnessN );
- WALBERLA_LOG_INFO_ON_ROOT(" - damping coeff t = " << dampingT );
- WALBERLA_LOG_INFO_ON_ROOT(" - stiffness coeff t= " << stiffnessT );
WALBERLA_LOG_INFO_ON_ROOT(" - coeff of restitution = " << restitutionCoeff );
WALBERLA_LOG_INFO_ON_ROOT(" - coeff of friction = " << frictionCoeff );
@@ -875,10 +863,7 @@ int main( int argc, char **argv )
// linear model
mesa_pd::kernel::LinearSpringDashpot linearCollisionResponse(1);
- linearCollisionResponse.setStiffnessN(0,0,stiffnessN);
- linearCollisionResponse.setDampingN(0,0,dampingN);
- linearCollisionResponse.setStiffnessT(0,0,stiffnessT);
- linearCollisionResponse.setDampingT(0,0,dampingT);
+ linearCollisionResponse.setStiffnessAndDamping(0,0,restitutionCoeff,collisionTime,kappa,Mij);
//linearCollisionResponse.setFrictionCoefficientStatic(0,0,frictionCoeff); // not used in this test case
linearCollisionResponse.setFrictionCoefficientDynamic(0,0,frictionCoeff);
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
index 8a339e6915f29dff6fcaea552b68923deb5540ce..0890c4ef497d9b840374a0cb4fe75bea2d4ff7c7 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
@@ -666,16 +666,6 @@ int main( int argc, char **argv )
const real_t particleMass = densityRatio * sphereVolume;
const real_t Mij = particleMass; // * particleMass / ( real_t(2) * particleMass ); // Mij = M for sphere-wall collision
- const real_t lnDryResCoeff = std::log(restitutionCoeff);
- const real_t stiffnessCoeff = Mij * ( math::pi * math::pi + lnDryResCoeff * lnDryResCoeff ) / (collisionTime * collisionTime); // Costa et al., Eq. 18
- const real_t dampingCoeff = - real_t(2) * Mij * lnDryResCoeff / collisionTime; // Costa et al., Eq. 18
- //const real_t contactDuration = real_t(2) * math::pi * Mij / ( std::sqrt( real_t(4) * Mij * stiffnessCoeff - dampingCoeff * dampingCoeff )); //formula from Uhlman
-
- const real_t var_a = std::sqrt( math::pi * math::pi + (std::log(restitutionCoeff) * std::log(restitutionCoeff) ) );
- const real_t TnLimit1 = var_a * diameter / uIn * std::exp(-std::asin(math::pi / var_a ) / math::pi); // Costa et al, Eq. 20
- const real_t TnLimit2 = std::sqrt(diameter / gravitationalAcceleration * var_a * var_a / std::abs(real_t(1) - densityFluid/densitySphere)); // Costa et al, Eq. 21
-
- const real_t overlapDueToWeight = std::abs(densitySphere - densityFluid) * gravitationalAcceleration * sphereVolume / stiffnessCoeff; // Costa et al., given in text after Eq. 34
const real_t uInCrit = real_t(9) * StCrit * viscosity / ( densityRatio * diameter);
@@ -741,12 +731,8 @@ int main( int argc, char **argv )
}
else
{
- WALBERLA_LOG_INFO_ON_ROOT(" - using linear collision model with fixed parameters:" );
- WALBERLA_LOG_INFO_ON_ROOT(" - damping coeff = " << dampingCoeff );
- WALBERLA_LOG_INFO_ON_ROOT(" - stiffness coeff = " << stiffnessCoeff );
- WALBERLA_LOG_INFO_ON_ROOT(" - overlap due to particle weight = " << overlapDueToWeight);
+ WALBERLA_LOG_INFO_ON_ROOT(" - using linear collision model with fixed parameters" );
}
- WALBERLA_LOG_INFO_ON_ROOT(" - TnCrit1 = " << TnLimit1 << ", TnCrit2 = " << TnLimit2 );
WALBERLA_LOG_INFO_ON_ROOT(" - coeff of restitution = " << restitutionCoeff );
WALBERLA_LOG_INFO_ON_ROOT(" - number of RPD sub cycles = " << numRPDSubCycles );
WALBERLA_LOG_INFO_ON_ROOT(" - lubrication correction = " << (useLubricationCorrection ? "yes" : "no") );
@@ -963,12 +949,11 @@ int main( int argc, char **argv )
// linear model
mesa_pd::kernel::LinearSpringDashpot linearCollisionResponse(1);
- linearCollisionResponse.setStiffnessN(0,0,stiffnessCoeff);
- linearCollisionResponse.setDampingN(0,0,dampingCoeff);
+ linearCollisionResponse.setStiffnessAndDamping(0,0,restitutionCoeff,collisionTime,real_t(0),Mij); // no response in tangential direction
// nonlinear model for ACTM
mesa_pd::kernel::NonLinearSpringDashpot nonLinearCollisionResponse(1, collisionTime);
- nonLinearCollisionResponse.setLnCORsqr(0,0,lnDryResCoeff * lnDryResCoeff);
+ nonLinearCollisionResponse.setLnCORsqr(0,0,std::log(restitutionCoeff) * std::log(restitutionCoeff));
nonLinearCollisionResponse.setMeff(0,0,Mij);
mesa_pd::mpi::ReduceProperty reduceProperty;
diff --git a/python/mesa_pd/templates/kernel/LinearSpringDashpot.templ.h b/python/mesa_pd/templates/kernel/LinearSpringDashpot.templ.h
index f5c1b8e224bdf5e61eb3bd23e545211e63f12f59..2e86c1c1d150ea53a4dc5344c3823a2e6a905aa0 100644
--- a/python/mesa_pd/templates/kernel/LinearSpringDashpot.templ.h
+++ b/python/mesa_pd/templates/kernel/LinearSpringDashpot.templ.h
@@ -81,6 +81,21 @@ public:
const real_t& penetrationDepth,
const real_t& dt) const;
+ real_t calcCoefficientOfRestitution(const size_t type1,
+ const size_t type2,
+ const real_t effectiveMass);
+
+ real_t calcCollisionTime(const size_t type1,
+ const size_t type2,
+ const real_t effectiveMass);
+
+ void setStiffnessAndDamping(const size_t type1,
+ const size_t type2,
+ const real_t coefficientOfRestitution,
+ const real_t collisionTime,
+ const real_t poissonsRatio,
+ const real_t effectiveMass);
+
{% for param in parameters %}
/// assumes this parameter is symmetric
void set{{param | capFirst}}(const size_t type1, const size_t type2, const real_t& val);
@@ -174,7 +189,7 @@ inline void LinearSpringDashpot::operator()(const size_t p_idx1,
impactVelocityMagnitude = relVel.length();
}
- //TODO: move to own tangential integration kernel?
+ //TODO: move to own tangential displacement integration kernel
Vec3 rotatedTangentialDisplacement = tangentialSpringDisplacement - contactNormal * (contactNormal * tangentialSpringDisplacement);
Vec3 newTangentialSpringDisplacement = rotatedTangentialDisplacement.sqrLength() <= real_t(0) ? // avoid division by zero
Vec3(real_t(0)) :
@@ -218,7 +233,6 @@ inline void LinearSpringDashpot::operator()(const size_t p_idx1,
// if tangential force falls below coulomb limit, we are back in sticking
if( fTLS.length() < fFrictionAbsDynamic )
{
- //TODO really?
isSticking = true;
}
}
@@ -248,6 +262,39 @@ inline void LinearSpringDashpot::operator()(const size_t p_idx1,
}
}
+inline real_t LinearSpringDashpot::calcCoefficientOfRestitution(const size_t type1,
+ const size_t type2,
+ const real_t effectiveMass)
+{
+ auto a = real_t(0.5) * getDampingN(type1, type2) / effectiveMass;
+ return std::exp(-a * math::pi / std::sqrt(getStiffnessN(type1, type2) / effectiveMass - a*a));
+}
+
+inline real_t LinearSpringDashpot::calcCollisionTime(const size_t type1,
+ const size_t type2,
+ const real_t effectiveMass)
+{
+ auto a = real_t(0.5) * getDampingN(type1, type2) / effectiveMass;
+ return math::pi / std::sqrt( getStiffnessN(type1, type2)/effectiveMass - a*a);
+}
+
+inline void LinearSpringDashpot::setStiffnessAndDamping(const size_t type1,
+ const size_t type2,
+ const real_t coefficientOfRestitution,
+ const real_t collisionTime,
+ const real_t poissonsRatio,
+ const real_t effectiveMass)
+{
+ const real_t kappa = real_t(2) * ( real_t(1) - poissonsRatio ) / ( real_t(2) - poissonsRatio ) ;
+ const real_t lnDryResCoeff = std::log(coefficientOfRestitution);
+
+ setStiffnessN(type1, type2, effectiveMass * ( math::pi * math::pi + lnDryResCoeff * lnDryResCoeff ) / (collisionTime * collisionTime) );
+ setDampingN (type1, type2, - real_t(2) * effectiveMass * lnDryResCoeff / collisionTime );
+
+ setStiffnessT(type1, type2, kappa * getStiffnessN(type1, type2) );
+ setDampingT (type1, type2, std::sqrt(kappa) * getDampingN(type1, type2) );
+}
+
} //namespace kernel
} //namespace mesa_pd
} //namespace walberla
diff --git a/src/mesa_pd/kernel/LinearSpringDashpot.h b/src/mesa_pd/kernel/LinearSpringDashpot.h
index 5239f29b6ef41879ee8ff913ff977fbfcc0f4d12..46db0270b614fed6e80c45f802400fff39370b5c 100644
--- a/src/mesa_pd/kernel/LinearSpringDashpot.h
+++ b/src/mesa_pd/kernel/LinearSpringDashpot.h
@@ -86,6 +86,21 @@ public:
const real_t& penetrationDepth,
const real_t& dt) const;
+ real_t calcCoefficientOfRestitution(const size_t type1,
+ const size_t type2,
+ const real_t effectiveMass);
+
+ real_t calcCollisionTime(const size_t type1,
+ const size_t type2,
+ const real_t effectiveMass);
+
+ void setStiffnessAndDamping(const size_t type1,
+ const size_t type2,
+ const real_t coefficientOfRestitution,
+ const real_t collisionTime,
+ const real_t poissonsRatio,
+ const real_t effectiveMass);
+
/// assumes this parameter is symmetric
void setStiffnessN(const size_t type1, const size_t type2, const real_t& val);
@@ -278,7 +293,7 @@ inline void LinearSpringDashpot::operator()(const size_t p_idx1,
impactVelocityMagnitude = relVel.length();
}
- //TODO: move to own tangential integration kernel?
+ //TODO: move to own tangential displacement integration kernel
Vec3 rotatedTangentialDisplacement = tangentialSpringDisplacement - contactNormal * (contactNormal * tangentialSpringDisplacement);
Vec3 newTangentialSpringDisplacement = rotatedTangentialDisplacement.sqrLength() <= real_t(0) ? // avoid division by zero
Vec3(real_t(0)) :
@@ -322,7 +337,6 @@ inline void LinearSpringDashpot::operator()(const size_t p_idx1,
// if tangential force falls below coulomb limit, we are back in sticking
if( fTLS.length() < fFrictionAbsDynamic )
{
- //TODO really?
isSticking = true;
}
}
@@ -352,6 +366,39 @@ inline void LinearSpringDashpot::operator()(const size_t p_idx1,
}
}
+inline real_t LinearSpringDashpot::calcCoefficientOfRestitution(const size_t type1,
+ const size_t type2,
+ const real_t effectiveMass)
+{
+ auto a = real_t(0.5) * getDampingN(type1, type2) / effectiveMass;
+ return std::exp(-a * math::pi / std::sqrt(getStiffnessN(type1, type2) / effectiveMass - a*a));
+}
+
+inline real_t LinearSpringDashpot::calcCollisionTime(const size_t type1,
+ const size_t type2,
+ const real_t effectiveMass)
+{
+ auto a = real_t(0.5) * getDampingN(type1, type2) / effectiveMass;
+ return math::pi / std::sqrt( getStiffnessN(type1, type2)/effectiveMass - a*a);
+}
+
+inline void LinearSpringDashpot::setStiffnessAndDamping(const size_t type1,
+ const size_t type2,
+ const real_t coefficientOfRestitution,
+ const real_t collisionTime,
+ const real_t poissonsRatio,
+ const real_t effectiveMass)
+{
+ const real_t kappa = real_t(2) * ( real_t(1) - poissonsRatio ) / ( real_t(2) - poissonsRatio ) ;
+ const real_t lnDryResCoeff = std::log(coefficientOfRestitution);
+
+ setStiffnessN(type1, type2, effectiveMass * ( math::pi * math::pi + lnDryResCoeff * lnDryResCoeff ) / (collisionTime * collisionTime) );
+ setDampingN (type1, type2, - real_t(2) * effectiveMass * lnDryResCoeff / collisionTime );
+
+ setStiffnessT(type1, type2, kappa * getStiffnessN(type1, type2) );
+ setDampingT (type1, type2, std::sqrt(kappa) * getDampingN(type1, type2) );
+}
+
} //namespace kernel
} //namespace mesa_pd
} //namespace walberla
\ No newline at end of file