WallLocalSpringSliderDashpot.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-2020 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 
28 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
29 
30 template<class CloudType>
32 (
33  scalar& rMin,
34  scalar& rhoMax,
35  scalar& UMagMax
36 ) const
37 {
38  rMin = vGreat;
39  rhoMax = -vGreat;
40  UMagMax = -vGreat;
41 
42  forAllConstIter(typename CloudType, this->owner(), iter)
43  {
44  const typename CloudType::parcelType& p = iter();
45 
46  // Finding minimum diameter to avoid excessive arithmetic
47 
48  scalar dEff = p.d();
49 
50  if (useEquivalentSize_)
51  {
52  dEff *= cbrt(p.nParticle()*volumeFactor_);
53  }
54 
55  rMin = min(dEff, rMin);
56 
57  rhoMax = max(p.rho(), rhoMax);
58 
59  UMagMax = max
60  (
61  mag(p.U()) + mag(p.omega())*dEff/2,
62  UMagMax
63  );
64  }
65 
66  // Transform the minimum diameter into minimum radius
67  // rMin = dMin/2
68 
69  rMin /= 2.0;
70 }
71 
72 
73 template<class CloudType>
75 (
76  typename CloudType::parcelType& p,
77  const point& site,
78  const WallSiteData<vector>& data,
79  scalar pREff,
80  bool cohesion
81 ) const
82 {
83  // wall patch index
84  label wPI = patchMap_[data.patchIndex()];
85 
86  // data for this patch
87  scalar Estar = Estar_[wPI];
88  scalar Gstar = Gstar_[wPI];
89  scalar alpha = alpha_[wPI];
90  scalar b = b_[wPI];
91  scalar mu = mu_[wPI];
92  scalar cohesionEnergyDensity = cohesionEnergyDensity_[wPI];
93  cohesion = cohesion && cohesion_[wPI];
94 
95  vector r_PW = p.position() - site;
96 
97  vector U_PW = p.U() - data.wallData();
98 
99  scalar r_PW_mag = mag(r_PW);
100 
101  scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);
102 
103  vector rHat_PW = r_PW/(r_PW_mag + vSmall);
104 
105  scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
106 
107  scalar etaN = alpha*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
108 
109  vector fN_PW =
110  rHat_PW
111  *(kN*pow(normalOverlapMag, b) - etaN*(U_PW & rHat_PW));
112 
113  // Cohesion force, energy density multiplied by the area of wall/particle
114  // overlap
115  if (cohesion)
116  {
117  fN_PW +=
118  -cohesionEnergyDensity
119  *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
120  *rHat_PW;
121  }
122 
123  p.f() += fN_PW;
124 
125  vector USlip_PW =
126  U_PW - (U_PW & rHat_PW)*rHat_PW
127  + (p.omega() ^ (pREff*-rHat_PW));
128 
129  scalar deltaT = this->owner().mesh().time().deltaTValue();
130 
131  vector& tangentialOverlap_PW =
132  p.collisionRecords().matchWallRecord(-r_PW, pREff).collisionData();
133 
134  tangentialOverlap_PW += USlip_PW*deltaT;
135 
136  scalar tangentialOverlapMag = mag(tangentialOverlap_PW);
137 
138  if (tangentialOverlapMag > vSmall)
139  {
140  scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar;
141 
142  scalar etaT = etaN;
143 
144  // Tangential force
145  vector fT_PW;
146 
147  if (kT*tangentialOverlapMag > mu*mag(fN_PW))
148  {
149  // Tangential force greater than sliding friction,
150  // particle slips
151 
152  fT_PW = -mu*mag(fN_PW)*USlip_PW/mag(USlip_PW);
153 
154  tangentialOverlap_PW = Zero;
155  }
156  else
157  {
158  fT_PW =
159  -kT*tangentialOverlapMag
160  *tangentialOverlap_PW/tangentialOverlapMag
161  - etaT*USlip_PW;
162  }
163 
164  p.f() += fT_PW;
165 
166  p.torque() += (pREff*-rHat_PW) ^ fT_PW;
167  }
168 }
169 
170 
171 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
172 
173 template<class CloudType>
175 (
176  const dictionary& dict,
177  CloudType& cloud
178 )
179 :
180  WallModel<CloudType>(dict, cloud, typeName),
181  Estar_(),
182  Gstar_(),
183  alpha_(),
184  b_(),
185  mu_(),
186  cohesionEnergyDensity_(),
187  cohesion_(),
188  patchMap_(),
189  maxEstarIndex_(-1),
190  collisionResolutionSteps_
191  (
192  this->coeffDict().template lookup<scalar>("collisionResolutionSteps")
193  ),
194  volumeFactor_(1.0),
195  useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
196 {
197  if (useEquivalentSize_)
198  {
199  volumeFactor_ =
200  this->coeffDict().template lookup<scalar>("volumeFactor");
201  }
202 
203  scalar pNu = this->owner().constProps().poissonsRatio();
204 
205  scalar pE = this->owner().constProps().youngsModulus();
206 
207  const polyMesh& mesh = cloud.mesh();
208 
209  const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
210 
211  patchMap_.setSize(bMesh.size(), -1);
212 
213  DynamicList<label> wallPatchIndices;
214 
215  forAll(bMesh, patchi)
216  {
217  if (isA<wallPolyPatch>(bMesh[patchi]))
218  {
219  wallPatchIndices.append(bMesh[patchi].index());
220  }
221  }
222 
223  label nWallPatches = wallPatchIndices.size();
224 
225  Estar_.setSize(nWallPatches);
226  Gstar_.setSize(nWallPatches);
227  alpha_.setSize(nWallPatches);
228  b_.setSize(nWallPatches);
229  mu_.setSize(nWallPatches);
230  cohesionEnergyDensity_.setSize(nWallPatches);
231  cohesion_.setSize(nWallPatches);
232 
233  scalar maxEstar = -great;
234 
235  forAll(wallPatchIndices, wPI)
236  {
237  const dictionary& patchCoeffDict
238  (
239  this->coeffDict().subDict(bMesh[wallPatchIndices[wPI]].name())
240  );
241 
242  patchMap_[wallPatchIndices[wPI]] = wPI;
243 
244  scalar nu = patchCoeffDict.template lookup<scalar>("poissonsRatio");
245 
246  scalar E = patchCoeffDict.template lookup<scalar>("youngsModulus");
247 
248  Estar_[wPI] = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
249 
250  Gstar_[wPI] = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
251 
252  alpha_[wPI] = patchCoeffDict.template lookup<scalar>("alpha");
253 
254  b_[wPI] = patchCoeffDict.template lookup<scalar>("b");
255 
256  mu_[wPI] = patchCoeffDict.template lookup<scalar>("mu");
257 
258  cohesionEnergyDensity_[wPI] =
259  patchCoeffDict.lookup<scalar>("cohesionEnergyDensity");
260 
261  cohesion_[wPI] = (mag(cohesionEnergyDensity_[wPI]) > vSmall);
262 
263  if (Estar_[wPI] > maxEstar)
264  {
265  maxEstarIndex_ = wPI;
266 
267  maxEstar = Estar_[wPI];
268  }
269  }
270 }
271 
272 
273 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
274 
275 template<class CloudType>
277 {}
278 
279 
280 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
281 
282 template<class CloudType>
284 (
285  const typename CloudType::parcelType& p
286 ) const
287 {
288  if (useEquivalentSize_)
289  {
290  return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
291  }
292  else
293  {
294  return p.d()/2;
295  }
296 }
297 
298 
299 template<class CloudType>
301 {
302  return true;
303 }
304 
305 
306 template<class CloudType>
308 {
309  if (!(this->owner().size()))
310  {
311  return 1;
312  }
313 
314  scalar rMin;
315  scalar rhoMax;
316  scalar UMagMax;
317 
318  findMinMaxProperties(rMin, rhoMax, UMagMax);
319 
320  // Note: pi^(7/5)*(5/4)^(2/5) = 5.429675
321  scalar minCollisionDeltaT =
322  5.429675
323  *rMin
324  *pow(rhoMax/(Estar_[maxEstarIndex_]*sqrt(UMagMax) + vSmall), 0.4)
325  /collisionResolutionSteps_;
326 
327  return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
328 }
329 
330 
331 template<class CloudType>
333 (
334  typename CloudType::parcelType& p,
335  const List<point>& flatSitePoints,
336  const List<WallSiteData<vector>>& flatSiteData,
337  const List<point>& sharpSitePoints,
338  const List<WallSiteData<vector>>& sharpSiteData
339 ) const
340 {
341  scalar pREff = this->pREff(p);
342 
343  forAll(flatSitePoints, siteI)
344  {
345  evaluateWall
346  (
347  p,
348  flatSitePoints[siteI],
349  flatSiteData[siteI],
350  pREff,
351  true
352  );
353  }
354 
355  forAll(sharpSitePoints, siteI)
356  {
357  // Treating sharp sites like flat sites, except suppress cohesion
358 
359  evaluateWall
360  (
361  p,
362  sharpSitePoints[siteI],
363  sharpSiteData[siteI],
364  pREff,
365  false
366  );
367  }
368 }
369 
370 
371 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
WallLocalSpringSliderDashpot(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
DSMCCloud< dsmcParcel > CloudType
virtual label nSubCycles() const
For WallModels that control the timestep, calculate the.
virtual scalar pREff(const typename CloudType::parcelType &p) const
Return the effective radius for a particle for the model.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
dimensionedScalar sqrt(const dimensionedScalar &ds)
Forces between particles and walls, interacting with a spring, slider, damper model.
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.
Definition: Switch.H:60
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
virtual bool controlsTimestep() const
Whether the WallModel has a timestep limit that will.
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
dimensionedScalar cbrt(const dimensionedScalar &ds)
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
static const zero Zero
Definition: zero.H:97
Templated wall interaction class.
Definition: PairCollision.H:51
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
Stores the patch ID and templated data to represent a collision with a wall to be passed to the wall ...
Definition: WallSiteData.H:50
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
Foam::polyBoundaryMesh.
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:221
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label patchi
vector point
Point is a vector.
Definition: point.H:41
const polyBoundaryMesh & bMesh
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844