diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H index 699d789adb..2cadf01317 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H @@ -37,14 +37,17 @@ if (mesh.changing()) pcorrTypes ); - // dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0); surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); + #ifndef divUCorr + #define divUCorr + #endif + while (pimple.correctNonOrthogonal()) { fvScalarMatrix pcorrEqn ( - fvm::laplacian(rAUf, pcorr) == fvc::div(phi) + fvm::laplacian(rAUf, pcorr) == fvc::div(phi) divUCorr ); pcorrEqn.setReference(pRefCell, pRefValue); @@ -56,5 +59,7 @@ if (mesh.changing()) } } + #undef divUCorr + #include "continuityErrs.H" } diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H index fb38fed37e..1172cf0707 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H @@ -82,7 +82,16 @@ 0 ); - tphiAlpha() += tphiAlphaCorr(); + // Under-relax the correction for all but the 1st corrector + if (aCorr == 0) + { + tphiAlpha() += tphiAlphaCorr(); + } + else + { + alpha1 = 0.5*alpha1 + 0.5*alpha10; + tphiAlpha() += 0.5*tphiAlphaCorr(); + } } else { diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C index 8bd1984539..e5ae566f6e 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C @@ -82,7 +82,6 @@ int main(int argc, char *argv[]) dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0) ); - #include "../interFoam/interDyMFoam/correctPhi.H" #include "createUf.H" #include "CourantNo.H" #include "setInitialDeltaT.H" @@ -106,6 +105,9 @@ int main(int argc, char *argv[]) { if (pimple.firstIter() || moveMeshOuterCorrectors) { + // Store divU from the previous mesh for the correctPhi + volScalarField divU(fvc::div(fvc::absolute(phi, U))); + scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime(); mesh.update(); @@ -125,6 +127,7 @@ int main(int argc, char *argv[]) // Calculate absolute flux from the mapped surface velocity phi = mesh.Sf() & Uf; + #define divUCorr -divU #include "../interFoam/interDyMFoam/correctPhi.H" // Make the flux relative to the mesh motion diff --git a/tutorials/multiphase/interPhaseChangeDyMFoam/propeller/system/fvSchemes b/tutorials/multiphase/interPhaseChangeDyMFoam/propeller/system/fvSchemes index b14f07a132..663657c047 100644 --- a/tutorials/multiphase/interPhaseChangeDyMFoam/propeller/system/fvSchemes +++ b/tutorials/multiphase/interPhaseChangeDyMFoam/propeller/system/fvSchemes @@ -32,7 +32,6 @@ divSchemes div(phi,alpha) Gauss vanLeer; div(phirb,alpha) Gauss linear; - UD Gauss upwind; div(rhoPhi,U) Gauss linearUpwind grad(U); div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind;