From db99252fc3b80b8ef518c4cb3a4e85167250562d Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Tue, 28 Feb 2012 15:51:02 +0000
Subject: [PATCH] MRFtwoPhaseEulerFoam: new solver; MRF version of
twoPhaseEulerFoam
---
.../MRFTwoPhaseEulerFoam.C | 119 +++++++++++++++++
.../MRFtwoPhaseEulerFoam/Make/files | 3 +
.../MRFtwoPhaseEulerFoam/Make/options | 17 +++
.../MRFtwoPhaseEulerFoam/UEqns.H | 99 ++++++++++++++
.../MRFtwoPhaseEulerFoam/createMRFZones.H | 3 +
.../MRFtwoPhaseEulerFoam/pEqn.H | 121 ++++++++++++++++++
6 files changed, 362 insertions(+)
create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C
create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files
create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options
create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H
create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H
create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C
new file mode 100644
index 00000000000..0620fecd67e
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C
@@ -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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files
new file mode 100644
index 00000000000..45960722aef
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files
@@ -0,0 +1,3 @@
+MRFTwoPhaseEulerFoam.C
+
+EXE = $(FOAM_APPBIN)/MRFTwoPhaseEulerFoam
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options
new file mode 100644
index 00000000000..b9b19059da6
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options
@@ -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
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H
new file mode 100644
index 00000000000..2aba5b9d3d0
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H
@@ -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();
+ }
+}
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H
new file mode 100644
index 00000000000..0d653536241
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H
@@ -0,0 +1,3 @@
+ MRFZones mrfZones(mesh);
+ mrfZones.correctBoundaryVelocity(Ua);
+ mrfZones.correctBoundaryVelocity(Ub);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H
new file mode 100644
index 00000000000..d89136b1721
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H
@@ -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"
--
GitLab