Mass conservation problem with compressibleInterDyMFoam
 

Hi Foamers,

I am working on a simple 2D case that flow passes by a channel with one side of the channel wall being concave. Flow goes in from the bottom and goes out from the top. The concave wall is fixed and the left flat wall is moving towards the concave wall.


I have specified:
Pressure: outlet fixed at 1 bar, zeroGradient for others;
Velocity: flowRateInletVelocity 1e-05 at the inlet, pressureInletOutletVelocity at the outlet, fixed value 0 at the concave wall, and movingWallVelocity at the moving wall;

I used compressibleInterFoam for my injection simulation. Therefore I do my moving mesh test based on this solver. Here the fluid is water (although I still call it gasoline...), and I have switched off the compressibility by setting psi = 0.

First test no mesh motion, pointMotionUx = 0, mass conservation is fine. Here is an extract of the log file
Courant Number mean: 0.000352001 max: 0.00046794
deltaT = 2e-05
Time = 0.0067

DICPCG: Solving for cellMotionUx, Initial residual = 0, Final residual = 0, No Iterations 0
Execution time for mesh.update() = 0.01 s
time step continuity errors : sum local = 1.14092e-10, global = 1.14034e-10, cumulative = -2.17556e-06
GAMGPCG: Solving for pcorr, Initial residual = 1, Final residual = 8.96484e-06, No Iterations 5
time step continuity errors : sum local = 1.4302e-15, global = -1.4698e-16, cumulative = -2.17556e-06
--> FOAM Warning :
From function cubeRootVolDelta::calcDelta()
in file cubeRootVolDelta/cubeRootVolDelta.C at line 52
Case is 2D, LES is not strictly applicable

PIMPLE: iteration 1
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -1.28786e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 1.2972e-05, Final residual = 1.71303e-14, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 0.00117604, Final residual = 2.75328e-06, No Iterations 2
max(U) 0.0499816
min(p_rgh) 99999.4
PIMPLE: iteration 2
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -1.28786e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 7.40712e-06, Final residual = 5.33731e-15, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 1.9535e-05, Final residual = 6.37423e-08, No Iterations 3
max(U) 0.0499816
min(p_rgh) 99999.4
PIMPLE: iteration 3
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -1.28786e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 7.26948e-06, Final residual = 5.15672e-15, No Iterations 1
GAMGPCG: Solving for p_rgh, Initial residual = 5.26857e-06, Final residual = 2.05796e-08, No Iterations 2
max(U) 0.0499816
min(p_rgh) 99999.4
ExecutionTime = 13.05 s


Courant Number mean: 0.000352 max: 0.000467942
deltaT = 2e-05
Time = 0.00672

DICPCG: Solving for cellMotionUx, Initial residual = 0, Final residual = 0, No Iterations 0
Execution time for mesh.update() = 0.01 s
time step continuity errors : sum local = 1.13852e-10, global = 1.13793e-10, cumulative = -2.17545e-06
GAMGPCG: Solving for pcorr, Initial residual = 1, Final residual = 8.96323e-06, No Iterations 5
time step continuity errors : sum local = 1.42864e-15, global = -1.46837e-16, cumulative = -2.17545e-06
--> FOAM Warning :
From function cubeRootVolDelta::calcDelta()
in file cubeRootVolDelta/cubeRootVolDelta.C at line 52
Case is 2D, LES is not strictly applicable

PIMPLE: iteration 1
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -1.28786e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 1.32896e-05, Final residual = 1.7448e-14, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 0.00117269, Final residual = 2.75305e-06, No Iterations 2
max(U) 0.0499818
min(p_rgh) 99999.4
PIMPLE: iteration 2
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -1.28786e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 7.50453e-06, Final residual = 5.21501e-15, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 1.94767e-05, Final residual = 6.45389e-08, No Iterations 3
max(U) 0.0499818
min(p_rgh) 99999.4
PIMPLE: iteration 3
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -1.28786e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 7.29273e-06, Final residual = 5.1861e-15, No Iterations 1
GAMGPCG: Solving for p_rgh, Initial residual = 5.25134e-06, Final residual = 2.06509e-08, No Iterations 2
max(U) 0.0499818
min(p_rgh) 99999.4
ExecutionTime = 13.08 s


Second test, pointMotionUx = 0.2, mass is not conserved.
Courant Number mean: 0.0128567 max: 0.0207858
deltaT = 2e-05
Time = 0.0067

DICPCG: Solving for cellMotionUx, Initial residual = 4.419e-07, Final residual = 9.74887e-09, No Iterations 12
Execution time for mesh.update() = 0.01 s
time step continuity errors : sum local = 0.000138757, global = 0.000138757, cumulative = 0.0905193
GAMGPCG: Solving for pcorr, Initial residual = 1, Final residual = 9.55594e-06, No Iterations 7
time step continuity errors : sum local = 0.000138771, global = 0.000138771, cumulative = 0.0906581
--> FOAM Warning :
From function cubeRootVolDelta::calcDelta()
in file cubeRootVolDelta/cubeRootVolDelta.C at line 52
Case is 2D, LES is not strictly applicable

PIMPLE: iteration 1
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -1.9762e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 0.00238355, Final residual = 2.36096e-12, No Iterations 2
GAMG: Solving for p_rgh, Initial residual = 0.000295135, Final residual = 1.29158e-06, No Iterations 2
max(U) 2.03828
min(p_rgh) 99515.3
PIMPLE: iteration 2
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -1.95399e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 2.86226e-07, Final residual = 4.05804e-11, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 7.16766e-05, Final residual = 5.46045e-07, No Iterations 3
max(U) 2.03828
min(p_rgh) 99514.4
PIMPLE: iteration 3
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -1.93179e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 9.52316e-09, Final residual = 9.52316e-09, No Iterations 0
GAMGPCG: Solving for p_rgh, Initial residual = 1.78146e-05, Final residual = 8.30045e-08, No Iterations 2
max(U) 2.03828
min(p_rgh) 99514.7
ExecutionTime = 14.46 s


Courant Number mean: 0.0128578 max: 0.0207894
deltaT = 2e-05
Time = 0.00672

DICPCG: Solving for cellMotionUx, Initial residual = 4.42075e-07, Final residual = 9.59846e-09, No Iterations 12
Execution time for mesh.update() = 0.01 s
time step continuity errors : sum local = 0.000138776, global = 0.000138776, cumulative = 0.0907968
GAMGPCG: Solving for pcorr, Initial residual = 1, Final residual = 8.71665e-06, No Iterations 7
time step continuity errors : sum local = 0.00013879, global = 0.00013879, cumulative = 0.0909356
--> FOAM Warning :
From function cubeRootVolDelta::calcDelta()
in file cubeRootVolDelta/cubeRootVolDelta.C at line 52
Case is 2D, LES is not strictly applicable

PIMPLE: iteration 1
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -2.13163e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 0.00237817, Final residual = 2.35453e-12, No Iterations 2
GAMG: Solving for p_rgh, Initial residual = 0.000361452, Final residual = 1.39167e-06, No Iterations 2
max(U) 2.03827
min(p_rgh) 99514.7
PIMPLE: iteration 2
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -2.10942e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 2.84647e-07, Final residual = 4.18919e-11, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 8.95334e-05, Final residual = 6.61969e-07, No Iterations 3
max(U) 2.03827
min(p_rgh) 99513.9
PIMPLE: iteration 3
MULES: Solving for alphaGasoline
Liquid phase volume fraction = 1 Min(alpha1) = 1 Min(alpha2) = -2.10942e-14
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for T, Initial residual = 7.85789e-09, Final residual = 7.85789e-09, No Iterations 0
GAMGPCG: Solving for p_rgh, Initial residual = 2.31585e-05, Final residual = 7.22024e-08, No Iterations 2
max(U) 2.03827
min(p_rgh) 99514.2
ExecutionTime = 14.5 s


When I run this setup with pimpleDyMFoam, mass is conserved too.

Is there a bug here? Your comments and help are mostly welcome!


Cheers,
Leslie

In fact it's the same problem with interDyMFoam. Is it a general problem of the multiphase solvers?

The bug locates in pEqn.H

surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho,U)
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
);
phi = phiHbyA;

Henry Weller mentioned on your bug report[0001110]: "Resolved by commit f3639a9c9d2baed720f2890e819f17af92a03a23"

What does it mean? Is the problem fixed in the last OF version?

I ran into the same problem by using openfoam extend 3.1.1. I tried to make the implementation of the mesh motion. The mass in the system is now constant so that is fine, but now i have instabilities due to high internal residuals. As mentioned that implementation is not sufficient, but is there an easy possibility to fix that problem in the pEqn in openfoam extend. Out of the architecture of openfoam 2.4. i fail to implement the changes in my solver in openfoam 3.1.1.