Artificial high velocities at the interface using interFoam
 

Hi all,

I'm simulation flow of water around an object using the interFoam solver and k-Omega SST turbulence model. At the moment, I'm struggling with some artificial high velocities which occur at the interface in the upper air zone. At the beginning of the simulation, some air interface cells have velocities of up to 20 times the normal maximum velocities in my water zone, which really reduce my time steps or let the simulation crash.

For the meshing, I used blockMesh and sHM including layers. checkMesh gives

...
Minumum face area = 3.3007033e-07. Maximum face area = 0.0026526524. Face area magnitudes OK.
Min volume = 4.735839e-10. Max volume = 5.2955194e-05. Total volume = 1.6946638. Cell volumes OK.
Mesh non-orthogonality Max: 74.576967 average: 6.9868146
*Number of severely non-orthogonal faces: 70.
Non-orthogonality check OK.
<<Writing 70 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 2.4098523 OK.

Mesh OK.


As I'm not interested in what is happening on the air part, I already changed the interFoam solver to ignore the convective term on the air side. This improved the simulation, but it's still not running very well.

Can someone give me a hint on how to modify my solver settings or improve my mesh?
A picture and the settings are attached/given below.


BTW: Is it somehow possible to 'sharpen' the interface by changing the settings? Or is this only possible by refining the mesh at the interface?

Arne


http://img833.imageshack.us/img833/5...mfoto1l.th.png http://img135.imageshack.us/img135/6...mfoto2d.th.png

Hi,

try using a limited scheme for div(rho*phi,U). You are currently using central differences (linear), which is not necessarily stable (is your cell Re < 2)?

Maybe try

div(rho*phi,U) Gauss limitedLinearV 1;

which should be quite accurate, if it works. If you need something more robust:

div(rho*phi,U) Gauss linearUpwindV cellLimited Gauss linear 1;

Best,

Thanks Alberto,

I guess you saved a lot of my time! I tried both schemes, and both simulations are working. The Gauss limitedLinearV 1 seems more reliable and accurate to me.

Nevertheless, artificial velocities at the interface still occure, but are much lower. They are "generated" directly at the simulation start, and never disappear (the latter maybe because the convective term on the air side is ignored?).
Do you have an idea where they are coming from and, what is more important, how I can completely avoid them? Any hints on which direction I could have a closer look?

Initial and BC: For the initial state, the whole water domain has a unique velocity given by setFields. Air has zero velocity. At the inlet, velocity is set as constant at the water part using groovyBC and zero for the air part during the simulation.
For the initial state, I also changed this to a smoother transition for alpha and U giving the interface cells half values (0.5 and half water velocity). This did not solve it as well...

Arne

More accurate, yes. More stable, I am not so sure :D

Try with the unmodified code to see if this still happens.

Best,
Alberto

Hi Alberto,

Although I changed and primarily tried different fvSolution settings for a comparison of PCG and GAMG solver, I also did calculations with modified and unmodified interFoam code. For the 'original' one, these artificial velocities occur as well, in the near-field of the mediums air, water and solid structure (with layers).
I noticed that some small cells in the boundary layer around the structure at the water/air transition have a very high value for k, in an order of 10 times the normal value for cells around the structure.
But I don't exactly know what to do with them...

Another question: Do you also have a suggestion for the 'right' fvSolution settings (PCG, GAMG) in terms of speed and accuracy for interFoam and around 2 million cells?

Arne

You might want to see if it still happens without turbulence model (laminar).

It is more a question of speed than accuracy. I tent to use GAMG for pressure, and CG methods for the rest, which is probably the default setting in the tutorials. If GAMG gives you troubles, use PCG for pressure too: it will take more iterations, but it seems to be a little more robust.

Best,
Alberto

Hello Arne,

I am interested in ignore the convective term on the air side, can you let me known how did you do?, is it improving your simulation? (i mean if you are avoiding very small time steps?

Advanced thanks

Pablo

Pablo, it did not really improve the simulation with regard to smaller time steps and high velocities in the air phase.

Depending on what you want to do and if the air phase is not of great importance for you, you could try modifying the solver (e.g. interFoam, copied and compiled to my_interFoam) and set the velocities in alpha1 (all cels smaller than a value of lets say 0.05) to zero every timestep. For my case, this saved computational time. Nevertheless, in case of surface waves, you have to be careful with regard to wave damping.

Ignoring the convective term would also be done in the solver within the U equation.

Arne

Hello Arne,

Set velocities 0 when alpha1<0.05 is clear for me, it is like a BC in every timestep, but modify the convertive term, i can not see too clear to implement.

We have
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
So we can modify U or rhoPhi, but how??, can you be more explicit?

Advanced thanks

Please have a look at the slides from Eric Paterson: http://www.google.de/url?sa=t&source...bWpkzA&cad=rja

Here, gamma/alpha1 is explicitly included in the equation and therefore the term becomes zero in air phase.

Now it is very clear and easy. Thank you very much.

Pablo

Hi Arne,

My version from interFoam is rotational and inviscid, and i modify with ;

fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ alpha1*fvm::div(rhoPhi, U)

but surprised that i got unstable solution. Any idea? if i need to modify another piece of code?

Pablo

Hello,
I am trying to set the velocities to 0 in the air phase as you have mentioned. could you please share how did you do that? I am new to OpenFoam.
Thanks

You can try something like this:

forAll(U, celli)
{
if (alpha1[celli] < 0.01) {
U[celli] = 0.0;
}

}

Info <<"/////////////////////// update U Air ///////////////////////////////// " << nl;

You can write in a file like mycorrectAir.H, and call from UEqn.H, after convective term.

Pablo

That's what I tried as well, but also got unstable results, depending on the mesh. Therefore I did not use it.

@ cfd_user2011: Concerning setting the velocity in air to zero: You have to modify the (interFoam) solver, like Pablo said. Please have a look at this tutorial http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam to see how a copy of the solver can be modified and compiled.

@ Pablo: I'm encountering that the deltaT significantly drops (by a factor 2) durin the whole simulation (using modified solver), compared to normal interFoam. Do you get the same?

Hi Arne,

I am not using exactly 0 for velocities at the air, i think that is too agressive. I am relaxing only the air phase and i am getting nice results, i mean i am running mesh that before was stopped because unstabilities at the air.

About alpha1*fvm::div(rhoPhi, U) the idea looks right but maybe the mesh , allways i got unstable.

Pablo

Yes, that might be true. Setting the velocity to 0.7 times computed velocity results in a speedup of around 15 percent in my case...

Arne

I'd like to add something to this thread, since I encoutered similar problems for my case. This problems disappeared with reducing the density ratio in my case from 1000 to 10.

Since I don't like the clipping in the air velocity I was able to get better results by applying the argument of Brackbill ( http://www.sciencedirect.com/science...2199919290240Y ) for surface tension, on the additional source term in the Navier-Stokes equations resutling for the reformulation of the pressure.

Basically what I did, was replacing this (in pEqn.H and UEqn.H of interFoam)
By
The argument by Brackbill is summarized as follows: as the interface thickness goes to zero, you can multiply the body force by a function that is 1 at the interface. Since \grad\rho is zero except at the interface and 2alpha1 is 1 there, I think this is justified (please correct me if I'm wrong)

For me it solved the issues with the high air velocities, and this implementation is at least more justified than just clipping U.

Thanks to share Bernhard, i will try to test ASAP.

Pablo

Hi, the solution
is not compiling for me... I got this answer :

It seems that the =0.0 is not working for a vector class... How should I write it better ?

thanks in advance !
Joris

Hi Joris

You could type

However, you ought to be careful as you are threading onto an extremely dangerous path. What you are doing here is equivalent to be removing energy from the system.

I do not have the solution to circumvent the high velocities, but I my case I have found the above approach to yield utterly wrong results.

Best regards,

Niels

Thanks niels for your answer.

I'll try that tomorrow.

My case is the following : Water is flowing onto a slope by its own weigth. The bc are cyclic in the slope direction. The top boundary is set to athmosphere so the air is free to enter and leave the system.

AT time 0 it gives :

https://picasaweb.google.com/lh/phot...eat=directlinkhttps://lh4.googleusercontent.com/-j...38/alphat0.png

In reality, only a thin layer of air is moved by the flowing water, but it is quickly slow down upwards by the athmosphere.

As I set up a athmosphere bc, I guess it is not too important to loose energy from the air field. Only the water field matter for me (the pressure and velocity of water would contribute almost totally to the water depth level). I think Shallow water equation are derived from Navier Stokes taking a free velocity bc for the water surface. So in my case, air should influate a minimum the water surface. (maybe I should decrease a lot air viscosity for that?)

However, since you are removing energy from a system, where the Navier-Stokes equations are solved simultaneously for both air and water, then you indirectly removes energy from the water, since energy from the water will be used to re-accelerate the air. Work put into the air, which is then removed in subsequent time steps.

The approach could be working for your case; you merely have to exhibit extreme caution.

/ Niels

I see what you mean, I will be carefull though.

I 'll post later my rersults,

cheers
joris

compile problem
 

hello joris.hey

I got same problem.

did you solve it ???

please advise me

To be sure
 

I am facing a similar problem and I have tried all the solutions mentioned :

alpha1*div(rhophi,U) => well then the air flow is not transported and saturates where it is generated.
U air to 0 or a percentage of the solution value (that should work but it does not to violent a create instability)
I tried increasing the density : it is a bit better
I tried also to do higher rho and Brackbill.

What was your final solution Bernhard?
are you still doing the alpha1*div and the increasing rho plus blackbill?

Thanks
Remi

Hi Remi,

I have been bothered by this problem for quite a long time. Which of those solutions mentioned in the earlier posts works best for you?
Cheers,

Yuanchuan Liu

Hi Yuanchuan,

I will increase the density of the air. But that depends on what you are trying to do?

My problem was when putting an object near the interface and moving it. Then I have vorticity (shedding) generated by the object (chose your frame it happens if you put a flow or move the object). And when the vorticity reaches the interface it kind of generates fast flow in the easy to transport air. Making it thicker (just rho_water / rho_air =100 is still big) really reduce this not necessary transport. Killing the flow in air is definitely not a good idea.

Now that I think about it, would making a viscosity gradient help? (high viscosity as you move up in the air to dissipate the energy or just higher air viscosity to damp the transport, it is kind of what people do just by making the grid more loose in the air? numerical viscosity is bigger!).

Let me know what is your final call.

Best

Remi

Hi Remi,

If I understand you correctly, you offer two ways to mitigate the effects of air convection:
1. Increase air density so that it is not so easy for air to transport.
2. Increase air viscosity so that even if high velocity occurs it will damp out quickly.
These two methods might work for my current case where a cylinder with a horinzontal thick plate oscillates vertically in xoz plane since air is not important here. Besides, I do not need to change the code, which is fine.
However, if I later want to take into account the aerodynamic force acting on the object, altering the properties of the air phase might be improper as air becomes important in this case. Of course, all the methods mentioned in this thread will change the real physics and are actually based on the assumption that air does not play a major part in my cases. I am wondering if there is any way to solve this problem without the need to alter the original problem. Really curious about how commercial solvers deal with this problem.
Anyway, thanks very much for sharing your thoughts with me. I will try them to see if they will make a difference.

Cheers,

Yuanchuan

a RANS in air only
 

Hi all,

I am still playing with this problem. I have recently tried a new idea. I activate turbulence in Air only. I created a new k-E RASModel where I use the scalar alpha1 to turn it on or off whether I am in air or water. This is quite convenient but might be wrong since I am in 2D. It however gives rather good results somehow in my case (at least qualitatively compares well with my experiments) I have a manually oscillating pincher near a free surface.

Any other comments on your side?


Best

Rmi

Hi,

I would like to make my modest contribution to this interesting subject. I agree that it is not the perfect solution but, nevertheless, it is better than getting kinetic turbulent energy divergence.

On the other hand, I have tested this loop and
1.- It is very slow
2.- It does not work with decomposed cases, since it only loops over the cells in the local processor

I suggest this solution for the velocity correction:

It is just an approximation to the Heaviside step function for the alpha1 scalar. The parameters can be changed in order to make it steeper and to modify the threshold value of alpha1

It is faster than the forAll loop, and it works fine in decomposed case.

Now I am thinking on how to modify this function so that it does not put to zero all the values of velocity where alpha1<0.1 but only the "peak" values, lowering them to an average velocity value. Any suggestion will be welcome.

Robert

Greetings to all!

I'm finding this issue with artificially high velocities very strange. I've taken a quick look into the whole thread and I couldn't find a single test case for this issue.

Can anyone provide a simple test case that demonstrates this?
Because I suspect that the problem for this has another origin, namely initialization problems, such as described here: http://www.cfd-online.com/Forums/ope...tml#post404292 - post #7

Best regards,
Bruno

Hi Remi!

How did you manage to activate the turbulence in air only? (by the way... shouldn't you activate it in water only?) I'm having troubles with using lookup for alpha1 in the turbulence model. When I use

the turbulence model compiles fine, but when running interFoam I get the error:

If I use alpha.water instead of alpha1 it doesn't compile... Did you had similar problems?

Thanks, Matej

Quick answer:
My guess is that you did this:
which is incorrect..

Try using this:

Dear wyldckat,

those interfacial-velocities occur e.g. in mutliphaseEulerFoam (2.3.x), too.

As a test case one can use the bubbleColumn tutorial.
To better study the behavior I modified U.air, an alpha.air/water in a way,
that there is no inflow of air.

inlet: alpha.air = 0 , alpha.water = 1
inlet: U.air = 0


As previously stated in the thread, when I reduce the density gradient of both phases,
the interface velocities disappear (e.g. 995 to 1000).

I already looked into the effects of e.g. surface tension&interface compression,
but the effect is only very little influenced by turning on VOF-interface compression with surface tension.

I also followed your advise and initialized the pressure field according to gravity/mass effects, but that did not make the velocities disappear, either.

//EDIT: My first funkySetField approach was based on the setFieldsDict setup. I now modified the blockMeshdict to create 100 cells in y-direction,
so it should be easier to initialize the correct pressure field. I'll update the post if I make any further progress

//EDIT2: With this initialization the pressure is dropping at the liquid-gas interface after one iteration, while liquid downflow and gas upward flow is indicated in the U-fields.
Case Files are attached.
It seems that the solver calculates the hydrostatic pressure, as if the first "liquid cell" is not fully acting with the liquid density onto the pressure, as the pressure drop,
is in the order of half a cell-height filled with water.

Comparing those results with a setup without pressure initialization, there is no significant improvement regarding the artificial velocities at the interface, thus here the pressure initialization might not help resolving the problem.




Another way to reduce the velocities is to use the implemented damping function of multiphaseEulerFoam,
as stated by Kent:
http://www.cfd-online.com/Forums/ope...tml#post396320

... but isn't this effecting the overall physics?

Might a damper only at the surface be a less powerful intervention?



//EDIT3:
http://openfoam.org/mantisbt/view.php?id=1379

Following this bug report, the currents at the interface seem to be a general issue in "stagnand" two phase flow?
Is there a general recommendation in cases were the interface currents are higher than the liquid velocities, e.g. induced by liquid natural convection?


Sincerely,
Dominik

I wonder is there any progress in addressing this issue?

portal: http://www.cfd-online.com/Forums/ope...tml#post594956

Here is a related topic with some info.
http://www.cfd-online.com/Forums/ope...tml#post576093

Hi!

I tried to apply this code but it doesn't seem to re-assign the values of U to 0 (vector(0,0,0,)) after I added this next to the convective term. Was this code correct in the first place?