alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField.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) 2015 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 
28 #include "fvPatchFieldMapper.H"
29 #include "volFields.H"
31 #include "twoPhaseSystem.H"
32 #include "phaseSystem.H"
35 #include "wallFvPatch.H"
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41 namespace compressible
42 {
43 
44 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
45 
46 scalar alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::maxExp_
47  = 50.0;
48 scalar alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::tolerance_
49  = 0.01;
50 label alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::maxIters_
51  = 10;
52 
53 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
54 
55 void alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::checkType()
56 {
57  if (!isA<wallFvPatch>(patch()))
58  {
60  (
61  "alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField"
62  "::checkType()"
63  )
64  << "Patch type for patch " << patch().name() << " must be wall\n"
65  << "Current patch type is " << patch().type() << nl
66  << exit(FatalError);
67  }
68 }
69 
70 
71 tmp<scalarField>
72 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::Psmooth
73 (
74  const scalarField& Prat
75 ) const
76 {
77  return 9.24*(pow(Prat, 0.75) - 1.0)*(1.0 + 0.28*exp(-0.007*Prat));
78 }
79 
80 
81 tmp<scalarField>
82 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::yPlusTherm
83 (
84  const scalarField& P,
85  const scalarField& Prat
86 ) const
87 {
88 
89  tmp<scalarField> typtf(new scalarField(this->size()));
90  scalarField& yptf = typtf();
91 
92  forAll(yptf, faceI)
93  {
94  scalar ypt = 11.0;
95 
96  for (int i=0; i<maxIters_; i++)
97  {
98  scalar f = ypt - (log(E_*ypt)/kappa_ + P[faceI])/Prat[faceI];
99  scalar df = 1.0 - 1.0/(ypt*kappa_*Prat[faceI]);
100  scalar yptNew = ypt - f/df;
101 
102  if (yptNew < VSMALL)
103  {
104  yptf[faceI] = 0;
105  }
106  else if (mag(yptNew - ypt) < tolerance_)
107  {
108  yptf[faceI] = yptNew;
109  }
110  else
111  {
112  ypt = yptNew;
113  }
114  }
115 
116  yptf[faceI] = ypt;
117  }
118 
119  return typtf;
120 }
121 
122 
123 
124 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
125 
128 (
129  const fvPatch& p,
130  const DimensionedField<scalar, volMesh>& iF
131 )
132 :
134  Prt_(0.85),
135  Cmu_(0.09),
136  kappa_(0.41),
137  E_(9.8),
138  fixedDmdt_(0.0)
139 {
140  checkType();
141 }
142 
143 
144 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
145 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
146 (
147  const fvPatch& p,
148  const DimensionedField<scalar, volMesh>& iF,
149  const dictionary& dict
150 )
151 :
152  alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF, dict),
153  Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85)),
154  Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
155  kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
156  E_(dict.lookupOrDefault<scalar>("E", 9.8)),
157  fixedDmdt_(dict.lookupOrDefault<scalar>("fixedDmdt", 0.0))
158 {}
159 
160 
161 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
162 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
163 (
164  const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField& ptf,
165  const fvPatch& p,
166  const DimensionedField<scalar, volMesh>& iF,
167  const fvPatchFieldMapper& mapper
168 )
169 :
170  alphatPhaseChangeWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
171  Prt_(ptf.Prt_),
172  Cmu_(ptf.Cmu_),
173  kappa_(ptf.kappa_),
174  E_(ptf.E_),
175  fixedDmdt_(ptf.fixedDmdt_)
176 {}
177 
178 
179 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
180 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
181 (
182  const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField& awfpsf
183 )
184 :
185  alphatPhaseChangeWallFunctionFvPatchScalarField(awfpsf),
186  Prt_(awfpsf.Prt_),
187  Cmu_(awfpsf.Cmu_),
188  kappa_(awfpsf.kappa_),
189  E_(awfpsf.E_),
190  fixedDmdt_(awfpsf.fixedDmdt_)
191 {}
192 
193 
194 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
195 alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
196 (
197  const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField& awfpsf,
198  const DimensionedField<scalar, volMesh>& iF
199 )
200 :
201  alphatPhaseChangeWallFunctionFvPatchScalarField(awfpsf, iF),
202  Prt_(awfpsf.Prt_),
203  Cmu_(awfpsf.Cmu_),
204  kappa_(awfpsf.kappa_),
205  E_(awfpsf.E_),
206  fixedDmdt_(awfpsf.fixedDmdt_)
207 {}
208 
209 
210 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
211 
212 void alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
213 {
214  if (updated())
215  {
216  return;
217  }
218 
219  // Lookup the fluid model
220  const ThermalPhaseChangePhaseSystem
221  <
222  MomentumTransferPhaseSystem<twoPhaseSystem>
223  >& fluid =
224  refCast
225  <
226  const ThermalPhaseChangePhaseSystem
227  <
228  MomentumTransferPhaseSystem<twoPhaseSystem>
229  >
230  >
231  (
232  db().lookupObject<phaseSystem>("phaseProperties")
233  );
234 
235  const phaseModel& liquid
236  (
237  fluid.phase1().name() == dimensionedInternalField().group()
238  ? fluid.phase1()
239  : fluid.phase2()
240  );
241 
242  const label patchi = patch().index();
243 
244  // Retrieve turbulence properties from model
245  const compressibleTurbulenceModel& turbModel =
246  db().lookupObject<compressibleTurbulenceModel>
247  (
248  IOobject::groupName
249  (
250  compressibleTurbulenceModel::propertiesName,
252  )
253  );
254 
255  const scalar Cmu25 = pow025(Cmu_);
256 
257  const scalarField& y = turbModel.y()[patchi];
258 
259  const tmp<scalarField> tmuw = turbModel.mu(patchi);
260  const scalarField& muw = tmuw();
261 
262  const tmp<scalarField> talphaw = liquid.thermo().alpha(patchi);
263  const scalarField& alphaw = talphaw();
264 
265  scalarField& alphatw = *this;
266 
267  const tmp<volScalarField> tk = turbModel.k();
268  const volScalarField& k = tk();
269  const fvPatchScalarField& kw = k.boundaryField()[patchi];
270 
271  const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
272  const scalarField magUp(mag(Uw.patchInternalField() - Uw));
273  const scalarField magGradUw(mag(Uw.snGrad()));
274 
275  const fvPatchScalarField& rhow = turbModel.rho().boundaryField()[patchi];
276  const fvPatchScalarField& hew =
277  liquid.thermo().he().boundaryField()[patchi];
278 
279  const fvPatchScalarField& Tw =
280  liquid.thermo().T().boundaryField()[patchi];
281 
282  scalarField Tp(Tw.patchInternalField());
283 
284  // Heat flux [W/m2] - lagging alphatw
285  const scalarField qDot
286  (
287  (alphatw + alphaw)*hew.snGrad()
288  );
289 
290  scalarField uTau(Cmu25*sqrt(kw));
291 
292  scalarField yPlus(uTau*y/(muw/rhow));
293 
294  scalarField Pr(muw/alphaw);
295 
296  // Molecular-to-turbulent Prandtl number ratio
297  scalarField Prat(Pr/Prt_);
298 
299  // Thermal sublayer thickness
300  scalarField P(this->Psmooth(Prat));
301 
302  scalarField yPlusTherm(this->yPlusTherm(P, Prat));
303 
304  scalarField alphatConv(this->size(), 0.0);
305 
306  // Populate boundary values
307  forAll(alphatw, faceI)
308  {
309  // Evaluate new effective thermal diffusivity
310  scalar alphaEff = 0.0;
311  if (yPlus[faceI] < yPlusTherm[faceI])
312  {
313  scalar A = qDot[faceI]*rhow[faceI]*uTau[faceI]*y[faceI];
314  scalar B = qDot[faceI]*Pr[faceI]*yPlus[faceI];
315  scalar C = Pr[faceI]*0.5*rhow[faceI]*uTau[faceI]*sqr(magUp[faceI]);
316  alphaEff = A/(B + C + VSMALL);
317  }
318  else
319  {
320  scalar A = qDot[faceI]*rhow[faceI]*uTau[faceI]*y[faceI];
321  scalar B =
322  qDot[faceI]*Prt_*(1.0/kappa_*log(E_*yPlus[faceI]) + P[faceI]);
323  scalar magUc =
324  uTau[faceI]/kappa_*log(E_*yPlusTherm[faceI]) - mag(Uw[faceI]);
325  scalar C =
326  0.5*rhow[faceI]*uTau[faceI]
327  *(Prt_*sqr(magUp[faceI]) + (Pr[faceI] - Prt_)*sqr(magUc));
328  alphaEff = A/(B + C + VSMALL);
329  }
330 
331  // Update convective heat transfer turbulent thermal diffusivity
332  alphatConv[faceI] = max(0.0, alphaEff - alphaw[faceI]);
333  }
334 
335  dmdt_ = fixedDmdt_;
336 
337  operator==(alphatConv);
338 
339  fixedValueFvPatchScalarField::updateCoeffs();
340 }
341 
342 
344 (
345  Ostream& os
346 ) const
347 {
349  os.writeKeyword("Prt") << Prt_ << token::END_STATEMENT << nl;
350  os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
351  os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
352  os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
353  os.writeKeyword("fixedDmdt") << fixedDmdt_ << token::END_STATEMENT << nl;
354  dmdt_.writeEntry("dmdt", os);
355  writeEntry("value", os);
356 }
357 
358 
359 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
360 
362 (
364  alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
365 );
366 
367 
368 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
369 
370 } // End namespace compressible
371 } // End namespace Foam
372 
373 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensionedScalar Pr("Pr", dimless, laminarTransport)
scalar magUp
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
multiphaseSystem & fluid
Definition: createFields.H:10
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > mag(const dimensioned< Type > &)
const char *const group
Group name for atomic constants.
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
dimensionedScalar exp(const dimensionedScalar &ds)
Namespace for OpenFOAM.
runTime write()
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField(const fvPatch &, const DimensionedField< scalar, volMesh > &)
Construct from patch and internal field.
dictionary dict
dimensionedScalar log(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:260
#define makePatchTypeField(PatchTypeField, typePatchTypeField)
Definition: fvPatchField.H:651
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
volScalarField & p
Definition: createFields.H:51
scalar uTau
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
label patchi
dimensionedScalar pow025(const dimensionedScalar &ds)
Macros for easy insertion into run-time selection tables.
fvPatchField< scalar > fvPatchScalarField
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
rDeltaT dimensionedInternalField()
scalar yPlus
bool compressible
Definition: pEqn.H:40
error FatalError
fvPatchField< vector > fvPatchVectorField
To & refCast(From &r)
Reference type cast template function.
Definition: typeInfo.H:106
dimensionedSymmTensor sqr(const dimensionedVector &dv)
alphatPhaseChangeWallFunctionFvPatchScalarField(const fvPatch &, const DimensionedField< scalar, volMesh > &)
Construct from patch and internal field.
volScalarField scalarField(fieldObject, mesh)
volScalarField & C