compressibleTwoPhaseEulerFoam: Further updates to include new thermo structure

This commit is contained in:
Henry
2013-01-17 18:02:55 +00:00
parent 4737661d6c
commit fe8bcb4905
37 changed files with 156 additions and 130 deletions

View File

@ -14,15 +14,13 @@
volScalarField& he1 = thermo1.he();
volScalarField& he2 = thermo2.he();
Info<< max(he1) << min(he1) << endl;
fvScalarMatrix he1Eqn
(
fvm::ddt(alpha1, he1)
+ fvm::div(alphaPhi1, he1)
- fvm::laplacian(k1, he1)
==
heatTransferCoeff*(he1/thermo1.Cp())/rho1
heatTransferCoeff*(thermo1.he(p, thermo2.T())/thermo1.Cp())/rho1
- fvm::Sp(heatTransferCoeff/thermo1.Cp()/rho1, he1)
+ alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))
);
@ -33,7 +31,7 @@
+ fvm::div(alphaPhi2, he2)
- fvm::laplacian(k2, he2)
==
heatTransferCoeff*(he2/thermo2.Cp())/rho2
heatTransferCoeff*(thermo2.he(p, thermo1.T())/thermo2.Cp())/rho2
- fvm::Sp(heatTransferCoeff/thermo2.Cp()/rho2, he2)
+ alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))
);

View File

@ -1,9 +1,9 @@
surfaceScalarField alphaPhi1("alphaPhi", phi1);
surfaceScalarField alphaPhi2("alphaPhi", phi2);
surfaceScalarField alphaPhi1("alphaPhi" + phase1Name, phi1);
surfaceScalarField alphaPhi2("alphaPhi" + phase2Name, phi2);
{
word scheme("div(phi,alpha)");
word schemer("div(phir,alpha)");
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phir,alpha)");
surfaceScalarField phic("phic", phi);
surfaceScalarField phir("phir", phi1 - phi2);
@ -56,15 +56,55 @@ surfaceScalarField alphaPhi2("alphaPhi", phi2);
}
}
dimensionedScalar totalDeltaT = runTime.deltaT();
if (nAlphaSubCycles > 1)
{
alphaPhi1 = dimensionedScalar("0", alphaPhi1.dimensions(), 0);
}
fvScalarMatrix alpha1Eqn
for
(
fvm::ddt(alpha1)
+ fvm::div(phic, alpha1, scheme)
+ fvm::div(-fvc::flux(-phir, alpha2, schemer), alpha1, schemer)
==
fvm::Sp(Sp, alpha1) + Su
);
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
surfaceScalarField alphaPhic1
(
fvc::flux
(
phic,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
alpha1,
alpharScheme
)
);
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phi,
alphaPhic1,
Sp,
Su,
(g0.value() > 0 ? alphaMax : 1),
0
);
if (nAlphaSubCycles > 1)
{
alphaPhi1 += (runTime.deltaT()/totalDeltaT)*alphaPhic1;
}
else
{
alphaPhi1 = alphaPhic1;
}
}
if (g0.value() > 0.0)
{
@ -74,30 +114,25 @@ surfaceScalarField alphaPhi2("alphaPhi", phi2);
fvc::interpolate((1.0/rho1)*rAU1)
*g0*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
// ppMagf =
// fvc::interpolate((1.0/rho1)*rAU1)
// *fvc::interpolate
// (
// g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
// );
alpha1Eqn -= fvm::laplacian
fvScalarMatrix alpha1Eqn
(
alpha1f*ppMagf,
alpha1,
"laplacian(alphaPpMag,alpha1)"
fvm::ddt(alpha1) - fvc::ddt(alpha1)
- fvm::laplacian
(
alpha1f*ppMagf,
alpha1,
"laplacian(alpha1PpMag,alpha1)"
)
);
alpha1Eqn.relax();
alpha1Eqn.solve();
#include "packingLimiter.H"
alphaPhi1 += alpha1Eqn.flux();
}
alpha1Eqn.relax();
alpha1Eqn.solve();
//***HGW temporary boundedness-fix pending the introduction of MULES
alpha1 = max(min(alpha1, scalar(1)), scalar(0));
#include "packingLimiter.H"
alphaPhi1 = alpha1Eqn.flux();
alphaPhi2 = phi - alphaPhi1;
alpha2 = scalar(1) - alpha1;

View File

@ -31,6 +31,8 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "rhoThermo.H"
#include "nearWallDist.H"
#include "wallFvPatch.H"

View File

@ -16,14 +16,14 @@
(
transportProperties.found("phases")
? wordList(transportProperties.lookup("phases"))[0]
: "phase1"
: "1"
);
word phase2Name
(
transportProperties.found("phases")
? wordList(transportProperties.lookup("phases"))[1]
: "phase2"
: "2"
);
autoPtr<phaseModel> phase1 = phaseModel::New

View File

@ -73,7 +73,7 @@ Foam::kineticTheoryModels::conductivityModels::Gidaspow::kappa
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const

View File

@ -77,7 +77,7 @@ public:
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const;

View File

@ -76,7 +76,7 @@ Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const

View File

@ -82,7 +82,7 @@ public:
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const;

View File

@ -73,7 +73,7 @@ Foam::kineticTheoryModels::conductivityModels::Syamlal::kappa
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const

View File

@ -77,7 +77,7 @@ public:
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const;

View File

@ -109,7 +109,7 @@ public:
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const = 0;

View File

@ -71,7 +71,7 @@ Foam::kineticTheoryModels::granularPressureModels::Lun::granularPressureCoeff
(
const volScalarField& alpha1,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const
{
@ -87,7 +87,7 @@ granularPressureCoeffPrime
const volScalarField& alpha1,
const volScalarField& g0,
const volScalarField& g0prime,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const
{

View File

@ -76,7 +76,7 @@ public:
(
const volScalarField& alpha1,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const;
@ -85,7 +85,7 @@ public:
const volScalarField& alpha1,
const volScalarField& g0,
const volScalarField& g0prime,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const;
};

View File

@ -74,7 +74,7 @@ granularPressureCoeff
(
const volScalarField& alpha1,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const
{
@ -90,7 +90,7 @@ granularPressureCoeffPrime
const volScalarField& alpha1,
const volScalarField& g0,
const volScalarField& g0prime,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const
{

View File

@ -76,7 +76,7 @@ public:
(
const volScalarField& alpha1,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const;
@ -85,7 +85,7 @@ public:
const volScalarField& alpha1,
const volScalarField& g0,
const volScalarField& g0prime,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const;
};

View File

@ -109,7 +109,7 @@ public:
(
const volScalarField& alpha1,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const = 0;
@ -119,7 +119,7 @@ public:
const volScalarField& alpha1,
const volScalarField& g0,
const volScalarField& g0prime,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const dimensionedScalar& e
) const = 0;
};

View File

@ -45,8 +45,7 @@ Foam::kineticTheoryModel::kineticTheoryModel
phi1_(phase1.phi()),
draga_(draga),
rho1_(phase1.rho()[0]), //***HGW
nu1_(phase1.nu()()[0]), //***HGW
rho1_(phase1.rho()),
kineticTheoryProperties_
(
@ -203,9 +202,9 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradU1t)
volScalarField da_(phase1_.d());
surfaceScalarField phi(1.5*rho1_*phi1_*fvc::interpolate(alpha1_));
surfaceScalarField phi(1.5*phi1_*fvc::interpolate(rho1_*alpha1_));
volTensorField dU(gradU1t.T()); //fvc::grad(U1_);
volTensorField dU(gradU1t.T());
volSymmTensorField D(symm(dU));
// NB, drag = K*alpha1*alpha2,

View File

@ -64,8 +64,7 @@ class kineticTheoryModel
const dragModel& draga_;
const dimensionedScalar& rho1_;
const dimensionedScalar& nu1_;
const volScalarField& rho1_;
//- dictionary holding the modeling info
IOdictionary kineticTheoryProperties_;

View File

@ -67,7 +67,7 @@ Foam::kineticTheoryModels::viscosityModels::Gidaspow::mu1
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const

View File

@ -76,7 +76,7 @@ public:
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const;

View File

@ -75,7 +75,7 @@ Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::mu1
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const

View File

@ -84,7 +84,7 @@ public:
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const;

View File

@ -67,7 +67,7 @@ Foam::kineticTheoryModels::viscosityModels::Syamlal::mu1
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const

View File

@ -77,7 +77,7 @@ public:
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const;

View File

@ -58,7 +58,7 @@ Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModels::noneViscosity::mu1
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const

View File

@ -75,7 +75,7 @@ public:
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const;

View File

@ -111,7 +111,7 @@ public:
const volScalarField& alpha1,
const volScalarField& Theta,
const volScalarField& g0,
const dimensionedScalar& rho1,
const volScalarField& rho1,
const volScalarField& da,
const dimensionedScalar& e
) const = 0;

View File

@ -51,7 +51,15 @@ Foam::phaseModel::phaseModel
dimensionedScalar("alpha", dimless, 0)
),
name_(phaseName),
phaseDict_(transportProperties.subDict(phaseName)),
phaseDict_
(
transportProperties.subDict
(
phaseName == "1" || phaseName == "2"
? "phase" + phaseName
: word(phaseName)
)
),
thermo_(rhoThermo::New(mesh, phaseName)),
U_
(

View File

@ -1,3 +1,4 @@
#include "readTimeControls.H"
int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));

View File

@ -29,8 +29,6 @@ air
water
{
//R 1e10;
diameterModel constant;
constantCoeffs
{

View File

@ -29,23 +29,23 @@ divSchemes
{
default none;
div(phi,alpha) Gauss limitedLinear01 1;
div(phir,alpha) Gauss limitedLinear01 1;
div(alphaPhi,Uair) Gauss limitedLinearV 1;
div(alphaPhi,Uwater) Gauss limitedLinearV 1;
div(phi,alpha) Gauss limitedLinear01 1;
div(phir,alpha) Gauss limitedLinear01 1;
div(alphaPhiair,Uair) Gauss limitedLinearV 1;
div(alphaPhiwater,Uwater) Gauss limitedLinearV 1;
div(phiair,Uair) Gauss limitedLinearV 1;
div(phiwater,Uwater) Gauss limitedLinearV 1;
div((alphaair*Rc)) Gauss linear;
div((alphawater*Rc)) Gauss linear;
div(alphaPhi,hair) Gauss limitedLinear 1;
div(alphaPhi,hwater) Gauss limitedLinear 1;
div(phiwater,Uwater) Gauss limitedLinearV 1;
div((alphaair*Rc)) Gauss linear;
div((alphawater*Rc)) Gauss linear;
div(alphaPhiair,hair) Gauss limitedLinear 1;
div(alphaPhiwater,hwater) Gauss limitedLinear 1;
div(alphaPhiwater,k) Gauss limitedLinear 1;
div(alphaPhiwater,epsilon) Gauss limitedLinear 1;
div(phi,Theta) Gauss limitedLinear 1;
div(phidair,p) Gauss upwind;
div(alphaPhiwater,epsilon) Gauss limitedLinear 1;
div(phi,Theta) Gauss limitedLinear 1;
div(phidair,p) Gauss upwind;
div(phidwater,p) Gauss upwind;
div(phiair,Kair) Gauss limitedLinear 1;
div(phiwater,Kwater) Gauss limitedLinear 1;
div(phiwater,Kwater) Gauss limitedLinear 1;
}
laplacianSchemes
@ -67,7 +67,6 @@ fluxRequired
{
default no;
p ;
alphaair ;
}

View File

@ -49,8 +49,8 @@ solvers
"h.*"
{
solver PCG; //PBiCG;
preconditioner DIC; //DILU;
solver PBiCG;
preconditioner DILU;
tolerance 1e-8;
relTol 0;
}
@ -93,8 +93,8 @@ PIMPLE
nOuterCorrectors 1;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
nAlphaCorr 2;
correctAlpha no;
nAlphaCorr 1;
nAlphaSubCycles 2;
}
relaxationFactors

View File

@ -22,18 +22,20 @@ boundaryField
{
inlet
{
type zeroGradient;
type fixedFluxPressure;
value $internalField;
}
outlet
{
type fixedValue;
value uniform 1e5;
value $internalField;
}
walls
{
type zeroGradient;
type fixedFluxPressure;
value $internalField;
}
frontAndBackPlanes

View File

@ -38,6 +38,7 @@ FoamFile
frontAndBackPlanes
{
type empty;
inGroups 1(empty);
nFaces 12000;
startFace 12230;
}

View File

@ -17,14 +17,6 @@ FoamFile
phase1
{
rho 2500;
rho0 2500;
R 1e10;
Cp 6000;
nu 1e-06;
kappa 0.613;
diameterModel constant;
constantCoeffs
{
@ -34,14 +26,6 @@ phase1
phase2
{
rho 1.2;
rho0 0;
R 287;
Cp 721;
nu 1.5e-05;
kappa 2.63e-2;
diameterModel constant;
constantCoeffs
{

View File

@ -29,23 +29,23 @@ divSchemes
{
default none;
div(alphaPhi1,U1) Gauss limitedLinearV 1;
div(alphaPhi2,U2) Gauss limitedLinearV 1;
div(phi1,U1) Gauss limitedLinearV 1;
div(phi2,U2) Gauss limitedLinearV 1;
div(alphaPhi1,T1) Gauss limitedLinear 1;
div(alphaPhi2,T2) Gauss limitedLinear 1;
div(alphaPhi2,k) Gauss limitedLinear 1;
div(alphaPhi2,epsilon) Gauss limitedLinear 1;
div(phi,alpha1) Gauss limitedLinear01 1;
div(phir,alpha1) Gauss limitedLinear01 1;
div(phi,Theta) Gauss limitedLinear 1;
div((alpha1*Rc1)) Gauss linear;
div((alpha2*Rc2)) Gauss linear;
div(phid1,p) Gauss upwind;
div(phid2,p) Gauss upwind;
div(phi1,K1) Gauss limitedLinear 1;
div(phi2,K2) Gauss limitedLinear 1;
div(alphaPhi1,U1) Gauss limitedLinearV 1;
div(alphaPhi2,U2) Gauss limitedLinearV 1;
div(phi1,U1) Gauss limitedLinearV 1;
div(phi2,U2) Gauss limitedLinearV 1;
div(alphaPhi1,h1) Gauss limitedLinear 1;
div(alphaPhi2,h2) Gauss limitedLinear 1;
div(alphaPhi2,k) Gauss limitedLinear 1;
div(alphaPhi2,epsilon) Gauss limitedLinear 1;
div(phi,alpha) Gauss limitedLinear01 1;
div(phir,alpha) Gauss limitedLinear01 1;
div(phi,Theta) Gauss limitedLinear 1;
div((alpha1*Rc)) Gauss linear;
div((alpha2*Rc)) Gauss linear;
div(phid1,p) Gauss upwind;
div(phid2,p) Gauss upwind;
div(phi1,K1) Gauss limitedLinear 1;
div(phi2,K2) Gauss limitedLinear 1;
}
laplacianSchemes

View File

@ -47,7 +47,7 @@ solvers
relTol 0;
}
"T.*"
"h.*"
{
solver PBiCG;
preconditioner DILU;
@ -93,8 +93,8 @@ PIMPLE
nOuterCorrectors 1;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
nAlphaCorr 2;
correctAlpha yes;
nAlphaCorr 1;
nAlphaSubCycles 2;
}
relaxationFactors
@ -105,7 +105,7 @@ relaxationFactors
equations
{
"U.*" 1;
"T.*" 1;
"h.*" 1;
"alpha.*" 1;
"Theta.*" 1;
"k.*" 1;