BrownianMotionForce.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 
26 #include "BrownianMotionForce.H"
27 #include "mathematicalConstants.H"
28 #include "fundamentalConstants.H"
29 #include "demandDrivenData.H"
30 #include "turbulenceModel.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  (
64  this->owner().U().group()
65  );
66 
67  if (obr.foundObject<turbulenceModel>(turbName))
68  {
69  const turbulenceModel& model =
70  obr.lookupObject<turbulenceModel>(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_(readScalar(this->coeffs().lookup("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 scalar dt,
164  const scalar mass,
165  const scalar Re,
166  const scalar muc
167 ) const
168 {
169  forceSuSp value(Zero, 0.0);
170 
171  const scalar dp = p.d();
172  const scalar Tc = p.Tc();
173 
174  const scalar alpha = 2.0*lambda_/dp;
175  const scalar cc = 1.0 + alpha*(1.257 + 0.4*exp(-1.1/alpha));
176 
177  // Boltzmann constant
178  const scalar kb = physicoChemical::k.value();
179 
180  scalar f = 0;
181  if (turbulence_)
182  {
183  const label celli = p.cell();
184  const volScalarField& k = *kPtr_;
185  const scalar kc = k[celli];
186  const scalar Dp = kb*Tc*cc/(3*mathematical::pi*muc*dp);
187  f = sqrt(2.0*sqr(kc)*sqr(Tc)/(Dp*dt));
188  }
189  else
190  {
191  const scalar s0 =
192  216*muc*kb*Tc/(sqr(mathematical::pi)*pow5(dp)*sqr(p.rho())*cc);
193  f = mass*sqrt(mathematical::pi*s0/dt);
194  }
195 
196 
197  // To generate a cubic distribution (3 independent directions) :
198  // const scalar sqrt2 = sqrt(2.0);
199  // for (direction dir = 0; dir < vector::nComponents; dir++)
200  // {
201  // const scalar x = rndGen_.sample01<scalar>();
202  // const scalar eta = sqrt2*erfInv(2*x - 1.0);
203  // value.Su()[dir] = f*eta;
204  // }
205 
206 
207  // To generate a spherical distribution:
208 
209  cachedRandom& rnd = this->owner().rndGen();
210 
211  const scalar theta = rnd.sample01<scalar>()*twoPi;
212  const scalar u = 2*rnd.sample01<scalar>() - 1;
213 
214  const scalar a = sqrt(1 - sqr(u));
215  const vector dir(a*cos(theta), a*sin(theta), u);
216 
217  value.Su() = f*mag(rnd.GaussNormal<scalar>())*dir;
218 
219  return value;
220 }
221 
222 
223 // ************************************************************************* //
Collection of constants.
const char *const group
Group name for atomic constants.
dictionary dict
U
Definition: pEqn.H:83
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)
const vector & Su() const
Return const access to the explicit contribution [kg.m/s2].
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:137
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
dimensionedSymmTensor sqr(const dimensionedVector &dv)
bool isTmp() const
Return true if this is really a temporary object.
Definition: tmpI.H:146
ThermalDiffusivity< CompressibleTurbulenceModel< fluidThermo > > turbulenceModel
dimensionedScalar sqrt(const dimensionedScalar &ds)
const fvMesh & mesh() const
Return the mesh database.
Abstract base class for particle forces.
Definition: ParticleForce.H:53
Random number generator.
Definition: cachedRandom.H:63
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
Type GaussNormal()
Return a sample whose components are normally distributed.
Type sample01()
Return a sample whose components lie in the range 0-1.
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar exp(const dimensionedScalar &ds)
static const word propertiesName
Default name of the turbulence properties dictionary.
static word groupName(Name name, const word &group)
const Type & value() const
Return const reference to value.
static const zero Zero
Definition: zero.H:91
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
errorManip< error > abort(error &err)
Definition: errorManip.H:131
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:218
dimensionedScalar sin(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:262
labelList f(nPoints)
virtual void cacheFields(const bool store)
Cache fields.
const dimensionedScalar k
Boltzmann constant.
Template functions to aid in the implementation of demand driven data.
virtual forceSuSp calcCoupled(const typename CloudType::parcelType &p, const scalar dt, const scalar mass, const scalar Re, const scalar muc) const
Calculate the coupled force.
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:198
A class for managing temporary objects.
Definition: PtrList.H:53
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
void deleteDemandDrivenData(DataPtr &dataPtr)
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69
Calculates particle Brownian motion force.