WallSpringSliderDashpot.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 "mathematicalConstants.H"
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 template<class CloudType>
34 (
35  scalar& rMin,
36  scalar& rhoMax,
37  scalar& UMagMax
38 ) const
39 {
40  rMin = vGreat;
41  rhoMax = -vGreat;
42  UMagMax = -vGreat;
43 
44  forAllConstIter(typename CloudType, this->owner(), iter)
45  {
46  const typename CloudType::parcelType& p = iter();
47 
48  // Finding minimum diameter to avoid excessive arithmetic
49 
50  scalar dEff = p.d();
51 
52  if (useEquivalentSize_)
53  {
54  dEff *= cbrt(p.nParticle()*volumeFactor_);
55  }
56 
57  rMin = min(dEff, rMin);
58 
59  rhoMax = max(p.rho(), rhoMax);
60 
61  UMagMax = max
62  (
63  mag(p.U()) + mag(p.omega())*dEff/2,
64  UMagMax
65  );
66  }
67 
68  // Transform the minimum diameter into minimum radius
69  // rMin = dMin/2
70 
71  rMin /= 2.0;
72 }
73 
74 
75 template<class CloudType>
77 (
78  typename CloudType::parcelType& p,
79  const point& site,
80  const WallSiteData<vector>& data,
81  scalar pREff,
82  scalar kN,
83  bool cohesion
84 ) const
85 {
86  const polyMesh& mesh = this->owner().mesh();
87 
88  vector r_PW = p.position(mesh) - site;
89 
90  vector U_PW = p.U() - data.wallData();
91 
92  scalar r_PW_mag = mag(r_PW);
93 
94  scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);
95 
96  vector rHat_PW = r_PW/(r_PW_mag + vSmall);
97 
98  scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
99 
100  vector fN_PW =
101  rHat_PW
102  *(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW));
103 
104  // Cohesion force, energy density multiplied by the area of wall/particle
105  // overlap
106  if (cohesion)
107  {
108  fN_PW +=
109  -cohesionEnergyDensity_
110  *constant::mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
111  *rHat_PW;
112  }
113 
114  p.f() += fN_PW;
115 
116  vector USlip_PW =
117  U_PW - (U_PW & rHat_PW)*rHat_PW
118  + (p.omega() ^ (pREff*-rHat_PW));
119 
120  scalar deltaT = this->owner().mesh().time().deltaTValue();
121 
122  vector& tangentialOverlap_PW =
123  p.collisionRecords().matchWallRecord(-r_PW, pREff).collisionData();
124 
125  tangentialOverlap_PW += USlip_PW*deltaT;
126 
127  scalar tangentialOverlapMag = mag(tangentialOverlap_PW);
128 
129  if (tangentialOverlapMag > vSmall)
130  {
131  scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar_;
132 
133  scalar etaT = etaN;
134 
135  // Tangential force
136  vector fT_PW;
137 
138  if (kT*tangentialOverlapMag > mu_*mag(fN_PW))
139  {
140  // Tangential force greater than sliding friction,
141  // particle slips
142 
143  fT_PW = -mu_*mag(fN_PW)*USlip_PW/mag(USlip_PW);
144 
145  tangentialOverlap_PW = Zero;
146  }
147  else
148  {
149  fT_PW = - kT*tangentialOverlap_PW - etaT*USlip_PW;
150  }
151 
152  p.f() += fT_PW;
153 
154  p.torque() += (pREff*-rHat_PW) ^ fT_PW;
155  }
156 }
157 
158 
159 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
160 
161 template<class CloudType>
163 (
164  const dictionary& dict,
166 )
167 :
168  WallModel<CloudType>(dict, cloud, typeName),
169  Estar_(),
170  Gstar_(),
171  alpha_(this->coeffDict().template lookup<scalar>("alpha")),
172  b_(this->coeffDict().template lookup<scalar>("b")),
173  mu_(this->coeffDict().template lookup<scalar>("mu")),
174  cohesionEnergyDensity_
175  (
176  this->coeffDict().template lookup<scalar>("cohesionEnergyDensity")
177  ),
178  cohesion_(false),
179  collisionResolutionSteps_
180  (
181  this->coeffDict().template lookup<scalar>("collisionResolutionSteps")
182  ),
183  volumeFactor_(1.0),
184  useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
185 {
186  if (useEquivalentSize_)
187  {
188  volumeFactor_ =
189  this->coeffDict().template lookup<scalar>("volumeFactor");
190  }
191 
192  scalar nu = this->coeffDict().template lookup<scalar>("poissonsRatio");
193 
194  scalar E = this->coeffDict().template lookup<scalar>("youngsModulus");
195 
196  scalar pNu = this->owner().constProps().poissonsRatio();
197 
198  scalar pE = this->owner().constProps().youngsModulus();
199 
200  Estar_ = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
201 
202  Gstar_ = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
203 
204  cohesion_ = (mag(cohesionEnergyDensity_) > vSmall);
205 }
206 
207 
208 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
209 
210 template<class CloudType>
212 {}
213 
214 
215 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
216 
217 template<class CloudType>
219 (
220  const typename CloudType::parcelType& p
221 ) const
222 {
223  if (useEquivalentSize_)
224  {
225  return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
226  }
227  else
228  {
229  return p.d()/2;
230  }
231 }
232 
233 
234 template<class CloudType>
236 {
237  return true;
238 }
239 
240 
241 template<class CloudType>
243 {
244  if (!this->owner().size())
245  {
246  return 1;
247  }
248 
249  scalar rMin, rhoMax, UMagMax;
250  findMinMaxProperties(rMin, rhoMax, UMagMax);
251 
252  // Note: pi^(7/5)*(5/4)^(2/5) = 5.429675
253  const scalar minCollisionDeltaT =
254  5.429675
255  *rMin
256  *pow(rhoMax/(Estar_*sqrt(UMagMax) + small), 0.4)
257  /collisionResolutionSteps_;
258 
259  return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
260 }
261 
262 
263 template<class CloudType>
265 (
266  typename CloudType::parcelType& p,
267  const List<point>& flatSitePoints,
268  const List<WallSiteData<vector>>& flatSiteData,
269  const List<point>& sharpSitePoints,
270  const List<WallSiteData<vector>>& sharpSiteData
271 ) const
272 {
273  scalar pREff = this->pREff(p);
274 
275  scalar kN = (4.0/3.0)*sqrt(pREff)*Estar_;
276 
277  forAll(flatSitePoints, siteI)
278  {
279  evaluateWall
280  (
281  p,
282  flatSitePoints[siteI],
283  flatSiteData[siteI],
284  pREff,
285  kN,
286  cohesion_
287  );
288  }
289 
290  forAll(sharpSitePoints, siteI)
291  {
292  // Treating sharp sites like flat sites, except suppress cohesion
293 
294  evaluateWall
295  (
296  p,
297  sharpSitePoints[siteI],
298  sharpSiteData[siteI],
299  pREff,
300  kN,
301  false
302  );
303  }
304 }
305 
306 
307 // ************************************************************************* //
#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
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
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
Definition: Switch.H:61
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
Forces between particles and walls, interacting with a spring, slider, damper model.
virtual ~WallSpringSliderDashpot()
Destructor.
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.
WallSpringSliderDashpot(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
virtual bool controlsTimestep() const
Whether the WallModel has a timestep limit that will.
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
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)
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