[ngj]

Hi Joris

You could type

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

[joris.hey]

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?)

[ngj]

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

[joris.hey]

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

I 'll post later my rersults,

cheers
joris

[wes1204]

hello joris.hey

I got same problem.

did you solve it ???

please advise me

[rcarmi]

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)

- ghf*fvc::snGrad(rho)

By

- 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

[SailorLiu]

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

[rcarmi]

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

[SailorLiu]

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

[rcarmi]

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

[rcastilla]

Originally Posted by ngj

Hi Joris

You could type

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:

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

[wyldckat]

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

[matejmuller]

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

const volScalarField& alpha1 = mesh_.lookupObject<volScalarField>("alpha1");

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

--> 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<Type>(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

[wyldckat]

Quick answer:

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:

const volScalarField& alpha.water = mesh_.lookupObject<volScalarField>("alpha.water");

which is incorrect..

Try using this:

const volScalarField& alpha1 = mesh_.lookupObject<volScalarField>("alpha.water");

[dschmidt]

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.

  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

[Ya_Squall2010]

I wonder is there any progress in addressing this issue?

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

[kwardle]

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

[cuikahlil]

Originally Posted by pablodecastillo

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

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?

[2538sukham]

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)

- ghf*fvc::snGrad(rho)

By

- 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.

Wow. This actually works like a charm. Atleast in damBreak case, the air velocity is reduced.