Commit e6b67f67 authored by Andrew Heather's avatar Andrew Heather
Browse files

ENH: Clean-up after latest Foundation integrations

parent 45381b10
{
//const volScalarField& psi = thermo.psi();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
......@@ -104,9 +102,6 @@
rho = thermo.rho();
thermo.rho() = max(thermo.rho(), rhoMin);
thermo.rho() = min(thermo.rho(), rhoMax);
if (!simple.transonic())
{
rho.relax();
......
......@@ -112,8 +112,6 @@ p.correctBoundaryConditions();
// Recalculate density from the relaxed pressure
rho = thermo.rho();
thermo.rho() = max(thermo.rho(), rhoMin);
thermo.rho() = min(thermo.rho(), rhoMax);
if (!simple.transonic())
{
......
......@@ -90,9 +90,5 @@
}
rho = thermo.rho();
thermo.rho() = max(thermo.rho(), rhoMin);
thermo.rho() = min(thermo.rho(), rhoMax);
rho.relax();
}
......@@ -12,9 +12,9 @@
- fvm::laplacian(turbulence->alphaEff(), he)
==
rho*(U&g)
+ Qdot
+ parcels.Sh(he)
+ radiation->Sh(thermo)
+ combustion->Sh()
+ fvOptions(rho, he)
);
......@@ -25,6 +25,7 @@
EEqn.solve();
fvOptions.correct(he);
thermo.correct();
radiation->correct();
......
......@@ -10,10 +10,7 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
......@@ -37,11 +34,6 @@ EXE_LIBS = \
-llagrangianTurbulence \
-lspecie \
-lfluidThermophysicalModels \
-lliquidProperties \
-lliquidMixtureProperties \
-lsolidProperties \
-lsolidMixtureProperties \
-lthermophysicalFunctions \
-lreactionThermophysicalModels \
-lSLGThermo \
-lchemistryModel \
......
......@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
{
combustion->correct();
dQ = combustion->dQ();
Qdot = combustion->Qdot();
volScalarField Yt(0.0*Y[0]);
forAll(Y, i)
......
......@@ -103,18 +103,18 @@ forAll(Y, i)
}
fields.add(thermo.he());
volScalarField dQ
volScalarField Qdot
(
IOobject
(
"dQ",
"Qdot",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
dimensionedScalar("Qdot", dimEnergy/dimTime, 0.0)
);
#include "createMRF.H"
......
......@@ -38,6 +38,9 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
......@@ -45,7 +48,7 @@ Description
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "fixedFluxPressureFvPatchScalarField.H"
#include "fvcSmooth.H"
#include "basicKinematicMPPICCloud.H"
......@@ -59,16 +62,22 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
#include "createTimeControls.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
turbulence->validate();
if (!LTS)
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
......@@ -76,9 +85,17 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
}
runTime++;
......@@ -133,6 +150,8 @@ int main(int argc, char *argv[])
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
mixture.correct();
#include "UEqn.H"
// --- Pressure corrector loop
......
......@@ -4,6 +4,7 @@ EXE_INC = \
-I. \
-I./IncompressibleTwoPhaseMixtureTurbulenceModels/lnInclude \
-I$(interFoamPath) \
-I../VoF \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
......@@ -31,7 +32,7 @@ EXE_LIBS = \
-lmeshTools \
-llagrangian \
-llagrangianIntermediate \
-lthermophysicalFunctions \
-lthermophysicalProperties \
-lspecie \
-lincompressibleTransportModels \
-limmiscibleIncompressibleTwoPhaseMixture \
......
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
// Standard face-flux compression coefficient
surfaceScalarField phic
(
mixture.cAlpha()*mag(alphaPhic/mesh.magSf())
);
// Add the optional isotropic compression contribution
if (icAlpha > 0)
{
phic *= (1.0 - icAlpha);
phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
}
// Do not compress interface at non-coupled boundary faces
// (inlets, outlets etc.)
surfaceScalarField::Boundary& phicBf = phic.boundaryFieldRef();
forAll(phic.boundaryField(), patchi)
{
fvsPatchScalarField& phicp = phicBf[patchi];
if (!phicp.coupled())
{
phicp == 0;
}
}
tmp<surfaceScalarField> tphiAlpha;
if (MULESCorr)
{
fvScalarMatrix alpha1Eqn
(
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alphac, alpha1)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
alphaPhic,
upwind<scalar>(mesh, alphaPhic)
).fvmDiv(alphaPhic, alpha1)
- fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), alpha1)
);
alpha1Eqn.solve();
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
alphaPhi = tphiAlphaUD();
if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
{
Info<< "Applying the previous iteration compression flux" << endl;
MULES::correct
(
alphac,
alpha1,
alphaPhi,
tphiAlphaCorr0.ref(),
zeroField(), zeroField(),
1, 0
);
alphaPhi += tphiAlphaCorr0();
}
// Cache the upwind-flux
tphiAlphaCorr0 = tphiAlphaUD;
alpha2 = 1.0 - alpha1;
mixture.correct();
}
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
surfaceScalarField phir(phic*mixture.nHatf());
tmp<surfaceScalarField> tphiAlphaUn
(
fvc::flux
(
alphaPhic,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
);
if (MULESCorr)
{
tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - alphaPhi);
volScalarField alpha10("alpha10", alpha1);
//MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
MULES::correct
(
alphac,
alpha1,
tphiAlphaUn(),
tphiAlphaCorr.ref(),
zeroField(), zeroField(),
1, 0
);
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
alphaPhi += tphiAlphaCorr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
alphaPhi += 0.5*tphiAlphaCorr();
}
}
else
{
alphaPhi = tphiAlphaUn;
MULES::explicitSolve
(
alphac,
alpha1,
alphaPhic,
alphaPhi,
zeroField(), zeroField(),
1, 0
);
}
alpha2 = 1.0 - alpha1;
mixture.correct();
}
rhoPhi = alphaPhi*(rho1 - rho2) + alphaPhic*rho2;
if (alphaApplyPrevCorr && MULESCorr)
{
tphiAlphaCorr0 = alphaPhi - tphiAlphaCorr0;
}
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum
(
IOobject
(
"rhoPhiSum",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", rhoPhi.dimensions(), 0)
);
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
#include "alphaEqn.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
#include "alphaEqn.H"
}
rho == alpha1*rho1 + alpha2*rho2;
mu = mixture.mu();
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createRDeltaT.H"
Info<< "Reading field U\n" << endl;
volVectorField U
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading transportProperties\n" << endl;
immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
#include "createPhi.H"
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
Info<< "Reading transportProperties\n" << endl;
immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store mu as incompressibleTwoPhaseMixture does not store it
volScalarField mu
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
IOobject
(
"mu",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
mixture.mu(),
calculatedFvPatchScalarField::typeName
);
// Mass flux
surfaceScalarField rhoPhi
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
alpha1*rho1 + alpha2*rho2
);
rho.oldTime();
// Need to store mu as incompressibleTwoPhaseMixture does not store it
volScalarField mu
(
IOobject
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
fvc::interpolate(rho)*phi
);
#include "readGravitationalAcceleration.H"
"mu",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
mixture.mu(),
calculatedFvPatchScalarField::typeName
);
#include "readhRef.H"
#include "gh.H"
volScalarField p
// Mass flux
surfaceScalarField rhoPhi
(
IOobject
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
"rhoPhi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
fvc::interpolate(rho)*phi
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p
(
IOobject
(
p,
p_rgh,
mesh.solutionDict().subDict("PIMPLE"),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),