[JunaidAhmad]

Hello,

I am trying to write a code on VOF Model in C++ using Finite Volume Method. Previously i used FVM and success fully write a code for Lid Driven Cavity and Natural Convection now i want to write it for 2D Dam Breaking like openFoam. Now i am trying it for some time but i am not getting the logic behind VOF.

Please give you suggestions, recommend books, or paper, and most importantly governing equations and algorithm.

My brain is not moving in any direction on that problem.

Please Help.

Junaid

[hilllike]

Dear Junaid,
I wrote a VOf code based on CICSAM (Ubbink, you could find his doctoral thesis online). This is an easy method based on original D-A VOF method, similar with FAC-VOF.

There are some other VOF methods like PLIC-VOF, which are different in surface Reconstruction. Maybe you can find something you need in a review (Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows) written by Denis Gueyffier.

Ren

[JunaidAhmad]

Hi, Thanks for your reply.

I have lots of questions for you.
First is,
as mentioned in the thesis that First solve "alpha" equation. and then momentum equation. As i have mentioned in my post that i have code the lid driven cavity. in that i did not use densities. kindly tell me how to handle density in momentum equation.
Problem i am working on is vertical column of water, so the source term is gravity, give me some tips to solve momentum eq and tips for the boundary conditions of volume fraction (alpha) equation.

Regards,
Junaid

[leflix]

Originally Posted by JunaidAhmad

As i have mentioned in my post that i have code the lid driven cavity. in that i did not use densities. kindly tell me how to handle density in momentum equation.

first of all you should change your solver in a variable density , and variable viscosity navier-stokes solver
Then you have to solve alpha equation, which is I guess the volume fraction.
One you got the distribution of alpha, ro=alpha*ro1 + (1-alpha)* ro2
where ro1 is the density of the first fluid, ro2 is the density of the second fluid
Do the same with the viscosity
Then solve your navier-stokes equation with new ro and new viscosity.
once you got your velocity distribution you can advect alpha for the next time step, and so on...

Problem i am working on is vertical column of water, so the source term is gravity, give me some tips to solve momentum eq and tips for the boundary conditions of volume fraction (alpha) equation.
Junaid

generally Neumann boundary condition for alpha should be ok.

[hilllike]

Originally Posted by JunaidAhmad

Hi, Thanks for your reply.

I have lots of questions for you.
First is,
as mentioned in the thesis that First solve "alpha" equation. and then momentum equation. As i have mentioned in my post that i have code the lid driven cavity. in that i did not use densities. kindly tell me how to handle density in momentum equation.
Problem i am working on is vertical column of water, so the source term is gravity, give me some tips to solve momentum eq and tips for the boundary conditions of volume fraction (alpha) equation.

Regards,
Junaid

Leflix gives a good answer of your questions.

In my code I use a smoothed matrix to smooth the VOF function to make sure the density calculated from VOF function as smoothed as possible according to a master thesis.That is a good thesis can give you a quick understand of numerical method for two phase flow named'' A navier stokes solver for single and two phase flow'' by Kim Motoyoshi Kalland and you could find a source code used by the author in that thesis.


hope that can help you.

Ren

[leflix]

Thanx Ren for this interesting reference ! It will definitely help

[JunaidAhmad]

Thanks guys,

how can we handle the boundary condition regarding Atmospheric pressure in a simple vertical column of water or dam break problem.

Domain is surrounded with three walls and an open surface.

at the wall can we use dp/dx or dp/dy equal to zero. by using p[0,j]=p[1,j] where j=0 to ny for zero wall gradiant. and similarly u and v both equal to zero at walls.

now the question is what would be the condition at the open surface for both pressure and velocity.

and most importantly what would be the boundary condition for the volumn fraction (alpha ) on walls and the open surface.

Kindly correct me if i am wrong. because when i run the simulation the volumn fraction or the phase contours shows that water phase disappears or sometimes it blows up when it disappears it moves downwards not sidewise.
please comment.

[leflix]

Originally Posted by JunaidAhmad

how can we handle the boundary condition regarding Atmospheric pressure in a simple vertical column of water or dam break problem.

Domain is surrounded with three walls and an open surface.

at the wall can we use dp/dx or dp/dy equal to zero. by using p[0,j]=p[1,j] where j=0 to ny for zero wall gradiant. and similarly u and v both equal to zero at walls.

now the question is what would be the condition at the open surface for both pressure and velocity.

it depends on you velocity-pressure algorithm. If you use SIMPLE like algorithm or projection method neumann boundary condition for the pressure is a good choice. Your implementaion saying p[0,j]=p[1,j] is right.
For velocity BC instead of no slip use rather free slip.

and most importantly what would be the boundary condition for the volumn fraction (alpha ) on walls and the open surface.

check my previous answer I already answered to you on that issue

Kindly correct me if i am wrong. because when i run the simulation the volumn fraction or the phase contours shows that water phase disappears or sometimes it blows up when it disappears it moves downwards not sidewise.
please comment.

To identify from where the problem could arise:
First verify your Navier-Stokes solver in a pure monophase case, but with a given distribution of density and viscosity.
Then check your solver for the transport equation of the mass fraction alone with a constant velocity field.
Then couple the both

[hilllike]

Originally Posted by JunaidAhmad

Thanks guys,
Kindly correct me if i am wrong. because when i run the simulation the volumn fraction or the phase contours shows that water phase disappears or sometimes it blows up when it disappears it moves downwards not sidewise.
please comment.

Leflix is really kind.

I don't know what kind of VOF you used.
When the water phase disappeared, you should check your volume convection equation.
for example in one dimensional problem, inside the water phase (not the surface) Ap*Fp(n+1)=b0*Fp(n)+Ae*Fe+Aw*Fw, if Ap>b0+Ae+Aw then the water phase will disappear quickly.

[JunaidAhmad]

Hi Guys,

As leflix suggested to use it for mono phase . So i use it on lid driven cavity.
Now the eq i used is
dF/dt+d(uF)/dx+d(vF)/dy=0

discritize Form

fn[col][row]=f[col][row]-(ue*(f[col+1][row]+f[col][row])-uw*(f[col-1][row]+f[col][row])+vnn*(f[col][row-1]+f[col][row])-vs*(f[col][row]+f[col][row+1]))/apo

where
apo=(dx*dy)/dt

i am using uniform grid 5 x 5

i started with initializing the fn field by
0 1 2 3 4 5
-------------
0| 0 0 0 0 0 0
1| 0|0 0 0 0| 0
2| 0|1 1 0 0| 0
3| 0|1 1 0 0| 0
4| 0|1 1 0 0| 0
5| 0 0 0 0 0 0

the outer nodes are the ghost nodes.

the above discretize form is solved using explict formulation only for the internal nodes.

Now the problem is that the total volume with in the cavity reduces and become negtive. for the above grid the total volume is 6. as i iteration increase the total volume decrease.

So any suggestion to solve this problem.
(The problem is that the mass with in the system is not conserved. is it due to some boundary condition or anything else. the code i am using for lid driven cavity is working fine and the velocity profiles are compared with ghia et al. for this problem i am using 81 x 81 grid).

[JunaidAhmad]

Hi Guys,

As leflix suggested to use it for mono phase . So i use it on lid driven cavity.
Now the eq i used is
dF/dt+d(uF)/dx+d(vF)/dy=0

discritize Form

fn[col][row]=f[col][row]-(ue*(f[col+1][row]+f[col][row])/2-uw*(f[col-1][row]+f[col][row])/2+vnn*(f[col][row-1]+f[col][row])/2-vs*(f[col][row]+f[col][row+1])/2)/apo

ue =u velocity at east face
uw =u velocity at west face
vnn =v velocity at north face
vs =v velocity at south face

where
apo=(dx*dy)/dt

i am using uniform grid

i started with initializing the "f" (volume fraction) field by let say we have grid 5 x 5
0 1 2 3 4 5
-------------
0| 0 0 0 0 0 0
1| 0|0 0 0 0| 0
2| 0|1 1 0 0| 0
3| 0|1 1 0 0| 0
4| 0|1 1 0 0| 0
5| 0 0 0 0 0 0

the outer nodes are the ghost nodes.

the above discretize form is solved using explict formulation only for the internal nodes i.e. from 1 to 4 in this case and calculate new f field and copy new field in the old field by copying the values.

Now the problem is that the total volume with in the cavity reduces and become negtive. for the above grid the total volume is 6. as iteration increases the total volume decreases.

So any suggestion to solve this problem.
(The problem is that the mass with in the system is not conserved. is it due to some boundary condition or anything else. the code i am using for lid driven cavity is working fine and the velocity profiles are compared with ghia et al. for this problem i am using 81 x 81 grid).

[JunaidAhmad]

Will some body give some input.

[hilllike]

Originally Posted by JunaidAhmad

Hi Guys,

As leflix suggested to use it for mono phase . So i use it on lid driven cavity.
Now the eq i used is
dF/dt+d(uF)/dx+d(vF)/dy=0

discritize Form

fn[col][row]=f[col][row]-(ue*(f[col+1][row]+f[col][row])/2-uw*(f[col-1][row]+f[col][row])/2+vnn*(f[col][row-1]+f[col][row])/2-vs*(f[col][row]+f[col][row+1])/2)/apo

ue =u velocity at east face
uw =u velocity at west face
vnn =v velocity at north face
vs =v velocity at south face

where
apo=(dx*dy)/dt


i am using uniform grid

i started with initializing the "f" (volume fraction) field by let say we have grid 5 x 5
0 1 2 3 4 5
-------------
0| 0 0 0 0 0 0
1| 0|0 0 0 0| 0
2| 0|1 1 0 0| 0
3| 0|1 1 0 0| 0
4| 0|1 1 0 0| 0
5| 0 0 0 0 0 0

the outer nodes are the ghost nodes.

the above discretize form is solved using explict formulation only for the internal nodes i.e. from 1 to 4 in this case and calculate new f field and copy new field in the old field by copying the values.

Now the problem is that the total volume with in the cavity reduces and become negtive. for the above grid the total volume is 6. as iteration increases the total volume decreases.

So any suggestion to solve this problem.
(The problem is that the mass with in the system is not conserved. is it due to some boundary condition or anything else. the code i am using for lid driven cavity is working fine and the velocity profiles are compared with ghia et al. for this problem i am using 81 x 81 grid).

I think you should read some papers about vof.
You can't discretize the volume convection in this form. VOF method is not easy like that.

[JunaidAhmad]

any recommendation for papers?