research/OpenFOAM-2.4.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'master' of github.com:OpenFOAM/OpenFOAM-2.3.x

Browse source
This commit is contained in:
sergio 2014-04-30 12:04:40 +01:00
commit cde7f593dc
620 changed files with 7488 additions and 10336 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
applications/solvers/multiphase/cavitatingFoam/createFields.H
View file

@ -57,12 +57,12 @@
Info<< "Reading transportProperties\n" << endl;
incompressibleTwoPhaseMixture twoPhaseProperties(U, phi);
incompressibleTwoPhaseMixture mixture(U, phi);
volScalarField& alphav(twoPhaseProperties.alpha1());
volScalarField& alphav(mixture.alpha1());
alphav.oldTime();
volScalarField& alphal(twoPhaseProperties.alpha2());
volScalarField& alphal(mixture.alpha2());
Info<< "Creating compressibilityModel\n" << endl;
autoPtr<barotropicCompressibilityModel> psiModel =
@ -85,5 +85,5 @@
// Create incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
incompressible::turbulenceModel::New(U, phi, mixture)
);

1
applications/solvers/multiphase/compressibleInterFoam/Make/options
View file

@ -1,5 +1,6 @@
EXE_INC = \
-ItwoPhaseMixtureThermo \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \

8
applications/solvers/multiphase/compressibleInterFoam/TEqn.H
View file

@ -3,19 +3,19 @@
(
fvm::ddt(rho, T)
+ fvm::div(rhoPhi, T)
- fvm::laplacian(twoPhaseProperties.alphaEff(turbulence->mut()), T)
- fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
+ (
fvc::div(fvc::absolute(phi, U), p)
+ fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
)
*(
alpha1/twoPhaseProperties.thermo1().Cv()
+ alpha2/twoPhaseProperties.thermo2().Cv()
alpha1/mixture.thermo1().Cv()
+ alpha2/mixture.thermo2().Cv()
)
);
TEqn.relax();
TEqn.solve();
twoPhaseProperties.correct();
mixture.correct();
}

2
applications/solvers/multiphase/compressibleInterFoam/UEqn.H
View file

@ -16,7 +16,7 @@
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
interface.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()

1
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options
View file

@ -1,6 +1,7 @@
EXE_INC = \
-I.. \
-I../twoPhaseMixtureThermo \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \

2
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
View file

@ -15,7 +15,7 @@
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
interface.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);

24
applications/solvers/multiphase/compressibleInterFoam/createFields.H
View file

@ -29,19 +29,19 @@
#include "createPhi.H"
Info<< "Constructing twoPhaseMixtureThermo\n" << endl;
twoPhaseMixtureThermo twoPhaseProperties(mesh);
twoPhaseMixtureThermo mixture(mesh);
volScalarField& alpha1(twoPhaseProperties.alpha1());
volScalarField& alpha2(twoPhaseProperties.alpha2());
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
Info<< "Reading thermophysical properties\n" << endl;
volScalarField& p = twoPhaseProperties.p();
volScalarField& T = twoPhaseProperties.T();
volScalarField& rho1 = twoPhaseProperties.thermo1().rho();
const volScalarField& psi1 = twoPhaseProperties.thermo1().psi();
volScalarField& rho2 = twoPhaseProperties.thermo2().rho();
const volScalarField& psi2 = twoPhaseProperties.thermo2().psi();
volScalarField& p = mixture.p();
volScalarField& T = mixture.T();
volScalarField& rho1 = mixture.thermo1().rho();
const volScalarField& psi1 = mixture.thermo1().psi();
volScalarField& rho2 = mixture.thermo2().rho();
const volScalarField& psi2 = mixture.thermo2().psi();
volScalarField rho
(
@ -57,7 +57,7 @@
);
dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
dimensionedScalar pMin(mixture.lookup("pMin"));
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
@ -85,12 +85,12 @@
);
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);
interfaceProperties interface(alpha1, U, mixture);
// Construct compressible turbulence model
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New(rho, U, rhoPhi, twoPhaseProperties)
compressible::turbulenceModel::New(rho, U, rhoPhi, mixture)
);
Info<< "Creating field kinetic energy K\n" << endl;

2
applications/solvers/multiphase/compressibleInterFoam/pEqn.H
View file

@ -15,7 +15,7 @@
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
interface.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);

1
applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/Make/options
View file

@ -1,4 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude

2
applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
View file

@ -1,5 +1,7 @@
EXE_INC = \
-I../interFoam \
-ImultiphaseMixtureThermo/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \

10
applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H
View file

@ -2,16 +2,16 @@
fvScalarMatrix TEqn
(
fvm::ddt(rho, T)
+ fvm::div(multiphaseProperties.rhoPhi(), T)
- fvm::laplacian(multiphaseProperties.alphaEff(turbulence->mut()), T)
+ fvm::div(mixture.rhoPhi(), T)
- fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
+ (
fvc::div(fvc::absolute(phi, U), p)
+ fvc::ddt(rho, K) + fvc::div(multiphaseProperties.rhoPhi(), K)
)*multiphaseProperties.rCv()
+ fvc::ddt(rho, K) + fvc::div(mixture.rhoPhi(), K)
)*mixture.rCv()
);
TEqn.relax();
TEqn.solve();
multiphaseProperties.correct();
mixture.correct();
}

4
applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H
View file

@ -1,7 +1,7 @@
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(multiphaseProperties.rhoPhi(), U)
+ fvm::div(mixture.rhoPhi(), U)
+ turbulence->divDevRhoReff(U)
);
@ -16,7 +16,7 @@
fvc::reconstruct
(
(
multiphaseProperties.surfaceTensionForce()
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()

6
applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -76,9 +76,9 @@ int main(int argc, char *argv[])
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
multiphaseProperties.solve();
mixture.solve();
solve(fvm::ddt(rho) + fvc::div(multiphaseProperties.rhoPhi()));
solve(fvm::ddt(rho) + fvc::div(mixture.rhoPhi()));
#include "UEqn.H"
#include "TEqn.H"

14
applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
View file

@ -29,10 +29,10 @@
#include "createPhi.H"
Info<< "Constructing multiphaseMixtureThermo\n" << endl;
multiphaseMixtureThermo multiphaseProperties(U, phi);
multiphaseMixtureThermo mixture(U, phi);
volScalarField& p = multiphaseProperties.p();
volScalarField& T = multiphaseProperties.T();
volScalarField& p = mixture.p();
volScalarField& T = mixture.T();
volScalarField rho
(
@ -43,12 +43,12 @@
mesh,
IOobject::READ_IF_PRESENT
),
multiphaseProperties.rho()
mixture.rho()
);
Info<< max(rho) << min(rho);
dimensionedScalar pMin(multiphaseProperties.lookup("pMin"));
dimensionedScalar pMin(mixture.lookup("pMin"));
Info<< "Calculating field g.h\n" << endl;
@ -62,8 +62,8 @@
(
rho,
U,
multiphaseProperties.rhoPhi(),
multiphaseProperties
mixture.rhoPhi(),
mixture
)
);

1
applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options
View file

@ -1,4 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude

20
applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
View file

@ -15,7 +15,7 @@
surfaceScalarField phig
(
(
multiphaseProperties.surfaceTensionForce()
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);
@ -32,13 +32,13 @@
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);
PtrList<fvScalarMatrix> p_rghEqnComps(multiphaseProperties.phases().size());
PtrList<fvScalarMatrix> p_rghEqnComps(mixture.phases().size());
label phasei = 0;
forAllConstIter
(
PtrDictionary<phaseModel>,
multiphaseProperties.phases(),
mixture.phases(),
phase
)
{
@ -74,7 +74,7 @@
forAllConstIter
(
PtrDictionary<phaseModel>,
multiphaseProperties.phases(),
mixture.phases(),
phase
)
{
@ -105,14 +105,14 @@
if (pimple.finalNonOrthogonalIter())
{
// p = max(p_rgh + multiphaseProperties.rho()*gh, pMin);
// p_rgh = p - multiphaseProperties.rho()*gh;
// p = max(p_rgh + mixture.rho()*gh, pMin);
// p_rgh = p - mixture.rho()*gh;
phasei = 0;
forAllIter
(
PtrDictionary<phaseModel>,
multiphaseProperties.phases(),
mixture.phases(),
phase
)
{
@ -129,11 +129,11 @@
}
}
p = max(p_rgh + multiphaseProperties.rho()*gh, pMin);
p = max(p_rgh + mixture.rho()*gh, pMin);
// Update densities from change in p_rgh
multiphaseProperties.correctRho(p_rgh - p_rgh_0);
rho = multiphaseProperties.rho();
mixture.correctRho(p_rgh - p_rgh_0);
rho = mixture.rho();
K = 0.5*magSqr(U);

1
applications/solvers/multiphase/driftFluxFoam/Make/files
View file

@ -1,4 +1,5 @@
incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.C
compressibleTurbulenceModels.C
driftFluxFoam.C
EXE = $(FOAM_APPBIN)/driftFluxFoam

8
applications/solvers/multiphase/driftFluxFoam/Make/options
View file

@ -9,13 +9,17 @@ EXE_INC = \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I./relativeVelocityModels/lnInclude
EXE_LIBS = \
-ldriftFluxTransportModels \
-ldriftFluxRelativeVelocityModels \
-lfiniteVolume \
-lmeshTools \
-lsampling \
-lfvOptions \
-lincompressibleTransportModels \
-ldriftFluxTransportModels \
-ldriftFluxRelativeVelocityModels
-lturbulenceModels \
-lcompressibleTurbulenceModels

3
applications/solvers/multiphase/driftFluxFoam/UEqn.H
View file

@ -5,8 +5,7 @@
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
+ fvc::div(UdmModel.tauDm())
- fvm::laplacian(muEff, U)
- fvc::div(muEff*dev2(T(fvc::grad(U))))
+ turbulence->divDevRhoReff(U)
==
fvOptions(rho, U)
);

56
applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
View file

@ -27,7 +27,7 @@
mesh,
linear<scalar>(mesh),
fv::uncorrectedSnGrad<scalar>(mesh)
).fvmLaplacian(fvc::interpolate(mut/rho), alpha1)
).fvmLaplacian(fvc::interpolate(turbulence->nut()), alpha1)
);
Info<< "Phase-1 volume fraction = "
@ -43,9 +43,23 @@
{
Info<< "Applying the previous iteration correction flux" << endl;
#ifdef LTSSOLVE
MULES::LTScorrect(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
MULES::LTScorrect
(
alpha1,
phiAlpha,
tphiAlphaCorr0(),
mixture.alphaMax(),
0
);
#else
MULES::correct(alpha1, phiAlpha, tphiAlphaCorr0(), 1, 0);
MULES::correct
(
alpha1,
phiAlpha,
tphiAlphaCorr0(),
mixture.alphaMax(),
0
);
#endif
phiAlpha += tphiAlphaCorr0();
@ -79,9 +93,23 @@
volScalarField alpha10(alpha1);
#ifdef LTSSOLVE
MULES::LTScorrect(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
MULES::LTScorrect
(
alpha1,
tphiAlphaUn(),
tphiAlphaCorr(),
mixture.alphaMax(),
0
);
#else
MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
MULES::correct
(
alpha1,
tphiAlphaUn(),
tphiAlphaCorr(),
mixture.alphaMax(),
0
);
#endif
// Under-relax the correction for all but the 1st corrector
@ -100,9 +128,23 @@
phiAlpha = tphiAlphaUn;
#ifdef LTSSOLVE
MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
MULES::explicitLTSSolve
(
alpha1,
phi,
phiAlpha,
mixture.alphaMax(),
0
);
#else
MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
MULES::explicitSolve
(
alpha1,
phi,
phiAlpha,
mixture.alphaMax(),
0
);
#endif
}
}

4
applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
View file

@ -55,7 +55,7 @@
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1) - fvc::ddt(alpha1)
- fvm::laplacian(mut/rho, alpha1)
- fvm::laplacian(turbulence->nut(), alpha1)
);
alpha1Eqn.solve(mesh.solver("alpha1Diffusion"));
@ -71,5 +71,5 @@
}
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
rho == alpha1*rho1 + alpha2*rho2;
rho = mixture.rho();
}

73
applications/solvers/multiphase/driftFluxFoam/compressibleTurbulenceModels.C Normal file
View file

@ -0,0 +1,73 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 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/>.
\*---------------------------------------------------------------------------*/
#include "CompressibleTurbulenceModel.H"
#include "incompressibleTwoPhaseInteractingMixture.H"
#include "addToRunTimeSelectionTable.H"
#include "makeTurbulenceModel.H"
#include "laminar.H"
#include "RASModel.H"
#include "LESModel.H"
makeBaseTurbulenceModel
(
geometricOneField,
volScalarField,
compressibleTurbulenceModel,
CompressibleTurbulenceModel,
incompressibleTwoPhaseInteractingMixture
);
#define makeRASModel(Type) \
makeTemplatedTurbulenceModel \
( \
incompressibleTwoPhaseInteractingMixtureCompressibleTurbulenceModel, \
RAS, \
Type \
)
#define makeLESModel(Type) \
makeTemplatedTurbulenceModel \
( \
incompressibleTwoPhaseInteractingMixtureCompressibleTurbulenceModel, \
LES, \
Type \
)
#include "kEpsilon.H"
makeRASModel(kEpsilon);
#include "buoyantKEpsilon.H"
makeRASModel(buoyantKEpsilon);
#include "Smagorinsky.H"
makeLESModel(Smagorinsky);
#include "kEqn.H"
makeLESModel(kEqn);
// ************************************************************************* //

196
applications/solvers/multiphase/driftFluxFoam/createFields.H
View file

@ -29,17 +29,14 @@
#include "createPhi.H"
// Transport
// ~~~~~~~~~
Info<< "Reading transportProperties\n" << endl;
incompressibleTwoPhaseInteractingMixture twoPhaseProperties(U, phi);
incompressibleTwoPhaseInteractingMixture mixture(U, phi);
volScalarField& alpha1(twoPhaseProperties.alpha1());
volScalarField& alpha2(twoPhaseProperties.alpha2());
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = twoPhaseProperties.rhod();
const dimensionedScalar& rho2 = twoPhaseProperties.rhoc();
const dimensionedScalar& rho1 = mixture.rhod();
const dimensionedScalar& rho2 = mixture.rhoc();
IOdictionary transportProperties
(
@ -64,9 +61,8 @@
IOobject::NO_READ,
IOobject::NO_WRITE
),
alpha1*rho1 + alpha2*rho2
mixture.rho()
);
rho.oldTime();
// Mass flux
surfaceScalarField rhoPhi
@ -84,190 +80,28 @@
// Relative Velocity
// ~~~~~~~~~~~~~~~~~
autoPtr<relativeVelocityModel> UdmModelPtr
(
relativeVelocityModel::New
(
transportProperties,
twoPhaseProperties
mixture
)
);
relativeVelocityModel& UdmModel(UdmModelPtr());
// Turbulence
// ~~~~~~~~~~
IOdictionary RASProperties
// Construct compressible turbulence model
autoPtr
<
CompressibleTurbulenceModel<incompressibleTwoPhaseInteractingMixture>
> turbulence
(
IOobject
(
"RASProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
CompressibleTurbulenceModel<incompressibleTwoPhaseInteractingMixture>
::New(rho, U, rhoPhi, mixture)
);
Switch turbulence(RASProperties.lookup("turbulence"));
dictionary kEpsilonDict(RASProperties.subDictPtr("kEpsilonCoeffs"));
dimensionedScalar Cmu
(
dimensionedScalar::lookupOrAddToDict
(
"Cmu",
kEpsilonDict,
0.09
)
);
dimensionedScalar C1
(
dimensionedScalar::lookupOrAddToDict
(
"C1",
kEpsilonDict,
1.44
)
);
dimensionedScalar C2
(
dimensionedScalar::lookupOrAddToDict
(
"C2",
kEpsilonDict,
1.92
)
);
dimensionedScalar C3
(
dimensionedScalar::lookupOrAddToDict
(
"C3",
kEpsilonDict,
0.85
)
);
dimensionedScalar sigmak
(
dimensionedScalar::lookupOrAddToDict
(
"sigmak",
kEpsilonDict,
1.0
)
);
dimensionedScalar sigmaEps
(
dimensionedScalar::lookupOrAddToDict
(
"sigmaEps",
kEpsilonDict,
1.3
)
);
dictionary wallFunctionDict(RASProperties.subDictPtr("wallFunctionCoeffs"));
dimensionedScalar kappa
(
dimensionedScalar::lookupOrAddToDict
(
"kappa",
wallFunctionDict,
0.41
)
);
dimensionedScalar E
(
dimensionedScalar::lookupOrAddToDict
(
"E",
wallFunctionDict,
9.8
)
);
if (RASProperties.lookupOrDefault("printCoeffs", false))
{
Info<< "kEpsilonCoeffs" << kEpsilonDict << nl
<< "wallFunctionCoeffs" << wallFunctionDict << endl;
}
nearWallDist y(mesh);
Info<< "Reading field k\n" << endl;
volScalarField k
(
IOobject
(
"k",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field epsilon\n" << endl;
volScalarField epsilon
(
IOobject
(
"epsilon",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Calculating field mut\n" << endl;
volScalarField mut
(
IOobject
(
"mut",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Cmu*rho*sqr(k)/epsilon
);
Info<< "Calculating field muEff\n" << endl;
volScalarField muEff
(
IOobject
(
"muEff",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mut + twoPhaseProperties.mu()
);
// Pressure
// ~~~~~~~~
Info<< "Calculating field (g.h)f\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
@ -309,6 +143,4 @@
// MULES Correction
// ~~~~~~~~~~~~~~~~
tmp<surfaceScalarField> tphiAlphaCorr0;

10
applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
View file

@ -38,10 +38,8 @@ Description
#include "subCycle.H"
#include "incompressibleTwoPhaseInteractingMixture.H"
#include "relativeVelocityModel.H"
#include "nearWallDist.H"
#include "wallFvPatch.H"
#include "bound.H"
#include "Switch.H"
#include "turbulenceModel.H"
#include "CompressibleTurbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "fixedFluxPressureFvPatchScalarField.H"
@ -86,7 +84,7 @@ int main(int argc, char *argv[])
#include "alphaEqnSubCycle.H"
twoPhaseProperties.correct();
mixture.correct();
#include "UEqn.H"
@ -98,7 +96,7 @@ int main(int argc, char *argv[])
if (pimple.turbCorr())
{
#include "kEpsilon.H"
turbulence->correct();
}
}

21
applications/solvers/multiphase/driftFluxFoam/incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.C
View file

@ -82,6 +82,13 @@ incompressibleTwoPhaseInteractingMixture
rhod_("rho", dimDensity, muModel_->viscosityProperties().lookup("rho")),
rhoc_("rho", dimDensity, nucModel_->viscosityProperties().lookup("rho")),
dd_
(
"d",
dimLength,
muModel_->viscosityProperties().lookupOrDefault("d", 0.0)
),
alphaMax_(muModel_->viscosityProperties().lookupOrDefault("alphaMax", 1.0)),
U_(U),
phi_(phi),
@ -118,6 +125,20 @@ bool Foam::incompressibleTwoPhaseInteractingMixture::read()
muModel_->viscosityProperties().lookup("rho") >> rhod_;
nucModel_->viscosityProperties().lookup("rho") >> rhoc_;
dd_ = dimensionedScalar
(
"d",
dimLength,
muModel_->viscosityProperties().lookupOrDefault("d", 0)
);
alphaMax_ =
muModel_->viscosityProperties().lookupOrDefault
(
"alphaMax",
1.0
);
return true;
}
else

48
applications/solvers/multiphase/driftFluxFoam/incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.H
View file

@ -69,6 +69,12 @@ protected:
dimensionedScalar rhod_;
dimensionedScalar rhoc_;
//- Optional diameter of the dispersed phase particles
dimensionedScalar dd_;
//- Optional maximum dispersed phase-fraction (e.g. packing limit)
scalar alphaMax_;
const volVectorField& U_;
const surfaceScalarField& phi_;
@ -121,6 +127,19 @@ public:
return rhoc_;
};
//- Return the diameter of the dispersed-phase particles
const dimensionedScalar& dd() const
{
return dd_;
}
//- Optional maximum phase-fraction (e.g. packing limit)
// Defaults to 1
scalar alphaMax() const
{
return alphaMax_;
}
//- Return const-access to the mixture velocity
const volVectorField& U() const
{
@ -133,21 +152,36 @@ public:
return mu_;
}
//- Return the dynamic mixture viscosity for patch
virtual tmp<scalarField> mu(const label patchi) const
{
return mu_.boundaryField()[patchi];
}
//- Return the mixture density
virtual tmp<volScalarField> rho() const
{
return alpha1_*rhod_ + alpha2_*rhoc_;
}
//- Return the mixture density for patch
virtual tmp<scalarField> rho(const label patchi) const
{
return
alpha1_.boundaryField()[patchi]*rhod_.value()
+ alpha2_.boundaryField()[patchi]*rhoc_.value();
}
//- Return the mixture viscosity
virtual tmp<volScalarField> nu() const
{
notImplemented("incompressibleTwoPhaseInteractingMixture::nu()");
return volScalarField::null();
return mu_/rho();
}
//- Return the mixture viscosity for patch
virtual tmp<scalarField> nu(const label patchi) const
{
notImplemented
(
"incompressibleTwoPhaseInteractingMixture::nu(const label)"
);
return scalarField::null();
return mu_.boundaryField()[patchi]/rho(patchi);
}
//- Correct the laminar viscosity

82
applications/solvers/multiphase/driftFluxFoam/kEpsilon.H
View file

@ -1,82 +0,0 @@
if (turbulence)
{
if (mesh.changing())
{
y.correct();
}
dimensionedScalar k0("k0", k.dimensions(), 0);
dimensionedScalar kMin("kMin", k.dimensions(), SMALL);
dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0);
dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField G(mut*(tgradU() && dev(twoSymm(tgradU()))));
tgradU.clear();
volScalarField Gcoef
(
Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin)
);
volScalarField mul(twoPhaseProperties.mu());
#include "wallFunctions.H"
// Dissipation equation
fvScalarMatrix epsEqn
(
fvm::ddt(rho, epsilon)
+ fvm::div(rhoPhi, epsilon)
- fvm::laplacian
(
mut/sigmaEps + mul, epsilon,
"laplacian(DepsilonEff,epsilon)"
)
==
C1*G*epsilon/(k + kMin)
- fvm::SuSp(C1*(1.0 - C3)*Gcoef, epsilon)
- fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon)
);
#include "wallDissipation.H"
epsEqn.relax();
epsEqn.solve();
bound(epsilon, epsilon0);
// Turbulent kinetic energy equation
fvScalarMatrix kEqn
(
fvm::ddt(rho, k)
+ fvm::div(rhoPhi, k)
- fvm::laplacian
(
mut/sigmak + mul, k,
"laplacian(DkEff,k)"
)
==
G
- fvm::SuSp(Gcoef, k)
- fvm::Sp(rho*epsilon/(k + kMin), k)
);
kEqn.relax();
kEqn.solve();
bound(k, k0);
//- Re-calculate viscosity
mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin);
#include "wallViscosity.H"
muEff = mut + mul;
}
else
{
muEff = mut + twoPhaseProperties.mu();
}

4
applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.C
View file

@ -54,7 +54,9 @@ Foam::relativeVelocityModel::relativeVelocityModel
(
"Udm",
alphac_.time().timeName(),
alphac_.mesh()
alphac_.mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
alphac_.mesh(),
dimensionedVector("Udm", dimVelocity, vector::zero),

85
applications/solvers/multiphase/driftFluxFoam/wallFunctions.H
View file

@ -1,85 +0,0 @@
{
labelList cellBoundaryFaceCount(epsilon.size(), 0);
const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
const scalar Cmu75 = ::pow(Cmu.value(), 0.75);
const scalar kappa_ = kappa.value();
const fvPatchList& patches = mesh.boundary();
//- Initialise the near-wall P field to zero
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
epsilon[faceCelli] = 0.0;
G[faceCelli] = 0.0;
}
}
}
//- Accumulate the wall face contributions to epsilon and G
// Increment cellBoundaryFaceCount for each face for averaging
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
const scalarField& mutw = mut.boundaryField()[patchi];
const scalarField& mulw = mul.boundaryField()[patchi];
scalarField magFaceGradU
(
mag(U.boundaryField()[patchi].snGrad())
);
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
// For corner cells (with two boundary or more faces),
// epsilon and G in the near-wall cell are calculated
// as an average
cellBoundaryFaceCount[faceCelli]++;
epsilon[faceCelli] +=
Cmu75*::pow(k[faceCelli], 1.5)
/(kappa_*y[patchi][facei]);
G[faceCelli] +=
(mutw[facei] + mulw[facei])
*magFaceGradU[facei]
*Cmu25*::sqrt(k[faceCelli])
/(kappa_*y[patchi][facei]);
}
}
}
// perform the averaging
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
cellBoundaryFaceCount[faceCelli] = 1;
}
}
}
}

38
applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
View file

@ -1,38 +0,0 @@
{
const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
const scalar kappa_ = kappa.value();
const scalar E_ = E.value();
const fvPatchList& patches = mesh.boundary();
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
scalarField& mutw = mut.boundaryField()[patchi];
const scalarField& mulw = mul.boundaryField()[patchi];
const scalarField& rhow = rho.boundaryField()[patchi];
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
scalar yPlus =
Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
/(mulw[facei]/rhow[facei]);
if (yPlus > 11.6)
{
mutw[facei] =
mulw[facei]*(yPlus*kappa_/::log(E_*yPlus) - 1);
}
else
{
mutw[facei] = 0.0;
}
}
}
}
}

1
applications/solvers/multiphase/interFoam/Allwclean
View file

@ -4,7 +4,6 @@ set -x
wclean
wclean interDyMFoam
wclean MRFInterFoam
wclean porousInterFoam
wclean LTSInterFoam
wclean interMixingFoam

1
applications/solvers/multiphase/interFoam/Allwmake
View file

@ -4,7 +4,6 @@ set -x
wmake
wmake interDyMFoam
wmake MRFInterFoam
wmake porousInterFoam
wmake LTSInterFoam
wmake interMixingFoam

7
applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
View file

@ -40,8 +40,7 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "incompressibleTwoPhaseMixture.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"
#include "fvcSmooth.H"
#include "pimpleControl.H"
@ -82,12 +81,12 @@ int main(int argc, char *argv[])
{
#include "alphaControls.H"
twoPhaseProperties.correct();
#define LTSSOLVE
#include "alphaEqnSubCycle.H"
#undef LTSSOLVE
mixture.correct();
turbulence->correct();
#include "UEqn.H"

6
applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options
View file

@ -5,16 +5,14 @@ EXE_INC = \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-ltwoPhaseMixture \
-ltwoPhaseProperties \
-lincompressibleTransportModels \
-limmiscibleIncompressibleTwoPhaseMixture \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \

122
applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
View file

@ -1,122 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 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
MRFInterFoam
Description
Solver for 2 incompressible, isothermal immiscible fluids using 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.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
For a two-fluid approach see twoPhaseEulerFoam.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "incompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"
#include "IOMRFZoneList.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createMRFZones.H"
#include "readTimeControls.H"
pimpleControl pimple(mesh);
#include "createPrghCorrTypes.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "alphaControls.H"
twoPhaseProperties.correct();
#include "alphaEqnSubCycle.H"
interface.correct();
#include "zonePhaseVolumes.H"
#include "UEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.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;
}
// ************************************************************************* //

3
applications/solvers/multiphase/interFoam/MRFInterFoam/Make/files
View file

@ -1,3 +0,0 @@
MRFInterFoam.C
EXE = $(FOAM_APPBIN)/MRFInterFoam

38
applications/solvers/multiphase/interFoam/MRFInterFoam/UEqn.H
View file

@ -1,38 +0,0 @@
surfaceScalarField muEff
(
"muEff",
twoPhaseProperties.muf()
+ fvc::interpolate(rho*turbulence->nut())
);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
==
fvOptions(rho, U)
);
mrfZones.addCoriolis(rho, UEqn);
UEqn.relax();
fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
}

2
applications/solvers/multiphase/interFoam/MRFInterFoam/createMRFZones.H
View file

@ -1,2 +0,0 @@
IOMRFZoneList mrfZones(mesh);
mrfZones.correctBoundaryVelocity(U);

72
applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
View file

@ -1,72 +0,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);
mrfZones.makeRelative(phiHbyA);
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);
phiHbyA += phig;
// Update the fixedFluxPressure BCs to ensure flux consistency
setSnGrad<fixedFluxPressureFvPatchScalarField>
(
p_rgh.boundaryField(),
(
phiHbyA.boundaryField()
- mrfZones.relative(mesh.Sf().boundaryField() & U.boundaryField())
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
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);
}
}
#include "continuityErrs.H"
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

21
applications/solvers/multiphase/interFoam/MRFInterFoam/zonePhaseVolumes.H
View file

@ -1,21 +0,0 @@
{
const scalarField& V = mesh.V();
forAll(mesh.cellZones(), czi)
{
const labelList& cellLabels = mesh.cellZones()[czi];
scalar phaseVolume = 0;
forAll(cellLabels, cli)
{
label celli = cellLabels[cli];
phaseVolume += alpha1[celli]*V[celli];
}
reduce(phaseVolume, sumOp<scalar>());
Info<< "Phase volume in zone " << mesh.cellZones()[czi].name()
<< " = " << phaseVolume*1e6 << " ml " << endl;
}
}

12
applications/solvers/multiphase/interFoam/Make/options
View file

@ -1,23 +1,21 @@
EXE_INC = -ggdb3 \
EXE_INC = \
-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/incompressible/turbulenceModel \
-I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-ltwoPhaseMixture \
-linterfaceProperties \
-ltwoPhaseProperties \
-lincompressibleTransportModels \
-limmiscibleIncompressibleTwoPhaseMixture \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \
-lmeshTools \
-lfvOptions \
-lmeshTools \
-lsampling

8
applications/solvers/multiphase/interFoam/UEqn.H
View file

@ -3,10 +3,14 @@
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
+ turbulence->divDevRhoReff(rho, U)
==
fvOptions(rho, U)
);
UEqn.relax();
fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
@ -16,10 +20,12 @@
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
fvOptions.correct(U);
}

4
applications/solvers/multiphase/interFoam/alphaCourantNo.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -41,7 +41,7 @@ if (mesh.nInternalFaces())
{
scalarField sumPhi
(
pos(alpha1 - 0.01)*pos(0.99 - alpha1)
mixture.nearInterface()().internalField()
*fvc::surfaceSum(mag(phi))().internalField()
);

10
applications/solvers/multiphase/interFoam/alphaEqn.H
View file

@ -3,13 +3,13 @@
word alpharScheme("div(phirb,alpha)");
// Standard face-flux compression coefficient
surfaceScalarField phic(interface.cAlpha()*mag(phi/mesh.magSf()));
surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
// Add the optional isotropic compression contribution
if (icAlpha > 0)
{
phic *= (1.0 - icAlpha);
phic += (interface.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
}
// Do not compress interface at non-coupled boundary faces
@ -74,12 +74,12 @@
alpha2 = 1.0 - alpha1;
interface.correct();
mixture.correct();
}
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
surfaceScalarField phir(phic*interface.nHatf());
surfaceScalarField phir(phic*mixture.nHatf());
tmp<surfaceScalarField> tphiAlphaUn
(
@ -132,7 +132,7 @@
alpha2 = 1.0 - alpha1;
interface.correct();
mixture.correct();
}
rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;

18
applications/solvers/multiphase/interFoam/createFields.H
View file

@ -30,13 +30,13 @@
Info<< "Reading transportProperties\n" << endl;
incompressibleTwoPhaseMixture twoPhaseProperties(U, phi);
immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
volScalarField& alpha1(twoPhaseProperties.alpha1());
volScalarField& alpha2(twoPhaseProperties.alpha2());
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
// Need to store rho for ddt(rho, U)
@ -70,14 +70,10 @@
);
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
incompressible::turbulenceModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
@ -131,7 +127,9 @@
p_rgh = p - rho*gh;
}
fv::IOoptionList fvOptions(mesh);
// MULES Correction
tmp<surfaceScalarField> tphiAlphaCorr0;

6
applications/solvers/multiphase/interFoam/interDyMFoam/Make/options
View file

@ -5,6 +5,7 @@ EXE_INC = \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
@ -13,10 +14,7 @@ EXE_INC = \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-ltwoPhaseMixture \
-linterfaceProperties \
-ltwoPhaseProperties \
-lincompressibleTransportModels \
-limmiscibleIncompressibleTwoPhaseMixture \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \

11
applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
View file

@ -36,8 +36,7 @@ Description
#include "dynamicFvMesh.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "incompressibleTwoPhaseMixture.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
@ -121,7 +120,7 @@ int main(int argc, char *argv[])
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
interface.correct();
mixture.correct();
}
if (mesh.changing() && checkMeshCourantNo)
@ -131,11 +130,9 @@ int main(int argc, char *argv[])
}
#include "alphaControls.H"
twoPhaseProperties.correct();
#include "alphaEqnSubCycle.H"
interface.correct();
mixture.correct();
#include "UEqn.H"

2
applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
View file

@ -22,7 +22,7 @@
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);

11
applications/solvers/multiphase/interFoam/interFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -40,8 +40,7 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "incompressibleTwoPhaseMixture.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
@ -84,11 +83,9 @@ int main(int argc, char *argv[])
while (pimple.loop())
{
#include "alphaControls.H"
twoPhaseProperties.correct();
#include "alphaEqnSubCycle.H"
interface.correct();
mixture.correct();
#include "UEqn.H"

3
applications/solvers/multiphase/interFoam/interMixingFoam/Make/files
View file

@ -1,5 +1,6 @@
incompressibleThreePhaseMixture/threePhaseMixture.C
incompressibleThreePhaseMixture/incompressibleThreePhaseMixture.C
threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
immiscibleIncompressibleThreePhaseMixture/immiscibleIncompressibleThreePhaseMixture.C
interMixingFoam.C
EXE = $(FOAM_APPBIN)/interMixingFoam

1
applications/solvers/multiphase/interFoam/interMixingFoam/Make/options
View file

@ -1,6 +1,7 @@
EXE_INC = \
-I.. \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-IimmiscibleIncompressibleThreePhaseMixture \
-IincompressibleThreePhaseMixture \
-IthreePhaseInterfaceProperties \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \

60
applications/solvers/multiphase/interFoam/interMixingFoam/alphaCourantNo.H
View file

@ -1,60 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 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/>.
Global
CourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar maxAlphaCo
(
readScalar(runTime.controlDict().lookup("maxAlphaCo"))
);
scalar alphaCoNum = 0.0;
scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
scalarField sumPhi
(
max
(
pos(alpha1 - 0.01)*pos(0.99 - alpha1),
pos(alpha2 - 0.01)*pos(0.99 - alpha2)
)*fvc::surfaceSum(mag(phi))().internalField()
);
alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
meanAlphaCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
}
Info<< "Interface Courant Number mean: " << meanAlphaCoNum
<< " max: " << alphaCoNum << endl;
// ************************************************************************* //

2
applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H
View file

@ -12,7 +12,7 @@
IOobject::NO_READ,
IOobject::NO_WRITE
),
interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf()
mixture.cAlpha()*mag(phi/mesh.magSf())*mixture.nHatf()
);
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)

22
applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
View file

@ -28,17 +28,17 @@
#include "createPhi.H"
threePhaseMixture threePhaseProperties(U, phi);
immiscibleIncompressibleThreePhaseMixture mixture(U, phi);
volScalarField& alpha1(threePhaseProperties.alpha1());
volScalarField& alpha2(threePhaseProperties.alpha2());
volScalarField& alpha3(threePhaseProperties.alpha3());
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
volScalarField& alpha3(mixture.alpha3());
const dimensionedScalar& rho1 = threePhaseProperties.rho1();
const dimensionedScalar& rho2 = threePhaseProperties.rho2();
const dimensionedScalar& rho3 = threePhaseProperties.rho3();
const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
const dimensionedScalar& rho3 = mixture.rho3();
dimensionedScalar D23(threePhaseProperties.lookup("D23"));
dimensionedScalar D23(mixture.lookup("D23"));
// Need to store rho for ddt(rho, U)
volScalarField rho
@ -73,14 +73,10 @@
);
// Construct interface from alpha distribution
threePhaseInterfaceProperties interface(threePhaseProperties);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, threePhaseProperties)
incompressible::turbulenceModel::New(U, phi, mixture)
);

46
applications/solvers/multiphase/interFoam/interMixingFoam/immiscibleIncompressibleThreePhaseMixture/immiscibleIncompressibleThreePhaseMixture.C Normal file
View file

@ -0,0 +1,46 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 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/>.
\*---------------------------------------------------------------------------*/
#include "immiscibleIncompressibleThreePhaseMixture.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::immiscibleIncompressibleThreePhaseMixture::
immiscibleIncompressibleThreePhaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi
)
:
incompressibleThreePhaseMixture(U, phi),
threePhaseInterfaceProperties
(
static_cast<incompressibleThreePhaseMixture&>(*this)
)
{}
// ************************************************************************* //

92
applications/solvers/multiphase/interFoam/interMixingFoam/immiscibleIncompressibleThreePhaseMixture/immiscibleIncompressibleThreePhaseMixture.H Normal file
View file

@ -0,0 +1,92 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 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/>.
Class
Foam::immiscibleIncompressibleThreePhaseMixture
Description
An immiscible incompressible two-phase mixture transport model
SourceFiles
immiscibleIncompressibleThreePhaseMixture.C
\*---------------------------------------------------------------------------*/
#ifndef immiscibleIncompressibleThreePhaseMixture_H
#define immiscibleIncompressibleThreePhaseMixture_H
#include "incompressibleThreePhaseMixture.H"
#include "threePhaseInterfaceProperties.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class immiscibleIncompressibleThreePhaseMixture Declaration
\*---------------------------------------------------------------------------*/
class immiscibleIncompressibleThreePhaseMixture
:
public incompressibleThreePhaseMixture,
public threePhaseInterfaceProperties
{
public:
// Constructors
//- Construct from components
immiscibleIncompressibleThreePhaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi
);
//- Destructor
virtual ~immiscibleIncompressibleThreePhaseMixture()
{}
// Member Functions
//- Correct the transport and interface properties
virtual void correct()
{
incompressibleThreePhaseMixture::correct();
threePhaseInterfaceProperties::correct();
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

19
applications/solvers/multiphase/interFoam/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C → applications/solvers/multiphase/interFoam/interMixingFoam/incompressibleThreePhaseMixture/incompressibleThreePhaseMixture.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,7 +23,7 @@ License
\*---------------------------------------------------------------------------*/
#include "threePhaseMixture.H"
#include "incompressibleThreePhaseMixture.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceFields.H"
#include "fvc.H"
@ -31,7 +31,7 @@ License
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
//- Calculate and return the laminar viscosity
void Foam::threePhaseMixture::calcNu()
void Foam::incompressibleThreePhaseMixture::calcNu()
{
nuModel1_->correct();
nuModel2_->correct();
@ -44,7 +44,7 @@ void Foam::threePhaseMixture::calcNu()
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::threePhaseMixture::threePhaseMixture
Foam::incompressibleThreePhaseMixture::incompressibleThreePhaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi
@ -163,7 +163,8 @@ Foam::threePhaseMixture::threePhaseMixture
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField> Foam::threePhaseMixture::mu() const
Foam::tmp<Foam::volScalarField>
Foam::incompressibleThreePhaseMixture::mu() const
{
return tmp<volScalarField>
(
@ -178,7 +179,8 @@ Foam::tmp<Foam::volScalarField> Foam::threePhaseMixture::mu() const
}
Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::muf() const
Foam::tmp<Foam::surfaceScalarField>
Foam::incompressibleThreePhaseMixture::muf() const
{
surfaceScalarField alpha1f(fvc::interpolate(alpha1_));
surfaceScalarField alpha2f(fvc::interpolate(alpha2_));
@ -197,7 +199,8 @@ Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::muf() const
}
Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::nuf() const
Foam::tmp<Foam::surfaceScalarField>
Foam::incompressibleThreePhaseMixture::nuf() const
{
surfaceScalarField alpha1f(fvc::interpolate(alpha1_));
surfaceScalarField alpha2f(fvc::interpolate(alpha2_));
@ -218,7 +221,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::nuf() const
}
bool Foam::threePhaseMixture::read()
bool Foam::incompressibleThreePhaseMixture::read()
{
if (transportModel::read())
{

18
applications/solvers/multiphase/interFoam/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H → applications/solvers/multiphase/interFoam/interMixingFoam/incompressibleThreePhaseMixture/incompressibleThreePhaseMixture.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,17 +22,17 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
threePhaseMixture
incompressibleThreePhaseMixture
Description
SourceFiles
threePhaseMixture.C
incompressibleThreePhaseMixture.C
\*---------------------------------------------------------------------------*/
#ifndef threePhaseMixture_H
#define threePhaseMixture_H
#ifndef incompressibleThreePhaseMixture_H
#define incompressibleThreePhaseMixture_H
#include "incompressible/transportModel/transportModel.H"
#include "IOdictionary.H"
@ -46,10 +46,10 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class threePhaseMixture Declaration
Class incompressibleThreePhaseMixture Declaration
\*---------------------------------------------------------------------------*/
class threePhaseMixture
class incompressibleThreePhaseMixture
:
public IOdictionary,
public transportModel
@ -89,7 +89,7 @@ public:
// Constructors
//- Construct from components
threePhaseMixture
incompressibleThreePhaseMixture
(
const volVectorField& U,
const surfaceScalarField& phi
@ -97,7 +97,7 @@ public:
//- Destructor
~threePhaseMixture()
~incompressibleThreePhaseMixture()
{}

9
applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
View file

@ -33,8 +33,7 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "subCycle.H"
#include "threePhaseMixture.H"
#include "threePhaseInterfaceProperties.H"
#include "immiscibleIncompressibleThreePhaseMixture.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
@ -78,13 +77,9 @@ int main(int argc, char *argv[])
while (pimple.loop())
{
#include "alphaControls.H"
threePhaseProperties.correct();
#include "alphaEqnsSubCycle.H"
interface.correct();
#define twoPhaseProperties threePhaseProperties
mixture.correct();
#include "UEqn.H"

26
applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,6 +37,7 @@ Description
#include "surfaceInterpolate.H"
#include "fvcDiv.H"
#include "fvcGrad.H"
#include "fvcSnGrad.H"
// * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
@ -166,7 +167,7 @@ void Foam::threePhaseInterfaceProperties::calculateK()
Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
(
const threePhaseMixture& mixture
const incompressibleThreePhaseMixture& mixture
)
:
mixture_(mixture),
@ -212,4 +213,25 @@ Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::tmp<Foam::surfaceScalarField>
Foam::threePhaseInterfaceProperties::surfaceTensionForce() const
{
return fvc::interpolate(sigmaK())*fvc::snGrad(mixture_.alpha1());
}
Foam::tmp<Foam::volScalarField>
Foam::threePhaseInterfaceProperties::nearInterface() const
{
return max
(
pos(mixture_.alpha1() - 0.01)*pos(0.99 - mixture_.alpha1()),
pos(mixture_.alpha2() - 0.01)*pos(0.99 - mixture_.alpha2())
);
}
// ************************************************************************* //

17
applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,7 +37,7 @@ SourceFiles
#ifndef threePhaseInterfaceProperties_H
#define threePhaseInterfaceProperties_H
#include "threePhaseMixture.H"
#include "incompressibleThreePhaseMixture.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -53,7 +53,7 @@ class threePhaseInterfaceProperties
{
// Private data
const threePhaseMixture& mixture_;
const incompressibleThreePhaseMixture& mixture_;
//- Compression coefficient
scalar cAlpha_;
@ -98,7 +98,10 @@ public:
// Constructors
//- Construct from volume fraction field alpha and IOdictionary
threePhaseInterfaceProperties(const threePhaseMixture& mixture);
threePhaseInterfaceProperties
(
const incompressibleThreePhaseMixture& mixture
);
// Member Functions
@ -138,6 +141,12 @@ public:
return sigma()*K_;
}
tmp<surfaceScalarField> surfaceTensionForce() const;
//- Indicator of the proximity of the interface
// Field values are 1 near and 0 away for the interface.
tmp<volScalarField> nearInterface() const;
void correct()
{
calculateK();

6
applications/solvers/multiphase/interFoam/pEqn.H
View file

@ -11,13 +11,13 @@
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
fvOptions.makeRelative(phiHbyA);
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);
@ -30,7 +30,7 @@
p_rgh.boundaryField(),
(
phiHbyA.boundaryField()
- (mesh.Sf().boundaryField() & U.boundaryField())
- fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);

5
applications/solvers/multiphase/interFoam/porousInterFoam/Make/options
View file

@ -4,6 +4,7 @@ EXE_INC = \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
@ -11,9 +12,7 @@ EXE_INC = \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-ltwoPhaseMixture \
-linterfaceProperties \
-ltwoPhaseProperties \
-limmiscibleIncompressibleTwoPhaseMixture \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \

6
applications/solvers/multiphase/interFoam/porousInterFoam/UEqn.H
View file

@ -1,12 +1,12 @@
surfaceScalarField muEff
(
"muEff",
twoPhaseProperties.muf()
mixture.muf()
+ fvc::interpolate(rho*turbulence->nut())
);
// Calculate and cache mu for the porous media
volScalarField mu(twoPhaseProperties.mu());
volScalarField mu(mixture.mu());
fvVectorMatrix UEqn
(
@ -33,7 +33,7 @@
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()

9
applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
View file

@ -41,8 +41,7 @@ Description
#include "fvCFD.H"
#include "CMULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "incompressibleTwoPhaseMixture.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulenceModel.H"
#include "IOporosityModelList.H"
#include "pimpleControl.H"
@ -87,11 +86,9 @@ int main(int argc, char *argv[])
while (pimple.loop())
{
#include "alphaControls.H"
twoPhaseProperties.correct();
#include "alphaEqnSubCycle.H"
interface.correct();
mixture.correct();
#include "UEqn.H"

2
applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
View file

@ -17,7 +17,7 @@
fvc::reconstruct
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
interface.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()

2
applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
View file

@ -5,7 +5,7 @@
surfaceScalarField phir("phir", phic*interface.nHatf());
Pair<tmp<volScalarField> > vDotAlphal =
twoPhaseProperties->vDotAlphal();
mixture->vDotAlphal();
const volScalarField& vDotcAlphal = vDotAlphal[0]();
const volScalarField& vDotvAlphal = vDotAlphal[1]();
const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);

16
applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
View file

@ -30,15 +30,15 @@
Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr<phaseChangeTwoPhaseMixture> twoPhaseProperties =
autoPtr<phaseChangeTwoPhaseMixture> mixture =
phaseChangeTwoPhaseMixture::New(U, phi);
volScalarField& alpha1(twoPhaseProperties->alpha1());
volScalarField& alpha2(twoPhaseProperties->alpha2());
volScalarField& alpha1(mixture->alpha1());
volScalarField& alpha2(mixture->alpha2());
const dimensionedScalar& rho1 = twoPhaseProperties->rho1();
const dimensionedScalar& rho2 = twoPhaseProperties->rho2();
const dimensionedScalar& pSat = twoPhaseProperties->pSat();
const dimensionedScalar& rho1 = mixture->rho1();
const dimensionedScalar& rho2 = mixture->rho2();
const dimensionedScalar& pSat = mixture->pSat();
// Need to store rho for ddt(rho, U)
@ -58,12 +58,12 @@
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties());
interfaceProperties interface(alpha1, U, mixture());
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties())
incompressible::turbulenceModel::New(U, phi, mixture())
);

2
applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
View file

@ -154,7 +154,7 @@ int main(int argc, char *argv[])
dimensionedScalar("0", dimMass/dimTime, 0)
);
twoPhaseProperties->correct();
mixture->correct();
#include "alphaEqnSubCycle.H"
interface.correct();

4
applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
View file

@ -22,7 +22,7 @@
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
interface.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);
@ -39,7 +39,7 @@
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);
Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
Pair<tmp<volScalarField> > vDotP = mixture->vDotP();
const volScalarField& vDotcP = vDotP[0]();
const volScalarField& vDotvP = vDotP[1]();

2
applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
View file

@ -100,7 +100,7 @@ int main(int argc, char *argv[])
dimensionedScalar("0", dimMass/dimTime, 0)
);
twoPhaseProperties->correct();
mixture->correct();
#include "alphaEqnSubCycle.H"
interface.correct();

4
applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
View file

@ -16,7 +16,7 @@
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
interface.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf()
);
@ -33,7 +33,7 @@
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);
Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
Pair<tmp<volScalarField> > vDotP = mixture->vDotP();
const volScalarField& vDotcP = vDotP[0]();
const volScalarField& vDotvP = vDotP[1]();

3
applications/solvers/multiphase/multiphaseInterFoam/Allwclean
View file

@ -2,7 +2,8 @@
cd ${0%/*} || exit 1 # run from this directory
set -x
wclean libso multiphaseMixture
wclean
wclean MRFMultiphaseInterFoam
wclean multiphaseInterDyMFoam
# ----------------------------------------------------------------- end-of-file

2
applications/solvers/multiphase/multiphaseInterFoam/Allwmake
View file

@ -4,6 +4,6 @@ set -x
wmake libso multiphaseMixture
wmake
wmake MRFMultiphaseInterFoam
wmake multiphaseInterDyMFoam
# ----------------------------------------------------------------- end-of-file

3
applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/Make/files
View file

@ -1,3 +0,0 @@
MRFMultiphaseInterFoam.C
EXE = $(FOAM_APPBIN)/MRFMultiphaseInterFoam

19
applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/Make/options
View file

@ -1,19 +0,0 @@
EXE_INC = \
-I.. \
-I../../interFoam \
-I../../interFoam/MRFInterFoam \
-I../multiphaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lmultiphaseInterFoam \
-linterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume

35
applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/UEqn.H
View file

@ -1,35 +0,0 @@
surfaceScalarField muEff
(
"muEff",
mixture.muf()
+ fvc::interpolate(rho*turbulence->nut())
);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(mixture.rhoPhi(), U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
);
mrfZones.addCoriolis(rho, UEqn);
UEqn.relax();
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
}

68
applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
View file

@ -1,68 +0,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);
mrfZones.makeRelative(phiHbyA);
surfaceScalarField phig
(
- ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
);
phiHbyA += phig;
// Update the fixedFluxPressure BCs to ensure flux consistency
setSnGrad<fixedFluxPressureFvPatchScalarField>
(
p_rgh.boundaryField(),
(
phiHbyA.boundaryField()
- mrfZones.relative(mesh.Sf().boundaryField() & U.boundaryField())
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
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();
}
}
#include "continuityErrs.H"
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

26
applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/zonePhaseVolumes.H
View file

@ -1,26 +0,0 @@
{
const scalarField& V = mesh.V();
forAll(mesh.cellZones(), czi)
{
const labelList& cellLabels = mesh.cellZones()[czi];
forAllConstIter(PtrDictionary<phase>, mixture.phases(), iter)
{
const volScalarField& alpha = iter();
scalar phaseVolume = 0;
forAll(cellLabels, cli)
{
label celli = cellLabels[cli];
phaseVolume += alpha[celli]*V[celli];
}
reduce(phaseVolume, sumOp<scalar>());
Info<< alpha.name()
<< " phase volume in zone " << mesh.cellZones()[czi].name()
<< " = " << phaseVolume*1e6 << " ml " << endl;
}
}
}

10
applications/solvers/multiphase/multiphaseInterFoam/Make/options
View file

@ -5,7 +5,10 @@ EXE_INC = \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lmultiphaseInterFoam \
@ -14,4 +17,7 @@ EXE_LIBS = \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume
-lfiniteVolume \
-lfvOptions \
-lmeshTools \
-lsampling

25
applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
View file

@ -1,25 +0,0 @@
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(mixture.rhoPhi(), U)
+ turbulence->divDevRhoReff(rho, U)
);
UEqn.relax();
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
}

3
applications/solvers/multiphase/multiphaseInterFoam/createFields.H
View file

@ -44,6 +44,7 @@
);
rho.oldTime();
const surfaceScalarField& rhoPhi(mixture.rhoPhi());
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
@ -90,3 +91,5 @@
pRefValue - getRefCellValue(p, pRefCell)
);
}
fv::IOoptionList fvOptions(mesh);

3
applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/Make/files Normal file
View file

@ -0,0 +1,3 @@
multiphaseInterDyMFoam.C
EXE = $(FOAM_APPBIN)/multiphaseInterDyMFoam

16
applications/solvers/multiphase/interFoam/MRFInterFoam/Make/options → applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/Make/options
View file

@ -1,24 +1,30 @@
EXE_INC = \
-I.. \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I../../interFoam/interDyMFoam \
-I../../interFoam \
-I../multiphaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-ltwoPhaseMixture \
-lmultiphaseInterFoam \
-linterfaceProperties \
-ltwoPhaseProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \
-lmeshTools \
-ldynamicMesh \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-lfvOptions \
-lmeshTools \
-lsampling

73
applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C → applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/multiphaseInterDyMFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
MRFMultiphaseInterFoam
multiphaseInterFoam
Description
Solver for n incompressible fluids which captures the interfaces and
@ -33,10 +33,11 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "multiphaseMixture.H"
#include "turbulenceModel.H"
#include "IOMRFZoneList.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -45,16 +46,31 @@ int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createMRFZones.H"
#include "readTimeControls.H"
#include "createPrghCorrTypes.H"
volScalarField rAU
(
IOobject
(
"rAU",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
);
#include "correctPhi.H"
#include "createUf.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
@ -64,22 +80,57 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readTimeControls.H"
#include "readControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
mixture.solve();
rho = mixture.rho();
#include "zonePhaseVolumes.H"
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
mesh.update();
if (mesh.changing())
{
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
gh = g & mesh.C();
ghf = g & mesh.Cf();
}
if (mesh.changing() && correctPhi)
{
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & Uf;
#include "correctPhi.H"
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
mixture.correct();
}
if (mesh.changing() && checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
mixture.solve();
rho = mixture.rho();
#include "UEqn.H"
// --- Pressure corrector loop

3
applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -36,6 +36,7 @@ Description
#include "multiphaseMixture.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

12
applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
View file

@ -298,10 +298,7 @@ Foam::multiphaseMixture::surfaceTensionForce() const
void Foam::multiphaseMixture::solve()
{
forAllIter(PtrDictionary<phase>, phases_, iter)
{
iter().correct();
}
correct();
const Time& runTime = mesh_.time();
@ -347,7 +344,12 @@ void Foam::multiphaseMixture::solve()
void Foam::multiphaseMixture::correct()
{}
{
forAllIter(PtrDictionary<phase>, phases_, iter)
{
iter().correct();
}
}
Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv

70
applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
View file

@ -1,70 +0,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
(
(
mixture.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())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
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();
}
}
#include "continuityErrs.H"
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

3
applications/solvers/multiphase/settlingFoam/Make/files
View file

@ -1,3 +0,0 @@
settlingFoam.C
EXE = $(FOAM_APPBIN)/settlingFoam

32
applications/solvers/multiphase/settlingFoam/UEqn.H
View file

@ -1,32 +0,0 @@
// Solve the Momentum equation
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ fvc::div
(
(Alpha/(scalar(1.001) - Alpha))*(sqr(rhoc)/rho)*Vdj*Vdj,
"div(phiVdj,Vdj)"
)
- fvm::laplacian(muEff, U, "laplacian(muEff,U)")
- (fvc::grad(U) & fvc::grad(muEff))
);
UEqn.relax();
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
}

34
applications/solvers/multiphase/settlingFoam/alphaEqn.H
View file

@ -1,34 +0,0 @@
{
surfaceScalarField phiAlpha
(
IOobject
(
"phiAlpha",
runTime.timeName(),
mesh
),
phi + rhoc*(mesh.Sf() & fvc::interpolate(Vdj))
);
fvScalarMatrix AlphaEqn
(
fvm::ddt(rho, Alpha)
+ fvm::div(phiAlpha, Alpha)
- fvm::laplacian(mut, Alpha)
);
AlphaEqn.relax();
AlphaEqn.solve();
Info<< "Solid phase fraction = "
<< Alpha.weightedAverage(mesh.V()).value()
<< " Min(Alpha) = " << min(Alpha).value()
<< " Max(Alpha) = " << max(Alpha).value()
<< endl;
Alpha.min(1.0);
Alpha.max(0.0);
rho == rhoc/(scalar(1) + (rhoc/rhod - 1.0)*Alpha);
alpha == rho*Alpha/rhod;
}

20
applications/solvers/multiphase/settlingFoam/calcVdj.H
View file

@ -1,20 +0,0 @@
if (VdjModel == "general")
{
Vdj = V0*
(
exp(-a*max(alpha - alphaMin, scalar(0)))
- exp(-a1*max(alpha - alphaMin, scalar(0)))
);
}
else if (VdjModel == "simple")
{
Vdj = V0*pow(10.0, -a*alpha);
}
else
{
FatalErrorIn(args.executable())
<< "Unknown VdjModel : " << VdjModel
<< abort(FatalError);
}
Vdj.correctBoundaryConditions();

21
applications/solvers/multiphase/settlingFoam/compressibleContinuityErrs.H
View file

@ -1,21 +0,0 @@
scalar sumLocalContErr =
runTime.deltaTValue()*
mag
(
fvc::ddt(rho)
+ fvc::div(phi)
)().weightedAverage(rho*mesh.V()).value();
scalar globalContErr =
runTime.deltaTValue()*
(
fvc::ddt(rho)
+ fvc::div(phi)
)().weightedAverage(rho*mesh.V()).value();
cumulativeContErr += globalContErr;
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;

39
applications/solvers/multiphase/settlingFoam/correctViscosity.H
View file

@ -1,39 +0,0 @@
{
mul = muc +
plasticViscosity
(
plasticViscosityCoeff,
plasticViscosityExponent,
alpha
);
if (BinghamPlastic)
{
volScalarField tauy = yieldStress
(
yieldStressCoeff,
yieldStressExponent,
yieldStressOffset,
alpha
);
mul =
tauy/
(
mag(fvc::grad(U))
+ 1.0e-4*
(
tauy
+ dimensionedScalar
(
"deltaTauy",
tauy.dimensions(),
1.0e-15
)
)/mul
)
+ mul;
}
mul = min(mul, muMax);
}

383
applications/solvers/multiphase/settlingFoam/createFields.H
View file

@ -1,383 +0,0 @@
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 Alpha\n" << endl;
volScalarField Alpha
(
IOobject
(
"Alpha",
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
);
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
dimensionedScalar rhoc(transportProperties.lookup("rhoc"));
dimensionedScalar rhod(transportProperties.lookup("rhod"));
dimensionedScalar muc(transportProperties.lookup("muc"));
dimensionedScalar muMax(transportProperties.lookup("muMax"));
dimensionedScalar plasticViscosityCoeff
(
transportProperties.lookup("plasticViscosityCoeff")
);
dimensionedScalar plasticViscosityExponent
(
transportProperties.lookup("plasticViscosityExponent")
);
dimensionedScalar yieldStressCoeff
(
transportProperties.lookup("yieldStressCoeff")
);
dimensionedScalar yieldStressExponent
(
transportProperties.lookup("yieldStressExponent")
);
dimensionedScalar yieldStressOffset
(
transportProperties.lookup("yieldStressOffset")
);
Switch BinghamPlastic(transportProperties.lookup("BinghamPlastic"));
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
rhoc/(scalar(1) + (rhoc/rhod - 1.0)*Alpha)
);
volScalarField alpha
(
IOobject
(
"alpha",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
rho*Alpha/rhod
);
#include "compressibleCreatePhi.H"
Info<< "Calculating field mul\n" << endl;
volScalarField mul
(
IOobject
(
"mul",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
muc
+ plasticViscosity
(
plasticViscosityCoeff,
plasticViscosityExponent,
Alpha
)
);
Info<< "Initialising field Vdj\n" << endl;
volVectorField Vdj
(
IOobject
(
"Vdj",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector("0.0", U.dimensions(), vector::zero),
U.boundaryField().types()
);
Info<< "Selecting Drift-Flux model " << endl;
const word VdjModel(transportProperties.lookup("VdjModel"));
Info<< tab << VdjModel << " selected\n" << endl;
const dictionary& VdjModelCoeffs
(
transportProperties.subDict(VdjModel + "Coeffs")
);
dimensionedVector V0(VdjModelCoeffs.lookup("V0"));
dimensionedScalar a(VdjModelCoeffs.lookup("a"));
dimensionedScalar a1(VdjModelCoeffs.lookup("a1"));
dimensionedScalar alphaMin(VdjModelCoeffs.lookup("alphaMin"));
IOdictionary RASProperties
(
IOobject
(
"RASProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
Switch turbulence(RASProperties.lookup("turbulence"));
dictionary kEpsilonDict(RASProperties.subDictPtr("kEpsilonCoeffs"));
dimensionedScalar Cmu
(
dimensionedScalar::lookupOrAddToDict
(
"Cmu",
kEpsilonDict,
0.09
)
);
dimensionedScalar C1
(
dimensionedScalar::lookupOrAddToDict
(
"C1",
kEpsilonDict,
1.44
)
);
dimensionedScalar C2
(
dimensionedScalar::lookupOrAddToDict
(
"C2",
kEpsilonDict,
1.92
)
);
dimensionedScalar C3
(
dimensionedScalar::lookupOrAddToDict
(
"C3",
kEpsilonDict,
0.85
)
);
dimensionedScalar sigmak
(
dimensionedScalar::lookupOrAddToDict
(
"sigmak",
kEpsilonDict,
1.0
)
);
dimensionedScalar sigmaEps
(
dimensionedScalar::lookupOrAddToDict
(
"sigmaEps",
kEpsilonDict,
1.3
)
);
dictionary wallFunctionDict(RASProperties.subDictPtr("wallFunctionCoeffs"));
dimensionedScalar kappa
(
dimensionedScalar::lookupOrAddToDict
(
"kappa",
wallFunctionDict,
0.41
)
);
dimensionedScalar E
(
dimensionedScalar::lookupOrAddToDict
(
"E",
wallFunctionDict,
9.8
)
);
if (RASProperties.lookupOrDefault("printCoeffs", false))
{
Info<< "kEpsilonCoeffs" << kEpsilonDict << nl
<< "wallFunctionCoeffs" << wallFunctionDict << endl;
}
nearWallDist y(mesh);
Info<< "Reading field k\n" << endl;
volScalarField k
(
IOobject
(
"k",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field epsilon\n" << endl;
volScalarField epsilon
(
IOobject
(
"epsilon",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Calculating field mut\n" << endl;
volScalarField mut
(
IOobject
(
"mut",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Cmu*rho*sqr(k)/epsilon
);
Info<< "Calculating field muEff\n" << endl;
volScalarField muEff
(
IOobject
(
"muEff",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mut + mul
);
Info<< "Calculating field (g.h)f\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("PIMPLE"),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}

78
applications/solvers/multiphase/settlingFoam/kEpsilon.H
View file

@ -1,78 +0,0 @@
if (turbulence)
{
if (mesh.changing())
{
y.correct();
}
dimensionedScalar k0("k0", k.dimensions(), 0);
dimensionedScalar kMin("kMin", k.dimensions(), SMALL);
dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0);
dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
volScalarField divU(fvc::div(phi/fvc::interpolate(rho)));
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField G(2*mut*(tgradU() && dev(symm(tgradU()))));
tgradU.clear();
volScalarField Gcoef
(
Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin)
);
#include "wallFunctions.H"
// Dissipation equation
fvScalarMatrix epsEqn
(
fvm::ddt(rho, epsilon)
+ fvm::div(phi, epsilon)
- fvm::laplacian
(
mut/sigmaEps + muc, epsilon,
"laplacian(DepsilonEff,epsilon)"
)
==
C1*G*epsilon/(k + kMin)
- fvm::SuSp(C1*(1.0 - C3)*Gcoef + (2.0/3.0*C1)*rho*divU, epsilon)
- fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon)
);
#include "wallDissipation.H"
epsEqn.relax();
epsEqn.solve();
bound(epsilon, epsilon0);
// Turbulent kinetic energy equation
fvScalarMatrix kEqn
(
fvm::ddt(rho, k)
+ fvm::div(phi, k)
- fvm::laplacian
(
mut/sigmak + muc, k,
"laplacian(DkEff,k)"
)
==
G
- fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k)
- fvm::Sp(rho*epsilon/(k + kMin), k)
);
kEqn.relax();
kEqn.solve();
bound(k, k0);
//- Re-calculate viscosity
mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin);
#include "wallViscosity.H"
}
muEff = mut + mul;

70
applications/solvers/multiphase/settlingFoam/pEqn.H
View file

@ -1,70 +0,0 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
surfaceScalarField phiHbyA
(
"phiHbyA",
(
(fvc::interpolate(rho*HbyA) & mesh.Sf())
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
)
);
surfaceScalarField phig
(
- ghf*fvc::snGrad(rho)*rhorAUf*mesh.magSf()
);
phiHbyA += phig;
// Update the fixedFluxPressure BCs to ensure flux consistency
setSnGrad<fixedFluxPressureFvPatchScalarField>
(
p_rgh.boundaryField(),
(
phiHbyA.boundaryField()
- (mesh.Sf().boundaryField() & U.boundaryField())
*rho.boundaryField()
)/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rhorAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
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())/rhorAUf);
U.correctBoundaryConditions();
}
}
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
}

17
applications/solvers/multiphase/settlingFoam/plasticViscosity.H
View file

@ -1,17 +0,0 @@
volScalarField plasticViscosity
(
const dimensionedScalar& plasticViscosityCoeff,
const dimensionedScalar& plasticViscosityExponent,
const volScalarField& alpha
)
{
tmp<volScalarField> tfld
(
plasticViscosityCoeff*
(
pow(10.0, plasticViscosityExponent*alpha + SMALL) - scalar(1)
)
);
return tfld();
}

109
applications/solvers/multiphase/settlingFoam/settlingFoam.C
View file

@ -1,109 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 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
settlingFoam
Description
Solver for 2 incompressible fluids for simulating the settling of the
dispersed phase.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "nearWallDist.H"
#include "wallFvPatch.H"
#include "bound.H"
#include "Switch.H"
#include "plasticViscosity.H"
#include "yieldStress.H"
#include "pimpleControl.H"
#include "fixedFluxPressureFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "rhoEqn.H"
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "calcVdj.H"
#include "UEqn.H"
#include "alphaEqn.H"
#include "correctViscosity.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
#include "kEpsilon.H"
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

85
applications/solvers/multiphase/settlingFoam/wallFunctions.H
View file

@ -1,85 +0,0 @@
{
labelList cellBoundaryFaceCount(epsilon.size(), 0);
const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
const scalar Cmu75 = ::pow(Cmu.value(), 0.75);
const scalar kappa_ = kappa.value();
const scalar muc_ = muc.value();
const fvPatchList& patches = mesh.boundary();
//- Initialise the near-wall P field to zero
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
epsilon[faceCelli] = 0.0;
G[faceCelli] = 0.0;
}
}
}
//- Accumulate the wall face contributions to epsilon and G
// Increment cellBoundaryFaceCount for each face for averaging
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
const scalarField& mutw = mut.boundaryField()[patchi];
scalarField magFaceGradU
(
mag(U.boundaryField()[patchi].snGrad())
);
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
// For corner cells (with two boundary or more faces),
// epsilon and G in the near-wall cell are calculated
// as an average
cellBoundaryFaceCount[faceCelli]++;
epsilon[faceCelli] +=
Cmu75*::pow(k[faceCelli], 1.5)
/(kappa_*y[patchi][facei]);
G[faceCelli] +=
(mutw[facei] + muc_)
*magFaceGradU[facei]
*Cmu25*::sqrt(k[faceCelli])
/(kappa_*y[patchi][facei]);
}
}
}
// perform the averaging
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
cellBoundaryFaceCount[faceCelli] = 1;
}
}
}
}

Some files were not shown because too many files have changed in this diff Show more