From 862fa9e10674d63383a03b2499f53d015443c97c Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Sat, 9 Dec 2017 21:03:59 +0000 Subject: [PATCH] compressibleInterFoam family: Corrected transonic option Resolves bug-report https://bugs.openfoam.org/view.php?id=2785 --- .../compressibleInterDyMFoam/pEqn.H | 42 +++++++++++++------ .../compressibleInterFilmFoam/pEqn.H | 42 +++++++++++++------ .../multiphase/compressibleInterFoam/pEqn.H | 42 +++++++++++++------ .../reactingMultiphaseEulerFoam/pU/pEqn.H | 4 -- .../reactingTwoPhaseEulerFoam/pU/pEqn.H | 2 - .../reactingTwoPhaseEulerFoam/pUf/pEqn.H | 2 - .../multiphase/twoPhaseEulerFoam/pU/pEqn.H | 2 - .../multiphase/twoPhaseEulerFoam/pUf/pEqn.H | 2 - .../fvMatrices/fvMatrix/fvMatrix.C | 25 +++-------- 9 files changed, 93 insertions(+), 70 deletions(-) diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H index b1eace14f..dd949f5da 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H @@ -31,27 +31,43 @@ if (pimple.transonic()) { + #include "rhofs.H" + surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi); surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi); p_rghEqnComp1 = - fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1) - + correction - ( - psi1*fvm::ddt(p_rgh) - + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) + pos(alpha1) + *( + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f) + - (fvOptions(alpha1, mixture.thermo1().rho())&rho1) + )/rho1 + - fvc::ddt(alpha1) - fvc::div(alphaPhi1) + + (alpha1/rho1) + *correction + ( + psi1*fvm::ddt(p_rgh) + + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) + ) ); - deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr()); p_rghEqnComp1.ref().relax(); p_rghEqnComp2 = - fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2) - + correction - ( - psi2*fvm::ddt(p_rgh) - + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) - ); - deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr()); + pos(alpha2) + *( + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f) + - (fvOptions(alpha2, mixture.thermo2().rho())&rho2) + )/rho2 + - fvc::ddt(alpha2) - fvc::div(alphaPhi2) + + (alpha2/rho2) + *correction + ( + psi2*fvm::ddt(p_rgh) + + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) + ) + ); p_rghEqnComp2.ref().relax(); } else diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H index 27196f0ee..7867727f8 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H @@ -28,27 +28,43 @@ if (pimple.transonic()) { + #include "rhofs.H" + surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi); surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi); p_rghEqnComp1 = - fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1) - + correction - ( - psi1*fvm::ddt(p_rgh) - + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) + pos(alpha1) + *( + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f) + - (fvOptions(alpha1, mixture.thermo1().rho())&rho1) + )/rho1 + - fvc::ddt(alpha1) - fvc::div(alphaPhi1) + + (alpha1/rho1) + *correction + ( + psi1*fvm::ddt(p_rgh) + + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) + ) ); - deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr()); p_rghEqnComp1.ref().relax(); p_rghEqnComp2 = - fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2) - + correction - ( - psi2*fvm::ddt(p_rgh) - + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) - ); - deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr()); + pos(alpha2) + *( + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f) + - (fvOptions(alpha2, mixture.thermo2().rho())&rho2) + )/rho2 + - fvc::ddt(alpha2) - fvc::div(alphaPhi2) + + (alpha2/rho2) + *correction + ( + psi2*fvm::ddt(p_rgh) + + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) + ) + ); p_rghEqnComp2.ref().relax(); } else diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H index 54bbcc5e4..e7f00b91a 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H @@ -28,27 +28,43 @@ if (pimple.transonic()) { + #include "rhofs.H" + surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi); surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi); p_rghEqnComp1 = - fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1) - + correction - ( - psi1*fvm::ddt(p_rgh) - + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) + pos(alpha1) + *( + ( + fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f) + - (fvOptions(alpha1, mixture.thermo1().rho())&rho1) + )/rho1 + - fvc::ddt(alpha1) - fvc::div(alphaPhi1) + + (alpha1/rho1) + *correction + ( + psi1*fvm::ddt(p_rgh) + + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) + ) ); - deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr()); p_rghEqnComp1.ref().relax(); p_rghEqnComp2 = - fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2) - + correction - ( - psi2*fvm::ddt(p_rgh) - + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) - ); - deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr()); + pos(alpha2) + *( + ( + fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f) + - (fvOptions(alpha2, mixture.thermo2().rho())&rho2) + )/rho2 + - fvc::ddt(alpha2) - fvc::div(alphaPhi2) + + (alpha2/rho2) + *correction + ( + psi2*fvm::ddt(p_rgh) + + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) + ) + ); p_rghEqnComp2.ref().relax(); } else diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H index a164b73f2..c6ec21aae 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H @@ -331,10 +331,6 @@ while (pimple.correct()) ).ptr() ); - deleteDemandDrivenData - ( - pEqnComps[phasei].faceFluxCorrectionPtr() - ); pEqnComps[phasei].relax(); } else diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H index c30d4f146..93c6e3b58 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H @@ -255,7 +255,6 @@ while (pimple.correct()) ) ); - deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr()); pEqnComp1.ref().relax(); } else @@ -297,7 +296,6 @@ while (pimple.correct()) ) ); - deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr()); pEqnComp2.ref().relax(); } else diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H index d2202eb89..b0c8e46b1 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H @@ -241,7 +241,6 @@ while (pimple.correct()) ) ); - deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr()); pEqnComp1.ref().relax(); } else @@ -283,7 +282,6 @@ while (pimple.correct()) ) ); - deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr()); pEqnComp2.ref().relax(); } else diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H index eabc3961d..71d75bd26 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H @@ -262,7 +262,6 @@ while (pimple.correct()) + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) ) ); - deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr()); pEqnComp1.ref().relax(); pEqnComp2 = @@ -278,7 +277,6 @@ while (pimple.correct()) + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) ) ); - deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr()); pEqnComp2.ref().relax(); } else diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H index f8388b2d4..0e1b3b7ff 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H @@ -243,7 +243,6 @@ while (pimple.correct()) + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh) ) ); - deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr()); pEqnComp1.ref().relax(); pEqnComp2 = @@ -259,7 +258,6 @@ while (pimple.correct()) + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh) ) ); - deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr()); pEqnComp2.ref().relax(); } else diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C index 111aecfb7..7cf0ab065 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C @@ -1384,14 +1384,9 @@ Foam::tmp> Foam::correction { tmp> tAcorr = A - (A & A.psi()); - if - ( - (A.hasUpper() || A.hasLower()) - && A.psi().mesh().fluxRequired(A.psi().name()) - ) - { - tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr(); - } + // Delete the faceFluxCorrection from the correction matrix + // as it does not have a clear meaning or purpose + deleteDemandDrivenData(tAcorr.ref().faceFluxCorrectionPtr()); return tAcorr; } @@ -1405,17 +1400,9 @@ Foam::tmp> Foam::correction { tmp> tAcorr = tA - (tA() & tA().psi()); - // Note the matrix coefficients are still that of matrix A - const fvMatrix& A = tAcorr(); - - if - ( - (A.hasUpper() || A.hasLower()) - && A.psi().mesh().fluxRequired(A.psi().name()) - ) - { - tAcorr.ref().faceFluxCorrectionPtr() = (-A.flux()).ptr(); - } + // Delete the faceFluxCorrection from the correction matrix + // as it does not have a clear meaning or purpose + deleteDemandDrivenData(tAcorr.ref().faceFluxCorrectionPtr()); return tAcorr; }