PairSpringSliderDashpot.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
27 #include "polyMesh.H"
28 
29 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
30 
31 template<class CloudType>
33 (
34  scalar& RMin,
35  scalar& rhoMax,
36  scalar& UMagMax
37 ) const
38 {
39  RMin = vGreat;
40  rhoMax = -vGreat;
41  UMagMax = -vGreat;
42 
43  forAllConstIter(typename CloudType, this->owner(), iter)
44  {
45  const typename CloudType::parcelType& p = iter();
46 
47  // Finding minimum diameter to avoid excessive arithmetic
48 
49  scalar dEff = p.d();
50 
51  if (useEquivalentSize_)
52  {
53  dEff *= cbrt(p.nParticle()*volumeFactor_);
54  }
55 
56  RMin = min(dEff, RMin);
57 
58  rhoMax = max(p.rho(), rhoMax);
59 
60  UMagMax = max
61  (
62  mag(p.U()) + mag(p.omega())*dEff/2,
63  UMagMax
64  );
65  }
66 
67  // Transform the minimum diameter into minimum radius
68  // rMin = dMin/2
69  // then rMin into minimum R,
70  // 1/RMin = 1/rMin + 1/rMin,
71  // RMin = rMin/2 = dMin/4
72 
73  RMin /= 4.0;
74 
75  // Multiply by two to create the worst-case relative velocity
76 
77  UMagMax = 2*UMagMax;
78 }
79 
80 
81 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
82 
83 template<class CloudType>
85 (
86  const dictionary& dict,
88 )
89 :
90  PairModel<CloudType>(dict, cloud, typeName),
91  Estar_(),
92  Gstar_(),
93  alpha_(this->coeffDict().template lookup<scalar>("alpha")),
94  b_(this->coeffDict().template lookup<scalar>("b")),
95  mu_(this->coeffDict().template lookup<scalar>("mu")),
96  cohesionEnergyDensity_
97  (
98  this->coeffDict().template lookup<scalar>("cohesionEnergyDensity")
99  ),
100  cohesion_(false),
101  collisionResolutionSteps_
102  (
103  this->coeffDict().template lookup<scalar>("collisionResolutionSteps")
104  ),
105  volumeFactor_(1.0),
106  useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
107 {
108  if (useEquivalentSize_)
109  {
110  volumeFactor_ =
111  this->coeffDict().template lookup<scalar>("volumeFactor");
112  }
113 
114  scalar nu = this->owner().constProps().poissonsRatio();
115 
116  scalar E = this->owner().constProps().youngsModulus();
117 
118  Estar_ = E/(2.0*(1.0 - sqr(nu)));
119 
120  scalar G = E/(2.0*(1.0 + nu));
121 
122  Gstar_ = G/(2.0*(2.0 - nu));
123 
124  cohesion_ = (mag(cohesionEnergyDensity_) > vSmall);
125 }
126 
127 
128 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
129 
130 template<class CloudType>
132 {}
133 
134 
135 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
136 
137 template<class CloudType>
139 {
140  return true;
141 }
142 
143 
144 template<class CloudType>
146 {
147  if (!this->owner().size())
148  {
149  return 1;
150  }
151 
152  scalar RMin, rhoMax, UMagMax;
153  findMinMaxProperties(RMin, rhoMax, UMagMax);
154 
155  // Note: pi^(7/5)*(5/4)^(2/5) = 5.429675
156  const scalar minCollisionDeltaT =
157  5.429675
158  *RMin
159  *pow(rhoMax/(Estar_*sqrt(UMagMax) + small), 0.4)
160  /collisionResolutionSteps_;
161 
162  return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
163 }
164 
165 
166 template<class CloudType>
168 (
169  typename CloudType::parcelType& pA,
170  typename CloudType::parcelType& pB
171 ) const
172 {
173  const polyMesh& mesh = this->owner().mesh();
174 
175  vector r_AB = pA.position(mesh) - pB.position(mesh);
176 
177  scalar dAEff = pA.d();
178 
179  if (useEquivalentSize_)
180  {
181  dAEff *= cbrt(pA.nParticle()*volumeFactor_);
182  }
183 
184  scalar dBEff = pB.d();
185 
186  if (useEquivalentSize_)
187  {
188  dBEff *= cbrt(pB.nParticle()*volumeFactor_);
189  }
190 
191  scalar r_AB_mag = mag(r_AB);
192 
193  scalar normalOverlapMag = 0.5*(dAEff + dBEff) - r_AB_mag;
194 
195  if (normalOverlapMag > 0)
196  {
197  // Particles in collision
198 
199  vector rHat_AB = r_AB/(r_AB_mag + vSmall);
200 
201  vector U_AB = pA.U() - pB.U();
202 
203  // Effective radius
204  scalar R = 0.5*dAEff*dBEff/(dAEff + dBEff);
205 
206  // Effective mass
207  scalar M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());
208 
209  scalar kN = (4.0/3.0)*sqrt(R)*Estar_;
210 
211  scalar etaN = alpha_*sqrt(M*kN)*pow025(normalOverlapMag);
212 
213  // Normal force
214  vector fN_AB =
215  rHat_AB
216  *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
217 
218  // Cohesion force, energy density multiplied by the area of
219  // particle-particle overlap
220  if (cohesion_)
221  {
222  fN_AB +=
223  -cohesionEnergyDensity_
224  *overlapArea(dAEff/2.0, dBEff/2.0, r_AB_mag)
225  *rHat_AB;
226  }
227 
228  pA.f() += fN_AB;
229  pB.f() += -fN_AB;
230 
231  vector USlip_AB =
232  U_AB - (U_AB & rHat_AB)*rHat_AB
233  - ((dAEff/2*pA.omega() + dBEff/2*pB.omega()) ^ rHat_AB);
234 
235  scalar deltaT = mesh.time().deltaTValue();
236 
237  vector& tangentialOverlap_AB =
238  pA.collisionRecords().matchPairRecord
239  (
240  pB.origProc(),
241  pB.origId()
242  ).collisionData();
243 
244  vector& tangentialOverlap_BA =
245  pB.collisionRecords().matchPairRecord
246  (
247  pA.origProc(),
248  pA.origId()
249  ).collisionData();
250 
251  vector deltaTangentialOverlap_AB = USlip_AB*deltaT;
252 
253  tangentialOverlap_AB += deltaTangentialOverlap_AB;
254  tangentialOverlap_BA += -deltaTangentialOverlap_AB;
255 
256  scalar tangentialOverlapMag = mag(tangentialOverlap_AB);
257 
258  if (tangentialOverlapMag > vSmall)
259  {
260  scalar kT = 8.0*sqrt(R*normalOverlapMag)*Gstar_;
261 
262  scalar etaT = etaN;
263 
264  // Tangential force
265  vector fT_AB;
266 
267  if (kT*tangentialOverlapMag > mu_*mag(fN_AB))
268  {
269  // Tangential force greater than sliding friction,
270  // particle slips
271 
272  fT_AB = -mu_*mag(fN_AB)*USlip_AB/mag(USlip_AB);
273 
274  tangentialOverlap_AB = Zero;
275  tangentialOverlap_BA = Zero;
276  }
277  else
278  {
279  fT_AB = - kT*tangentialOverlap_AB - etaT*USlip_AB;
280  }
281 
282  pA.f() += fT_AB;
283  pB.f() += -fT_AB;
284 
285  pA.torque() += (dAEff/2*-rHat_AB) ^ fT_AB;
286  pB.torque() += (dBEff/2*rHat_AB) ^ -fT_AB;
287  }
288  }
289 }
290 
291 
292 // ************************************************************************* //
#define M(I)
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:80
const List< typename ParcelType::constantProperties > & constProps() const
Return all of the constant properties.
Definition: DSMCCloudI.H:95
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:225
Templated pair interaction class.
Definition: PairModel.H:53
const dictionary & coeffDict() const
Return the coefficients dictionary.
Definition: PairModel.C:70
const CloudType & owner() const
Return the owner cloud object.
Definition: PairModel.C:55
Pair forces between particles colliding with a spring, slider, damper model.
virtual label nSubCycles() const
For PairModels that control the timestep, calculate the.
virtual ~PairSpringSliderDashpot()
Destructor.
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.
virtual bool controlsTimestep() const
Whether the PairModel has a timestep limit that will.
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
Definition: Switch.H:61
scalar deltaTValue() const
Return time step value.
Definition: TimeStateI.H:34
A cloud is a collection of lagrangian particles.
Definition: cloud.H:55
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
const Time & time() const
Return time.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
const dimensionedScalar G
Newtonian constant of gravitation.
static const zero Zero
Definition: zero.H:97
DSMCCloud< dsmcParcel > CloudType
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
static scalar R(const scalar a, const scalar x)
Definition: invIncGamma.C:102
dimensionedScalar cbrt(const dimensionedScalar &ds)
dimensionedScalar pow025(const dimensionedScalar &ds)
dictionary dict
volScalarField & p