mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
MRFtwoPhaseEulerFoam: new solver; MRF version of twoPhaseEulerFoam
This commit is contained in:
Henry
@ -0,0 +1,119 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011 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
|
||||
twoPhaseEulerFoam
|
||||
|
||||
Description
|
||||
Solver for a system of 2 incompressible fluid phases with one phase
|
||||
dispersed, e.g. gas bubbles in a liquid or solid particles in a gas.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
#include "nearWallDist.H"
|
||||
#include "wallFvPatch.H"
|
||||
#include "Switch.H"
|
||||
|
||||
#include "IFstream.H"
|
||||
#include "OFstream.H"
|
||||
|
||||
#include "dragModel.H"
|
||||
#include "phaseModel.H"
|
||||
#include "kineticTheoryModel.H"
|
||||
|
||||
#include "pimpleControl.H"
|
||||
#include "MRFZones.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
#include "setRootCase.H"
|
||||
|
||||
#include "createTime.H"
|
||||
#include "createMesh.H"
|
||||
#include "readGravitationalAcceleration.H"
|
||||
#include "createFields.H"
|
||||
#include "readPPProperties.H"
|
||||
#include "initContinuityErrs.H"
|
||||
#include "createMRFZones.H"
|
||||
#include "readTimeControls.H"
|
||||
#include "CourantNo.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 "liftDragCoeffs.H"
|
||||
#include "UEqns.H"
|
||||
|
||||
// --- Pressure corrector loop
|
||||
while (pimple.correct())
|
||||
{
|
||||
#include "pEqn.H"
|
||||
|
||||
if (correctAlpha && !pimple.finalIter())
|
||||
{
|
||||
#include "alphaEqn.H"
|
||||
}
|
||||
}
|
||||
|
||||
#include "DDtU.H"
|
||||
|
||||
if (pimple.turbCorr())
|
||||
{
|
||||
#include "kEpsilon.H"
|
||||
}
|
||||
}
|
||||
|
||||
#include "write.H"
|
||||
|
||||
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
|
||||
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
|
||||
<< nl << endl;
|
||||
}
|
||||
|
||||
Info<< "End\n" << endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,3 @@
|
||||
MRFTwoPhaseEulerFoam.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/MRFTwoPhaseEulerFoam
|
||||
@ -0,0 +1,17 @@
|
||||
EXE_INC = \
|
||||
-I.. \
|
||||
-I../../bubbleFoam \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude \
|
||||
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
|
||||
-IturbulenceModel \
|
||||
-I../kineticTheoryModels/lnInclude \
|
||||
-I../interfacialModels/lnInclude \
|
||||
-I../phaseModel/lnInclude \
|
||||
|
||||
EXE_LIBS = \
|
||||
-lEulerianInterfacialModels \
|
||||
-lfiniteVolume \
|
||||
-lmeshTools \
|
||||
-lincompressibleTransportModels \
|
||||
-lphaseModel \
|
||||
-lkineticTheoryModel
|
||||
@ -0,0 +1,99 @@
|
||||
fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime);
|
||||
fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
|
||||
|
||||
{
|
||||
{
|
||||
volTensorField gradUaT(T(fvc::grad(Ua)));
|
||||
|
||||
if (kineticTheory.on())
|
||||
{
|
||||
kineticTheory.solve(gradUaT);
|
||||
nuEffa = kineticTheory.mua()/rhoa;
|
||||
}
|
||||
else // If not using kinetic theory is using Ct model
|
||||
{
|
||||
nuEffa = sqr(Ct)*nutb + nua;
|
||||
}
|
||||
|
||||
volTensorField Rca
|
||||
(
|
||||
"Rca",
|
||||
(((2.0/3.0)*I)*nuEffa)*tr(gradUaT) - nuEffa*gradUaT
|
||||
);
|
||||
|
||||
if (kineticTheory.on())
|
||||
{
|
||||
Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);
|
||||
}
|
||||
|
||||
surfaceScalarField phiRa
|
||||
(
|
||||
-fvc::interpolate(nuEffa)*mesh.magSf()*fvc::snGrad(alpha)
|
||||
/fvc::interpolate(alpha + scalar(0.001))
|
||||
);
|
||||
|
||||
UaEqn =
|
||||
(
|
||||
(scalar(1) + Cvm*rhob*beta/rhoa)*
|
||||
(
|
||||
fvm::ddt(Ua)
|
||||
+ fvm::div(phia, Ua, "div(phia,Ua)")
|
||||
- fvm::Sp(fvc::div(phia), Ua)
|
||||
)
|
||||
|
||||
- fvm::laplacian(nuEffa, Ua)
|
||||
+ fvc::div(Rca)
|
||||
|
||||
+ fvm::div(phiRa, Ua, "div(phia,Ua)")
|
||||
- fvm::Sp(fvc::div(phiRa), Ua)
|
||||
+ (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)
|
||||
==
|
||||
// g // Buoyancy term transfered to p-equation
|
||||
- fvm::Sp(beta/rhoa*K, Ua)
|
||||
//+ beta/rhoa*K*Ub // Explicit drag transfered to p-equation
|
||||
- beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)
|
||||
);
|
||||
mrfZones.addCoriolis(UaEqn);
|
||||
UaEqn.relax();
|
||||
}
|
||||
|
||||
{
|
||||
volTensorField gradUbT(T(fvc::grad(Ub)));
|
||||
volTensorField Rcb
|
||||
(
|
||||
"Rcb",
|
||||
(((2.0/3.0)*I)*nuEffb)*tr(gradUbT) - nuEffb*gradUbT
|
||||
);
|
||||
|
||||
surfaceScalarField phiRb
|
||||
(
|
||||
-fvc::interpolate(nuEffb)*mesh.magSf()*fvc::snGrad(beta)
|
||||
/fvc::interpolate(beta + scalar(0.001))
|
||||
);
|
||||
|
||||
UbEqn =
|
||||
(
|
||||
(scalar(1) + Cvm*rhob*alpha/rhob)*
|
||||
(
|
||||
fvm::ddt(Ub)
|
||||
+ fvm::div(phib, Ub, "div(phib,Ub)")
|
||||
- fvm::Sp(fvc::div(phib), Ub)
|
||||
)
|
||||
|
||||
- fvm::laplacian(nuEffb, Ub)
|
||||
+ fvc::div(Rcb)
|
||||
|
||||
+ fvm::div(phiRb, Ub, "div(phib,Ub)")
|
||||
- fvm::Sp(fvc::div(phiRb), Ub)
|
||||
|
||||
+ (fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb)
|
||||
==
|
||||
// g // Buoyancy term transfered to p-equation
|
||||
- fvm::Sp(alpha/rhob*K, Ub)
|
||||
//+ alpha/rhob*K*Ua // Explicit drag transfered to p-equation
|
||||
+ alpha/rhob*(liftCoeff + Cvm*rhob*DDtUa)
|
||||
);
|
||||
mrfZones.addCoriolis(UbEqn);
|
||||
UbEqn.relax();
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,3 @@
|
||||
MRFZones mrfZones(mesh);
|
||||
mrfZones.correctBoundaryVelocity(Ua);
|
||||
mrfZones.correctBoundaryVelocity(Ub);
|
||||
@ -0,0 +1,121 @@
|
||||
{
|
||||
surfaceScalarField alphaf(fvc::interpolate(alpha));
|
||||
surfaceScalarField betaf(scalar(1) - alphaf);
|
||||
|
||||
volScalarField rUaA(1.0/UaEqn.A());
|
||||
volScalarField rUbA(1.0/UbEqn.A());
|
||||
|
||||
rUaAf = fvc::interpolate(rUaA);
|
||||
surfaceScalarField rUbAf(fvc::interpolate(rUbA));
|
||||
|
||||
volVectorField HbyAa("HbyAa", Ua);
|
||||
HbyAa = rUaA*UaEqn.H();
|
||||
|
||||
volVectorField HbyAb("HbyAb", Ub);
|
||||
HbyAb = rUbA*UbEqn.H();
|
||||
|
||||
mrfZones.absoluteFlux(phia.oldTime());
|
||||
mrfZones.absoluteFlux(phia);
|
||||
|
||||
mrfZones.absoluteFlux(phib.oldTime());
|
||||
mrfZones.absoluteFlux(phib);
|
||||
|
||||
surfaceScalarField phiDraga
|
||||
(
|
||||
fvc::interpolate(beta/rhoa*K*rUaA)*phib + rUaAf*(g & mesh.Sf())
|
||||
);
|
||||
|
||||
if (g0.value() > 0.0)
|
||||
{
|
||||
phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
|
||||
}
|
||||
|
||||
if (kineticTheory.on())
|
||||
{
|
||||
phiDraga -= rUaAf*fvc::snGrad(kineticTheory.pa()/rhoa)*mesh.magSf();
|
||||
}
|
||||
|
||||
|
||||
surfaceScalarField phiDragb
|
||||
(
|
||||
fvc::interpolate(alpha/rhob*K*rUbA)*phia + rUbAf*(g & mesh.Sf())
|
||||
);
|
||||
|
||||
// Fix for gravity on outlet boundary.
|
||||
forAll(p.boundaryField(), patchi)
|
||||
{
|
||||
if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
|
||||
{
|
||||
phiDraga.boundaryField()[patchi] = 0.0;
|
||||
phiDragb.boundaryField()[patchi] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
surfaceScalarField phiHbyAa
|
||||
(
|
||||
"phiHbyAa",
|
||||
(fvc::interpolate(HbyAa) & mesh.Sf())
|
||||
+ fvc::ddtPhiCorr(rUaA, Ua, phia)
|
||||
+ phiDraga
|
||||
);
|
||||
mrfZones.relativeFlux(phiHbyAa);
|
||||
|
||||
surfaceScalarField phiHbyAb
|
||||
(
|
||||
"phiHbyAb",
|
||||
(fvc::interpolate(HbyAb) & mesh.Sf())
|
||||
+ fvc::ddtPhiCorr(rUbA, Ub, phib)
|
||||
+ phiDragb
|
||||
);
|
||||
mrfZones.relativeFlux(phiHbyAb);
|
||||
|
||||
surfaceScalarField phiHbyA("phiHbyA", alphaf*phiHbyAa + betaf*phiHbyAb);
|
||||
|
||||
surfaceScalarField Dp
|
||||
(
|
||||
"Dp",
|
||||
alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
|
||||
);
|
||||
|
||||
while (pimple.correctNonOrthogonal())
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
||||
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
|
||||
|
||||
if (pimple.finalNonOrthogonalIter())
|
||||
{
|
||||
surfaceScalarField SfGradp(pEqn.flux()/Dp);
|
||||
|
||||
phia.boundaryField() ==
|
||||
(fvc::interpolate(Ua) & mesh.Sf())().boundaryField();
|
||||
mrfZones.relativeFlux(phia);
|
||||
phia = phiHbyAa - rUaAf*SfGradp/rhoa;
|
||||
|
||||
phib.boundaryField() ==
|
||||
(fvc::interpolate(Ub) & mesh.Sf())().boundaryField();
|
||||
mrfZones.relativeFlux(phib);
|
||||
phib = phiHbyAb - rUbAf*SfGradp/rhob;
|
||||
|
||||
phi = alphaf*phia + betaf*phib;
|
||||
|
||||
p.relax();
|
||||
SfGradp = pEqn.flux()/Dp;
|
||||
|
||||
Ua = HbyAa + fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa);
|
||||
Ua.correctBoundaryConditions();
|
||||
|
||||
Ub = HbyAb + fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
|
||||
Ub.correctBoundaryConditions();
|
||||
|
||||
U = alpha*Ua + beta*Ub;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#include "continuityErrs.H"
|
||||
Reference in New Issue
Block a user