Commit d141b3c9 authored by sergio's avatar sergio
Browse files

ENH: Adding interCondensingEvaporatingFoam and tutorial

parent e8d73e55
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
set -x
wclean libso temperaturePhaseChangeTwoPhaseMixtures
wclean
#------------------------------------------------------------------------------
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
# Parse arguments for library compilation
targetType=libso
. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
set -x
wmake $targetType temperaturePhaseChangeTwoPhaseMixtures
wmake
#------------------------------------------------------------------------------
interCondensatingEvaporatingFoam.C
EXE = $(FOAM_APPBIN)/interCondensatingEvaporatingFoam
interPhaseChangePath = $(FOAM_SOLVERS)/multiphase/interPhaseChangeFoam
EXE_INC = \
-ItemperaturePhaseChangeTwoPhaseMixtures/lnInclude \
-I$(interPhaseChangePath) \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-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)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude\
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lphaseTemperatureChangeTwoPhaseMixtures \
-ltwoPhaseMixture \
-linterfaceProperties \
-ltwoPhaseProperties \
-lincompressibleTransportModels \
-lturbulenceModels \
-lincompressibleTurbulenceModels \
-lfiniteVolume \
-lmeshTools \
-lfvOptions \
-lsampling \
-lfluidThermophysicalModels
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
+ turbulence->divDevRhoReff(rho, U)
);
UEqn.relax();
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
interface.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
}
volScalarField contErr(fvc::ddt(rho) + fvc::div(rhoPhi));
scalar sumLocalContErr = runTime.deltaTValue()*
mag(contErr)().weightedAverage(mesh.V()).value();
scalar globalContErr = runTime.deltaTValue()*
contErr.weightedAverage(mesh.V()).value();
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< endl;
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"
// Create p before the thermo
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh
);
// Creating e based thermo
autoPtr<twoPhaseMixtureEThermo> thermo;
thermo.set(new twoPhaseMixtureEThermo(U, phi));
// Create mixture and
Info<< "Creating temperaturePhaseChangeTwoPhaseMixture\n" << endl;
autoPtr<temperaturePhaseChangeTwoPhaseMixture> mixture =
temperaturePhaseChangeTwoPhaseMixture::New(thermo(), mesh);
volScalarField& T = thermo->T();
volScalarField& e = thermo->he();
// Correct e from T and alpha
thermo->correct();
volScalarField& alpha1(thermo->alpha1());
volScalarField& alpha2(thermo->alpha2());
const dimensionedScalar& rho1 = thermo->rho1();
const dimensionedScalar& rho2 = thermo->rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
alpha1*rho1 + alpha2*rho2,
alpha1.boundaryField().types()
);
rho.oldTime();
// Construct interface from alpha1 distribution
interfaceProperties interface
(
alpha1,
U,
thermo->transportPropertiesDict()
);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, thermo())
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
//Update p with rho
p = p_rgh + rho*gh;
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
// Turbulent Prandtl number
dimensionedScalar Prt("Prt", dimless, thermo->transportPropertiesDict());
volScalarField kappaEff
(
IOobject
(
"kappaEff",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
thermo->kappa()
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
Info<< "Creating field pDivU\n" << endl;
volScalarField pDivU
(
IOobject
(
"pDivU",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("pDivU", p.dimensions()/dimTime, 0)
);
{
tmp<volScalarField> tcp(thermo->Cp());
const volScalarField& cp = tcp();
kappaEff = thermo->kappa() + rho*cp*turbulence->nut()/Prt;
pDivU = dimensionedScalar("pDivU", p.dimensions()/dimTime, 0.0);
if (thermo->pDivU())
{
pDivU = (p*fvc::div(rhoPhi/fvc::interpolate(rho)));
}
fvScalarMatrix eEqn
(
fvm::ddt(rho, e)
+ fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
+ fvm::div(rhoPhi, e)
- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), e)
- fvm::laplacian(kappaEff/cp, e)
+ pDivU
);
eEqn.relax();
eEqn.solve();
thermo->correct();
Info<< "min/max(T) = " << min(T).value() << ", "
<< max(T).value() <<endl;
}
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 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 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
interCondensatingEvaporatingFoam
Group
grpMultiphaseSolvers
Description
Solver for 2 incompressible, non-isothermal immiscible fluids with
phase-change (evaporation-condensation) between a fluid and its vapour.
Uses a VOF (volume of fluid) phase-fraction based interface capturing
approach.
The momentum, energy 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 "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "twoPhaseMixtureEThermo.H"
#include "temperaturePhaseChangeTwoPhaseMixture.H"
#include "turbulentTransportModel.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
pimpleControl pimple(mesh);
#include "readGravitationalAcceleration.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "createTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "alphaControls.H"
surfaceScalarField rhoPhi
(
IOobject
(
"rhoPhi",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", dimMass/dimTime, 0)
);
#include "alphaEqnSubCycle.H"
mixture->correct();
#include "UEqn.H"
#include "eEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
#include "continuityError.H"
if (pimple.turbCorr())
{
turbulence->correct();
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
{
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
surfaceScalarField phig
(
(
interface.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);
phiHbyA += phig;
// Update the fixedFluxPressure BCs to ensure flux consistency
setSnGrad<fixedFluxPressureFvPatchScalarField>
(
p_rgh.boundaryField(),
(
phiHbyA.boundaryField()
- (mesh.Sf().boundaryField() & U.boundaryField())
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);
Pair<tmp<volScalarField> > vDot = mixture->vDot();
const volScalarField& vDotc = vDot[0]();
const volScalarField& vDotv = vDot[1]();
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rAUf, p_rgh)
- (vDotc - vDotv)
);
p_rghEqn.setReference(pRefCell, pRefValue);
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + p_rghEqn.flux();
U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
}
}
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}
temperaturePhaseChangeTwoPhaseMixtures/newtemperaturePhaseChangeTwoPhaseMixture.C
temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.C
thermoIncompressibleTwoPhaseMixture/thermoIncompressibleTwoPhaseMixture.C
twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.C
constant/constant.C
LIB = $(FOAM_LIBBIN)/libphaseTemperatureChangeTwoPhaseMixtures
EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
LIB_LIBS = \
-ltwoPhaseMixture \
-linterfaceProperties \
-ltwoPhaseProperties \
-lincompressibleTransportModels \
-lfiniteVolume \
-lfluidThermophysicalModels
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 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 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/>.
\*---------------------------------------------------------------------------*/
#include "constant.H"
#include "addToRunTimeSelectionTable.H"
#include "fvcGrad.H"
#include "twoPhaseMixtureEThermo.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace temperaturePhaseChangeTwoPhaseMixtures
{
defineTypeNameAndDebug(constant, 0);
addToRunTimeSelectionTable
(
temperaturePhaseChangeTwoPhaseMixture,
constant,
components
);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::temperaturePhaseChangeTwoPhaseMixtures::