From 8c4571a050e07e84d92f81fb3682b6b7d7cd7435 Mon Sep 17 00:00:00 2001 From: Henry Date: Fri, 24 Feb 2012 17:12:33 +0000 Subject: [PATCH] twoPhaseEulerFoam and compressibleTwoPhaseEulerFoam: Improved flux formulation for BCs --- .../compressibleTwoPhaseEulerFoam/pEqn.H | 115 ++++++++++-------- .../multiphase/twoPhaseEulerFoam/pEqn.H | 13 +- 2 files changed, 71 insertions(+), 57 deletions(-) diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H index cd2c6b7ad8..d2ef15776d 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H @@ -2,6 +2,51 @@ rho1 = rho10 + psi1*p; rho2 = rho20 + psi2*p; + surfaceScalarField alpha1f(fvc::interpolate(alpha1)); + surfaceScalarField alpha2f(scalar(1) - alpha1f); + + volScalarField rAU1(1.0/U1Eqn.A()); + volScalarField rAU2(1.0/U2Eqn.A()); + + surfaceScalarField rAlphaAU1f = fvc::interpolate(alpha1*rAU1); + surfaceScalarField rAlphaAU2f = fvc::interpolate(alpha2*rAU2); + + volVectorField HbyA1("HbyA1", U1); + HbyA1 = rAU1*U1Eqn.H(); + + volVectorField HbyA2("HbyA2", U2); + HbyA2 = rAU2*U2Eqn.H(); + + surfaceScalarField phiHbyA1 + ( + "phiHbyA1", + (fvc::interpolate(HbyA1) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1) + + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2 + + rAlphaAU1f*(g & mesh.Sf()) + ); + + surfaceScalarField phiHbyA2 + ( + "phiHbyA2", + (fvc::interpolate(HbyA2) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2) + + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1 + + rAlphaAU2f*(g & mesh.Sf()) + ); + + surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2); + + surfaceScalarField Dp + ( + "Dp", + mag + ( + alpha1f*rAlphaAU1f/fvc::interpolate(rho1) + + alpha2f*rAlphaAU2f/fvc::interpolate(rho2) + ) + ); + tmp pEqnComp1; tmp pEqnComp2; @@ -24,57 +69,11 @@ - fvc::Sp(fvc::div(phid2), p); } - surfaceScalarField alpha1f(fvc::interpolate(alpha1)); - surfaceScalarField alpha2f(scalar(1) - alpha1f); - - volVectorField U10 = U1; - volVectorField U20 = U2; - - volScalarField rAU1(1.0/U1Eqn.A()); - volScalarField rAU2(1.0/U2Eqn.A()); - - surfaceScalarField rAlphaAU1f = fvc::interpolate(alpha1*rAU1); - surfaceScalarField rAlphaAU2f = fvc::interpolate(alpha2*rAU2); - - U1 = rAU1*U1Eqn.H(); - U2 = rAU2*U2Eqn.H(); - - { - surfaceScalarField phi10 = phi1; - - phi1 = - ( - (fvc::interpolate(U1) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1) - + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2 - + rAlphaAU1f*(g & mesh.Sf()) - ); - - phi2 = - ( - (fvc::interpolate(U2) & mesh.Sf()) - + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2) - + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi10 - + rAlphaAU2f*(g & mesh.Sf()) - ); - } - - surfaceScalarField phi0("phi0", alpha1f*phi1 + alpha2f*phi2); - phi = phi0; - adjustPhi(phi, U, p); - - surfaceScalarField Dp - ( - "Dp", - mag(alpha1f*rAlphaAU1f/fvc::interpolate(rho1) - + alpha2f*rAlphaAU2f/fvc::interpolate(rho2)) - ); - while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqnIncomp ( - fvc::div(phi) + fvc::div(phiHbyA) - fvm::laplacian(Dp, p) ); @@ -85,13 +84,21 @@ + (alpha2/rho2)*pEqnComp2() ) + pEqnIncomp, - mesh.solver(p.select(pimple.finalInnerIter()))); + mesh.solver(p.select(pimple.finalInnerIter())) + ); if (pimple.finalNonOrthogonalIter()) { surfaceScalarField mSfGradp = pEqnIncomp.flux()/Dp; - phi1 += rAlphaAU1f*mSfGradp/fvc::interpolate(rho1); - phi2 += rAlphaAU2f*mSfGradp/fvc::interpolate(rho2); + + phi1.boundaryField() == + (fvc::interpolate(U1) & mesh.Sf())().boundaryField(); + phi1 = phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1); + + phi2.boundaryField() == + (fvc::interpolate(U2) & mesh.Sf())().boundaryField(); + phi2 = phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2); + phi = alpha1f*phi1 + alpha2f*phi2; dgdt = @@ -103,7 +110,10 @@ p.relax(); mSfGradp = pEqnIncomp.flux()/Dp; - U1 += (1.0/rho1)*rAU1*dragCoeff*U20 + volVectorField U10 = U1; + + U1 = HbyA1 + + (1.0/rho1)*rAU1*dragCoeff*U2 + fvc::reconstruct ( rAlphaAU1f*(g & mesh.Sf()) @@ -111,7 +121,8 @@ ); U1.correctBoundaryConditions(); - U2 += (1.0/rho2)*rAU2*dragCoeff*U10 + U2 = HbyA2 + + (1.0/rho2)*rAU2*dragCoeff*U10 + fvc::reconstruct ( rAlphaAU2f*(g & mesh.Sf()) diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H index b5d57d9c61..78e5587a3a 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H @@ -5,11 +5,6 @@ volScalarField rUaA(1.0/UaEqn.A()); volScalarField rUbA(1.0/UbEqn.A()); - phia.boundaryField() == - (fvc::interpolate(Ua) & mesh.Sf())().boundaryField(); - phib.boundaryField() == - (fvc::interpolate(Ub) & mesh.Sf())().boundaryField(); - rUaAf = fvc::interpolate(rUaA); surfaceScalarField rUbAf(fvc::interpolate(rUbA)); @@ -51,6 +46,7 @@ surfaceScalarField phiHbyAa ( + "phiHbyAa", (fvc::interpolate(HbyAa) & mesh.Sf()) + fvc::ddtPhiCorr(rUaA, Ua, phia) + phiDraga @@ -58,6 +54,7 @@ surfaceScalarField phiHbyAb ( + "phiHbyAb", (fvc::interpolate(HbyAb) & mesh.Sf()) + fvc::ddtPhiCorr(rUbA, Ub, phib) + phiDragb @@ -86,8 +83,14 @@ { surfaceScalarField SfGradp(pEqn.flux()/Dp); + phia.boundaryField() == + (fvc::interpolate(Ua) & mesh.Sf())().boundaryField(); phia = phiHbyAa - rUaAf*SfGradp/rhoa; + + phib.boundaryField() == + (fvc::interpolate(Ub) & mesh.Sf())().boundaryField(); phib = phiHbyAb - rUbAf*SfGradp/rhob; + phi = alphaf*phia + betaf*phib; p.relax();