solvers::XiFluid: New solver module for compressible premixed/partially-premixed combustion
executed with foamRun for single region simulations of foamMultiRun for
multi-region simulations. Replaces XiFoam and all the corresponding
tutorials have been updated and moved to tutorials/modules/XiFluid.
Class
Foam::solvers::XiFluid
Description
Solver module for compressible premixed/partially-premixed combustion with
turbulence modelling.
Combusting RANS code using the b-Xi two-equation model.
Xi may be obtained by either the solution of the Xi transport
equation or from an algebraic expression. Both approaches are
based on Gulder's flame speed correlation which has been shown
to be appropriate by comparison with the results from the
spectral model.
Strain effects are encorporated directly into the Xi equation
but not in the algebraic approximation. Further work need to be
done on this issue, particularly regarding the enhanced removal rate
caused by flame compression. Analysis using results of the spectral
model will be required.
For cases involving very lean Propane flames or other flames which are
very strain-sensitive, a transport equation for the laminar flame
speed is present. This equation is derived using heuristic arguments
involving the strain time scale and the strain-rate at extinction.
the transport velocity is the same as that for the Xi equation.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient and steady simulations.
Optional fvModels and fvConstraints are provided to enhance the simulation
in many ways including adding various sources, chemical reactions,
combustion, Lagrangian particles, radiation, surface film etc. and
constraining or limiting the solution.
Reference:
\verbatim
Greenshields, C. J., & Weller, H. G. (2022).
Notes on Computational Fluid Dynamics: General Principles.
CFD Direct Ltd.: Reading, UK.
\endverbatim
SourceFiles
XiFluid.C
See also
Foam::solvers::fluidSolver
Foam::solvers::isothermalFluid
This commit is contained in:
Henry Weller
@ -1,6 +1,4 @@
|
||||
EXE_INC = \
|
||||
-I.. \
|
||||
-I../ignition/lnInclude \
|
||||
-IXiModels/transport \
|
||||
-IXiModels/fixed \
|
||||
-IXiModels/algebraic \
|
||||
@ -26,6 +24,7 @@ EXE_INC = \
|
||||
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/ignition/lnInclude \
|
||||
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
|
||||
-I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
|
||||
-I$(LIB_SRC)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \
|
||||
@ -36,7 +35,6 @@ EXE_INC = \
|
||||
-I$(LIB_SRC)/triSurface/lnInclude \
|
||||
|
||||
EXE_LIBS = \
|
||||
-lXiIgnition \
|
||||
-lmeshTools \
|
||||
-lfluidThermophysicalModels \
|
||||
-lmulticomponentThermophysicalModels \
|
||||
@ -45,6 +43,7 @@ EXE_LIBS = \
|
||||
-lcompressibleMomentumTransportModels \
|
||||
-lcoupledThermophysicalTransportModels \
|
||||
-llaminarFlameSpeedModels \
|
||||
-lXiIgnition \
|
||||
-lfiniteVolume \
|
||||
-lfvModels \
|
||||
-lfvConstraints
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -1,9 +0,0 @@
|
||||
#!/bin/sh
|
||||
cd ${0%/*} || exit 1 # Run from this directory
|
||||
|
||||
wclean libso ignition
|
||||
|
||||
wclean
|
||||
wclean PDRFoam
|
||||
|
||||
#------------------------------------------------------------------------------
|
||||
@ -1,12 +0,0 @@
|
||||
#!/bin/sh
|
||||
cd ${0%/*} || exit 1 # Run from this directory
|
||||
|
||||
# Parse arguments for library compilation
|
||||
. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
|
||||
|
||||
wmake $targetType ignition
|
||||
|
||||
wmake $targetType
|
||||
wmake $targetType PDRFoam
|
||||
|
||||
#------------------------------------------------------------------------------
|
||||
@ -1,26 +0,0 @@
|
||||
{
|
||||
volScalarField& hea = thermo.he();
|
||||
|
||||
fvScalarMatrix EaEqn
|
||||
(
|
||||
fvm::ddt(rho, hea) + mvConvection->fvmDiv(phi, hea)
|
||||
+ fvc::ddt(rho, K) + fvc::div(phi, K)
|
||||
+ (
|
||||
hea.name() == "ea"
|
||||
? fvc::div(fvc::absolute(phi, rho, U), p/rho)
|
||||
: -dpdt
|
||||
)
|
||||
+ thermophysicalTransport.divq(hea)
|
||||
+ fvModels.source(rho, hea)
|
||||
);
|
||||
|
||||
EaEqn.relax();
|
||||
|
||||
fvConstraints.constrain(EaEqn);
|
||||
|
||||
EaEqn.solve();
|
||||
|
||||
fvConstraints.constrain(hea);
|
||||
|
||||
thermo.correct();
|
||||
}
|
||||
@ -1,32 +0,0 @@
|
||||
if (ign.ignited())
|
||||
{
|
||||
volScalarField& heau = thermo.heu();
|
||||
|
||||
fvScalarMatrix heauEqn
|
||||
(
|
||||
fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau)
|
||||
+ (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
|
||||
+ (
|
||||
heau.name() == "eau"
|
||||
? fvc::div(fvc::absolute(phi, rho, U), p/rho)
|
||||
*rho/thermo.rhou()
|
||||
: -dpdt*rho/thermo.rhou()
|
||||
)
|
||||
+ thermophysicalTransport.divq(heau)
|
||||
|
||||
// These terms cannot be used in partially-premixed combustion due to
|
||||
// the resultant inconsistency between ft and heau transport.
|
||||
// A possible solution would be to solve for ftu as well as ft.
|
||||
//- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
|
||||
//+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
|
||||
|
||||
==
|
||||
fvModels.source(rho, heau)
|
||||
);
|
||||
|
||||
fvConstraints.constrain(heauEqn);
|
||||
|
||||
heauEqn.solve();
|
||||
|
||||
fvConstraints.constrain(heau);
|
||||
}
|
||||
@ -1,3 +0,0 @@
|
||||
XiFoam.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/XiFoam
|
||||
@ -1,20 +0,0 @@
|
||||
fvVectorMatrix UEqn
|
||||
(
|
||||
fvm::ddt(rho, U) + fvm::div(phi, U)
|
||||
+ MRF.DDt(rho, U)
|
||||
+ turbulence->divDevTau(U)
|
||||
==
|
||||
fvModels.source(rho, U)
|
||||
);
|
||||
|
||||
UEqn.relax();
|
||||
|
||||
fvConstraints.constrain(UEqn);
|
||||
|
||||
if (pimple.momentumPredictor())
|
||||
{
|
||||
solve(UEqn == -fvc::grad(p));
|
||||
|
||||
fvConstraints.constrain(U);
|
||||
K = 0.5*magSqr(U);
|
||||
}
|
||||
@ -1,206 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
|
||||
\\/ 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/>.
|
||||
|
||||
Application
|
||||
XiFoam
|
||||
|
||||
Description
|
||||
Solver for compressible premixed/partially-premixed combustion with
|
||||
turbulence modelling.
|
||||
|
||||
Combusting RANS code using the b-Xi two-equation model.
|
||||
Xi may be obtained by either the solution of the Xi transport
|
||||
equation or from an algebraic expression. Both approaches are
|
||||
based on Gulder's flame speed correlation which has been shown
|
||||
to be appropriate by comparison with the results from the
|
||||
spectral model.
|
||||
|
||||
Strain effects are encorporated directly into the Xi equation
|
||||
but not in the algebraic approximation. Further work need to be
|
||||
done on this issue, particularly regarding the enhanced removal rate
|
||||
caused by flame compression. Analysis using results of the spectral
|
||||
model will be required.
|
||||
|
||||
For cases involving very lean Propane flames or other flames which are
|
||||
very strain-sensitive, a transport equation for the laminar flame
|
||||
speed is present. This equation is derived using heuristic arguments
|
||||
involving the strain time scale and the strain-rate at extinction.
|
||||
the transport velocity is the same as that for the Xi equation.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
#include "psiuMulticomponentThermo.H"
|
||||
#include "compressibleMomentumTransportModels.H"
|
||||
#include "RASThermophysicalTransportModel.H"
|
||||
#include "unityLewisEddyDiffusivity.H"
|
||||
#include "laminarFlameSpeed.H"
|
||||
#include "ignition.H"
|
||||
#include "Switch.H"
|
||||
#include "pimpleControl.H"
|
||||
#include "CorrectPhi.H"
|
||||
#include "fvModels.H"
|
||||
#include "fvConstraints.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
#include "postProcess.H"
|
||||
|
||||
#include "setRootCaseLists.H"
|
||||
#include "createTime.H"
|
||||
#include "createMesh.H"
|
||||
#include "createDyMControls.H"
|
||||
#include "initContinuityErrs.H"
|
||||
#include "readCombustionProperties.H"
|
||||
#include "createFields.H"
|
||||
#include "createFieldRefs.H"
|
||||
#include "createRhoUfIfPresent.H"
|
||||
|
||||
turbulence->validate();
|
||||
|
||||
#include "compressibleCourantNo.H"
|
||||
#include "setInitialDeltaT.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
Info<< "\nStarting time loop\n" << endl;
|
||||
|
||||
while (pimple.run(runTime))
|
||||
{
|
||||
#include "readDyMControls.H"
|
||||
|
||||
// Store divrhoU from the previous mesh so that it can be mapped
|
||||
// and used in correctPhi to ensure the corrected phi has the
|
||||
// same divergence
|
||||
autoPtr<volScalarField> divrhoU;
|
||||
if (correctPhi)
|
||||
{
|
||||
divrhoU = new volScalarField
|
||||
(
|
||||
"divrhoU",
|
||||
fvc::div(fvc::absolute(phi, rho, U))
|
||||
);
|
||||
}
|
||||
|
||||
#include "compressibleCourantNo.H"
|
||||
#include "setDeltaT.H"
|
||||
|
||||
fvModels.preUpdateMesh();
|
||||
|
||||
// Store momentum to set rhoUf for introduced faces.
|
||||
autoPtr<volVectorField> rhoU;
|
||||
if (rhoUf.valid())
|
||||
{
|
||||
rhoU = new volVectorField("rhoU", rho*U);
|
||||
}
|
||||
|
||||
// Update the mesh for topology change, mesh to mesh mapping
|
||||
mesh.update();
|
||||
|
||||
runTime++;
|
||||
|
||||
Info<< "Time = " << runTime.userTimeName() << nl << endl;
|
||||
|
||||
// --- Pressure-velocity PIMPLE corrector loop
|
||||
while (pimple.loop())
|
||||
{
|
||||
if (pimple.firstPimpleIter() || pimple.moveMeshOuterCorrectors())
|
||||
{
|
||||
// Move the mesh
|
||||
mesh.move();
|
||||
|
||||
if (mesh.changing())
|
||||
{
|
||||
MRF.update();
|
||||
|
||||
if (correctPhi)
|
||||
{
|
||||
#include "correctPhi.H"
|
||||
}
|
||||
|
||||
if (checkMeshCourantNo)
|
||||
{
|
||||
#include "meshCourantNo.H"
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if
|
||||
(
|
||||
!pimple.simpleRho()
|
||||
&& pimple.firstPimpleIter()
|
||||
)
|
||||
{
|
||||
#include "rhoEqn.H"
|
||||
}
|
||||
|
||||
fvModels.correct();
|
||||
|
||||
if (pimple.predictTransport())
|
||||
{
|
||||
turbulence->predict();
|
||||
thermophysicalTransport.predict();
|
||||
}
|
||||
|
||||
#include "UEqn.H"
|
||||
#include "ftEqn.H"
|
||||
#include "bEqn.H"
|
||||
#include "EauEqn.H"
|
||||
#include "EaEqn.H"
|
||||
|
||||
if (!ign.ignited())
|
||||
{
|
||||
thermo.heu() == thermo.he();
|
||||
}
|
||||
|
||||
// --- Pressure corrector loop
|
||||
while (pimple.correct())
|
||||
{
|
||||
#include "pEqn.H"
|
||||
}
|
||||
|
||||
if (pimple.correctTransport())
|
||||
{
|
||||
turbulence->correct();
|
||||
thermophysicalTransport.correct();
|
||||
}
|
||||
}
|
||||
|
||||
rho = thermo.rho();
|
||||
|
||||
runTime.write();
|
||||
|
||||
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
|
||||
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
|
||||
<< nl << endl;
|
||||
}
|
||||
|
||||
Info<< "End\n" << endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,282 +0,0 @@
|
||||
if (ign.ignited())
|
||||
{
|
||||
// progress variable
|
||||
// ~~~~~~~~~~~~~~~~~
|
||||
volScalarField c("c", scalar(1) - b);
|
||||
|
||||
// Unburnt gas density
|
||||
// ~~~~~~~~~~~~~~~~~~~
|
||||
volScalarField rhou(thermo.rhou());
|
||||
|
||||
// Calculate flame normal etc.
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
volVectorField n("n", fvc::grad(b));
|
||||
|
||||
volScalarField mgb(mag(n));
|
||||
|
||||
dimensionedScalar dMgb = 1.0e-3*
|
||||
(b*c*mgb)().weightedAverage(mesh.V())
|
||||
/((b*c)().weightedAverage(mesh.V()) + small)
|
||||
+ dimensionedScalar(mgb.dimensions(), small);
|
||||
|
||||
mgb += dMgb;
|
||||
|
||||
surfaceVectorField SfHat(mesh.Sf()/mesh.magSf());
|
||||
surfaceVectorField nfVec(fvc::interpolate(n));
|
||||
nfVec += SfHat*(fvc::snGrad(b) - (SfHat & nfVec));
|
||||
nfVec /= (mag(nfVec) + dMgb);
|
||||
surfaceScalarField nf((mesh.Sf() & nfVec));
|
||||
n /= mgb;
|
||||
|
||||
|
||||
#include "StCorr.H"
|
||||
|
||||
// Calculate turbulent flame speed flux
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
surfaceScalarField phiSt("phiSt", fvc::interpolate(rhou*StCorr*Su*Xi)*nf);
|
||||
|
||||
scalar StCoNum = max
|
||||
(
|
||||
mesh.surfaceInterpolation::deltaCoeffs()
|
||||
*mag(phiSt)/(fvc::interpolate(rho)*mesh.magSf())
|
||||
).value()*runTime.deltaTValue();
|
||||
|
||||
Info<< "Max St-Courant Number = " << StCoNum << endl;
|
||||
|
||||
// Create b equation
|
||||
// ~~~~~~~~~~~~~~~~~
|
||||
fvScalarMatrix bEqn
|
||||
(
|
||||
fvm::ddt(rho, b)
|
||||
+ mvConvection->fvmDiv(phi, b)
|
||||
+ fvm::div(phiSt, b)
|
||||
- fvm::Sp(fvc::div(phiSt), b)
|
||||
- fvm::laplacian(thermophysicalTransport.DEff(b), b)
|
||||
==
|
||||
fvModels.source(rho, b)
|
||||
);
|
||||
|
||||
|
||||
// Add ignition cell contribution to b-equation
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
#include "ignite.H"
|
||||
|
||||
|
||||
// Solve for b
|
||||
// ~~~~~~~~~~~
|
||||
bEqn.relax();
|
||||
|
||||
fvConstraints.constrain(bEqn);
|
||||
|
||||
bEqn.solve();
|
||||
|
||||
fvConstraints.constrain(b);
|
||||
|
||||
Info<< "min(b) = " << min(b).value() << endl;
|
||||
|
||||
|
||||
// Calculate coefficients for Gulder's flame speed correlation
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
volScalarField up(uPrimeCoef*sqrt((2.0/3.0)*turbulence->k()));
|
||||
// volScalarField up(sqrt(mag(diag(n * n) & diag(turbulence->r()))));
|
||||
|
||||
volScalarField epsilon(pow(uPrimeCoef, 3)*turbulence->epsilon());
|
||||
|
||||
volScalarField tauEta(sqrt(thermo.muu()/(rhou*epsilon)));
|
||||
|
||||
volScalarField Reta
|
||||
(
|
||||
up
|
||||
/ (
|
||||
sqrt(epsilon*tauEta)
|
||||
+ dimensionedScalar(up.dimensions(), 1e-8)
|
||||
)
|
||||
);
|
||||
|
||||
// volScalarField l = 0.337*k*sqrt(k)/epsilon;
|
||||
// Reta *= max((l - dimensionedScalar(dimLength, 1.5e-3))/l, 0);
|
||||
|
||||
// Calculate Xi flux
|
||||
// ~~~~~~~~~~~~~~~~~
|
||||
surfaceScalarField phiXi
|
||||
(
|
||||
phiSt
|
||||
- fvc::interpolate
|
||||
(
|
||||
fvc::laplacian(thermophysicalTransport.DEff(b), b)/mgb
|
||||
)*nf
|
||||
+ fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
|
||||
);
|
||||
|
||||
|
||||
// Calculate mean and turbulent strain rates
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
volVectorField Ut(U + Su*Xi*n);
|
||||
volScalarField sigmat((n & n)*fvc::div(Ut) - (n & fvc::grad(Ut) & n));
|
||||
|
||||
volScalarField sigmas
|
||||
(
|
||||
((n & n)*fvc::div(U) - (n & fvc::grad(U) & n))/Xi
|
||||
+ (
|
||||
(n & n)*fvc::div(Su*n)
|
||||
- (n & fvc::grad(Su*n) & n)
|
||||
)*(Xi + scalar(1))/(2*Xi)
|
||||
);
|
||||
|
||||
|
||||
// Calculate the unstrained laminar flame speed
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
volScalarField Su0(unstrainedLaminarFlameSpeed()());
|
||||
|
||||
|
||||
// Calculate the laminar flame speed in equilibrium with the applied strain
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
volScalarField SuInf(Su0*max(scalar(1) - sigmas/sigmaExt, scalar(0.01)));
|
||||
|
||||
if (SuModel == "unstrained")
|
||||
{
|
||||
Su == Su0;
|
||||
}
|
||||
else if (SuModel == "equilibrium")
|
||||
{
|
||||
Su == SuInf;
|
||||
}
|
||||
else if (SuModel == "transport")
|
||||
{
|
||||
// Solve for the strained laminar flame speed
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
volScalarField Rc
|
||||
(
|
||||
(sigmas*SuInf*(Su0 - SuInf) + sqr(SuMin)*sigmaExt)
|
||||
/(sqr(Su0 - SuInf) + sqr(SuMin))
|
||||
);
|
||||
|
||||
fvScalarMatrix SuEqn
|
||||
(
|
||||
fvm::ddt(rho, Su)
|
||||
+ fvm::div(phi + phiXi, Su, "div(phiXi,Su)")
|
||||
- fvm::Sp(fvc::div(phiXi), Su)
|
||||
==
|
||||
- fvm::SuSp(-rho*Rc*Su0/Su, Su)
|
||||
- fvm::SuSp(rho*(sigmas + Rc), Su)
|
||||
+ fvModels.source(rho, Su)
|
||||
);
|
||||
|
||||
SuEqn.relax();
|
||||
|
||||
fvConstraints.constrain(SuEqn);
|
||||
|
||||
SuEqn.solve();
|
||||
|
||||
fvConstraints.constrain(Su);
|
||||
|
||||
// Limit the maximum Su
|
||||
// ~~~~~~~~~~~~~~~~~~~~
|
||||
Su.min(SuMax);
|
||||
Su.max(SuMin);
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalError
|
||||
<< args.executable() << " : Unknown Su model " << SuModel
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
|
||||
// Calculate Xi according to the selected flame wrinkling model
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
if (XiModel == "fixed")
|
||||
{
|
||||
// Do nothing, Xi is fixed!
|
||||
}
|
||||
else if (XiModel == "algebraic")
|
||||
{
|
||||
// Simple algebraic model for Xi based on Gulders correlation
|
||||
// with a linear correction function to give a plausible profile for Xi
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
Xi == scalar(1) +
|
||||
(scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
|
||||
*XiCoef*sqrt(up/(Su + SuMin))*Reta;
|
||||
}
|
||||
else if (XiModel == "transport")
|
||||
{
|
||||
// Calculate Xi transport coefficients based on Gulders correlation
|
||||
// and DNS data for the rate of generation
|
||||
// with a linear correction function to give a plausible profile for Xi
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
volScalarField XiEqStar
|
||||
(
|
||||
scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta
|
||||
);
|
||||
|
||||
volScalarField XiEq
|
||||
(
|
||||
scalar(1.001)
|
||||
+ (
|
||||
scalar(1)
|
||||
+ (2*XiShapeCoef)
|
||||
*(scalar(0.5) - min(max(b, scalar(0)), scalar(1)))
|
||||
)*(XiEqStar - scalar(1.001))
|
||||
);
|
||||
|
||||
volScalarField Gstar(0.28/tauEta);
|
||||
volScalarField R(Gstar*XiEqStar/(XiEqStar - scalar(1)));
|
||||
volScalarField G(R*(XiEq - scalar(1.001))/XiEq);
|
||||
|
||||
// R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;
|
||||
|
||||
// Solve for the flame wrinkling
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
fvScalarMatrix XiEqn
|
||||
(
|
||||
fvm::ddt(rho, Xi)
|
||||
+ fvm::div(phi + phiXi, Xi, "div(phiXi,Xi)")
|
||||
- fvm::Sp(fvc::div(phiXi), Xi)
|
||||
==
|
||||
rho*R
|
||||
- fvm::Sp(rho*(R - G), Xi)
|
||||
- fvm::Sp
|
||||
(
|
||||
rho*max
|
||||
(
|
||||
sigmat - sigmas,
|
||||
dimensionedScalar(sigmat.dimensions(), 0)
|
||||
),
|
||||
Xi
|
||||
)
|
||||
+ fvModels.source(rho, Xi)
|
||||
);
|
||||
|
||||
XiEqn.relax();
|
||||
|
||||
fvConstraints.constrain(XiEqn);
|
||||
|
||||
XiEqn.solve();
|
||||
|
||||
fvConstraints.constrain(Xi);
|
||||
|
||||
// Correct boundedness of Xi
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
Xi.max(1.0);
|
||||
Info<< "max(Xi) = " << max(Xi).value() << endl;
|
||||
Info<< "max(XiEq) = " << max(XiEq).value() << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalError
|
||||
<< args.executable() << " : Unknown Xi model " << XiModel
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
Info<< "Combustion progress = "
|
||||
<< 100*(scalar(1) - b)().weightedAverage(mesh.V()).value() << "%"
|
||||
<< endl;
|
||||
|
||||
St = Xi*Su;
|
||||
}
|
||||
@ -1,19 +0,0 @@
|
||||
// Calculate absolute flux
|
||||
// from the mapped surface velocity
|
||||
phi = mesh.Sf() & rhoUf();
|
||||
|
||||
correctUphiBCs(rho, U, phi, true);
|
||||
|
||||
CorrectPhi
|
||||
(
|
||||
phi,
|
||||
p,
|
||||
rho,
|
||||
psi,
|
||||
dimensionedScalar("rAUf", dimTime, 1),
|
||||
divrhoU(),
|
||||
pimple
|
||||
);
|
||||
|
||||
// Make the fluxes relative to the mesh-motion
|
||||
fvc::makeRelative(phi, rho, U);
|
||||
@ -1 +0,0 @@
|
||||
const volScalarField& psi = thermo.psi();
|
||||
@ -1,157 +0,0 @@
|
||||
Info<< "Reading thermophysical properties\n" << endl;
|
||||
|
||||
autoPtr<psiuMulticomponentThermo> pThermo
|
||||
(
|
||||
psiuMulticomponentThermo::New(mesh)
|
||||
);
|
||||
psiuMulticomponentThermo& thermo = pThermo();
|
||||
thermo.validate(args.executable(), "ha", "ea");
|
||||
|
||||
basicCombustionMixture& composition = thermo.composition();
|
||||
|
||||
volScalarField rho
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"rho",
|
||||
runTime.name(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
thermo.renameRho()
|
||||
);
|
||||
|
||||
volScalarField& p = thermo.p();
|
||||
|
||||
volScalarField& b = composition.Y("b");
|
||||
Info<< "min(b) = " << min(b).value() << endl;
|
||||
|
||||
|
||||
Info<< "\nReading field U\n" << endl;
|
||||
volVectorField U
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"U",
|
||||
runTime.name(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
#include "compressibleCreatePhi.H"
|
||||
|
||||
mesh.schemes().setFluxRequired(p.name());
|
||||
|
||||
Info<< "Creating turbulence model\n" << endl;
|
||||
autoPtr<compressible::momentumTransportModel> turbulence
|
||||
(
|
||||
compressible::momentumTransportModel::New
|
||||
(
|
||||
rho,
|
||||
U,
|
||||
phi,
|
||||
thermo
|
||||
)
|
||||
);
|
||||
|
||||
Info<< "Creating thermophysical transport model\n" << endl;
|
||||
turbulenceThermophysicalTransportModels::unityLewisEddyDiffusivity
|
||||
<
|
||||
RASThermophysicalTransportModel
|
||||
<
|
||||
ThermophysicalTransportModel
|
||||
<
|
||||
compressibleMomentumTransportModel,
|
||||
fluidThermo
|
||||
>
|
||||
>
|
||||
> thermophysicalTransport(turbulence(), thermo, true);
|
||||
|
||||
Info<< "Creating field dpdt\n" << endl;
|
||||
volScalarField dpdt
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"dpdt",
|
||||
runTime.name(),
|
||||
mesh
|
||||
),
|
||||
mesh,
|
||||
dimensionedScalar(p.dimensions()/dimTime, 0)
|
||||
);
|
||||
|
||||
Info<< "Creating field kinetic energy K\n" << endl;
|
||||
volScalarField K("K", 0.5*magSqr(U));
|
||||
|
||||
Info<< "Creating field Xi\n" << endl;
|
||||
volScalarField Xi
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"Xi",
|
||||
runTime.name(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
|
||||
Info<< "Creating the unstrained laminar flame speed\n" << endl;
|
||||
autoPtr<laminarFlameSpeed> unstrainedLaminarFlameSpeed
|
||||
(
|
||||
laminarFlameSpeed::New(thermo)
|
||||
);
|
||||
|
||||
|
||||
Info<< "Reading strained laminar flame speed field Su\n" << endl;
|
||||
volScalarField Su
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"Su",
|
||||
runTime.name(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
dimensionedScalar SuMin = 0.01*Su.average();
|
||||
dimensionedScalar SuMax = 4*Su.average();
|
||||
|
||||
Info<< "Calculating turbulent flame speed field St\n" << endl;
|
||||
volScalarField St
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"St",
|
||||
runTime.name(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
Xi*Su
|
||||
);
|
||||
|
||||
|
||||
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
|
||||
|
||||
if (composition.contains("ft"))
|
||||
{
|
||||
fields.add(composition.Y("ft"));
|
||||
}
|
||||
|
||||
fields.add(b);
|
||||
fields.add(thermo.he());
|
||||
fields.add(thermo.heu());
|
||||
|
||||
#include "createMRF.H"
|
||||
#include "createFvModels.H"
|
||||
#include "createFvConstraints.H"
|
||||
@ -1,30 +0,0 @@
|
||||
tmp<fv::convectionScheme<scalar>> mvConvection
|
||||
(
|
||||
fv::convectionScheme<scalar>::New
|
||||
(
|
||||
mesh,
|
||||
fields,
|
||||
phi,
|
||||
mesh.schemes().div("div(phi,ft_b_ha_hau)")
|
||||
)
|
||||
);
|
||||
|
||||
if (composition.contains("ft"))
|
||||
{
|
||||
volScalarField& ft = composition.Y("ft");
|
||||
|
||||
fvScalarMatrix ftEqn
|
||||
(
|
||||
fvm::ddt(rho, ft)
|
||||
+ mvConvection->fvmDiv(phi, ft)
|
||||
- fvm::laplacian(thermophysicalTransport.DEff(ft), ft)
|
||||
==
|
||||
fvModels.source(rho, ft)
|
||||
);
|
||||
|
||||
fvConstraints.constrain(ftEqn);
|
||||
|
||||
ftEqn.solve();
|
||||
|
||||
fvConstraints.constrain(ft);
|
||||
}
|
||||
@ -1,11 +0,0 @@
|
||||
EXE_INC = \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
-I$(LIB_SRC)/dynamicMesh/lnInclude \
|
||||
-I$(LIB_SRC)/fvMotionSolver/lnInclude
|
||||
|
||||
LIB_LIBS = \
|
||||
-lfiniteVolume \
|
||||
-lmeshTools \
|
||||
-ldynamicMesh \
|
||||
-lfvMotionSolvers
|
||||
@ -1,89 +0,0 @@
|
||||
rho = thermo.rho();
|
||||
|
||||
volScalarField rAU(1.0/UEqn.A());
|
||||
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
|
||||
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
|
||||
|
||||
surfaceScalarField phiHbyA
|
||||
(
|
||||
"phiHbyA",
|
||||
fvc::interpolate(rho)*fvc::flux(HbyA)
|
||||
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf))
|
||||
);
|
||||
|
||||
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
|
||||
|
||||
// Update the pressure BCs to ensure flux consistency
|
||||
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
|
||||
|
||||
fvc::makeRelative(phiHbyA, rho, U);
|
||||
|
||||
if (pimple.transonic())
|
||||
{
|
||||
const surfaceScalarField phid
|
||||
(
|
||||
"phid",
|
||||
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
|
||||
);
|
||||
|
||||
while (pimple.correctNonOrthogonal())
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::ddt(psi, p)
|
||||
+ fvm::div(phid, p)
|
||||
- fvm::laplacian(rhorAUf, p)
|
||||
==
|
||||
fvModels.source(psi, p, rho.name())
|
||||
);
|
||||
|
||||
pEqn.solve();
|
||||
|
||||
if (pimple.finalNonOrthogonalIter())
|
||||
{
|
||||
phi == pEqn.flux();
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
while (pimple.correctNonOrthogonal())
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::ddt(psi, p)
|
||||
+ fvc::div(phiHbyA)
|
||||
- fvm::laplacian(rhorAUf, p)
|
||||
==
|
||||
fvModels.source(psi, p, rho.name())
|
||||
);
|
||||
|
||||
pEqn.solve();
|
||||
|
||||
if (pimple.finalNonOrthogonalIter())
|
||||
{
|
||||
phi = phiHbyA + pEqn.flux();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#include "rhoEqn.H"
|
||||
#include "compressibleContinuityErrs.H"
|
||||
|
||||
U = HbyA - rAU*fvc::grad(p);
|
||||
U.correctBoundaryConditions();
|
||||
fvConstraints.constrain(U);
|
||||
K = 0.5*magSqr(U);
|
||||
|
||||
// Correct rhoUf if the mesh is moving
|
||||
fvc::correctRhoUf(rhoUf, rho, U, phi, MRF);
|
||||
|
||||
if (thermo.dpdt())
|
||||
{
|
||||
dpdt = fvc::ddt(p);
|
||||
|
||||
if (mesh.moving())
|
||||
{
|
||||
dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
|
||||
}
|
||||
}
|
||||
@ -1,45 +0,0 @@
|
||||
Info<< "Reading combustion properties\n" << endl;
|
||||
|
||||
IOdictionary combustionProperties
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"combustionProperties",
|
||||
runTime.constant(),
|
||||
mesh,
|
||||
IOobject::MUST_READ_IF_MODIFIED,
|
||||
IOobject::NO_WRITE
|
||||
)
|
||||
);
|
||||
|
||||
word SuModel
|
||||
(
|
||||
combustionProperties.lookup("SuModel")
|
||||
);
|
||||
|
||||
dimensionedScalar sigmaExt
|
||||
(
|
||||
combustionProperties.lookup("sigmaExt")
|
||||
);
|
||||
|
||||
word XiModel
|
||||
(
|
||||
combustionProperties.lookup("XiModel")
|
||||
);
|
||||
|
||||
dimensionedScalar XiCoef
|
||||
(
|
||||
combustionProperties.lookup("XiCoef")
|
||||
);
|
||||
|
||||
dimensionedScalar XiShapeCoef
|
||||
(
|
||||
combustionProperties.lookup("XiShapeCoef")
|
||||
);
|
||||
|
||||
dimensionedScalar uPrimeCoef
|
||||
(
|
||||
combustionProperties.lookup("uPrimeCoef")
|
||||
);
|
||||
|
||||
ignition ign(combustionProperties, runTime, mesh);
|
||||
@ -9,6 +9,7 @@ wmake $targetType incompressibleFluid
|
||||
wmake $targetType isothermalFluid
|
||||
wmake $targetType fluid
|
||||
wmake $targetType multicomponentFluid
|
||||
wmake $targetType XiFluid
|
||||
wmake $targetType VoFSolver
|
||||
incompressibleVoF/Allwmake $targetType $*
|
||||
compressibleVoF/Allwmake $targetType $*
|
||||
|
||||
4
applications/solvers/modules/XiFluid/Make/files
Normal file
4
applications/solvers/modules/XiFluid/Make/files
Normal file
@ -0,0 +1,4 @@
|
||||
XiFluid.C
|
||||
thermophysicalPredictor.C
|
||||
|
||||
LIB = $(FOAM_LIBBIN)/libXiFluid
|
||||
@ -1,21 +1,25 @@
|
||||
EXE_INC = \
|
||||
-Iignition/lnInclude \
|
||||
-I$(LIB_SRC)/physicalProperties/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \
|
||||
-I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
|
||||
-I$(FOAM_SOLVERS)/modules/isothermalFluid/lnInclude \
|
||||
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
|
||||
-I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
|
||||
-I$(LIB_SRC)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \
|
||||
-I$(LIB_SRC)/ThermophysicalTransportModels/fluid/lnInclude \
|
||||
-I$(LIB_SRC)/ThermophysicalTransportModels/fluidThermo/lnInclude \
|
||||
-I$(LIB_SRC)/physicalProperties/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/ignition/lnInclude \
|
||||
-I$(LIB_SRC)/finiteVolume/cfdTools \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||
-I$(LIB_SRC)/sampling/lnInclude \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
-I$(LIB_SRC)/sampling/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lXiIgnition \
|
||||
LIB_LIBS = \
|
||||
-lfluidSolver \
|
||||
-lisothermalFluid \
|
||||
-lmomentumTransportModels \
|
||||
-lcompressibleMomentumTransportModels \
|
||||
-lfluidThermoThermophysicalTransportModels \
|
||||
@ -23,6 +27,7 @@ EXE_LIBS = \
|
||||
-lmulticomponentThermophysicalModels \
|
||||
-lspecie \
|
||||
-llaminarFlameSpeedModels \
|
||||
-lXiIgnition \
|
||||
-lfiniteVolume \
|
||||
-lfvModels \
|
||||
-lfvConstraints \
|
||||
199
applications/solvers/modules/XiFluid/XiFluid.C
Normal file
199
applications/solvers/modules/XiFluid/XiFluid.C
Normal file
@ -0,0 +1,199 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
|
||||
\\/ 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 "XiFluid.H"
|
||||
#include "localEulerDdtScheme.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
namespace solvers
|
||||
{
|
||||
defineTypeNameAndDebug(XiFluid, 0);
|
||||
addToRunTimeSelectionTable(solver, XiFluid, fvMesh);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::solvers::XiFluid::XiFluid(fvMesh& mesh)
|
||||
:
|
||||
isothermalFluid
|
||||
(
|
||||
mesh,
|
||||
autoPtr<fluidThermo>(psiuMulticomponentThermo::New(mesh).ptr())
|
||||
),
|
||||
|
||||
thermo(refCast<psiuMulticomponentThermo>(isothermalFluid::thermo)),
|
||||
|
||||
composition(thermo.composition()),
|
||||
|
||||
b(composition.Y("b")),
|
||||
|
||||
unstrainedLaminarFlameSpeed
|
||||
(
|
||||
laminarFlameSpeed::New(thermo)
|
||||
),
|
||||
|
||||
Su
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"Su",
|
||||
runTime.name(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
),
|
||||
|
||||
SuMin(0.01*Su.average()),
|
||||
SuMax(4*Su.average()),
|
||||
|
||||
Xi
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"Xi",
|
||||
runTime.name(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
),
|
||||
|
||||
St
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"St",
|
||||
runTime.name(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
Xi*Su
|
||||
),
|
||||
|
||||
combustionProperties
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"combustionProperties",
|
||||
runTime.constant(),
|
||||
mesh,
|
||||
IOobject::MUST_READ_IF_MODIFIED,
|
||||
IOobject::NO_WRITE
|
||||
)
|
||||
),
|
||||
|
||||
SuModel
|
||||
(
|
||||
combustionProperties.lookup("SuModel")
|
||||
),
|
||||
|
||||
sigmaExt
|
||||
(
|
||||
combustionProperties.lookup("sigmaExt")
|
||||
),
|
||||
|
||||
XiModel
|
||||
(
|
||||
combustionProperties.lookup("XiModel")
|
||||
),
|
||||
|
||||
XiCoef
|
||||
(
|
||||
combustionProperties.lookup("XiCoef")
|
||||
),
|
||||
|
||||
XiShapeCoef
|
||||
(
|
||||
combustionProperties.lookup("XiShapeCoef")
|
||||
),
|
||||
|
||||
uPrimeCoef
|
||||
(
|
||||
combustionProperties.lookup("uPrimeCoef")
|
||||
),
|
||||
|
||||
ign(combustionProperties, runTime, mesh),
|
||||
|
||||
thermophysicalTransport
|
||||
(
|
||||
momentumTransport(),
|
||||
thermo,
|
||||
true
|
||||
)
|
||||
{
|
||||
thermo.validate(type(), "ha", "ea");
|
||||
|
||||
if (composition.contains("ft"))
|
||||
{
|
||||
fields.add(composition.Y("ft"));
|
||||
}
|
||||
|
||||
fields.add(b);
|
||||
fields.add(thermo.he());
|
||||
fields.add(thermo.heu());
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::solvers::XiFluid::~XiFluid()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
void Foam::solvers::XiFluid::prePredictor()
|
||||
{
|
||||
isothermalFluid::prePredictor();
|
||||
|
||||
if (pimple.predictTransport())
|
||||
{
|
||||
thermophysicalTransport.predict();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::solvers::XiFluid::postCorrector()
|
||||
{
|
||||
isothermalFluid::postCorrector();
|
||||
|
||||
if (pimple.correctTransport())
|
||||
{
|
||||
thermophysicalTransport.correct();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
264
applications/solvers/modules/XiFluid/XiFluid.H
Normal file
264
applications/solvers/modules/XiFluid/XiFluid.H
Normal file
@ -0,0 +1,264 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
|
||||
\\/ 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/>.
|
||||
|
||||
Class
|
||||
Foam::solvers::XiFluid
|
||||
|
||||
Description
|
||||
Solver module for compressible premixed/partially-premixed combustion with
|
||||
turbulence modelling.
|
||||
|
||||
Combusting RANS code using the b-Xi two-equation model.
|
||||
Xi may be obtained by either the solution of the Xi transport
|
||||
equation or from an algebraic expression. Both approaches are
|
||||
based on Gulder's flame speed correlation which has been shown
|
||||
to be appropriate by comparison with the results from the
|
||||
spectral model.
|
||||
|
||||
Strain effects are encorporated directly into the Xi equation
|
||||
but not in the algebraic approximation. Further work need to be
|
||||
done on this issue, particularly regarding the enhanced removal rate
|
||||
caused by flame compression. Analysis using results of the spectral
|
||||
model will be required.
|
||||
|
||||
For cases involving very lean Propane flames or other flames which are
|
||||
very strain-sensitive, a transport equation for the laminar flame
|
||||
speed is present. This equation is derived using heuristic arguments
|
||||
involving the strain time scale and the strain-rate at extinction.
|
||||
the transport velocity is the same as that for the Xi equation.
|
||||
|
||||
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
|
||||
pseudo-transient and steady simulations.
|
||||
|
||||
Optional fvModels and fvConstraints are provided to enhance the simulation
|
||||
in many ways including adding various sources, chemical reactions,
|
||||
combustion, Lagrangian particles, radiation, surface film etc. and
|
||||
constraining or limiting the solution.
|
||||
|
||||
Reference:
|
||||
\verbatim
|
||||
Greenshields, C. J., & Weller, H. G. (2022).
|
||||
Notes on Computational Fluid Dynamics: General Principles.
|
||||
CFD Direct Ltd.: Reading, UK.
|
||||
\endverbatim
|
||||
|
||||
SourceFiles
|
||||
XiFluid.C
|
||||
|
||||
See also
|
||||
Foam::solvers::fluidSolver
|
||||
Foam::solvers::isothermalFluid
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef XiFluid_H
|
||||
#define XiFluid_H
|
||||
|
||||
#include "isothermalFluid.H"
|
||||
#include "psiuMulticomponentThermo.H"
|
||||
#include "RASThermophysicalTransportModel.H"
|
||||
#include "unityLewisEddyDiffusivity.H"
|
||||
#include "laminarFlameSpeed.H"
|
||||
#include "ignition.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
namespace solvers
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class XiFluid Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class XiFluid
|
||||
:
|
||||
public isothermalFluid
|
||||
{
|
||||
|
||||
protected:
|
||||
|
||||
// Thermophysical properties
|
||||
|
||||
psiuMulticomponentThermo& thermo;
|
||||
|
||||
|
||||
// Composition
|
||||
|
||||
//- Reference to the combustion mixture
|
||||
basicCombustionMixture& composition;
|
||||
|
||||
//- Reference to the combustion regress variable
|
||||
// obtained from the combustion mixture
|
||||
volScalarField& b;
|
||||
|
||||
//- Set of fields used for the multivariate convection scheme
|
||||
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
|
||||
|
||||
|
||||
// Reactions
|
||||
|
||||
//- Laminar flame-speed model
|
||||
autoPtr<laminarFlameSpeed> unstrainedLaminarFlameSpeed;
|
||||
|
||||
//- Laminar flame-speed field
|
||||
volScalarField Su;
|
||||
|
||||
//- Minimum laminar flame-speed allowed for numerical stability
|
||||
dimensionedScalar SuMin;
|
||||
|
||||
//- Maximum laminar flame-speed allowed for numerical stability
|
||||
dimensionedScalar SuMax;
|
||||
|
||||
//- Flame wrinkling coefficient field
|
||||
volScalarField Xi;
|
||||
|
||||
//- Turbulent flame-speed field
|
||||
volScalarField St;
|
||||
|
||||
//- Dictionary of combustion model coefficients
|
||||
IOdictionary combustionProperties;
|
||||
|
||||
//- Name of the strained laminar flame-speed model
|
||||
word SuModel;
|
||||
|
||||
//- Laminar flame extinction strain-rate
|
||||
dimensionedScalar sigmaExt;
|
||||
|
||||
//- Name of the flame wrinkling model
|
||||
word XiModel;
|
||||
|
||||
//- Flame wrinkling model coefficient
|
||||
dimensionedScalar XiCoef;
|
||||
|
||||
//- Flame wrinkling model shape coefficient
|
||||
dimensionedScalar XiShapeCoef;
|
||||
|
||||
//- Flame wrinkling model u' coefficient
|
||||
dimensionedScalar uPrimeCoef;
|
||||
|
||||
//- Ignition model
|
||||
ignition ign;
|
||||
|
||||
|
||||
// Thermophysical transport
|
||||
|
||||
turbulenceThermophysicalTransportModels::unityLewisEddyDiffusivity
|
||||
<
|
||||
RASThermophysicalTransportModel
|
||||
<
|
||||
ThermophysicalTransportModel
|
||||
<
|
||||
compressibleMomentumTransportModel,
|
||||
fluidThermo
|
||||
>
|
||||
>
|
||||
> thermophysicalTransport;
|
||||
|
||||
|
||||
// Protected member functions
|
||||
|
||||
//- Solve the ft equation for partially-premixed mixtures
|
||||
void ftSolve
|
||||
(
|
||||
const fv::convectionScheme<scalar>& mvConvection
|
||||
);
|
||||
|
||||
//- Calculate and return the turbulent flame-speed kernel correction
|
||||
dimensionedScalar StCorr
|
||||
(
|
||||
const volScalarField& c,
|
||||
const surfaceScalarField& nf,
|
||||
const dimensionedScalar& dMgb
|
||||
) const;
|
||||
|
||||
//- Solve the Xi and regress variable equations
|
||||
void bSolve
|
||||
(
|
||||
const fv::convectionScheme<scalar>& mvConvection
|
||||
);
|
||||
|
||||
//- Solve the unburnt energy equation
|
||||
void EauSolve
|
||||
(
|
||||
const fv::convectionScheme<scalar>& mvConvection
|
||||
);
|
||||
|
||||
//- Solve the energy equation
|
||||
void EaSolve
|
||||
(
|
||||
const fv::convectionScheme<scalar>& mvConvection
|
||||
);
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("XiFluid");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from region mesh
|
||||
XiFluid(fvMesh& mesh);
|
||||
|
||||
//- Disallow default bitwise copy construction
|
||||
XiFluid(const XiFluid&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~XiFluid();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Called at the start of the PIMPLE loop
|
||||
virtual void prePredictor();
|
||||
|
||||
//- Construct and solve the energy equation,
|
||||
// convert to temperature
|
||||
// and update thermophysical and transport properties
|
||||
virtual void thermophysicalPredictor();
|
||||
|
||||
//- Correct the momentum and thermophysical transport modelling
|
||||
virtual void postCorrector();
|
||||
|
||||
|
||||
// Member Operators
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
void operator=(const XiFluid&) = delete;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace solvers
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
591
applications/solvers/modules/XiFluid/thermophysicalPredictor.C
Normal file
591
applications/solvers/modules/XiFluid/thermophysicalPredictor.C
Normal file
@ -0,0 +1,591 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration | Website: https://openfoam.org
|
||||
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
|
||||
\\/ 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 "XiFluid.H"
|
||||
#include "fvcSnGrad.H"
|
||||
#include "fvcLaplacian.H"
|
||||
#include "fvcDdt.H"
|
||||
#include "fvcMeshPhi.H"
|
||||
#include "fvmDiv.H"
|
||||
#include "fvmSup.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
void Foam::solvers::XiFluid::ftSolve
|
||||
(
|
||||
const fv::convectionScheme<scalar>& mvConvection
|
||||
)
|
||||
{
|
||||
volScalarField& ft = composition.Y("ft");
|
||||
|
||||
fvScalarMatrix ftEqn
|
||||
(
|
||||
fvm::ddt(rho, ft)
|
||||
+ mvConvection.fvmDiv(phi, ft)
|
||||
- fvm::laplacian(thermophysicalTransport.DEff(ft), ft)
|
||||
==
|
||||
fvModels().source(rho, ft)
|
||||
);
|
||||
|
||||
fvConstraints().constrain(ftEqn);
|
||||
|
||||
ftEqn.solve();
|
||||
|
||||
fvConstraints().constrain(ft);
|
||||
}
|
||||
|
||||
|
||||
Foam::dimensionedScalar Foam::solvers::XiFluid::StCorr
|
||||
(
|
||||
const volScalarField& c,
|
||||
const surfaceScalarField& nf,
|
||||
const dimensionedScalar& dMgb
|
||||
) const
|
||||
{
|
||||
dimensionedScalar StCorr("StCorr", dimless, 1.0);
|
||||
|
||||
if (ign.igniting())
|
||||
{
|
||||
// Calculate volume of ignition kernel
|
||||
const dimensionedScalar Vk("Vk", dimVolume, gSum(c*mesh.V().field()));
|
||||
dimensionedScalar Ak("Ak", dimArea, 0.0);
|
||||
|
||||
if (Vk.value() > small)
|
||||
{
|
||||
// Calculate kernel area from its volume
|
||||
// and the dimensionality of the case
|
||||
|
||||
switch(mesh.nGeometricD())
|
||||
{
|
||||
case 3:
|
||||
{
|
||||
// Assume it is part-spherical
|
||||
const scalar sphereFraction
|
||||
(
|
||||
combustionProperties.lookup<scalar>
|
||||
(
|
||||
"ignitionSphereFraction"
|
||||
)
|
||||
);
|
||||
|
||||
Ak = sphereFraction*4.0*constant::mathematical::pi
|
||||
*pow
|
||||
(
|
||||
3.0*Vk
|
||||
/(sphereFraction*4.0*constant::mathematical::pi),
|
||||
2.0/3.0
|
||||
);
|
||||
}
|
||||
break;
|
||||
|
||||
case 2:
|
||||
{
|
||||
// Assume it is part-circular
|
||||
const dimensionedScalar thickness
|
||||
(
|
||||
combustionProperties.lookup("ignitionThickness")
|
||||
);
|
||||
|
||||
const scalar circleFraction
|
||||
(
|
||||
combustionProperties.lookup<scalar>
|
||||
(
|
||||
"ignitionCircleFraction"
|
||||
)
|
||||
);
|
||||
|
||||
Ak = circleFraction*constant::mathematical::pi*thickness
|
||||
*sqrt
|
||||
(
|
||||
4.0*Vk
|
||||
/(
|
||||
circleFraction
|
||||
*thickness
|
||||
*constant::mathematical::pi
|
||||
)
|
||||
);
|
||||
}
|
||||
break;
|
||||
|
||||
case 1:
|
||||
// Assume it is plane or two planes
|
||||
Ak = dimensionedScalar
|
||||
(
|
||||
combustionProperties.lookup("ignitionKernelArea")
|
||||
);
|
||||
break;
|
||||
}
|
||||
|
||||
// Calculate kernel area from b field consistent with the
|
||||
// discretisation of the b equation.
|
||||
const volScalarField mgb
|
||||
(
|
||||
fvc::div(nf, b, "div(phiSt,b)") - b*fvc::div(nf) + dMgb
|
||||
);
|
||||
const dimensionedScalar AkEst = gSum(mgb*mesh.V().field());
|
||||
|
||||
StCorr.value() = max(min((Ak/AkEst).value(), 10.0), 1.0);
|
||||
|
||||
Info<< "StCorr = " << StCorr.value() << endl;
|
||||
}
|
||||
}
|
||||
|
||||
return StCorr;
|
||||
}
|
||||
|
||||
|
||||
void Foam::solvers::XiFluid::bSolve
|
||||
(
|
||||
const fv::convectionScheme<scalar>& mvConvection
|
||||
)
|
||||
{
|
||||
// progress variable
|
||||
// ~~~~~~~~~~~~~~~~~
|
||||
const volScalarField c("c", scalar(1) - b);
|
||||
|
||||
// Unburnt gas density
|
||||
// ~~~~~~~~~~~~~~~~~~~
|
||||
const volScalarField rhou(thermo.rhou());
|
||||
|
||||
// Calculate flame normal etc.
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
volVectorField n("n", fvc::grad(b));
|
||||
|
||||
volScalarField mgb(mag(n));
|
||||
|
||||
const dimensionedScalar dMgb = 1.0e-3*
|
||||
(b*c*mgb)().weightedAverage(mesh.V())
|
||||
/((b*c)().weightedAverage(mesh.V()) + small)
|
||||
+ dimensionedScalar(mgb.dimensions(), small);
|
||||
|
||||
mgb += dMgb;
|
||||
|
||||
const surfaceVectorField SfHat(mesh.Sf()/mesh.magSf());
|
||||
surfaceVectorField nfVec(fvc::interpolate(n));
|
||||
nfVec += SfHat*(fvc::snGrad(b) - (SfHat & nfVec));
|
||||
nfVec /= (mag(nfVec) + dMgb);
|
||||
surfaceScalarField nf((mesh.Sf() & nfVec));
|
||||
n /= mgb;
|
||||
|
||||
// Calculate turbulent flame speed flux
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
const surfaceScalarField phiSt
|
||||
(
|
||||
"phiSt",
|
||||
fvc::interpolate(rhou*StCorr(c, nf, dMgb)*Su*Xi)*nf
|
||||
);
|
||||
|
||||
const scalar StCoNum = max
|
||||
(
|
||||
mesh.surfaceInterpolation::deltaCoeffs()
|
||||
*mag(phiSt)/(fvc::interpolate(rho)*mesh.magSf())
|
||||
).value()*runTime.deltaTValue();
|
||||
|
||||
Info<< "Max St-Courant Number = " << StCoNum << endl;
|
||||
|
||||
// Create b equation
|
||||
// ~~~~~~~~~~~~~~~~~
|
||||
fvScalarMatrix bEqn
|
||||
(
|
||||
fvm::ddt(rho, b)
|
||||
+ mvConvection.fvmDiv(phi, b)
|
||||
+ fvm::div(phiSt, b)
|
||||
- fvm::Sp(fvc::div(phiSt), b)
|
||||
- fvm::laplacian(thermophysicalTransport.DEff(b), b)
|
||||
==
|
||||
fvModels().source(rho, b)
|
||||
);
|
||||
|
||||
|
||||
// Add ignition cell contribution to b-equation
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
forAll(ign.sites(), i)
|
||||
{
|
||||
const ignitionSite& ignSite = ign.sites()[i];
|
||||
|
||||
if (ignSite.igniting())
|
||||
{
|
||||
forAll(ignSite.cells(), icelli)
|
||||
{
|
||||
label ignCell = ignSite.cells()[icelli];
|
||||
Info<< "Igniting cell " << ignCell;
|
||||
|
||||
Info<< " state :"
|
||||
<< ' ' << b[ignCell]
|
||||
<< ' ' << Xi[ignCell]
|
||||
<< ' ' << Su[ignCell]
|
||||
<< ' ' << mgb[ignCell]
|
||||
<< endl;
|
||||
|
||||
bEqn.diag()[ignSite.cells()[icelli]] +=
|
||||
(
|
||||
ignSite.strength()*ignSite.cellVolumes()[icelli]
|
||||
*rhou[ignSite.cells()[icelli]]/ignSite.duration()
|
||||
)/(b[ignSite.cells()[icelli]] + 0.001);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Solve for b
|
||||
// ~~~~~~~~~~~
|
||||
bEqn.relax();
|
||||
|
||||
fvConstraints().constrain(bEqn);
|
||||
|
||||
bEqn.solve();
|
||||
|
||||
fvConstraints().constrain(b);
|
||||
|
||||
Info<< "min(b) = " << min(b).value() << endl;
|
||||
|
||||
|
||||
// Calculate coefficients for Gulder's flame speed correlation
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
const volScalarField up(uPrimeCoef*sqrt((2.0/3.0)*momentumTransport->k()));
|
||||
// volScalarField up(sqrt(mag(diag(n * n) & diag(momentumTransport->r()))));
|
||||
|
||||
const volScalarField epsilon
|
||||
(
|
||||
pow(uPrimeCoef, 3)*momentumTransport->epsilon()
|
||||
);
|
||||
|
||||
const volScalarField tauEta(sqrt(thermo.muu()/(rhou*epsilon)));
|
||||
|
||||
const volScalarField Reta
|
||||
(
|
||||
up
|
||||
/ (
|
||||
sqrt(epsilon*tauEta)
|
||||
+ dimensionedScalar(up.dimensions(), 1e-8)
|
||||
)
|
||||
);
|
||||
|
||||
// volScalarField l = 0.337*k*sqrt(k)/epsilon;
|
||||
// Reta *= max((l - dimensionedScalar(dimLength, 1.5e-3))/l, 0);
|
||||
|
||||
// Calculate Xi flux
|
||||
// ~~~~~~~~~~~~~~~~~
|
||||
const surfaceScalarField phiXi
|
||||
(
|
||||
phiSt
|
||||
- fvc::interpolate
|
||||
(
|
||||
fvc::laplacian(thermophysicalTransport.DEff(b), b)/mgb
|
||||
)*nf
|
||||
+ fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
|
||||
);
|
||||
|
||||
|
||||
// Calculate mean and turbulent strain rates
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
const volVectorField Ut(U + Su*Xi*n);
|
||||
const volScalarField sigmat((n & n)*fvc::div(Ut) - (n & fvc::grad(Ut) & n));
|
||||
|
||||
const volScalarField sigmas
|
||||
(
|
||||
((n & n)*fvc::div(U) - (n & fvc::grad(U) & n))/Xi
|
||||
+ (
|
||||
(n & n)*fvc::div(Su*n)
|
||||
- (n & fvc::grad(Su*n) & n)
|
||||
)*(Xi + scalar(1))/(2*Xi)
|
||||
);
|
||||
|
||||
|
||||
// Calculate the unstrained laminar flame speed
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
const volScalarField Su0(unstrainedLaminarFlameSpeed()());
|
||||
|
||||
|
||||
// Calculate the laminar flame speed in equilibrium with the applied strain
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
const volScalarField SuInf
|
||||
(
|
||||
Su0*max(scalar(1) - sigmas/sigmaExt, scalar(0.01))
|
||||
);
|
||||
|
||||
if (SuModel == "unstrained")
|
||||
{
|
||||
Su == Su0;
|
||||
}
|
||||
else if (SuModel == "equilibrium")
|
||||
{
|
||||
Su == SuInf;
|
||||
}
|
||||
else if (SuModel == "transport")
|
||||
{
|
||||
// Solve for the strained laminar flame speed
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
const volScalarField Rc
|
||||
(
|
||||
(sigmas*SuInf*(Su0 - SuInf) + sqr(SuMin)*sigmaExt)
|
||||
/(sqr(Su0 - SuInf) + sqr(SuMin))
|
||||
);
|
||||
|
||||
fvScalarMatrix SuEqn
|
||||
(
|
||||
fvm::ddt(rho, Su)
|
||||
+ fvm::div(phi + phiXi, Su, "div(phiXi,Su)")
|
||||
- fvm::Sp(fvc::div(phiXi), Su)
|
||||
==
|
||||
- fvm::SuSp(-rho*Rc*Su0/Su, Su)
|
||||
- fvm::SuSp(rho*(sigmas + Rc), Su)
|
||||
+ fvModels().source(rho, Su)
|
||||
);
|
||||
|
||||
SuEqn.relax();
|
||||
|
||||
fvConstraints().constrain(SuEqn);
|
||||
|
||||
SuEqn.solve();
|
||||
|
||||
fvConstraints().constrain(Su);
|
||||
|
||||
// Limit the maximum Su
|
||||
// ~~~~~~~~~~~~~~~~~~~~
|
||||
Su.min(SuMax);
|
||||
Su.max(SuMin);
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Unknown Su model " << SuModel
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
|
||||
// Calculate Xi according to the selected flame wrinkling model
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
if (XiModel == "fixed")
|
||||
{
|
||||
// Do nothing, Xi is fixed!
|
||||
}
|
||||
else if (XiModel == "algebraic")
|
||||
{
|
||||
// Simple algebraic model for Xi based on Gulders correlation
|
||||
// with a linear correction function to give a plausible profile for Xi
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
Xi == scalar(1) +
|
||||
(scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
|
||||
*XiCoef*sqrt(up/(Su + SuMin))*Reta;
|
||||
}
|
||||
else if (XiModel == "transport")
|
||||
{
|
||||
// Calculate Xi transport coefficients based on Gulders correlation
|
||||
// and DNS data for the rate of generation
|
||||
// with a linear correction function to give a plausible profile for Xi
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
const volScalarField XiEqStar
|
||||
(
|
||||
scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta
|
||||
);
|
||||
|
||||
const volScalarField XiEq
|
||||
(
|
||||
scalar(1.001)
|
||||
+ (
|
||||
scalar(1)
|
||||
+ (2*XiShapeCoef)
|
||||
*(scalar(0.5) - min(max(b, scalar(0)), scalar(1)))
|
||||
)*(XiEqStar - scalar(1.001))
|
||||
);
|
||||
|
||||
const volScalarField Gstar(0.28/tauEta);
|
||||
const volScalarField R(Gstar*XiEqStar/(XiEqStar - scalar(1)));
|
||||
const volScalarField G(R*(XiEq - scalar(1.001))/XiEq);
|
||||
|
||||
// R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;
|
||||
|
||||
// Solve for the flame wrinkling
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
fvScalarMatrix XiEqn
|
||||
(
|
||||
fvm::ddt(rho, Xi)
|
||||
+ fvm::div(phi + phiXi, Xi, "div(phiXi,Xi)")
|
||||
- fvm::Sp(fvc::div(phiXi), Xi)
|
||||
==
|
||||
rho*R
|
||||
- fvm::Sp(rho*(R - G), Xi)
|
||||
- fvm::Sp
|
||||
(
|
||||
rho*max
|
||||
(
|
||||
sigmat - sigmas,
|
||||
dimensionedScalar(sigmat.dimensions(), 0)
|
||||
),
|
||||
Xi
|
||||
)
|
||||
+ fvModels().source(rho, Xi)
|
||||
);
|
||||
|
||||
XiEqn.relax();
|
||||
|
||||
fvConstraints().constrain(XiEqn);
|
||||
|
||||
XiEqn.solve();
|
||||
|
||||
fvConstraints().constrain(Xi);
|
||||
|
||||
// Correct boundedness of Xi
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
Xi.max(1.0);
|
||||
Info<< "max(Xi) = " << max(Xi).value() << endl;
|
||||
Info<< "max(XiEq) = " << max(XiEq).value() << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Unknown Xi model " << XiModel
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
Info<< "Combustion progress = "
|
||||
<< 100*(scalar(1) - b)().weightedAverage(mesh.V()).value() << "%"
|
||||
<< endl;
|
||||
|
||||
St = Xi*Su;
|
||||
}
|
||||
|
||||
|
||||
void Foam::solvers::XiFluid::EauSolve
|
||||
(
|
||||
const fv::convectionScheme<scalar>& mvConvection
|
||||
)
|
||||
{
|
||||
volScalarField& heau = thermo.heu();
|
||||
|
||||
const volScalarField::Internal rhoByRhou(rho()/thermo.rhou()());
|
||||
|
||||
fvScalarMatrix heauEqn
|
||||
(
|
||||
fvm::ddt(rho, heau) + mvConvection.fvmDiv(phi, heau)
|
||||
+ rhoByRhou
|
||||
*(
|
||||
(fvc::ddt(rho, K) + fvc::div(phi, K))()
|
||||
+ pressureWork
|
||||
(
|
||||
heau.name() == "eau"
|
||||
? mvConvection.fvcDiv(phi, p/rho)()
|
||||
: -dpdt
|
||||
)
|
||||
)
|
||||
+ thermophysicalTransport.divq(heau)
|
||||
|
||||
// These terms cannot be used in partially-premixed combustion due to
|
||||
// the resultant inconsistency between ft and heau transport.
|
||||
// A possible solution would be to solve for ftu as well as ft.
|
||||
//- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
|
||||
//+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
|
||||
|
||||
==
|
||||
fvModels().source(rho, heau)
|
||||
);
|
||||
|
||||
fvConstraints().constrain(heauEqn);
|
||||
|
||||
heauEqn.solve();
|
||||
|
||||
fvConstraints().constrain(heau);
|
||||
}
|
||||
|
||||
|
||||
void Foam::solvers::XiFluid::EaSolve
|
||||
(
|
||||
const fv::convectionScheme<scalar>& mvConvection
|
||||
)
|
||||
{
|
||||
volScalarField& hea = thermo.he();
|
||||
|
||||
fvScalarMatrix EaEqn
|
||||
(
|
||||
fvm::ddt(rho, hea) + mvConvection.fvmDiv(phi, hea)
|
||||
+ fvc::ddt(rho, K) + fvc::div(phi, K)
|
||||
+ pressureWork
|
||||
(
|
||||
hea.name() == "ea"
|
||||
? mvConvection.fvcDiv(phi, p/rho)()
|
||||
: -dpdt
|
||||
)
|
||||
+ thermophysicalTransport.divq(hea)
|
||||
==
|
||||
(
|
||||
buoyancy.valid()
|
||||
? fvModels().source(rho, hea) + rho*(U & buoyancy->g)
|
||||
: fvModels().source(rho, hea)
|
||||
)
|
||||
);
|
||||
|
||||
EaEqn.relax();
|
||||
|
||||
fvConstraints().constrain(EaEqn);
|
||||
|
||||
EaEqn.solve();
|
||||
|
||||
fvConstraints().constrain(hea);
|
||||
}
|
||||
|
||||
|
||||
void Foam::solvers::XiFluid::thermophysicalPredictor()
|
||||
{
|
||||
tmp<fv::convectionScheme<scalar>> mvConvection
|
||||
(
|
||||
fv::convectionScheme<scalar>::New
|
||||
(
|
||||
mesh,
|
||||
fields,
|
||||
phi,
|
||||
mesh.schemes().div("div(phi,ft_b_ha_hau)")
|
||||
)
|
||||
);
|
||||
|
||||
if (composition.contains("ft"))
|
||||
{
|
||||
ftSolve(mvConvection());
|
||||
}
|
||||
|
||||
if (ign.ignited())
|
||||
{
|
||||
bSolve(mvConvection());
|
||||
EauSolve(mvConvection());
|
||||
}
|
||||
|
||||
EaSolve(mvConvection());
|
||||
|
||||
if (!ign.ignited())
|
||||
{
|
||||
thermo.heu() == thermo.he();
|
||||
}
|
||||
|
||||
thermo.correct();
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -25,6 +25,7 @@ License
|
||||
|
||||
#include "compressibleInterPhaseThermophysicalTransportModel.H"
|
||||
#include "compressibleInterPhaseTransportModel.H"
|
||||
#include "fvcSnGrad.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -25,6 +25,7 @@ License
|
||||
|
||||
#include "compressibleVoF.H"
|
||||
#include "localEulerDdtScheme.H"
|
||||
#include "fvcDdt.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
@ -63,6 +63,7 @@ See also
|
||||
#include "compressibleInterPhaseThermophysicalTransportModel.H"
|
||||
#include "buoyancy.H"
|
||||
#include "pressureReference.H"
|
||||
#include "fvmSup.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -27,6 +27,12 @@ License
|
||||
#include "constrainHbyA.H"
|
||||
#include "constrainPressure.H"
|
||||
#include "adjustPhi.H"
|
||||
#include "fvcMeshPhi.H"
|
||||
#include "fvcFlux.H"
|
||||
#include "fvcDdt.H"
|
||||
#include "fvcSnGrad.H"
|
||||
#include "fvcReconstruct.H"
|
||||
#include "fvmDiv.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -24,6 +24,9 @@ License
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "compressibleVoF.H"
|
||||
#include "fvcMeshPhi.H"
|
||||
#include "fvcDdt.H"
|
||||
#include "fvmDiv.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -52,7 +52,9 @@ Foam::solvers::fluid::fluid(fvMesh& mesh)
|
||||
thermo
|
||||
)
|
||||
)
|
||||
{}
|
||||
{
|
||||
thermo.validate(type(), "h", "e");
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
@ -24,6 +24,8 @@ License
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fluid.H"
|
||||
#include "fvcDdt.H"
|
||||
#include "fvmDiv.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
@ -38,7 +40,7 @@ void Foam::solvers::fluid::thermophysicalPredictor()
|
||||
+ pressureWork
|
||||
(
|
||||
he.name() == "e"
|
||||
? fvc::div(phi, p/rho)
|
||||
? fvc::div(phi, p/rho)()
|
||||
: -dpdt
|
||||
)
|
||||
+ thermophysicalTransport->divq(he)
|
||||
|
||||
@ -27,6 +27,10 @@ License
|
||||
#include "constrainHbyA.H"
|
||||
#include "constrainPressure.H"
|
||||
#include "adjustPhi.H"
|
||||
#include "fvcMeshPhi.H"
|
||||
#include "fvcFlux.H"
|
||||
#include "fvcDdt.H"
|
||||
#include "fvcSnGrad.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -58,6 +58,7 @@ See also
|
||||
#include "viscosityModel.H"
|
||||
#include "incompressibleMomentumTransportModels.H"
|
||||
#include "pressureReference.H"
|
||||
#include "IOMRFZoneList.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -24,6 +24,7 @@ License
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "incompressibleFluid.H"
|
||||
#include "fvmDiv.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -25,6 +25,7 @@ License
|
||||
|
||||
#include "incompressibleFluid.H"
|
||||
#include "CorrectPhi.H"
|
||||
#include "geometricZeroField.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -26,6 +26,7 @@ License
|
||||
#include "incompressibleVoF.H"
|
||||
#include "localEulerDdtScheme.H"
|
||||
#include "CorrectPhi.H"
|
||||
#include "geometricZeroField.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
@ -28,6 +28,11 @@ License
|
||||
#include "constrainPressure.H"
|
||||
#include "adjustPhi.H"
|
||||
#include "findRefCell.H"
|
||||
#include "fvcMeshPhi.H"
|
||||
#include "fvcFlux.H"
|
||||
#include "fvcDdt.H"
|
||||
#include "fvcSnGrad.H"
|
||||
#include "fvcReconstruct.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -27,6 +27,13 @@ License
|
||||
#include "constrainHbyA.H"
|
||||
#include "constrainPressure.H"
|
||||
#include "adjustPhi.H"
|
||||
#include "fvcMeshPhi.H"
|
||||
#include "fvcFlux.H"
|
||||
#include "fvcDdt.H"
|
||||
#include "fvcSnGrad.H"
|
||||
#include "fvcReconstruct.H"
|
||||
#include "fvcVolumeIntegrate.H"
|
||||
#include "fvmDiv.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -27,6 +27,13 @@ License
|
||||
#include "constrainHbyA.H"
|
||||
#include "constrainPressure.H"
|
||||
#include "adjustPhi.H"
|
||||
#include "fvcMeshPhi.H"
|
||||
#include "fvcFlux.H"
|
||||
#include "fvcDdt.H"
|
||||
#include "fvcSnGrad.H"
|
||||
#include "fvcReconstruct.H"
|
||||
#include "fvcVolumeIntegrate.H"
|
||||
#include "fvmDiv.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -26,6 +26,10 @@ License
|
||||
#include "isothermalFluid.H"
|
||||
#include "localEulerDdtScheme.H"
|
||||
#include "hydrostaticInitialisation.H"
|
||||
#include "fvcMeshPhi.H"
|
||||
#include "fvcVolumeIntegrate.H"
|
||||
#include "fvcReconstruct.H"
|
||||
#include "fvcSnGrad.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
@ -56,10 +60,10 @@ void Foam::solvers::isothermalFluid::continuityErrors()
|
||||
|
||||
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
|
||||
|
||||
Foam::tmp<Foam::volScalarField>
|
||||
Foam::tmp<Foam::volScalarField::Internal>
|
||||
Foam::solvers::isothermalFluid::pressureWork
|
||||
(
|
||||
const tmp<volScalarField>& work
|
||||
const tmp<volScalarField::Internal>& work
|
||||
) const
|
||||
{
|
||||
if (mesh.moving())
|
||||
@ -71,11 +75,11 @@ Foam::solvers::isothermalFluid::pressureWork
|
||||
fvc::interpolate(rho)*fvc::meshPhi(rho, U),
|
||||
p/rho,
|
||||
"div(phi,(p|rho))"
|
||||
);
|
||||
)();
|
||||
}
|
||||
else
|
||||
{
|
||||
return work;
|
||||
return move(work);
|
||||
}
|
||||
}
|
||||
|
||||
@ -178,7 +182,6 @@ Foam::solvers::isothermalFluid::isothermalFluid
|
||||
// Read the controls
|
||||
read();
|
||||
|
||||
thermo.validate("isothermalFluid", "h", "e");
|
||||
mesh.schemes().setFluxRequired(p.name());
|
||||
momentumTransport->validate();
|
||||
|
||||
|
||||
@ -59,6 +59,7 @@ See also
|
||||
#include "compressibleMomentumTransportModels.H"
|
||||
#include "buoyancy.H"
|
||||
#include "pressureReference.H"
|
||||
#include "IOMRFZoneList.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
@ -94,7 +95,7 @@ protected:
|
||||
|
||||
//- Rate of change of the pressure
|
||||
// Used in the enthalpy equation
|
||||
volScalarField dpdt;
|
||||
volScalarField::Internal dpdt;
|
||||
|
||||
|
||||
// Optional buoyancy
|
||||
@ -204,7 +205,10 @@ private:
|
||||
protected:
|
||||
|
||||
//- Adds the mesh-motion work to the pressure work term provided
|
||||
tmp<volScalarField> pressureWork(const tmp<volScalarField>&) const;
|
||||
tmp<volScalarField::Internal> pressureWork
|
||||
(
|
||||
const tmp<volScalarField::Internal>&
|
||||
) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
@ -24,6 +24,7 @@ License
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "isothermalFluid.H"
|
||||
#include "fvmDiv.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
@ -24,6 +24,7 @@ License
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "isothermalFluid.H"
|
||||
#include "fvcFlux.H"
|
||||
#include "fvcSmooth.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
@ -66,6 +66,8 @@ Foam::solvers::multicomponentFluid::multicomponentFluid(fvMesh& mesh)
|
||||
)
|
||||
)
|
||||
{
|
||||
thermo.validate(type(), "h", "e");
|
||||
|
||||
forAll(Y, i)
|
||||
{
|
||||
fields.add(Y[i]);
|
||||
|
||||
@ -24,6 +24,7 @@ License
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "multicomponentFluid.H"
|
||||
#include "fvcDdt.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
||||
|
||||
@ -80,7 +81,7 @@ void Foam::solvers::multicomponentFluid::thermophysicalPredictor()
|
||||
+ pressureWork
|
||||
(
|
||||
he.name() == "e"
|
||||
? mvConvection->fvcDiv(phi, p/rho)
|
||||
? mvConvection->fvcDiv(phi, p/rho)()
|
||||
: -dpdt
|
||||
)
|
||||
+ thermophysicalTransport->divq(he)
|
||||
|
||||
@ -25,8 +25,8 @@ License
|
||||
|
||||
#include "cavitation.H"
|
||||
#include "phaseSystem.H"
|
||||
#include "zeroGradientFvPatchFields.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
#include "fvCFD.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
Some files were not shown because too many files have changed in this diff Show More
Reference in New Issue
Block a user