# Artificial high velocities at the interface using interFoam

 Register Blogs Members List Search Today's Posts Mark Forums Read

 February 15, 2012, 12:32 #21 Senior Member   Niels Gjoel Jacobsen Join Date: Mar 2009 Location: Deltares, Delft, The Netherlands Posts: 1,754 Rep Power: 29 Hi Joris You could type Code: forAll(U, celli) { if (alpha1[celli] < 0.01) U[celli] = vector::zero; } 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

 February 15, 2012, 14:10 #22 New Member   Join Date: May 2011 Posts: 15 Rep Power: 8 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 : 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?)

 February 16, 2012, 03:51 #23 Senior Member   Niels Gjoel Jacobsen Join Date: Mar 2009 Location: Deltares, Delft, The Netherlands Posts: 1,754 Rep Power: 29 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 vsammartano likes this.

 February 16, 2012, 04:14 #24 New Member   Join Date: May 2011 Posts: 15 Rep Power: 8 I see what you mean, I will be carefull though. I 'll post later my rersults, cheers joris

 May 5, 2014, 11:59 compile problem #25 New Member   Baek, Donghae Join Date: Jan 2013 Location: Seoul Posts: 24 Rep Power: 6 hello joris.hey I got same problem. did you solve it ??? please advise me

October 16, 2014, 17:18
To be sure
#26
New Member

Remi Carmi
Join Date: Jul 2014
Posts: 15
Rep Power: 5
Quote:
 Originally Posted by Bernhard 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) Code: - ghf*fvc::snGrad(rho) By Code: - ghf*fvc::snGrad(rho)*2.0*fvc::interpolate(alpha1) 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.

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

 November 7, 2014, 07:52 #27 Member   Yuanchuan Liu Join Date: Oct 2012 Posts: 30 Rep Power: 6 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

 November 7, 2014, 15:13 #28 New Member   Remi Carmi Join Date: Jul 2014 Posts: 15 Rep Power: 5 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

 November 7, 2014, 15:44 #29 Member   Yuanchuan Liu Join Date: Oct 2012 Posts: 30 Rep Power: 6 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

 October 2, 2015, 04:49 a RANS in air only #30 New Member   Remi Carmi Join Date: Jul 2014 Posts: 15 Rep Power: 5 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

October 23, 2015, 04:42
#31
Member

Robert Castilla
Join Date: Apr 2009
Location: Spain
Posts: 95
Rep Power: 10
Quote:
 Originally Posted by ngj Hi Joris You could type Code: forAll(U, celli) { if (alpha1[celli] < 0.01) U[celli] = vector::zero; } 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
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:

Code:
Info<< "Umax : " << max(mag(U)).value() << " U average: " << average(mag(U)).value() << endl;

U *= 1/(1+exp(-1000*(alpha1-scalar(0.1))));

Info<< "Update U" << endl;
Info<< "Umax : " << max(mag(U)).value() << " U average: " << average(mag(U)).value() << endl;
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

 October 24, 2015, 13:12 #32 Super Moderator   Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,020 Blog Entries: 39 Rep Power: 109 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 vsammartano likes this. __________________ OpenFOAM: FAQ | Getting started Forum: How to get help, to post code/output and forum guide What am I doing/planning: blog/wiki Read this before sending me PM

 October 29, 2015, 21:33 #33 New Member   Matej Muller Join Date: Oct 2011 Location: Slovenia Posts: 23 Rep Power: 7 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 Code: const volScalarField& alpha1 = mesh_.lookupObject("alpha1"); the turbulence model compiles fine, but when running interFoam I get the error: Code: --> FOAM FATAL ERROR: request for volScalarField alpha1_ from objectRegistry region0 failed available objects of type volScalarField are 15 ( alpha.water_0 interfaceProperties:K nut alpha.water rho k p_rgh nu gh nu1 p rho_0 nu2 alpha.air epsilon ) From function objectRegistry::lookupObject(const word&) const in file /opt/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 198. FOAM aborting If I use alpha.water instead of alpha1 it doesn't compile... Did you had similar problems? Thanks, Matej Last edited by matejmuller; October 29, 2015 at 22:37.

October 31, 2015, 10:56
#34
Super Moderator

Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,020
Blog Entries: 39
Rep Power: 109
Quote:
 Originally Posted by matejmuller If I use alpha.water instead of alpha1 it doesn't compile... Did you had similar problems?
My guess is that you did this:
Code:
const volScalarField& alpha.water = mesh_.lookupObject<volScalarField>("alpha.water");
which is incorrect..

Try using this:
Code:
const volScalarField& alpha1 = mesh_.lookupObject<volScalarField>("alpha.water");

November 6, 2015, 07:30
#35
New Member

Dominik Schmidt
Join Date: Mar 2014
Posts: 11
Rep Power: 5
Quote:
 I've taken a quick look into the whole thread and I couldn't find a single test case for this issue

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.

Code:
  expressions
(
pressureWater
{
field p;
expression "1000*9.81*(0.7-pos().y)+(1-0.7)*1*9.81";
condition  "(pos().y<0.7) && (pos().y>=0)";
keepPatches 1;
}

pressureAir
{
field p;
expression "1*9.81*(1-pos().y)";
condition  "(pos().y<1) && (pos().y>=0.7)";
keepPatches 1;
}

);

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
Attached Files
 bubbleColumn.23x.artificialVelocities.zip (88.2 KB, 5 views)

Last edited by dschmidt; November 9, 2015 at 12:22.

 April 14, 2016, 03:24 #36 Member   YS Join Date: Jan 2010 Posts: 80 Rep Power: 9 I wonder is there any progress in addressing this issue? portal: http://www.cfd-online.com/Forums/ope...tml#post594956

 April 14, 2016, 12:01 #37 Senior Member   Kent Wardle Join Date: Mar 2009 Location: Illinois, USA Posts: 218 Rep Power: 14 Here is a related topic with some info. http://www.cfd-online.com/Forums/ope...tml#post576093 wyldckat likes this.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Saturn CFX 45 February 8, 2016 05:42 andersson.j OpenFOAM Running, Solving & CFD 0 February 8, 2011 11:43 Pankaj CFX 9 November 23, 2009 05:05 Andrea CFX 2 October 11, 2004 05:12 lego CFX 3 November 5, 2002 21:09

All times are GMT -4. The time now is 00:49.