research/OpenFOAM-5.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

compressibleInterFoam: More consistent with interFoam and added partial support for LTS

Browse source
This commit is contained in:
Henry Weller 2017-02-08 20:50:07 +00:00
parent 11e298382e
commit 1aec2695f1
17 changed files with 165 additions and 127 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
applications/solvers/multiphase/compressibleInterFoam/Make/options
View file

@ -1,4 +1,6 @@
EXE_INC = \
-I. \
-I../interFoam \
-ItwoPhaseMixtureThermo \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \

1
applications/solvers/multiphase/compressibleInterFoam/TEqn.H
View file

@ -24,5 +24,6 @@
fvOptions.correct(T);
mixture.correctThermo();
mixture.correct();
}

2
applications/solvers/multiphase/compressibleInterFoam/UEqn.H
View file

@ -20,7 +20,7 @@
fvc::reconstruct
(
(
interface.surfaceTensionForce()
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()

41
applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H → applications/solvers/multiphase/compressibleInterFoam/alphaEqn.H
View file

@ -2,48 +2,11 @@
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir(phic*interface.nHatf());
surfaceScalarField phir(phic*mixture.nHatf());
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
volScalarField::Internal Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
);
volScalarField::Internal Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
divU*min(alpha1, scalar(1))
);
forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]*alpha1[celli];
Su[celli] += dgdt[celli]*alpha1[celli];
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
}
}
#include "alphaSuSp.H"
surfaceScalarField alphaPhi1
(

8
applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H → applications/solvers/multiphase/compressibleInterFoam/alphaEqnSubCycle.H
View file

@ -1,8 +1,6 @@
{
#include "alphaControls.H"
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
phic = min(mixture.cAlpha()*phic, max(phic));
volScalarField divU(fvc::div(fvc::absolute(phi, U)));
@ -27,7 +25,7 @@
!(++alphaSubCycle).end();
)
{
#include "alphaEqns.H"
#include "alphaEqn.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
@ -35,6 +33,6 @@
}
else
{
#include "alphaEqns.H"
#include "alphaEqn.H"
}
}

37
applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H Normal file
View file

@ -0,0 +1,37 @@
volScalarField::Internal Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
);
volScalarField::Internal Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
divU*min(alpha1, scalar(1))
);
forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]*alpha1[celli];
Su[celli] += dgdt[celli]*alpha1[celli];
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
}
}

1
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options
View file

@ -1,5 +1,6 @@
EXE_INC = \
-I.. \
-I../../interFoam \
-I../twoPhaseMixtureThermo \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \

107
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
View file

@ -41,13 +41,14 @@ Description
#include "dynamicFvMesh.H"
#include "MULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "twoPhaseMixtureThermo.H"
#include "turbulentFluidThermoModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -81,64 +82,67 @@ int main(int argc, char *argv[])
{
#include "readControls.H"
{
// Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
// Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
// Do any mesh changes
mesh.update();
if (mesh.changing())
{
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
}
if (mesh.changing() && correctPhi)
{
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & Uf;
#include "correctPhi.H"
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
}
}
if (mesh.changing() && checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
runTime++;
turbulence->correct();
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "alphaEqnsSubCycle.H"
// correct interface on first PIMPLE corrector
if (pimple.corr() == 1)
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
interface.correct();
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
mesh.update();
if (mesh.changing())
{
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
}
if (mesh.changing() && correctPhi)
{
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & Uf;
#include "correctPhi.H"
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
mixture.correct();
}
if (mesh.changing() && checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
solve(fvm::ddt(rho) + fvc::div(rhoPhi));
#include "UEqn.H"
@ -149,14 +153,13 @@ int main(int argc, char *argv[])
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
rho = alpha1*rho1 + alpha2*rho2;
// Correct p_rgh for consistency with p and the updated densities
p_rgh = p - rho*gh;
p_rgh.correctBoundaryConditions();
runTime.write();
Info<< "ExecutionTime = "

5
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/createControls.H
View file

@ -9,3 +9,8 @@ bool checkMeshCourantNo
(
pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false)
);
bool moveMeshOuterCorrectors
(
pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false)
);

2
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
View file

@ -12,7 +12,7 @@
surfaceScalarField phig
(
(
interface.surfaceTensionForce()
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);

3
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
View file

@ -4,3 +4,6 @@ correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
checkMeshCourantNo =
pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
moveMeshOuterCorrectors =
pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false);

34
applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
View file

@ -39,12 +39,13 @@ Description
#include "MULES.H"
#include "subCycle.H"
#include "rhoThermo.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "twoPhaseMixtureThermo.H"
#include "turbulentFluidThermoModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -59,8 +60,6 @@ int main(int argc, char *argv[])
#include "createTimeControls.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
volScalarField& p = mixture.p();
volScalarField& T = mixture.T();
@ -69,6 +68,13 @@ int main(int argc, char *argv[])
turbulence->validate();
if (!LTS)
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
@ -76,8 +82,17 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
}
runTime++;
@ -86,13 +101,8 @@ int main(int argc, char *argv[])
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "alphaEqnsSubCycle.H"
// correct interface on first PIMPLE corrector
if (pimple.corr() == 1)
{
interface.correct();
}
#include "alphaControls.H"
#include "alphaEqnSubCycle.H"
solve(fvm::ddt(rho) + fvc::div(rhoPhi));

7
applications/solvers/multiphase/compressibleInterFoam/createFields.H
View file

@ -1,3 +1,5 @@
#include "createRDeltaT.H"
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
@ -29,7 +31,7 @@ volVectorField U
#include "createPhi.H"
Info<< "Constructing twoPhaseMixtureThermo\n" << endl;
twoPhaseMixtureThermo mixture(mesh);
twoPhaseMixtureThermo mixture(U, phi);
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
@ -89,9 +91,6 @@ volScalarField dgdt
pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
);
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, mixture);
// Construct compressible turbulence model
autoPtr<compressible::turbulenceModel> turbulence
(

2
applications/solvers/multiphase/compressibleInterFoam/pEqn.H
View file

@ -12,7 +12,7 @@
surfaceScalarField phig
(
(
interface.surfaceTensionForce()
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);

2
applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/Make/options
View file

@ -2,6 +2,7 @@ EXE_INC = \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
LIB_LIBS = \
@ -9,4 +10,5 @@ LIB_LIBS = \
-lfluidThermophysicalModels \
-lspecie \
-ltwoPhaseMixture \
-linterfaceProperties \
-lfiniteVolume

22
applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -40,11 +40,13 @@ namespace Foam
Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo
(
const fvMesh& mesh
const volVectorField& U,
const surfaceScalarField& phi
)
:
psiThermo(mesh, word::null),
twoPhaseMixture(mesh, *this),
psiThermo(U.mesh(), word::null),
twoPhaseMixture(U.mesh(), *this),
interfaceProperties(alpha1(), U, *this),
thermo1_(nullptr),
thermo2_(nullptr)
{
@ -58,8 +60,8 @@ Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo
T2.write();
}
thermo1_ = rhoThermo::New(mesh, phase1Name());
thermo2_ = rhoThermo::New(mesh, phase2Name());
thermo1_ = rhoThermo::New(U.mesh(), phase1Name());
thermo2_ = rhoThermo::New(U.mesh(), phase2Name());
thermo1_->validate(phase1Name(), "e");
thermo2_->validate(phase2Name(), "e");
@ -76,17 +78,23 @@ Foam::twoPhaseMixtureThermo::~twoPhaseMixtureThermo()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::twoPhaseMixtureThermo::correct()
void Foam::twoPhaseMixtureThermo::correctThermo()
{
thermo1_->he() = thermo1_->he(p_, T_);
thermo1_->correct();
thermo2_->he() = thermo2_->he(p_, T_);
thermo2_->correct();
}
void Foam::twoPhaseMixtureThermo::correct()
{
psi_ = alpha1()*thermo1_->psi() + alpha2()*thermo2_->psi();
mu_ = alpha1()*thermo1_->mu() + alpha2()*thermo2_->mu();
alpha_ = alpha1()*thermo1_->alpha() + alpha2()*thermo2_->alpha();
interfaceProperties::correct();
}

16
applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -39,6 +39,7 @@ SourceFiles
#include "rhoThermo.H"
#include "psiThermo.H"
#include "twoPhaseMixture.H"
#include "interfaceProperties.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -52,7 +53,8 @@ namespace Foam
class twoPhaseMixtureThermo
:
public psiThermo,
public twoPhaseMixture
public twoPhaseMixture,
public interfaceProperties
{
// Private data
@ -71,10 +73,11 @@ public:
// Constructors
//- Construct from mesh
//- Construct from components
twoPhaseMixtureThermo
(
const fvMesh& mesh
const volVectorField& U,
const surfaceScalarField& phi
);
@ -104,7 +107,10 @@ public:
return thermo2_();
}
//- Update properties
//- Correct the thermodynamics of each phase
virtual void correctThermo();
//- Update mixture properties
virtual void correct();
//- Return true if the equation of state is incompressible