From b94d46d0e7bd8c2a85fca313758419214c41e6c4 Mon Sep 17 00:00:00 2001 From: Henry Date: Thu, 2 Aug 2012 15:29:08 +0100 Subject: [PATCH] multiphase and fireFoam: set phi based on the interpolated U before the pEqn construction for the p BCs --- .../solvers/combustion/fireFoam/pEqn.H | 4 +-- .../solvers/multiphase/bubbleFoam/pEqn.H | 12 +++---- .../multiphase/compressibleInterFoam/pEqn.H | 3 +- .../compressibleTwoPhaseEulerFoam/pEqn.H | 35 +++++++++++-------- .../multiphase/interFoam/MRFInterFoam/pEqn.H | 9 ++--- .../multiphase/interFoam/interDyMFoam/pEqn.H | 4 +-- .../solvers/multiphase/interFoam/pEqn.H | 3 +- .../multiphase/interPhaseChangeFoam/pEqn.H | 4 +-- .../multiphase/multiphaseEulerFoam/pEqn.H | 15 +++++--- .../MRFMultiphaseInterFoam/pEqn.H | 9 ++--- .../multiphase/multiphaseInterFoam/pEqn.H | 4 +-- .../solvers/multiphase/settlingFoam/pEqn.H | 4 +-- .../multiphase/twoLiquidMixingFoam/pEqn.H | 4 +-- .../multiphase/twoPhaseEulerFoam/pEqn.H | 29 ++++++++------- 14 files changed, 64 insertions(+), 75 deletions(-) diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H index eb263b8a50..22c3f0139a 100644 --- a/applications/solvers/combustion/fireFoam/pEqn.H +++ b/applications/solvers/combustion/fireFoam/pEqn.H @@ -5,8 +5,8 @@ surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU)); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); phi.boundaryField() = - fvc::interpolate(rho.boundaryField()) - *(fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField()); + fvc::interpolate(rho.boundaryField()*U.boundaryField()) + & mesh.Sf().boundaryField(); surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); diff --git a/applications/solvers/multiphase/bubbleFoam/pEqn.H b/applications/solvers/multiphase/bubbleFoam/pEqn.H index af3f6e2822..d438b1124f 100644 --- a/applications/solvers/multiphase/bubbleFoam/pEqn.H +++ b/applications/solvers/multiphase/bubbleFoam/pEqn.H @@ -14,12 +14,6 @@ volVectorField HbyA2("HbyA2", U2); HbyA2 = rAU2*U2Eqn.H(); - U.boundaryField() = - alpha1.boundaryField()*U1.boundaryField() - + alpha2.boundaryField()*U2.boundaryField(); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); - surfaceScalarField phiDrag1 ( fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2 @@ -44,16 +38,18 @@ ( (fvc::interpolate(HbyA1) & mesh.Sf()) + fvc::ddtPhiCorr(rAU1, U1, phi1) - + phiDrag1 ); surfaceScalarField phiHbyA2 ( (fvc::interpolate(HbyA2) & mesh.Sf()) + fvc::ddtPhiCorr(rAU2, U2, phi2) - + phiDrag2 ); + phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2; + + phiHbyA1 += phiDrag1; + phiHbyA2 += phiDrag2; surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2); surfaceScalarField Dp diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H index e0c13b16e0..80011d1ab5 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H @@ -29,8 +29,6 @@ volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); surfaceScalarField phiHbyA ( @@ -38,6 +36,7 @@ (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); + phi = phiHbyA; surfaceScalarField phig ( diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H index 915e714082..926abe2a09 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H @@ -2,11 +2,6 @@ rho1 = rho10 + psi1*p; rho2 = rho20 + psi2*p; - mrfZones.absoluteFlux(phi1.oldTime()); - mrfZones.absoluteFlux(phi1); - mrfZones.absoluteFlux(phi2.oldTime()); - mrfZones.absoluteFlux(phi2); - surfaceScalarField alpha1f(fvc::interpolate(alpha1)); surfaceScalarField alpha2f(scalar(1) - alpha1f); @@ -22,11 +17,10 @@ volVectorField HbyA2("HbyA2", U2); HbyA2 = rAU2*U2Eqn.H(); - U.boundaryField() = - alpha1.boundaryField()*U1.boundaryField() - + alpha2.boundaryField()*U2.boundaryField(); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); + mrfZones.absoluteFlux(phi1.oldTime()); + mrfZones.absoluteFlux(phi1); + mrfZones.absoluteFlux(phi2.oldTime()); + mrfZones.absoluteFlux(phi2); surfaceScalarField ppDrag("ppDrag", 0.0*phi1); @@ -47,18 +41,29 @@ "phiHbyA1", (fvc::interpolate(HbyA1) & mesh.Sf()) + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1) - + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2 - + ppDrag - + rAlphaAU1f*(g & mesh.Sf()) ); - mrfZones.relativeFlux(phiHbyA1); surfaceScalarField phiHbyA2 ( "phiHbyA2", (fvc::interpolate(HbyA2) & mesh.Sf()) + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2) - + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1 + ); + + phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2; + mrfZones.relativeFlux(phi); + + phiHbyA1 += + ( + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2 + + ppDrag + + rAlphaAU1f*(g & mesh.Sf()) + ); + mrfZones.relativeFlux(phiHbyA1); + + phiHbyA2 += + ( + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1 + rAlphaAU2f*(g & mesh.Sf()) ); mrfZones.relativeFlux(phiHbyA2); diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H index 28de546ce0..1f4f0bbbeb 100644 --- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H +++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H @@ -5,20 +5,15 @@ volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - mrfZones.absoluteFlux(phi); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); - mrfZones.relativeFlux(phi); - surfaceScalarField phiHbyA ( "phiHbyA", (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - mrfZones.relativeFlux(phiHbyA); - adjustPhi(phiHbyA, U, p_rgh); + mrfZones.relativeFlux(phiHbyA); + phi = phiHbyA; surfaceScalarField phig ( diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H index 42721c0018..5b089550c5 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H @@ -4,8 +4,6 @@ volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); surfaceScalarField phiHbyA ( @@ -21,6 +19,8 @@ fvc::makeAbsolute(phiHbyA, U); } + phiAbs = phiHbyA; + surfaceScalarField phig ( ( diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H index a21d5b7bc8..53c2161f97 100644 --- a/applications/solvers/multiphase/interFoam/pEqn.H +++ b/applications/solvers/multiphase/interFoam/pEqn.H @@ -4,8 +4,6 @@ volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); surfaceScalarField phiHbyA ( @@ -15,6 +13,7 @@ ); adjustPhi(phiHbyA, U, p_rgh); + phi = phiHbyA; surfaceScalarField phig ( diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H index aa3ed1fe85..7ff971ed7f 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H @@ -4,8 +4,6 @@ volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); surfaceScalarField phiHbyA ( @@ -13,8 +11,8 @@ (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - adjustPhi(phiHbyA, U, p_rgh); + phi = phiHbyA; surfaceScalarField phig ( diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H index 9ed3fecce9..07d867b1e7 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H @@ -58,6 +58,8 @@ dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0) ); + phi = dimensionedScalar("phi", phi.dimensions(), 0); + phasei = 0; forAllIter(PtrDictionary, fluid.phases(), iter) { @@ -74,14 +76,17 @@ ( (fvc::interpolate(HbyAs[phasei]) & mesh.Sf()) + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi()) - + rAlphaAUfs[phasei] + ); + mrfZones.relativeFlux(phiHbyAs[phasei]); + + phi += alphafs[phasei]*phiHbyAs[phasei]; + + phiHbyAs[phasei] += + rAlphaAUfs[phasei] *( fluid.surfaceTension(phase)*mesh.magSf()/phase.rho() + (g & mesh.Sf()) - ) - ); - - mrfZones.relativeFlux(phiHbyAs[phasei]); + ); multiphaseSystem::dragModelTable::const_iterator dmIter = fluid.dragModels().begin(); diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H index 5e87fa0319..823216897c 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H @@ -5,20 +5,15 @@ volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - mrfZones.absoluteFlux(phi); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); - mrfZones.relativeFlux(phi); - surfaceScalarField phiHbyA ( "phiHbyA", (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - mrfZones.relativeFlux(phiHbyA); - adjustPhi(phiHbyA, U, p_rgh); + mrfZones.relativeFlux(phiHbyA); + phi = phiHbyA; surfaceScalarField phig ( diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H index 2e26fb538a..82bf3dbd59 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H +++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H @@ -4,8 +4,6 @@ volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); surfaceScalarField phiHbyA ( @@ -13,8 +11,8 @@ (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - adjustPhi(phiHbyA, U, p_rgh); + phi = phiHbyA; surfaceScalarField phig ( diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H index 99c15656dc..6d61ec8730 100644 --- a/applications/solvers/multiphase/settlingFoam/pEqn.H +++ b/applications/solvers/multiphase/settlingFoam/pEqn.H @@ -5,9 +5,6 @@ volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - phi.boundaryField() = - fvc::interpolate(rho.boundaryField()*U.boundaryField()) - & mesh.Sf().boundaryField(); surfaceScalarField phiHbyA ( @@ -18,6 +15,7 @@ + fvc::ddtPhiCorr(rAU, rho, U, phi) ) ); + phi = phiHbyA; surfaceScalarField phig ( diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H index 810a753c9d..e090522de6 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H @@ -4,8 +4,6 @@ volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H(); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); surfaceScalarField phiHbyA ( @@ -13,8 +11,8 @@ (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::ddtPhiCorr(rAU, rho, U, phi) ); - adjustPhi(phiHbyA, U, p_rgh); + phi = phiHbyA; surfaceScalarField phig ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H index d6bde3812b..546e1a3e41 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H @@ -18,13 +18,6 @@ mrfZones.absoluteFlux(phi1); mrfZones.absoluteFlux(phi2.oldTime()); mrfZones.absoluteFlux(phi2); - mrfZones.absoluteFlux(phi); - - U.boundaryField() = - alpha1.boundaryField()*U1.boundaryField() - + alpha2.boundaryField()*U2.boundaryField(); - phi.boundaryField() = - fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField(); surfaceScalarField ppDrag("ppDrag", 0.0*phi1); @@ -43,18 +36,29 @@ "phiHbyA1", (fvc::interpolate(HbyA1) & mesh.Sf()) + fvc::ddtPhiCorr(rAU1, U1, phi1) - + fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 - + ppDrag - + rAU1f*(g & mesh.Sf()) ); - mrfZones.relativeFlux(phiHbyA1); surfaceScalarField phiHbyA2 ( "phiHbyA2", (fvc::interpolate(HbyA2) & mesh.Sf()) + fvc::ddtPhiCorr(rAU2, U2, phi2) - + fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + ); + + phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2; + mrfZones.relativeFlux(phi); + + phiHbyA1 += + ( + fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + + ppDrag + + rAU1f*(g & mesh.Sf()) + ); + mrfZones.relativeFlux(phiHbyA1); + + phiHbyA2 += + ( + fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + rAU2f*(g & mesh.Sf()) ); mrfZones.relativeFlux(phiHbyA2); @@ -63,7 +67,6 @@ mrfZones.relativeFlux(phi1); mrfZones.relativeFlux(phi2.oldTime()); mrfZones.relativeFlux(phi2); - mrfZones.relativeFlux(phi); surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);