30 template<
class CloudType>
44 const typename CloudType::parcelType& p = iter();
50 if (useEquivalentSize_)
52 dEff *=
cbrt(p.nParticle()*volumeFactor_);
55 RMin =
min(dEff, RMin);
57 rhoMax =
max(p.rho(), rhoMax);
61 mag(p.U()) +
mag(p.omega())*dEff/2,
82 template<
class CloudType>
92 alpha_(this->coeffDict().
template lookup<scalar>(
"alpha")),
93 b_(this->coeffDict().
template lookup<scalar>(
"b")),
94 mu_(this->coeffDict().
template lookup<scalar>(
"mu")),
95 cohesionEnergyDensity_
97 this->coeffDict().
template lookup<scalar>(
"cohesionEnergyDensity")
100 collisionResolutionSteps_
102 this->coeffDict().
template lookup<scalar>(
"collisionResolutionSteps")
105 useEquivalentSize_(
Switch(this->coeffDict().
lookup(
"useEquivalentSize")))
107 if (useEquivalentSize_)
110 this->coeffDict().template lookup<scalar>(
"volumeFactor");
113 scalar nu = this->owner().constProps().poissonsRatio();
115 scalar E = this->owner().constProps().youngsModulus();
117 Estar_ = E/(2.0*(1.0 -
sqr(nu)));
119 scalar
G = E/(2.0*(1.0 + nu));
121 Gstar_ = G/(2.0*(2.0 - nu));
123 cohesion_ = (
mag(cohesionEnergyDensity_) > vSmall);
129 template<
class CloudType>
136 template<
class CloudType>
143 template<
class CloudType>
146 if (!(this->owner().size()))
155 findMinMaxProperties(RMin, rhoMax, UMagMax);
158 scalar minCollisionDeltaT =
161 *
pow(rhoMax/(Estar_*
sqrt(UMagMax) + vSmall), 0.4)
162 /collisionResolutionSteps_;
164 return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
168 template<
class CloudType>
175 vector r_AB = (pA.position() - pB.position());
177 scalar dAEff = pA.d();
179 if (useEquivalentSize_)
181 dAEff *=
cbrt(pA.nParticle()*volumeFactor_);
184 scalar dBEff = pB.d();
186 if (useEquivalentSize_)
188 dBEff *=
cbrt(pB.nParticle()*volumeFactor_);
191 scalar r_AB_mag =
mag(r_AB);
193 scalar normalOverlapMag = 0.5*(dAEff + dBEff) - r_AB_mag;
195 if (normalOverlapMag > 0)
199 vector rHat_AB = r_AB/(r_AB_mag + vSmall);
201 vector U_AB = pA.U() - pB.U();
204 scalar
R = 0.5*dAEff*dBEff/(dAEff + dBEff);
207 scalar
M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());
209 scalar kN = (4.0/3.0)*
sqrt(R)*Estar_;
211 scalar etaN = alpha_*
sqrt(M*kN)*
pow025(normalOverlapMag);
216 *(kN*
pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
223 -cohesionEnergyDensity_
224 *overlapArea(dAEff/2.0, dBEff/2.0, r_AB_mag)
232 U_AB - (U_AB & rHat_AB)*rHat_AB
233 - ((dAEff/2*pA.omega() + dBEff/2*pB.omega()) ^ rHat_AB);
235 scalar deltaT = this->owner().mesh().time().deltaTValue();
237 vector& tangentialOverlap_AB =
238 pA.collisionRecords().matchPairRecord
244 vector& tangentialOverlap_BA =
245 pB.collisionRecords().matchPairRecord
251 vector deltaTangentialOverlap_AB = USlip_AB*deltaT;
253 tangentialOverlap_AB += deltaTangentialOverlap_AB;
254 tangentialOverlap_BA += -deltaTangentialOverlap_AB;
256 scalar tangentialOverlapMag =
mag(tangentialOverlap_AB);
258 if (tangentialOverlapMag > vSmall)
260 scalar kT = 8.0*
sqrt(R*normalOverlapMag)*Gstar_;
267 if (kT*tangentialOverlapMag > mu_*
mag(fN_AB))
272 fT_AB = -mu_*
mag(fN_AB)*USlip_AB/
mag(USlip_AB);
274 tangentialOverlap_AB =
Zero;
275 tangentialOverlap_BA =
Zero;
279 fT_AB = - kT*tangentialOverlap_AB - etaT*USlip_AB;
285 pA.torque() += (dAEff/2*-rHat_AB) ^ fT_AB;
286 pB.torque() += (dBEff/2*rHat_AB) ^ -fT_AB;
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
DSMCCloud< dsmcParcel > CloudType
A list of keyword definitions, which are a keyword followed by any number of values (e...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensionedScalar pow025(const dimensionedScalar &ds)
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none/any.
Pair forces between particles colliding with a spring, slider, damper model.
PairSpringSliderDashpot(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
virtual void evaluatePair(typename CloudType::parcelType &pA, typename CloudType::parcelType &pB) const
Calculate the pair interaction between parcels.
stressControl lookup("compactNormalStress") >> compactNormalStress
virtual bool controlsTimestep() const
Whether the PairModel has a timestep limit that will.
dimensionedScalar cbrt(const dimensionedScalar &ds)
Templated pair interaction class.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
ParcelType parcelType
Type of parcel the cloud was instantiated for.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
#define R(A, B, C, D, E, F, K, M)
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual ~PairSpringSliderDashpot()
Destructor.
Templated base class for dsmc cloud.
virtual label nSubCycles() const
For PairModels that control the timestep, calculate the.
const dimensionedScalar & G
Newtonian constant of gravitation.