Commit 91e7cab0 authored by henry's avatar henry
Browse files

Add interPhaseChangeFoam solver

parent 23fb7670
interPhaseChangeFoam.C
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/newPhaseChangeTwoPhaseMixture.C
phaseChangeTwoPhaseMixtures/Kunz/Kunz.C
phaseChangeTwoPhaseMixtures/Merkle/Merkle.C
phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C
EXE = $(FOAM_APPBIN)/interPhaseChangeFoam
EXE_INC = \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/LESmodels \
-I$(LIB_SRC)/LESmodels/LESdeltas/lnInclude \
-IphaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleLESmodels \
-lfiniteVolume
surfaceScalarField muf =
twoPhaseProperties->muf()
+ fvc::interpolate(rho*turbulence->nuSgs());
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
- fvm::laplacian(muf, U)
- (fvc::grad(U) & fvc::grad(muf))
//- fvc::div(muf*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))
);
UEqn.relax();
if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(pd)
) * mesh.magSf()
)
);
}
{
# include "continuityErrs.H"
wordList pcorrTypes(pd.boundaryField().types());
for (label i=0; i<pd.boundaryField().size(); i++)
{
if (pd.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", pd.dimensions(), 0.0),
pcorrTypes
);
dimensionedScalar rUAf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
adjustPhi(phi, U, pcorr);
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pdRefCell, pdRefValue);
pcorrEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pcorrEqn.flux();
}
}
# include "continuityErrs.H"
}
Info<< "Reading field pd\n" << endl;
volScalarField pd
(
IOobject
(
"pd",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field gamma\n" << endl;
volScalarField gamma
(
IOobject
(
"gamma",
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<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr<phaseChangeTwoPhaseMixture> twoPhaseProperties =
phaseChangeTwoPhaseMixture::New(U, phi, "gamma");
const dimensionedScalar& rho1 = twoPhaseProperties->rho1();
const dimensionedScalar& rho2 = twoPhaseProperties->rho2();
const dimensionedScalar& pSat = twoPhaseProperties->pSat();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
gamma*rho1 + (scalar(1) - gamma)*rho2,
gamma.boundaryField().types()
);
rho.oldTime();
label pdRefCell = 0;
scalar pdRefValue = 0.0;
setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
Info<< "Calculating field g.h" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pd + rho*gh
);
// Construct interface from gamma distribution
interfaceProperties interface(gamma, U, twoPhaseProperties());
// Construct LES model
autoPtr<LESmodel> turbulence
(
LESmodel::New(U, phi, twoPhaseProperties())
);
{
word gammaScheme("div(phi,gamma)");
word gammarScheme("div(phirb,gamma)");
surfaceScalarField phir("phir", phic*interface.nHatf());
for (int gCorr=0; gCorr<nGammaCorr; gCorr++)
{
surfaceScalarField phiGamma =
fvc::flux
(
phi,
gamma,
gammaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, scalar(1) - gamma, gammarScheme),
gamma,
gammarScheme
);
Pair<tmp<volScalarField> > vDotAlphal =
twoPhaseProperties->vDotAlphal();
const volScalarField& vDotcAlphal = vDotAlphal[0]();
const volScalarField& vDotvAlphal = vDotAlphal[1]();
volScalarField Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
vDotvAlphal - vDotcAlphal
);
volScalarField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
divU*gamma
+ vDotcAlphal
);
//MULES::explicitSolve(gamma, phi, phiGamma, 1, 0);
//MULES::explicitSolve(oneField(), gamma, phi, phiGamma, Sp, Su, 1, 0);
MULES::implicitSolve(oneField(), gamma, phi, phiGamma, Sp, Su, 1, 0);
rhoPhi +=
(runTime.deltaT()/totalDeltaT)
*(phiGamma*(rho1 - rho2) + phi*rho2);
}
Info<< "Liquid phase volume fraction = "
<< gamma.weightedAverage(mesh.V()).value()
<< " Min(gamma) = " << min(gamma).value()
<< " Max(gamma) = " << max(gamma).value()
<< endl;
}
surfaceScalarField rhoPhi
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0)
);
{
label nGammaCorr
(
readLabel(piso.lookup("nGammaCorr"))
);
label nGammaSubCycles
(
readLabel(piso.lookup("nGammaSubCycles"))
);
surfaceScalarField phic = mag(phi/mesh.magSf());
phic = min(interface.cGamma()*phic, max(phic));
volScalarField divU = fvc::div(phi);
dimensionedScalar totalDeltaT = runTime.deltaT();
if (nGammaSubCycles > 1)
{
for
(
subCycle<volScalarField> gammaSubCycle(gamma, nGammaSubCycles);
!(++gammaSubCycle).end();
)
{
# include "gammaEqn.H"
}
}
else
{
# include "gammaEqn.H"
}
if (nOuterCorr == 1)
{
interface.correct();
}
rho == gamma*rho1 + (scalar(1) - gamma)*rho2;
}
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
interPhaseChangeFoam
Description
Solver for 2 incompressible, isothermal immiscible fluids with phase-change.
Uses a VOF (volume of fluid) phase-fraction based interface capturing
approach. The momentum and other fluid properties are of the "mixture" and
a single momentum equation is solved.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "phaseChangeTwoPhaseMixture.H"
#include "incompressible/LESmodel/LESmodel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "readEnvironmentalProperties.H"
# include "readPISOControls.H"
# include "initContinuityErrs.H"
# include "createFields.H"
# include "readTimeControls.H"
# include "correctPhi.H"
# include "CourantNo.H"
# include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
# include "readPISOControls.H"
# include "readTimeControls.H"
# include "CourantNo.H"
# include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
twoPhaseProperties->correct();
# include "gammaEqnSubCycle.H"
turbulence->correct();
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}
#include "continuityErrs.H"
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);
U = rUA*UEqn.H();
surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);
phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();
adjustPhi(phi, U, pd);
Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
const volScalarField& vDotcP = vDotP[0]();
const volScalarField& vDotvP = vDotP[1]();
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvc::div(phi) - fvm::laplacian(rUAf, pd)
+ (vDotvP - vDotcP)*(rho*gh - pSat) + fvm::Sp(vDotvP - vDotcP, pd)
);
pdEqn.setReference(pdRefCell, pdRefValue);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pdEqn.solve(mesh.solver(pd.name() + "Final"));
}
else
{
pdEqn.solve(mesh.solver(pd.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi += pdEqn.flux();
}
}
p = pd + rho*gh;
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
}
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "Kunz.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace phaseChangeTwoPhaseMixtures
{
defineTypeNameAndDebug(Kunz, 0);
addToRunTimeSelectionTable(phaseChangeTwoPhaseMixture, Kunz, components);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::phaseChangeTwoPhaseMixtures::Kunz::Kunz
(
const volVectorField& U,
const surfaceScalarField& phi,
const word& alpha1Name
)
:
phaseChangeTwoPhaseMixture(typeName, U, phi, alpha1Name),
UInf_(phaseChangeTwoPhaseMixtureCoeffs_.lookup("UInf")),
tInf_(phaseChangeTwoPhaseMixtureCoeffs_.lookup("tInf")),
Cc_(phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cc")),
Cv_(phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cv")),
p0_("0", pSat().dimensions(), 0.0),
mcCoeff_(Cc_*rho2()/tInf_),
mvCoeff_(Cv_*rho2()/(0.5*rho1()*sqr(UInf_)*tInf_))
{
correct();
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::Pair<Foam::tmp<Foam::volScalarField> >
Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotAlphal() const
{
const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
volScalarField limitedAlpha1 = min(max(alpha1_, scalar(0)), scalar(1));
return Pair<tmp<volScalarField> >
(
mcCoeff_*sqr(limitedAlpha1)
*max(p - pSat(), p0_)/max(p - pSat(), 0.01*pSat()),
mvCoeff_*min(p - pSat(), p0_)
);
}
Foam::Pair<Foam::tmp<Foam::volScalarField> >
Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotP() const
{
const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
volScalarField limitedAlpha1 = min(max(alpha1_, scalar(0)), scalar(1));
return Pair<tmp<volScalarField> >
(
mcCoeff_*sqr(limitedAlpha1)*(1.0 - limitedAlpha1)
*pos(p - pSat())/max(p - pSat(), 0.01*pSat()),