diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H index 08282e78b3..a882a58249 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H @@ -98,60 +98,3 @@ Info<< "Creating field kinetic energy K\n" << endl; volScalarField K1(IOobject::groupName("K", phase1.name()), 0.5*magSqr(U1)); volScalarField K2(IOobject::groupName("K", phase2.name()), 0.5*magSqr(U2)); - - - volScalarField rAU1 - ( - IOobject - ( - IOobject::groupName("rAU", phase1.name()), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0) - ); - - volScalarField rAU2 - ( - IOobject - ( - IOobject::groupName("rAU", phase2.name()), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0) - ); - - surfaceScalarField rAU1f - ( - IOobject - ( - IOobject::groupName("rAUf", phase1.name()), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0) - ); - - surfaceScalarField rAU2f - ( - IOobject - ( - IOobject::groupName("rAUf", phase2.name()), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0) - ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H index fe3ba4350a..1f84d197f5 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H @@ -1,11 +1,65 @@ -surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1)); -surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f); +surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1)); +surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1); -rAU1 = 1.0/U1Eqn.A(); -rAU2 = 1.0/U2Eqn.A(); +volScalarField rAU1(IOobject::groupName("rAU", phase1.name()), 1.0/U1Eqn.A()); +volScalarField rAU2(IOobject::groupName("rAU", phase2.name()), 1.0/U2Eqn.A()); + +surfaceScalarField alpharAUf1(fvc::interpolate(alpha1*rAU1)); +surfaceScalarField alpharAUf2(fvc::interpolate(alpha2*rAU2)); + +// Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes +tmp phiF1; +tmp phiF2; +{ + // Turbulent-dispersion diffusivity + volScalarField D(fluid.D()); + + // Phase-1 turbulent dispersion and particle-pressure flux + surfaceScalarField Df1 + ( + fvc::interpolate + ( + rAU1*(D + phase1.turbulence().pPrime()) + ) + ); + + // Phase-2 turbulent dispersion and particle-pressure flux + surfaceScalarField Df2 + ( + fvc::interpolate + ( + rAU2*(D + phase2.turbulence().pPrime()) + ) + ); + + // Cache the net diffusivity for implicit diffusion treatment in the + // phase-fraction equation + if (implicitPhasePressure) + { + fluid.pPrimeByA() = Df1 + Df2; + } + + // Lift and wall-lubrication forces + volVectorField F(fluid.F()); + + // Phase-fraction face-gradient + surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf()); + + // Phase-1 dispersion, lift and wall-lubrication flux + phiF1 = + ( + Df1*snGradAlpha1 + + (fvc::interpolate(rAU1*F) & mesh.Sf()) + ); + + // Phase-1 dispersion, lift and wall-lubrication flux + phiF2 = + ( + - Df2*snGradAlpha1 + - (fvc::interpolate(rAU2*F) & mesh.Sf()) + ); +} -surfaceScalarField alpharAU1f(fvc::interpolate(alpha1*rAU1)); -surfaceScalarField alpharAU2f(fvc::interpolate(alpha2*rAU2)); // --- Pressure corrector loop while (pimple.correct()) @@ -33,61 +87,38 @@ while (pimple.correct()) ); HbyA2 = rAU2*U2Eqn.H(); - // Face force fluxes - tmp phiF1; - tmp phiF2; - - // Turbulent diffusion, particle-pressure lift and wall-lubrication fluxes - { - volScalarField D(fluid.D()); - volVectorField F(fluid.F()); - surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf()); - - phiF1 = - ( - fvc::interpolate - ( - rAU1*(D + phase1.turbulence().pPrime()) - )*snGradAlpha1 - + (fvc::interpolate(rAU1*F) & mesh.Sf()) - ); - - phiF2 = - ( - - fvc::interpolate - ( - rAU2*(D + phase2.turbulence().pPrime()) - )*snGradAlpha1 - - (fvc::interpolate(rAU2*F) & mesh.Sf()) - ); - } - // Mean density for buoyancy force and p_rgh -> p volScalarField rho("rho", fluid.rho()); - // Add Buoyancy force - { - surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf()); + surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + ghf*fvc::snGrad(rho)*mesh.magSf() + ); - phiF1() += - alpharAU1f - *( - ghSnGradRho - - alpha2f*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf()) - ); + surfaceScalarField phig1 + ( + alpharAUf1 + *( + ghSnGradRho + - alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf()) + ) + ); + + surfaceScalarField phig2 + ( + alpharAUf2 + *( + ghSnGradRho + - alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf()) + ) + ); - phiF2() += - alpharAU2f - *( - ghSnGradRho - - alpha1f*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf()) - ); - } // ddtPhiCorr filter -- only apply in pure(ish) phases - surfaceScalarField alpha1fBar(fvc::interpolate(fvc::average(alpha1f))); - surfaceScalarField phiCorrCoeff1(pos(alpha1fBar - 0.99)); - surfaceScalarField phiCorrCoeff2(pos(0.01 - alpha1fBar)); + surfaceScalarField alphaf1Bar(fvc::interpolate(fvc::average(alphaf1))); + surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99)); + surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar)); forAll(mesh.boundary(), patchi) { @@ -114,6 +145,7 @@ while (pimple.correct()) - (fvc::interpolate(U1.oldTime()) & mesh.Sf()) )/runTime.deltaT() - phiF1() + - phig1 ); // Phase-2 predicted flux @@ -127,19 +159,20 @@ while (pimple.correct()) - (fvc::interpolate(U2.oldTime()) & mesh.Sf()) )/runTime.deltaT() - phiF2() + - phig2 ); // Face-drag coefficients - surfaceScalarField D1f(fvc::interpolate(rAU1*Kd)); - surfaceScalarField D2f(fvc::interpolate(rAU2*Kd)); + surfaceScalarField rAUKd1(fvc::interpolate(rAU1*Kd)); + surfaceScalarField rAUKd2(fvc::interpolate(rAU2*Kd)); // Construct the mean predicted flux // including explicit drag contributions based on absolute fluxes surfaceScalarField phiHbyA ( "phiHbyA", - alpha1f*(phiHbyA1 + D1f*mrfZones.absolute(phi2)) - + alpha2f*(phiHbyA2 + D2f*mrfZones.absolute(phi1)) + alphaf1*(phiHbyA1 + rAUKd1*mrfZones.absolute(phi2)) + + alphaf2*(phiHbyA2 + rAUKd2*mrfZones.absolute(phi1)) ); mrfZones.makeRelative(phiHbyA); @@ -147,7 +180,7 @@ while (pimple.correct()) surfaceScalarField rAUf ( "rAUf", - mag(alpha1f*alpharAU1f + alpha2f*alpharAU2f) + mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2) ); // Update the fixedFluxPressure BCs to ensure flux consistency @@ -158,9 +191,9 @@ while (pimple.correct()) phiHbyA.boundaryField() - mrfZones.relative ( - alpha1f.boundaryField() + alphaf1.boundaryField() *(mesh.Sf().boundaryField() & U1.boundaryField()) - + alpha2f.boundaryField() + + alphaf2.boundaryField() *(mesh.Sf().boundaryField() & U2.boundaryField()) ) )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) @@ -256,21 +289,22 @@ while (pimple.correct()) { surfaceScalarField phi1s ( - phiHbyA1 + alpharAU1f*mSfGradp + phiHbyA1 + alpharAUf1*mSfGradp ); surfaceScalarField phi2s ( - phiHbyA2 + alpharAU2f*mSfGradp + phiHbyA2 + alpharAUf2*mSfGradp ); surfaceScalarField phir ( - ((phi1s + D1f*phi2s) - (phi2s + D2f*phi1s))/(1 - D1f*D2f) + ((phi1s + rAUKd1*phi2s) - (phi2s + rAUKd2*phi1s)) + /(1 - rAUKd1*rAUKd2) ); - phi1 = phi + alpha2f*phir; - phi2 = phi - alpha1f*phir; + phi1 = phi + alphaf2*phir; + phi2 = phi - alphaf1*phir; } // Compressibility correction for phase-fraction equations @@ -287,21 +321,23 @@ while (pimple.correct()) // Partial-elimination phase-velocity corrector { - volVectorField U1s + volVectorField Us1 ( - HbyA1 + fvc::reconstruct(alpharAU1f*mSfGradp - phiF1) + HbyA1 + + fvc::reconstruct(alpharAUf1*mSfGradp - phiF1() - phig1) ); - volVectorField U2s + volVectorField Us2 ( - HbyA2 + fvc::reconstruct(alpharAU2f*mSfGradp - phiF2) + HbyA2 + + fvc::reconstruct(alpharAUf2*mSfGradp - phiF2() - phig2) ); volScalarField D1(rAU1*Kd); volScalarField D2(rAU2*Kd); - U = alpha1*(U1s + D1*U2) + alpha2*(U2s + D2*U1); - volVectorField Ur(((1 - D2)*U1s - (1 - D1)*U2s)/(1 - D1*D2)); + U = alpha1*(Us1 + D1*U2) + alpha2*(Us2 + D2*U1); + volVectorField Ur(((1 - D2)*Us1 - (1 - D1)*Us2)/(1 - D1*D2)); U1 = U + alpha2*Ur; U1.correctBoundaryConditions(); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H index 0e5a47a8ff..5d39ce0d7c 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H @@ -1,9 +1,9 @@ -surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1)); -surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f); +surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1)); +surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1); -surfaceScalarField alphaRho1f0 +surfaceScalarField alphaRhof10 ( - "alphaRho1f0", + "alphaRhof10", fvc::interpolate ( max(alpha1.oldTime(), fluid.residualAlpha(phase1)) @@ -11,9 +11,9 @@ surfaceScalarField alphaRho1f0 ) ); -surfaceScalarField alphaRho2f0 +surfaceScalarField alphaRhof20 ( - "alphaRho2f0", + "alphaRhof20", fvc::interpolate ( max(alpha2.oldTime(), fluid.residualAlpha(phase2)) @@ -27,34 +27,76 @@ surfaceScalarField Kdf("Kdf", fluid.Kdf()); // Virtual-mass coefficient surfaceScalarField Vmf("Vmf", fluid.Vmf()); -// Lift and wall-lubrication forces -surfaceScalarField Ff("Ff", fluid.Ff()); - -tmp snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf()); -tmp D(fluid.D()); - -// Turbulent-dispersion force for phase 1 -surfaceScalarField Ftdf1 +surfaceScalarField rAUf1 ( - IOobject::groupName("Ftdf", phase1.name()), - fvc::interpolate(D() + phase1.turbulence().pPrime())*snGradAlpha1() + IOobject::groupName("rAUf", phase1.name()), + 1.0 + /( + (alphaRhof10 + Vmf)/runTime.deltaT() + + fvc::interpolate(U1Eqn.A()) + + Kdf + ) ); -// Turbulent-dispersion force for phase 2 -surfaceScalarField Ftdf2 +surfaceScalarField rAUf2 ( - IOobject::groupName("Ftdf", phase2.name()), - fvc::interpolate(D + phase2.turbulence().pPrime())*snGradAlpha1 + IOobject::groupName("rAUf", phase2.name()), + 1.0 + /( + (alphaRhof20 + Vmf)/runTime.deltaT() + + fvc::interpolate(U2Eqn.A()) + + Kdf + ) ); +// Turbulent dispersion, particle-pressure, lift and wall-lubrication forces +tmp Ff1; +tmp Ff2; +{ + // Turbulent-dispersion diffusivity + volScalarField D(fluid.D()); + + // Phase-1 turbulent dispersion and particle-pressure diffusivity + surfaceScalarField Df1 + ( + fvc::interpolate(D + phase1.turbulence().pPrime()) + ); + + // Phase-2 turbulent dispersion and particle-pressure diffusivity + surfaceScalarField Df2 + ( + fvc::interpolate(D + phase2.turbulence().pPrime()) + ); + + // Cache the net diffusivity for implicit diffusion treatment in the + // phase-fraction equation + if (implicitPhasePressure) + { + fluid.pPrimeByA() = rAUf1*Df1 + rAUf2*Df2; + } + + // Lift and wall-lubrication forces + surfaceScalarField Ff(fluid.Ff()); + + // Phase-fraction face-gradient + surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf()); + + // Phase-1 dispersion, lift and wall-lubrication force + Ff1 = Df1*snGradAlpha1 + Ff; + + // Phase-2 dispersion, lift and wall-lubrication force + Ff2 = -Df2*snGradAlpha1 - Ff; +} + + while (pimple.correct()) { // Update continuity errors due to temperature changes #include "correctContErrs.H" - surfaceScalarField rho1f(fvc::interpolate(rho1)); - surfaceScalarField rho2f(fvc::interpolate(rho2)); + surfaceScalarField rhof1(fvc::interpolate(rho1)); + surfaceScalarField rhof2(fvc::interpolate(rho2)); // Correct flux BCs to be consistent with the velocity BCs phi1.boundaryField() == @@ -62,41 +104,18 @@ while (pimple.correct()) phi2.boundaryField() == mrfZones.relative(mesh.Sf().boundaryField() & U2.boundaryField()); - rAU1f = + surfaceScalarField alpharAUf1 ( - IOobject::groupName("rAUf", phase1.name()), - 1.0 - /( - (alphaRho1f0 + Vmf)/runTime.deltaT() - + fvc::interpolate(U1Eqn.A()) - + Kdf - ) + IOobject::groupName("alpharAUf", phase1.name()), + max(alphaf1, fluid.residualAlpha(phase1))*rAUf1 ); - rAU2f = + surfaceScalarField alpharAUf2 ( - IOobject::groupName("rAUf", phase2.name()), - 1.0 - /( - (alphaRho2f0 + Vmf)/runTime.deltaT() - + fvc::interpolate(U2Eqn.A()) - + Kdf - ) + IOobject::groupName("alpharAUf", phase2.name()), + max(alphaf2, fluid.residualAlpha(phase2))*rAUf2 ); - surfaceScalarField rAlphaAU1f - ( - IOobject::groupName("rAlphaAUf", phase1.name()), - max(alpha1f, fluid.residualAlpha(phase1))*rAU1f - ); - - surfaceScalarField rAlphaAU2f - ( - IOobject::groupName("rAlphaAUf", phase2.name()), - max(alpha2f, fluid.residualAlpha(phase2))*rAU2f - ); - - volScalarField rho("rho", fluid.rho()); surfaceScalarField ghSnGradRho ( @@ -104,25 +123,25 @@ while (pimple.correct()) ghf*fvc::snGrad(rho)*mesh.magSf() ); - // Add the phase-1 buoyancy force - surfaceScalarField phiF1 + // Phase-1 buoyancy flux + surfaceScalarField phig1 ( - IOobject::groupName("phiF", phase1.name()), - rAlphaAU1f + IOobject::groupName("phig", phase1.name()), + alpharAUf1 *( ghSnGradRho - - alpha2f*(rho1f - rho2f)*(g & mesh.Sf()) + - alphaf2*(rhof1 - rhof2)*(g & mesh.Sf()) ) ); - // Add the phase-2 buoyancy force - surfaceScalarField phiF2 + // Phase-2 buoyancy flux + surfaceScalarField phig2 ( - IOobject::groupName("phiF", phase2.name()), - rAlphaAU2f + IOobject::groupName("phig", phase2.name()), + alpharAUf2 *( ghSnGradRho - - alpha1f*(rho2f - rho1f)*(g & mesh.Sf()) + - alphaf1*(rhof2 - rhof1)*(g & mesh.Sf()) ) ); @@ -135,15 +154,14 @@ while (pimple.correct()) ); phiHbyA1 = - rAU1f + rAUf1 *( - (alphaRho1f0 + Vmf) + (alphaRhof10 + Vmf) *mrfZones.absolute(phi1.oldTime())/runTime.deltaT() + (fvc::interpolate(U1Eqn.H()) & mesh.Sf()) + Vmf*ddtPhi2 + Kdf*mrfZones.absolute(phi2) - - Ff - - Ftdf1 + - Ff1() ); // Phase-2 predicted flux @@ -154,32 +172,31 @@ while (pimple.correct()) ); phiHbyA2 = - rAU2f + rAUf2 *( - (alphaRho2f0 + Vmf) + (alphaRhof20 + Vmf) *mrfZones.absolute(phi2.oldTime())/runTime.deltaT() + (fvc::interpolate(U2Eqn.H()) & mesh.Sf()) + Vmf*ddtPhi1 + Kdf*mrfZones.absolute(phi1) - + Ff - + Ftdf2 + - Ff2() ); surfaceScalarField phiHbyA ( "phiHbyA", - alpha1f*(phiHbyA1 - phiF1) + alpha2f*(phiHbyA2 - phiF2) + alphaf1*(phiHbyA1 - phig1) + alphaf2*(phiHbyA2 - phig2) ); mrfZones.makeRelative(phiHbyA); - phiHbyA1 -= phiF1; - phiHbyA2 -= phiF2; + phiHbyA1 -= phig1; + phiHbyA2 -= phig2; surfaceScalarField rAUf ( "rAUf", - mag(alpha1f*rAlphaAU1f + alpha2f*rAlphaAU2f) + mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2) ); // Update the fixedFluxPressure BCs to ensure flux consistency @@ -189,8 +206,8 @@ while (pimple.correct()) ( phiHbyA.boundaryField() - ( - alpha1f.boundaryField()*phi1.boundaryField() - + alpha2f.boundaryField()*phi2.boundaryField() + alphaf1.boundaryField()*phi1.boundaryField() + + alphaf2.boundaryField()*phi2.boundaryField() ) )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) ); @@ -280,25 +297,25 @@ while (pimple.correct()) surfaceScalarField phi1s ( phiHbyA1 - + rAlphaAU1f*mSfGradp - - rAU1f*Kdf*mrfZones.absolute(phi2) + + alpharAUf1*mSfGradp + - rAUf1*Kdf*mrfZones.absolute(phi2) ); surfaceScalarField phi2s ( phiHbyA2 - + rAlphaAU2f*mSfGradp - - rAU2f*Kdf*mrfZones.absolute(phi1) + + alpharAUf2*mSfGradp + - rAUf2*Kdf*mrfZones.absolute(phi1) ); surfaceScalarField phir ( - ((phi2s + rAU2f*Kdf*phi1s) - (phi1s + rAU1f*Kdf*phi2s)) - /(1.0 - rAU1f*rAU2f*sqr(Kdf)) + ((phi2s + rAUf2*Kdf*phi1s) - (phi1s + rAUf1*Kdf*phi2s)) + /(1.0 - rAUf1*rAUf2*sqr(Kdf)) ); - phi1 = phi - alpha2f*phir; - phi2 = phi + alpha1f*phir; + phi1 = phi - alphaf2*phir; + phi2 = phi + alphaf1*phir; U1 = fvc::reconstruct(mrfZones.absolute(phi1)); U1.correctBoundaryConditions(); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C index c9679c3439..adb37bc0b3 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C @@ -50,11 +50,6 @@ int main(int argc, char *argv[]) pimpleControl pimple(mesh); - Switch faceMomentum - ( - pimple.dict().lookupOrDefault("faceMomentum", false) - ); - #include "createFields.H" #include "createMRFZones.H" #include "createFvOptions.H" @@ -63,6 +58,19 @@ int main(int argc, char *argv[]) #include "CourantNos.H" #include "setInitialDeltaT.H" + Switch faceMomentum + ( + pimple.dict().lookupOrDefault("faceMomentum", false) + ); + + Switch implicitPhasePressure + ( + mesh.solverDict(alpha1.name()).lookupOrDefault + ( + "implicitPhasePressure", false + ) + ); + #include "pUf/createDDtU.H" #include "pU/createDDtU.H" diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C index 391b38a76a..59b727dad5 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C @@ -365,13 +365,6 @@ void Foam::twoPhaseSystem::solve() const surfaceScalarField& phi1 = phase1_.phi(); const surfaceScalarField& phi2 = phase2_.phi(); - Switch faceMomentum - ( - //pimple.dict().lookupOrDefault("faceMomentum", false) - mesh_.solutionDict().subDict("PIMPLE") - .lookupOrDefault("faceMomentum", false) - ); - const dictionary& alphaControls = mesh_.solverDict ( alpha1.name() @@ -379,10 +372,6 @@ void Foam::twoPhaseSystem::solve() label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles"))); label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr"))); - Switch implicitPhasePressure - ( - alphaControls.lookupOrDefault("implicitPhasePressure", false) - ); word alphaScheme("div(phi," + alpha1.name() + ')'); word alpharScheme("div(phir," + alpha1.name() + ')'); @@ -395,48 +384,11 @@ void Foam::twoPhaseSystem::solve() surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0)))); - tmp pPrimeByA; - - if (implicitPhasePressure) + if (pPrimeByA_.valid()) { - if (faceMomentum) - { - const surfaceScalarField& rAU1f = - mesh_.lookupObject - ( - IOobject::groupName("rAUf", phase1_.name()) - ); - const surfaceScalarField& rAU2f = - mesh_.lookupObject - ( - IOobject::groupName("rAUf", phase2_.name()) - ); - - volScalarField D(this->D()); - - pPrimeByA = - rAU1f*fvc::interpolate(D + phase1_.turbulence().pPrime()) - + rAU2f*fvc::interpolate(D + phase2_.turbulence().pPrime()); - } - else - { - const volScalarField& rAU1 = mesh_.lookupObject - ( - IOobject::groupName("rAU", phase1_.name()) - ); - const volScalarField& rAU2 = mesh_.lookupObject - ( - IOobject::groupName("rAU", phase2_.name()) - ); - - pPrimeByA = - fvc::interpolate(rAU1*phase1_.turbulence().pPrime()) - + fvc::interpolate(rAU2*phase2_.turbulence().pPrime()); - } - surfaceScalarField phiP ( - pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf() + pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf() ); phic += alpha1f*phiP; @@ -571,12 +523,12 @@ void Foam::twoPhaseSystem::solve() phase1_.alphaPhi() = alphaPhic1; } - if (implicitPhasePressure) + if (pPrimeByA_.valid()) { fvScalarMatrix alpha1Eqn ( fvm::ddt(alpha1) - fvc::ddt(alpha1) - - fvm::laplacian(alpha1f*pPrimeByA(), alpha1, "bounded") + - fvm::laplacian(alpha1f*pPrimeByA_(), alpha1, "bounded") ); alpha1Eqn.relax(); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H index 974af111c1..c6c6a95125 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H @@ -78,9 +78,12 @@ class twoPhaseSystem //- Total volumetric flux surfaceScalarField phi_; - //- Dilatation term + //- Dilatation term volScalarField dgdt_; + //- Optional dispersion diffusivity + tmp pPrimeByA_; + //- Unordered phase pair autoPtr pair_; @@ -113,7 +116,7 @@ class twoPhaseSystem autoPtr > turbulentDispersion_; - //- + //- Default residual alpha (0) static dimensionedScalar zeroResidualAlpha_; @@ -227,6 +230,9 @@ public: //- Return non-const access to the dilatation parameter inline volScalarField& dgdt(); + + //- Return non-const access to the dispersion diffusivity + inline tmp& pPrimeByA(); }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystemI.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystemI.H index e832b0fa1e..94bc4d7f92 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystemI.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystemI.H @@ -95,4 +95,10 @@ inline Foam::volScalarField& Foam::twoPhaseSystem::dgdt() } +inline Foam::tmp& Foam::twoPhaseSystem::pPrimeByA() +{ + return pPrimeByA_; +} + + // ************************************************************************* //