mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
ENH: centralize more libraries in src/phaseSystemModels
- prelude to code refactoring
NOTE
no source code change in this commit, only relocation,
renaming and adjustment of Make/{files,options}
This commit is contained in:
Mark Olesen
committed by
Sergio Ferraris
Sergio Ferraris
parent
42ce617b43
commit
03526e2097
@ -0,0 +1,12 @@
|
||||
#include "CourantNo.H"
|
||||
|
||||
{
|
||||
scalar UrCoNum = 0.5*gMax
|
||||
(
|
||||
fvc::surfaceSum(mag(phi1 - phi2))().primitiveField()/mesh.V().field()
|
||||
)*runTime.deltaTValue();
|
||||
|
||||
Info<< "Max Ur Courant Number = " << UrCoNum << endl;
|
||||
|
||||
CoNum = max(CoNum, UrCoNum);
|
||||
}
|
||||
@ -0,0 +1,47 @@
|
||||
for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
|
||||
{
|
||||
fluid.correctEnergyTransport();
|
||||
|
||||
autoPtr<phaseSystem::heatTransferTable>
|
||||
heatTransferPtr(fluid.heatTransfer());
|
||||
|
||||
phaseSystem::heatTransferTable&
|
||||
heatTransfer = heatTransferPtr();
|
||||
|
||||
if (!phase1.isothermal())
|
||||
{
|
||||
fvScalarMatrix E1Eqn
|
||||
(
|
||||
phase1.heEqn()
|
||||
==
|
||||
*heatTransfer[phase1.name()]
|
||||
+ alpha1*rho1*(U1&g)
|
||||
+ fvOptions(alpha1, rho1, thermo1.he())
|
||||
);
|
||||
|
||||
E1Eqn.relax();
|
||||
fvOptions.constrain(E1Eqn);
|
||||
E1Eqn.solve();
|
||||
fvOptions.correct(thermo1.he());
|
||||
}
|
||||
|
||||
if (!phase2.isothermal())
|
||||
{
|
||||
fvScalarMatrix E2Eqn
|
||||
(
|
||||
phase2.heEqn()
|
||||
==
|
||||
*heatTransfer[phase2.name()]
|
||||
+ alpha2*rho2*(U2&g)
|
||||
+ fvOptions(alpha2, rho2, phase2.thermoRef().he())
|
||||
);
|
||||
|
||||
E2Eqn.relax();
|
||||
fvOptions.constrain(E2Eqn);
|
||||
E2Eqn.solve();
|
||||
fvOptions.correct(thermo2.he());
|
||||
}
|
||||
|
||||
fluid.correctThermo();
|
||||
fluid.correct();
|
||||
}
|
||||
@ -0,0 +1,3 @@
|
||||
reactingTwoPhaseEulerFoam.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/reactingTwoPhaseEulerFoam
|
||||
@ -0,0 +1,27 @@
|
||||
phaseSystem = $(LIB_SRC)/phaseSystemModels/reactingEuler
|
||||
|
||||
EXE_INC = \
|
||||
-I${phaseSystem}/twoPhaseSystem/lnInclude \
|
||||
-I${phaseSystem}/twoPhaseCompressibleTurbulenceModels/lnInclude \
|
||||
-I${phaseSystem}/phaseSystems/lnInclude \
|
||||
-I${phaseSystem}/interfacialCompositionModels/lnInclude \
|
||||
-I${phaseSystem}/interfacialModels/lnInclude \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
-I$(LIB_SRC)/sampling/lnInclude \
|
||||
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
|
||||
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
|
||||
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
|
||||
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
|
||||
-I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lfiniteVolume \
|
||||
-lfvOptions \
|
||||
-lmeshTools \
|
||||
-lsampling \
|
||||
-lreactingPhaseSystem \
|
||||
-lreactingTwoPhaseSystem \
|
||||
-lreactingEulerianInterfacialModels \
|
||||
-lreactingEulerianInterfacialCompositionModels \
|
||||
-ltwoPhaseReactingTurbulenceModels
|
||||
@ -0,0 +1,45 @@
|
||||
{
|
||||
autoPtr<phaseSystem::massTransferTable>
|
||||
massTransferPtr(fluid.massTransfer());
|
||||
|
||||
phaseSystem::massTransferTable&
|
||||
massTransfer(massTransferPtr());
|
||||
|
||||
if (!phase1.pure())
|
||||
{
|
||||
UPtrList<volScalarField>& Y1 = phase1.YActiveRef();
|
||||
|
||||
forAll(Y1, i)
|
||||
{
|
||||
fvScalarMatrix Y1iEqn
|
||||
(
|
||||
phase1.YiEqn(Y1[i])
|
||||
==
|
||||
*massTransfer[Y1[i].name()]
|
||||
+ fvOptions(alpha1, rho1, Y1[i])
|
||||
);
|
||||
|
||||
Y1iEqn.relax();
|
||||
Y1iEqn.solve(mesh.solver("Yi"));
|
||||
}
|
||||
}
|
||||
|
||||
if (!phase2.pure())
|
||||
{
|
||||
UPtrList<volScalarField>& Y2 = phase2.YActiveRef();
|
||||
|
||||
forAll(Y2, i)
|
||||
{
|
||||
fvScalarMatrix Y2iEqn
|
||||
(
|
||||
phase2.YiEqn(Y2[i])
|
||||
==
|
||||
*massTransfer[Y2[i].name()]
|
||||
+ fvOptions(alpha2, rho2, Y2[i])
|
||||
);
|
||||
|
||||
Y2iEqn.relax();
|
||||
Y2iEqn.solve(mesh.solver("Yi"));
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,27 @@
|
||||
phaseModel& phase1 = fluid.phase1();
|
||||
phaseModel& phase2 = fluid.phase2();
|
||||
|
||||
const volScalarField& alpha1 = phase1;
|
||||
const volScalarField& alpha2 = phase2;
|
||||
|
||||
volVectorField& U1 = phase1.URef();
|
||||
surfaceScalarField& phi1 = phase1.phiRef();
|
||||
const surfaceScalarField& alphaPhi1 = phase1.alphaPhi();
|
||||
|
||||
volVectorField& U2 = phase2.URef();
|
||||
surfaceScalarField& phi2 = phase2.phiRef();
|
||||
const surfaceScalarField& alphaPhi2 = phase2.alphaPhi();
|
||||
|
||||
surfaceScalarField& phi = fluid.phi();
|
||||
|
||||
rhoThermo& thermo1 = phase1.thermoRef();
|
||||
rhoThermo& thermo2 = phase2.thermoRef();
|
||||
|
||||
volScalarField& rho1 = thermo1.rho();
|
||||
const volScalarField& psi1 = thermo1.psi();
|
||||
|
||||
volScalarField& rho2 = thermo2.rho();
|
||||
const volScalarField& psi2 = thermo2.psi();
|
||||
|
||||
const IOMRFZoneList& MRF = fluid.MRF();
|
||||
fv::options& fvOptions = fluid.fvOptions();
|
||||
@ -0,0 +1,48 @@
|
||||
#include "createRDeltaT.H"
|
||||
#include "readGravitationalAcceleration.H"
|
||||
#include "readhRef.H"
|
||||
|
||||
Info<< "Creating phaseSystem\n" << endl;
|
||||
|
||||
autoPtr<twoPhaseSystem> fluidPtr
|
||||
(
|
||||
twoPhaseSystem::New(mesh)
|
||||
);
|
||||
twoPhaseSystem& fluid = fluidPtr();
|
||||
|
||||
dimensionedScalar pMin
|
||||
(
|
||||
"pMin",
|
||||
dimPressure,
|
||||
fluid
|
||||
);
|
||||
|
||||
#include "gh.H"
|
||||
|
||||
volScalarField& p = fluid.phase1().thermoRef().p();
|
||||
|
||||
Info<< "Reading field p_rgh\n" << endl;
|
||||
volScalarField p_rgh
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"p_rgh",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
label pRefCell = 0;
|
||||
scalar pRefValue = 0.0;
|
||||
setRefCell
|
||||
(
|
||||
p,
|
||||
p_rgh,
|
||||
pimple.dict(),
|
||||
pRefCell,
|
||||
pRefValue
|
||||
);
|
||||
mesh.setFluxRequired(p_rgh.name());
|
||||
@ -0,0 +1,38 @@
|
||||
Info<< "Constructing momentum equations" << endl;
|
||||
|
||||
fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVolume/dimTime);
|
||||
fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVolume/dimTime);
|
||||
|
||||
{
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
momentumTransferPtr(fluid.momentumTransfer());
|
||||
|
||||
phaseSystem::momentumTransferTable&
|
||||
momentumTransfer(momentumTransferPtr());
|
||||
|
||||
{
|
||||
U1Eqn =
|
||||
(
|
||||
phase1.UEqn()
|
||||
==
|
||||
*momentumTransfer[phase1.name()]
|
||||
+ fvOptions(alpha1, rho1, U1)
|
||||
);
|
||||
U1Eqn.relax();
|
||||
fvOptions.constrain(U1Eqn);
|
||||
fvOptions.correct(U1);
|
||||
}
|
||||
|
||||
{
|
||||
U2Eqn =
|
||||
(
|
||||
phase2.UEqn()
|
||||
==
|
||||
*momentumTransfer[phase2.name()]
|
||||
+ fvOptions(alpha2, rho2, U2)
|
||||
);
|
||||
U2Eqn.relax();
|
||||
fvOptions.constrain(U2Eqn);
|
||||
fvOptions.correct(U2);
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,429 @@
|
||||
const surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1));
|
||||
const surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
|
||||
|
||||
PtrList<volScalarField> rAUs;
|
||||
rAUs.append
|
||||
(
|
||||
new volScalarField
|
||||
(
|
||||
IOobject::groupName("rAU", phase1.name()),
|
||||
1.0
|
||||
/(
|
||||
U1Eqn.A()
|
||||
+ byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
|
||||
)
|
||||
)
|
||||
);
|
||||
rAUs.append
|
||||
(
|
||||
new volScalarField
|
||||
(
|
||||
IOobject::groupName("rAU", phase2.name()),
|
||||
1.0
|
||||
/(
|
||||
U2Eqn.A()
|
||||
+ byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
|
||||
)
|
||||
)
|
||||
);
|
||||
const volScalarField& rAU1 = rAUs[0];
|
||||
const volScalarField& rAU2 = rAUs[1];
|
||||
|
||||
const surfaceScalarField alpharAUf1
|
||||
(
|
||||
fvc::interpolate(max(alpha1, phase1.residualAlpha())*rAU1)
|
||||
);
|
||||
const surfaceScalarField alpharAUf2
|
||||
(
|
||||
fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2)
|
||||
);
|
||||
|
||||
// Drag coefficients
|
||||
const volScalarField Kd(fluid.Kd());
|
||||
const volScalarField rAUKd1(rAU1*Kd);
|
||||
const volScalarField rAUKd2(rAU2*Kd);
|
||||
const surfaceScalarField phiKd1(fvc::interpolate(rAUKd1));
|
||||
const surfaceScalarField phiKd2(fvc::interpolate(rAUKd2));
|
||||
|
||||
// Explicit force fluxes
|
||||
PtrList<surfaceScalarField> phiFs(fluid.phiFs(rAUs));
|
||||
const surfaceScalarField& phiF1 = phiFs[0];
|
||||
const surfaceScalarField& phiF2 = phiFs[1];
|
||||
|
||||
// --- Pressure corrector loop
|
||||
while (pimple.correct())
|
||||
{
|
||||
volScalarField rho("rho", fluid.rho());
|
||||
|
||||
// Correct p_rgh for consistency with p and the updated densities
|
||||
p_rgh = p - rho*gh;
|
||||
|
||||
// Correct fixed-flux BCs to be consistent with the velocity BCs
|
||||
MRF.correctBoundaryFlux(U1, phi1);
|
||||
MRF.correctBoundaryFlux(U2, phi2);
|
||||
|
||||
// Combined buoyancy and force fluxes
|
||||
const surfaceScalarField ghSnGradRho
|
||||
(
|
||||
"ghSnGradRho",
|
||||
ghf*fvc::snGrad(rho)*mesh.magSf()
|
||||
);
|
||||
|
||||
const surfaceScalarField phigF1
|
||||
(
|
||||
alpharAUf1
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
|
||||
)
|
||||
+ phiF1
|
||||
);
|
||||
|
||||
const surfaceScalarField phigF2
|
||||
(
|
||||
alpharAUf2
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
|
||||
)
|
||||
+ phiF2
|
||||
);
|
||||
|
||||
// Predicted velocities
|
||||
volVectorField HbyA1
|
||||
(
|
||||
IOobject::groupName("HbyA", phase1.name()),
|
||||
U1
|
||||
);
|
||||
HbyA1 =
|
||||
rAU1
|
||||
*(
|
||||
U1Eqn.H()
|
||||
+ byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
|
||||
*U1.oldTime()
|
||||
);
|
||||
|
||||
volVectorField HbyA2
|
||||
(
|
||||
IOobject::groupName("HbyA", phase2.name()),
|
||||
U2
|
||||
);
|
||||
HbyA2 =
|
||||
rAU2
|
||||
*(
|
||||
U2Eqn.H()
|
||||
+ byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
|
||||
*U2.oldTime()
|
||||
);
|
||||
|
||||
// Correction force fluxes
|
||||
PtrList<surfaceScalarField> ddtCorrByAs(fluid.ddtCorrByAs(rAUs));
|
||||
|
||||
// Predicted fluxes
|
||||
const surfaceScalarField phiHbyA1
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase1.name()),
|
||||
fvc::flux(HbyA1) - phigF1 - ddtCorrByAs[0]
|
||||
);
|
||||
|
||||
const surfaceScalarField phiHbyA2
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase2.name()),
|
||||
fvc::flux(HbyA2) - phigF2 - ddtCorrByAs[1]
|
||||
);
|
||||
|
||||
ddtCorrByAs.clear();
|
||||
|
||||
// Drag fluxes
|
||||
PtrList<surfaceScalarField> phiKdPhis(fluid.phiKdPhis(rAUs));
|
||||
const surfaceScalarField& phiKdPhi1 = phiKdPhis[0];
|
||||
const surfaceScalarField& phiKdPhi2 = phiKdPhis[1];
|
||||
|
||||
// Total predicted flux
|
||||
surfaceScalarField phiHbyA
|
||||
(
|
||||
"phiHbyA",
|
||||
alphaf1*(phiHbyA1 - phiKdPhi1) + alphaf2*(phiHbyA2 - phiKdPhi2)
|
||||
);
|
||||
|
||||
MRF.makeRelative(phiHbyA);
|
||||
|
||||
phiKdPhis.clear();
|
||||
|
||||
// Construct pressure "diffusivity"
|
||||
const surfaceScalarField rAUf
|
||||
(
|
||||
"rAUf",
|
||||
mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2)
|
||||
);
|
||||
|
||||
// Update the fixedFluxPressure BCs to ensure flux consistency
|
||||
setSnGrad<fixedFluxPressureFvPatchScalarField>
|
||||
(
|
||||
p_rgh.boundaryFieldRef(),
|
||||
(
|
||||
phiHbyA.boundaryField()
|
||||
- (
|
||||
alphaf1.boundaryField()*phi1.boundaryField()
|
||||
+ alphaf2.boundaryField()*phi2.boundaryField()
|
||||
)
|
||||
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
|
||||
);
|
||||
|
||||
// Construct the compressibility parts of the pressure equation
|
||||
tmp<fvScalarMatrix> pEqnComp1, pEqnComp2;
|
||||
if (phase1.compressible())
|
||||
{
|
||||
if (pimple.transonic())
|
||||
{
|
||||
const surfaceScalarField phid1
|
||||
(
|
||||
IOobject::groupName("phid", phase1.name()),
|
||||
fvc::interpolate(psi1)*phi1
|
||||
);
|
||||
|
||||
pEqnComp1 =
|
||||
(
|
||||
fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
|
||||
)/rho1
|
||||
+ correction
|
||||
(
|
||||
(alpha1/rho1)*
|
||||
(
|
||||
psi1*fvm::ddt(p_rgh)
|
||||
+ fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
|
||||
)
|
||||
);
|
||||
|
||||
pEqnComp1.ref().relax();
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp1 =
|
||||
(
|
||||
fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
|
||||
)/rho1
|
||||
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
|
||||
}
|
||||
}
|
||||
if (phase2.compressible())
|
||||
{
|
||||
if (pimple.transonic())
|
||||
{
|
||||
const surfaceScalarField phid2
|
||||
(
|
||||
IOobject::groupName("phid", phase2.name()),
|
||||
fvc::interpolate(psi2)*phi2
|
||||
);
|
||||
|
||||
pEqnComp2 =
|
||||
(
|
||||
fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
|
||||
)/rho2
|
||||
+ correction
|
||||
(
|
||||
(alpha2/rho2)*
|
||||
(
|
||||
psi2*fvm::ddt(p_rgh)
|
||||
+ fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
|
||||
)
|
||||
);
|
||||
|
||||
pEqnComp2.ref().relax();
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp2 =
|
||||
(
|
||||
fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
|
||||
)/rho2
|
||||
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
|
||||
}
|
||||
}
|
||||
|
||||
// Add option sources
|
||||
{
|
||||
if (fvOptions.appliesToField(rho1.name()))
|
||||
{
|
||||
tmp<fvScalarMatrix> optEqn1 = fvOptions(alpha1, rho1);
|
||||
if (pEqnComp1.valid())
|
||||
{
|
||||
pEqnComp1.ref() -= (optEqn1 & rho1)/rho1;
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp1 = fvm::Su(- (optEqn1 & rho1)/rho1, p_rgh);
|
||||
}
|
||||
}
|
||||
if (fvOptions.appliesToField(rho2.name()))
|
||||
{
|
||||
tmp<fvScalarMatrix> optEqn2 = fvOptions(alpha2, rho2);
|
||||
if (pEqnComp2.valid())
|
||||
{
|
||||
pEqnComp2.ref() -= (optEqn2 & rho2)/rho2;
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp2 = fvm::Su(- (optEqn2 & rho2)/rho2, p_rgh);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Add mass transfer
|
||||
{
|
||||
PtrList<volScalarField> dmdts(fluid.dmdts());
|
||||
if (dmdts.set(0))
|
||||
{
|
||||
if (pEqnComp1.valid())
|
||||
{
|
||||
pEqnComp1.ref() -= dmdts[0]/rho1;
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp1 = fvm::Su(- dmdts[0]/rho1, p_rgh);
|
||||
}
|
||||
}
|
||||
if (dmdts.set(1))
|
||||
{
|
||||
if (pEqnComp2.valid())
|
||||
{
|
||||
pEqnComp2.ref() -= dmdts[1]/rho2;
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp2 = fvm::Su(- dmdts[1]/rho2, p_rgh);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Cache p prior to solve for density update
|
||||
const volScalarField p_rgh_0(p_rgh);
|
||||
|
||||
// Iterate over the pressure equation to correct for non-orthogonality
|
||||
while (pimple.correctNonOrthogonal())
|
||||
{
|
||||
// Construct the transport part of the pressure equation
|
||||
fvScalarMatrix pEqnIncomp
|
||||
(
|
||||
fvc::div(phiHbyA)
|
||||
- fvm::laplacian(rAUf, p_rgh)
|
||||
);
|
||||
|
||||
{
|
||||
fvScalarMatrix pEqn(pEqnIncomp);
|
||||
|
||||
if (pEqnComp1.valid())
|
||||
{
|
||||
pEqn += pEqnComp1();
|
||||
}
|
||||
|
||||
if (pEqnComp2.valid())
|
||||
{
|
||||
pEqn += pEqnComp2();
|
||||
}
|
||||
|
||||
pEqn.solve();
|
||||
}
|
||||
|
||||
// Correct fluxes and velocities on last non-orthogonal iteration
|
||||
if (pimple.finalNonOrthogonalIter())
|
||||
{
|
||||
phi = phiHbyA + pEqnIncomp.flux();
|
||||
|
||||
surfaceScalarField mSfGradp
|
||||
(
|
||||
"mSfGradp",
|
||||
pEqnIncomp.flux()/rAUf
|
||||
);
|
||||
|
||||
// Partial-elimination phase-flux corrector
|
||||
{
|
||||
const surfaceScalarField phi1s
|
||||
(
|
||||
phiHbyA1 + alpharAUf1*mSfGradp
|
||||
);
|
||||
|
||||
const surfaceScalarField phi2s
|
||||
(
|
||||
phiHbyA2 + alpharAUf2*mSfGradp
|
||||
);
|
||||
|
||||
surfaceScalarField phir
|
||||
(
|
||||
((phi1s + phiKd1*phi2s) - (phi2s + phiKd2*phi1s))
|
||||
/(1 - phiKd1*phiKd2)
|
||||
);
|
||||
|
||||
phi1 = phi + alphaf2*phir;
|
||||
phi2 = phi - alphaf1*phir;
|
||||
}
|
||||
|
||||
// Set the phase dilatation rates
|
||||
if (pEqnComp1.valid())
|
||||
{
|
||||
phase1.divU(-pEqnComp1 & p_rgh);
|
||||
}
|
||||
if (pEqnComp2.valid())
|
||||
{
|
||||
phase2.divU(-pEqnComp2 & p_rgh);
|
||||
}
|
||||
|
||||
// Optionally relax pressure for velocity correction
|
||||
p_rgh.relax();
|
||||
|
||||
mSfGradp = pEqnIncomp.flux()/rAUf;
|
||||
|
||||
// Partial-elimination phase-velocity corrector
|
||||
{
|
||||
const volVectorField Us1
|
||||
(
|
||||
HbyA1
|
||||
+ fvc::reconstruct(alpharAUf1*mSfGradp - phigF1)
|
||||
);
|
||||
|
||||
const volVectorField Us2
|
||||
(
|
||||
HbyA2
|
||||
+ fvc::reconstruct(alpharAUf2*mSfGradp - phigF2)
|
||||
);
|
||||
|
||||
const volVectorField U
|
||||
(
|
||||
alpha1*(Us1 + rAUKd1*U2) + alpha2*(Us2 + rAUKd2*U1)
|
||||
);
|
||||
|
||||
const volVectorField Ur
|
||||
(
|
||||
((1 - rAUKd2)*Us1 - (1 - rAUKd1)*Us2)/(1 - rAUKd1*rAUKd2)
|
||||
);
|
||||
|
||||
U1 = U + alpha2*Ur;
|
||||
U1.correctBoundaryConditions();
|
||||
fvOptions.correct(U1);
|
||||
|
||||
U2 = U - alpha1*Ur;
|
||||
U2.correctBoundaryConditions();
|
||||
fvOptions.correct(U2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Update and limit the static pressure
|
||||
p = max(p_rgh + rho*gh, pMin);
|
||||
|
||||
// Limit p_rgh
|
||||
p_rgh = p - rho*gh;
|
||||
|
||||
// Update densities from change in p_rgh
|
||||
rho1 += psi1*(p_rgh - p_rgh_0);
|
||||
rho2 += psi2*(p_rgh - p_rgh_0);
|
||||
|
||||
// Correct p_rgh for consistency with p and the updated densities
|
||||
rho = fluid.rho();
|
||||
p_rgh = p - rho*gh;
|
||||
p_rgh.correctBoundaryConditions();
|
||||
}
|
||||
@ -0,0 +1,43 @@
|
||||
Info<< "Constructing face momentum equations" << endl;
|
||||
|
||||
fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVolume/dimTime);
|
||||
fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVolume/dimTime);
|
||||
|
||||
{
|
||||
fluid.momentumTransfer(); // !!! Update coefficients shouldn't be necessary
|
||||
// This should be done on demand
|
||||
|
||||
autoPtr<phaseSystem::momentumTransferTable>
|
||||
momentumTransferPtr(fluid.momentumTransferf());
|
||||
|
||||
phaseSystem::momentumTransferTable&
|
||||
momentumTransfer(momentumTransferPtr());
|
||||
|
||||
{
|
||||
U1Eqn =
|
||||
(
|
||||
phase1.UfEqn()
|
||||
==
|
||||
*momentumTransfer[phase1.name()]
|
||||
+ fvOptions(alpha1, rho1, U1)
|
||||
);
|
||||
U1Eqn.relax();
|
||||
fvOptions.constrain(U1Eqn);
|
||||
U1.correctBoundaryConditions();
|
||||
fvOptions.correct(U1);
|
||||
}
|
||||
|
||||
{
|
||||
U2Eqn =
|
||||
(
|
||||
phase2.UfEqn()
|
||||
==
|
||||
*momentumTransfer[phase2.name()]
|
||||
+ fvOptions(alpha2, rho2, U2)
|
||||
);
|
||||
U2Eqn.relax();
|
||||
fvOptions.constrain(U2Eqn);
|
||||
U2.correctBoundaryConditions();
|
||||
fvOptions.correct(U2);
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,21 @@
|
||||
tmp<surfaceScalarField> trDeltaTf;
|
||||
|
||||
if (LTS && faceMomentum)
|
||||
{
|
||||
trDeltaTf = tmp<surfaceScalarField>
|
||||
(
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
fv::localEulerDdt::rDeltaTfName,
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::READ_IF_PRESENT,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh,
|
||||
dimensionedScalar("one", dimless/dimTime, 1)
|
||||
)
|
||||
);
|
||||
}
|
||||
@ -0,0 +1,410 @@
|
||||
surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1));
|
||||
surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
|
||||
|
||||
surfaceScalarField alphaRhof10
|
||||
(
|
||||
"alphaRhof10",
|
||||
fvc::interpolate
|
||||
(
|
||||
max(alpha1.oldTime(), phase1.residualAlpha())
|
||||
*rho1.oldTime()
|
||||
)
|
||||
);
|
||||
|
||||
surfaceScalarField alphaRhof20
|
||||
(
|
||||
"alphaRhof20",
|
||||
fvc::interpolate
|
||||
(
|
||||
max(alpha2.oldTime(), phase2.residualAlpha())
|
||||
*rho2.oldTime()
|
||||
)
|
||||
);
|
||||
|
||||
// Drag coefficient
|
||||
const surfaceScalarField Kdf("Kdf", fluid.Kdf());
|
||||
|
||||
// Diagonal coefficients
|
||||
PtrList<surfaceScalarField> AFfs(fluid.AFfs());
|
||||
|
||||
PtrList<surfaceScalarField> rAUfs;
|
||||
rAUfs.append
|
||||
(
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("rAUf", phase1.name()),
|
||||
1.0
|
||||
/(
|
||||
byDt(alphaRhof10)
|
||||
+ fvc::interpolate(U1Eqn.A())
|
||||
+ AFfs[0]
|
||||
)
|
||||
)
|
||||
);
|
||||
rAUfs.append
|
||||
(
|
||||
new surfaceScalarField
|
||||
(
|
||||
IOobject::groupName("rAUf", phase2.name()),
|
||||
1.0
|
||||
/(
|
||||
byDt(alphaRhof20)
|
||||
+ fvc::interpolate(U2Eqn.A())
|
||||
+ AFfs[0]
|
||||
)
|
||||
)
|
||||
);
|
||||
const surfaceScalarField& rAUf1 = rAUfs[0];
|
||||
const surfaceScalarField& rAUf2 = rAUfs[1];
|
||||
|
||||
AFfs.clear();
|
||||
|
||||
// Explicit force fluxes
|
||||
PtrList<surfaceScalarField> phiFfs(fluid.phiFfs(rAUfs));
|
||||
const surfaceScalarField& phiFf1 = phiFfs[0];
|
||||
const surfaceScalarField& phiFf2 = phiFfs[1];
|
||||
|
||||
|
||||
while (pimple.correct())
|
||||
{
|
||||
volScalarField rho("rho", fluid.rho());
|
||||
|
||||
// Correct p_rgh for consistency with p and the updated densities
|
||||
p_rgh = p - rho*gh;
|
||||
|
||||
surfaceScalarField rhof1(fvc::interpolate(rho1));
|
||||
surfaceScalarField rhof2(fvc::interpolate(rho2));
|
||||
|
||||
// Correct fixed-flux BCs to be consistent with the velocity BCs
|
||||
MRF.correctBoundaryFlux(U1, phi1);
|
||||
MRF.correctBoundaryFlux(U2, phi2);
|
||||
|
||||
const surfaceScalarField alpharAUf1
|
||||
(
|
||||
IOobject::groupName("alpharAUf", phase1.name()),
|
||||
max(alphaf1, phase1.residualAlpha())*rAUf1
|
||||
);
|
||||
|
||||
const surfaceScalarField alpharAUf2
|
||||
(
|
||||
IOobject::groupName("alpharAUf", phase2.name()),
|
||||
max(alphaf2, phase2.residualAlpha())*rAUf2
|
||||
);
|
||||
|
||||
// Combined buoyancy and force fluxes
|
||||
const surfaceScalarField ghSnGradRho
|
||||
(
|
||||
"ghSnGradRho",
|
||||
ghf*fvc::snGrad(rho)*mesh.magSf()
|
||||
);
|
||||
|
||||
const surfaceScalarField phigF1
|
||||
(
|
||||
alpharAUf1
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf2*(rhof1 - rhof2)*(g & mesh.Sf())
|
||||
)
|
||||
+ phiFf1
|
||||
);
|
||||
|
||||
const surfaceScalarField phigF2
|
||||
(
|
||||
alpharAUf2
|
||||
*(
|
||||
ghSnGradRho
|
||||
- alphaf1*(rhof2 - rhof1)*(g & mesh.Sf())
|
||||
)
|
||||
+ phiFf2
|
||||
);
|
||||
|
||||
|
||||
// Predicted fluxes
|
||||
surfaceScalarField phiHbyA1
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase1.name()),
|
||||
phi1
|
||||
);
|
||||
|
||||
phiHbyA1 =
|
||||
rAUf1
|
||||
*(
|
||||
alphaRhof10*byDt(MRF.absolute(phi1.oldTime()))
|
||||
+ fvc::flux(U1Eqn.H())
|
||||
)
|
||||
- phigF1;
|
||||
|
||||
surfaceScalarField phiHbyA2
|
||||
(
|
||||
IOobject::groupName("phiHbyA", phase2.name()),
|
||||
phi2
|
||||
);
|
||||
|
||||
phiHbyA2 =
|
||||
rAUf2
|
||||
*(
|
||||
alphaRhof20*byDt(MRF.absolute(phi2.oldTime()))
|
||||
+ fvc::flux(U2Eqn.H())
|
||||
)
|
||||
- phigF2;
|
||||
|
||||
// Drag fluxes
|
||||
PtrList<surfaceScalarField> phiKdPhifs(fluid.phiKdPhifs(rAUfs));
|
||||
const surfaceScalarField& phiKdPhif1 = phiKdPhifs[0];
|
||||
const surfaceScalarField& phiKdPhif2 = phiKdPhifs[1];
|
||||
|
||||
// Total predicted flux
|
||||
surfaceScalarField phiHbyA
|
||||
(
|
||||
"phiHbyA",
|
||||
alphaf1*(phiHbyA1 - phiKdPhif1) + alphaf2*(phiHbyA2 - phiKdPhif2)
|
||||
);
|
||||
MRF.makeRelative(phiHbyA);
|
||||
|
||||
phiKdPhifs.clear();
|
||||
|
||||
// Construct pressure "diffusivity"
|
||||
const surfaceScalarField rAUf
|
||||
(
|
||||
"rAUf",
|
||||
mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2)
|
||||
);
|
||||
|
||||
|
||||
// Update the fixedFluxPressure BCs to ensure flux consistency
|
||||
setSnGrad<fixedFluxPressureFvPatchScalarField>
|
||||
(
|
||||
p_rgh.boundaryFieldRef(),
|
||||
(
|
||||
phiHbyA.boundaryField()
|
||||
- (
|
||||
alphaf1.boundaryField()*phi1.boundaryField()
|
||||
+ alphaf2.boundaryField()*phi2.boundaryField()
|
||||
)
|
||||
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
|
||||
);
|
||||
|
||||
// Construct the compressibility parts of the pressure equation
|
||||
tmp<fvScalarMatrix> pEqnComp1, pEqnComp2;
|
||||
|
||||
if (phase1.compressible())
|
||||
{
|
||||
if (pimple.transonic())
|
||||
{
|
||||
surfaceScalarField phid1
|
||||
(
|
||||
IOobject::groupName("phid", phase1.name()),
|
||||
fvc::interpolate(psi1)*phi1
|
||||
);
|
||||
|
||||
pEqnComp1 =
|
||||
(
|
||||
fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
|
||||
)/rho1
|
||||
+ correction
|
||||
(
|
||||
(alpha1/rho1)*
|
||||
(
|
||||
psi1*fvm::ddt(p_rgh)
|
||||
+ fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
|
||||
)
|
||||
);
|
||||
|
||||
pEqnComp1.ref().relax();
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp1 =
|
||||
(
|
||||
fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
|
||||
)/rho1
|
||||
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
|
||||
}
|
||||
}
|
||||
|
||||
if (phase2.compressible())
|
||||
{
|
||||
if (pimple.transonic())
|
||||
{
|
||||
surfaceScalarField phid2
|
||||
(
|
||||
IOobject::groupName("phid", phase2.name()),
|
||||
fvc::interpolate(psi2)*phi2
|
||||
);
|
||||
|
||||
pEqnComp2 =
|
||||
(
|
||||
fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
|
||||
)/rho2
|
||||
+ correction
|
||||
(
|
||||
(alpha2/rho2)*
|
||||
(
|
||||
psi2*fvm::ddt(p_rgh)
|
||||
+ fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
|
||||
)
|
||||
);
|
||||
|
||||
pEqnComp2.ref().relax();
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp2 =
|
||||
(
|
||||
fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
|
||||
)/rho2
|
||||
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
|
||||
}
|
||||
}
|
||||
|
||||
// Add option sources
|
||||
{
|
||||
if (fvOptions.appliesToField(rho1.name()))
|
||||
{
|
||||
tmp<fvScalarMatrix> optEqn1 = fvOptions(alpha1, rho1);
|
||||
if (pEqnComp1.valid())
|
||||
{
|
||||
pEqnComp1.ref() -= (optEqn1 & rho1)/rho1;
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp1 = fvm::Su(- (optEqn1 & rho1)/rho1, p_rgh);
|
||||
}
|
||||
}
|
||||
if (fvOptions.appliesToField(rho2.name()))
|
||||
{
|
||||
tmp<fvScalarMatrix> optEqn2 = fvOptions(alpha2, rho2);
|
||||
if (pEqnComp2.valid())
|
||||
{
|
||||
pEqnComp2.ref() -= (optEqn2 & rho2)/rho2;
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp2 = fvm::Su(- (optEqn2 & rho2)/rho2, p_rgh);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Add mass transfer
|
||||
{
|
||||
PtrList<volScalarField> dmdts(fluid.dmdts());
|
||||
if (dmdts.set(0))
|
||||
{
|
||||
if (pEqnComp1.valid())
|
||||
{
|
||||
pEqnComp1.ref() -= dmdts[0]/rho1;
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp1 = fvm::Su(- dmdts[0]/rho1, p_rgh);
|
||||
}
|
||||
}
|
||||
if (dmdts.set(1))
|
||||
{
|
||||
if (pEqnComp2.valid())
|
||||
{
|
||||
pEqnComp2.ref() -= dmdts[1]/rho2;
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp2 = fvm::Su(- dmdts[1]/rho2, p_rgh);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Cache p prior to solve for density update
|
||||
volScalarField p_rgh_0("p_rgh_0", p_rgh);
|
||||
|
||||
while (pimple.correctNonOrthogonal())
|
||||
{
|
||||
fvScalarMatrix pEqnIncomp
|
||||
(
|
||||
fvc::div(phiHbyA)
|
||||
- fvm::laplacian(rAUf, p_rgh)
|
||||
);
|
||||
|
||||
{
|
||||
fvScalarMatrix pEqn(pEqnIncomp);
|
||||
|
||||
if (pEqnComp1.valid())
|
||||
{
|
||||
pEqn += pEqnComp1();
|
||||
}
|
||||
|
||||
if (pEqnComp2.valid())
|
||||
{
|
||||
pEqn += pEqnComp2();
|
||||
}
|
||||
|
||||
solve
|
||||
(
|
||||
pEqn,
|
||||
mesh.solver(p_rgh.select(pimple.finalInnerIter()))
|
||||
);
|
||||
}
|
||||
|
||||
if (pimple.finalNonOrthogonalIter())
|
||||
{
|
||||
surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
|
||||
|
||||
phi = phiHbyA + pEqnIncomp.flux();
|
||||
|
||||
const surfaceScalarField phi1s
|
||||
(
|
||||
phiHbyA1 + alpharAUf1*mSfGradp
|
||||
);
|
||||
|
||||
const surfaceScalarField phi2s
|
||||
(
|
||||
phiHbyA2 + alpharAUf2*mSfGradp
|
||||
);
|
||||
|
||||
surfaceScalarField phir
|
||||
(
|
||||
((phi2s + rAUf2*Kdf*phi1s) - (phi1s + rAUf1*Kdf*phi2s))
|
||||
/(1.0 - rAUf1*rAUf2*sqr(Kdf))
|
||||
);
|
||||
|
||||
phi1 = phi - alphaf2*phir;
|
||||
phi2 = phi + alphaf1*phir;
|
||||
|
||||
U1 = fvc::reconstruct(MRF.absolute(phi1));
|
||||
U1.correctBoundaryConditions();
|
||||
fvOptions.correct(U1);
|
||||
|
||||
U2 = fvc::reconstruct(MRF.absolute(phi2));
|
||||
U2.correctBoundaryConditions();
|
||||
fvOptions.correct(U2);
|
||||
|
||||
// Set the phase dilatation rates
|
||||
if (pEqnComp1.valid())
|
||||
{
|
||||
phase1.divU(-pEqnComp1 & p_rgh);
|
||||
}
|
||||
if (pEqnComp2.valid())
|
||||
{
|
||||
phase2.divU(-pEqnComp2 & p_rgh);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Update and limit the static pressure
|
||||
p = max(p_rgh + rho*gh, pMin);
|
||||
|
||||
// Limit p_rgh
|
||||
p_rgh = p - rho*gh;
|
||||
|
||||
// Update densities from change in p_rgh
|
||||
rho1 += psi1*(p_rgh - p_rgh_0);
|
||||
rho2 += psi2*(p_rgh - p_rgh_0);
|
||||
|
||||
// Correct p_rgh for consistency with p and the updated densities
|
||||
rho = fluid.rho();
|
||||
p_rgh = p - rho*gh;
|
||||
p_rgh.correctBoundaryConditions();
|
||||
}
|
||||
@ -0,0 +1,154 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
Copyright (C) 2020 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
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
|
||||
reactingTwoPhaseEulerFoam
|
||||
|
||||
Group
|
||||
grpMultiphaseSolvers
|
||||
|
||||
Description
|
||||
Solver for a system of two compressible fluid phases with a common pressure,
|
||||
but otherwise separate properties. The type of phase model is run time
|
||||
selectable and can optionally represent multiple species and in-phase
|
||||
reactions. The phase system is also run time selectable and can optionally
|
||||
represent different types of momentum, heat and mass transfer.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
#include "twoPhaseSystem.H"
|
||||
#include "phaseCompressibleTurbulenceModel.H"
|
||||
#include "pimpleControl.H"
|
||||
#include "localEulerDdtScheme.H"
|
||||
#include "fvcSmooth.H"
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
argList::addNote
|
||||
(
|
||||
"Solver for a system of two compressible fluid phases with a"
|
||||
" common pressure, but otherwise separate properties."
|
||||
);
|
||||
|
||||
#include "postProcess.H"
|
||||
|
||||
#include "addCheckCaseOptions.H"
|
||||
#include "setRootCaseLists.H"
|
||||
#include "createTime.H"
|
||||
#include "createMesh.H"
|
||||
#include "createControl.H"
|
||||
#include "createTimeControls.H"
|
||||
#include "createFields.H"
|
||||
#include "createFieldRefs.H"
|
||||
|
||||
if (!LTS)
|
||||
{
|
||||
#include "CourantNo.H"
|
||||
#include "setInitialDeltaT.H"
|
||||
}
|
||||
|
||||
bool faceMomentum
|
||||
(
|
||||
pimple.dict().getOrDefault("faceMomentum", false)
|
||||
);
|
||||
|
||||
#include "pUf/createRDeltaTf.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
Info<< "\nStarting time loop\n" << endl;
|
||||
|
||||
while (runTime.run())
|
||||
{
|
||||
#include "readTimeControls.H"
|
||||
|
||||
int nEnergyCorrectors
|
||||
(
|
||||
pimple.dict().getOrDefault<int>("nEnergyCorrectors", 1)
|
||||
);
|
||||
|
||||
if (LTS)
|
||||
{
|
||||
#include "setRDeltaT.H"
|
||||
if (faceMomentum)
|
||||
{
|
||||
#include "setRDeltaTf.H"
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
#include "CourantNos.H"
|
||||
#include "setDeltaT.H"
|
||||
}
|
||||
|
||||
++runTime;
|
||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
||||
|
||||
// --- Pressure-velocity PIMPLE corrector loop
|
||||
while (pimple.loop())
|
||||
{
|
||||
fluid.solve();
|
||||
fluid.correct();
|
||||
|
||||
#include "YEqns.H"
|
||||
|
||||
if (faceMomentum)
|
||||
{
|
||||
#include "pUf/UEqns.H"
|
||||
#include "EEqns.H"
|
||||
#include "pUf/pEqn.H"
|
||||
}
|
||||
else
|
||||
{
|
||||
#include "pU/UEqns.H"
|
||||
#include "EEqns.H"
|
||||
#include "pU/pEqn.H"
|
||||
}
|
||||
|
||||
fluid.correctKinematics();
|
||||
|
||||
if (pimple.turbCorr())
|
||||
{
|
||||
fluid.correctTurbulence();
|
||||
}
|
||||
}
|
||||
|
||||
runTime.write();
|
||||
|
||||
runTime.printExecutionTime(Info);
|
||||
}
|
||||
|
||||
Info<< "End\n" << endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,37 @@
|
||||
{
|
||||
volScalarField& rDeltaT = trDeltaT.ref();
|
||||
|
||||
const dictionary& pimpleDict = pimple.dict();
|
||||
|
||||
scalar maxCo
|
||||
(
|
||||
pimpleDict.getOrDefault<scalar>("maxCo", 0.2)
|
||||
);
|
||||
|
||||
scalar maxDeltaT
|
||||
(
|
||||
pimpleDict.getOrDefault<scalar>("maxDeltaT", GREAT)
|
||||
);
|
||||
|
||||
scalar rDeltaTSmoothingCoeff
|
||||
(
|
||||
pimpleDict.getOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.02)
|
||||
);
|
||||
|
||||
// Set the reciprocal time-step from the local Courant number
|
||||
rDeltaT.ref() = max
|
||||
(
|
||||
1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
|
||||
fvc::surfaceSum(max(mag(phi1), mag(phi2)))()()
|
||||
/((2*maxCo)*mesh.V())
|
||||
);
|
||||
|
||||
// Update tho boundary values of the reciprocal time-step
|
||||
rDeltaT.correctBoundaryConditions();
|
||||
|
||||
fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
|
||||
|
||||
Info<< "Flow time scale min/max = "
|
||||
<< gMin(1/rDeltaT.primitiveField())
|
||||
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;
|
||||
}
|
||||
@ -0,0 +1 @@
|
||||
trDeltaTf.ref() = fvc::interpolate(fv::localEulerDdt::localRDeltaT(mesh));
|
||||
Reference in New Issue
Block a user