[Lydia]

Hello Foamers,

being fairly new to the twoPhaseEulerFoam solver I'm trying to simulate a simple box - closed on all sides and the bottom, open to atmosphere on top.
(Boundary conditions see below).

The box is completely filled with water. There is no inflow, so I would actually expect, that nothing happens.
BUT unfortunately velocities (up to 89 m/s after 4.8s ) and upward flow develops within the box, which gets bigger with time. (See attached pictures.)

Form the pictures of the first written result after 0.1s you can see that the problem seems to evolve from the bottom.

Has anyone an explanation for that?

Looking forward to your answers.
Kind regards,
Lydia

My boundary conditions:

            walls                                 top    
Theta    zeroGradient;                     inletOutlet; uniform 1.0e-8;
Ua        fixedValue; uniform (0 0 0);   zeroGradient;
Ub        fixedValue; uniform (0 0 0);   zeroGradient;
alpha    zeroGradient;                      inletOutlet; uniform 1.0;
epsilon  zeroGradient;                      inletOutlet; uniform 10.0;
k          zeroGradient;                      inletOutlet; uniform 1.0;
p          buoyant pressure; uniform 0; fixedValue; uniform 0;

As initialization: 

the whole box is filled with alpha=1
via the setFieldsDict:

regions
(boxToCell   {   box (-999 -999 -999) (999 999 999);
 fieldValues (volScalarFieldValue alpha 1);}

[alberto]

Could you attach your case?

[Lydia]

For a better understanding , I attached the 0/ system/ and constant/ folder.

Hoping that this gives a better insight into what I did, I'm looking forward to any hint about possible reasons for the aforementioned behavior.

Best regards,

Lydia

[GerhardHolzinger]

Originally Posted by Lydia

The box is completely filled with water.

As initialization: 

the whole box is filled with alpha=1
via the setFieldsDict:

Please correct me if I am wrong, but alpha = 1 means the box is filled with air. Alpha is the volume fraction of the dispersed phase, which is defined if I am not wrong, in the transportProperties with the keyword dragPhase, but I am not 100% sure at this.

[Lydia]

Thanks Gerhard. You're correct with this:

alpha = 1 means the box is filled with air.

But the results still don't make sense to me...? No matter if it's water or air in the box - I still would expect nothing to happen. Why do i get these high velocities?

I also tried the same with alpha=0. Here the velocities are much smaller.
Does the solver maybe have problems with only air in the domain?

[Lydia]

Hello again,

I just performed a few tests, where I modified the bubble column tutorial and I get similar results to my previous tests.

In the first test I set alpha to 0 in the whole domain.
In the second test alpha was set to 1 in the whole domain.
All other settings were taken over from the tutorial setup.

From the attached pictures it is visible, that velocities evolve although there is no inflow or other sources that could causes flow within the domain. In the second test it even seems that air is flowing in from the bottom. Is there maybe a mistake in my modifications? (I attached the case files to this post).


Or has anyone an explanation for this?

Wishing a you all a nice weekend,

Lydia

[alberto]

Hi Lydia,

I checked your alpha0 case. I have the following comments:

1. There is a problem in how you define the bottom condition. Since it is a wall, the boundary condition in blockMeshDict should be changed to "wall" instead than of "patch".
2. The twoPhaseEulerFoam solver uses the "phase-intensive" form of the momentum equation. This form is more robust when the phase fraction approaches zero, however it will return non-zero velocities in such a limit (you can show that for alpha = 0, in absence of acceleration terms, it will return the terminal velocity). This is not a problem being alpha = 0 in those regions, but you need to know this to interpret your results.

I modified the boundary condition as I wrote above, and results look fine (see attachments for pictures and cases). You will see a non-zero velocity for phase a in the cells before the outlet for the reason at point 2. Pictures are taken at t = 20s.

[Lydia]

Hello Alberto,

Thank you for your comments and the attached pictures +case!

I didn't know that the definition of the bottom as a 'wall' instead of a 'patch' would make such a difference, since no special wall-functions were set and the velocities were set to 'fixedValue (0 0 0)' at the 'bottom' patch. But from your pictures I really seems to make a difference, so I'll keep that in mind.

Concerning your second comment I don't understand what you mean with the "terminal velocity"? I only started reading through H. Rusche's Thesis and hope to understand more soon.

Best regards,
Lydia

[alberto]

Hi Lydia,

I'll try to explain the terminal velocity in simple words with one example. Think to a single particle falling into a fluid. We only consider two forces: drag and external force due to the acceleration (gravity, for example).

If the particle initially has zero velocity, the acceleration will cause an increase in its velocity. However, the drag force exerted by the fluid onto the particle will increase (it depends on the relative velocity), and at some point it will balance the force caused by the acceleration. At that point the velocity of the particle won't change anymore, and such a velocity is the terminal velocity.

Best,

[Lydia]

Hello Alberto and all other interested Foamers,

being back at the office today, I just tested your modified case and found out that the modifications actually didn't make the difference.
The difference between your and my results come from a different visualization:
I assume, you visualized the cell-to-point filtered data of Ua and Ub
whereas I visualized the cell data.

The attached pictures show the different visualizations of the same result. I don't understand how Paraview averages these values, since it seems that after the averaging the velocity in the whole box (except at the very top) is 0, whereas without the averaging all values range between 0.3 and 0.31 (according to the 'Data Range' in Paraview). Is there an explanation for that?


Additionally, I tested another case which also produces incomprehensible (for me at least...) results. I took the same box and filled it partially with water (see attached case; filling via setFields). Somehow the water level increases with time although there is no inflow defined. Did I make a mistake in my setup? Or why does the water level rise?

I'm looking forward to your answers.

Best regards,
Lydia

[alberto]

Originally Posted by Lydia

Hello Alberto and all other interested Foamers,

being back at the office today, I just tested your modified case and found out that the modifications actually didn't make the difference.
The difference between your and my results come from a different visualization:
I assume, you visualized the cell-to-point filtered data of Ua and Ub
whereas I visualized the cell data.

The attached pictures show the different visualizations of the same result. I don't understand how Paraview averages these values, since it seems that after the averaging the velocity in the whole box (except at the very top) is 0, whereas without the averaging all values range between 0.3 and 0.31 (according to the 'Data Range' in Paraview). Is there an explanation for that?

As I said before, you can disregard any non-zero value of the velocity where alpha is zero. This is due to the formulation of the momentum equation, which you can find described in the paper of Oliveira, P.J., Issa, R.I., 2003. Numerical aspects of an algorithm for the Eulerian simulation of two-phase flows. International Journal for Numerical Methods in Fluids 43, 1177–1198.

You can post-process your ouput data by checking where alpha < alphaSmall (for example 1.0e-6), and set Ua to zero there, if you prefer to have a cleaner view.

Additionally, I tested another case which also produces incomprehensible (for me at least...) results. I took the same box and filled it partially with water (see attached case; filling via setFields). Somehow the water level increases with time although there is no inflow defined. Did I make a mistake in my setup? Or why does the water level rise?

Please, try setting dragPhase to a, since the drag phase is air in your problem.

Best,

[Lydia]

Hello Alberto,

Thanks for the reference to Oliveira's paper!

Originally Posted by alberto

As I said before, you can disregard any non-zero value of the velocity where alpha is zero. This is due to the formulation of the momentum equation, which you can find described in the paper of Oliveira, P.J., Issa, R.I., 2003. Numerical aspects of an algorithm for the Eulerian simulation of two-phase flows. International Journal for Numerical Methods in Fluids 43, 1177–1198.

It is now understandable for me that the solver has "problems" when alpha approaches 0.

Can I conclude from that, that it is not recommendable to model flows where the phases are completely separated in most regions and only penetrate each other in a small region of the domain?

In specific, I'm thinking of water (free-surface) flowing over an edge, plunging into a pool. The water air mixture (bubble development) in the pool is of special interest.


Best regards,

Lydia

P.S. Setting the drag phase to a in the earlier described case with water and air didn't change the results. The domain is still filling without reason.

[alberto]

Originally Posted by Lydia

Hello Alberto,

Thanks for the reference to Oliveira's paper!

It is now understandable for me that the solver has "problems" when alpha approaches 0.

Can I conclude from that, that it is not recommendable to model flows where the phases are completely separated in most regions and only penetrate each other in a small region of the domain?

The conclusion is not correct. The implementation of the code is exactly done to address the problem when alpha approaches to zero. The fact you see non-zero velocities where alpha -> 0 does not represent a problem at all, since the momentum is anyways zero there. If you do a literature search, a similar approach is used in many codes (CFDLib, NEPTUNE, ...).

You might want to see http://openfoamwiki.net/index.php/BubbleFoam and notice that the momentum equation is treated similarly to what suggested by Oliveira.

In specific, I'm thinking of water (free-surface) flowing over an edge, plunging into a pool. The water air mixture (bubble development) in the pool is of special interest.

That seems to me a different problem, where you might need surface-tracking capabilities,in addition to surface tracking. OpenFOAM has multiphaseEulerFoam to do this.

P.S. Setting the drag phase to a in the earlier described case with water and air didn't change the results. The domain is still filling without reason.

I attach a case which does not present this problem. Note that wall conditions for p were changed to zeroGradient. Please, let me know if you have questions.

Best,

[kwardle]

Lydia,
The momentum equation formulation in multiphaseEulerFoam is a little different from the one in bubbleFoam in terms of phase conservation, so it might behave a little different (hopefully better!) than what you have seen. Indeed if you are interested in a problem with both a free surface and dispersed phases also this might be useful to you. That said, the development of this solver was targeted at a three-phase flow (liquid-liquid-air) in which you have a free-surface for the liquid-air pairs and dispersion for the liquid-liquid part--that is, interface sharpening is either on or off everywhere for a given phase pair. I did tinker with an earlier two-phase version of the same solver in which I implemented 'dynamic' switching of interface sharpening (a la Cerne et al. 2001) which would essentially allow you to capture regions with a sharp interface and regions dispersed for the same phase pair. I may look at adding this into the multiphase version but so far I haven't needed it and have gone on to working on other things.
Hope this is helpful.
-Kent

[alberto]

One note on this, since I receive often emails on the topic

The fact that in twoPhaseEulerFoam the velocity of the dispersed phase (Ua) is not zero when the phase fraction alpha is zero is expected, since it is a property of the phase intensive form of the momentum equation, which is defined when alpha -> 0, and provides a finite non-zero value of Ua in such a case. This also does not represent a problem in general, since the sub-models are typically function of alpha*(1-alpha). (See Oliveira's paper for details).

As Kent said, multiphaseEulerFoam (and compressibleTwoPhaseEulerFoam) uses the conservative form of the momentum equation (and MULES to solve for the phase fractions and ensure they are bounded). To avoid the singularity when the phase fraction is small, the drag term is manipulated so that the values of alpha and beta have a strictly positive inferior bound (not zero). Similarly, the magnitude of the relative velocity is kept strictly positive. This treatment has the important advantage of keeping the conservative form, which is key in certain situations. You will still see however the non-zero velocity for a phase whose phase fraction is zero, and, as before this is fine (remove it in the post-processing).