WallSpringSliderDashpot.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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  scalar kN,
81  bool cohesion
82 ) const
83 {
84  vector r_PW = p.position() - site;
85 
86  vector U_PW = p.U() - data.wallData();
87 
88  scalar r_PW_mag = mag(r_PW);
89 
90  scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);
91 
92  vector rHat_PW = r_PW/(r_PW_mag + VSMALL);
93 
94  scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
95 
96  vector fN_PW =
97  rHat_PW
98  *(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW));
99 
100  // Cohesion force, energy density multiplied by the area of wall/particle
101  // overlap
102  if (cohesion)
103  {
104  fN_PW +=
105  -cohesionEnergyDensity_
106  *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
107  *rHat_PW;
108  }
109 
110  p.f() += fN_PW;
111 
112  vector USlip_PW =
113  U_PW - (U_PW & rHat_PW)*rHat_PW
114  + (p.omega() ^ (pREff*-rHat_PW));
115 
116  scalar deltaT = this->owner().mesh().time().deltaTValue();
117 
118  vector& tangentialOverlap_PW =
119  p.collisionRecords().matchWallRecord(-r_PW, pREff).collisionData();
120 
121  tangentialOverlap_PW += USlip_PW*deltaT;
122 
123  scalar tangentialOverlapMag = mag(tangentialOverlap_PW);
124 
125  if (tangentialOverlapMag > VSMALL)
126  {
127  scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar_;
128 
129  scalar etaT = etaN;
130 
131  // Tangential force
132  vector fT_PW;
133 
134  if (kT*tangentialOverlapMag > mu_*mag(fN_PW))
135  {
136  // Tangential force greater than sliding friction,
137  // particle slips
138 
139  fT_PW = -mu_*mag(fN_PW)*USlip_PW/mag(USlip_PW);
140 
141  tangentialOverlap_PW = Zero;
142  }
143  else
144  {
145  fT_PW =
146  -kT*tangentialOverlapMag
147  *tangentialOverlap_PW/tangentialOverlapMag
148  - etaT*USlip_PW;
149  }
150 
151  p.f() += fT_PW;
152 
153  p.torque() += (pREff*-rHat_PW) ^ fT_PW;
154  }
155 }
156 
157 
158 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
159 
160 template<class CloudType>
162 (
163  const dictionary& dict,
164  CloudType& cloud
165 )
166 :
167  WallModel<CloudType>(dict, cloud, typeName),
168  Estar_(),
169  Gstar_(),
170  alpha_(readScalar(this->coeffDict().lookup("alpha"))),
171  b_(readScalar(this->coeffDict().lookup("b"))),
172  mu_(readScalar(this->coeffDict().lookup("mu"))),
173  cohesionEnergyDensity_
174  (
175  readScalar(this->coeffDict().lookup("cohesionEnergyDensity"))
176  ),
177  cohesion_(false),
178  collisionResolutionSteps_
179  (
180  readScalar
181  (
182  this->coeffDict().lookup("collisionResolutionSteps")
183  )
184  ),
185  volumeFactor_(1.0),
186  useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
187 {
188  if (useEquivalentSize_)
189  {
190  volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
191  }
192 
193  scalar nu = readScalar(this->coeffDict().lookup("poissonsRatio"));
194 
195  scalar E = readScalar(this->coeffDict().lookup("youngsModulus"));
196 
197  scalar pNu = this->owner().constProps().poissonsRatio();
198 
199  scalar pE = this->owner().constProps().youngsModulus();
200 
201  Estar_ = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
202 
203  Gstar_ = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
204 
205  cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL);
206 }
207 
208 
209 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
210 
211 template<class CloudType>
213 {}
214 
215 
216 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
217 
218 template<class CloudType>
220 (
221  const typename CloudType::parcelType& p
222 ) const
223 {
224  if (useEquivalentSize_)
225  {
226  return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
227  }
228  else
229  {
230  return p.d()/2;
231  }
232 }
233 
234 
235 template<class CloudType>
237 {
238  return true;
239 }
240 
241 
242 template<class CloudType>
244 {
245  if (!(this->owner().size()))
246  {
247  return 1;
248  }
249 
250  scalar rMin;
251  scalar rhoMax;
252  scalar UMagMax;
253 
254  findMinMaxProperties(rMin, rhoMax, UMagMax);
255 
256  // Note: pi^(7/5)*(5/4)^(2/5) = 5.429675
257  scalar minCollisionDeltaT =
258  5.429675
259  *rMin
260  *pow(rhoMax/(Estar_*sqrt(UMagMax) + VSMALL), 0.4)
261  /collisionResolutionSteps_;
262 
263  return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
264 }
265 
266 
267 template<class CloudType>
269 (
270  typename CloudType::parcelType& p,
271  const List<point>& flatSitePoints,
272  const List<WallSiteData<vector>>& flatSiteData,
273  const List<point>& sharpSitePoints,
274  const List<WallSiteData<vector>>& sharpSiteData
275 ) const
276 {
277  scalar pREff = this->pREff(p);
278 
279  scalar kN = (4.0/3.0)*sqrt(pREff)*Estar_;
280 
281  forAll(flatSitePoints, siteI)
282  {
283  evaluateWall
284  (
285  p,
286  flatSitePoints[siteI],
287  flatSiteData[siteI],
288  pREff,
289  kN,
290  cohesion_
291  );
292  }
293 
294  forAll(sharpSitePoints, siteI)
295  {
296  // Treating sharp sites like flat sites, except suppress cohesion
297 
298  evaluateWall
299  (
300  p,
301  sharpSitePoints[siteI],
302  sharpSiteData[siteI],
303  pREff,
304  kN,
305  false
306  );
307  }
308 }
309 
310 
311 // ************************************************************************* //
#define readScalar
Definition: doubleScalar.C:38
PtrList< dimensionedScalar > rhoMax(fluidRegions.size())
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
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)
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, y/n, t/f, or none.
Definition: Switch.H:60
Forces between particles and walls, interacting with a spring, slider, damper model.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
virtual ~WallSpringSliderDashpot()
Destructor.
virtual scalar pREff(const typename CloudType::parcelType &p) const
Return the effective radius for a particle for the model.
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensionedScalar cbrt(const dimensionedScalar &ds)
static const zero Zero
Definition: zero.H:91
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
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:217
virtual label nSubCycles() const
For WallModels that control the timestep, calculate the.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
WallSpringSliderDashpot(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual bool controlsTimestep() const
Whether the WallModel has a timestep limit that will.
vector point
Point is a vector.
Definition: point.H:41
dimensioned< scalar > mag(const dimensioned< Type > &)
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
volScalarField & nu