Commit 175f03f8 authored by Henry's avatar Henry
Browse files

Removed twoPhaseEulerFoam and renamed compressibleTwoPhaseEulerFoam -> twoPhaseEulerFoam

parent d54c9c73
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wclean libso twoPhaseSystem
wclean libso interfacialModels
wclean libso phaseIncompressibleTurbulenceModels
wclean
# ----------------------------------------------------------------- end-of-file
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wmakeLnInclude interfacialModels
wmake libso twoPhaseSystem
wmake libso interfacialModels
wmake libso phaseIncompressibleTurbulenceModels
wmake
# ----------------------------------------------------------------- end-of-file
# include "CourantNo.H"
{
scalar UrCoNum = 0.5*gMax
(
fvc::surfaceSum(mag(phi1 - phi2))().internalField()/mesh.V().field()
)*runTime.deltaTValue();
Info<< "Max Ur Courant Number = " << UrCoNum << endl;
CoNum = max(CoNum, UrCoNum);
}
{
DDtU1 =
fvc::ddt(U1)
+ fvc::div(phi1, U1)
- fvc::div(phi1)*U1;
mrfZones.addCoriolis(U1, DDtU1);
DDtU2 =
fvc::ddt(U2)
+ fvc::div(phi2, U2)
- fvc::div(phi2)*U2;
mrfZones.addCoriolis(U2, DDtU2);
}
compressibleTwoPhaseEulerFoam.C
EXE = $(FOAM_APPBIN)/compressibleTwoPhaseEulerFoam
EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
-IphaseIncompressibleTurbulenceModels/lnInclude \
-IinterfacialModels/lnInclude \
-ItwoPhaseSystem/lnInclude \
-Iaveraging
EXE_LIBS = \
-lfluidThermophysicalModels \
-lspecie \
-lturbulenceModels \
-lincompressibleTurbulenceModels \
-lphaseIncompressibleTurbulenceModels \
-lincompressibleTransportModels \
-lcompressibleTwoPhaseSystem \
-lcompressibleEulerianInterfacialModels \
-lfiniteVolume \
-lmeshTools
mrfZones.correctBoundaryVelocity(U1);
mrfZones.correctBoundaryVelocity(U2);
mrfZones.correctBoundaryVelocity(U);
fvVectorMatrix U1Eqn(U1, U1.dimensions()*dimVol/dimTime);
fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
volScalarField dragCoeff(fluid.dragCoeff());
{
volVectorField liftForce(fluid.liftForce(U));
{
U1Eqn =
(
fvm::ddt(alpha1, U1)
+ fvm::div(alphaPhi1, U1)
// Compressibity correction
- fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1)
+ fluid.Cvm()*rho2*alpha1*alpha2/rho1*
(
fvm::ddt(U1)
+ fvm::div(phi1, U1)
- fvm::Sp(fvc::div(phi1), U1)
)
+ turbulence1->divDevReff(U1)
==
- fvm::Sp(dragCoeff/rho1, U1)
- alpha1*alpha2/rho1*(liftForce - fluid.Cvm()*rho2*DDtU2)
);
mrfZones.addCoriolis(alpha1*(1 + fluid.Cvm()*rho2*alpha2/rho1), U1Eqn);
U1Eqn.relax();
}
{
U2Eqn =
(
fvm::ddt(alpha2, U2)
+ fvm::div(alphaPhi2, U2)
// Compressibity correction
- fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2)
+ fluid.Cvm()*rho2*alpha1*alpha2/rho2*
(
fvm::ddt(U2)
+ fvm::div(phi2, U2)
- fvm::Sp(fvc::div(phi2), U2)
)
+ turbulence2->divDevReff(U2)
==
- fvm::Sp(dragCoeff/rho2, U2)
+ alpha1*alpha2/rho2*(liftForce + fluid.Cvm()*rho2*DDtU1)
);
mrfZones.addCoriolis(alpha2*(1 + fluid.Cvm()*rho2*alpha1/rho2), U2Eqn);
U2Eqn.relax();
}
}
{
word alphaScheme("div(phi," + alpha1.name() + ')');
word alpharScheme("div(phir," + alpha1.name() + ')');
surfaceScalarField phic("phic", phi);
surfaceScalarField phir("phir", phi1 - phi2);
surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
tmp<surfaceScalarField> pPrimeByA;
if (implicitPhasePressure)
{
pPrimeByA =
fvc::interpolate((1.0/rho1)*rAU1*turbulence1().pPrime())
+ fvc::interpolate((1.0/rho2)*rAU2*turbulence2().pPrime());
surfaceScalarField phiP
(
pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh.magSf()
);
phic += alpha1f*phiP;
phir += phiP;
}
for (int acorr=0; acorr<nAlphaCorr; acorr++)
{
volScalarField::DimensionedInternalField Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
);
volScalarField::DimensionedInternalField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
fvc::div(phi)*min(alpha1, scalar(1))
);
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]);
}
}
dimensionedScalar totalDeltaT = runTime.deltaT();
if (nAlphaSubCycles > 1)
{
alphaPhi1 = dimensionedScalar("0", alphaPhi1.dimensions(), 0);
}
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
surfaceScalarField alphaPhic1
(
fvc::flux
(
phic,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
alpha1,
alpharScheme
)
);
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phi,
alphaPhic1,
Sp,
Su,
1,
0
);
if (nAlphaSubCycles > 1)
{
alphaPhi1 += (runTime.deltaT()/totalDeltaT)*alphaPhic1;
}
else
{
alphaPhi1 = alphaPhic1;
}
}
if (implicitPhasePressure)
{
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1) - fvc::ddt(alpha1)
- fvm::laplacian(alpha1f*pPrimeByA, alpha1, "bounded")
);
alpha1Eqn.relax();
alpha1Eqn.solve();
alphaPhi1 += alpha1Eqn.flux();
}
alphaPhi2 = phi - alphaPhi1;
alpha2 = scalar(1) - alpha1;
Info<< "Dispersed phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}
}
rho = fluid.rho();
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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
compressibleTwoPhaseEulerFoam
Description
Solver for a system of 2 compressible fluid phases with one phase
dispersed, e.g. gas bubbles in a liquid including heat-transfer.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "rhoThermo.H"
#include "twoPhaseSystem.H"
#include "dragModel.H"
#include "heatTransferModel.H"
#include "PhaseIncompressibleTurbulenceModel.H"
#include "pimpleControl.H"
#include "IOMRFZoneList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
#include "createMRFZones.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
#include "CourantNos.H"
#include "setInitialDeltaT.H"
pimpleControl pimple(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTwoPhaseEulerFoamControls.H"
#include "CourantNos.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "alphaEqn.H"
#include "EEqns.H"
#include "UEqns.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
#include "DDtU.H"
if (pimple.turbCorr())
{
turbulence1->correct();
turbulence2->correct();
}
}
#include "write.H"
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
Info<< "Creating twoPhaseSystem\n" << endl;
twoPhaseSystem fluid(mesh);
phaseModel& phase1 = fluid.phase1();
phaseModel& phase2 = fluid.phase2();
volScalarField& alpha1 = phase1;
volScalarField& alpha2 = phase2;
volVectorField& U1 = phase1.U();
surfaceScalarField& phi1 = phase1.phi();
surfaceScalarField alphaPhi1
(
IOobject::groupName("alphaPhi", phase1.name()),
fvc::interpolate(alpha1)*phi1
);
volVectorField& U2 = phase2.U();
surfaceScalarField& phi2 = phase2.phi();
surfaceScalarField alphaPhi2
(
IOobject::groupName("alphaPhi", phase2.name()),
fvc::interpolate(alpha2)*phi2
);
dimensionedScalar pMin
(
"pMin",
dimPressure,
fluid.lookup("pMin")
);
rhoThermo& thermo1 = phase1.thermo();
rhoThermo& thermo2 = phase2.thermo();
volScalarField& p = thermo1.p();
volScalarField& rho1 = thermo1.rho();
const volScalarField& psi1 = thermo1.psi();
volScalarField& rho2 = thermo2.rho();
const volScalarField& psi2 = thermo2.psi();
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fluid.U()
);
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fluid.phi()
);
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fluid.rho()
);
Info<< "Calculating field DDtU1 and DDtU2\n" << endl;
volVectorField DDtU1
(
fvc::ddt(U1)
+ fvc::div(phi1, U1)
- fvc::div(phi1)*U1
);
volVectorField DDtU2
(
fvc::ddt(U2)
+ fvc::div(phi2, U2)
- fvc::div(phi2)*U2
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
volScalarField rAU1
(
IOobject
(
IOobject::groupName("rAU", phase1.name()),
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("zero", dimensionSet(0, 0, 1, 0, 0), 0.0)
);
volScalarField rAU2
(
IOobject
(
IOobject::groupName("rAU", phase2.name()),
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("zero", dimensionSet(0, 0, 1, 0, 0), 0.0)
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
volScalarField dgdt
(
pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
);
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K1(IOobject::groupName("K", phase1.name()), 0.5*magSqr(U1));
volScalarField K2(IOobject::groupName("K", phase2.name()), 0.5*magSqr(U2));
autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> >
turbulence1
(
PhaseIncompressibleTurbulenceModel<phaseModel>::New
(
alpha1,
U1,
alphaPhi1,
phi1,
phase1
)
);
autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> >
turbulence2
(
PhaseIncompressibleTurbulenceModel<phaseModel>::New
(
alpha2,
U2,
alphaPhi2,
phi2,
phase2
)
);