Compare commits

...

1 Commits

Author SHA1 Message Date
Kutalmis Bercin
54dcde5f9a INT: compressibleInterFoam: assimilate DyM solver
Based on changes from openfoam.org
2024-05-01 16:44:42 +01:00
58 changed files with 96 additions and 452 deletions

View File

@ -7,7 +7,6 @@ wclean libso surfaceTensionModels
wclean libso VoFphaseCompressibleTurbulenceModels
wclean
wclean compressibleInterDyMFoam
wclean compressibleInterFilmFoam
wclean compressibleInterIsoFoam

View File

@ -8,7 +8,6 @@ wmake $targetType surfaceTensionModels
wmake $targetType VoFphaseCompressibleTurbulenceModels
wmake $targetType
wmake $targetType compressibleInterDyMFoam
wmake $targetType compressibleInterFilmFoam
wmake $targetType compressibleInterIsoFoam
wmake $targetType overCompressibleInterDyMFoam

View File

@ -1,21 +1,25 @@
EXE_INC = \
-I../VoF \
-ItwoPhaseMixtureThermo \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-IVoFphaseCompressibleTurbulenceModels/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)/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 \
-IVoFphaseCompressibleTurbulenceModels/lnInclude
-I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lfvOptions \
-lmeshTools \
-ldynamicMesh \
-ldynamicFvMesh \
-ltwoPhaseMixtureThermo \
-ltwoPhaseSurfaceTension \
-lcompressibleTransportModels \

View File

@ -1,3 +0,0 @@
compressibleInterDyMFoam.C
EXE = $(FOAM_APPBIN)/compressibleInterDyMFoam

View File

@ -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

View File

@ -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())
);

View File

@ -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;
}
// ************************************************************************* //

View File

@ -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);
}

View File

@ -5,8 +5,8 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) OpenCFD OpenCFD Ltd.
Copyright (C) 2011-2018 OpenFOAM Foundation
Copyright (C) 2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -32,7 +32,9 @@ Group
Description
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
momentum equation is solved.
@ -45,6 +47,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
@ -53,6 +56,7 @@ Description
#include "compressibleInterPhaseTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -70,30 +74,31 @@ int main(int argc, char *argv[])
#include "addCheckCaseOptions.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
#include "createUfIfPresent.H"
volScalarField& p = mixture.p();
volScalarField& T = mixture.T();
const volScalarField& psi1 = mixture.thermo1().psi();
const volScalarField& psi2 = mixture.thermo2().psi();
if (!LTS)
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
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)
{
@ -113,6 +118,44 @@ int main(int argc, char *argv[])
// --- 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 (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 "compressibleAlphaEqnSubCycle.H"

View File

@ -8,6 +8,4 @@ CorrectPhi
pimple
);
//***HGW phi.oldTime() = phi;
#include "continuityErrs.H"

View File

@ -6,7 +6,7 @@
(
"phiHbyA",
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);
@ -23,6 +23,9 @@
// 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;
@ -34,8 +37,7 @@
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)
@ -48,11 +50,9 @@
+ 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)
@ -65,7 +65,6 @@
+ fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
)
);
p_rghEqnComp2.ref().relax();
}
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
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
mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));

View File

@ -14,7 +14,7 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application compressibleInterDyMFoam;
application compressibleInterFoam;
startFrom startTime;

View File

@ -73,8 +73,8 @@ solvers
U
{
solver smoothSolver;
smoother GaussSeidel;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-06;
relTol 0;
nSweeps 1;
@ -82,8 +82,8 @@ solvers
"(T|k|B|nuTilda).*"
{
solver smoothSolver;
smoother symGaussSeidel;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-08;
relTol 0;
}

View File

@ -15,7 +15,7 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application compressibleInterDyMFoam;
application compressibleInterFoam;
startFrom startTime;

View File

@ -67,8 +67,8 @@ solvers
"(U|k|epsilon|T)"
{
solver smoothSolver;
smoother symGaussSeidel;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-08;
relTol 0.1;
}