InterFoam viscosity problem / possible bug
 

The InterFoam solver seems to be giving the wrong answer when the viscosity is large. To demonstrate this I generated a test case called waterDrop which I uploaded as a zip file. It is essentially identical to the damBreak tutorial except that the liquid is now in a rectangular blob suspended above the ground by about 0.218 meters. At t=0 it is dropped and falls down. The expected fall time can be obtained from the usual physics formula: t = sqrt(2*h/g), which in this case is about 0.21 seconds. It works as expected for low viscosity, but surprisingly the fall time is found to increase when the viscosity is increased to a large value. As uploaded, the viscosity of the "water" has been increased from 1e-6 to 1e0 and the fall time has increased to about 0.78 seconds. When further increased to 1e2 the water actually floats upwards and disappears out the atmosphere patch! (The viscosity nu for phase1 is set in the file constant/transportProperties. To run and generate log file, enter the commands: blockMesh, checkMesh, setFields, interFoam | tee log1) I have tried adjusting most of the parameters in the fvSolution file with little success. Setting momentumPredictor to yes helps somewhat as does increasing nCorrectors but drastically increases the run time and doesn't completely fix the problem. Decreasing Courant number in the controlDict file also has some beneficial effect. I'm hoping that someone will know how to solve this problem, or else can confirm my suspicion that this is a symptom of some kind of bug in the software.
Thanks - Paul

Hi, I can reproduce this on OpenFOAM 1.7.x too. Did you report it as a bug?

Link to the (closed) bug report: http://openfoam.com/mantisbt/view.php?id=156

Indeed, increasing the viscosity ratio so much makes the whole idea of VOF, based on a mixture approach, questionable...

From a numerical point of view, implementing the PIMPLE procedure to allow under-relaxation is straightforward, and I am not so sure I agree it is much more expensive than performing hundreds PISO correctors. In multiphase codes you often end up needing an unsteady SIMPLE or a PIMPLE procedure to couple the different equations.

Best,

The key point of course is the large viscosity ratio is generating the problem. This can be countered by simultaneously increasing the viscosity of the light phase. Depending upon the nature of the problem, this may not introduce terribly significant errors.

Yes, but with such high viscosity ratio, it would be better to adopt models where each phase has its own equations instead of relying on a mixture approach which is valid when the phase properties are similar.

A couple of comments:

The severity of the problem seems to be tied to the value of the "water" viscosity, rather than the viscosity ratio. I ran a variety of cases which indicate this. For example, with the "water" viscosity set to 0.1, there is a moderate effect where the fall time is increased to 0.33 seconds from the expected value of about 0.20. Keeping 0.1 for the "water" and increasing the viscosity of the "air" by a factor of 100 from 1.48e-5 to 1.48e-3 there is only a barely perceptible change in the fall time, increasing to about 0.335. The viscosity of the "air" can also be decreased all the way to zero with negligible effect. When increasing the density of the "air" from 1 to 10 the effect is also quite small.

I am also interested using non-Newtonian viscosity such as the Herschel-Bulkley model. When the shear rate decreases, such a fluid can transition from moderate viscosity to nearly infinite viscosity. It would be very desirable for the simulation to perform reliably in both regimes. Incidentally, I recently reported a bug in that model http://openfoam.com/mantisbt/view.php?id=153.

Paul

Effective "g" value is scaled down
 

For moderately high viscosity (up to nu=1.0 at least) the block of "water" does not behave erratically, but rather it moves in the downwards direction as expected from gravity. The only problem is that it accelerates more slowly as though the value of "g" was scaled down by some factor that depends on the water's viscosity. This remains approximately true even when "g" is drastically increased. For example, in the nu=1.0 case, if g is increased by a factor of 100 from 9.81 to 981, the fall time of the block is reduced by roughly a factor of 10 from 0.78 to about 0.07 sec, which is as expected since the fall time should be proportional to the inverse square root of g. The effective value of g based on the fall time t is roughly 0.073 times the actual value of g: [geff/g = 2*h/(g*t^2)]. There also appears to be some size dependence to the effect. If the block dimensions are reduced by a factor of 2, the fall time increases from 0.76 to 0.99 sec.

The consistent nature of the error in these cases suggest that there is something about the algorithm or the way it is implemented that is resulting in the reduced effectiveness of gravity in the problem. If anyone has any comments or suggestions of how to possibly fix the problem without increasing the run time by orders of magnitude I would greatly appreciate it.

Thanks, Paul

Hi Paul,

I believe the problem is of numerical origin. Did you check the evolution of the residuals and the convergence of your solution?

Alberto

Reducing all tolerances in fvSolution by 10x has almost no effect. Here is a typical section of output from the log file:
__________________________________________________ ___
Courant Number mean: 0.0359508 max: 0.397045
Interface Courant Number mean: 0.00321944 max: 0.246874
deltaT = 0.00111111
Time = 0.0411111

MULES: Solving for alpha1
Liquid phase volume fraction = 0.01716 Min(alpha1) = -8.72931e-40 Max(alpha1) = 0.999951
MULES: Solving for alpha1
Liquid phase volume fraction = 0.01716 Min(alpha1) = -4.82315e-40 Max(alpha1) = 0.999951
DICPCG: Solving for p_rgh, Initial residual = 0.00259989, Final residual = 1.1547e-05, No Iterations 22
time step continuity errors : sum local = 2.37286e-05, global = 6.0088e-08, cumulative = -1.85707e-05
DICPCG: Solving for p_rgh, Initial residual = 6.39632e-05, Final residual = 2.712e-07, No Iterations 36
time step continuity errors : sum local = 5.47611e-07, global = -6.02456e-08, cumulative = -1.86309e-05
DICPCG: Solving for p_rgh, Initial residual = 2.33225e-05, Final residual = 7.87965e-09, No Iterations 63
time step continuity errors : sum local = 1.59131e-08, global = 1.33036e-10, cumulative = -1.86308e-05
ExecutionTime = 4.05 s ClockTime = 5 s
__________________________________________________ _____

Forcing a greatly smaller time step by say 10x or 100x does lead to improvement, but this is undesirable since it results in an increase in run time by a similar factor.

Again, it seems very strange that the motion of the block of water seems as though it is responding to a gravitational force of greatly reduced magnitude. It seems to accelerate normally, with no sign of reaching a terminal velocity as it would if the sluggish response was a result of the air viscosity.

Paul

You should plot the initial residuals of the equations and study their trend to see if the solution of the coupled system converges.

The fact that improvements appear reducing the time step seem to confirm what Henry wrote in the bugzilla. You might need an unsteady SIMPLE/PIMPLE algorithm with under-relaxation of the solution. This would allow larger time steps, since you will perform outer correctors.

Best,
Alberto

momentumPredictor eliminates gravity problem
 

The "momentum predictor" step is normally turned off in interFoam. I'm not sure of the reasoning behind this, but I have found that turning it on leads to a very dramatic improvement in the dynamics for high viscosity cases (much more than I had originally thought). If we consider a 10% error in the fall rate of the block of water to be the acceptable limit, then turning on the momentum predictor increases the maximum usable viscosity for the water from about 0.02 to at least 4.0, i.e. at least a factor of 200 times! I say "at least" because above 4.0 other instabilities take over that cause the block of water to move about erratically, which is quite different from the "reduced gravity" effect that is observed below this level.

I checked again on the suggested possibility of a "viscosity ratio" problem. For nu=4.0 the viscosity of the air can be reduced to zero with no significant effect on the motion of the block of water. For nu=10.0, where the behavior is erratic, increasing the viscosity of the air by a factor of 100 does change the behavior somewhat, but it remains very erratic, i.e. reducing the viscosity ratio by 100x has no beneficial effect.

For low viscosity, the momentum predictor seems to effect the run time by only about 10%. So it might be appropriate to have it on by default and set to "yes" in the file fvSolution for examples such as damBreak. It might also be useful to include some discussion of this in the user manual when discussing the damBreak tutorial.

Any comments or suggestions will be greatly appreciated,
Thanks, Paul

Thank you for your updates.

Yes, solving the momentum predictor is expected to play an important role when viscous stresses are important or dominant.

What kind of instability? Altered interface (high velocities there?), or complete crash?

Alberto,

The block of water has a tendency in these very high viscosity cases (nu>4) to rotate while moving, sometimes going sideways or even upwards and not necessarily in a straight line. Sometimes it will collide with one of the sides or go out the atmosphere patch, usually within about 0.4 seconds or so. It usually remains intact and block shaped but may show some deformations in the surface.

You mentioned that it is known that the momentum predictor is important when the viscosity is high (or low Re) -- do you know of a reference that discusses this? (I am curious to see what other people have said about it)

Thanks, Paul

I would say that what you observe (the deformation and non-vertical movement) is a consequence of the high velocities you should observe at the interface.

There was a discussion on this forum about that. The summary is simply the following: if you do not solve for the momentum predictor, you estimate U with a Jacobi sweep. This approximately accounts for the viscous stresses, which is acceptable if they are not dominant. In cases like yours, where viscous stresses are very strong, solving the momentum equation leads to a more accurate prediction of U before solving the pressure equation, and it should have stabilizing effects of the iterative procedure.

If you increase the viscosity as much as you do, you are basically solving for a blob that tends to behave like a solid in one block of cells, and something fluid on the other side of the interface. This corresponds to a sharp viscosity gradient, which will destabilise your solution.

Did you try to put a limiter on the gradients of the viscous term and use a completely limited scheme? This is very common in commercial codes.

Best,

High-viscosity problem with interFoam
 

Hi,

As I seem to face a similar problem like Paul, I would like to contribute to this thread. I'm currently using a modified version of compressibleInterFoam to model time-dependent physical foam expansion (mostly polystyrene and silicone oil with carbon dioxide as blowing gas). The viscosity of the foam phase can reach values well above 10 Pas.

Now, as Paul already stated, the problem seems to be mostly tied to the maximum viscosity value and not the gradient at the interface. I ran several simple simulations with a vertical rectangular channel containing a certain amount of polymer that is sitting at the bottom. No density change whatsoever was imposed, so one would expect the interface to remain immobile. However, this is only true for low viscosity values. For values above 1 Pas however, the polymer immediately starts to move and shift upwards in the center due to high velocities at the interface. Raising the air viscosity to values similar to the polymer has very little effect, so the viscosity gradient does not seem to be the only problem.

The case runs perfectly with calm interface if gravity is set to a very low value (eg -0.01 m^2/s), so I would assume that the core of the problem seems to be some kind of (numerical?) interference of the gravitational and viscous forces for high viscosity values.

I would be very grateful If anyone could shed some more light on this issue, as for now my only solution to this problem is to drastically reduce gravity in my calculations, which of course is not always appropriate.

Florian

Hi,

I tried to reproduce your problem without success. I took the depthCharge2D case, removed the sphere from setFields, so to fill the bottom of the container with phase1, and I set nu = 0.01 (10 Pa s / (1000 kg/m^3)), since the density of phase 1 is 1000 kg/m^3. Gravity is the correct value along y.

I ran the case, and the field of alpha1 seems fine (picture attached).

Maybe you could attach a small test case reproducing the problem.

Best,

Hi Alberto,

thanks for your fast reply. You are right, it seems to work with the official release of compressibleInterFoam. The solver I'm currently using was initially developed by a former colleague of me and is based on an early version of compressibleInterFoam. One thing that just caught my eye is that the solver expects dynamic viscosity values from transportproperties.dict, and uses an altered twoPhaseMixture-class to calculate muf in Ueqn.h as follows:

"muf",
alpha1f*fvc::interpolate(muModel1_->mu()) + (scalar(1) - alpha1f)*fvc::interpolate(muModel2_->mu())

"muf" is called in Ueqn.h right before the velocity equation. Other than that, I think twoPhaseMixture is only used to update the viscosity for each time step (.correct).
Maybe this is the root of my interface instability, although I don't see yet where problems could arise by specifying mu instead of nu.

To be sure, I will try and change parts of the code to be consistent with the official compressibleInterFoam and see if my problems persist. Hopefully I can provide some more information in the future.

Florian

Problem solved
 

Hi,
I finally managed to locate the error in my calculations. Thankfully it's not related to the solver itself. but rather to the discretization and solution procedure. I changed parts of my case to match fvSchemes and fvSolutions from the depthCharge-tutorial, and now it's working like a charm. Thank you very much Alberto for pointing me in the right direction.

If anyone else is having trouble dealing with multiphase simulations involving high viscosity gradients: one measure that tremendously improved my velocity field was switching off the momentum predictor. With predictor steps, I get gradually increasing velocities at the interface that eventually get out of hand.

Thanks for the hint with the momentum predictor. Has anyone experience with perssure-dependent viscosity models in interFoam? I get high pressure fluctuations / viscosity fluctuations depending on whether I deal with p or with p_rgh, and I wonder wether I have to switch to compressible two phase flow to deal with it.