Compare commits
1 Commits
develop.me
...
feature-co
| Author | SHA1 | Date | |
|---|---|---|---|
| 54dcde5f9a |
@ -7,7 +7,6 @@ wclean libso surfaceTensionModels
|
|||||||
wclean libso VoFphaseCompressibleTurbulenceModels
|
wclean libso VoFphaseCompressibleTurbulenceModels
|
||||||
|
|
||||||
wclean
|
wclean
|
||||||
wclean compressibleInterDyMFoam
|
|
||||||
wclean compressibleInterFilmFoam
|
wclean compressibleInterFilmFoam
|
||||||
wclean compressibleInterIsoFoam
|
wclean compressibleInterIsoFoam
|
||||||
|
|
||||||
|
|||||||
@ -8,7 +8,6 @@ wmake $targetType surfaceTensionModels
|
|||||||
wmake $targetType VoFphaseCompressibleTurbulenceModels
|
wmake $targetType VoFphaseCompressibleTurbulenceModels
|
||||||
|
|
||||||
wmake $targetType
|
wmake $targetType
|
||||||
wmake $targetType compressibleInterDyMFoam
|
|
||||||
wmake $targetType compressibleInterFilmFoam
|
wmake $targetType compressibleInterFilmFoam
|
||||||
wmake $targetType compressibleInterIsoFoam
|
wmake $targetType compressibleInterIsoFoam
|
||||||
wmake $targetType overCompressibleInterDyMFoam
|
wmake $targetType overCompressibleInterDyMFoam
|
||||||
|
|||||||
@ -1,21 +1,25 @@
|
|||||||
EXE_INC = \
|
EXE_INC = \
|
||||||
-I../VoF \
|
-I../VoF \
|
||||||
-ItwoPhaseMixtureThermo \
|
-ItwoPhaseMixtureThermo \
|
||||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
-IVoFphaseCompressibleTurbulenceModels/lnInclude \
|
||||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||||
|
-I$(LIB_SRC)/dynamicMesh/lnInclude \
|
||||||
|
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
|
||||||
|
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||||
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
|
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
|
||||||
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
|
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
|
||||||
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
|
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
|
||||||
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
|
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
|
||||||
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
|
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
|
||||||
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
|
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
|
||||||
-I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
|
-I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude
|
||||||
-IVoFphaseCompressibleTurbulenceModels/lnInclude
|
|
||||||
|
|
||||||
EXE_LIBS = \
|
EXE_LIBS = \
|
||||||
-lfiniteVolume \
|
-lfiniteVolume \
|
||||||
-lfvOptions \
|
-lfvOptions \
|
||||||
-lmeshTools \
|
-lmeshTools \
|
||||||
|
-ldynamicMesh \
|
||||||
|
-ldynamicFvMesh \
|
||||||
-ltwoPhaseMixtureThermo \
|
-ltwoPhaseMixtureThermo \
|
||||||
-ltwoPhaseSurfaceTension \
|
-ltwoPhaseSurfaceTension \
|
||||||
-lcompressibleTransportModels \
|
-lcompressibleTransportModels \
|
||||||
|
|||||||
@ -1,3 +0,0 @@
|
|||||||
compressibleInterDyMFoam.C
|
|
||||||
|
|
||||||
EXE = $(FOAM_APPBIN)/compressibleInterDyMFoam
|
|
||||||
@ -1,35 +0,0 @@
|
|||||||
EXE_INC = \
|
|
||||||
-I.. \
|
|
||||||
-I../../VoF \
|
|
||||||
-I../twoPhaseMixtureThermo \
|
|
||||||
-I../VoFphaseCompressibleTurbulenceModels/lnInclude \
|
|
||||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
|
||||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
|
||||||
-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)/TurbulenceModels/turbulenceModels/lnInclude \
|
|
||||||
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
|
|
||||||
-I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
|
|
||||||
-I$(LIB_SRC)/dynamicMesh/lnInclude \
|
|
||||||
-I$(LIB_SRC)/dynamicFvMesh/lnInclude
|
|
||||||
|
|
||||||
EXE_LIBS = \
|
|
||||||
-lfiniteVolume \
|
|
||||||
-lfvOptions \
|
|
||||||
-lmeshTools \
|
|
||||||
-ltwoPhaseMixtureThermo \
|
|
||||||
-ltwoPhaseSurfaceTension \
|
|
||||||
-lcompressibleTransportModels \
|
|
||||||
-lfluidThermophysicalModels \
|
|
||||||
-lspecie \
|
|
||||||
-ltwoPhaseMixture \
|
|
||||||
-ltwoPhaseProperties \
|
|
||||||
-linterfaceProperties \
|
|
||||||
-lturbulenceModels \
|
|
||||||
-lcompressibleTurbulenceModels \
|
|
||||||
-lthermoTools \
|
|
||||||
-lVoFphaseCompressibleTurbulenceModels \
|
|
||||||
-ldynamicMesh \
|
|
||||||
-ldynamicFvMesh
|
|
||||||
@ -1,43 +0,0 @@
|
|||||||
volScalarField::Internal Sp
|
|
||||||
(
|
|
||||||
IOobject
|
|
||||||
(
|
|
||||||
"Sp",
|
|
||||||
runTime.timeName(),
|
|
||||||
mesh
|
|
||||||
),
|
|
||||||
mesh,
|
|
||||||
dimensionedScalar(dgdt.dimensions(), Zero)
|
|
||||||
);
|
|
||||||
|
|
||||||
volScalarField::Internal Su
|
|
||||||
(
|
|
||||||
IOobject
|
|
||||||
(
|
|
||||||
"Su",
|
|
||||||
runTime.timeName(),
|
|
||||||
mesh
|
|
||||||
),
|
|
||||||
mesh,
|
|
||||||
dimensionedScalar(dgdt.dimensions(), Zero)
|
|
||||||
);
|
|
||||||
|
|
||||||
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]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
volScalarField::Internal divU
|
|
||||||
(
|
|
||||||
mesh.moving()
|
|
||||||
? fvc::div(phiCN() + mesh.phi())
|
|
||||||
: fvc::div(phiCN())
|
|
||||||
);
|
|
||||||
@ -1,190 +0,0 @@
|
|||||||
/*---------------------------------------------------------------------------*\
|
|
||||||
========= |
|
|
||||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
||||||
\\ / O peration |
|
|
||||||
\\ / A nd | www.openfoam.com
|
|
||||||
\\/ M anipulation |
|
|
||||||
-------------------------------------------------------------------------------
|
|
||||||
Copyright (C) 2011-2017 OpenFOAM Foundation
|
|
||||||
Copyright (C) 2019 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
|
|
||||||
compressibleInterDyMFoam
|
|
||||||
|
|
||||||
Description
|
|
||||||
Solver for two compressible, non-isothermal immiscible fluids using a VOF
|
|
||||||
(volume of fluid) phase-fraction based interface capturing approach,
|
|
||||||
with optional mesh motion and mesh topology changes including adaptive
|
|
||||||
re-meshing.
|
|
||||||
|
|
||||||
The momentum and other fluid properties are of the "mixture" and a single
|
|
||||||
momentum equation is solved.
|
|
||||||
|
|
||||||
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
|
|
||||||
|
|
||||||
\*---------------------------------------------------------------------------*/
|
|
||||||
|
|
||||||
#include "fvCFD.H"
|
|
||||||
#include "dynamicFvMesh.H"
|
|
||||||
#include "CMULES.H"
|
|
||||||
#include "EulerDdtScheme.H"
|
|
||||||
#include "localEulerDdtScheme.H"
|
|
||||||
#include "CrankNicolsonDdtScheme.H"
|
|
||||||
#include "subCycle.H"
|
|
||||||
#include "compressibleInterPhaseTransportModel.H"
|
|
||||||
#include "pimpleControl.H"
|
|
||||||
#include "fvOptions.H"
|
|
||||||
#include "CorrectPhi.H"
|
|
||||||
#include "fvcSmooth.H"
|
|
||||||
|
|
||||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
||||||
|
|
||||||
int main(int argc, char *argv[])
|
|
||||||
{
|
|
||||||
argList::addNote
|
|
||||||
(
|
|
||||||
"Solver for two compressible, non-isothermal immiscible fluids"
|
|
||||||
" using VOF phase-fraction based interface capturing.\n"
|
|
||||||
"With optional mesh motion and mesh topology changes including"
|
|
||||||
" adaptive re-meshing."
|
|
||||||
);
|
|
||||||
|
|
||||||
#include "postProcess.H"
|
|
||||||
|
|
||||||
#include "setRootCaseLists.H"
|
|
||||||
#include "createTime.H"
|
|
||||||
#include "createDynamicFvMesh.H"
|
|
||||||
#include "initContinuityErrs.H"
|
|
||||||
#include "createDyMControls.H"
|
|
||||||
#include "createFields.H"
|
|
||||||
#include "createUf.H"
|
|
||||||
#include "CourantNo.H"
|
|
||||||
#include "setInitialDeltaT.H"
|
|
||||||
|
|
||||||
volScalarField& p = mixture.p();
|
|
||||||
volScalarField& T = mixture.T();
|
|
||||||
const volScalarField& psi1 = mixture.thermo1().psi();
|
|
||||||
const volScalarField& psi2 = mixture.thermo2().psi();
|
|
||||||
|
|
||||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
||||||
Info<< "\nStarting time loop\n" << endl;
|
|
||||||
|
|
||||||
while (runTime.run())
|
|
||||||
{
|
|
||||||
#include "readDyMControls.H"
|
|
||||||
|
|
||||||
// Store divU and divUp 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)));
|
|
||||||
volScalarField divUp("divUp", fvc::div(fvc::absolute(phi, U), p));
|
|
||||||
|
|
||||||
if (LTS)
|
|
||||||
{
|
|
||||||
#include "setRDeltaT.H"
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
#include "CourantNo.H"
|
|
||||||
#include "alphaCourantNo.H"
|
|
||||||
#include "setDeltaT.H"
|
|
||||||
}
|
|
||||||
|
|
||||||
++runTime;
|
|
||||||
|
|
||||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
|
||||||
|
|
||||||
// --- Pressure-velocity PIMPLE corrector loop
|
|
||||||
while (pimple.loop())
|
|
||||||
{
|
|
||||||
if (pimple.firstIter() || moveMeshOuterCorrectors)
|
|
||||||
{
|
|
||||||
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
|
|
||||||
|
|
||||||
mesh.update();
|
|
||||||
|
|
||||||
if (mesh.changing())
|
|
||||||
{
|
|
||||||
MRF.update();
|
|
||||||
|
|
||||||
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 "compressibleAlphaEqnSubCycle.H"
|
|
||||||
|
|
||||||
turbulence.correctPhasePhi();
|
|
||||||
|
|
||||||
#include "UEqn.H"
|
|
||||||
#include "TEqn.H"
|
|
||||||
|
|
||||||
// --- Pressure corrector loop
|
|
||||||
while (pimple.correct())
|
|
||||||
{
|
|
||||||
#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();
|
|
||||||
|
|
||||||
runTime.printExecutionTime(Info);
|
|
||||||
}
|
|
||||||
|
|
||||||
Info<< "End\n" << endl;
|
|
||||||
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
// ************************************************************************* //
|
|
||||||
@ -1,145 +0,0 @@
|
|||||||
{
|
|
||||||
volScalarField rAU("rAU", 1.0/UEqn.A());
|
|
||||||
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
|
|
||||||
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
|
|
||||||
surfaceScalarField phiHbyA
|
|
||||||
(
|
|
||||||
"phiHbyA",
|
|
||||||
fvc::flux(HbyA)
|
|
||||||
+ MRF.zeroFilter(fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf))
|
|
||||||
);
|
|
||||||
MRF.makeRelative(phiHbyA);
|
|
||||||
|
|
||||||
surfaceScalarField phig
|
|
||||||
(
|
|
||||||
(
|
|
||||||
mixture.surfaceTensionForce()
|
|
||||||
- ghf*fvc::snGrad(rho)
|
|
||||||
)*rAUf*mesh.magSf()
|
|
||||||
);
|
|
||||||
|
|
||||||
phiHbyA += phig;
|
|
||||||
|
|
||||||
// Update the pressure BCs to ensure flux consistency
|
|
||||||
constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
|
|
||||||
|
|
||||||
// Make the fluxes relative to the mesh motion
|
|
||||||
fvc::makeRelative(phiHbyA, U);
|
|
||||||
|
|
||||||
tmp<fvScalarMatrix> p_rghEqnComp1;
|
|
||||||
tmp<fvScalarMatrix> p_rghEqnComp2;
|
|
||||||
|
|
||||||
if (pimple.transonic())
|
|
||||||
{
|
|
||||||
#include "rhofs.H"
|
|
||||||
|
|
||||||
surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
|
|
||||||
surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
|
|
||||||
|
|
||||||
p_rghEqnComp1 =
|
|
||||||
pos(alpha1)
|
|
||||||
*(
|
|
||||||
(
|
|
||||||
fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
|
|
||||||
- (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
|
|
||||||
)/rho1
|
|
||||||
- fvc::ddt(alpha1) - fvc::div(alphaPhi1)
|
|
||||||
+ (alpha1/rho1)
|
|
||||||
*correction
|
|
||||||
(
|
|
||||||
psi1*fvm::ddt(p_rgh)
|
|
||||||
+ fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
|
|
||||||
)
|
|
||||||
);
|
|
||||||
p_rghEqnComp1.ref().relax();
|
|
||||||
|
|
||||||
p_rghEqnComp2 =
|
|
||||||
pos(alpha2)
|
|
||||||
*(
|
|
||||||
(
|
|
||||||
fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
|
|
||||||
- (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
|
|
||||||
)/rho2
|
|
||||||
- fvc::ddt(alpha2) - fvc::div(alphaPhi2)
|
|
||||||
+ (alpha2/rho2)
|
|
||||||
*correction
|
|
||||||
(
|
|
||||||
psi2*fvm::ddt(p_rgh)
|
|
||||||
+ fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
|
|
||||||
)
|
|
||||||
);
|
|
||||||
p_rghEqnComp2.ref().relax();
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
p_rghEqnComp1 =
|
|
||||||
fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
|
|
||||||
+ fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
|
|
||||||
|
|
||||||
p_rghEqnComp2 =
|
|
||||||
fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
|
|
||||||
+ fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Cache p_rgh prior to solve for density update
|
|
||||||
volScalarField p_rgh_0(p_rgh);
|
|
||||||
|
|
||||||
while (pimple.correctNonOrthogonal())
|
|
||||||
{
|
|
||||||
fvScalarMatrix p_rghEqnIncomp
|
|
||||||
(
|
|
||||||
fvc::div(phiHbyA)
|
|
||||||
- fvm::laplacian(rAUf, p_rgh)
|
|
||||||
);
|
|
||||||
|
|
||||||
solve
|
|
||||||
(
|
|
||||||
(
|
|
||||||
(max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
|
|
||||||
+ (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
|
|
||||||
)
|
|
||||||
+ p_rghEqnIncomp,
|
|
||||||
p_rgh.select(pimple.finalInnerIter())
|
|
||||||
);
|
|
||||||
|
|
||||||
if (pimple.finalNonOrthogonalIter())
|
|
||||||
{
|
|
||||||
p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
|
|
||||||
p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
|
|
||||||
|
|
||||||
dgdt =
|
|
||||||
(
|
|
||||||
pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
|
|
||||||
- pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
|
|
||||||
);
|
|
||||||
|
|
||||||
phi = phiHbyA + p_rghEqnIncomp.flux();
|
|
||||||
|
|
||||||
U = HbyA
|
|
||||||
+ rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
|
|
||||||
U.correctBoundaryConditions();
|
|
||||||
fvOptions.correct(U);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
{
|
|
||||||
Uf = fvc::interpolate(U);
|
|
||||||
surfaceVectorField n(mesh.Sf()/mesh.magSf());
|
|
||||||
Uf += n*(fvc::absolute(phi, U)/mesh.magSf() - (n & Uf));
|
|
||||||
}
|
|
||||||
|
|
||||||
// Update densities from change in p_rgh
|
|
||||||
mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
|
|
||||||
mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
|
|
||||||
|
|
||||||
rho = alpha1*rho1 + alpha2*rho2;
|
|
||||||
|
|
||||||
// Correct p_rgh for consistency with p and the updated densities
|
|
||||||
p = max(p_rgh + rho*gh, pMin);
|
|
||||||
p_rgh = p - rho*gh;
|
|
||||||
p_rgh.correctBoundaryConditions();
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
K = 0.5*magSqr(U);
|
|
||||||
}
|
|
||||||
@ -5,8 +5,8 @@
|
|||||||
\\ / A nd | www.openfoam.com
|
\\ / A nd | www.openfoam.com
|
||||||
\\/ M anipulation |
|
\\/ M anipulation |
|
||||||
-------------------------------------------------------------------------------
|
-------------------------------------------------------------------------------
|
||||||
Copyright (C) 2011-2017 OpenFOAM Foundation
|
Copyright (C) 2011-2018 OpenFOAM Foundation
|
||||||
Copyright (C) OpenCFD OpenCFD Ltd.
|
Copyright (C) 2024 OpenCFD Ltd.
|
||||||
-------------------------------------------------------------------------------
|
-------------------------------------------------------------------------------
|
||||||
License
|
License
|
||||||
This file is part of OpenFOAM.
|
This file is part of OpenFOAM.
|
||||||
@ -32,7 +32,9 @@ Group
|
|||||||
|
|
||||||
Description
|
Description
|
||||||
Solver for two compressible, non-isothermal immiscible fluids using a VOF
|
Solver for two compressible, non-isothermal immiscible fluids using a VOF
|
||||||
(volume of fluid) phase-fraction based interface capturing approach.
|
(volume of fluid) phase-fraction based interface capturing approach,
|
||||||
|
with optional mesh motion and mesh topology changes including adaptive
|
||||||
|
re-meshing.
|
||||||
|
|
||||||
The momentum and other fluid properties are of the "mixture" and a single
|
The momentum and other fluid properties are of the "mixture" and a single
|
||||||
momentum equation is solved.
|
momentum equation is solved.
|
||||||
@ -45,6 +47,7 @@ Description
|
|||||||
\*---------------------------------------------------------------------------*/
|
\*---------------------------------------------------------------------------*/
|
||||||
|
|
||||||
#include "fvCFD.H"
|
#include "fvCFD.H"
|
||||||
|
#include "dynamicFvMesh.H"
|
||||||
#include "CMULES.H"
|
#include "CMULES.H"
|
||||||
#include "EulerDdtScheme.H"
|
#include "EulerDdtScheme.H"
|
||||||
#include "localEulerDdtScheme.H"
|
#include "localEulerDdtScheme.H"
|
||||||
@ -53,6 +56,7 @@ Description
|
|||||||
#include "compressibleInterPhaseTransportModel.H"
|
#include "compressibleInterPhaseTransportModel.H"
|
||||||
#include "pimpleControl.H"
|
#include "pimpleControl.H"
|
||||||
#include "fvOptions.H"
|
#include "fvOptions.H"
|
||||||
|
#include "CorrectPhi.H"
|
||||||
#include "fvcSmooth.H"
|
#include "fvcSmooth.H"
|
||||||
|
|
||||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||||
@ -70,30 +74,31 @@ int main(int argc, char *argv[])
|
|||||||
#include "addCheckCaseOptions.H"
|
#include "addCheckCaseOptions.H"
|
||||||
#include "setRootCaseLists.H"
|
#include "setRootCaseLists.H"
|
||||||
#include "createTime.H"
|
#include "createTime.H"
|
||||||
#include "createMesh.H"
|
#include "createDynamicFvMesh.H"
|
||||||
#include "createControl.H"
|
#include "initContinuityErrs.H"
|
||||||
#include "createTimeControls.H"
|
#include "createDyMControls.H"
|
||||||
#include "createFields.H"
|
#include "createFields.H"
|
||||||
|
#include "CourantNo.H"
|
||||||
|
#include "setInitialDeltaT.H"
|
||||||
|
#include "createUfIfPresent.H"
|
||||||
|
|
||||||
volScalarField& p = mixture.p();
|
volScalarField& p = mixture.p();
|
||||||
volScalarField& T = mixture.T();
|
volScalarField& T = mixture.T();
|
||||||
const volScalarField& psi1 = mixture.thermo1().psi();
|
const volScalarField& psi1 = mixture.thermo1().psi();
|
||||||
const volScalarField& psi2 = mixture.thermo2().psi();
|
const volScalarField& psi2 = mixture.thermo2().psi();
|
||||||
|
|
||||||
if (!LTS)
|
|
||||||
{
|
|
||||||
#include "readTimeControls.H"
|
|
||||||
#include "CourantNo.H"
|
|
||||||
#include "setInitialDeltaT.H"
|
|
||||||
}
|
|
||||||
|
|
||||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||||
|
|
||||||
Info<< "\nStarting time loop\n" << endl;
|
Info<< "\nStarting time loop\n" << endl;
|
||||||
|
|
||||||
while (runTime.run())
|
while (runTime.run())
|
||||||
{
|
{
|
||||||
#include "readTimeControls.H"
|
#include "readDyMControls.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)));
|
||||||
|
|
||||||
if (LTS)
|
if (LTS)
|
||||||
{
|
{
|
||||||
@ -113,6 +118,44 @@ int main(int argc, char *argv[])
|
|||||||
// --- Pressure-velocity PIMPLE corrector loop
|
// --- Pressure-velocity PIMPLE corrector loop
|
||||||
while (pimple.loop())
|
while (pimple.loop())
|
||||||
{
|
{
|
||||||
|
if (pimple.firstIter() || moveMeshOuterCorrectors)
|
||||||
|
{
|
||||||
|
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
|
||||||
|
|
||||||
|
mesh.update();
|
||||||
|
|
||||||
|
if (mesh.changing())
|
||||||
|
{
|
||||||
|
MRF.update();
|
||||||
|
|
||||||
|
Info<< "Execution time for mesh.update() = "
|
||||||
|
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
|
||||||
|
<< " s" << endl;
|
||||||
|
|
||||||
|
gh = (g & mesh.C()) - ghRef;
|
||||||
|
ghf = (g & mesh.Cf()) - ghRef;
|
||||||
|
|
||||||
|
if (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 (checkMeshCourantNo)
|
||||||
|
{
|
||||||
|
#include "meshCourantNo.H"
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
#include "alphaControls.H"
|
#include "alphaControls.H"
|
||||||
#include "compressibleAlphaEqnSubCycle.H"
|
#include "compressibleAlphaEqnSubCycle.H"
|
||||||
|
|
||||||
|
|||||||
@ -8,6 +8,4 @@ CorrectPhi
|
|||||||
pimple
|
pimple
|
||||||
);
|
);
|
||||||
|
|
||||||
//***HGW phi.oldTime() = phi;
|
|
||||||
|
|
||||||
#include "continuityErrs.H"
|
#include "continuityErrs.H"
|
||||||
@ -6,7 +6,7 @@
|
|||||||
(
|
(
|
||||||
"phiHbyA",
|
"phiHbyA",
|
||||||
fvc::flux(HbyA)
|
fvc::flux(HbyA)
|
||||||
+ MRF.zeroFilter(fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi))
|
+ MRF.zeroFilter(fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi, Uf))
|
||||||
);
|
);
|
||||||
MRF.makeRelative(phiHbyA);
|
MRF.makeRelative(phiHbyA);
|
||||||
|
|
||||||
@ -23,6 +23,9 @@
|
|||||||
// Update the pressure BCs to ensure flux consistency
|
// Update the pressure BCs to ensure flux consistency
|
||||||
constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
|
constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
|
||||||
|
|
||||||
|
// Make the fluxes relative to the mesh motion
|
||||||
|
fvc::makeRelative(phiHbyA, U);
|
||||||
|
|
||||||
tmp<fvScalarMatrix> p_rghEqnComp1;
|
tmp<fvScalarMatrix> p_rghEqnComp1;
|
||||||
tmp<fvScalarMatrix> p_rghEqnComp2;
|
tmp<fvScalarMatrix> p_rghEqnComp2;
|
||||||
|
|
||||||
@ -34,8 +37,7 @@
|
|||||||
surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
|
surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
|
||||||
|
|
||||||
p_rghEqnComp1 =
|
p_rghEqnComp1 =
|
||||||
pos(alpha1)
|
(
|
||||||
*(
|
|
||||||
(
|
(
|
||||||
fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
|
fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
|
||||||
- (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
|
- (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
|
||||||
@ -48,11 +50,9 @@
|
|||||||
+ fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
|
+ fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
|
||||||
)
|
)
|
||||||
);
|
);
|
||||||
p_rghEqnComp1.ref().relax();
|
|
||||||
|
|
||||||
p_rghEqnComp2 =
|
p_rghEqnComp2 =
|
||||||
pos(alpha2)
|
(
|
||||||
*(
|
|
||||||
(
|
(
|
||||||
fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
|
fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
|
||||||
- (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
|
- (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
|
||||||
@ -65,7 +65,6 @@
|
|||||||
+ fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
|
+ fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
|
||||||
)
|
)
|
||||||
);
|
);
|
||||||
p_rghEqnComp2.ref().relax();
|
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
@ -94,6 +93,21 @@
|
|||||||
);
|
);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
if (mesh.moving())
|
||||||
|
{
|
||||||
|
p_rghEqnComp1.ref() += fvc::div(mesh.phi())*alpha1;
|
||||||
|
p_rghEqnComp2.ref() += fvc::div(mesh.phi())*alpha2;
|
||||||
|
}
|
||||||
|
|
||||||
|
p_rghEqnComp1.ref() *= pos(alpha1);
|
||||||
|
p_rghEqnComp2.ref() *= pos(alpha2);
|
||||||
|
|
||||||
|
if (pimple.transonic())
|
||||||
|
{
|
||||||
|
p_rghEqnComp1.ref().relax();
|
||||||
|
p_rghEqnComp2.ref().relax();
|
||||||
|
}
|
||||||
|
|
||||||
// Cache p_rgh prior to solve for density update
|
// Cache p_rgh prior to solve for density update
|
||||||
volScalarField p_rgh_0(p_rgh);
|
volScalarField p_rgh_0(p_rgh);
|
||||||
|
|
||||||
@ -131,6 +145,9 @@
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// Correct Uf if the mesh is moving
|
||||||
|
fvc::correctUf(Uf, U, fvc::absolute(phi, U));
|
||||||
|
|
||||||
// Update densities from change in p_rgh
|
// Update densities from change in p_rgh
|
||||||
mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
|
mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
|
||||||
mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
|
mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
|
||||||
|
|||||||
@ -14,7 +14,7 @@ FoamFile
|
|||||||
}
|
}
|
||||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||||
|
|
||||||
application compressibleInterDyMFoam;
|
application compressibleInterFoam;
|
||||||
|
|
||||||
startFrom startTime;
|
startFrom startTime;
|
||||||
|
|
||||||
@ -73,8 +73,8 @@ solvers
|
|||||||
|
|
||||||
U
|
U
|
||||||
{
|
{
|
||||||
solver smoothSolver;
|
solver PBiCGStab;
|
||||||
smoother GaussSeidel;
|
preconditioner DILU;
|
||||||
tolerance 1e-06;
|
tolerance 1e-06;
|
||||||
relTol 0;
|
relTol 0;
|
||||||
nSweeps 1;
|
nSweeps 1;
|
||||||
@ -82,8 +82,8 @@ solvers
|
|||||||
|
|
||||||
"(T|k|B|nuTilda).*"
|
"(T|k|B|nuTilda).*"
|
||||||
{
|
{
|
||||||
solver smoothSolver;
|
solver PBiCGStab;
|
||||||
smoother symGaussSeidel;
|
preconditioner DILU;
|
||||||
tolerance 1e-08;
|
tolerance 1e-08;
|
||||||
relTol 0;
|
relTol 0;
|
||||||
}
|
}
|
||||||
@ -15,7 +15,7 @@ FoamFile
|
|||||||
}
|
}
|
||||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||||
|
|
||||||
application compressibleInterDyMFoam;
|
application compressibleInterFoam;
|
||||||
|
|
||||||
startFrom startTime;
|
startFrom startTime;
|
||||||
|
|
||||||
@ -67,8 +67,8 @@ solvers
|
|||||||
|
|
||||||
"(U|k|epsilon|T)"
|
"(U|k|epsilon|T)"
|
||||||
{
|
{
|
||||||
solver smoothSolver;
|
solver PBiCGStab;
|
||||||
smoother symGaussSeidel;
|
preconditioner DILU;
|
||||||
tolerance 1e-08;
|
tolerance 1e-08;
|
||||||
relTol 0.1;
|
relTol 0.1;
|
||||||
}
|
}
|
||||||
Reference in New Issue
Block a user