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