diff --git a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
index 13e1feb546e1e9d9efdf2965b99b5f8ff5d8b193..970e97157f4c3875cc807031d5daf5b1d9be764c 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/TEqn.H
@@ -5,8 +5,8 @@
+ fvm::div(rhoPhi, T)
- fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
+ (
- fvc::div(fvc::absolute(phi, U), p)
- + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
+ divU*p
+ + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
)
*(
alpha1/mixture.thermo1().Cv()
@@ -20,4 +20,4 @@
mixture.correct();
Info<< "min(T) " << min(T).value() << endl;
-}
+}
\ No newline at end of file
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 7cbbd26b141ce1df7f8ddc552ce3d80a5d900999..daef88206b26a52cdbe563b3d0253709614fffab 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -70,6 +70,21 @@ int main(int argc, char *argv[])
turbulence->validate();
+
+ volScalarField rAU
+ (
+ IOobject
+ (
+ "rAU",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("rAU", dimTime/rho.dimensions(), 1.0)
+ );
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
@@ -77,45 +92,17 @@ int main(int argc, char *argv[])
{
#include "readControls.H"
- {
- // Store divU from the previous mesh so that it can be mapped
- // and used in correctPhi to ensure the corrected phi has the
- // same divergence
- volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
-
- #include "CourantNo.H"
- #include "setDeltaT.H"
-
- runTime++;
-
- Info<< "Time = " << runTime.timeName() << nl << endl;
-
- scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
-
- // Do any mesh changes
- mesh.update();
+ // Store divU from the previous mesh so that it can be mapped
+ // and used in correctPhi to ensure the corrected phi has the
+ // same divergence
+ volScalarField divU("divU0",fvc::div(fvc::absolute(phi, U)));
- if (mesh.changing())
- {
- Info<< "Execution time for mesh.update() = "
- << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
- << " s" << endl;
-
- gh = (g & mesh.C()) - ghRef;
- ghf = (g & mesh.Cf()) - ghRef;
- }
-
- if (mesh.changing() && correctPhi)
- {
- // Calculate absolute flux from the mapped surface velocity
- phi = mesh.Sf() & Uf;
+ #include "CourantNo.H"
+ #include "setDeltaT.H"
- #include "correctPhi.H"
+ runTime++;
- // Make the fluxes relative to the mesh motion
- fvc::makeRelative(phi, U);
- }
- }
+ Info<< "Time = " << runTime.timeName() << nl << endl;
if (mesh.changing() && checkMeshCourantNo)
{
@@ -127,6 +114,46 @@ int main(int argc, char *argv[])
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
+ if (pimple.firstIter())
+ {
+ scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+ mesh.update();
+
+ if (mesh.changing())
+ {
+ Info<< "Execution time for mesh.update() = "
+ << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+ << " s" << endl;
+
+ gh = (g & mesh.C()) - ghRef;
+ ghf = (g & mesh.Cf()) - ghRef;
+ }
+
+ if ((correctPhi && mesh.changing()) || mesh.topoChanging())
+ {
+ // Calculate absolute flux from the mapped surface velocity
+ // SAF: temporary fix until mapped Uf is assessed
+ Uf = fvc::interpolate(U);
+
+ phi = mesh.Sf() & Uf;
+
+ #include "correctPhi.H"
+
+ fvc::makeRelative(phi, U);
+
+ mixture.correct();
+
+ mesh.topoChanging(false);
+ }
+
+ if (mesh.changing() && checkMeshCourantNo)
+ {
+ #include "meshCourantNo.H"
+ }
+ }
+
+
#include "alphaEqnsSubCycle.H"
// correct interface on first PIMPLE corrector
@@ -145,6 +172,11 @@ int main(int argc, char *argv[])
{
#include "pEqn.H"
}
+
+ if (pimple.turbCorr())
+ {
+ turbulence->correct();
+ }
}
rho = alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index 859fc9cc4e8aa240f7d280467dfa42017c67449a..34e39faddef662ed991bc30b8981e55c0704342e 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -87,15 +87,6 @@
if (pimple.finalNonOrthogonalIter())
{
- p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
- p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
-
- dgdt =
- (
- pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
- - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
- );
-
phi = phiHbyA + p_rghEqnIncomp.flux();
U = HbyA
@@ -116,8 +107,19 @@
rho = alpha1*rho1 + alpha2*rho2;
+ p = max(p_rgh + rho*gh, pMin);
+ p_rgh = p - rho*gh;
+
+ dgdt =
+ (
+ pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
+ - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
+ );
+
K = 0.5*magSqr(U);
Info<< "max(U) " << max(mag(U)).value() << endl;
Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
}
+
+
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index 2a113880dee0d5f7b783e8cdc0446075b7c03e2e..e4d413964d0f13b6fa7e0d9981d2978681c916ae 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -93,6 +93,7 @@ int main(int argc, char *argv[])
solve(fvm::ddt(rho) + fvc::div(rhoPhi));
#include "UEqn.H"
+ volScalarField divU(fvc::div(fvc::absolute(phi, U)));
#include "TEqn.H"
// --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index a4a3fb2e2914150398b4ea3f6da31ee5e58e8b4e..f24b72eb684413803a0345fb75f692827803784d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -84,15 +84,6 @@
if (pimple.finalNonOrthogonalIter())
{
- p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
- p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
-
- dgdt =
- (
- pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
- - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
- );
-
phi = phiHbyA + p_rghEqnIncomp.flux();
U = HbyA
@@ -101,14 +92,21 @@
}
}
- // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-
// Update densities from change in p_rgh
rho1 += psi1*(p_rgh - p_rgh_0);
rho2 += psi2*(p_rgh - p_rgh_0);
rho = alpha1*rho1 + alpha2*rho2;
+ p = max(p_rgh + rho*gh, pMin);
+ p_rgh = p - rho*gh;
+
+ dgdt =
+ (
+ pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
+ - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
+ );
+
K = 0.5*magSqr(U);
Info<< "max(U) " << max(mag(U)).value() << endl;
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 340f2c948edaf07f481c05e3203dd58b65b85389..401ece3fde71622172420515756259decca179e9 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -133,9 +133,12 @@ int main(int argc, char *argv[])
ghf = (g & mesh.Cf()) - ghRef;
}
- if (mesh.changing() && correctPhi)
+ if ((mesh.changing() && correctPhi) || mesh.topoChanging())
{
// Calculate absolute flux from the mapped surface velocity
+ // SAF: temporary fix until mapped Uf is assessed
+ Uf = fvc::interpolate(U);
+
phi = mesh.Sf() & Uf;
#include "correctPhi.H"
@@ -144,6 +147,8 @@ int main(int argc, char *argv[])
fvc::makeRelative(phi, U);
mixture.correct();
+
+ mesh.topoChanging(false);
}
if (mesh.changing() && checkMeshCourantNo)