diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C index 8e3500c9d6ad7f80d38c7e360c2dabf5ec2f7ae5..1face6ddd57a9ad664f9a7811403fc8cc49d9a68 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C @@ -58,6 +58,9 @@ int main(int argc, char *argv[]) #include "CourantNo.H" #include "setInitialDeltaT.H" + surfaceScalarField phiAbs("phiAbs", phi); + fvc::makeAbsolute(phiAbs, U); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -67,9 +70,6 @@ int main(int argc, char *argv[]) #include "alphaCourantNo.H" #include "CourantNo.H" - // Make the fluxes absolute - fvc::makeAbsolute(phi, U); - #include "setDeltaT.H" runTime++; @@ -78,8 +78,18 @@ int main(int argc, char *argv[]) scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime(); - // Do any mesh changes - mesh.update(); + { + // Calculate the relative velocity used to map the relative flux phi + volVectorField Urel("Urel", U); + + if (mesh.moving()) + { + Urel -= fvc::reconstruct(fvc::meshPhi(U)); + } + + // Do any mesh changes + mesh.update(); + } if (mesh.changing()) { @@ -96,9 +106,6 @@ int main(int argc, char *argv[]) #include "correctPhi.H" } - // Make the fluxes relative to the mesh motion - fvc::makeRelative(phi, U); - if (mesh.changing() && checkMeshCourantNo) { #include "meshCourantNo.H" diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H index 84c1fe45b6a40666d5a7a1f184da5a1700eab522..4b90d9949929bc3a5b7d761acc8294604c9d4b51 100644 --- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H @@ -3,16 +3,19 @@ surfaceScalarField rAUf(fvc::interpolate(rAU)); U = rAU*UEqn.H(); - surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf())); + + phiAbs = + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phiAbs); if (p_rgh.needReference()) { - fvc::makeRelative(phiU, U); - adjustPhi(phiU, U, p_rgh); - fvc::makeAbsolute(phiU, U); + fvc::makeRelative(phiAbs, U); + adjustPhi(phiAbs, U, p_rgh); + fvc::makeAbsolute(phiAbs, U); } - phi = phiU + + phi = phiAbs + ( fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - ghf*fvc::snGrad(rho) @@ -38,11 +41,13 @@ } } - U += rAU*fvc::reconstruct((phi - phiU)/rAUf); + U += rAU*fvc::reconstruct((phi - phiAbs)/rAUf); U.correctBoundaryConditions(); #include "continuityErrs.H" + phiAbs = phi; + // Make the fluxes relative to the mesh motion fvc::makeRelative(phi, U);