In case of the two methods HOL and SEM implemented in PHOENICS code for simulating free surface flows, the continuity equation was formulated in term of volumetric conservation as

D(Ln(RHO))/Dt+du/dx+dv/dy+dw/dz=0.

REF: Free surface flow and heat transfer in cavities: the sea algorithm, K.A Percilious & al 1995 ,Numerical heat transfer Part B, 27:487-507, 1995

1- what are the steps followed to discretise this equation by finite volume method ?

2- What is the discretized form of this equation ?

I don't found answers in Cham site.

Thank you.

Dear Medhi

The steps followed to discretize this (continuity) equation in FV is the same general principle followed by PHOENICS in any general transport equation, ie. take the volume integral over a control volume (scalar or pressure control volume in this case) and use Gauss theorem to substitute div(V(u,v,w)) for fluxes at cell's faces. The way in which PHOENICS determines the fluxes is by using first order upwind scheme.

Nevertheless, it is a very delicate equation which is coupled to pressure equation by SIMPLEST algorithm. It is used as source term into pressure correction equation, so it should be treated carefully if you want to "change" anything on it. I do not suggest you this.

Depending on the type of flow you have the final discrete equation for a general node p with neighbours e(east), w(west), n, s, h(high) and l(low); is like this

Ae*Ue+(Aw*Uw)+(An*Vn)+(As*Vs)+(Ah*Wh)+(Al*Wl)+(Ap* rho_p)_n+1+(Ap*rho_p)_n = 0

If you flow is incompressible you should elliminate the last two terms corresponding to compressibility effect; if not, the indices n+1 and n in those terms correspond to the actual and previous time steps respectivelly.

The way in which the equations are discretized in PHOENICS is partially summarized in TR99 manual. It is old (mostly refered to shareware version), but instructive. You can buy it in CHAM's web store, or ask to friends in our forum.

Regards

Kike

The single phase continuity equation can be written as:

d(rho)/dt + d( rho*ui )/dxi = 0

This can also be written as:

D(ln(rho))/Dt + d(ui)/dxi = 0

Assuming incompressible flow:

d(ui)/dxi = 0

This implies that knowledge of the density is immaterial for the solution of the continuity equation if the flow is incompressible.

1- I don't undrestand how the continuity condition in terms of volumetric conservation is valid even when the density changes from point to point at the interface.

2- Why this form (in terms of volumetric conservation) is more easy to use and to program after discretisation.

Dear Medhi

VOF method is an approximate method to track the interfaces' dynamics in two phase flows. You are using a single flow (variable properties) to model a system of two sets of motion equations and interface boundary conditions, so interface's dynamic is treated in an approximate way. There is no sharp line (or surface) between two well defined regions; instead, you have (at least) one cell with a "mixture" of the two phases (without to mention surface tension).

Originally (ie. rigurously) you should have a continuity eq. for each phase which should give you the final eq. you had writen (for incompressible phases). Therefore, you have that for both phases, the mass conservation eq. is the same! ie. div(v)=0

(1) That is exactly why you can use the volumetric form for continuity eq.. You only have to guarantee a divergence free velocity field in the whole domain, no matter the phase you are, even in interface cells.

(2) The answer for the second question is closely related. Now, let suppose you include the density term into the equations. Formally there is no inconsistency, because you are adding a zero for each continuity eq.

As there is no eq. for pressure (becasue incompressibility) and it is derived from continuity eq., to determine the pressure you should satisfy mass conservation in cells with interfaces considering the mass balance for both phases at a time.

Now, lets go for the discrete aspect of it. To do this balance you should know where (cell faces) they come from and where each phase go to (whit its density). To write correctly the continuity eq. for interface cells, it forces you to know: (a)interface orientation with respect to each cell face and (b)the partial area of each cell face occupied by each phase.

In 2D simulations it could be hard, but try to imagine all this in 3D. Moreover, put the usual density ratio you have in two phase problems (1/1000) and you will have a partial ansuwer to your second question.

Regards

Kike