33 template<
class CloudType>
45 for (
const typename CloudType::parcelType&
p : this->owner())
51 if (useEquivalentSize_)
53 dEff *=
cbrt(
p.nParticle()*volumeFactor_);
56 RMin =
min(dEff, RMin);
83 template<
class CloudType>
93 alpha_(this->coeffDict().getScalar(
"alpha")),
94 b_(this->coeffDict().getScalar(
"b")),
95 mu_(this->coeffDict().getScalar(
"mu")),
96 cohesionEnergyDensity_
98 this->coeffDict().getScalar(
"cohesionEnergyDensity")
101 collisionResolutionSteps_
103 this->coeffDict().getScalar(
"collisionResolutionSteps")
106 useEquivalentSize_(
Switch(this->coeffDict().
lookup(
"useEquivalentSize")))
108 if (useEquivalentSize_)
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>
139 if (!(this->owner().size()))
148 findMinMaxProperties(RMin,
rhoMax, UMagMax);
151 scalar minCollisionDeltaT =
155 /collisionResolutionSteps_;
157 return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
161 template<
class CloudType>
168 vector r_AB = (pA.position() - pB.position());
170 scalar dAEff = pA.d();
172 if (useEquivalentSize_)
174 dAEff *=
cbrt(pA.nParticle()*volumeFactor_);
177 scalar dBEff = pB.d();
179 if (useEquivalentSize_)
181 dBEff *=
cbrt(pB.nParticle()*volumeFactor_);
184 scalar r_AB_mag =
mag(r_AB);
186 scalar normalOverlapMag = 0.5*(dAEff + dBEff) - r_AB_mag;
188 if (normalOverlapMag > 0)
193 scalar forceCoeff = this->forceCoeff(pA, pB);
195 vector rHat_AB = r_AB/(r_AB_mag + VSMALL);
197 vector U_AB = pA.U() - pB.U();
200 scalar
R = 0.5*dAEff*dBEff/(dAEff + dBEff);
203 scalar
M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());
205 scalar kN = (4.0/3.0)*
sqrt(
R)*Estar_;
207 scalar etaN = alpha_*
sqrt(
M*kN)*
pow025(normalOverlapMag);
212 *(kN*
pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
219 -cohesionEnergyDensity_
220 *overlapArea(dAEff/2.0, dBEff/2.0, r_AB_mag)
230 U_AB - (U_AB & rHat_AB)*rHat_AB
231 - ((dAEff/2*pA.omega() + dBEff/2*pB.omega()) ^ rHat_AB);
233 scalar deltaT = this->owner().mesh().time().deltaTValue();
235 vector& tangentialOverlap_AB =
236 pA.collisionRecords().matchPairRecord
242 vector& tangentialOverlap_BA =
243 pB.collisionRecords().matchPairRecord
249 vector deltaTangentialOverlap_AB = USlip_AB*deltaT;
251 tangentialOverlap_AB += deltaTangentialOverlap_AB;
252 tangentialOverlap_BA += -deltaTangentialOverlap_AB;
254 scalar tangentialOverlapMag =
mag(tangentialOverlap_AB);
256 if (tangentialOverlapMag > VSMALL)
258 scalar kT = 8.0*
sqrt(
R*normalOverlapMag)*Gstar_;
265 if (kT*tangentialOverlapMag > mu_*
mag(fN_AB))
270 fT_AB = -mu_*
mag(fN_AB)*USlip_AB/
mag(USlip_AB);
272 tangentialOverlap_AB =
Zero;
273 tangentialOverlap_BA =
Zero;
277 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;
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A list of keyword definitions, which are a keyword followed by a number of values (eg...
const dimensionedScalar G
Newtonian constant of gravitation.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
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, any/none. Also accepts 0/1 as a string and shortcuts t/f, y/n.
Pair forces between particles colliding with a spring, slider, damper model.
Lookup type of boundary radiation properties.
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.
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, IOobjectOption::readOption readOpt=IOobjectOption::MUST_READ) const
Find entry and assign to T val. FatalIOError if it is found and the number of tokens is incorrect...
virtual bool controlsTimestep() const
Whether the PairModel has a timestep limit that will.
A cloud is a registry collection of lagrangian particles.
const CloudType & owner() const
Return the owner cloud object.
dimensionedScalar cbrt(const dimensionedScalar &ds)
Templated pair interaction class.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
ParcelType parcelType
Type of parcel the cloud was instantiated for.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const dictionary & coeffDict() const
Return the coefficients dictionary.
#define R(A, B, C, D, E, F, K, M)
Templated base class for dsmc cloud.
virtual label nSubCycles() const
For PairModels that control the timestep, calculate the.
const List< typename ParcelType::constantProperties > & constProps() const
Return all of the constant properties.
const dimensionedScalar rhoMax
static constexpr const zero Zero
Global zero (0)