Newer
Older
Henry Weller
committed
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
Henry Weller
committed
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "alphatWallBoilingWallFunctionFvPatchScalarField.H"
#include "fvPatchFieldMapper.H"
#include "addToRunTimeSelectionTable.H"
#include "twoPhaseSystem.H"
#include "phaseSystem.H"
#include "ThermalPhaseChangePhaseSystem.H"
#include "MomentumTransferPhaseSystem.H"
#include "compressibleTurbulenceModel.H"
#include "ThermalDiffusivity.H"
#include "PhaseCompressibleTurbulenceModel.H"
#include "saturationModel.H"
#include "wallFvPatch.H"
#include "uniformDimensionedFields.H"
#include "mathematicalConstants.H"
Henry Weller
committed
using namespace Foam::constant::mathematical;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
<
Foam::compressible::
alphatWallBoilingWallFunctionFvPatchScalarField::phaseType
>
Foam::compressible::
alphatWallBoilingWallFunctionFvPatchScalarField::phaseTypeNames_
{
{ phaseType::vaporPhase, "vapor" },
{ phaseType::liquidPhase, "liquid" },
};
Henry Weller
committed
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace compressible
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
alphatWallBoilingWallFunctionFvPatchScalarField::
alphatWallBoilingWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF),
phaseType_(liquidPhase),
Henry Weller
committed
relax_(0.5),
AbyV_(p.size(), 0),
alphatConv_(p.size(), 0),
dDep_(p.size(), 1e-5),
qq_(p.size(), 0),
partitioningModel_(nullptr),
nucleationSiteModel_(nullptr),
departureDiamModel_(nullptr),
departureFreqModel_(nullptr)
Henry Weller
committed
{
AbyV_ = this->patch().magSf();
forAll(AbyV_, facei)
Henry Weller
committed
{
const label faceCelli = this->patch().faceCells()[facei];
Henry Weller
committed
AbyV_[facei] /= iF.mesh().V()[faceCelli];
}
}
alphatWallBoilingWallFunctionFvPatchScalarField::
alphatWallBoilingWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF, dict),
phaseType_(phaseTypeNames_.lookup("phaseType", dict)),
Henry Weller
committed
relax_(dict.lookupOrDefault<scalar>("relax", 0.5)),
AbyV_(p.size(), 0),
alphatConv_(p.size(), 0),
dDep_(p.size(), 1e-5),
qq_(p.size(), 0),
partitioningModel_(nullptr),
nucleationSiteModel_(nullptr),
departureDiamModel_(nullptr),
departureFreqModel_(nullptr)
Henry Weller
committed
{
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
switch (phaseType_)
{
case vaporPhase:
{
partitioningModel_ =
wallBoilingModels::partitioningModel::New
(
dict.subDict("partitioningModel")
);
dmdt_ = 0;
break;
}
case liquidPhase:
{
partitioningModel_ =
wallBoilingModels::partitioningModel::New
(
dict.subDict("partitioningModel")
);
nucleationSiteModel_ =
wallBoilingModels::nucleationSiteModel::New
(
dict.subDict("nucleationSiteModel")
);
departureDiamModel_ =
wallBoilingModels::departureDiameterModel::New
(
dict.subDict("departureDiamModel")
);
departureFreqModel_ =
wallBoilingModels::departureFrequencyModel::New
(
dict.subDict("departureFreqModel")
);
if (dict.found("dDep"))
{
dDep_ = scalarField("dDep", dict, p.size());
}
if (dict.found("qQuenching"))
{
qq_ = scalarField("qQuenching", dict, p.size());
}
break;
}
}
Henry Weller
committed
if (dict.found("alphatConv"))
{
alphatConv_ = scalarField("alphatConv", dict, p.size());
}
AbyV_ = this->patch().magSf();
forAll(AbyV_, facei)
Henry Weller
committed
{
const label faceCelli = this->patch().faceCells()[facei];
Henry Weller
committed
AbyV_[facei] /= iF.mesh().V()[faceCelli];
}
}
alphatWallBoilingWallFunctionFvPatchScalarField::
alphatWallBoilingWallFunctionFvPatchScalarField
(
const alphatWallBoilingWallFunctionFvPatchScalarField& psf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
(
psf,
p,
iF,
mapper
),
relax_(psf.relax_),
AbyV_(psf.AbyV_),
alphatConv_(psf.alphatConv_, mapper),
dDep_(psf.dDep_, mapper),
qq_(psf.qq_, mapper),
partitioningModel_(psf.partitioningModel_),
nucleationSiteModel_(psf.nucleationSiteModel_),
departureDiamModel_(psf.departureDiamModel_),
departureFreqModel_(psf.departureFreqModel_)
Henry Weller
committed
{}
alphatWallBoilingWallFunctionFvPatchScalarField::
alphatWallBoilingWallFunctionFvPatchScalarField
(
const alphatWallBoilingWallFunctionFvPatchScalarField& psf
)
:
alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf),
relax_(psf.relax_),
AbyV_(psf.AbyV_),
alphatConv_(psf.alphatConv_),
dDep_(psf.dDep_),
partitioningModel_(psf.partitioningModel_),
nucleationSiteModel_(psf.nucleationSiteModel_),
departureDiamModel_(psf.departureDiamModel_),
departureFreqModel_(psf.departureFreqModel_)
Henry Weller
committed
{}
alphatWallBoilingWallFunctionFvPatchScalarField::
alphatWallBoilingWallFunctionFvPatchScalarField
(
const alphatWallBoilingWallFunctionFvPatchScalarField& psf,
const DimensionedField<scalar, volMesh>& iF
)
:
alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf, iF),
relax_(psf.relax_),
AbyV_(psf.AbyV_),
alphatConv_(psf.alphatConv_),
dDep_(psf.dDep_),
partitioningModel_(psf.partitioningModel_),
nucleationSiteModel_(psf.nucleationSiteModel_),
departureDiamModel_(psf.departureDiamModel_),
departureFreqModel_(psf.departureFreqModel_)
Henry Weller
committed
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
// Check that partitioningModel has been constructed
if (!partitioningModel_.valid())
{
FatalErrorInFunction
<< "partitioningModel has not been constructed!"
<< abort(FatalError);
}
Henry Weller
committed
// Lookup the fluid model
const ThermalPhaseChangePhaseSystem
Henry Weller
committed
<
MomentumTransferPhaseSystem<twoPhaseSystem>
>& fluid =
refCast
Henry Weller
committed
<
const ThermalPhaseChangePhaseSystem
<
MomentumTransferPhaseSystem<twoPhaseSystem>
>
>
(
db().lookupObject<phaseSystem>("phaseProperties")
);
const label patchi = patch().index();
switch (phaseType_)
Henry Weller
committed
{
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
case vaporPhase:
{
const phaseModel& vapor
(
fluid.phase1().name() == internalField().group()
? fluid.phase1()
: fluid.phase2()
);
const phaseModel& liquid(fluid.otherPhase(vapor));
// Liquid phase fraction at the wall
const scalarField liquidw(liquid.boundaryField()[patchi]);
// Vapor Liquid phase fraction at the wall
const scalarField vaporw(vapor.boundaryField()[patchi]);
const scalarField fLiquid
(
partitioningModel_->fLiquid(liquidw)
);
operator==
(
calcAlphat(*this)*(1 - fLiquid)/max(vaporw, scalar(1e-8))
);
break;
}
case liquidPhase:
{
// Check that nucleationSiteModel has been constructed
if (!nucleationSiteModel_.valid())
{
FatalErrorInFunction
<< "nucleationSiteModel has not been constructed!"
<< abort(FatalError);
}
// Check that departureDiameterModel has been constructed
if (!departureDiamModel_.valid())
{
FatalErrorInFunction
<< "departureDiameterModel has not been constructed!"
<< abort(FatalError);
}
// Check that nucleationSiteModel has been constructed
if (!departureFreqModel_.valid())
{
FatalErrorInFunction
<< "departureFrequencyModel has not been constructed!"
<< abort(FatalError);
}
const phaseModel& liquid
(
fluid.phase1().name() == internalField().group()
? fluid.phase1()
: fluid.phase2()
);
const phaseModel& vapor(fluid.otherPhase(liquid));
// Retrieve turbulence properties from model
const phaseCompressibleTurbulenceModel& turbModel =
liquid.turbulence();
const tmp<scalarField> tnutw = turbModel.nut(patchi);
const scalar Cmu25(pow025(Cmu_));
const scalarField& y = turbModel.y()[patchi];
const tmp<scalarField> tmuw = turbModel.mu(patchi);
const scalarField& muw = tmuw();
const tmp<scalarField> talphaw = liquid.thermo().alpha(patchi);
const scalarField& alphaw = talphaw();
const tmp<volScalarField> tk = turbModel.k();
const volScalarField& k = tk();
const fvPatchScalarField& kw = k.boundaryField()[patchi];
const fvPatchVectorField& Uw =
turbModel.U().boundaryField()[patchi];
const scalarField magUp(mag(Uw.patchInternalField() - Uw));
const scalarField magGradUw(mag(Uw.snGrad()));
const fvPatchScalarField& rhow =
turbModel.rho().boundaryField()[patchi];
const fvPatchScalarField& hew =
liquid.thermo().he().boundaryField()[patchi];
const fvPatchScalarField& Tw =
liquid.thermo().T().boundaryField()[patchi];
const scalarField Tc(Tw.patchInternalField());
const scalarField uTau(Cmu25*sqrt(kw));
const scalarField yPlus(uTau*y/(muw/rhow));
const scalarField Pr(muw/alphaw);
// Molecular-to-turbulent Prandtl number ratio
const scalarField Prat(Pr/Prt_);
// Thermal sublayer thickness
const scalarField P(this->Psmooth(Prat));
const scalarField yPlusTherm(this->yPlusTherm(P, Prat));
const fvPatchScalarField& rhoLiquidw =
liquid.turbulence().rho().boundaryField()[patchi];
const fvPatchScalarField& rhoVaporw =
vapor.turbulence().rho().boundaryField()[patchi];
tmp<volScalarField> tCp = liquid.thermo().Cp();
const volScalarField& Cp = tCp();
const fvPatchScalarField& Cpw = Cp.boundaryField()[patchi];
// Saturation temperature
const tmp<volScalarField> tTsat =
fluid.saturation().Tsat(liquid.thermo().p());
const volScalarField& Tsat = tTsat();
const fvPatchScalarField& Tsatw(Tsat.boundaryField()[patchi]);
const scalarField Tsatc(Tsatw.patchInternalField());
const fvPatchScalarField& pw =
liquid.thermo().p().boundaryField()[patchi];
const scalarField L
(
vapor.thermo().he(pw,Tsatc,patchi)-hew.patchInternalField()
);
// Liquid phase fraction at the wall
const scalarField liquidw(liquid.boundaryField()[patchi]);
const scalarField fLiquid(partitioningModel_->fLiquid(liquidw));
Henry Weller
committed
// Convective thermal diffusivity
alphatConv_ = calcAlphat(alphatConv_);
for (label i=0; i<10; i++)
{
// Liquid temperature at y+=250 is estimated from logarithmic
// thermal wall function (Koncar, Krepper & Egorov, 2005)
const scalarField Tplus_y250(Prt_*(log(E_*250)/kappa_ + P));
const scalarField Tplus(Prt_*(log(E_*yPlus)/kappa_ + P));
scalarField Tl(Tw - (Tplus_y250/Tplus)*(Tw - Tc));
Tl = max(Tc - 40, Tl);
// Nucleation site density:
const scalarField N
(
nucleationSiteModel_->N
(
liquid,
vapor,
patchi,
Tl,
Tsatw,
L
)
);
// Bubble departure diameter:
dDep_ = departureDiamModel_->dDeparture
liquid,
vapor,
patchi,
Tl,
Tsatw,
L
);
// Bubble departure frequency:
const scalarField fDep
(
departureFreqModel_->fDeparture
(
liquid,
vapor,
patchi,
)
);
// Area fractions:
// Del Valle & Kenning (1985)
const scalarField Ja
(
rhoLiquidw*Cpw*(Tsatw - Tl)/(rhoVaporw*L)
);
const scalarField Al
(
fLiquid*4.8*exp( min(-Ja/80,log(VGREAT)))
);
const scalarField A2(min(pi*sqr(dDep_)*N*Al/4, scalar(1)));
const scalarField A1(max(1 - A2, scalar(1e-4)));
const scalarField A2E(min(pi*sqr(dDep_)*N*Al/4, scalar(5)));
// Volumetric mass source in the near wall cell due to the
// wall boiling
dmdt_ =
(1 - relax_)*dmdt_
+ relax_*(1.0/6.0)*A2E*dDep_*rhoVaporw*fDep*AbyV_;
// Volumetric source in the near wall cell due to the wall
// boiling
mDotL_ = dmdt_*L;
// Quenching heat transfer coefficient
const scalarField hQ
(
2*(alphaw*Cpw)*fDep*sqrt((0.8/fDep)/(pi*alphaw/rhow))
);
// Quenching heat flux
qq_ = (A2*hQ*max(Tw - Tl, scalar(0)));
// Effective thermal diffusivity that corresponds to the
// calculated convective, quenching and evaporative heat fluxes
operator==
(
(
A1*alphatConv_
+ (qq_ + qe())/max(hew.snGrad(), scalar(1e-16))
)
/max(liquidw, scalar(1e-8))
);
scalarField TsupPrev(max((Tw - Tsatw),scalar(0)));
const_cast<fvPatchScalarField&>(Tw).evaluate();
scalarField TsupNew(max((Tw - Tsatw),scalar(0)));
scalar maxErr(max(mag(TsupPrev - TsupNew)));
if (maxErr < 1e-1)
{
if (i > 0)
{
Info<< "Wall boiling wall function iterations: "
<< i + 1 << endl;
}
break;
}
if (debug)
{
const scalarField qc
(
fLiquid*A1*(alphatConv_ + alphaw)*hew.snGrad()
);
const scalarField qEff
(
liquidw*(*this + alphaw)*hew.snGrad()
);
Info<< " L: " << gMin(L) << " - " << gMax(L) << endl;
Info<< " Tl: " << gMin(Tl) << " - " << gMax(Tl) << endl;
Info<< " N: " << gMin(N) << " - " << gMax(N) << endl;
Info<< " dDep_: " << gMin(dDep_) << " - "
<< gMax(dDep_) << endl;
Info<< " fDep: " << gMin(fDep) << " - "
<< gMax(fDep) << endl;
Info<< " Al: " << gMin(Al) << " - " << gMax(Al) << endl;
Info<< " A1: " << gMin(A1) << " - " << gMax(A1) << endl;
Info<< " A2: " << gMin(A2) << " - " << gMax(A2) << endl;
Info<< " A2E: " << gMin(A2E) << " - "
<< gMax(A2E) << endl;
Info<< " dmdtW: " << gMin(dmdt_) << " - "
<< gMax(dmdt_) << endl;
Info<< " qc: " << gMin(qc) << " - " << gMax(qc) << endl;
Info<< " qq: " << gMin(fLiquid*qq()) << " - "
<< gMax(fLiquid*qq()) << endl;
Info<< " qe: " << gMin(fLiquid*qe()) << " - "
<< gMax(fLiquid*qe()) << endl;
Info<< " qEff: " << gMin(qEff) << " - "
<< gMax(qEff) << endl;
Info<< " alphat: " << gMin(*this) << " - "
<< gMax(*this) << endl;
Info<< " alphatConv: " << gMin(alphatConv_)
<< " - " << gMax(alphatConv_) << endl;
}
}
break;
}
default:
{
FatalErrorInFunction
<< "Unknown phase type. Valid types are: "
<< phaseTypeNames_ << nl << exit(FatalError);
}
Henry Weller
committed
}
fixedValueFvPatchScalarField::updateCoeffs();
}
void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fvPatchField<scalar>::write(os);
os.writeEntry("phaseType", phaseTypeNames_[phaseType_]);
os.writeEntry("relax", relax_);
switch (phaseType_)
{
case vaporPhase:
{
os.beginBlock("partitioningModel");
partitioningModel_->write(os);
os.endBlock();
break;
}
case liquidPhase:
{
os.beginBlock("partitioningModel");
partitioningModel_->write(os);
os.endBlock();
os.beginBlock("nucleationSiteModel");
nucleationSiteModel_->write(os);
os.endBlock();
os.beginBlock("departureDiamModel");
departureDiamModel_->write(os);
os.endBlock();
os.beginBlock("departureFreqModel");
departureFreqModel_->write(os);
os.endBlock();
break;
}
}
Henry Weller
committed
dmdt_.writeEntry("dmdt", os);
dDep_.writeEntry("dDep", os);
qq_.writeEntry("qQuenching", os);
Henry Weller
committed
alphatConv_.writeEntry("alphatConv", os);
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField
(
fvPatchScalarField,
alphatWallBoilingWallFunctionFvPatchScalarField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace compressible
} // End namespace Foam
// ************************************************************************* //