Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
mattijs
2012-03-01 12:52:06 +00:00
255 changed files with 1981 additions and 1950 deletions

View File

@ -49,6 +49,11 @@
singlePhaseTransportModel laminarTransport(U, phi); singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
const volScalarField nu(laminarTransport.nu()); const volScalarField nu(laminarTransport.nu());
volScalarField mu volScalarField mu

View File

@ -1,11 +1,11 @@
{ {
DDtUa = DDtU1 =
fvc::ddt(Ua) fvc::ddt(U1)
+ fvc::div(phia, Ua) + fvc::div(phi1, U1)
- fvc::div(phia)*Ua; - fvc::div(phi1)*U1;
DDtUb = DDtU2 =
fvc::ddt(Ub) fvc::ddt(U2)
+ fvc::div(phib, Ub) + fvc::div(phi2, U2)
- fvc::div(phib)*Ub; - fvc::div(phi2)*U2;
} }

View File

@ -1,72 +1,74 @@
fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime); fvVectorMatrix U1Eqn(U1, U1.dimensions()*dimVol/dimTime);
fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime); fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
{ {
volTensorField Rca(-nuEffa*(T(fvc::grad(Ua)))); volTensorField Rc1(-nuEff1*(T(fvc::grad(U1))));
Rca = Rca + (2.0/3.0)*sqr(Ct)*I*k - (2.0/3.0)*I*tr(Rca); Rc1 = Rc1 + (2.0/3.0)*sqr(Ct)*I*k - (2.0/3.0)*I*tr(Rc1);
surfaceScalarField phiRa surfaceScalarField phiR1
( (
- fvc::interpolate(nuEffa) - fvc::interpolate(nuEff1)
*mesh.magSf()*fvc::snGrad(alpha)/fvc::interpolate(alpha + scalar(0.001)) *mesh.magSf()*fvc::snGrad(alpha1)
/fvc::interpolate(alpha1 + scalar(0.001))
); );
UaEqn = U1Eqn =
( (
(scalar(1) + Cvm*rhob*beta/rhoa)* (scalar(1) + Cvm*rho2*alpha2/rho1)*
( (
fvm::ddt(Ua) fvm::ddt(U1)
+ fvm::div(phia, Ua, "div(phia,Ua)") + fvm::div(phi1, U1, "div(phi1,U1)")
- fvm::Sp(fvc::div(phia), Ua) - fvm::Sp(fvc::div(phi1), U1)
) )
- fvm::laplacian(nuEffa, Ua) - fvm::laplacian(nuEff1, U1)
+ fvc::div(Rca) + fvc::div(Rc1)
+ fvm::div(phiRa, Ua, "div(phia,Ua)") + fvm::div(phiR1, U1, "div(phi1,U1)")
- fvm::Sp(fvc::div(phiRa), Ua) - fvm::Sp(fvc::div(phiR1), U1)
+ (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca) + (fvc::grad(alpha1)/(fvc::average(alpha1) + scalar(0.001)) & Rc1)
== ==
// g // Buoyancy term transfered to p-equation // g // Buoyancy term transfered to p-equation
- fvm::Sp(beta/rhoa*dragCoef, Ua) - fvm::Sp(alpha2/rho1*dragCoef, U1)
//+ beta/rhoa*dragCoef*Ub // Explicit drag transfered to p-equation //+ alpha2/rho1*dragCoef*U2 // Explicit drag transfered to p-equation
- beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb) - alpha2/rho1*(liftCoeff - Cvm*rho2*DDtU2)
); );
UaEqn.relax(); U1Eqn.relax();
volTensorField Rcb(-nuEffb*T(fvc::grad(Ub))); volTensorField Rc2(-nuEff2*T(fvc::grad(U2)));
Rcb = Rcb + (2.0/3.0)*I*k - (2.0/3.0)*I*tr(Rcb); Rc2 = Rc2 + (2.0/3.0)*I*k - (2.0/3.0)*I*tr(Rc2);
surfaceScalarField phiRb surfaceScalarField phiR2
( (
- fvc::interpolate(nuEffb) - fvc::interpolate(nuEff2)
*mesh.magSf()*fvc::snGrad(beta)/fvc::interpolate(beta + scalar(0.001)) *mesh.magSf()*fvc::snGrad(alpha2)
/fvc::interpolate(alpha2 + scalar(0.001))
); );
UbEqn = U2Eqn =
( (
(scalar(1) + Cvm*rhob*alpha/rhob)* (scalar(1) + Cvm*rho2*alpha1/rho2)*
( (
fvm::ddt(Ub) fvm::ddt(U2)
+ fvm::div(phib, Ub, "div(phib,Ub)") + fvm::div(phi2, U2, "div(phi2,U2)")
- fvm::Sp(fvc::div(phib), Ub) - fvm::Sp(fvc::div(phi2), U2)
) )
- fvm::laplacian(nuEffb, Ub) - fvm::laplacian(nuEff2, U2)
+ fvc::div(Rcb) + fvc::div(Rc2)
+ fvm::div(phiRb, Ub, "div(phib,Ub)") + fvm::div(phiR2, U2, "div(phi2,U2)")
- fvm::Sp(fvc::div(phiRb), Ub) - fvm::Sp(fvc::div(phiR2), U2)
+ (fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb) + (fvc::grad(alpha2)/(fvc::average(alpha2) + scalar(0.001)) & Rc2)
== ==
// g // Buoyancy term transfered to p-equation // g // Buoyancy term transfered to p-equation
- fvm::Sp(alpha/rhob*dragCoef, Ub) - fvm::Sp(alpha1/rho2*dragCoef, U2)
//+ alpha/rhob*dragCoef*Ua // Explicit drag transfered to p-equation //+ alpha1/rho2*dragCoef*U1 // Explicit drag transfered to p-equation
+ alpha/rhob*(liftCoeff + Cvm*rhob*DDtUa) + alpha1/rho2*(liftCoeff + Cvm*rho2*DDtU1)
); );
UbEqn.relax(); U2Eqn.relax();
} }

View File

@ -1,7 +1,7 @@
{ {
word scheme("div(phi,alpha)"); word scheme("div(phi,alpha1)");
surfaceScalarField phir(phia - phib); surfaceScalarField phir(phi1 - phi2);
Info<< "Max Ur Courant Number = " Info<< "Max Ur Courant Number = "
<< ( << (
@ -15,42 +15,47 @@
for (int acorr=0; acorr<nAlphaCorr; acorr++) for (int acorr=0; acorr<nAlphaCorr; acorr++)
{ {
fvScalarMatrix alphaEqn fvScalarMatrix alpha1Eqn
( (
fvm::ddt(alpha) fvm::ddt(alpha1)
+ fvm::div(phi, alpha, scheme) + fvm::div(phi, alpha1, scheme)
+ fvm::div(-fvc::flux(-phir, beta, scheme), alpha, scheme) + fvm::div(-fvc::flux(-phir, alpha2, scheme), alpha1, scheme)
); );
alphaEqn.relax(); alpha1Eqn.relax();
alphaEqn.solve(); alpha1Eqn.solve();
/* /*
fvScalarMatrix betaEqn fvScalarMatrix alpha2Eqn
( (
fvm::ddt(beta) fvm::ddt(alpha2)
+ fvm::div(phi, beta, scheme) + fvm::div(phi, alpha2, scheme)
+ fvm::div(-fvc::flux(phir, scalar(1) - beta, scheme), beta, scheme) + fvm::div
(
-fvc::flux(phir, scalar(1) - alpha2, scheme),
alpha2,
scheme
)
); );
betaEqn.relax(); alpha2Eqn.relax();
betaEqn.solve(); alpha2Eqn.solve();
alpha = alpha1 =
0.5 0.5
*( *(
scalar(1) scalar(1)
+ sqr(scalar(1) - beta) + sqr(scalar(1) - alpha2)
- sqr(scalar(1) - alpha) - sqr(scalar(1) - alpha1)
); );
*/ */
beta = scalar(1) - alpha; alpha2 = scalar(1) - alpha1;
} }
Info<< "Dispersed phase volume fraction = " Info<< "Dispersed phase volume fraction = "
<< alpha.weightedAverage(mesh.V()).value() << alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha) = " << min(alpha).value() << " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha) = " << max(alpha).value() << " Max(alpha1) = " << max(alpha1).value()
<< endl; << endl;
} }
rho = alpha*rhoa + beta*rhob; rho = alpha1*rho1 + alpha2*rho2;

View File

@ -90,7 +90,7 @@ int main(int argc, char *argv[])
if (pimple.turbCorr()) if (pimple.turbCorr())
{ {
#include "kEpsilon.H" #include "kEpsilon.H"
nuEffa = sqr(Ct)*nutb + nua; nuEff1 = sqr(Ct)*nut2 + nu1;
} }
} }

View File

@ -1,9 +1,9 @@
Info<< "Reading field alpha\n" << endl; Info<< "Reading field alpha1\n" << endl;
volScalarField alpha volScalarField alpha1
( (
IOobject IOobject
( (
"alpha", "alpha1",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::MUST_READ, IOobject::MUST_READ,
@ -12,18 +12,18 @@
mesh mesh
); );
volScalarField beta volScalarField alpha2
( (
IOobject IOobject
( (
"beta", "alpha2",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
scalar(1) - alpha scalar(1) - alpha1
//,alpha.boundaryField().types() //,alpha1.boundaryField().types()
); );
Info<< "Reading field p\n" << endl; Info<< "Reading field p\n" << endl;
@ -40,12 +40,12 @@
mesh mesh
); );
Info<< "Reading field Ua\n" << endl; Info<< "Reading field U1\n" << endl;
volVectorField Ua volVectorField U1
( (
IOobject IOobject
( (
"Ua", "U1",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::MUST_READ, IOobject::MUST_READ,
@ -54,12 +54,12 @@
mesh mesh
); );
Info<< "Reading field Ub\n" << endl; Info<< "Reading field U2\n" << endl;
volVectorField Ub volVectorField U2
( (
IOobject IOobject
( (
"Ub", "U2",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::MUST_READ, IOobject::MUST_READ,
@ -78,7 +78,7 @@
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
alpha*Ua + beta*Ub alpha1*U1 + alpha2*U2
); );
@ -96,34 +96,34 @@
) )
); );
dimensionedScalar rhoa dimensionedScalar rho1
( (
transportProperties.lookup("rhoa") transportProperties.lookup("rho1")
); );
dimensionedScalar rhob dimensionedScalar rho2
( (
transportProperties.lookup("rhob") transportProperties.lookup("rho2")
); );
dimensionedScalar nua dimensionedScalar nu1
( (
transportProperties.lookup("nua") transportProperties.lookup("nu1")
); );
dimensionedScalar nub dimensionedScalar nu2
( (
transportProperties.lookup("nub") transportProperties.lookup("nu2")
); );
dimensionedScalar da dimensionedScalar d1
( (
transportProperties.lookup("da") transportProperties.lookup("d1")
); );
dimensionedScalar db dimensionedScalar d2
( (
transportProperties.lookup("db") transportProperties.lookup("d2")
); );
dimensionedScalar Cvm dimensionedScalar Cvm
@ -141,8 +141,8 @@
transportProperties.lookup("Ct") transportProperties.lookup("Ct")
); );
#include "createPhia.H" #include "createPhi1.H"
#include "createPhib.H" #include "createPhi2.H"
surfaceScalarField phi surfaceScalarField phi
( (
@ -152,8 +152,8 @@
runTime.timeName(), runTime.timeName(),
mesh mesh
), ),
fvc::interpolate(alpha)*phia fvc::interpolate(alpha1)*phi1
+ fvc::interpolate(beta)*phib + fvc::interpolate(alpha2)*phi2
); );
volScalarField rho volScalarField rho
@ -164,25 +164,25 @@
runTime.timeName(), runTime.timeName(),
mesh mesh
), ),
alpha*rhoa + beta*rhob alpha1*rho1 + alpha2*rho2
); );
#include "createRASTurbulence.H" #include "createRASTurbulence.H"
Info<< "Calculating field DDtUa and DDtUb\n" << endl; Info<< "Calculating field DDtU1 and DDtU2\n" << endl;
volVectorField DDtUa volVectorField DDtU1
( (
fvc::ddt(Ua) fvc::ddt(U1)
+ fvc::div(phia, Ua) + fvc::div(phi1, U1)
- fvc::div(phia)*Ua - fvc::div(phi1)*U1
); );
volVectorField DDtUb volVectorField DDtU2
( (
fvc::ddt(Ub) fvc::ddt(U2)
+ fvc::div(phib, Ub) + fvc::div(phi2, U2)
- fvc::div(phib)*Ub - fvc::div(phi2)*U2
); );

View File

@ -1,24 +1,24 @@
IOobject phiaHeader IOobject phi1Header
( (
"phia", "phi1",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ IOobject::NO_READ
); );
autoPtr<surfaceScalarField> phiaPtr(NULL); autoPtr<surfaceScalarField> phi1Ptr(NULL);
if (phiaHeader.headerOk()) if (phi1Header.headerOk())
{ {
Info<< "Reading face flux field phia\n" << endl; Info<< "Reading face flux field phi1\n" << endl;
phiaPtr.reset phi1Ptr.reset
( (
new surfaceScalarField new surfaceScalarField
( (
IOobject IOobject
( (
"phia", "phi1",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::MUST_READ, IOobject::MUST_READ,
@ -30,38 +30,38 @@
} }
else else
{ {
Info<< "Calculating face flux field phia\n" << endl; Info<< "Calculating face flux field phi1\n" << endl;
wordList phiTypes wordList phiTypes
( (
Ua.boundaryField().size(), U1.boundaryField().size(),
calculatedFvPatchScalarField::typeName calculatedFvPatchScalarField::typeName
); );
forAll(Ua.boundaryField(), i) forAll(U1.boundaryField(), i)
{ {
if (isA<fixedValueFvPatchVectorField>(Ua.boundaryField()[i])) if (isA<fixedValueFvPatchVectorField>(U1.boundaryField()[i]))
{ {
phiTypes[i] = fixedValueFvPatchScalarField::typeName; phiTypes[i] = fixedValueFvPatchScalarField::typeName;
} }
} }
phiaPtr.reset phi1Ptr.reset
( (
new surfaceScalarField new surfaceScalarField
( (
IOobject IOobject
( (
"phia", "phi1",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
fvc::interpolate(Ua) & mesh.Sf(), fvc::interpolate(U1) & mesh.Sf(),
phiTypes phiTypes
) )
); );
} }
surfaceScalarField& phia = phiaPtr(); surfaceScalarField& phi1 = phi1Ptr();

View File

@ -1,24 +1,24 @@
IOobject phibHeader IOobject phi2Header
( (
"phib", "phi2",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ IOobject::NO_READ
); );
autoPtr<surfaceScalarField> phibPtr(NULL); autoPtr<surfaceScalarField> phi2Ptr(NULL);
if (phibHeader.headerOk()) if (phi2Header.headerOk())
{ {
Info<< "Reading face flux field phib\n" << endl; Info<< "Reading face flux field phi2\n" << endl;
phibPtr.reset phi2Ptr.reset
( (
new surfaceScalarField new surfaceScalarField
( (
IOobject IOobject
( (
"phib", "phi2",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::MUST_READ, IOobject::MUST_READ,
@ -30,38 +30,38 @@
} }
else else
{ {
Info<< "Calculating face flux field phib\n" << endl; Info<< "Calculating face flux field phi2\n" << endl;
wordList phiTypes wordList phiTypes
( (
Ub.boundaryField().size(), U2.boundaryField().size(),
calculatedFvPatchScalarField::typeName calculatedFvPatchScalarField::typeName
); );
forAll(Ub.boundaryField(), i) forAll(U2.boundaryField(), i)
{ {
if (isA<fixedValueFvPatchVectorField>(Ub.boundaryField()[i])) if (isA<fixedValueFvPatchVectorField>(U2.boundaryField()[i]))
{ {
phiTypes[i] = fixedValueFvPatchScalarField::typeName; phiTypes[i] = fixedValueFvPatchScalarField::typeName;
} }
} }
phibPtr.reset phi2Ptr.reset
( (
new surfaceScalarField new surfaceScalarField
( (
IOobject IOobject
( (
"phib", "phi2",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
fvc::interpolate(Ub) & mesh.Sf(), fvc::interpolate(U2) & mesh.Sf(),
phiTypes phiTypes
) )
); );
} }
surfaceScalarField& phib = phibPtr(); surfaceScalarField& phi2 = phi2Ptr();

View File

@ -51,21 +51,21 @@
) )
); );
dimensionedScalar alphak dimensionedScalar alpha1k
( (
dimensionedScalar::lookupOrAddToDict dimensionedScalar::lookupOrAddToDict
( (
"alphak", "alpha1k",
kEpsilonDict, kEpsilonDict,
1.0 1.0
) )
); );
dimensionedScalar alphaEps dimensionedScalar alpha1Eps
( (
dimensionedScalar::lookupOrAddToDict dimensionedScalar::lookupOrAddToDict
( (
"alphaEps", "alpha1Eps",
kEpsilonDict, kEpsilonDict,
0.76923 0.76923
) )
@ -135,12 +135,12 @@
); );
Info<< "Calculating field nutb\n" << endl; Info<< "Calculating field nut2\n" << endl;
volScalarField nutb volScalarField nut2
( (
IOobject IOobject
( (
"nutb", "nut2",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
@ -149,30 +149,30 @@
Cmu*sqr(k)/epsilon Cmu*sqr(k)/epsilon
); );
Info<< "Calculating field nuEffa\n" << endl; Info<< "Calculating field nuEff1\n" << endl;
volScalarField nuEffa volScalarField nuEff1
( (
IOobject IOobject
( (
"nuEffa", "nuEff1",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
sqr(Ct)*nutb + nua sqr(Ct)*nut2 + nu1
); );
Info<< "Calculating field nuEffb\n" << endl; Info<< "Calculating field nuEff2\n" << endl;
volScalarField nuEffb volScalarField nuEff2
( (
IOobject IOobject
( (
"nuEffb", "nuEff2",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
nutb + nub nut2 + nu2
); );

View File

@ -5,9 +5,9 @@ if (turbulence)
y.correct(); y.correct();
} }
tmp<volTensorField> tgradUb = fvc::grad(Ub); tmp<volTensorField> tgradU2 = fvc::grad(U2);
volScalarField G(2*nutb*(tgradUb() && dev(symm(tgradUb())))); volScalarField G(2*nut2*(tgradU2() && dev(symm(tgradU2()))));
tgradUb.clear(); tgradU2.clear();
#include "wallFunctions.H" #include "wallFunctions.H"
@ -15,11 +15,11 @@ if (turbulence)
fvScalarMatrix epsEqn fvScalarMatrix epsEqn
( (
fvm::ddt(epsilon) fvm::ddt(epsilon)
+ fvm::div(phib, epsilon) + fvm::div(phi2, epsilon)
- fvm::Sp(fvc::div(phib), epsilon) - fvm::Sp(fvc::div(phi2), epsilon)
- fvm::laplacian - fvm::laplacian
( (
alphaEps*nuEffb, epsilon, alpha1Eps*nuEff2, epsilon,
"laplacian(DepsilonEff,epsilon)" "laplacian(DepsilonEff,epsilon)"
) )
== ==
@ -39,11 +39,11 @@ if (turbulence)
fvScalarMatrix kEqn fvScalarMatrix kEqn
( (
fvm::ddt(k) fvm::ddt(k)
+ fvm::div(phib, k) + fvm::div(phi2, k)
- fvm::Sp(fvc::div(phib), k) - fvm::Sp(fvc::div(phi2), k)
- fvm::laplacian - fvm::laplacian
( (
alphak*nuEffb, k, alpha1k*nuEff2, k,
"laplacian(DkEff,k)" "laplacian(DkEff,k)"
) )
== ==
@ -56,9 +56,9 @@ if (turbulence)
k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8)); k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8));
//- Re-calculate turbulence viscosity //- Re-calculate turbulence viscosity
nutb = Cmu*sqr(k)/epsilon; nut2 = Cmu*sqr(k)/epsilon;
#include "wallViscosity.H" #include "wallViscosity.H"
} }
nuEffb = nutb + nub; nuEff2 = nut2 + nu2;

View File

@ -1,23 +1,23 @@
volVectorField Ur(Ua - Ub); volVectorField Ur(U1 - U2);
volScalarField magUr(mag(Ur)); volScalarField magUr(mag(Ur));
volScalarField CdaMagUr volScalarField Cd1MagUr
( (
(24.0*nub/da)*(scalar(1) + 0.15*pow(da*magUr/nub, 0.687)) (24.0*nu2/d1)*(scalar(1) + 0.15*pow(d1*magUr/nu2, 0.687))
); );
volScalarField CdbMagUr volScalarField Cd2MagUr
( (
(24.0*nua/db)*(scalar(1) + 0.15*pow(db*magUr/nua, 0.687)) (24.0*nu1/d2)*(scalar(1) + 0.15*pow(d2*magUr/nu1, 0.687))
); );
volScalarField dragCoef volScalarField dragCoef
( (
"Cd", "Cd",
0.75*(beta*rhob*CdaMagUr/da + alpha*rhoa*CdbMagUr/db) 0.75*(alpha2*rho2*Cd1MagUr/d1 + alpha1*rho1*Cd2MagUr/d2)
); );
volVectorField liftCoeff volVectorField liftCoeff
( (
Cl*(beta*rhob + alpha*rhoa)*(Ur ^ fvc::curl(U)) Cl*(alpha2*rho2 + alpha1*rho1)*(Ur ^ fvc::curl(U))
); );

View File

@ -1,45 +1,45 @@
{ {
surfaceScalarField alphaf(fvc::interpolate(alpha)); surfaceScalarField alpha1f(fvc::interpolate(alpha1));
surfaceScalarField betaf(scalar(1) - alphaf); surfaceScalarField alpha2f(scalar(1) - alpha1f);
volScalarField rUaA(1.0/UaEqn.A()); volScalarField rAU1(1.0/U1Eqn.A());
volScalarField rUbA(1.0/UbEqn.A()); volScalarField rAU2(1.0/U2Eqn.A());
surfaceScalarField rUaAf(fvc::interpolate(rUaA)); surfaceScalarField rAU1f(fvc::interpolate(rAU1));
surfaceScalarField rUbAf(fvc::interpolate(rUbA)); surfaceScalarField rAU2f(fvc::interpolate(rAU2));
Ua = rUaA*UaEqn.H(); U1 = rAU1*U1Eqn.H();
Ub = rUbA*UbEqn.H(); U2 = rAU2*U2Eqn.H();
surfaceScalarField phiDraga surfaceScalarField phiDrag1
( (
fvc::interpolate(beta/rhoa*dragCoef*rUaA)*phib + rUaAf*(g & mesh.Sf()) fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2 + rAU1f*(g & mesh.Sf())
); );
surfaceScalarField phiDragb surfaceScalarField phiDrag2
( (
fvc::interpolate(alpha/rhob*dragCoef*rUbA)*phia + rUbAf*(g & mesh.Sf()) fvc::interpolate(alpha1/rho2*dragCoef*rAU2)*phi1 + rAU2f*(g & mesh.Sf())
); );
forAll(p.boundaryField(), patchi) forAll(p.boundaryField(), patchi)
{ {
if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi])) if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
{ {
phiDraga.boundaryField()[patchi] = 0.0; phiDrag1.boundaryField()[patchi] = 0.0;
phiDragb.boundaryField()[patchi] = 0.0; phiDrag2.boundaryField()[patchi] = 0.0;
} }
} }
phia = (fvc::interpolate(Ua) & mesh.Sf()) + fvc::ddtPhiCorr(rUaA, Ua, phia) phi1 = (fvc::interpolate(U1) & mesh.Sf()) + fvc::ddtPhiCorr(rAU1, U1, phi1)
+ phiDraga; + phiDrag1;
phib = (fvc::interpolate(Ub) & mesh.Sf()) + fvc::ddtPhiCorr(rUbA, Ub, phib) phi2 = (fvc::interpolate(U2) & mesh.Sf()) + fvc::ddtPhiCorr(rAU2, U2, phi2)
+ phiDragb; + phiDrag2;
phi = alphaf*phia + betaf*phib; phi = alpha1f*phi1 + alpha2f*phi2;
surfaceScalarField Dp surfaceScalarField Dp
( (
"(rho*(1|A(U)))", "(rho*(1|A(U)))",
alphaf*rUaAf/rhoa + betaf*rUbAf/rhob alpha1f*rAU1f/rho1 + alpha2f*rAU2f/rho2
); );
while (pimple.correctNonOrthogonal()) while (pimple.correctNonOrthogonal())
@ -57,22 +57,22 @@
{ {
surfaceScalarField SfGradp(pEqn.flux()/Dp); surfaceScalarField SfGradp(pEqn.flux()/Dp);
phia -= rUaAf*SfGradp/rhoa; phi1 -= rAU1f*SfGradp/rho1;
phib -= rUbAf*SfGradp/rhob; phi2 -= rAU2f*SfGradp/rho2;
phi = alphaf*phia + betaf*phib; phi = alpha1f*phi1 + alpha2f*phi2;
p.relax(); p.relax();
SfGradp = pEqn.flux()/Dp; SfGradp = pEqn.flux()/Dp;
Ua += (fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa)); U1 += (fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1));
//Ua += rUaA*(fvc::reconstruct(phiDraga/rUaAf - SfGradp/rhoa)); //U1 += rAU1*(fvc::reconstruct(phiDrag1/rAU1f - SfGradp/rho1));
Ua.correctBoundaryConditions(); U1.correctBoundaryConditions();
Ub += (fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob)); U2 += (fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2));
//Ub += rUbA*(fvc::reconstruct(phiDragb/rUbAf - SfGradp/rhob)); //U2 += rAU2*(fvc::reconstruct(phiDrag2/rAU2f - SfGradp/rho2));
Ub.correctBoundaryConditions(); U2.correctBoundaryConditions();
U = alpha*Ua + beta*Ub; U = alpha1*U1 + alpha2*U2;
} }
} }
} }

View File

@ -4,7 +4,7 @@
scalar Cmu25 = ::pow(Cmu.value(), 0.25); scalar Cmu25 = ::pow(Cmu.value(), 0.25);
scalar Cmu75 = ::pow(Cmu.value(), 0.75); scalar Cmu75 = ::pow(Cmu.value(), 0.75);
scalar kappa_ = kappa.value(); scalar kappa_ = kappa.value();
scalar nub_ = nub.value(); scalar nu2_ = nu2.value();
const fvPatchList& patches = mesh.boundary(); const fvPatchList& patches = mesh.boundary();
@ -33,9 +33,9 @@
if (isA<wallFvPatch>(currPatch)) if (isA<wallFvPatch>(currPatch))
{ {
const scalarField& nutbw = nutb.boundaryField()[patchi]; const scalarField& nut2w = nut2.boundaryField()[patchi];
scalarField magFaceGradU(mag(Ub.boundaryField()[patchi].snGrad())); scalarField magFaceGradU(mag(U2.boundaryField()[patchi].snGrad()));
forAll(currPatch, facei) forAll(currPatch, facei)
{ {
@ -52,7 +52,7 @@
/(kappa_*y[patchi][facei]); /(kappa_*y[patchi][facei]);
G[faceCelli] += G[faceCelli] +=
(nutbw[facei] + nub_)*magFaceGradU[facei] (nut2w[facei] + nu2_)*magFaceGradU[facei]
*Cmu25*::sqrt(k[faceCelli]) *Cmu25*::sqrt(k[faceCelli])
/(kappa_*y[patchi][facei]); /(kappa_*y[patchi][facei]);
} }

View File

@ -2,7 +2,7 @@
scalar Cmu25 = ::pow(Cmu.value(), 0.25); scalar Cmu25 = ::pow(Cmu.value(), 0.25);
scalar kappa_ = kappa.value(); scalar kappa_ = kappa.value();
scalar E_ = E.value(); scalar E_ = E.value();
scalar nub_ = nub.value(); scalar nu2_ = nu2.value();
const fvPatchList& patches = mesh.boundary(); const fvPatchList& patches = mesh.boundary();
@ -12,7 +12,7 @@
if (isA<wallFvPatch>(currPatch)) if (isA<wallFvPatch>(currPatch))
{ {
scalarField& nutw = nutb.boundaryField()[patchi]; scalarField& nutw = nut2.boundaryField()[patchi];
forAll(currPatch, facei) forAll(currPatch, facei)
{ {
@ -20,11 +20,11 @@
// calculate yPlus // calculate yPlus
scalar yPlus = scalar yPlus =
Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nub_; Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nu2_;
if (yPlus > 11.6) if (yPlus > 11.6)
{ {
nutw[facei] = nub_*(yPlus*kappa_/::log(E_*yPlus) -1); nutw[facei] = nu2_*(yPlus*kappa_/::log(E_*yPlus) -1);
} }
else else
{ {

View File

@ -10,7 +10,7 @@
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
Ua - Ub U1 - U2
); );
runTime.write(); runTime.write();

View File

@ -5,8 +5,8 @@
fvScalarMatrix T1Eqn fvScalarMatrix T1Eqn
( (
fvm::ddt(alpha1, T1) fvm::ddt(alpha1, T1)
+ fvm::div(alphaPhi1, T1) + fvm::div(alpha1Phi1, T1)
- fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), T1) - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alpha1Phi1), T1)
- fvm::laplacian(kByCp1, T1) - fvm::laplacian(kByCp1, T1)
== ==
heatTransferCoeff*T2/Cp1/rho1 heatTransferCoeff*T2/Cp1/rho1
@ -17,8 +17,8 @@
fvScalarMatrix T2Eqn fvScalarMatrix T2Eqn
( (
fvm::ddt(alpha2, T2) fvm::ddt(alpha2, T2)
+ fvm::div(alphaPhi2, T2) + fvm::div(alpha1Phi2, T2)
- fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), T2) - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alpha1Phi2), T2)
- fvm::laplacian(kByCp2, T2) - fvm::laplacian(kByCp2, T2)
== ==
heatTransferCoeff*T1/Cp2/rho2 heatTransferCoeff*T1/Cp2/rho2

View File

@ -8,7 +8,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
if (kineticTheory.on()) if (kineticTheory.on())
{ {
kineticTheory.solve(gradU1T); kineticTheory.solve(gradU1T);
nuEff1 = kineticTheory.mua()/rho1; nuEff1 = kineticTheory.mu1()/rho1;
} }
else // If not using kinetic theory is using Ct model else // If not using kinetic theory is using Ct model
{ {
@ -31,8 +31,8 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
(scalar(1) + Cvm*rho2*alpha2/rho1)* (scalar(1) + Cvm*rho2*alpha2/rho1)*
( (
fvm::ddt(alpha1, U1) fvm::ddt(alpha1, U1)
+ fvm::div(alphaPhi1, U1) + fvm::div(alpha1Phi1, U1)
- fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1) - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alpha1Phi1), U1)
) )
- fvm::laplacian(alpha1*nuEff1, U1) - fvm::laplacian(alpha1*nuEff1, U1)
+ fvc::div(alpha1*Rc1) + fvc::div(alpha1*Rc1)
@ -57,8 +57,8 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
(scalar(1) + Cvm*rho2*alpha1/rho2)* (scalar(1) + Cvm*rho2*alpha1/rho2)*
( (
fvm::ddt(alpha2, U2) fvm::ddt(alpha2, U2)
+ fvm::div(alphaPhi2, U2) + fvm::div(alpha1Phi2, U2)
- fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2) - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alpha1Phi2), U2)
) )
- fvm::laplacian(alpha2*nuEff2, U2) - fvm::laplacian(alpha2*nuEff2, U2)
+ fvc::div(alpha2*Rc2) + fvc::div(alpha2*Rc2)

View File

@ -1,5 +1,5 @@
surfaceScalarField alphaPhi1("alphaPhi1", phi1); surfaceScalarField alpha1Phi1("alpha1Phi1", phi1);
surfaceScalarField alphaPhi2("alphaPhi2", phi2); surfaceScalarField alpha1Phi2("alpha1Phi2", phi2);
{ {
word scheme("div(phi,alpha1)"); word scheme("div(phi,alpha1)");
@ -68,17 +68,17 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
if (g0.value() > 0.0) if (g0.value() > 0.0)
{ {
ppMagf = rU1Af*fvc::interpolate ppMagf = rAU1f*fvc::interpolate
( (
(1.0/(rho1*(alpha1 + scalar(0.0001)))) (1.0/(rho1*(alpha1 + scalar(0.0001))))
*g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax) *g0*min(exp(preAlphaExp*(alpha1 - alpha1Max)), expMax)
); );
alpha1Eqn -= fvm::laplacian alpha1Eqn -= fvm::laplacian
( (
(fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf, (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
alpha1, alpha1,
"laplacian(alphaPpMag,alpha1)" "laplacian(alpha1PpMag,alpha1)"
); );
} }
@ -90,8 +90,8 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
#include "packingLimiter.H" #include "packingLimiter.H"
alphaPhi1 = alpha1Eqn.flux(); alpha1Phi1 = alpha1Eqn.flux();
alphaPhi2 = phi - alphaPhi1; alpha1Phi2 = phi - alpha1Phi1;
alpha2 = scalar(1) - alpha1; alpha2 = scalar(1) - alpha1;
Info<< "Dispersed phase volume fraction = " Info<< "Dispersed phase volume fraction = "

View File

@ -324,11 +324,11 @@
drag1 drag1
); );
surfaceScalarField rU1Af surfaceScalarField rAU1f
( (
IOobject IOobject
( (
"rU1Af", "rAU1f",
runTime.timeName(), runTime.timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,

View File

@ -52,21 +52,21 @@
) )
); );
dimensionedScalar alphak dimensionedScalar alpha1k
( (
dimensionedScalar::lookupOrAddToDict dimensionedScalar::lookupOrAddToDict
( (
"alphak", "alpha1k",
kEpsilonDict, kEpsilonDict,
1.0 1.0
) )
); );
dimensionedScalar alphaEps dimensionedScalar alpha1Eps
( (
dimensionedScalar::lookupOrAddToDict dimensionedScalar::lookupOrAddToDict
( (
"alphaEps", "alpha1Eps",
kEpsilonDict, kEpsilonDict,
0.76923 0.76923
) )

View File

@ -69,12 +69,12 @@ volScalarField heatTransferCoeff
<< exit(FatalError); << exit(FatalError);
} }
volScalarField alphaCoeff volScalarField alpha1Coeff
( (
(alpha1 + minInterfaceAlpha)*(alpha2 + minInterfaceAlpha) (alpha1 + minInterfaceAlpha)*(alpha2 + minInterfaceAlpha)
); );
dragCoeff *= alphaCoeff; dragCoeff *= alpha1Coeff;
heatTransferCoeff *= alphaCoeff; heatTransferCoeff *= alpha1Coeff;
liftForce = Cl*(alpha1*rho1 + alpha2*rho2)*(Ur ^ fvc::curl(U)); liftForce = Cl*(alpha1*rho1 + alpha2*rho2)*(Ur ^ fvc::curl(U));
} }

View File

@ -49,12 +49,12 @@ namespace dragModels
Foam::dragModels::Ergun::Ergun Foam::dragModels::Ergun::Ergun
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
dragModel(interfaceDict, alpha, phase1, phase2) dragModel(interfaceDict, alpha1, phase1, phase2)
{} {}
@ -71,12 +71,12 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::Ergun::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6))); volScalarField alpha2(max(scalar(1) - alpha1_, scalar(1.0e-6)));
return return
150.0*alpha_*phase2_.nu()*phase2_.rho() 150.0*alpha1_*phase2_.nu()*phase2_.rho()
/sqr(beta*phase1_.d()) /sqr(alpha2*phase1_.d())
+ 1.75*phase2_.rho()*Ur/(beta*phase1_.d()); + 1.75*phase2_.rho()*Ur/(alpha2*phase1_.d());
} }

View File

@ -68,7 +68,7 @@ public:
Ergun Ergun
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );

View File

@ -49,12 +49,12 @@ namespace dragModels
Foam::dragModels::Gibilaro::Gibilaro Foam::dragModels::Gibilaro::Gibilaro
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
dragModel(interfaceDict, alpha, phase1, phase2) dragModel(interfaceDict, alpha1, phase1, phase2)
{} {}
@ -71,9 +71,9 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::Gibilaro::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6))); volScalarField alpha2(max(scalar(1) - alpha1_, scalar(1.0e-6)));
volScalarField bp(pow(beta, -2.8)); volScalarField bp(pow(alpha2, -2.8));
volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(alpha2*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
return (17.3/Re + scalar(0.336))*phase2_.rho()*Ur*bp/phase1_.d(); return (17.3/Re + scalar(0.336))*phase2_.rho()*Ur*bp/phase1_.d();
} }

View File

@ -68,7 +68,7 @@ public:
Gibilaro Gibilaro
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );

View File

@ -49,12 +49,12 @@ namespace dragModels
Foam::dragModels::GidaspowErgunWenYu::GidaspowErgunWenYu Foam::dragModels::GidaspowErgunWenYu::GidaspowErgunWenYu
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
dragModel(interfaceDict, alpha, phase1, phase2) dragModel(interfaceDict, alpha1, phase1, phase2)
{} {}
@ -71,9 +71,9 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6))); volScalarField alpha2(max(scalar(1) - alpha1_, scalar(1.0e-6)));
volScalarField d(phase1_.d()); volScalarField d(phase1_.d());
volScalarField bp(pow(beta, -2.65)); volScalarField bp(pow(alpha2, -2.65));
volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds volScalarField Cds
@ -85,12 +85,12 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
// Wen and Yu (1966) // Wen and Yu (1966)
return return
( (
pos(beta - 0.8) pos(alpha2 - 0.8)
*(0.75*Cds*phase2_.rho()*Ur*bp/d) *(0.75*Cds*phase2_.rho()*Ur*bp/d)
+ neg(beta - 0.8) + neg(alpha2 - 0.8)
*( *(
150.0*alpha_*phase2_.nu()*phase2_.rho()/(sqr(beta*d)) 150.0*alpha1_*phase2_.nu()*phase2_.rho()/(sqr(alpha2*d))
+ 1.75*phase2_.rho()*Ur/(beta*d) + 1.75*phase2_.rho()*Ur/(alpha2*d)
) )
); );
} }

View File

@ -66,7 +66,7 @@ public:
GidaspowErgunWenYu GidaspowErgunWenYu
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );

View File

@ -49,12 +49,12 @@ namespace dragModels
Foam::dragModels::GidaspowSchillerNaumann::GidaspowSchillerNaumann Foam::dragModels::GidaspowSchillerNaumann::GidaspowSchillerNaumann
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
dragModel(interfaceDict, alpha, phase1, phase2) dragModel(interfaceDict, alpha1, phase1, phase2)
{} {}
@ -71,10 +71,10 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(scalar(1) - alpha_, scalar(1e-6))); volScalarField alpha2(max(scalar(1) - alpha1_, scalar(1e-6)));
volScalarField bp(pow(beta, -2.65)); volScalarField bp(pow(alpha2, -2.65));
volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(alpha2*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds volScalarField Cds
( (
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)

View File

@ -75,7 +75,7 @@ public:
GidaspowSchillerNaumann GidaspowSchillerNaumann
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );

View File

@ -49,12 +49,12 @@ namespace dragModels
Foam::dragModels::SchillerNaumann::SchillerNaumann Foam::dragModels::SchillerNaumann::SchillerNaumann
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
dragModel(interfaceDict, alpha, phase1, phase2) dragModel(interfaceDict, alpha1, phase1, phase2)
{} {}

View File

@ -64,7 +64,7 @@ public:
SchillerNaumann SchillerNaumann
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );

View File

@ -49,12 +49,12 @@ namespace dragModels
Foam::dragModels::SyamlalOBrien::SyamlalOBrien Foam::dragModels::SyamlalOBrien::SyamlalOBrien
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
dragModel(interfaceDict, alpha, phase1, phase2) dragModel(interfaceDict, alpha1, phase1, phase2)
{} {}
@ -71,12 +71,12 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6))); volScalarField alpha2(max(scalar(1) - alpha1_, scalar(1.0e-6)));
volScalarField A(pow(beta, 4.14)); volScalarField A(pow(alpha2, 4.14));
volScalarField B volScalarField B
( (
neg(beta - 0.85)*(0.8*pow(beta, 1.28)) neg(alpha2 - 0.85)*(0.8*pow(alpha2, 1.28))
+ pos(beta - 0.85)*(pow(beta, 2.65)) + pos(alpha2 - 0.85)*(pow(alpha2, 2.65))
); );
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));

View File

@ -67,7 +67,7 @@ public:
SyamlalOBrien SyamlalOBrien
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );

View File

@ -49,12 +49,12 @@ namespace dragModels
Foam::dragModels::WenYu::WenYu Foam::dragModels::WenYu::WenYu
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
dragModel(interfaceDict, alpha, phase1, phase2) dragModel(interfaceDict, alpha1, phase1, phase2)
{} {}
@ -71,8 +71,8 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6))); volScalarField alpha2(max(scalar(1) - alpha1_, scalar(1.0e-6)));
volScalarField bp(pow(beta, -2.65)); volScalarField bp(pow(alpha2, -2.65));
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds volScalarField Cds

View File

@ -78,7 +78,7 @@ public:
WenYu WenYu
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );

View File

@ -39,13 +39,13 @@ namespace Foam
Foam::dragModel::dragModel Foam::dragModel::dragModel
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
interfaceDict_(interfaceDict), interfaceDict_(interfaceDict),
alpha_(alpha), alpha1_(alpha1),
phase1_(phase1), phase1_(phase1),
phase2_(phase2) phase2_(phase2)
{} {}

View File

@ -55,7 +55,7 @@ protected:
// Protected data // Protected data
const dictionary& interfaceDict_; const dictionary& interfaceDict_;
const volScalarField& alpha_; const volScalarField& alpha1_;
const phaseModel& phase1_; const phaseModel& phase1_;
const phaseModel& phase2_; const phaseModel& phase2_;
@ -75,11 +75,11 @@ public:
dictionary, dictionary,
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
), ),
(interfaceDict, alpha, phase1, phase2) (interfaceDict, alpha1, phase1, phase2)
); );
@ -88,7 +88,7 @@ public:
dragModel dragModel
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );
@ -103,7 +103,7 @@ public:
static autoPtr<dragModel> New static autoPtr<dragModel> New
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );
@ -112,12 +112,13 @@ public:
// Member Functions // Member Functions
//- the dragfunction K used in the momentum eq. //- the dragfunction K used in the momentum eq.
// ddt(alpha*rhoa*Ua) + ... = ... alpha*beta*K*(Ua-Ub) // ddt(alpha1*rho1*U1) + ... = ... alpha1*alpha2*K*(U1-U2)
// ddt(beta*rhob*Ub) + ... = ... alpha*beta*K*(Ub-Ua) // ddt(alpha2*rho2*U2) + ... = ... alpha1*alpha2*K*(U2-U1)
// ********************************** NB! ***************************** // ********************************** NB! *****************************
// for numerical reasons alpha and beta has been // for numerical reasons alpha1 and alpha2 has been
// extracted from the dragFunction K, // extracted from the dragFunction K,
// so you MUST divide K by alpha*beta when implemnting the drag function // so you MUST divide K by alpha1*alpha2 when implemnting the drag
// function
// ********************************** NB! ***************************** // ********************************** NB! *****************************
virtual tmp<volScalarField> K(const volScalarField& Ur) const = 0; virtual tmp<volScalarField> K(const volScalarField& Ur) const = 0;
}; };

View File

@ -30,7 +30,7 @@ License
Foam::autoPtr<Foam::dragModel> Foam::dragModel::New Foam::autoPtr<Foam::dragModel> Foam::dragModel::New
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
@ -58,7 +58,7 @@ Foam::autoPtr<Foam::dragModel> Foam::dragModel::New
<< exit(FatalError); << exit(FatalError);
} }
return cstrIter()(interfaceDict, alpha, phase1, phase2); return cstrIter()(interfaceDict, alpha1, phase1, phase2);
} }

View File

@ -49,12 +49,12 @@ namespace heatTransferModels
Foam::heatTransferModels::RanzMarshall::RanzMarshall Foam::heatTransferModels::RanzMarshall::RanzMarshall
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
heatTransferModel(interfaceDict, alpha, phase1, phase2) heatTransferModel(interfaceDict, alpha1, phase1, phase2)
{} {}

View File

@ -64,7 +64,7 @@ public:
RanzMarshall RanzMarshall
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );

View File

@ -39,13 +39,13 @@ namespace Foam
Foam::heatTransferModel::heatTransferModel Foam::heatTransferModel::heatTransferModel
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
interfaceDict_(interfaceDict), interfaceDict_(interfaceDict),
alpha_(alpha), alpha1_(alpha1),
phase1_(phase1), phase1_(phase1),
phase2_(phase2) phase2_(phase2)
{} {}

View File

@ -55,7 +55,7 @@ protected:
// Protected data // Protected data
const dictionary& interfaceDict_; const dictionary& interfaceDict_;
const volScalarField& alpha_; const volScalarField& alpha1_;
const phaseModel& phase1_; const phaseModel& phase1_;
const phaseModel& phase2_; const phaseModel& phase2_;
@ -75,11 +75,11 @@ public:
dictionary, dictionary,
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
), ),
(interfaceDict, alpha, phase1, phase2) (interfaceDict, alpha1, phase1, phase2)
); );
@ -88,7 +88,7 @@ public:
heatTransferModel heatTransferModel
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );
@ -103,7 +103,7 @@ public:
static autoPtr<heatTransferModel> New static autoPtr<heatTransferModel> New
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );
@ -112,12 +112,12 @@ public:
// Member Functions // Member Functions
//- the heat-transfer function K used in the enthalpy eq. //- the heat-transfer function K used in the enthalpy eq.
// ddt(alpha*rhoa*ha) + ... = ... alpha*beta*K*(Ta - Tb) // ddt(alpha1*rho1*ha) + ... = ... alpha1*alpha2*K*(Ta - Tb)
// ddt(beta*rhob*hb) + ... = ... alpha*beta*K*(Tb - Ta) // ddt(alpha2*rho2*hb) + ... = ... alpha1*alpha2*K*(Tb - Ta)
// ********************************** NB!***************************** // ********************************** NB!*****************************
// for numerical reasons alpha and beta has been // for numerical reasons alpha1 and alpha2 has been
// extracted from the heat-transfer function K, // extracted from the heat-transfer function K,
// so you MUST divide K by alpha*beta when implementing the // so you MUST divide K by alpha1*alpha2 when implementing the
// heat-transfer function // heat-transfer function
// ********************************** NB!***************************** // ********************************** NB!*****************************
virtual tmp<volScalarField> K(const volScalarField& Ur) const = 0; virtual tmp<volScalarField> K(const volScalarField& Ur) const = 0;

View File

@ -30,7 +30,7 @@ License
Foam::autoPtr<Foam::heatTransferModel> Foam::heatTransferModel::New Foam::autoPtr<Foam::heatTransferModel> Foam::heatTransferModel::New
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
@ -58,7 +58,7 @@ Foam::autoPtr<Foam::heatTransferModel> Foam::heatTransferModel::New
<< exit(FatalError); << exit(FatalError);
} }
return cstrIter()(interfaceDict, alpha, phase1, phase2); return cstrIter()(interfaceDict, alpha1, phase1, phase2);
} }

View File

@ -70,21 +70,21 @@ Foam::kineticTheoryModels::conductivityModels::Gidaspow::~Gidaspow()
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::conductivityModels::Gidaspow::kappa Foam::kineticTheoryModels::conductivityModels::Gidaspow::kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
{ {
const scalar sqrtPi = sqrt(constant::mathematical::pi); const scalar sqrtPi = sqrt(constant::mathematical::pi);
return rhoa*da*sqrt(Theta)* return rho1*da*sqrt(Theta)*
( (
2.0*sqr(alpha)*g0*(1.0 + e)/sqrtPi 2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
+ (9.0/8.0)*sqrtPi*g0*0.5*(1.0 + e)*sqr(alpha) + (9.0/8.0)*sqrtPi*g0*0.5*(1.0 + e)*sqr(alpha1)
+ (15.0/16.0)*sqrtPi*alpha + (15.0/16.0)*sqrtPi*alpha1
+ (25.0/64.0)*sqrtPi/((1.0 + e)*g0) + (25.0/64.0)*sqrtPi/((1.0 + e)*g0)
); );
} }

View File

@ -74,10 +74,10 @@ public:
tmp<volScalarField> kappa tmp<volScalarField> kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;

View File

@ -73,10 +73,10 @@ Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
@ -85,15 +85,15 @@ Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa
volScalarField lamda volScalarField lamda
( (
scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_ scalar(1) + da/(6.0*sqrt(2.0)*(alpha1 + scalar(1.0e-5)))/L_
); );
return rhoa*da*sqrt(Theta)* return rho1*da*sqrt(Theta)*
( (
2.0*sqr(alpha)*g0*(1.0 + e)/sqrtPi 2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
+ (9.0/8.0)*sqrtPi*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha) + (9.0/8.0)*sqrtPi*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha1)
/(49.0/16.0 - 33.0*e/16.0) /(49.0/16.0 - 33.0*e/16.0)
+ (15.0/16.0)*sqrtPi*alpha*(0.5*sqr(e) + 0.25*e - 0.75 + lamda) + (15.0/16.0)*sqrtPi*alpha1*(0.5*sqr(e) + 0.25*e - 0.75 + lamda)
/((49.0/16.0 - 33.0*e/16.0)*lamda) /((49.0/16.0 - 33.0*e/16.0)*lamda)
+ (25.0/64.0)*sqrtPi + (25.0/64.0)*sqrtPi
/((1.0 + e)*(49.0/16.0 - 33.0*e/16.0)*lamda*g0) /((1.0 + e)*(49.0/16.0 - 33.0*e/16.0)*lamda*g0)

View File

@ -79,10 +79,10 @@ public:
tmp<volScalarField> kappa tmp<volScalarField> kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;

View File

@ -70,22 +70,22 @@ Foam::kineticTheoryModels::conductivityModels::Syamlal::~Syamlal()
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::conductivityModels::Syamlal::kappa Foam::kineticTheoryModels::conductivityModels::Syamlal::kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
{ {
const scalar sqrtPi = sqrt(constant::mathematical::pi); const scalar sqrtPi = sqrt(constant::mathematical::pi);
return rhoa*da*sqrt(Theta)* return rho1*da*sqrt(Theta)*
( (
2.0*sqr(alpha)*g0*(1.0 + e)/sqrtPi 2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
+ (9.0/8.0)*sqrtPi*g0*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha) + (9.0/8.0)*sqrtPi*g0*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha1)
/(49.0/16.0 - 33.0*e/16.0) /(49.0/16.0 - 33.0*e/16.0)
+ (15.0/32.0)*sqrtPi*alpha/(49.0/16.0 - 33.0*e/16.0) + (15.0/32.0)*sqrtPi*alpha1/(49.0/16.0 - 33.0*e/16.0)
); );
} }

View File

@ -74,10 +74,10 @@ public:
tmp<volScalarField> kappa tmp<volScalarField> kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;

View File

@ -106,10 +106,10 @@ public:
virtual tmp<volScalarField> kappa virtual tmp<volScalarField> kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const = 0; ) const = 0;

View File

@ -72,9 +72,9 @@ Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson:: Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::
frictionalPressure frictionalPressure
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& eta, const dimensionedScalar& eta,
const dimensionedScalar& p const dimensionedScalar& p
@ -82,8 +82,8 @@ frictionalPressure
{ {
return return
Fr*pow(max(alpha - alphaMinFriction, scalar(0)), eta) Fr*pow(max(alpha1 - alpha1MinFriction, scalar(0)), eta)
/pow(max(alphaMax - alpha, scalar(5.0e-2)), p); /pow(max(alpha1Max - alpha1, scalar(5.0e-2)), p);
} }
@ -91,9 +91,9 @@ Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson:: Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::
frictionalPressurePrime frictionalPressurePrime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& eta, const dimensionedScalar& eta,
const dimensionedScalar& p const dimensionedScalar& p
@ -101,17 +101,18 @@ frictionalPressurePrime
{ {
return Fr* return Fr*
( (
eta*pow(max(alpha - alphaMinFriction, scalar(0)), eta - 1.0) eta*pow(max(alpha1 - alpha1MinFriction, scalar(0)), eta - 1.0)
*(alphaMax-alpha) + p*pow(max(alpha - alphaMinFriction, scalar(0)), eta) *(alpha1Max-alpha1)
)/pow(max(alphaMax - alpha, scalar(5.0e-2)), p + 1.0); + p*pow(max(alpha1 - alpha1MinFriction, scalar(0)), eta)
)/pow(max(alpha1Max - alpha1, scalar(5.0e-2)), p + 1.0);
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::muf Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::muf
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const volScalarField& pf, const volScalarField& pf,
const volSymmTensorField& D, const volSymmTensorField& D,
const dimensionedScalar& phi const dimensionedScalar& phi

View File

@ -74,9 +74,9 @@ public:
virtual tmp<volScalarField> frictionalPressure virtual tmp<volScalarField> frictionalPressure
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& eta, const dimensionedScalar& eta,
const dimensionedScalar& p const dimensionedScalar& p
@ -84,9 +84,9 @@ public:
virtual tmp<volScalarField> frictionalPressurePrime virtual tmp<volScalarField> frictionalPressurePrime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& eta, const dimensionedScalar& eta,
const dimensionedScalar& p const dimensionedScalar& p
@ -94,8 +94,8 @@ public:
virtual tmp<volScalarField> muf virtual tmp<volScalarField> muf
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const volScalarField& pf, const volScalarField& pf,
const volSymmTensorField& D, const volSymmTensorField& D,
const dimensionedScalar& phi const dimensionedScalar& phi

View File

@ -70,9 +70,9 @@ Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer:: Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::
frictionalPressure frictionalPressure
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& eta, const dimensionedScalar& eta,
const dimensionedScalar& p const dimensionedScalar& p
@ -80,7 +80,7 @@ frictionalPressure
{ {
return return
dimensionedScalar("1e24", dimensionSet(1, -1, -2, 0, 0), 1e24) dimensionedScalar("1e24", dimensionSet(1, -1, -2, 0, 0), 1e24)
*pow(Foam::max(alpha - alphaMinFriction, scalar(0)), 10.0); *pow(Foam::max(alpha1 - alpha1MinFriction, scalar(0)), 10.0);
} }
@ -88,9 +88,9 @@ Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer:: Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::
frictionalPressurePrime frictionalPressurePrime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& eta, const dimensionedScalar& eta,
const dimensionedScalar& p const dimensionedScalar& p
@ -98,15 +98,15 @@ frictionalPressurePrime
{ {
return return
dimensionedScalar("1e25", dimensionSet(1, -1, -2, 0, 0), 1e25) dimensionedScalar("1e25", dimensionSet(1, -1, -2, 0, 0), 1e25)
*pow(Foam::max(alpha - alphaMinFriction, scalar(0)), 9.0); *pow(Foam::max(alpha1 - alpha1MinFriction, scalar(0)), 9.0);
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const volScalarField& pf, const volScalarField& pf,
const volSymmTensorField& D, const volSymmTensorField& D,
const dimensionedScalar& phi const dimensionedScalar& phi
@ -123,10 +123,10 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
IOobject IOobject
( (
"muf", "muf",
alpha.mesh().time().timeName(), alpha1.mesh().time().timeName(),
alpha.mesh() alpha1.mesh()
), ),
alpha.mesh(), alpha1.mesh(),
dimensionedScalar("muf", dimensionSet(1, -1, -1, 0, 0), 0.0) dimensionedScalar("muf", dimensionSet(1, -1, -1, 0, 0), 0.0)
) )
); );
@ -135,7 +135,7 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
forAll (D, celli) forAll (D, celli)
{ {
if (alpha[celli] > alphaMax.value() - 5e-2) if (alpha1[celli] > alpha1Max.value() - 5e-2)
{ {
muff[celli] = muff[celli] =
0.5*pf[celli]*sin(phi.value()) 0.5*pf[celli]*sin(phi.value())

View File

@ -74,9 +74,9 @@ public:
virtual tmp<volScalarField> frictionalPressure virtual tmp<volScalarField> frictionalPressure
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& eta, const dimensionedScalar& eta,
const dimensionedScalar& p const dimensionedScalar& p
@ -84,9 +84,9 @@ public:
virtual tmp<volScalarField> frictionalPressurePrime virtual tmp<volScalarField> frictionalPressurePrime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& n, const dimensionedScalar& n,
const dimensionedScalar& p const dimensionedScalar& p
@ -94,8 +94,8 @@ public:
virtual tmp<volScalarField> muf virtual tmp<volScalarField> muf
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const volScalarField& pf, const volScalarField& pf,
const volSymmTensorField& D, const volSymmTensorField& D,
const dimensionedScalar& phi const dimensionedScalar& phi

View File

@ -106,9 +106,9 @@ public:
virtual tmp<volScalarField> frictionalPressure virtual tmp<volScalarField> frictionalPressure
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& eta, const dimensionedScalar& eta,
const dimensionedScalar& p const dimensionedScalar& p
@ -116,9 +116,9 @@ public:
virtual tmp<volScalarField> frictionalPressurePrime virtual tmp<volScalarField> frictionalPressurePrime
( (
const volScalarField& alphaf, const volScalarField& alpha1f,
const dimensionedScalar& alphaMinFriction, const dimensionedScalar& alpha1MinFriction,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const dimensionedScalar& Fr, const dimensionedScalar& Fr,
const dimensionedScalar& eta, const dimensionedScalar& eta,
const dimensionedScalar& p const dimensionedScalar& p
@ -126,8 +126,8 @@ public:
virtual tmp<volScalarField> muf virtual tmp<volScalarField> muf
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax, const dimensionedScalar& alpha1Max,
const volScalarField& pf, const volScalarField& pf,
const volSymmTensorField& D, const volSymmTensorField& D,
const dimensionedScalar& phi const dimensionedScalar& phi

View File

@ -69,14 +69,14 @@ Foam::kineticTheoryModels::granularPressureModels::Lun::~Lun()
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::granularPressureModels::Lun::granularPressureCoeff Foam::kineticTheoryModels::granularPressureModels::Lun::granularPressureCoeff
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
{ {
return rhoa*alpha*(1.0 + 2.0*(1.0 + e)*alpha*g0); return rho1*alpha1*(1.0 + 2.0*(1.0 + e)*alpha1*g0);
} }
@ -84,14 +84,14 @@ Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::granularPressureModels::Lun:: Foam::kineticTheoryModels::granularPressureModels::Lun::
granularPressureCoeffPrime granularPressureCoeffPrime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const volScalarField& g0prime, const volScalarField& g0prime,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
{ {
return rhoa*(1.0 + alpha*(1.0 + e)*(4.0*g0 + 2.0*g0prime*alpha)); return rho1*(1.0 + alpha1*(1.0 + e)*(4.0*g0 + 2.0*g0prime*alpha1));
} }

View File

@ -74,18 +74,18 @@ public:
tmp<volScalarField> granularPressureCoeff tmp<volScalarField> granularPressureCoeff
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;
tmp<volScalarField> granularPressureCoeffPrime tmp<volScalarField> granularPressureCoeffPrime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const volScalarField& g0prime, const volScalarField& g0prime,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;
}; };

View File

@ -72,14 +72,14 @@ Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::granularPressureModels::SyamlalRogersOBrien:: Foam::kineticTheoryModels::granularPressureModels::SyamlalRogersOBrien::
granularPressureCoeff granularPressureCoeff
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
{ {
return 2.0*rhoa*(1.0 + e)*sqr(alpha)*g0; return 2.0*rho1*(1.0 + e)*sqr(alpha1)*g0;
} }
@ -87,14 +87,14 @@ Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::granularPressureModels::SyamlalRogersOBrien:: Foam::kineticTheoryModels::granularPressureModels::SyamlalRogersOBrien::
granularPressureCoeffPrime granularPressureCoeffPrime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const volScalarField& g0prime, const volScalarField& g0prime,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
{ {
return rhoa*alpha*(1.0 + e)*(4.0*g0 + 2.0*g0prime*alpha); return rho1*alpha1*(1.0 + e)*(4.0*g0 + 2.0*g0prime*alpha1);
} }

View File

@ -74,18 +74,18 @@ public:
tmp<volScalarField> granularPressureCoeff tmp<volScalarField> granularPressureCoeff
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;
tmp<volScalarField> granularPressureCoeffPrime tmp<volScalarField> granularPressureCoeffPrime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const volScalarField& g0prime, const volScalarField& g0prime,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;
}; };

View File

@ -107,19 +107,19 @@ public:
//- Granular pressure coefficient //- Granular pressure coefficient
virtual tmp<volScalarField> granularPressureCoeff virtual tmp<volScalarField> granularPressureCoeff
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const = 0; ) const = 0;
//- Derivative of the granular pressure coefficient //- Derivative of the granular pressure coefficient
virtual tmp<volScalarField> granularPressureCoeffPrime virtual tmp<volScalarField> granularPressureCoeffPrime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& g0, const volScalarField& g0,
const volScalarField& g0prime, const volScalarField& g0prime,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const dimensionedScalar& e const dimensionedScalar& e
) const = 0; ) const = 0;
}; };

View File

@ -33,28 +33,28 @@ License
Foam::kineticTheoryModel::kineticTheoryModel Foam::kineticTheoryModel::kineticTheoryModel
( (
const Foam::phaseModel& phase1, const Foam::phaseModel& phase1,
const Foam::volVectorField& Ub, const Foam::volVectorField& U2,
const Foam::volScalarField& alpha, const Foam::volScalarField& alpha1,
const Foam::dragModel& draga const Foam::dragModel& draga
) )
: :
phase1_(phase1), phase1_(phase1),
Ua_(phase1.U()), U1_(phase1.U()),
Ub_(Ub), U2_(U2),
alpha_(alpha), alpha1_(alpha1),
phia_(phase1.phi()), phi1_(phase1.phi()),
draga_(draga), draga_(draga),
rhoa_(phase1.rho()), rho1_(phase1.rho()),
nua_(phase1.nu()), nu1_(phase1.nu()),
kineticTheoryProperties_ kineticTheoryProperties_
( (
IOobject IOobject
( (
"kineticTheoryProperties", "kineticTheoryProperties",
Ua_.time().constant(), U1_.time().constant(),
Ua_.mesh(), U1_.mesh(),
IOobject::MUST_READ, IOobject::MUST_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
) )
@ -98,8 +98,8 @@ Foam::kineticTheoryModel::kineticTheoryModel
) )
), ),
e_(kineticTheoryProperties_.lookup("e")), e_(kineticTheoryProperties_.lookup("e")),
alphaMax_(kineticTheoryProperties_.lookup("alphaMax")), alpha1Max_(kineticTheoryProperties_.lookup("alpha1Max")),
alphaMinFriction_(kineticTheoryProperties_.lookup("alphaMinFriction")), alpha1MinFriction_(kineticTheoryProperties_.lookup("alpha1MinFriction")),
Fr_(kineticTheoryProperties_.lookup("Fr")), Fr_(kineticTheoryProperties_.lookup("Fr")),
eta_(kineticTheoryProperties_.lookup("eta")), eta_(kineticTheoryProperties_.lookup("eta")),
p_(kineticTheoryProperties_.lookup("p")), p_(kineticTheoryProperties_.lookup("p")),
@ -109,24 +109,24 @@ Foam::kineticTheoryModel::kineticTheoryModel
IOobject IOobject
( (
"Theta", "Theta",
Ua_.time().timeName(), U1_.time().timeName(),
Ua_.mesh(), U1_.mesh(),
IOobject::MUST_READ, IOobject::MUST_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
Ua_.mesh() U1_.mesh()
), ),
mua_ mu1_
( (
IOobject IOobject
( (
"mua", "mu1",
Ua_.time().timeName(), U1_.time().timeName(),
Ua_.mesh(), U1_.mesh(),
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
Ua_.mesh(), U1_.mesh(),
dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0) dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
), ),
lambda_ lambda_
@ -134,12 +134,12 @@ Foam::kineticTheoryModel::kineticTheoryModel
IOobject IOobject
( (
"lambda", "lambda",
Ua_.time().timeName(), U1_.time().timeName(),
Ua_.mesh(), U1_.mesh(),
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
Ua_.mesh(), U1_.mesh(),
dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0) dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
), ),
pa_ pa_
@ -147,12 +147,12 @@ Foam::kineticTheoryModel::kineticTheoryModel
IOobject IOobject
( (
"pa", "pa",
Ua_.time().timeName(), U1_.time().timeName(),
Ua_.mesh(), U1_.mesh(),
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
Ua_.mesh(), U1_.mesh(),
dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0) dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
), ),
kappa_ kappa_
@ -160,12 +160,12 @@ Foam::kineticTheoryModel::kineticTheoryModel
IOobject IOobject
( (
"kappa", "kappa",
Ua_.time().timeName(), U1_.time().timeName(),
Ua_.mesh(), U1_.mesh(),
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
Ua_.mesh(), U1_.mesh(),
dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0) dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
), ),
gs0_ gs0_
@ -173,12 +173,12 @@ Foam::kineticTheoryModel::kineticTheoryModel
IOobject IOobject
( (
"gs0", "gs0",
Ua_.time().timeName(), U1_.time().timeName(),
Ua_.mesh(), U1_.mesh(),
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
Ua_.mesh(), U1_.mesh(),
dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 1.0) dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 1.0)
) )
{} {}
@ -192,7 +192,7 @@ Foam::kineticTheoryModel::~kineticTheoryModel()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::kineticTheoryModel::solve(const volTensorField& gradUat) void Foam::kineticTheoryModel::solve(const volTensorField& gradU1t)
{ {
if (!kineticTheory_) if (!kineticTheory_)
{ {
@ -203,24 +203,24 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
volScalarField da_(phase1_.d()); volScalarField da_(phase1_.d());
surfaceScalarField phi(1.5*rhoa_*phia_*fvc::interpolate(alpha_)); surfaceScalarField phi(1.5*rho1_*phi1_*fvc::interpolate(alpha1_));
volTensorField dU(gradUat.T()); //fvc::grad(Ua_); volTensorField dU(gradU1t.T()); //fvc::grad(U1_);
volSymmTensorField D(symm(dU)); volSymmTensorField D(symm(dU));
// NB, drag = K*alpha*beta, // NB, drag = K*alpha1*alpha2,
// (the alpha and beta has been extracted from the drag function for // (the alpha1 and alpha2 has been extracted from the drag function for
// numerical reasons) // numerical reasons)
volScalarField Ur(mag(Ua_ - Ub_)); volScalarField Ur(mag(U1_ - U2_));
volScalarField betaPrim(alpha_*(1.0 - alpha_)*draga_.K(Ur)); volScalarField alpha2Prim(alpha1_*(1.0 - alpha1_)*draga_.K(Ur));
// Calculating the radial distribution function (solid volume fraction is // Calculating the radial distribution function (solid volume fraction is
// limited close to the packing limit, but this needs improvements) // limited close to the packing limit, but this needs improvements)
// The solution is higly unstable close to the packing limit. // The solution is higly unstable close to the packing limit.
gs0_ = radialModel_->g0 gs0_ = radialModel_->g0
( (
min(max(alpha_, scalar(1e-6)), alphaMax_ - 0.01), min(max(alpha1_, scalar(1e-6)), alpha1Max_ - 0.01),
alphaMax_ alpha1Max_
); );
// particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45) // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
@ -228,18 +228,18 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
( (
granularPressureModel_->granularPressureCoeff granularPressureModel_->granularPressureCoeff
( (
alpha_, alpha1_,
gs0_, gs0_,
rhoa_, rho1_,
e_ e_
) )
); );
// 'thermal' conductivity (Table 3.3, p. 49) // 'thermal' conductivity (Table 3.3, p. 49)
kappa_ = conductivityModel_->kappa(alpha_, Theta_, gs0_, rhoa_, da_, e_); kappa_ = conductivityModel_->kappa(alpha1_, Theta_, gs0_, rho1_, da_, e_);
// particle viscosity (Table 3.2, p.47) // particle viscosity (Table 3.2, p.47)
mua_ = viscosityModel_->mua(alpha_, Theta_, gs0_, rhoa_, da_, e_); mu1_ = viscosityModel_->mu1(alpha1_, Theta_, gs0_, rho1_, da_, e_);
dimensionedScalar Tsmall dimensionedScalar Tsmall
( (
@ -254,22 +254,22 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
// dissipation (Eq. 3.24, p.50) // dissipation (Eq. 3.24, p.50)
volScalarField gammaCoeff volScalarField gammaCoeff
( (
12.0*(1.0 - sqr(e_))*sqr(alpha_)*rhoa_*gs0_*(1.0/da_)*ThetaSqrt/sqrtPi 12.0*(1.0 - sqr(e_))*sqr(alpha1_)*rho1_*gs0_*(1.0/da_)*ThetaSqrt/sqrtPi
); );
// Eq. 3.25, p. 50 Js = J1 - J2 // Eq. 3.25, p. 50 Js = J1 - J2
volScalarField J1(3.0*betaPrim); volScalarField J1(3.0*alpha2Prim);
volScalarField J2 volScalarField J2
( (
0.25*sqr(betaPrim)*da_*sqr(Ur) 0.25*sqr(alpha2Prim)*da_*sqr(Ur)
/(max(alpha_, scalar(1e-6))*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt)) /(max(alpha1_, scalar(1e-6))*rho1_*sqrtPi*(ThetaSqrt + TsmallSqrt))
); );
// bulk viscosity p. 45 (Lun et al. 1984). // bulk viscosity p. 45 (Lun et al. 1984).
lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi; lambda_ = (4.0/3.0)*sqr(alpha1_)*rho1_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi;
// stress tensor, Definitions, Table 3.1, p. 43 // stress tensor, Definitions, Table 3.1, p. 43
volSymmTensorField tau(2.0*mua_*D + (lambda_ - (2.0/3.0)*mua_)*tr(D)*I); volSymmTensorField tau(2.0*mu1_*D + (lambda_ - (2.0/3.0)*mu1_)*tr(D)*I);
if (!equilibrium_) if (!equilibrium_)
{ {
@ -279,7 +279,7 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
// wrong sign infront of laplacian // wrong sign infront of laplacian
fvScalarMatrix ThetaEqn fvScalarMatrix ThetaEqn
( (
fvm::ddt(1.5*alpha_*rhoa_, Theta_) fvm::ddt(1.5*alpha1_*rho1_, Theta_)
+ fvm::div(phi, Theta_, "div(phi,Theta)") + fvm::div(phi, Theta_, "div(phi,Theta)")
== ==
fvm::SuSp(-((PsCoeff*I) && dU), Theta_) fvm::SuSp(-((PsCoeff*I) && dU), Theta_)
@ -297,40 +297,40 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
{ {
// equilibrium => dissipation == production // equilibrium => dissipation == production
// Eq. 4.14, p.82 // Eq. 4.14, p.82
volScalarField K1(2.0*(1.0 + e_)*rhoa_*gs0_); volScalarField K1(2.0*(1.0 + e_)*rho1_*gs0_);
volScalarField K3 volScalarField K3
( (
0.5*da_*rhoa_* 0.5*da_*rho1_*
( (
(sqrtPi/(3.0*(3.0-e_))) (sqrtPi/(3.0*(3.0-e_)))
*(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha_*gs0_) *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha1_*gs0_)
+1.6*alpha_*gs0_*(1.0 + e_)/sqrtPi +1.6*alpha1_*gs0_*(1.0 + e_)/sqrtPi
) )
); );
volScalarField K2 volScalarField K2
( (
4.0*da_*rhoa_*(1.0 + e_)*alpha_*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0 4.0*da_*rho1_*(1.0 + e_)*alpha1_*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0
); );
volScalarField K4(12.0*(1.0 - sqr(e_))*rhoa_*gs0_/(da_*sqrtPi)); volScalarField K4(12.0*(1.0 - sqr(e_))*rho1_*gs0_/(da_*sqrtPi));
volScalarField trD(tr(D)); volScalarField trD(tr(D));
volScalarField tr2D(sqr(trD)); volScalarField tr2D(sqr(trD));
volScalarField trD2(tr(D & D)); volScalarField trD2(tr(D & D));
volScalarField t1(K1*alpha_ + rhoa_); volScalarField t1(K1*alpha1_ + rho1_);
volScalarField l1(-t1*trD); volScalarField l1(-t1*trD);
volScalarField l2(sqr(t1)*tr2D); volScalarField l2(sqr(t1)*tr2D);
volScalarField l3 volScalarField l3
( (
4.0 4.0
*K4 *K4
*max(alpha_, scalar(1e-6)) *max(alpha1_, scalar(1e-6))
*(2.0*K3*trD2 + K2*tr2D) *(2.0*K3*trD2 + K2*tr2D)
); );
Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha_ + 1.0e-4)*K4)); Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha1_ + 1.0e-4)*K4));
} }
Theta_.max(1.0e-15); Theta_.max(1.0e-15);
@ -340,9 +340,9 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
( (
frictionalStressModel_->frictionalPressure frictionalStressModel_->frictionalPressure
( (
alpha_, alpha1_,
alphaMinFriction_, alpha1MinFriction_,
alphaMax_, alpha1Max_,
Fr_, Fr_,
eta_, eta_,
p_ p_
@ -362,8 +362,8 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
( (
frictionalStressModel_->muf frictionalStressModel_->muf
( (
alpha_, alpha1_,
alphaMax_, alpha1Max_,
pf, pf,
D, D,
phi_ phi_
@ -371,16 +371,16 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
); );
// add frictional stress // add frictional stress
mua_ += muf; mu1_ += muf;
mua_.min(1.0e+2); mu1_.min(1.0e+2);
mua_.max(0.0); mu1_.max(0.0);
Info<< "kinTheory: max(Theta) = " << max(Theta_).value() << endl; Info<< "kinTheory: max(Theta) = " << max(Theta_).value() << endl;
volScalarField ktn(mua_/rhoa_); volScalarField ktn(mu1_/rho1_);
Info<< "kinTheory: min(nua) = " << min(ktn).value() Info<< "kinTheory: min(nu1) = " << min(ktn).value()
<< ", max(nua) = " << max(ktn).value() << endl; << ", max(nu1) = " << max(ktn).value() << endl;
Info<< "kinTheory: min(pa) = " << min(pa_).value() Info<< "kinTheory: min(pa) = " << min(pa_).value()
<< ", max(pa) = " << max(pa_).value() << endl; << ", max(pa) = " << max(pa_).value() << endl;

View File

@ -57,15 +57,15 @@ class kineticTheoryModel
// Private data // Private data
const phaseModel& phase1_; const phaseModel& phase1_;
const volVectorField& Ua_; const volVectorField& U1_;
const volVectorField& Ub_; const volVectorField& U2_;
const volScalarField& alpha_; const volScalarField& alpha1_;
const surfaceScalarField& phia_; const surfaceScalarField& phi1_;
const dragModel& draga_; const dragModel& draga_;
const dimensionedScalar& rhoa_; const dimensionedScalar& rho1_;
const dimensionedScalar& nua_; const dimensionedScalar& nu1_;
//- dictionary holding the modeling info //- dictionary holding the modeling info
IOdictionary kineticTheoryProperties_; IOdictionary kineticTheoryProperties_;
@ -92,10 +92,10 @@ class kineticTheoryModel
const dimensionedScalar e_; const dimensionedScalar e_;
//- maximum packing //- maximum packing
const dimensionedScalar alphaMax_; const dimensionedScalar alpha1Max_;
//- min value for which the frictional stresses are zero //- min value for which the frictional stresses are zero
const dimensionedScalar alphaMinFriction_; const dimensionedScalar alpha1MinFriction_;
//- material constant for frictional normal stress //- material constant for frictional normal stress
const dimensionedScalar Fr_; const dimensionedScalar Fr_;
@ -113,7 +113,7 @@ class kineticTheoryModel
volScalarField Theta_; volScalarField Theta_;
//- The granular viscosity //- The granular viscosity
volScalarField mua_; volScalarField mu1_;
//- The granular bulk viscosity //- The granular bulk viscosity
volScalarField lambda_; volScalarField lambda_;
@ -145,8 +145,8 @@ public:
kineticTheoryModel kineticTheoryModel
( (
const phaseModel& phase1, const phaseModel& phase1,
const volVectorField& Ub, const volVectorField& U2,
const volScalarField& alpha, const volScalarField& alpha1,
const dragModel& draga const dragModel& draga
); );
@ -157,16 +157,16 @@ public:
// Member Functions // Member Functions
void solve(const volTensorField& gradUat); void solve(const volTensorField& gradU1t);
bool on() const bool on() const
{ {
return kineticTheory_; return kineticTheory_;
} }
const volScalarField& mua() const const volScalarField& mu1() const
{ {
return mua_; return mu1_;
} }
const volScalarField& pa() const const volScalarField& pa() const

View File

@ -69,30 +69,30 @@ Foam::kineticTheoryModels::radialModels::CarnahanStarling::~CarnahanStarling()
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0 Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const ) const
{ {
return return
1.0/(1.0 - alpha) 1.0/(1.0 - alpha1)
+ 3.0*alpha/(2.0*sqr(1.0 - alpha)) + 3.0*alpha1/(2.0*sqr(1.0 - alpha1))
+ sqr(alpha)/(2.0*pow(1.0 - alpha, 3)); + sqr(alpha1)/(2.0*pow(1.0 - alpha1, 3));
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0prime Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0prime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const ) const
{ {
return return
- alpha/sqr(1.0 - alpha) - alpha1/sqr(1.0 - alpha1)
+ (3.0*(1.0 - alpha) + 6.0*sqr(alpha))/(2.0*(1.0 - alpha)) + (3.0*(1.0 - alpha1) + 6.0*sqr(alpha1))/(2.0*(1.0 - alpha1))
+ (2.0*alpha*(1.0 - alpha) + 3.0*pow(alpha, 3)) + (2.0*alpha1*(1.0 - alpha1) + 3.0*pow(alpha1, 3))
/(2.0*pow(1.0 - alpha, 4)); /(2.0*pow(1.0 - alpha1, 4));
} }

View File

@ -75,14 +75,14 @@ public:
tmp<volScalarField> g0 tmp<volScalarField> g0
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const; ) const;
tmp<volScalarField> g0prime tmp<volScalarField> g0prime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const; ) const;
}; };

View File

@ -69,24 +69,24 @@ Foam::kineticTheoryModels::radialModels::Gidaspow::~Gidaspow()
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::Gidaspow::g0 Foam::kineticTheoryModels::radialModels::Gidaspow::g0
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const ) const
{ {
return 0.6/(1.0 - pow(alpha/alphaMax, 1.0/3.0)); return 0.6/(1.0 - pow(alpha1/alpha1Max, 1.0/3.0));
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::Gidaspow::g0prime Foam::kineticTheoryModels::radialModels::Gidaspow::g0prime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const ) const
{ {
return return
(-1.0/5.0)*pow(alpha/alphaMax, -2.0/3.0) (-1.0/5.0)*pow(alpha1/alpha1Max, -2.0/3.0)
/(alphaMax*sqr(1.0 - pow(alpha/alphaMax, 1.0/3.0))); /(alpha1Max*sqr(1.0 - pow(alpha1/alpha1Max, 1.0/3.0)));
} }

View File

@ -74,14 +74,14 @@ public:
tmp<volScalarField> g0 tmp<volScalarField> g0
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const; ) const;
tmp<volScalarField> g0prime tmp<volScalarField> g0prime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const; ) const;
}; };

View File

@ -69,23 +69,23 @@ Foam::kineticTheoryModels::radialModels::LunSavage::~LunSavage()
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::LunSavage::g0 Foam::kineticTheoryModels::radialModels::LunSavage::g0
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const ) const
{ {
return pow(1.0 - alpha/alphaMax, -2.5*alphaMax); return pow(1.0 - alpha1/alpha1Max, -2.5*alpha1Max);
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::LunSavage::g0prime Foam::kineticTheoryModels::radialModels::LunSavage::g0prime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const ) const
{ {
return 2.5*alphaMax*alpha*pow(1.0 - alpha, -1.0 - 2.5*alphaMax); return 2.5*alpha1Max*alpha1*pow(1.0 - alpha1, -1.0 - 2.5*alpha1Max);
} }

View File

@ -74,14 +74,14 @@ public:
tmp<volScalarField> g0 tmp<volScalarField> g0
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const; ) const;
tmp<volScalarField> g0prime tmp<volScalarField> g0prime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const; ) const;
}; };

View File

@ -69,24 +69,24 @@ Foam::kineticTheoryModels::radialModels::SinclairJackson::~SinclairJackson()
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::SinclairJackson::g0 Foam::kineticTheoryModels::radialModels::SinclairJackson::g0
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const ) const
{ {
return 1.0/(1.0 - pow(alpha/alphaMax, 1.0/3.0)); return 1.0/(1.0 - pow(alpha1/alpha1Max, 1.0/3.0));
} }
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const ) const
{ {
return return
(-1.0/3.0)*pow(alpha/alphaMax, -2.0/3.0) (-1.0/3.0)*pow(alpha1/alpha1Max, -2.0/3.0)
/(alphaMax*sqr(1.0 - pow(alpha/alphaMax, 1.0/3.0))); /(alpha1Max*sqr(1.0 - pow(alpha1/alpha1Max, 1.0/3.0)));
} }

View File

@ -74,14 +74,14 @@ public:
tmp<volScalarField> g0 tmp<volScalarField> g0
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const; ) const;
tmp<volScalarField> g0prime tmp<volScalarField> g0prime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const; ) const;
}; };

View File

@ -107,15 +107,15 @@ public:
//- Radial distribution function //- Radial distribution function
virtual tmp<volScalarField> g0 virtual tmp<volScalarField> g0
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const = 0; ) const = 0;
//- Derivative of the radial distribution function //- Derivative of the radial distribution function
virtual tmp<volScalarField> g0prime virtual tmp<volScalarField> g0prime
( (
const volScalarField& alpha, const volScalarField& alpha1,
const dimensionedScalar& alphaMax const dimensionedScalar& alpha1Max
) const = 0; ) const = 0;
}; };

View File

@ -62,23 +62,23 @@ Foam::kineticTheoryModels::viscosityModels::Gidaspow::~Gidaspow()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::viscosityModels::Gidaspow::mua Foam::kineticTheoryModels::viscosityModels::Gidaspow::mu1
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
{ {
const scalar sqrtPi = sqrt(constant::mathematical::pi); const scalar sqrtPi = sqrt(constant::mathematical::pi);
return rhoa*da*sqrt(Theta)* return rho1*da*sqrt(Theta)*
( (
(4.0/5.0)*sqr(alpha)*g0*(1.0 + e)/sqrtPi (4.0/5.0)*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
+ (1.0/15.0)*sqrtPi*g0*(1.0 + e)*sqr(alpha) + (1.0/15.0)*sqrtPi*g0*(1.0 + e)*sqr(alpha1)
+ (1.0/6.0)*sqrtPi*alpha + (1.0/6.0)*sqrtPi*alpha1
+ (10.0/96.0)*sqrtPi/((1.0 + e)*g0) + (10.0/96.0)*sqrtPi/((1.0 + e)*g0)
); );
} }

View File

@ -71,12 +71,12 @@ public:
// Member functions // Member functions
tmp<volScalarField> mua tmp<volScalarField> mu1
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;

View File

@ -70,12 +70,12 @@ Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::~HrenyaSinclair()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::mua Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::mu1
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
@ -84,14 +84,14 @@ Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::mua
volScalarField lamda volScalarField lamda
( (
scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_ scalar(1) + da/(6.0*sqrt(2.0)*(alpha1 + scalar(1.0e-5)))/L_
); );
return rhoa*da*sqrt(Theta)* return rho1*da*sqrt(Theta)*
( (
(4.0/5.0)*sqr(alpha)*g0*(1.0 + e)/sqrtPi (4.0/5.0)*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
+ (1.0/15.0)*sqrtPi*g0*(1.0 + e)*(3.0*e - 1)*sqr(alpha)/(3.0-e) + (1.0/15.0)*sqrtPi*g0*(1.0 + e)*(3.0*e - 1)*sqr(alpha1)/(3.0-e)
+ (1.0/6.0)*sqrtPi*alpha*(0.5*lamda + 0.25*(3.0*e - 1.0)) + (1.0/6.0)*sqrtPi*alpha1*(0.5*lamda + 0.25*(3.0*e - 1.0))
/(0.5*(3.0 - e)*lamda) /(0.5*(3.0 - e)*lamda)
+ (10/96.0)*sqrtPi/((1.0 + e)*0.5*(3.0 - e)*g0*lamda) + (10/96.0)*sqrtPi/((1.0 + e)*0.5*(3.0 - e)*g0*lamda)
); );

View File

@ -79,12 +79,12 @@ public:
// Member functions // Member functions
tmp<volScalarField> mua tmp<volScalarField> mu1
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;

View File

@ -62,23 +62,23 @@ Foam::kineticTheoryModels::viscosityModels::Syamlal::~Syamlal()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::viscosityModels::Syamlal::mua Foam::kineticTheoryModels::viscosityModels::Syamlal::mu1
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
{ {
const scalar sqrtPi = sqrt(constant::mathematical::pi); const scalar sqrtPi = sqrt(constant::mathematical::pi);
return rhoa*da*sqrt(Theta)* return rho1*da*sqrt(Theta)*
( (
(4.0/5.0)*sqr(alpha)*g0*(1.0 + e)/sqrtPi (4.0/5.0)*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
+ (1.0/15.0)*sqrtPi*g0*(1.0 + e)*(3.0*e - 1.0)*sqr(alpha)/(3.0 - e) + (1.0/15.0)*sqrtPi*g0*(1.0 + e)*(3.0*e - 1.0)*sqr(alpha1)/(3.0 - e)
+ (1.0/6.0)*alpha*sqrtPi/(3.0 - e) + (1.0/6.0)*alpha1*sqrtPi/(3.0 - e)
); );
} }

View File

@ -72,12 +72,12 @@ public:
// Member functions // Member functions
tmp<volScalarField> mua tmp<volScalarField> mu1
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;

View File

@ -53,12 +53,12 @@ Foam::kineticTheoryModels::noneViscosity::~noneViscosity()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModels::noneViscosity::mua Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModels::noneViscosity::mu1
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
@ -68,7 +68,7 @@ Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModels::noneViscosity::mua
"0", "0",
dimensionSet(1, -1, -1, 0, 0, 0, 0), dimensionSet(1, -1, -1, 0, 0, 0, 0),
0.0 0.0
)*alpha; )*alpha1;
} }

View File

@ -70,12 +70,12 @@ public:
// Member functions // Member functions
tmp<volScalarField> mua tmp<volScalarField> mu1
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;

View File

@ -106,12 +106,12 @@ public:
// Member Functions // Member Functions
virtual tmp<volScalarField> mua virtual tmp<volScalarField> mu1
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const = 0; ) const = 0;

View File

@ -1,13 +1,13 @@
if (packingLimiter) if (packingLimiter)
{ {
// Calculating exceeding volume fractions // Calculating exceeding volume fractions
volScalarField alphaEx(max(alpha1 - alphaMax, scalar(0))); volScalarField alpha1Ex(max(alpha1 - alpha1Max, scalar(0)));
// Finding neighbouring cells of the whole domain // Finding neighbouring cells of the whole domain
labelListList neighbour = mesh.cellCells(); labelListList neighbour = mesh.cellCells();
scalarField cellVolumes(mesh.cellVolumes()); scalarField cellVolumes(mesh.cellVolumes());
forAll (alphaEx, celli) forAll (alpha1Ex, celli)
{ {
// Finding the labels of the neighbouring cells // Finding the labels of the neighbouring cells
labelList neighbourCell = neighbour[celli]; labelList neighbourCell = neighbour[celli];
@ -27,10 +27,10 @@
} }
neighboursEx += neighboursEx +=
alphaEx[neighbourCell[cellj]]*cellVolumes[celli] alpha1Ex[neighbourCell[cellj]]*cellVolumes[celli]
/neighboursNeighbourCellVolumes; /neighboursNeighbourCellVolumes;
} }
alpha1[celli] += neighboursEx - alphaEx[celli]; alpha1[celli] += neighboursEx - alpha1Ex[celli];
} }
} }

View File

@ -25,7 +25,7 @@ Class
Foam::diameterModel Foam::diameterModel
Description Description
Abstract base-class for dispersed-phase particle diameter models. A2stract base-class for dispersed-phase particle diameter models.
SourceFiles SourceFiles
diameterModel.C diameterModel.C

View File

@ -15,9 +15,9 @@
readScalar(ppProperties.lookup("preAlphaExp")) readScalar(ppProperties.lookup("preAlphaExp"))
); );
scalar alphaMax scalar alpha1Max
( (
readScalar(ppProperties.lookup("alphaMax")) readScalar(ppProperties.lookup("alpha1Max"))
); );
scalar expMax scalar expMax

View File

@ -15,10 +15,10 @@ if (turbulence)
fvScalarMatrix epsEqn fvScalarMatrix epsEqn
( (
fvm::ddt(alpha2, epsilon) fvm::ddt(alpha2, epsilon)
+ fvm::div(alphaPhi2, epsilon) + fvm::div(alpha1Phi2, epsilon)
- fvm::laplacian - fvm::laplacian
( (
alphaEps*nuEff2, epsilon, alpha1Eps*nuEff2, epsilon,
"laplacian(DepsilonEff,epsilon)" "laplacian(DepsilonEff,epsilon)"
) )
== ==
@ -40,10 +40,10 @@ if (turbulence)
fvScalarMatrix kEqn fvScalarMatrix kEqn
( (
fvm::ddt(alpha2, k) fvm::ddt(alpha2, k)
+ fvm::div(alphaPhi2, k) + fvm::div(alpha1Phi2, k)
- fvm::laplacian - fvm::laplacian
( (
alphak*nuEff2, k, alpha1k*nuEff2, k,
"laplacian(DkEff,k)" "laplacian(DkEff,k)"
) )
== ==

View File

@ -70,12 +70,12 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::Ergun::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(phase2_, scalar(1.0e-6))); volScalarField alpha2(max(phase2_, scalar(1.0e-6)));
return return
150.0*phase1_*phase2_.nu()*phase2_.rho() 150.0*phase1_*phase2_.nu()*phase2_.rho()
/sqr(beta*phase1_.d()) /sqr(alpha2*phase1_.d())
+ 1.75*phase2_.rho()*Ur/(beta*phase1_.d()); + 1.75*phase2_.rho()*Ur/(alpha2*phase1_.d());
} }

View File

@ -70,9 +70,9 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::Gibilaro::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(phase2_, scalar(1.0e-6))); volScalarField alpha2(max(phase2_, scalar(1.0e-6)));
volScalarField bp(pow(beta, -2.8)); volScalarField bp(pow(alpha2, -2.8));
volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(alpha2*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
return (17.3/Re + scalar(0.336))*phase2_.rho()*Ur*bp/phase1_.d(); return (17.3/Re + scalar(0.336))*phase2_.rho()*Ur*bp/phase1_.d();
} }

View File

@ -70,9 +70,9 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(phase2_, scalar(1.0e-6))); volScalarField alpha2(max(phase2_, scalar(1.0e-6)));
volScalarField d = phase1_.d(); volScalarField d = phase1_.d();
volScalarField bp(pow(beta, -2.65)); volScalarField bp(pow(alpha2, -2.65));
volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds volScalarField Cds
@ -84,12 +84,12 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
// Wen and Yu (1966) // Wen and Yu (1966)
return return
( (
pos(beta - 0.8) pos(alpha2 - 0.8)
*(0.75*Cds*phase2_.rho()*Ur*bp/d) *(0.75*Cds*phase2_.rho()*Ur*bp/d)
+ neg(beta - 0.8) + neg(alpha2 - 0.8)
*( *(
150.0*phase1_*phase2_.nu()*phase2_.rho()/(sqr(beta*d)) 150.0*phase1_*phase2_.nu()*phase2_.rho()/(sqr(alpha2*d))
+ 1.75*phase2_.rho()*Ur/(beta*d) + 1.75*phase2_.rho()*Ur/(alpha2*d)
) )
); );
} }

View File

@ -70,10 +70,10 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(phase2_, scalar(1e-6))); volScalarField alpha2(max(phase2_, scalar(1e-6)));
volScalarField bp(pow(beta, -2.65)); volScalarField bp(pow(alpha2, -2.65));
volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(alpha2*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds volScalarField Cds
( (
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re) neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)

View File

@ -70,12 +70,12 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(phase2_, scalar(1.0e-6))); volScalarField alpha2(max(phase2_, scalar(1.0e-6)));
volScalarField A(pow(beta, 4.14)); volScalarField A(pow(alpha2, 4.14));
volScalarField B volScalarField B
( (
neg(beta - 0.85)*(0.8*pow(beta, 1.28)) neg(alpha2 - 0.85)*(0.8*pow(alpha2, 1.28))
+ pos(beta - 0.85)*(pow(beta, 2.65)) + pos(alpha2 - 0.85)*(pow(alpha2, 2.65))
); );
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));

View File

@ -70,8 +70,8 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
const volScalarField& Ur const volScalarField& Ur
) const ) const
{ {
volScalarField beta(max(phase2_, scalar(1.0e-6))); volScalarField alpha2(max(phase2_, scalar(1.0e-6)));
volScalarField bp(pow(beta, -2.65)); volScalarField bp(pow(alpha2, -2.65));
volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
volScalarField Cds volScalarField Cds

View File

@ -129,12 +129,13 @@ public:
} }
//- the dragfunction K used in the momentum eq. //- the dragfunction K used in the momentum eq.
// ddt(alpha*rhoa*Ua) + ... = ... alpha*beta*K*(Ua-Ub) // ddt(alpha1*rho1*U1) + ... = ... alpha1*alpha2*K*(U1-U2)
// ddt(beta*rhob*Ub) + ... = ... alpha*beta*K*(Ub-Ua) // ddt(alpha2*rho2*U2) + ... = ... alpha1*alpha2*K*(U2-U1)
// ********************************** NB! ***************************** // ********************************** NB! *****************************
// for numerical reasons alpha and beta has been // for numerical reasons alpha1 and alpha2 has been
// extracted from the dragFunction K, // extracted from the dragFunction K,
// so you MUST divide K by alpha*beta when implemnting the drag function // so you MUST divide K by alpha1*alpha2 when implemnting the drag
// function
// ********************************** NB! ***************************** // ********************************** NB! *****************************
virtual tmp<volScalarField> K(const volScalarField& Ur) const = 0; virtual tmp<volScalarField> K(const volScalarField& Ur) const = 0;
}; };

View File

@ -49,12 +49,12 @@ namespace heatTransferModels
Foam::heatTransferModels::RanzMarshall::RanzMarshall Foam::heatTransferModels::RanzMarshall::RanzMarshall
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
heatTransferModel(interfaceDict, alpha, phase1, phase2) heatTransferModel(interfaceDict, alpha1, phase1, phase2)
{} {}

View File

@ -64,7 +64,7 @@ public:
RanzMarshall RanzMarshall
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );

View File

@ -39,13 +39,13 @@ namespace Foam
Foam::heatTransferModel::heatTransferModel Foam::heatTransferModel::heatTransferModel
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
: :
interfaceDict_(interfaceDict), interfaceDict_(interfaceDict),
alpha_(alpha), alpha1_(alpha1),
phase1_(phase1), phase1_(phase1),
phase2_(phase2) phase2_(phase2)
{} {}

View File

@ -55,7 +55,7 @@ protected:
// Protected data // Protected data
const dictionary& interfaceDict_; const dictionary& interfaceDict_;
const volScalarField& alpha_; const volScalarField& alpha1_;
const phaseModel& phase1_; const phaseModel& phase1_;
const phaseModel& phase2_; const phaseModel& phase2_;
@ -75,11 +75,11 @@ public:
dictionary, dictionary,
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
), ),
(interfaceDict, alpha, phase1, phase2) (interfaceDict, alpha1, phase1, phase2)
); );
@ -88,7 +88,7 @@ public:
heatTransferModel heatTransferModel
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );
@ -103,7 +103,7 @@ public:
static autoPtr<heatTransferModel> New static autoPtr<heatTransferModel> New
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
); );
@ -112,12 +112,12 @@ public:
// Member Functions // Member Functions
//- the heat-transfer function K used in the enthalpy eq. //- the heat-transfer function K used in the enthalpy eq.
// ddt(alpha*rhoa*ha) + ... = ... alpha*beta*K*(Ta - Tb) // ddt(alpha1*rho1*ha) + ... = ... alpha1*alpha2*K*(Ta - Tb)
// ddt(beta*rhob*hb) + ... = ... alpha*beta*K*(Tb - Ta) // ddt(alpha2*rho2*hb) + ... = ... alpha1*alpha2*K*(Tb - Ta)
// ********************************** NB! ***************************** // ********************************** NB! *****************************
// for numerical reasons alpha and beta has been // for numerical reasons alpha1 and alpha2 has been
// extracted from the heat-transfer function K, // extracted from the heat-transfer function K,
// so you MUST divide K by alpha*beta when implementing the // so you MUST divide K by alpha1*alpha2 when implementing the
// heat-transfer function // heat-transfer function
// ********************************** NB! ***************************** // ********************************** NB! *****************************
virtual tmp<volScalarField> K(const volScalarField& Ur) const = 0; virtual tmp<volScalarField> K(const volScalarField& Ur) const = 0;

View File

@ -30,7 +30,7 @@ License
Foam::autoPtr<Foam::heatTransferModel> Foam::heatTransferModel::New Foam::autoPtr<Foam::heatTransferModel> Foam::heatTransferModel::New
( (
const dictionary& interfaceDict, const dictionary& interfaceDict,
const volScalarField& alpha, const volScalarField& alpha1,
const phaseModel& phase1, const phaseModel& phase1,
const phaseModel& phase2 const phaseModel& phase2
) )
@ -58,7 +58,7 @@ Foam::autoPtr<Foam::heatTransferModel> Foam::heatTransferModel::New
<< exit(FatalError); << exit(FatalError);
} }
return cstrIter()(interfaceDict, alpha, phase1, phase2); return cstrIter()(interfaceDict, alpha1, phase1, phase2);
} }

View File

@ -70,21 +70,21 @@ Foam::kineticTheoryModels::conductivityModels::Gidaspow::~Gidaspow()
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::conductivityModels::Gidaspow::kappa Foam::kineticTheoryModels::conductivityModels::Gidaspow::kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
{ {
const scalar sqrtPi = sqrt(constant::mathematical::pi); const scalar sqrtPi = sqrt(constant::mathematical::pi);
return rhoa*da*sqrt(Theta)* return rho1*da*sqrt(Theta)*
( (
2.0*sqr(alpha)*g0*(1.0 + e)/sqrtPi 2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
+ (9.0/8.0)*sqrtPi*g0*0.5*(1.0 + e)*sqr(alpha) + (9.0/8.0)*sqrtPi*g0*0.5*(1.0 + e)*sqr(alpha1)
+ (15.0/16.0)*sqrtPi*alpha + (15.0/16.0)*sqrtPi*alpha1
+ (25.0/64.0)*sqrtPi/((1.0 + e)*g0) + (25.0/64.0)*sqrtPi/((1.0 + e)*g0)
); );
} }

View File

@ -74,10 +74,10 @@ public:
tmp<volScalarField> kappa tmp<volScalarField> kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;

View File

@ -73,10 +73,10 @@ Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::
Foam::tmp<Foam::volScalarField> Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const ) const
@ -85,15 +85,15 @@ Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa
volScalarField lamda volScalarField lamda
( (
scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_ scalar(1) + da/(6.0*sqrt(2.0)*(alpha1 + scalar(1.0e-5)))/L_
); );
return rhoa*da*sqrt(Theta)* return rho1*da*sqrt(Theta)*
( (
2.0*sqr(alpha)*g0*(1.0 + e)/sqrtPi 2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
+ (9.0/8.0)*sqrtPi*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha) + (9.0/8.0)*sqrtPi*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha1)
/(49.0/16.0 - 33.0*e/16.0) /(49.0/16.0 - 33.0*e/16.0)
+ (15.0/16.0)*sqrtPi*alpha*(0.5*sqr(e) + 0.25*e - 0.75 + lamda) + (15.0/16.0)*sqrtPi*alpha1*(0.5*sqr(e) + 0.25*e - 0.75 + lamda)
/((49.0/16.0 - 33.0*e/16.0)*lamda) /((49.0/16.0 - 33.0*e/16.0)*lamda)
+ (25.0/64.0)*sqrtPi + (25.0/64.0)*sqrtPi
/((1.0 + e)*(49.0/16.0 - 33.0*e/16.0)*lamda*g0) /((1.0 + e)*(49.0/16.0 - 33.0*e/16.0)*lamda*g0)

View File

@ -79,10 +79,10 @@ public:
tmp<volScalarField> kappa tmp<volScalarField> kappa
( (
const volScalarField& alpha, const volScalarField& alpha1,
const volScalarField& Theta, const volScalarField& Theta,
const volScalarField& g0, const volScalarField& g0,
const dimensionedScalar& rhoa, const dimensionedScalar& rho1,
const volScalarField& da, const volScalarField& da,
const dimensionedScalar& e const dimensionedScalar& e
) const; ) const;

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