zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

COMP: avoid ambiguous construct from tmp - solvers/ multiphase

Browse Source
This commit is contained in:
Mark Olesen
2010-12-21 09:53:19 +01:00
parent 5a8f925221
commit 078c427594
68 changed files with 399 additions and 286 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
applications/solvers/multiphase/bubbleFoam/UEqns.H
View File

@ -2,12 +2,14 @@ fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime);
fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
{
volTensorField Rca = -nuEffa*(fvc::grad(Ua)().T());
volTensorField Rca(-nuEffa*(fvc::grad(Ua)().T()));
Rca = Rca + (2.0/3.0)*sqr(Ct)*I*k - (2.0/3.0)*I*tr(Rca);
surfaceScalarField phiRa =
surfaceScalarField phiRa
(
- fvc::interpolate(nuEffa)
*mesh.magSf()*fvc::snGrad(alpha)/fvc::interpolate(alpha + scalar(0.001));
*mesh.magSf()*fvc::snGrad(alpha)/fvc::interpolate(alpha + scalar(0.001))
);
UaEqn =
(
@ -34,12 +36,14 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
UaEqn.relax();
volTensorField Rcb = -nuEffb*fvc::grad(Ub)().T();
volTensorField Rcb(-nuEffb*fvc::grad(Ub)().T());
Rcb = Rcb + (2.0/3.0)*I*k - (2.0/3.0)*I*tr(Rcb);
surfaceScalarField phiRb =
surfaceScalarField phiRb
(
- fvc::interpolate(nuEffb)
*mesh.magSf()*fvc::snGrad(beta)/fvc::interpolate(beta + scalar(0.001));
*mesh.magSf()*fvc::snGrad(beta)/fvc::interpolate(beta + scalar(0.001))
);
UbEqn =
(

2
applications/solvers/multiphase/bubbleFoam/alphaEqn.H
View File

@ -1,7 +1,7 @@
{
word scheme("div(phi,alpha)");
surfaceScalarField phir = phia - phib;
surfaceScalarField phir(phia - phib);
Info<< "Max Ur Courant Number = "
<< (

12
applications/solvers/multiphase/bubbleFoam/createFields.H
View File

@ -171,15 +171,19 @@
Info<< "Calculating field DDtUa and DDtUb\n" << endl;
volVectorField DDtUa =
volVectorField DDtUa
(
fvc::ddt(Ua)
+ fvc::div(phia, Ua)
- fvc::div(phia)*Ua;
- fvc::div(phia)*Ua
);
volVectorField DDtUb =
volVectorField DDtUb
(
fvc::ddt(Ub)
+ fvc::div(phib, Ub)
- fvc::div(phib)*Ub;
- fvc::div(phib)*Ub
);
Info<< "Calculating field g.h\n" << endl;

2
applications/solvers/multiphase/bubbleFoam/kEpsilon.H
View File

@ -6,7 +6,7 @@ if (turbulence)
}
tmp<volTensorField> tgradUb = fvc::grad(Ub);
volScalarField G = 2*nutb*(tgradUb() && dev(symm(tgradUb())));
volScalarField G(2*nutb*(tgradUb() && dev(symm(tgradUb()))));
tgradUb.clear();
#include "wallFunctions.H"

21
applications/solvers/multiphase/bubbleFoam/liftDragCoeffs.H
View File

@ -1,11 +1,15 @@
volVectorField Ur = Ua - Ub;
volScalarField magUr = mag(Ur);
volVectorField Ur(Ua - Ub);
volScalarField magUr(mag(Ur));
volScalarField CdaMagUr =
(24.0*nub/da)*(scalar(1) + 0.15*pow(da*magUr/nub, 0.687));
volScalarField CdaMagUr
(
(24.0*nub/da)*(scalar(1) + 0.15*pow(da*magUr/nub, 0.687))
);
volScalarField CdbMagUr =
(24.0*nua/db)*(scalar(1) + 0.15*pow(db*magUr/nua, 0.687));
volScalarField CdbMagUr
(
(24.0*nua/db)*(scalar(1) + 0.15*pow(db*magUr/nua, 0.687))
);
volScalarField dragCoef
(
@ -13,4 +17,7 @@
0.75*(beta*rhob*CdaMagUr/da + alpha*rhoa*CdbMagUr/db)
);
volVectorField liftCoeff = Cl*(beta*rhob + alpha*rhoa)*(Ur ^ fvc::curl(U));
volVectorField liftCoeff
(
Cl*(beta*rhob + alpha*rhoa)*(Ur ^ fvc::curl(U))
);

26
applications/solvers/multiphase/bubbleFoam/pEqn.H
View File

@ -1,20 +1,24 @@
{
surfaceScalarField alphaf = fvc::interpolate(alpha);
surfaceScalarField betaf = scalar(1) - alphaf;
surfaceScalarField alphaf(fvc::interpolate(alpha));
surfaceScalarField betaf(scalar(1) - alphaf);
volScalarField rUaA = 1.0/UaEqn.A();
volScalarField rUbA = 1.0/UbEqn.A();
volScalarField rUaA(1.0/UaEqn.A());
volScalarField rUbA(1.0/UbEqn.A());
surfaceScalarField rUaAf = fvc::interpolate(rUaA);
surfaceScalarField rUbAf = fvc::interpolate(rUbA);
surfaceScalarField rUaAf(fvc::interpolate(rUaA));
surfaceScalarField rUbAf(fvc::interpolate(rUbA));
Ua = rUaA*UaEqn.H();
Ub = rUbA*UbEqn.H();
surfaceScalarField phiDraga =
fvc::interpolate(beta/rhoa*dragCoef*rUaA)*phib + rUaAf*(g & mesh.Sf());
surfaceScalarField phiDragb =
fvc::interpolate(alpha/rhob*dragCoef*rUbA)*phia + rUbAf*(g & mesh.Sf());
surfaceScalarField phiDraga
(
fvc::interpolate(beta/rhoa*dragCoef*rUaA)*phib + rUaAf*(g & mesh.Sf())
);
surfaceScalarField phiDragb
(
fvc::interpolate(alpha/rhob*dragCoef*rUbA)*phia + rUbAf*(g & mesh.Sf())
);
forAll(p.boundaryField(), patchi)
{
@ -50,7 +54,7 @@
if (nonOrth == nNonOrthCorr)
{
surfaceScalarField SfGradp = pEqn.flux()/Dp;
surfaceScalarField SfGradp(pEqn.flux()/Dp);
phia -= rUaAf*SfGradp/rhoa;
phib -= rUbAf*SfGradp/rhob;

2
applications/solvers/multiphase/bubbleFoam/wallFunctions.H
View File

@ -34,7 +34,7 @@
{
const scalarField& nuw = nutb.boundaryField()[patchi];
scalarField magFaceGradU = mag(U.boundaryField()[patchi].snGrad());
scalarField magFaceGradU(mag(U.boundaryField()[patchi].snGrad()));
forAll(currPatch, facei)
{

6
applications/solvers/multiphase/cavitatingFoam/CourantNo.H
View File

@ -35,8 +35,10 @@ scalar acousticCoNum = 0.0;
if (mesh.nInternalFaces())
{
scalarField sumPhi =
fvc::surfaceSum(mag(phiv))().internalField();
scalarField sumPhi
(
fvc::surfaceSum(mag(phiv))().internalField()
);
CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

8
applications/solvers/multiphase/cavitatingFoam/pEqn.H
View File

@ -9,18 +9,18 @@
)/psi;
}
surfaceScalarField rhof = fvc::interpolate(rho, "rhof");
surfaceScalarField rhof(fvc::interpolate(rho, "rhof"));
volScalarField rAU = 1.0/UEqn.A();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU));
volVectorField HbyA = rAU*UEqn.H();
volVectorField HbyA(rAU*UEqn.H());
phiv = (fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phiv);
p.boundaryField().updateCoeffs();
surfaceScalarField phiGradp = rAUf*mesh.magSf()*fvc::snGrad(p);
surfaceScalarField phiGradp(rAUf*mesh.magSf()*fvc::snGrad(p));
phiv -= phiGradp/rhof;

12
applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H
View File

@ -2,7 +2,7 @@
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir = phic*interface.nHatf();
surfaceScalarField phir(phic*interface.nHatf());
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
@ -45,7 +45,8 @@
}
surfaceScalarField phiAlpha1 =
surfaceScalarField phiAlpha1
(
fvc::flux
(
phi,
@ -57,7 +58,8 @@
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
);
)
);
MULES::explicitSolve
(
@ -71,8 +73,8 @@
0
);
surfaceScalarField rho1f = fvc::interpolate(rho1);
surfaceScalarField rho2f = fvc::interpolate(rho2);
surfaceScalarField rho1f(fvc::interpolate(rho1));
surfaceScalarField rho2f(fvc::interpolate(rho2));
rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
alpha2 = scalar(1) - alpha1;

6
applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
View File

@ -9,15 +9,15 @@
readLabel(piso.lookup("nAlphaSubCycles"))
);
surfaceScalarField phic = mag(phi/mesh.magSf());
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
volScalarField divU = fvc::div(phi);
volScalarField divU(fvc::div(phi));
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
surfaceScalarField rhoPhiSum(0.0*rhoPhi);
for
(

6
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H
View File

@ -9,17 +9,17 @@
readLabel(piso.lookup("nAlphaSubCycles"))
);
surfaceScalarField phic = mag(phi/mesh.magSf());
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
fvc::makeAbsolute(phi, U);
volScalarField divU = fvc::div(phi);
volScalarField divU(fvc::div(phi));
fvc::makeRelative(phi, U);
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
surfaceScalarField rhoPhiSum(0.0*rhoPhi);
for
(

2
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
View File

@ -79,7 +79,7 @@ int main(int argc, char *argv[])
{
// Store divU from the previous mesh for the correctPhi
volScalarField divU = fvc::div(phi);
volScalarField divU(fvc::div(phi));
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();

4
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
View File

@ -1,6 +1,6 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf(fvc::interpolate(rAU));
tmp<fvScalarMatrix> p_rghEqnComp;

10
applications/solvers/multiphase/compressibleInterFoam/createFields.H
View File

@ -105,8 +105,8 @@
)
);
volScalarField rho1 = rho10 + psi1*p;
volScalarField rho2 = rho20 + psi2*p;
volScalarField rho1(rho10 + psi1*p);
volScalarField rho2(rho20 + psi2*p);
volScalarField rho
(
@ -138,8 +138,10 @@
fvc::interpolate(rho)*phi
);
volScalarField dgdt =
pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001));
volScalarField dgdt
(
pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
);
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);

4
applications/solvers/multiphase/compressibleInterFoam/pEqn.H
View File

@ -1,6 +1,6 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf(fvc::interpolate(rAU));
tmp<fvScalarMatrix> p_rghEqnComp;

12
applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
View File

@ -56,7 +56,7 @@ void Foam::MULES::explicitLTSSolve
const fvMesh& mesh = psi.mesh();
psi.correctBoundaryConditions();
surfaceScalarField phiBD = upwind<scalar>(psi.mesh(), phi).flux(psi);
surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
surfaceScalarField& phiCorr = phiPsi;
phiCorr -= phiBD;
@ -168,9 +168,11 @@ void Foam::MULES::implicitSolve
scalarField allCoLambda(mesh.nFaces());
{
surfaceScalarField Cof =
surfaceScalarField Cof
(
mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
*mag(phi)/mesh.magSf();
*mag(phi)/mesh.magSf()
);
slicedSurfaceScalarField CoLambda
(
@ -226,7 +228,7 @@ void Foam::MULES::implicitSolve
- Su
);
surfaceScalarField phiBD = psiConvectionDiffusion.flux();
surfaceScalarField phiBD(psiConvectionDiffusion.flux());
surfaceScalarField& phiCorr = phiPsi;
phiCorr -= phiBD;
@ -408,7 +410,7 @@ void Foam::MULES::limiter
if (psiPf.coupled())
{
scalarField psiPNf = psiPf.patchNeighbourField();
scalarField psiPNf(psiPf.patchNeighbourField());
forAll(phiCorrPf, pFacei)
{

10
applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
View File

@ -2,13 +2,14 @@
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phic = mag(phi/mesh.magSf());
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
surfaceScalarField phir = phic*interface.nHatf();
surfaceScalarField phir(phic*interface.nHatf());
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
surfaceScalarField phiAlpha =
surfaceScalarField phiAlpha
(
fvc::flux
(
phi,
@ -20,7 +21,8 @@
-fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
alpha1,
alpharScheme
);
)
);
MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
//MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);

2
applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
View File

@ -11,7 +11,7 @@ label nAlphaSubCycles
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
surfaceScalarField rhoPhiSum(0.0*rhoPhi);
for
(

4
applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
View File

@ -1,6 +1,6 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf(fvc::interpolate(rAU));
U = rAU*UEqn.H();
surfaceScalarField phiU

6
applications/solvers/multiphase/interFoam/alphaCourantNo.H
View File

@ -39,9 +39,11 @@ scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
scalarField sumPhi =
scalarField sumPhi
(
pos(alpha1 - 0.01)*pos(0.99 - alpha1)
*fvc::surfaceSum(mag(phi))().internalField();
*fvc::surfaceSum(mag(phi))().internalField()
);
alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

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

@ -2,13 +2,14 @@
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phic = mag(phi/mesh.magSf());
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
surfaceScalarField phir = phic*interface.nHatf();
surfaceScalarField phir(phic*interface.nHatf());
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
surfaceScalarField phiAlpha =
surfaceScalarField phiAlpha
(
fvc::flux
(
phi,
@ -20,7 +21,8 @@
-fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
alpha1,
alpharScheme
);
)
);
MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);

2
applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
View File

@ -11,7 +11,7 @@ label nAlphaSubCycles
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
surfaceScalarField rhoPhiSum(0.0*rhoPhi);
for
(

4
applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
View File

@ -1,6 +1,6 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf(fvc::interpolate(rAU));
U = rAU*UEqn.H();
surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf()));

4
applications/solvers/multiphase/interFoam/pEqn.H
View File

@ -1,6 +1,6 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf(fvc::interpolate(rAU));
U = rAU*UEqn.H();
surfaceScalarField phiU

11
applications/solvers/multiphase/interMixingFoam/alphaCourantNo.H
View File

@ -39,11 +39,14 @@ scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
scalarField sumPhi = max
scalarField sumPhi
(
pos(alpha1 - 0.01)*pos(0.99 - alpha1),
pos(alpha2 - 0.01)*pos(0.99 - alpha2)
)*fvc::surfaceSum(mag(phi))().internalField();
max
(
pos(alpha1 - 0.01)*pos(0.99 - alpha1),
pos(alpha2 - 0.01)*pos(0.99 - alpha2)
)*fvc::surfaceSum(mag(phi))().internalField()
);
alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

28
applications/solvers/multiphase/interMixingFoam/alphaEqns.H
View File

@ -40,7 +40,8 @@
// Create the complete convection flux for alpha1
surfaceScalarField phiAlpha1 =
surfaceScalarField phiAlpha1
(
fvc::flux
(
phi,
@ -58,11 +59,14 @@
-fvc::flux(-phir, alpha3, alpharScheme),
alpha1,
alpharScheme
);
)
);
// Create the bounded (upwind) flux for alpha1
surfaceScalarField phiAlpha1BD =
upwind<scalar>(mesh, phi).flux(alpha1);
surfaceScalarField phiAlpha1BD
(
upwind<scalar>(mesh, phi).flux(alpha1)
);
// Calculate the flux correction for alpha1
phiAlpha1 -= phiAlpha1BD;
@ -83,7 +87,8 @@
);
// Create the complete flux for alpha2
surfaceScalarField phiAlpha2 =
surfaceScalarField phiAlpha2
(
fvc::flux
(
phi,
@ -95,11 +100,14 @@
-fvc::flux(phir, alpha1, alpharScheme),
alpha2,
alpharScheme
);
)
);
// Create the bounded (upwind) flux for alpha2
surfaceScalarField phiAlpha2BD =
upwind<scalar>(mesh, phi).flux(alpha2);
surfaceScalarField phiAlpha2BD
(
upwind<scalar>(mesh, phi).flux(alpha2)
);
// Calculate the flux correction for alpha2
phiAlpha2 -= phiAlpha2BD;
@ -127,8 +135,8 @@
solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
// Create the diffusion coefficients for alpha2<->alpha3
volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2);
volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3);
volScalarField Dc23(D23*max(alpha3, scalar(0))*pos(alpha2));
volScalarField Dc32(D23*max(alpha2, scalar(0))*pos(alpha3));
// Add the diffusive flux for alpha3->alpha2
phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);

4
applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H
View File

@ -10,7 +10,7 @@ label nAlphaSubCycles
if (nAlphaSubCycles > 1)
{
surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
surfaceScalarField rhoPhiSum(0.0*rhoPhi);
dimensionedScalar totalDeltaT = runTime.deltaT();
for
@ -33,7 +33,7 @@ else
interface.correct();
{
volScalarField rhoNew = alpha1*rho1 + alpha2*rho2 + alpha3*rho3;
volScalarField rhoNew(alpha1*rho1 + alpha2*rho2 + alpha3*rho3);
//solve(fvm::ddt(rho) + fvc::div(rhoPhi));
//Info<< "density error = "

12
applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C
View File

@ -133,9 +133,9 @@ Foam::tmp<Foam::volScalarField> Foam::threePhaseMixture::mu() const
Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::muf() const
{
surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
surfaceScalarField alpha1f(fvc::interpolate(alpha1_));
surfaceScalarField alpha2f(fvc::interpolate(alpha2_));
surfaceScalarField alpha3f(fvc::interpolate(alpha3_));
return tmp<surfaceScalarField>
(
@ -152,9 +152,9 @@ Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::muf() const
Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::nuf() const
{
surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
surfaceScalarField alpha1f(fvc::interpolate(alpha1_));
surfaceScalarField alpha2f(fvc::interpolate(alpha2_));
surfaceScalarField alpha3f(fvc::interpolate(alpha3_));
return tmp<surfaceScalarField>
(

36
applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
View File

@ -81,31 +81,35 @@ void Foam::threePhaseInterfaceProperties::correctContactAngle
refCast<const alphaContactAngleFvPatchScalarField>
(alpha3[patchi]);
scalarField twoPhaseAlpha2 = max(a2cap, scalar(0));
scalarField twoPhaseAlpha3 = max(a3cap, scalar(0));
scalarField twoPhaseAlpha2(max(a2cap, scalar(0)));
scalarField twoPhaseAlpha3(max(a3cap, scalar(0)));
scalarField sumTwoPhaseAlpha =
twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL;
scalarField sumTwoPhaseAlpha
(
twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL
);
twoPhaseAlpha2 /= sumTwoPhaseAlpha;
twoPhaseAlpha3 /= sumTwoPhaseAlpha;
fvsPatchVectorField& nHatp = nHatb[patchi];
scalarField theta =
scalarField theta
(
convertToRad
*(
* (
twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
+ twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
);
)
);
vectorField nf = boundary[patchi].nf();
vectorField nf(boundary[patchi].nf());
// Reset nHatPatch to correspond to the contact angle
scalarField a12 = nHatp & nf;
scalarField a12(nHatp & nf);
scalarField b1 = cos(theta);
scalarField b1(cos(theta));
scalarField b2(nHatp.size());
@ -114,10 +118,10 @@ void Foam::threePhaseInterfaceProperties::correctContactAngle
b2[facei] = cos(acos(a12[facei]) - theta[facei]);
}
scalarField det = 1.0 - a12*a12;
scalarField det(1.0 - a12*a12);
scalarField a = (b1 - a12*b2)/det;
scalarField b = (b2 - a12*b1)/det;
scalarField a((b1 - a12*b2)/det);
scalarField b((b2 - a12*b1)/det);
nHatp = a*nf + b*nHatp;
@ -135,13 +139,13 @@ void Foam::threePhaseInterfaceProperties::calculateK()
const surfaceVectorField& Sf = mesh.Sf();
// Cell gradient of alpha
volVectorField gradAlpha = fvc::grad(alpha1);
volVectorField gradAlpha(fvc::grad(alpha1));
// Interpolated face-gradient of alpha
surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
// Face unit interface normal
surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
correctContactAngle(nHatfv.boundaryField());
// Face unit interface normal flux

4
applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H
View File

@ -125,8 +125,8 @@ public:
tmp<volScalarField> sigma() const
{
volScalarField limitedAlpha2 = max(mixture_.alpha2(), scalar(0));
volScalarField limitedAlpha3 = max(mixture_.alpha3(), scalar(0));
volScalarField limitedAlpha2(max(mixture_.alpha2(), scalar(0)));
volScalarField limitedAlpha3(max(mixture_.alpha3(), scalar(0)));
return
(limitedAlpha2*sigma12_ + limitedAlpha3*sigma13_)

6
applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
View File

@ -6,7 +6,8 @@
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
surfaceScalarField phiAlpha =
surfaceScalarField phiAlpha
(
fvc::flux
(
phi,
@ -18,7 +19,8 @@
-fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
alpha1,
alpharScheme
);
)
);
Pair<tmp<volScalarField> > vDotAlphal =
twoPhaseProperties->vDotAlphal();

4
applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
View File

@ -21,10 +21,10 @@ surfaceScalarField rhoPhi
readLabel(piso.lookup("nAlphaSubCycles"))
);
surfaceScalarField phic = mag(phi/mesh.magSf());
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
volScalarField divU = fvc::div(phi);
volScalarField divU(fvc::div(phi));
dimensionedScalar totalDeltaT = runTime.deltaT();

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

@ -1,6 +1,6 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf(fvc::interpolate(rAU));
U = rAU*UEqn.H();

4
applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C
View File

@ -68,7 +68,7 @@ Foam::Pair<Foam::tmp<Foam::volScalarField> >
Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotAlphal() const
{
const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
volScalarField limitedAlpha1 = min(max(alpha1_, scalar(0)), scalar(1));
volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
return Pair<tmp<volScalarField> >
(
@ -83,7 +83,7 @@ Foam::Pair<Foam::tmp<Foam::volScalarField> >
Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotP() const
{
const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
volScalarField limitedAlpha1 = min(max(alpha1_, scalar(0)), scalar(1));
volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
return Pair<tmp<volScalarField> >
(

2
applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C
View File

@ -80,7 +80,7 @@ Foam::Pair<Foam::tmp<Foam::volScalarField> >
Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotP() const
{
const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
volScalarField limitedAlpha1 = min(max(alpha1_, scalar(0)), scalar(1));
volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
return Pair<tmp<volScalarField> >
(

16
applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C
View File

@ -97,9 +97,11 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::pCoeff
const volScalarField& p
) const
{
volScalarField limitedAlpha1 = min(max(alpha1_, scalar(0)), scalar(1));
volScalarField rho =
(limitedAlpha1*rho1() + (scalar(1) - limitedAlpha1)*rho2());
volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
volScalarField rho
(
limitedAlpha1*rho1() + (scalar(1) - limitedAlpha1)*rho2()
);
return
(3*rho1()*rho2())*sqrt(2/(3*rho1()))
@ -111,9 +113,9 @@ Foam::Pair<Foam::tmp<Foam::volScalarField> >
Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotAlphal() const
{
const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
volScalarField limitedAlpha1 = min(max(alpha1_, scalar(0)), scalar(1));
volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
volScalarField pCoeff = this->pCoeff(p);
volScalarField pCoeff(this->pCoeff(p));
return Pair<tmp<volScalarField> >
(
@ -128,9 +130,9 @@ Foam::Pair<Foam::tmp<Foam::volScalarField> >
Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotP() const
{
const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
volScalarField limitedAlpha1 = min(max(alpha1_, scalar(0)), scalar(1));
volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
volScalarField apCoeff = limitedAlpha1*pCoeff(p);
volScalarField apCoeff(limitedAlpha1*pCoeff(p));
return Pair<tmp<volScalarField> >
(

2
applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
View File

@ -54,7 +54,7 @@ Foam::phaseChangeTwoPhaseMixture::phaseChangeTwoPhaseMixture
Foam::Pair<Foam::tmp<Foam::volScalarField> >
Foam::phaseChangeTwoPhaseMixture::vDotAlphal() const
{
volScalarField alphalCoeff = 1.0/rho1() - alpha1_*(1.0/rho1() - 1.0/rho2());
volScalarField alphalCoeff(1.0/rho1() - alpha1_*(1.0/rho1() - 1.0/rho2()));
Pair<tmp<volScalarField> > mDotAlphal = this->mDotAlphal();
return Pair<tmp<volScalarField> >

4
applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
View File

@ -1,6 +1,6 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf(fvc::interpolate(rAU));
U = rAU*UEqn.H();

6
applications/solvers/multiphase/multiphaseInterFoam/alphaCourantNo.H
View File

@ -39,9 +39,11 @@ scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
scalarField sumPhi =
scalarField sumPhi
(
mixture.nearInterface()().internalField()
*fvc::surfaceSum(mag(phi))().internalField();
* fvc::surfaceSum(mag(phi))().internalField()
);
alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

6
applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C
View File

@ -133,7 +133,11 @@ void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField(fvPatchScalarField, alphaContactAngleFvPatchScalarField);
makeNonTemplatedPatchTypeField
(
fvPatchScalarField,
alphaContactAngleFvPatchScalarField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

50
applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
View File

@ -273,7 +273,7 @@ void Foam::multiphaseMixture::solve()
if (nAlphaSubCycles > 1)
{
surfaceScalarField rhoPhiSum = 0.0*rhoPhi_;
surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
dimensionedScalar totalDeltaT = runTime.deltaT();
for
@ -314,9 +314,11 @@ Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
*/
surfaceVectorField gradAlphaf =
surfaceVectorField gradAlphaf
(
fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
- fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2));
- fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
);
// Face unit interface normal
return gradAlphaf/(mag(gradAlphaf) + deltaN_);
@ -361,9 +363,11 @@ void Foam::multiphaseMixture::correctContactAngle
vectorField& nHatPatch = nHatb[patchi];
vectorField AfHatPatch =
vectorField AfHatPatch
(
mesh_.Sf().boundaryField()[patchi]
/mesh_.magSf().boundaryField()[patchi];
/mesh_.magSf().boundaryField()[patchi]
);
alphaContactAngleFvPatchScalarField::thetaPropsTable::
const_iterator tp =
@ -396,21 +400,25 @@ void Foam::multiphaseMixture::correctContactAngle
scalar thetaR = convertToRad*tp().thetaR(matched);
// Calculated the component of the velocity parallel to the wall
vectorField Uwall =
vectorField Uwall
(
U_.boundaryField()[patchi].patchInternalField()
- U_.boundaryField()[patchi];
- U_.boundaryField()[patchi]
);
Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
// Find the direction of the interface parallel to the wall
vectorField nWall =
nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch;
vectorField nWall
(
nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
);
// Normalise nWall
nWall /= (mag(nWall) + SMALL);
// Calculate Uwall resolved normal to the interface parallel to
// the interface
scalarField uwall = nWall & Uwall;
scalarField uwall(nWall & Uwall);
theta += (thetaA - thetaR)*tanh(uwall/uTheta);
}
@ -418,9 +426,9 @@ void Foam::multiphaseMixture::correctContactAngle
// Reset nHatPatch to correspond to the contact angle
scalarField a12 = nHatPatch & AfHatPatch;
scalarField a12(nHatPatch & AfHatPatch);
scalarField b1 = cos(theta);
scalarField b1(cos(theta));
scalarField b2(nHatPatch.size());
@ -429,10 +437,10 @@ void Foam::multiphaseMixture::correctContactAngle
b2[facei] = cos(acos(a12[facei]) - theta[facei]);
}
scalarField det = 1.0 - a12*a12;
scalarField det(1.0 - a12*a12);
scalarField a = (b1 - a12*b2)/det;
scalarField b = (b2 - a12*b1)/det;
scalarField a((b1 - a12*b2)/det);
scalarField b((b2 - a12*b1)/det);
nHatPatch = a*AfHatPatch + b*nHatPatch;
@ -508,7 +516,7 @@ void Foam::multiphaseMixture::solveAlphas
)
);
surfaceScalarField phic = mag(phi_/mesh_.magSf());
surfaceScalarField phic(mag(phi_/mesh_.magSf()));
phic = min(cAlpha*phic, max(phic));
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
@ -527,7 +535,7 @@ void Foam::multiphaseMixture::solveAlphas
phase& refPhase = *refPhasePtr;
volScalarField refPhaseNew = refPhase;
volScalarField refPhaseNew(refPhase);
refPhaseNew == 1.0;
rhoPhi_ = phi_*refPhase.rho();
@ -550,9 +558,11 @@ void Foam::multiphaseMixture::solveAlphas
if (&alpha2 == &alpha) continue;
surfaceScalarField phir = phic*nHatf(alpha, alpha2);
surfaceScalarField phirb12 =
-fvc::flux(-phir, alpha2, alphacScheme);
surfaceScalarField phir(phic*nHatf(alpha, alpha2));
surfaceScalarField phirb12
(
-fvc::flux(-phir, alpha2, alphacScheme)
);
alphaEqn += fvm::div(phirb12, alpha, alphacScheme);
}

4
applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
View File

@ -1,6 +1,6 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf(fvc::interpolate(rAU));
U = rAU*UEqn.H();

2
applications/solvers/multiphase/settlingFoam/createFields.H
View File

@ -345,7 +345,7 @@
Info<< "Calculating field (g.h)f\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf = surfaceScalarField("gh", g & mesh.Cf());
surfaceScalarField ghf(surfaceScalarField("gh", g & mesh.Cf()));
volScalarField p
(

10
applications/solvers/multiphase/settlingFoam/kEpsilon.H
View File

@ -8,14 +8,16 @@ if (turbulence)
dimensionedScalar kMin("kMin", k.dimensions(), SMALL);
dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
volScalarField divU = fvc::div(phi/fvc::interpolate(rho));
volScalarField divU(fvc::div(phi/fvc::interpolate(rho)));
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField G = 2*mut*(tgradU() && dev(symm(tgradU())));
volScalarField G(2*mut*(tgradU() && dev(symm(tgradU()))));
tgradU.clear();
volScalarField Gcoef =
Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin);
volScalarField Gcoef
(
Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin)
);
#include "wallFunctions.H"

2
applications/solvers/multiphase/settlingFoam/pEqn.H
View File

@ -1,4 +1,4 @@
volScalarField rAU = 1.0/UEqn.A();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf
(

4
applications/solvers/multiphase/settlingFoam/plasticViscosity.H
View File

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

8
applications/solvers/multiphase/settlingFoam/wallFunctions.H
View File

@ -34,11 +34,13 @@
{
const scalarField& rhow = rho.boundaryField()[patchi];
const scalarField muw = mul.boundaryField()[patchi];
const scalarField muw(mul.boundaryField()[patchi]);
const scalarField& mutw = mut.boundaryField()[patchi];
scalarField magFaceGradU =
mag(U.boundaryField()[patchi].snGrad());
scalarField magFaceGradU
(
mag(U.boundaryField()[patchi].snGrad())
);
forAll(curPatch, facei)
{

2
applications/solvers/multiphase/settlingFoam/wallViscosity.H
View File

@ -13,7 +13,7 @@
{
const scalarField& rhow = rho.boundaryField()[patchi];
const scalarField muw = mul.boundaryField()[patchi];
const scalarField muw(mul.boundaryField()[patchi]);
scalarField& mutw = mut.boundaryField()[patchi];
forAll(curPatch, facei)

4
applications/solvers/multiphase/settlingFoam/yieldStress.H
View File

@ -6,7 +6,7 @@ volScalarField yieldStress
const volScalarField& alpha
)
{
return
tmp<volScalarField> tfld
(
yieldStressCoeff*
(
@ -14,4 +14,6 @@ volScalarField yieldStress
- pow(10.0, yieldStressExponent*yieldStressOffset)
)
);
return tfld();
}

4
applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
View File

@ -1,6 +1,6 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf(fvc::interpolate(rAU));
U = rAU*UEqn.H();
surfaceScalarField phiU

16
applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
View File

@ -3,7 +3,7 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
{
{
volTensorField gradUaT = fvc::grad(Ua)().T();
volTensorField gradUaT(fvc::grad(Ua)().T());
if (kineticTheory.on())
{
@ -26,9 +26,11 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);
}
surfaceScalarField phiRa =
surfaceScalarField phiRa
(
-fvc::interpolate(nuEffa)*mesh.magSf()*fvc::snGrad(alpha)
/fvc::interpolate(alpha + scalar(0.001));
/fvc::interpolate(alpha + scalar(0.001))
);
UaEqn =
(
@ -56,16 +58,18 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
}
{
volTensorField gradUbT = fvc::grad(Ub)().T();
volTensorField gradUbT(fvc::grad(Ub)().T());
volTensorField Rcb
(
"Rcb",
((2.0/3.0)*I)*(k + nuEffb*tr(gradUbT)) - nuEffb*gradUbT
);
surfaceScalarField phiRb =
surfaceScalarField phiRb
(
-fvc::interpolate(nuEffb)*mesh.magSf()*fvc::snGrad(beta)
/fvc::interpolate(beta + scalar(0.001));
/fvc::interpolate(beta + scalar(0.001))
);
UbEqn =
(

8
applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
View File

@ -2,13 +2,13 @@
word scheme("div(phi,alpha)");
word schemer("div(phir,alpha)");
surfaceScalarField phic = phi;
surfaceScalarField phir = phia - phib;
surfaceScalarField phic(phi);
surfaceScalarField phir(phia - phib);
if (g0.value() > 0.0)
{
surfaceScalarField alphaf = fvc::interpolate(alpha);
surfaceScalarField phipp = ppMagf*fvc::snGrad(alpha)*mesh.magSf();
surfaceScalarField alphaf(fvc::interpolate(alpha));
surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha)*mesh.magSf());
phir += phipp;
phic += fvc::interpolate(alpha)*phipp;
}

12
applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
View File

@ -132,15 +132,19 @@
Info<< "Calculating field DDtUa and DDtUb\n" << endl;
volVectorField DDtUa =
volVectorField DDtUa
(
fvc::ddt(Ua)
+ fvc::div(phia, Ua)
- fvc::div(phia)*Ua;
- fvc::div(phia)*Ua
);
volVectorField DDtUb =
volVectorField DDtUb
(
fvc::ddt(Ub)
+ fvc::div(phib, Ub)
- fvc::div(phib)*Ub;
- fvc::div(phib)*Ub
);
Info<< "Calculating field g.h\n" << endl;

2
applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/Ergun/Ergun.C
View File

@ -68,7 +68,7 @@ Foam::tmp<Foam::volScalarField> Foam::Ergun::K
const volScalarField& Ur
) const
{
volScalarField beta = max(scalar(1) - alpha_, scalar(1.0e-6));
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
return
150.0*alpha_*phaseb_.nu()*phaseb_.rho()

6
applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/Gibilaro/Gibilaro.C
View File

@ -68,9 +68,9 @@ Foam::tmp<Foam::volScalarField> Foam::Gibilaro::K
const volScalarField& Ur
) const
{
volScalarField beta = max(scalar(1) - alpha_, scalar(1.0e-6));
volScalarField bp = pow(beta, -2.8);
volScalarField Re = max(beta*Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
volScalarField bp(pow(beta, -2.8));
volScalarField Re(max(beta*Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
return (17.3/Re + scalar(0.336))*phaseb_.rho()*Ur*bp/phasea_.d();
}

8
applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
View File

@ -68,12 +68,12 @@ Foam::tmp<Foam::volScalarField> Foam::GidaspowErgunWenYu::K
const volScalarField& Ur
) const
{
volScalarField beta = max(scalar(1) - alpha_, scalar(1.0e-6));
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
volScalarField bp = pow(beta, -2.65);
volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
volScalarField Cds = 24.0*(1.0 + 0.15*pow(Re, 0.687))/Re;
volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{

8
applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
View File

@ -68,11 +68,11 @@ Foam::tmp<Foam::volScalarField> Foam::GidaspowSchillerNaumann::K
const volScalarField& Ur
) const
{
volScalarField beta = max(scalar(1) - alpha_, scalar(1e-6));
volScalarField bp = pow(beta, -2.65);
volScalarField beta(max(scalar(1) - alpha_, scalar(1e-6)));
volScalarField bp(pow(beta, -2.65));
volScalarField Re = max(beta*Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
volScalarField Cds = 24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re;
volScalarField Re(max(beta*Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{

4
applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
View File

@ -68,8 +68,8 @@ Foam::tmp<Foam::volScalarField> Foam::SchillerNaumann::K
const volScalarField& Ur
) const
{
volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
volScalarField Cds = 24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re;
volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{

17
applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
View File

@ -68,9 +68,9 @@ Foam::tmp<Foam::volScalarField> Foam::SyamlalOBrien::K
const volScalarField& Ur
) const
{
volScalarField beta = max(scalar(1) - alpha_, scalar(1.0e-6));
volScalarField A = pow(beta, 4.14);
volScalarField B = 0.8*pow(beta, 1.28);
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
volScalarField A(pow(beta, 4.14));
volScalarField B(0.8*pow(beta, 1.28));
forAll (beta, celli)
{
@ -80,14 +80,17 @@ Foam::tmp<Foam::volScalarField> Foam::SyamlalOBrien::K
}
}
volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
volScalarField Vr = 0.5*
volScalarField Vr
(
A - 0.06*Re + sqrt(sqr(0.06*Re) + 0.12*Re*(2.0*B - A) + sqr(A))
0.5*
(
A - 0.06*Re + sqrt(sqr(0.06*Re) + 0.12*Re*(2.0*B - A) + sqr(A))
)
);
volScalarField Cds = sqr(0.63 + 4.8*sqrt(Vr/Re));
volScalarField Cds(sqr(0.63 + 4.8*sqrt(Vr/Re)));
return 0.75*Cds*phaseb_.rho()*Ur/(phasea_.d()*sqr(Vr));
}

8
applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
View File

@ -68,11 +68,11 @@ Foam::tmp<Foam::volScalarField> Foam::WenYu::K
const volScalarField& Ur
) const
{
volScalarField beta = max(scalar(1) - alpha_, scalar(1.0e-6));
volScalarField bp = pow(beta, -2.65);
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
volScalarField bp(pow(beta, -2.65));
volScalarField Re = max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3));
volScalarField Cds = 24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re;
volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{

2
applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H
View File

@ -6,7 +6,7 @@ if (turbulence)
}
tmp<volTensorField> tgradUb = fvc::grad(Ub);
volScalarField G = 2*nutb*(tgradUb() && dev(symm(tgradUb())));
volScalarField G(2*nutb*(tgradUb() && dev(symm(tgradUb()))));
tgradUb.clear();
#include "wallFunctions.H"

6
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C
View File

@ -75,8 +75,10 @@ Foam::tmp<Foam::volScalarField> Foam::HrenyaSinclairConductivity::kappa
{
const scalar sqrtPi = sqrt(constant::mathematical::pi);
volScalarField lamda =
scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_;
volScalarField lamda
(
scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_
);
return rhoa*da*sqrt(Theta)*
(

108
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
View File

@ -202,16 +202,16 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
const scalar sqrtPi = sqrt(constant::mathematical::pi);
surfaceScalarField phi = 1.5*rhoa_*phia_*fvc::interpolate(alpha_);
surfaceScalarField phi(1.5*rhoa_*phia_*fvc::interpolate(alpha_));
volTensorField dU = gradUat.T();//fvc::grad(Ua_);
volSymmTensorField D = symm(dU);
volTensorField dU(gradUat.T()); //fvc::grad(Ua_);
volSymmTensorField D(symm(dU));
// NB, drag = K*alpha*beta,
// (the alpha and beta has been extracted from the drag function for
// numerical reasons)
volScalarField Ur = mag(Ua_ - Ub_);
volScalarField betaPrim = alpha_*(1.0 - alpha_)*draga_.K(Ur);
volScalarField Ur(mag(Ua_ - Ub_));
volScalarField betaPrim(alpha_*(1.0 - alpha_)*draga_.K(Ur));
// Calculating the radial distribution function (solid volume fraction is
// limited close to the packing limit, but this needs improvements)
@ -223,12 +223,15 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
);
// particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
volScalarField PsCoeff = granularPressureModel_->granularPressureCoeff
volScalarField PsCoeff
(
alpha_,
gs0_,
rhoa_,
e_
granularPressureModel_->granularPressureCoeff
(
alpha_,
gs0_,
rhoa_,
e_
)
);
// 'thermal' conductivity (Table 3.3, p. 49)
@ -245,23 +248,27 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
);
dimensionedScalar TsmallSqrt = sqrt(Tsmall);
volScalarField ThetaSqrt = sqrt(Theta_);
volScalarField ThetaSqrt(sqrt(Theta_));
// dissipation (Eq. 3.24, p.50)
volScalarField gammaCoeff =
12.0*(1.0 - sqr(e_))*sqr(alpha_)*rhoa_*gs0_*(1.0/da_)*ThetaSqrt/sqrtPi;
volScalarField gammaCoeff
(
12.0*(1.0 - sqr(e_))*sqr(alpha_)*rhoa_*gs0_*(1.0/da_)*ThetaSqrt/sqrtPi
);
// Eq. 3.25, p. 50 Js = J1 - J2
volScalarField J1 = 3.0*betaPrim;
volScalarField J2 =
volScalarField J1(3.0*betaPrim);
volScalarField J2
(
0.25*sqr(betaPrim)*da_*sqr(Ur)
/(max(alpha_, 1e-6)*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt));
/ (max(alpha_, 1e-6)*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt))
);
// bulk viscosity p. 45 (Lun et al. 1984).
lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi;
// 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*mua_*D + (lambda_ - (2.0/3.0)*mua_)*tr(D)*I);
if (!equilibrium_)
{
@ -289,27 +296,32 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
{
// equilibrium => dissipation == production
// Eq. 4.14, p.82
volScalarField K1 = 2.0*(1.0 + e_)*rhoa_*gs0_;
volScalarField K3 = 0.5*da_*rhoa_*
volScalarField K1(2.0*(1.0 + e_)*rhoa_*gs0_);
volScalarField K3
(
0.5*da_*rhoa_*
(
(sqrtPi/(3.0*(3.0-e_)))
*(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha_*gs0_)
+1.6*alpha_*gs0_*(1.0 + e_)/sqrtPi
);
)
);
volScalarField K2 =
4.0*da_*rhoa_*(1.0 + e_)*alpha_*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0;
volScalarField K2
(
4.0*da_*rhoa_*(1.0 + e_)*alpha_*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_))*rhoa_*gs0_/(da_*sqrtPi));
volScalarField trD = tr(D);
volScalarField tr2D = sqr(trD);
volScalarField trD2 = tr(D & D);
volScalarField trD(tr(D));
volScalarField tr2D(sqr(trD));
volScalarField trD2(tr(D & D));
volScalarField t1 = K1*alpha_ + rhoa_;
volScalarField l1 = -t1*trD;
volScalarField l2 = sqr(t1)*tr2D;
volScalarField l3 = 4.0*K4*max(alpha_, 1e-6)*(2.0*K3*trD2 + K2*tr2D);
volScalarField t1(K1*alpha_ + rhoa_);
volScalarField l1(-t1*trD);
volScalarField l2(sqr(t1)*tr2D);
volScalarField l3(4.0*K4*max(alpha_, 1e-6)*(2.0*K3*trD2 + K2*tr2D));
Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha_ + 1.0e-4)*K4));
}
@ -317,14 +329,17 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
Theta_.max(1.0e-15);
Theta_.min(1.0e+3);
volScalarField pf = frictionalStressModel_->frictionalPressure
volScalarField pf
(
alpha_,
alphaMinFriction_,
alphaMax_,
Fr_,
eta_,
p_
frictionalStressModel_->frictionalPressure
(
alpha_,
alphaMinFriction_,
alphaMax_,
Fr_,
eta_,
p_
)
);
PsCoeff += pf/(Theta_+Tsmall);
@ -336,23 +351,26 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
pa_ = PsCoeff*Theta_;
// frictional shear stress, Eq. 3.30, p. 52
volScalarField muf = frictionalStressModel_->muf
volScalarField muf
(
alpha_,
alphaMax_,
pf,
D,
phi_
frictionalStressModel_->muf
(
alpha_,
alphaMax_,
pf,
D,
phi_
)
);
// add frictional stress
// add frictional stress
mua_ += muf;
mua_.min(1.0e+2);
mua_.max(0.0);
Info<< "kinTheory: max(Theta) = " << max(Theta_).value() << endl;
volScalarField ktn = mua_/rhoa_;
volScalarField ktn(mua_/rhoa_);
Info<< "kinTheory: min(nua) = " << min(ktn).value()
<< ", max(nua) = " << max(ktn).value() << endl;

6
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.C
View File

@ -78,8 +78,10 @@ Foam::kineticTheoryModels::HrenyaSinclairViscosity::mua
{
const scalar sqrtPi = sqrt(constant::mathematical::pi);
volScalarField lamda =
scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_;
volScalarField lamda
(
scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_
);
return rhoa*da*sqrt(Theta)*
(

14
applications/solvers/multiphase/twoPhaseEulerFoam/liftDragCoeffs.H
View File

@ -1,18 +1,18 @@
volVectorField Ur = Ua - Ub;
volScalarField magUr = mag(Ur);
volVectorField Ur(Ua - Ub);
volScalarField magUr(mag(Ur));
volScalarField Ka = draga->K(magUr);
volScalarField K = Ka;
volScalarField Ka(draga->K(magUr));
volScalarField K(Ka);
if (dragPhase == "b")
{
volScalarField Kb = dragb->K(magUr);
volScalarField Kb(dragb->K(magUr));
K = Kb;
}
else if (dragPhase == "blended")
{
volScalarField Kb = dragb->K(magUr);
volScalarField Kb(dragb->K(magUr));
K = (beta*Ka + alpha*Kb);
}
volVectorField liftCoeff = Cl*(beta*rhob + alpha*rhoa)*(Ur ^ fvc::curl(U));
volVectorField liftCoeff(Cl*(beta*rhob + alpha*rhoa)*(Ur ^ fvc::curl(U)));

24
applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
View File

@ -1,21 +1,23 @@
{
surfaceScalarField alphaf = fvc::interpolate(alpha);
surfaceScalarField betaf = scalar(1) - alphaf;
surfaceScalarField alphaf(fvc::interpolate(alpha));
surfaceScalarField betaf(scalar(1) - alphaf);
volScalarField rUaA = 1.0/UaEqn.A();
volScalarField rUbA = 1.0/UbEqn.A();
volScalarField rUaA(1.0/UaEqn.A());
volScalarField rUbA(1.0/UbEqn.A());
phia == (fvc::interpolate(Ua) & mesh.Sf());
phib == (fvc::interpolate(Ub) & mesh.Sf());
rUaAf = fvc::interpolate(rUaA);
surfaceScalarField rUbAf = fvc::interpolate(rUbA);
surfaceScalarField rUbAf(fvc::interpolate(rUbA));
Ua = rUaA*UaEqn.H();
Ub = rUbA*UbEqn.H();
surfaceScalarField phiDraga =
fvc::interpolate(beta/rhoa*K*rUaA)*phib + rUaAf*(g & mesh.Sf());
surfaceScalarField phiDraga
(
fvc::interpolate(beta/rhoa*K*rUaA)*phib + rUaAf*(g & mesh.Sf())
);
if (g0.value() > 0.0)
{
@ -27,8 +29,10 @@
phiDraga -= rUaAf*fvc::snGrad(kineticTheory.pa()/rhoa)*mesh.magSf();
}
surfaceScalarField phiDragb =
fvc::interpolate(alpha/rhob*K*rUbA)*phia + rUbAf*(g & mesh.Sf());
surfaceScalarField phiDragb
(
fvc::interpolate(alpha/rhob*K*rUbA)*phia + rUbAf*(g & mesh.Sf())
);
// Fix for gravity on outlet boundary.
forAll(p.boundaryField(), patchi)
@ -66,7 +70,7 @@
if (nonOrth == nNonOrthCorr)
{
surfaceScalarField SfGradp = pEqn.flux()/Dp;
surfaceScalarField SfGradp(pEqn.flux()/Dp);
phia -= rUaAf*SfGradp/rhoa;
phib -= rUbAf*SfGradp/rhob;

4
applications/solvers/multiphase/twoPhaseEulerFoam/packingLimiter.H
View File

@ -1,11 +1,11 @@
if (packingLimiter)
{
// Calculating exceeding volume fractions
volScalarField alphaEx = max(alpha - alphaMax, scalar(0));
volScalarField alphaEx(max(alpha - alphaMax, scalar(0)));
// Finding neighbouring cells of the whole domain
labelListList neighbour = mesh.cellCells();
scalarField cellVolumes = mesh.cellVolumes();
scalarField cellVolumes(mesh.cellVolumes());
forAll (alphaEx, celli)
{