zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

CONT: Addition of compressibleIsoInterFOam and PLIC

Browse Source 
1) Implementation of the compressibleIsoInterFOam solver
   2) Implementation of a new PLIC interpolation scheme.
   3) New tutorials associated with the solvers

This implementation was carried out by Henning Scheufler (DLR) and Johan
Roenby (DHI), following :

\verbatim

Henning Scheufler, Johan Roenby,
Accurate and efficient surface reconstruction from volume fraction data
on general meshes, Journal of Computational Physics, 2019, doi
10.1016/j.jcp.2019.01.009

\endverbatim

The integration of the code was carried out by Andy Heather and Sergio
Ferraris from OpenCFD Ltd.
This commit is contained in:
Henning Scheufler
2020-06-05 14:31:15 +02:00
committed by Andrew Heather
parent 237f2e1076
commit 44a84d4778
166 changed files with 16335 additions and 3372 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/solvers/multiphase/compressibleInterIsoFoam/Make/files Normal file
View File

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

39
applications/solvers/multiphase/compressibleInterIsoFoam/Make/options Normal file
View File

@ -0,0 +1,39 @@
EXE_INC = \
-I. \
-I../VoF \
-I../compressibleInterFoam \
-I$(FOAM_SOLVERS)/multiphase/compressibleInterFoam/twoPhaseMixtureThermo \
-I$(FOAM_SOLVERS)/multiphase/compressibleInterFoam/VoFphaseCompressibleTurbulenceModels \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/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)/transportModels/geometricVoF/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lfvOptions \
-lsurfMesh \
-lmeshTools \
-ldynamicMesh \
-ldynamicFvMesh \
-ltwoPhaseMixtureThermo \
-ltwoPhaseSurfaceTension \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lspecie \
-ltwoPhaseMixture \
-ltwoPhaseProperties \
-linterfaceProperties \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-lVoFphaseCompressibleTurbulenceModels \
-lgeometricVoF

3
applications/solvers/multiphase/compressibleInterIsoFoam/alphaControls.H Normal file
View File

@ -0,0 +1,3 @@
const dictionary& alphaControls = mesh.solverDict(alpha1.name());
label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));

15
applications/solvers/multiphase/compressibleInterIsoFoam/alphaEqn.H Normal file
View File

@ -0,0 +1,15 @@
// Update alpha1
#include "alphaSuSp.H"
advector.advect(Sp,(Su + divU*min(alpha1(), scalar(1)))());
// Update rhoPhi
rhoPhi = advector.getRhoPhi(rho1, rho2);
alphaPhi10 = advector.alphaPhi();
alpha2 = 1.0 - alpha1;
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1
<< endl;

43
applications/solvers/multiphase/compressibleInterIsoFoam/alphaSuSp.H Normal file
View File

@ -0,0 +1,43 @@
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)
{
Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
}
else if (dgdt[celli] < 0.0)
{
Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
}
}
volScalarField::Internal divU
(
mesh.moving()
? fvc::div(phi + mesh.phi())
: fvc::div(phi)
);

70
applications/solvers/multiphase/compressibleInterIsoFoam/compressibleAlphaEqnSubCycle.H Normal file
View File

@ -0,0 +1,70 @@
if (pimple.nCorrPIMPLE() > 1)
{
if (!pimple.firstIter())
{
// Resetting alpha1 to value before advection in first PIMPLE
// iteration.
alpha1 = alpha1.oldTime();
}
}
tmp<surfaceScalarField> talphaPhi1(alphaPhi10);
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
talphaPhi1 = new surfaceScalarField
(
IOobject
(
"alphaPhi1",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(alphaPhi10.dimensions(), Zero)
);
surfaceScalarField rhoPhiSum
(
IOobject
(
"rhoPhiSum",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(rhoPhi.dimensions(), Zero)
);
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
#include "alphaEqn.H"
talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi10;
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
rhoPhi = rhoPhiSum;
}
else
{
#include "alphaEqn.H"
}
rho == alpha1*rho1 + alpha2*rho2;
const surfaceScalarField& alphaPhi1 = talphaPhi1();
surfaceScalarField alphaPhi2("alphaPhi2", phi - alphaPhi1);
volScalarField::Internal contErr
(
(
fvc::ddt(rho) + fvc::div(rhoPhi)
- (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
- (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
)()
);

203
applications/solvers/multiphase/compressibleInterIsoFoam/compressibleInterIsoFoam.C Normal file
View File

@ -0,0 +1,203 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020 Johan Roenby
Copyright (C) 2020 DLR
-------------------------------------------------------------------------------
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
compressibleInterFlow
Description
Solver derived from interFoam for two compressible, immiscible
fluids using the isoAdvector phase-fraction based interface capturing
approach, with optional mesh motion and mesh topology changes including
adaptive re-meshing.
Reference:
\verbatim
Roenby, J., Bredmose, H. and Jasak, H. (2016).
A computational method for sharp interface advection
Royal Society Open Science, 3
doi 10.1098/rsos.160405
\endverbatim
\*---------------------------------------------------------------------------*/
#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"
#include "dynamicRefineFvMesh.H"
#include "isoAdvection.H"
#include "twoPhaseMixtureThermo.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)));
#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();
if (isA<dynamicRefineFvMesh>(mesh))
{
advector.surf().reconstruct();
}
mesh.update();
if (mesh.changing())
{
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
if (isA<dynamicRefineFvMesh>(mesh))
{
advector.surf().mapAlphaField();
alpha2 = 1.0 - alpha1;
alpha2.correctBoundaryConditions();
rho == alpha1*rho1 + alpha2*rho2;
rho.correctBoundaryConditions();
rho.oldTime() = rho;
alpha2.oldTime() = alpha2;
}
MRF.update();
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
}
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"
volScalarField divUp("divUp", fvc::div(fvc::absolute(phi, U), p));
#include "TEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence.correct();
}
}
runTime.write();
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

13
applications/solvers/multiphase/compressibleInterIsoFoam/correctPhi.H Normal file
View File

@ -0,0 +1,13 @@
CorrectPhi
(
U,
phi,
p,
dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
divU,
pimple
);
//***HGW phi.oldTime() = phi;
#include "continuityErrs.H"

110
applications/solvers/multiphase/compressibleInterIsoFoam/createFields.H Normal file
View File

@ -0,0 +1,110 @@
#include "createRDeltaT.H"
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"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
);
#include "createPhi.H"
Info<< "Constructing twoPhaseMixtureThermo\n" << endl;
twoPhaseMixtureThermo mixture(U, phi);
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
Info<< "Reading thermophysical properties\n" << endl;
const volScalarField& rho1 = mixture.thermo1().rho();
const volScalarField& rho2 = mixture.thermo2().rho();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
alpha1*rho1 + alpha2*rho2
);
dimensionedScalar pMin
(
"pMin",
dimPressure,
mixture
);
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha1 solution before it is used in the U equation.
surfaceScalarField rhoPhi
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
fvc::interpolate(rho)*phi
);
volScalarField dgdt(alpha1*fvc::div(phi));
#include "createAlphaFluxes.H"
Foam::isoAdvection advector(alpha1,phi,U);
// Construct compressible turbulence model
compressibleInterPhaseTransportModel turbulence
(
rho,
U,
phi,
rhoPhi,
alphaPhi10,
mixture
);
#include "createK.H"
#include "createMRF.H"
#include "createFvOptions.H"

169
applications/solvers/multiphase/compressibleInterIsoFoam/pEqn.H Normal file
View File

@ -0,0 +1,169 @@
{
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 =
(
(
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 =
(
(
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
{
#include "rhofs.H"
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*psi1/rho1)*correction(fvm::ddt(p_rgh))
);
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*psi2/rho2)*correction(fvm::ddt(p_rgh))
);
}
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);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqnIncomp
(
fvc::div(phiHbyA)
- fvm::laplacian(rAUf, p_rgh)
);
solve
(
p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
mesh.solver(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 =
(
alpha1*(p_rghEqnComp2 & p_rgh)
- alpha2*(p_rghEqnComp1 & p_rgh)
);
phi = phiHbyA + p_rghEqnIncomp.flux();
U = HbyA
+ rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}
// Correct Uf if the mesh is moving
{
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_rgh = p - rho*gh;
p_rgh.correctBoundaryConditions();
K = 0.5*magSqr(U);
}

6
applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Make/options
View File

@ -7,7 +7,8 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/transportModels/geometricVoF/lnInclude
LIB_LIBS = \
-lfiniteVolume \
@ -15,4 +16,5 @@ LIB_LIBS = \
-lfluidThermophysicalModels \
-lreactionThermophysicalModels \
-lsolidThermo \
-lsolidSpecie
-lsolidSpecie \
-lgeometricVoF

6
applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceHeatResistance/interfaceHeatResistance.C
View File

@ -28,7 +28,7 @@ License
#include "interfaceHeatResistance.H"
#include "constants.H"
#include "isoCutCell.H"
#include "cutCellIso.H"
#include "volPointInterpolation.H"
#include "wallPolyPatch.H"
#include "fvcSmooth.H"
@ -52,7 +52,7 @@ interfaceHeatResistance<Thermo, OtherThermo>
volPointInterpolation::New(mesh).interpolate(alpha)
);
isoCutCell cutCell(mesh, ap);
cutCellIso cutCell(mesh, ap);
forAll(interfaceArea_, celli)
{
@ -61,7 +61,7 @@ interfaceHeatResistance<Thermo, OtherThermo>
if (status == 0) // cell is cut
{
interfaceArea_[celli] =
mag(cutCell.isoFaceArea())/mesh.V()[celli];
mag(cutCell.faceArea())/mesh.V()[celli];
}
}

6
applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C
View File

@ -27,7 +27,7 @@ License
#include "kineticGasEvaporation.H"
#include "constants.H"
#include "isoCutCell.H"
#include "cutCellIso.H"
#include "volPointInterpolation.H"
#include "wallPolyPatch.H"
#include "fvcSmooth.H"
@ -50,7 +50,7 @@ void Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
volPointInterpolation::New(mesh).interpolate(alpha)
);
isoCutCell cutCell(mesh, ap);
cutCellIso cutCell(mesh, ap);
forAll(interfaceArea_, celli)
{
@ -59,7 +59,7 @@ void Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
if (status == 0) // cell is cut
{
interfaceArea_[celli] =
mag(cutCell.isoFaceArea())/mesh.V()[celli];
mag(cutCell.faceArea())/mesh.V()[celli];
}
}

11
applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/Make/options
View File

@ -1,16 +1,17 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/geometricVoF/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
-I$(LIB_SRC)/transportModels/incompressible/lnInclude
LIB_LIBS = \
-lfiniteVolume \
-lgeometricVoF \
-ltwoPhaseMixture \
-linterfaceProperties \
-ltwoPhaseProperties \
-lincompressibleTransportModels \
-lfiniteVolume \
-lfluidThermophysicalModels

6
applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C
View File

@ -29,7 +29,7 @@ License
#include "interfaceHeatResistance.H"
#include "addToRunTimeSelectionTable.H"
#include "twoPhaseMixtureEThermo.H"
#include "isoCutCell.H"
#include "cutCellIso.H"
#include "volPointInterpolation.H"
#include "calculatedFvPatchFields.H"
#include "wallPolyPatch.H"
@ -342,7 +342,7 @@ updateInterface()
volPointInterpolation::New(mesh_).interpolate(mixture_.alpha1())
);
isoCutCell cutCell(mesh_, ap);
cutCellIso cutCell(mesh_, ap);
forAll(interfaceArea_, celli)
{
@ -351,7 +351,7 @@ updateInterface()
if (status == 0) // cell is cut
{
interfaceArea_[celli] =
mag(cutCell.isoFaceArea())/mesh_.V()[celli];
mag(cutCell.faceArea())/mesh_.V()[celli];
}
}

12
applications/solvers/multiphase/interIsoFoam/Make/options
View File

@ -1,19 +1,26 @@
EXE_INC = \
-I../VoF \
-I../interFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/geometricVoF/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude
-I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/geometricVoF/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lfvOptions \
-lsurfMesh \
-lmeshTools \
-lsampling \
-ldynamicFvMesh \
@ -22,4 +29,5 @@ EXE_LIBS = \
-limmiscibleIncompressibleTwoPhaseMixture \
-lturbulenceModels \
-lincompressibleTurbulenceModels \
-lwaveModels
-lwaveModels \
-lgeometricVoF

33
applications/solvers/multiphase/interIsoFoam/UEqn.H
View File

@ -1,33 +0,0 @@
MRF.correctBoundaryVelocity(U);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(rho, U)
==
fvOptions(rho, U)
);
UEqn.relax();
fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
fvOptions.correct(U);
}

59
applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H
View File

@ -1,59 +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
-------------------------------------------------------------------------------
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/>.
Global
alphaCourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar maxAlphaCo
(
runTime.controlDict().get<scalar>("maxAlphaCo")
);
scalar alphaCoNum = 0.0;
scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
scalarField sumPhi
(
mixture.nearInterface()().primitiveField()
*fvc::surfaceSum(mag(phi))().primitiveField()
);
alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
meanAlphaCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
}
Info<< "Interface Courant Number mean: " << meanAlphaCoNum
<< " max: " << alphaCoNum << endl;
// ************************************************************************* //

3
applications/solvers/multiphase/interIsoFoam/alphaEqn.H
View File

@ -6,7 +6,8 @@
}
// Updating alpha1
advector.advect();
#include "alphaSuSp.H"
advector.advect(Sp, Su);
// Making U absolute again after advection step
if (mesh.moving())

11
applications/solvers/multiphase/interIsoFoam/correctPhi.H
View File

@ -1,11 +0,0 @@
CorrectPhi
(
U,
phi,
p_rgh,
surfaceScalarField("rAUf", fvc::interpolate(rAU())),
geometricZeroField(),
pimple
);
#include "continuityErrs.H"

34
applications/solvers/multiphase/interIsoFoam/initCorrectPhi.H
View File

@ -1,34 +0,0 @@
tmp<volScalarField> rAU;
if (correctPhi)
{
rAU = new volScalarField
(
IOobject
(
"rAU",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("rAU", dimTime/dimDensity, 1)
);
#include "correctPhi.H"
}
else
{
CorrectPhi
(
U,
phi,
p_rgh,
dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
geometricZeroField(),
pimple
);
#include "continuityErrs.H"
}

19
applications/solvers/multiphase/interIsoFoam/interIsoFoam.C
View File

@ -9,6 +9,7 @@
Copyright (C) 2016 DHI
Copyright (C) 2017 OpenCFD Ltd.
Copyright (C) 2018 Johan Roenby
Copyright (C) 2019-2020 DLR
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -63,6 +64,7 @@ Description
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "fvcSmooth.H"
#include "dynamicRefineFvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -113,14 +115,29 @@ int main(int argc, char *argv[])
{
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
if (isA<dynamicRefineFvMesh>(mesh))
{
advector.surf().reconstruct();
}
mesh.update();
if (mesh.changing())
{
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
if (isA<dynamicRefineFvMesh>(mesh))
{
advector.surf().mapAlphaField();
alpha2 = 1.0 - alpha1;
alpha2.correctBoundaryConditions();
rho == alpha1*rho1 + alpha2*rho2;
rho.correctBoundaryConditions();
rho.oldTime() = rho;
alpha2.oldTime() = alpha2;
}
MRF.update();
if (correctPhi)

89
applications/solvers/multiphase/interIsoFoam/pEqn.H
View File

@ -1,89 +0,0 @@
{
if (correctPhi)
{
rAU.ref() = 1.0/UEqn.A();
}
else
{
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, phi, Uf))
);
MRF.makeRelative(phiHbyA);
if (p_rgh.needReference())
{
fvc::makeRelative(phiHbyA, U);
adjustPhi(phiHbyA, U, p_rgh);
fvc::makeAbsolute(phiHbyA, U);
}
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);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA - p_rghEqn.flux();
p_rgh.relax();
U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}
#include "continuityErrs.H"
// Correct Uf if the mesh is moving
fvc::correctUf(Uf, U, phi);
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
if (!correctPhi)
{
rAU.clear();
}
}

2
applications/solvers/multiphase/interIsoFoam/rhofs.H
View File

@ -1,2 +0,0 @@
const dimensionedScalar& rho1f(rho1);
const dimensionedScalar& rho2f(rho2);

55
applications/solvers/multiphase/interIsoFoam/setDeltaT.H
View File

@ -1,55 +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
-------------------------------------------------------------------------------
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/>.
Global
setDeltaT
Description
Reset the timestep to maintain a constant maximum courant Number.
Reduction of time-step is immediate, but increase is damped to avoid
unstable oscillations.
\*---------------------------------------------------------------------------*/
if (adjustTimeStep)
{
scalar maxDeltaTFact =
min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
runTime.setDeltaT
(
min
(
deltaTFact*runTime.deltaTValue(),
maxDeltaT
)
);
Info<< "deltaT = " << runTime.deltaTValue() << endl;
}
// ************************************************************************* //