InterFoam letting one phase leave the domain

Shibi · March 20, 2021, 17:30

Hello to all,

I would like to simulate the flow of a liquid and gas with InterFoam. I would like to have one of the phases (the gas) to leave the domain while the other (the liquid) fills the domain.

Does OpenFOAM have any boundaries for this?

What do you recommend to achieve this? Would it be achievable with a custom boundary condition that dynamically changes from outlet to wall based on the value of alpha?

Best regards

Shibi · March 28, 2021, 12:12

So, apparently this question has been going around in the forum for several years without a proper answer.

Some questions I still have:

how sudden should the outlet switch from a zeroGradient to a fixedValue ( 0 0 0)?

Is it necessary to have a similar dynamic boundary condition for the pressure (to switch from fixedValue to a zeroGradient?)

I am using the codedMixed boundary condition to write this. So far I have:

Code:

    myOutlet
    {
        type            codedMixed;

        refValue        uniform (0 0 0);
        refGradient     uniform (0 0 0);
        valueFraction   uniform 0;

        name    dynamicOutlet;   

        code
        #{
            // Code goes here

            // Gets the alpha field

           scalarField alphap =
                    patch().lookupPatchField<volScalarField, scalar>("alpha.liquid");

            alphap = max(alphap, scalar(0));
            alphap = min(alphap, scalar(1));

          // write value as a function of alpha
            this->valueFraction() = alphap;

        #};
    }

The code above is switching from zeroGradient to fixedValue based the value of alpha.
if alpha = 1, we get a pure fixedValue
if alpha =0 , we get a pure zeroGradient

Any contribution would be very welcome!

mAlletto · March 28, 2021, 12:35

For a rising air bubble in water I applied zeruGradient for alpha and this https://www.openfoam.com/documentati...alarField.html

For prgh at the upper boundary. In this way the air bubble left the upper boundary

Shibi · March 29, 2021, 08:05

Hi,

After some reading, I need to implement something like this:

$\begin{matrix} \vec{\mathbf{v}} = 0 \; \; \; \; if \; \alpha \geq 0.5 \\ \mathbf{\sigma}(\vec{\mathbf{v}}) \cdot \vec{\mathbf{n}} = 0 \; \; \; \; if \; \alpha < 0.5 \end{matrix}$

with:

$\mathbf{\sigma}(\vec{\mathbf{v}}) = -p \mathbf{I} + 2\mu \mathbf{D}$

$\mathbf{D} = \frac{1}{2}\left ( \nabla \vec{\mathbf{v}} + \nabla \vec{\mathbf{v}}^T\right )$

So, the first one is a no-slip boundary which is straightforward to implement.

Just need to loop over the patch faces and check the value of alpha is higher than 0.5 (or any specified value for alpha). If so, I set

Code:

this->refValue()[faceI] = vector::zero; 

this->valueFraction()[faceI] = 1.0;

The second one I do not know how to implement.

Any help be appreciated!

trailer · March 29, 2021, 09:09

Hello,

Although I do not have any input to give you, I am also interested in this matter!

Best of luck!

mAlletto · March 29, 2021, 11:04

Quote:

Originally Posted by Shibi

Hi,

After some reading, I need to implement something like this:

$\begin{matrix} \vec{\mathbf{v}} = 0 \; \; \; \; if \; \alpha \geq 0.5 \\ \mathbf{\sigma}(\vec{\mathbf{v}}) \cdot \vec{\mathbf{n}} = 0 \; \; \; \; if \; \alpha < 0.5 \end{matrix}$

with:

$\mathbf{\sigma}(\vec{\mathbf{v}}) = -p \mathbf{I} + 2\mu \mathbf{D}$

$\mathbf{D} = \frac{1}{2}\left ( \nabla \vec{\mathbf{v}} + \nabla \vec{\mathbf{v}}^T\right )$

So, the first one is a no-slip boundary which is straightforward to implement.

Just need to loop over the patch faces and check the value of alpha is higher than 0.5 (or any specified value for alpha). If so, I set

Code:

this->refValue()[faceI] = vector::zero; 

this->valueFraction()[faceI] = 1.0;

The second one I do not know how to implement.

Any help be appreciated!

You have three equations and four unknowns. If you fix the pressure than you have three equations and three unknowns from where you can obtain a relation for the velocity gradients

Shibi · March 30, 2021, 13:12

Quote:

Originally Posted by mAlletto

You have three equations and four unknowns. If you fix the pressure than you have three equations and three unknowns from where you can obtain a relation for the velocity gradients

Wouldn't it be better to create two boundaries. One that switches from zero velocity gradient to zero FixedValue and another for pressure that switches from fixedValue to zeroGradient?

Shibi · April 1, 2021, 08:58

Can anyone confirm this approach?

I have implemented this dynamic boundary, but now I am facing several convergence issues because the initial state of the problem has imposed value of pressure everywhere expect at the inlet ...

mAlletto · April 5, 2021, 12:51

Can you sketch your case? From the description above it's hard to guess what the exact problem is

Shibi · April 6, 2021, 06:17

Quote:

Originally Posted by mAlletto

Can you sketch your case? From the description above it's hard to guess what the exact problem is

Sure,

Filling of a plate with water (other kind of fluid) with air inside the cavity.

The water enters through an inlet and all other patches in the geometry should behave like a dynamic wall, letting the air out and keeping the water in.

mAlletto · April 6, 2021, 06:48

Where is the outlet and where the walls. The direction of g is missing

Shibi · April 6, 2021, 07:12

Quote:

Originally Posted by mAlletto

Where is the outlet and where the walls. The direction of g is missing

There is no outlet (at least explicitly specified). The walls should behave like an outlet up to and alpha of 0.5 and as a wall for values higher than 0.5 .

Gravity is not being considered in this problem.

mAlletto · April 6, 2021, 08:57

Quote:

Originally Posted by Shibi

Wouldn't it be better to create two boundaries. One that switches from zero velocity gradient to zero FixedValue and another for pressure that switches from fixedValue to zeroGradient?

This seem a good solution. How did you implement this. Seem you wrote a new boundary conditions

brunomramoa · April 6, 2021, 13:10

Hi,

I would also be very interested in this kind of boundary!

Based on your previous input, I wrote for U:

Code:

  wall
    {
        type            codedMixed;

        refValue        uniform (0 0 0);
        refGradient     uniform (0 0 0);
        valueFraction   uniform 0;

        name    dynamicOutlet;   

        code
        #{

       scalarField alphap =
         patch().lookupPatchField<volScalarField, scalar>("alpha.liquid");

        alphap = max(alphap, scalar(0));
        alphap = min(alphap, scalar(1));



        const scalar alphaValue = 0.5;

        // Loops over the patch
       forAll(patch(), faceI)
        {   

              if(alphap[faceI] > alphaValue)
              {
                // Imposes fixedValue
                  this->refValue()[faceI] = vector::zero; 
                  this->valueFraction()[faceI] = scalar(1); 
              }
              else
              {
                // Imposes zeroGradient
                  this->refGrad()[faceI] = vector::zero;
                  this->valueFraction()[faceI] = scalar(0);
              }
          
        }
      #};
    }

and for P

Code:


wall
    {
        type            codedMixed;
        refValue        uniform 0;
        refGradient     uniform 0;
        valueFraction   uniform 1;

        name    dynamicOutlet;   

        code
        #{
               scalarField alphap =
                     patch().lookupPatchField<volScalarField, scalar>("alpha.liquid");
                    
                alphap = max(alphap, scalar(0));
                alphap = min(alphap, scalar(1));

                 const scalar  alphaValue = 0.5;


                // Loops over the patch
                forAll(this->patch(), faceI)
                {

                        if (alphap[faceI] > alphaValue)
                        {
                            // Imposes zeroGradient
                            this->refGrad()[faceI] = scalar(0);
                            this->valueFraction()[faceI] = scalar(0);
                        }
                        else
                        {
                            // Imposes fixedValue
                            this->refValue()[faceI] = scalar(0); 
                            this->valueFraction()[faceI] = scalar(1);
                        }
                }
        #};
    }

I did some tests in a plate with 3 patches: inlet, outlet, walls with "standard boundaries ".

For:
U: Fixed velocity on the inlet, no-slip for the walls and zeroGradient for the outlet.
P: ZeroGradient for the inlet and the walls and fixedValue for the outlet

The simulation worked fine.

However, with the codedMixed boundary the convergence is not good...

Hopefully the code will get you started and someone else can fill in with a little CFD magic to make the boundary more well behaved

Shibi · April 21, 2021, 05:34

Any suggestions on a possible strategy to improve the convergence of this boundary condition?

Best Regards!

JuRe09 · August 17, 2022, 09:05

Quote:

Originally Posted by Shibi

Any suggestions on a possible strategy to improve the convergence of this boundary condition?

Best Regards!

Hey Shibi,
you made any progress imroving convergence for the coded BC?

I found this custom BC in the web; within the paper the code of the BC is also linked!: https://www.emerald.com/insight/cont...0190/full/html

However, in my case this BC blows up my simulation within the first iterations...

March 20, 2021, 17:30	InterFoam letting one phase leave the domain	#1
Shibi Member Join Date: Feb 2020 Posts: 90 Rep Power: 6	Hello to all, I would like to simulate the flow of a liquid and gas with InterFoam. I would like to have one of the phases (the gas) to leave the domain while the other (the liquid) fills the domain. Does OpenFOAM have any boundaries for this? What do you recommend to achieve this? Would it be achievable with a custom boundary condition that dynamically changes from outlet to wall based on the value of alpha? Best regards trailer likes this.

March 28, 2021, 12:12		#2
Shibi Member Join Date: Feb 2020 Posts: 90 Rep Power: 6	So, apparently this question has been going around in the forum for several years without a proper answer. Some questions I still have: how sudden should the outlet switch from a zeroGradient to a fixedValue ( 0 0 0)? Is it necessary to have a similar dynamic boundary condition for the pressure (to switch from fixedValue to a zeroGradient?) I am using the codedMixed boundary condition to write this. So far I have: Code: myOutlet { type codedMixed; refValue uniform (0 0 0); refGradient uniform (0 0 0); valueFraction uniform 0; name dynamicOutlet; code #{ // Code goes here // Gets the alpha field scalarField alphap = patch().lookupPatchField<volScalarField, scalar>("alpha.liquid"); alphap = max(alphap, scalar(0)); alphap = min(alphap, scalar(1)); // write value as a function of alpha this->valueFraction() = alphap; #}; } The code above is switching from zeroGradient to fixedValue based the value of alpha. if alpha = 1, we get a pure fixedValue if alpha =0 , we get a pure zeroGradient Any contribution would be very welcome! JuRe09 likes this.

March 29, 2021, 08:05		#4
Shibi Member Join Date: Feb 2020 Posts: 90 Rep Power: 6	Hi, After some reading, I need to implement something like this: $\begin{matrix} \vec{\mathbf{v}} = 0 \; \; \; \; if \; \alpha \geq 0.5 \\ \mathbf{\sigma}(\vec{\mathbf{v}}) \cdot \vec{\mathbf{n}} = 0 \; \; \; \; if \; \alpha < 0.5 \end{matrix}$ with: $\mathbf{\sigma}(\vec{\mathbf{v}}) = -p \mathbf{I} + 2\mu \mathbf{D}$ $\mathbf{D} = \frac{1}{2}\left ( \nabla \vec{\mathbf{v}} + \nabla \vec{\mathbf{v}}^T\right )$ So, the first one is a no-slip boundary which is straightforward to implement. Just need to loop over the patch faces and check the value of alpha is higher than 0.5 (or any specified value for alpha). If so, I set Code: this->refValue()[faceI] = vector::zero; this->valueFraction()[faceI] = 1.0; The second one I do not know how to implement. Any help be appreciated!

April 21, 2021, 05:34		#15
Shibi Member Join Date: Feb 2020 Posts: 90 Rep Power: 6	Any suggestions on a possible strategy to improve the convergence of this boundary condition? Best Regards! JuRe09 likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
[ICEM] How to generate an unstructured mesh for a 2 fluid domain rotor in ICEM CFD	Hazem24	ANSYS Meshing & Geometry	0	December 21, 2020 07:52
UDF to constrain one of the phase to stay in the domain	Koromajor	Fluent UDF and Scheme Programming	23	April 15, 2020 20:21
interFoam high air phase velocities	indy07cz	OpenFOAM Running, Solving & CFD	1	November 8, 2017 06:00
Sudden increase the residual of Maxwell's equations	hsezsz	CFX	2	October 13, 2016 07:58
phase fraction part of domain	RBJ	OpenFOAM Post-Processing	0	July 6, 2010 10:33

March 28, 2021, 12:35		#3
mAlletto Senior Member Michael Alletto Join Date: Jun 2018 Location: Bremen Posts: 616 Rep Power: 16	For a rising air bubble in water I applied zeruGradient for alpha and this https://www.openfoam.com/documentati...alarField.html For prgh at the upper boundary. In this way the air bubble left the upper boundary

March 29, 2021, 09:09		#5
trailer Member Join Date: Mar 2021 Posts: 39 Rep Power: 5	Hello, Although I do not have any input to give you, I am also interested in this matter! Best of luck!

April 1, 2021, 08:58		#8
Shibi Member Join Date: Feb 2020 Posts: 90 Rep Power: 6	Can anyone confirm this approach? I have implemented this dynamic boundary, but now I am facing several convergence issues because the initial state of the problem has imposed value of pressure everywhere expect at the inlet ...

April 5, 2021, 12:51		#9
mAlletto Senior Member Michael Alletto Join Date: Jun 2018 Location: Bremen Posts: 616 Rep Power: 16	Can you sketch your case? From the description above it's hard to guess what the exact problem is

April 6, 2021, 06:48		#11
mAlletto Senior Member Michael Alletto Join Date: Jun 2018 Location: Bremen Posts: 616 Rep Power: 16	Where is the outlet and where the walls. The direction of g is missing