diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H index 5e916dce2..3445c81fb 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H @@ -294,12 +294,11 @@ while (pimple.correct()) forAll(phases, phasei) { phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; + volScalarField& rho = phase.thermo().rho(); if (phase.compressible()) { - const volScalarField& alpha = phase; - const volScalarField& rho = phase.rho(); - if (pimple.transonic()) { surfaceScalarField phid @@ -357,21 +356,29 @@ while (pimple.correct()) ).ptr() ); } + } + else + { + pEqnComps.set + ( + phasei, + fvm::Su(-(fvOptions(alpha, rho)&rho)/rho, p_rgh).ptr() + ); + } - if (fluid.transfersMass(phase)) + if (fluid.transfersMass(phase)) + { + if (pEqnComps.set(phasei)) { - if (pEqnComps.set(phasei)) - { - pEqnComps[phasei] -= fluid.dmdt(phase)/rho; - } - else - { - pEqnComps.set - ( - phasei, - fvm::Su(-fluid.dmdt(phase)/rho, p_rgh) - ); - } + pEqnComps[phasei] -= fluid.dmdt(phase)/rho; + } + else + { + pEqnComps.set + ( + phasei, + fvm::Su(-fluid.dmdt(phase)/rho, p_rgh) + ); } } } @@ -421,7 +428,7 @@ while (pimple.correct()) phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp; // Set the phase dilatation rates - if (phase.compressible()) + if (pEqnComps.set(phasei)) { phase.divU(-pEqnComps[phasei] & p_rgh); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H index f268beae6..3365febdf 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H @@ -230,9 +230,10 @@ while (pimple.correct()) tmp pEqnComp2; // Construct the compressibility parts of the pressure equation - if (pimple.transonic()) + + if (phase1.compressible()) { - if (phase1.compressible()) + if (pimple.transonic()) { surfaceScalarField phid1 ( @@ -257,8 +258,24 @@ while (pimple.correct()) deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr()); pEqnComp1.ref().relax(); } + else + { + pEqnComp1 = + ( + phase1.continuityError() + - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) + )/rho1 + + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)); + } + } + else + { + pEqnComp1 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh); + } - if (phase2.compressible()) + if (phase2.compressible()) + { + if (pimple.transonic()) { surfaceScalarField phid2 ( @@ -279,23 +296,11 @@ while (pimple.correct()) + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) ) ); + deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr()); pEqnComp2.ref().relax(); } - } - else - { - if (phase1.compressible()) - { - pEqnComp1 = - ( - phase1.continuityError() - - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) - )/rho1 - + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)); - } - - if (phase2.compressible()) + else { pEqnComp2 = ( @@ -305,6 +310,10 @@ while (pimple.correct()) + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)); } } + else + { + pEqnComp2 = fvm::Su(-(fvOptions(alpha2, rho2)&rho2)/rho2, p_rgh); + } if (fluid.transfersMass()) { @@ -390,11 +399,11 @@ while (pimple.correct()) } // Set the phase dilatation rates - if (phase1.compressible()) + if (pEqnComp1.valid()) { phase1.divU(-pEqnComp1 & p_rgh); } - if (phase2.compressible()) + if (pEqnComp2.valid()) { phase2.divU(-pEqnComp2 & p_rgh); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H index 05bca24c4..d2202eb89 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H @@ -215,79 +215,91 @@ while (pimple.correct()) tmp pEqnComp1; tmp pEqnComp2; - if (pimple.transonic()) - { - surfaceScalarField phid1 - ( - IOobject::groupName("phid", phase1.name()), - fvc::interpolate(psi1)*phi1 - ); - surfaceScalarField phid2 - ( - IOobject::groupName("phid", phase2.name()), - fvc::interpolate(psi2)*phi2 - ); + // Construct the compressibility parts of the pressure equation - if (phase1.compressible()) + if (phase1.compressible()) + { + if (pimple.transonic()) { - pEqnComp1 = + surfaceScalarField phid1 ( - phase1.continuityError() - - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) - )/rho1 - + correction - ( - (alpha1/rho1)* - ( - psi1*fvm::ddt(p_rgh) - + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) - ) + IOobject::groupName("phid", phase1.name()), + fvc::interpolate(psi1)*phi1 ); + + pEqnComp1 = + ( + phase1.continuityError() + - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) + )/rho1 + + correction + ( + (alpha1/rho1)* + ( + psi1*fvm::ddt(p_rgh) + + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) + ) + ); + deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr()); pEqnComp1.ref().relax(); } - - if (phase2.compressible()) + else { - pEqnComp2 = - ( - phase2.continuityError() - - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) - )/rho2 - + correction - ( - (alpha2/rho2)* + pEqnComp1 = ( - psi2*fvm::ddt(p_rgh) - + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) - ) - ); - deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr()); - pEqnComp2.ref().relax(); + phase1.continuityError() + - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) + )/rho1 + + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)); } } else { - if (phase1.compressible()) - { - pEqnComp1 = - ( - phase1.continuityError() - - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1) - )/rho1 - + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh)); - } + pEqnComp1 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh); + } - if (phase2.compressible()) + if (phase2.compressible()) + { + if (pimple.transonic()) + { + surfaceScalarField phid2 + ( + IOobject::groupName("phid", phase2.name()), + fvc::interpolate(psi2)*phi2 + ); + + pEqnComp2 = + ( + phase2.continuityError() + - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) + )/rho2 + + correction + ( + (alpha2/rho2)* + ( + psi2*fvm::ddt(p_rgh) + + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) + ) + ); + + deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr()); + pEqnComp2.ref().relax(); + } + else { pEqnComp2 = - ( - phase2.continuityError() - - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) - )/rho2 - + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)); + ( + phase2.continuityError() + - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2) + )/rho2 + + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh)); } } + else + { + pEqnComp2 = fvm::Su(-(fvOptions(alpha2, rho2)&rho2)/rho2, p_rgh); + } if (fluid.transfersMass()) {