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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2018-2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
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  for (const typename CloudType::parcelType& p : this->owner())
46  {
47  // Finding minimum diameter to avoid excessive arithmetic
48 
49  scalar dEff = p.d();
50 
51  if (useEquivalentSize_)
52  {
53  dEff *= cbrt(p.nParticle()*volumeFactor_);
54  }
55 
56  RMin = min(dEff, RMin);
57 
58  rhoMax = max(p.rho(), rhoMax);
59 
60  UMagMax = max
61  (
62  mag(p.U()) + mag(p.omega())*dEff/2,
63  UMagMax
64  );
65  }
66 
67  // Transform the minimum diameter into minimum radius
68  // rMin = dMin/2
69  // then rMin into minimum R,
70  // 1/RMin = 1/rMin + 1/rMin,
71  // RMin = rMin/2 = dMin/4
72 
73  RMin /= 4.0;
74 
75  // Multiply by two to create the worst-case relative velocity
76 
77  UMagMax = 2*UMagMax;
78 }
79 
80 
81 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
82 
83 template<class CloudType>
85 (
86  const dictionary& dict,
88 )
89 :
90  PairModel<CloudType>(dict, cloud, typeName),
91  Estar_(),
92  Gstar_(),
93  alpha_(this->coeffDict().getScalar("alpha")),
94  b_(this->coeffDict().getScalar("b")),
95  mu_(this->coeffDict().getScalar("mu")),
96  cohesionEnergyDensity_
97  (
98  this->coeffDict().getScalar("cohesionEnergyDensity")
99  ),
100  cohesion_(false),
101  collisionResolutionSteps_
102  (
103  this->coeffDict().getScalar("collisionResolutionSteps")
104  ),
105  volumeFactor_(1.0),
106  useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
107 {
108  if (useEquivalentSize_)
109  {
110  this->coeffDict().readEntry("volumeFactor", volumeFactor_);
111  }
112 
113  scalar nu = this->owner().constProps().poissonsRatio();
114 
115  scalar E = this->owner().constProps().youngsModulus();
116 
117  Estar_ = E/(2.0*(1.0 - sqr(nu)));
118 
119  scalar G = E/(2.0*(1.0 + nu));
120 
121  Gstar_ = G/(2.0*(2.0 - nu));
122 
123  cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL);
124 }
125 
126 
127 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
128 
129 template<class CloudType>
131 {
132  return true;
133 }
134 
135 
136 template<class CloudType>
138 {
139  if (!(this->owner().size()))
140  {
141  return 1;
142  }
143 
144  scalar RMin;
145  scalar rhoMax;
146  scalar UMagMax;
147 
148  findMinMaxProperties(RMin, rhoMax, UMagMax);
149 
150  // Note: pi^(7/5)*(5/4)^(2/5) = 5.429675
151  scalar minCollisionDeltaT =
152  5.429675
153  *RMin
154  *pow(rhoMax/(Estar_*sqrt(UMagMax) + VSMALL), 0.4)
155  /collisionResolutionSteps_;
156 
157  return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
158 }
159 
160 
161 template<class CloudType>
163 (
164  typename CloudType::parcelType& pA,
165  typename CloudType::parcelType& pB
166 ) const
167 {
168  vector r_AB = (pA.position() - pB.position());
169 
170  scalar dAEff = pA.d();
171 
172  if (useEquivalentSize_)
173  {
174  dAEff *= cbrt(pA.nParticle()*volumeFactor_);
175  }
176 
177  scalar dBEff = pB.d();
178 
179  if (useEquivalentSize_)
180  {
181  dBEff *= cbrt(pB.nParticle()*volumeFactor_);
182  }
183 
184  scalar r_AB_mag = mag(r_AB);
185 
186  scalar normalOverlapMag = 0.5*(dAEff + dBEff) - r_AB_mag;
187 
188  if (normalOverlapMag > 0)
189  {
190  // Particles in collision
191 
192  // Force coefficient
193  scalar forceCoeff = this->forceCoeff(pA, pB);
194 
195  vector rHat_AB = r_AB/(r_AB_mag + VSMALL);
196 
197  vector U_AB = pA.U() - pB.U();
198 
199  // Effective radius
200  scalar R = 0.5*dAEff*dBEff/(dAEff + dBEff);
201 
202  // Effective mass
203  scalar M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());
204 
205  scalar kN = (4.0/3.0)*sqrt(R)*Estar_;
206 
207  scalar etaN = alpha_*sqrt(M*kN)*pow025(normalOverlapMag);
208 
209  // Normal force
210  vector fN_AB =
211  rHat_AB
212  *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
213 
214  // Cohesion force, energy density multiplied by the area of
215  // particle-particle overlap
216  if (cohesion_)
217  {
218  fN_AB +=
219  -cohesionEnergyDensity_
220  *overlapArea(dAEff/2.0, dBEff/2.0, r_AB_mag)
221  *rHat_AB;
222  }
223 
224  fN_AB *= forceCoeff;
225 
226  pA.f() += fN_AB;
227  pB.f() += -fN_AB;
228 
229  vector USlip_AB =
230  U_AB - (U_AB & rHat_AB)*rHat_AB
231  - ((dAEff/2*pA.omega() + dBEff/2*pB.omega()) ^ rHat_AB);
232 
233  scalar deltaT = this->owner().mesh().time().deltaTValue();
234 
235  vector& tangentialOverlap_AB =
236  pA.collisionRecords().matchPairRecord
237  (
238  pB.origProc(),
239  pB.origId()
240  ).collisionData();
241 
242  vector& tangentialOverlap_BA =
243  pB.collisionRecords().matchPairRecord
244  (
245  pA.origProc(),
246  pA.origId()
247  ).collisionData();
248 
249  vector deltaTangentialOverlap_AB = USlip_AB*deltaT;
250 
251  tangentialOverlap_AB += deltaTangentialOverlap_AB;
252  tangentialOverlap_BA += -deltaTangentialOverlap_AB;
253 
254  scalar tangentialOverlapMag = mag(tangentialOverlap_AB);
255 
256  if (tangentialOverlapMag > VSMALL)
257  {
258  scalar kT = 8.0*sqrt(R*normalOverlapMag)*Gstar_;
259 
260  scalar etaT = etaN;
261 
262  // Tangential force
263  vector fT_AB;
264 
265  if (kT*tangentialOverlapMag > mu_*mag(fN_AB))
266  {
267  // Tangential force greater than sliding friction,
268  // particle slips
269 
270  fT_AB = -mu_*mag(fN_AB)*USlip_AB/mag(USlip_AB);
271 
272  tangentialOverlap_AB = Zero;
273  tangentialOverlap_BA = Zero;
274  }
275  else
276  {
277  fT_AB = - kT*tangentialOverlap_AB - etaT*USlip_AB;
278  }
279 
280  fT_AB *= forceCoeff;
281 
282  pA.f() += fT_AB;
283  pB.f() += -fT_AB;
284 
285  pA.torque() += (dAEff/2*-rHat_AB) ^ fT_AB;
286  pB.torque() += (dBEff/2*rHat_AB) ^ -fT_AB;
287  }
288  }
289 }
290 
291 
292 // ************************************************************************* //
Foam::PairSpringSliderDashpot::nSubCycles
virtual label nSubCycles() const
For PairModels that control the timestep, calculate the.
Definition: PairSpringSliderDashpot.C:137
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::Switch
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
Definition: Switch.H:77
Foam::constant::universal::G
const dimensionedScalar G
Newtonian constant of gravitation.
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::pow025
dimensionedScalar pow025(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:133
nu
volScalarField & nu
Definition: readMechanicalProperties.H:176
R
#define R(A, B, C, D, E, F, K, M)
PairSpringSliderDashpot.H
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::radiation::lookup
Lookup type of boundary radiation properties.
Definition: lookup.H:63
Foam::DSMCCloud
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:71
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::PairSpringSliderDashpot::PairSpringSliderDashpot
PairSpringSliderDashpot(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
Definition: PairSpringSliderDashpot.C:85
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::cloud
A cloud is a registry collection of lagrangian particles.
Definition: cloud.H:57
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
rhoMax
const dimensionedScalar rhoMax
Definition: setRegionFluidFields.H:66
Foam::Vector< scalar >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::PairModel
Templated pair interaction class.
Definition: PairCollision.H:53
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
M
#define M(I)
Foam::PairSpringSliderDashpot::controlsTimestep
virtual bool controlsTimestep() const
Whether the PairModel has a timestep limit that will.
Definition: PairSpringSliderDashpot.C:130
Foam::PairSpringSliderDashpot::evaluatePair
virtual void evaluatePair(typename CloudType::parcelType &pA, typename CloudType::parcelType &pB) const
Calculate the pair interaction between parcels.
Definition: PairSpringSliderDashpot.C:163
Foam::DSMCCloud::parcelType
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:220
Foam::cbrt
dimensionedScalar cbrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:155
Foam::PairSpringSliderDashpot
Pair forces between particles colliding with a spring, slider, damper model.
Definition: PairSpringSliderDashpot.H:59