BrownianMotionForce.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-2020 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 
26 #include "BrownianMotionForce.H"
27 #include "mathematicalConstants.H"
28 #include "fundamentalConstants.H"
29 #include "demandDrivenData.H"
30 #include "momentumTransportModel.H"
31 
32 using namespace Foam::constant;
33 
34 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
35 
36 template<class CloudType>
37 Foam::scalar Foam::BrownianMotionForce<CloudType>::erfInv(const scalar y) const
38 {
39  const scalar a = 0.147;
40  scalar k = 2.0/(mathematical::pi*a) + 0.5*log(1.0 - y*y);
41  scalar h = log(1.0 - y*y)/a;
42  scalar x = sqrt(-k + sqrt(k*k - h));
43 
44  if (y < 0.0)
45  {
46  return -x;
47  }
48  else
49  {
50  return x;
51  }
52 }
53 
54 
55 template<class CloudType>
58 {
59  const objectRegistry& obr = this->owner().mesh();
60  const word turbName =
62  (
63  momentumTransportModel::typeName,
64  this->owner().U().group()
65  );
66 
67  if (obr.foundObject<momentumTransportModel>(turbName))
68  {
69  const momentumTransportModel& model =
70  obr.lookupObject<momentumTransportModel>(turbName);
71  return model.k();
72  }
73  else
74  {
76  << "Turbulence model not found in mesh database" << nl
77  << "Database objects include: " << obr.sortedToc()
78  << abort(FatalError);
79 
80  return tmp<volScalarField>(nullptr);
81  }
82 }
83 
84 
85 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
86 
87 template<class CloudType>
89 (
90  CloudType& owner,
91  const fvMesh& mesh,
92  const dictionary& dict
93 )
94 :
95  ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
96  rndGen_(owner.rndGen()),
97  lambda_(this->coeffs().template lookup<scalar>("lambda")),
98  turbulence_(readBool(this->coeffs().lookup("turbulence"))),
99  kPtr_(nullptr),
100  ownK_(false)
101 {}
102 
103 
104 template<class CloudType>
106 (
107  const BrownianMotionForce& bmf
108 )
109 :
111  rndGen_(bmf.rndGen_),
112  lambda_(bmf.lambda_),
113  turbulence_(bmf.turbulence_),
114  kPtr_(nullptr),
115  ownK_(false)
116 {}
117 
118 
119 // * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * //
120 
121 template<class CloudType>
123 {}
124 
125 
126 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
127 
128 template<class CloudType>
130 {
131  if (turbulence_)
132  {
133  if (store)
134  {
135  tmp<volScalarField> tk = kModel();
136  if (tk.isTmp())
137  {
138  kPtr_ = tk.ptr();
139  ownK_ = true;
140  }
141  else
142  {
143  kPtr_ = &tk();
144  ownK_ = false;
145  }
146  }
147  else
148  {
149  if (ownK_ && kPtr_)
150  {
151  deleteDemandDrivenData(kPtr_);
152  ownK_ = false;
153  }
154  }
155  }
156 }
157 
158 
159 template<class CloudType>
161 (
162  const typename CloudType::parcelType& p,
163  const typename CloudType::parcelType::trackingData& td,
164  const scalar dt,
165  const scalar mass,
166  const scalar Re,
167  const scalar muc
168 ) const
169 {
170  forceSuSp value(Zero, 0.0);
171 
172  const scalar dp = p.d();
173  const scalar Tc = td.Tc();
174 
175  const scalar alpha = 2.0*lambda_/dp;
176  const scalar cc = 1.0 + alpha*(1.257 + 0.4*exp(-1.1/alpha));
177 
178  // Boltzmann constant
179  const scalar kb = physicoChemical::k.value();
180 
181  scalar f = 0;
182  if (turbulence_)
183  {
184  const label celli = p.cell();
185  const volScalarField& k = *kPtr_;
186  const scalar kc = k[celli];
187  const scalar Dp = kb*Tc*cc/(3*mathematical::pi*muc*dp);
188  f = sqrt(2.0*sqr(kc)*sqr(Tc)/(Dp*dt));
189  }
190  else
191  {
192  const scalar s0 =
193  216*muc*kb*Tc/(sqr(mathematical::pi)*pow5(dp)*sqr(p.rho())*cc);
194  f = mass*sqrt(mathematical::pi*s0/dt);
195  }
196 
197 
198  // To generate a cubic distribution (3 independent directions) :
199  // const scalar sqrt2 = sqrt(2.0);
200  // for (direction dir = 0; dir < vector::nComponents; dir++)
201  // {
202  // const scalar x = rndGen_.sample01<scalar>();
203  // const scalar eta = sqrt2*erfInv(2*x - 1.0);
204  // value.Su()[dir] = f*eta;
205  // }
206 
207 
208  // To generate a spherical distribution:
209 
210  Random& rnd = this->owner().rndGen();
211 
212  const scalar theta = rnd.scalar01()*twoPi;
213  const scalar u = 2*rnd.scalar01() - 1;
214 
215  const scalar a = sqrt(1 - sqr(u));
216  const vector dir(a*cos(theta), a*sin(theta), u);
217 
218  value.Su() = f*mag(rnd.scalarNormal())*dir;
219 
220  return value;
221 }
222 
223 
224 // ************************************************************************* //
Collection of constants.
const char *const group
Group name for atomic constants.
dictionary dict
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
dimensionedScalar log(const dimensionedScalar &ds)
scalar scalarNormal()
Advance the state and return a scalar sample from a normal.
Definition: Random.C:31
const vector & Su() const
Return const access to the explicit contribution [kg m/s^2].
Definition: forceSuSpI.H:56
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
dimensionedSymmTensor sqr(const dimensionedVector &dv)
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
bool isTmp() const
Return true if this is really a temporary object.
Definition: tmpI.H:153
dimensionedScalar sqrt(const dimensionedScalar &ds)
const fvMesh & mesh() const
Return the mesh database.
Abstract base class for particle forces.
Definition: ParticleForce.H:53
bool readBool(Istream &)
Definition: boolIO.C:60
dimensionedScalar pow5(const dimensionedScalar &ds)
Fundamental dimensioned constants.
Helper container for force Su and Sp terms.
Definition: forceSuSp.H:61
CompressibleMomentumTransportModel< dynamicTransportModel > momentumTransportModel
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar exp(const dimensionedScalar &ds)
virtual forceSuSp calcCoupled(const typename CloudType::parcelType &p, const typename CloudType::parcelType::trackingData &td, const scalar dt, const scalar mass, const scalar Re, const scalar muc) const
Calculate the coupled force.
static word groupName(Name name, const word &group)
const Type & value() const
Return const reference to value.
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Random number generator.
Definition: Random.H:57
BrownianMotionForce(CloudType &owner, const fvMesh &mesh, const dictionary &dict)
Construct from mesh.
const scalar twoPi(2 *pi)
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:221
dimensionedScalar sin(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:260
labelList f(nPoints)
virtual void cacheFields(const bool store)
Cache fields.
Template functions to aid in the implementation of demand driven data.
U
Definition: pEqn.H:72
const dimensionedScalar k
Boltzmann constant.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual ~BrownianMotionForce()
Destructor.
T * ptr() const
Return tmp pointer for reuse.
Definition: tmpI.H:205
A class for managing temporary objects.
Definition: PtrList.H:53
void deleteDemandDrivenData(DataPtr &dataPtr)
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
scalar scalar01()
Advance the state and return a scalar sample from a uniform.
Definition: RandomI.H:57
Calculates particle Brownian motion force.