PairSpringSliderDashpot.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  // then rMin into minimum R,
69  // 1/RMin = 1/rMin + 1/rMin,
70  // RMin = rMin/2 = dMin/4
71 
72  RMin /= 4.0;
73 
74  // Multiply by two to create the worst-case relative velocity
75 
76  UMagMax = 2*UMagMax;
77 }
78 
79 
80 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
81 
82 template<class CloudType>
84 (
85  const dictionary& dict,
86  CloudType& cloud
87 )
88 :
89  PairModel<CloudType>(dict, cloud, typeName),
90  Estar_(),
91  Gstar_(),
92  alpha_(readScalar(this->coeffDict().lookup("alpha"))),
93  b_(readScalar(this->coeffDict().lookup("b"))),
94  mu_(readScalar(this->coeffDict().lookup("mu"))),
95  cohesionEnergyDensity_
96  (
97  readScalar(this->coeffDict().lookup("cohesionEnergyDensity"))
98  ),
99  cohesion_(false),
100  collisionResolutionSteps_
101  (
102  readScalar(this->coeffDict().lookup("collisionResolutionSteps"))
103  ),
104  volumeFactor_(1.0),
105  useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
106 {
107  if (useEquivalentSize_)
108  {
109  volumeFactor_ = readScalar(this->coeffDict().lookup("volumeFactor"));
110  }
111 
112  scalar nu = this->owner().constProps().poissonsRatio();
113 
114  scalar E = this->owner().constProps().youngsModulus();
115 
116  Estar_ = E/(2.0*(1.0 - sqr(nu)));
117 
118  scalar G = E/(2.0*(1.0 + nu));
119 
120  Gstar_ = G/(2.0*(2.0 - nu));
121 
122  cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL);
123 }
124 
125 
126 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
127 
128 template<class CloudType>
130 {}
131 
132 
133 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
134 
135 template<class CloudType>
137 {
138  return true;
139 }
140 
141 
142 template<class CloudType>
144 {
145  if (!(this->owner().size()))
146  {
147  return 1;
148  }
149 
150  scalar RMin;
151  scalar rhoMax;
152  scalar UMagMax;
153 
154  findMinMaxProperties(RMin, rhoMax, UMagMax);
155 
156  // Note: pi^(7/5)*(5/4)^(2/5) = 5.429675
157  scalar minCollisionDeltaT =
158  5.429675
159  *RMin
160  *pow(rhoMax/(Estar_*sqrt(UMagMax) + VSMALL), 0.4)
161  /collisionResolutionSteps_;
162 
163  return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
164 }
165 
166 
167 template<class CloudType>
169 (
170  typename CloudType::parcelType& pA,
171  typename CloudType::parcelType& pB
172 ) const
173 {
174  vector r_AB = (pA.position() - pB.position());
175 
176  scalar dAEff = pA.d();
177 
178  if (useEquivalentSize_)
179  {
180  dAEff *= cbrt(pA.nParticle()*volumeFactor_);
181  }
182 
183  scalar dBEff = pB.d();
184 
185  if (useEquivalentSize_)
186  {
187  dBEff *= cbrt(pB.nParticle()*volumeFactor_);
188  }
189 
190  scalar r_AB_mag = mag(r_AB);
191 
192  scalar normalOverlapMag = 0.5*(dAEff + dBEff) - r_AB_mag;
193 
194  if (normalOverlapMag > 0)
195  {
196  //Particles in collision
197 
198  vector rHat_AB = r_AB/(r_AB_mag + VSMALL);
199 
200  vector U_AB = pA.U() - pB.U();
201 
202  // Effective radius
203  scalar R = 0.5*dAEff*dBEff/(dAEff + dBEff);
204 
205  // Effective mass
206  scalar M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());
207 
208  scalar kN = (4.0/3.0)*sqrt(R)*Estar_;
209 
210  scalar etaN = alpha_*sqrt(M*kN)*pow025(normalOverlapMag);
211 
212  // Normal force
213  vector fN_AB =
214  rHat_AB
215  *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
216 
217  // Cohesion force, energy density multiplied by the area of
218  // particle-particle overlap
219  if (cohesion_)
220  {
221  fN_AB +=
222  -cohesionEnergyDensity_
223  *overlapArea(dAEff/2.0, dBEff/2.0, r_AB_mag)
224  *rHat_AB;
225  }
226 
227  pA.f() += fN_AB;
228  pB.f() += -fN_AB;
229 
230  vector USlip_AB =
231  U_AB - (U_AB & rHat_AB)*rHat_AB
232  + (pA.omega() ^ (dAEff/2*-rHat_AB))
233  - (pB.omega() ^ (dBEff/2*rHat_AB));
234 
235  scalar deltaT = this->owner().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 =
280  -kT*tangentialOverlapMag
281  *tangentialOverlap_AB/tangentialOverlapMag
282  - etaT*USlip_AB;
283  }
284 
285  pA.f() += fT_AB;
286  pB.f() += -fT_AB;
287 
288  pA.torque() += (dAEff/2*-rHat_AB) ^ fT_AB;
289  pB.torque() += (dBEff/2*rHat_AB) ^ -fT_AB;
290  }
291  }
292 }
293 
294 
295 // ************************************************************************* //
#define readScalar
Definition: doubleScalar.C:38
PtrList< dimensionedScalar > rhoMax(fluidRegions.size())
dictionary dict
virtual bool controlsTimestep() const
Whether the PairModel has a timestep limit that will.
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 > &)
const dimensionedScalar G
Newtonian constant of gravitation.
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
Pair forces between particles colliding with a spring, slider, damper model.
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.
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensionedScalar cbrt(const dimensionedScalar &ds)
Templated pair interaction class.
Definition: PairCollision.H:48
static const zero Zero
Definition: zero.H:91
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:217
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
#define R(A, B, C, D, E, F, K, M)
virtual label nSubCycles() const
For PairModels that control the timestep, calculate the.
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual ~PairSpringSliderDashpot()
Destructor.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
volScalarField & nu
#define M(I)