CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   VOF method (https://www.cfd-online.com/Forums/openfoam-solving/58063-vof-method.html)

Vinay Ramohalli Gopala (Gopala) March 7, 2005 13:41

Hello, I am working on mod
 
Hello,

I am working on modeling of mass-transfer between two immiscible liquids. Recently I started using OpenFoam.

I have tried some simple simulations by manipulating the dambreak problem,(for example a denser liquid drop falling in a less denser liquid and also the rayleigh-taylor instability problem).

I intend to use VOF for my project and it would be very helpful to know more about the method implemented in the present version of OpenFoam, i.e. is it CICSAM or a different one ?

Also I would like to know if the option in interFoam - movingMesh is the same as Adaptive grid refinement around the interFace ?

Thanks in advance.

Henry Weller (Henry) March 7, 2005 13:49

> i.e. is it CICSAM or a diff
 
> i.e. is it CICSAM or a different one

It's a new technique I developed a few years ago to resolve some of the fundamental problems of CICSAM and other traditional VOF interface compression methods. The differences have already been debated at length, have a look through previous threads on the subject.

> Also I would like to know if the option in interFoam -
movingMesh is the same as Adaptive grid refinement around
the interFace ?

No.

Ali (Ali) March 7, 2005 21:24

I was wondering what's the ac
 
I was wondering what's the actual and preferred expression for the compressive term in gamma equation and a very brief explanation of how it is derived. The formula for 'phir' in the code differs from what Henrik has mentioned in (Eq. (4.15)). In the code:

phir=cGamma()*mag(phi/mesh.magSf())*interface.nHatf()

while Henrik's suggests (if I have formulated it correctly):

phir=cGamma*nhatf*max(nhat * phi / mesh.magSf()**2 )

Are they really different (am I missing something) and if yes, which one is better?


2) where is the smoothing function for 'gamma' in interFoam? Does it have a significant effect on VOF performance or just affects surface tension prediction?

thanks

Henry Weller (Henry) March 8, 2005 03:58

I have tried various options
 
I have tried various options for the compression term in the gamma equation and the one I recommend is the one that is currently in interFoam. It is not derived, it is selected and there are other choices. If you would like to find out which option is best for you try them out on your case.

There is no "smoothing" function for gamma in interFoam, it was found to be detrimental to the overall performance of the interface capturing.

gopala March 9, 2005 15:49

Hi Henry, Can I get any ref
 
Hi Henry,

Can I get any references for the VOF method you have implemented ?

Thanks

henry March 9, 2005 15:56

Sorry, I haven't written many
 
Sorry, I haven't written many papers; it's too time consuming and stops me doing the interesting work.

hjasak March 9, 2005 16:55

There is a decent description
 
There is a decent description of the implemented Interface-Capturing Methodology (probably a bit out of date now) in a PhD Thesis by dr. Henrik Rusche:

Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions, Henrik Rusche, Imperial College of Science, Technology & Medicine, December 2002.

Hrv

henry March 9, 2005 17:00

That describes an old version
 
That describes an old version of my method, not the one currently in interFoam.

sega December 15, 2008 09:44

I'm currently having a look in
 
I'm currently having a look into the gammaEqn.H from OF 1.4.1

Well, some parts look like the equation (4.15) from Rusches PhD thesis.

But still I'm not sure about some issues:

- Where is the time derivative of gamma?
- What is interface.nHatf()? Is it the normal vector to the cell face?
- What doese MULES::explicitSolve01(gamma, phi, phiGamma) mean?

Thanks so far!

caw December 15, 2008 10:30

Dear Sebastian, the time de
 
Dear Sebastian,

the time derivative is not there because of the usage of MULES. Look at OpenFOAM-1.5.x\src\finiteVolume\fvMatrices\solvers\MULES\MU LES.H for details:

" MULES: Multidimensional universal limiter with explicit solution.
Solve a convective-only transport equation using an explicit universal
multi-dimensional limiter.
Parameters are the variable to solve, the normal convective flux and the
actual explicit flux of the variable which is also used to return limited
flux used in the bounded-solution. "

Best regards
Christian

eberberovic December 15, 2008 12:19

The time derivative is account
 
The time derivative is accounted for within the MULES solver. In the gammaEqn.H only explicit fluxes are calculated, which are needed in MULES.

The interface.nHatf() represents a cell face unit interface normal flux. It is evaluated from the dot product of the cell face surface vector and the interface unit normal calculated at the cell face:

nHatf_ = nHatfv & Sf,

where

nHatfv = gradGammaf/(mag(gradGammaf) + deltaN_).

and deltaN is a stabilization factor for the case of gradGammaf = 0.

For the full implementation look in src/transportModels/interfaceProperties/interfaceProperties.C

Regards.

mer December 16, 2008 05:15

Hi, I have two questions rela
 
Hi,
I have two questions related to the interFoam solver:

1) I used the utility "barycenter" posted by Sebastian. It works fine with axisymmetric cases. In 2D or 3D cases and using the utility, it doesn't work (the position of the bubble is in decrease when time increase). Is there a signification for that?
I find also that the position of the bubble increase when the value of Min(gamma) is negative. Is there a relation with Min(gamma) and the barycenter of the bubble?!

MULES: Solving for gamma
Liquid phase volume fraction = 0.99921 Min(gamma) = -2.20718e-11 Max(gamma) = 1

2) I want to switch to OF-1.5. In the file interFoam.C there is an additional line comparing to OF-1.4.1:

p = pd + rho*gh;

Is it a construction of the pressure p field or what? if yes, how can I get values of this field for different times.

Best regards

mer December 16, 2008 05:47

Hi, I have two questions rela
 
Hi,
I have two questions related to the interFoam solver:

1) I used the utility "barycenter" posted by Sebastian. It works fine with axisymmetric cases. In 2D or 3D cases and using the utility, it doesn't work (the position of the bubble is in decrease when time increase). Is there a signification for that?
I find also that the position of the bubble increase when the value of Min(gamma) is negative. Is there a relation with Min(gamma) and the barycenter of the bubble?!

MULES: Solving for gamma
Liquid phase volume fraction = 0.99921 Min(gamma) = -2.20718e-11 Max(gamma) = 1

2) I want to switch to OF-1.5. In the file interFoam.C there is an additional line comparing to OF-1.4.1:

p = pd + rho*gh;

Is it a construction of the pressure p field or what? if yes, how can I get values of this field for different times.

Best regards

sega December 16, 2008 12:54

Thanks so far for your respons
 
Thanks so far for your responses.
I will try to get them step by step (and respond if there are any questions).

First of all:
Is there a reference for MULES in some kind of citable form?

sega December 16, 2008 14:57

Dear Edin. Thanks for your
 
Dear Edin.

Thanks for your response.
I have read your answer carefully and looked the code up.

But I still have some questions related to the code.
Is this the correct analytical representation of the code?
http://www.cfd-online.com/OpenFOAM_D...es/1/10363.png

I'm not sure what to do with the deltaN.
In the code it's calculated like this:
1e-8/average(gamma.mesh().V()), 1/3)

What does gamma.mesh().V() mean?

See you (maybe at SLA ...)

eberberovic December 17, 2008 06:19

Sebastian, I will come up t
 
Sebastian,

I will come up to you. The coefficient cGamma is not in the analytical transport equation, but it is used in modeling the relative velocity.

gamma.mesh().V(), for a non-moving mesh simply takes the volumes of control cells throughout the domain, so the average of it to the power of 1/3 gives a representative cell dimension (lenght).

Regards.

sega December 17, 2008 10:38

Dear Edin. Thanks for your
 
Dear Edin.

Thanks for your answere.
So, the equation will take this form?

http://www.cfd-online.com/OpenFOAM_D...es/1/10372.png

bobatpurdue December 18, 2008 23:28

On a side comment, I plan on u
 
On a side comment, I plan on using interFoam to do 2-D flows. Does anyone know how to implement proper BC?

markc December 24, 2008 04:59

Hello All, One new question
 
Hello All,

One new question related to VOF solvers (interFoam).
I modelled a ship sailing in a basin with water and air. The watersurface at the Inlet and internal has been set at Z=4 m. After a few seconds of physical time (so not CPU time) the solver crashes on very low deltaT. Looking the results in paraview one can see what happens:
http://www.cfd-online.com/OpenFOAM_D...es/1/10436.jpg

The picture shows the outlet of the domain. It looks like the watersurface is attracted towards Z=0 at the outlet patch. Maybe it might eventually go down further, provided that the solver should not crash.
However, after that I translated the geometry such that the initial waterlevel is at Z=0 and kept the rest unchanged. Than the solution runs smooth as expected.
BC's on the outlet:
gamma and U: zeroGradient
pd: fixedValue uniform 0

My questions:
- is there any mechanism that forces a surface towards Z=0?
- Has the problem described here something to do with continuity and/or mass conservation?
- I set BC's on the top of the domain at fixedValue (vesselspeed) for U and zeroGradient for pd, this to reduce strong air vortices which influence solution stability. Should I better set the BC's on the top at atmosphere?

Hope someone can shine his/hers light, not because I got stuck (I found a work around) but out of curiousity,

Brgds,

Mark

Mark

sega December 24, 2008 06:22

If this is an outlet than mayb
 
If this is an outlet than maybe the water is simply "leaking" out of it.
Maybe you can try fixedValue 0 for U at the out- and inlet?

markc December 24, 2008 07:54

Hello Sebastian, The inlet
 
Hello Sebastian,

The inlet has fixedValue U uniform (5 0 0) (approx). I can try to set this BC on outlet as well but I think this BC is too strong, it may influence the results too much. Isn't it very common to set fixedValue U on inlet and zeroGradient on outlet and for p the reverse?

Brgds and nice days,

albcem December 25, 2008 02:02

Has anyone any comment on the
 
Has anyone any comment on the "bug" side of things, as opposed to the work arounds?

I have been working on similar flow setups, with the added complication that the air/water interface does not line up with the principle coordinate system. So the air/water interface is not at z=constant.

So far I have not been able to achieve a stable simulation, which blows up, about after 1.5s of physically simulated time. Of course, I may have other issues with the case setup, but elimination of other possibilities is always a help.

So if the interfoam or rasinterfoam solver buggy, if the air/water interface is not at z=0???

Thanks.

albcem

markc December 29, 2008 06:35

Albcem, One thing which com
 
Albcem,

One thing which comes to my mind: is your gravity vector set correctly?

Brgds,

Mark

albcem December 29, 2008 15:46

Hello Mark, Thanks for the
 
Hello Mark,

Thanks for the tip. However I double checked the consistency of the gravity direction. It checks out correctly.

In my latest boundary condition settings, here is what I impose:

U:

atmospheric_wind_profile-boat_velocity at the air side,

-boat_velocity at the water side.

These inlet and outlet conditions are imposed at all vertical walls.

Zero-gradient conditions at the top and bottom faces of the rectangular domain should allow for the flow to preserve conservation laws. (I get 1e-14 continuity error in each time-step)

pd:

Impose the pressure only at the top of the atmosphere following the Wigley hull example.

k-epsilon:

Impose inlet & inletOutlet types with non-zero k & epsilon values calculated according to Fluent manual recommendations using reference lengths typical of tested boats. These are imposed at vertical walls.

At top and bottom patches, I have zero-gradient.

----

For all variables, I have also tested zeroGradient, inletOutlet types at several of the vertical and horizontal walls as sense allowed so. All have been fruitless.

Possibilities on where things may be going wrong:

The wind profile I impose, obeying a basic 1/7 power law, does not go to zero exactly at the air/water interface which may be causing high gradients there. However, I observe descent convergence on the simpleFOAM solution. So I have doubts on this argument.

I use tetrahedral mesh, with seemingly fine enough resolution at the air/water interface (40mm->320mm variation). Can this be an issue? I have read papers in which tetrahedral mesh with such resolution seems to do fine...

I may be introducing too much or incorrect turbulence into the flow messing things up. I will try a laminar solution to see if the blow up is delayed or eliminated - As the simpleFLOW solution converges I doubt this one as well...

Out of ideas, anything else I can check? Would it help to upload the case somewhere?

Albcem

markc December 30, 2008 08:00

Hello Albcem, I do not know
 
Hello Albcem,

I do not know your BC's. However the air side of boat-like case setups often show high air velocities, causing instabilities. In my cases I am not interested in the air side, I just set the same uniform fixedValue on U on the air side as on the water side. I only distinghuish these two by gamma.



Attached a 0 directory containing my BC's. Only the gamma file is to be filled in for the actual case setup. I wrote a small utility for this (see elsewhere on this forum). It sets gamma 1 on internal field and Inlet for Z<=0 and gamma=0 elsewhere.
I also experienced instable solutions using tet meshes. Since OF-1.5 I use snappyHexMesh, which I strongly suggest.
However, we started with a question about Z=0 reference in VOF calculations. What is actually your problem and case setup?

Brgds,

Mark

markc December 30, 2008 08:02

Hello Albcem, I do not know
 
Hello Albcem,

I do not know your BC's. However the air side of boat-like case setups often show high air velocities, causing instabilities. In my cases I am not interested in the air side, I just set the same uniform fixedValue on U on the air side as on the water side. I only distinghuish these two by gamma.

http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif pd
http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif U
http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif gamma

Attached a 0 directory containing my BC's. Only the gamma file is to be filled in for the actual case setup. I wrote a small utility for this (see elsewhere on this forum). It sets gamma 1 on internal field and Inlet for Z<=0 and gamma=0 elsewhere.
I also experienced instable solutions using tet meshes. Since OF-1.5 I use snappyHexMesh, which I strongly suggest.
However, we started with a question about Z=0 reference in VOF calculations. What is actually your problem and case setup?

Brgds,

Mark

sega January 4, 2009 08:56

Hello. Happy New Year to al
 
Hello.

Happy New Year to all of you!

I have some further questions regarding the VOF methode in interFoam.

Concerning the gamma-equation (gammaEqn.H):
phir is limited to max(phic)!
Is this the maximum phic in the cell or in the whole computational domain?
Rusche wrote about limiting the phir in the transition area. This looks more like limiting it tho the whole domain ...

Concerning the U-equation (UEqn.H):
Am I right that

fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muf, U)
- (fvc::grad(U) & fvc::grad(muf))
)

is an equation for face values (before interpolation to the cell centers)?

So

(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(g amma)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(pd)
) * mesh.magSf()

has to be of this kind as well to satisfy the momentum predictor?

Since the snGrad-Schemes are face-interpolation-Schemes I think fvc::interpolate is one as well?
So what is fvc::reconstruct used for?
The Programmers and User Guide keep quiet about these two fvc-Schemes.

Greetings. S.

sega January 6, 2009 04:56

Let me put my question in anot
 
Let me put my question in another way:
How is the interface.sigmaK() discretisized?

Is this a source term? Then I suppose it is calculated explicitly, since it is called within fvc::interpolate?
But then it is a volume field?
Shouldn't it be a surface field?

amghsu January 12, 2009 05:31

Hello Sebastian, I am also
 
Hello Sebastian,

I am also studying the code of interfoam.

As I know of the interfoam solver. The sigmaK() is the product of the curvature and the surface tension coefficient. Since the surface tension effect is treated by CSF model, interfoam discretisized it as body force term.

You can read the PhD Thesis by dr. Henrik Rusche P.154~P.155 for further infomation.

Jim

sega January 12, 2009 11:45

Thank you Jim. Then, how is
 
Thank you Jim.

Then, how is this body force discretisized?
As far as I got it, sigmaK() is explicitly calculated (fvc) at the cell faces.
Does it need further treatment?

Greetings. S.

michi January 13, 2009 08:07

Hallo Sebastian, As far as
 
Hallo Sebastian,

As far as I see interface.sigmaK() (product of surface tension and curvature) returns a volScalarField. These values are then interpolated to the cell faces with fvc::interpolate(). Different schemes are here possible for example linear (= central differencing).

Michael

sega January 13, 2009 15:11

Thank you Michael. Ok, sigm
 
Thank you Michael.

Ok, sigmaK() returns a volScalarField following interfaceProperties.H:

tmp<volscalarfield> sigmaK() const
{
return sigma_*K_;

}

And K_ is called with calculateK() in interfaceProperties.C:


void interfaceProperties::calculateK()
{
const fvMesh& mesh = gamma_.mesh();
const surfaceVectorField& Sf = mesh.Sf();

// Cell gradient of gamma
volVectorField gradGamma = fvc::grad(gamma_);

// Interpolated face-gradient of gamma
surfaceVectorField gradGammaf = fvc::interpolate(gradGamma);
//gradGammaf -=
// (mesh.Sf()/mesh.magSf())
// *(fvc::snGrad(gamma_) - (mesh.Sf() & gradGammaf)/mesh.magSf());

// Face unit interface normal
surfaceVectorField nHatfv = gradGammaf/(mag(gradGammaf) + deltaN_);
correctContactAngle(nHatfv.boundaryField());

// Face unit interface normal flux
nHatf_ = nHatfv & Sf;

// Simple expression for curvature
K_ = -fvc::div(nHatf_);

}

As far is I see the gradGamma at the cell centers is interpolated to the cellfaces and named gradGammaf.
This gradGammaf is normalized to unit length and named nHatfv.
nHatfv is "dot producted" with the cell face vector Sf and named nHatf_.
The curvature K_ is than declared as -div(nHatf_).

If nHatf_ is calculated at the cell face how can sigma_*K_ be a volScalarField?

I'm still not getting it. Where am I thinking wrong?

michi January 13, 2009 18:29

Hallo Sebastian, there are di
 
Hallo Sebastian,
there are different versions of fvc::div(). It is possible to have surfaceFields or volFields inside the brackets (). But as far as I see all have volFields as return value. Here we put in the surfaceScalarField nHatf_ and get back the volScalarField K_. It seems right that K_ is a volScalarField because most other variables (like pressure) are also definded as volFields. The mathematics behind the fvc::div() operation can be seen in the Programmer's Guide 2.4.5 / 2.4.10 (look at surfaceIntegrate). Maybe a better reference for what is done here is the PhD thesis of Onno Ubbink (equation 3.47).

Greetings
Michael

sega January 14, 2009 15:41

Ok, Michael, Thank you. It
 
Ok, Michael, Thank you.

It dawns me slowly - Ubbink (3.47) makes it some sense to me.
Since all the gradients in in natfv are explicitly interpolated the fvc::div(nHatf_) is just a tensor operation and does not need any approximation?
And the "outer" fvc::interpolation() is meant for the curvation at the cell center.

But why interpolate it from the neighbours, if can be calculated from -sigma*div(nHatf_) ?

sega January 14, 2009 15:43

Ok, Michael, Thank you. It
 
Ok, Michael, Thank you.

It dawns me slowly - Ubbink (3.47) makes it some sense to me.
Since all the gradients in in natfv are explicitly interpolated the fvc::div(nHatf_) is just a tensor operation and does not need any approximation?
And the "outer" fvc::interpolation() is meant for the curvation at the cell center.

But why interpolate it from the neighbours, if can be calculated from -sigma*div(nHatf_) ?

What happens to the cell volume if the force due to surface tension is treated like a source term?

sega January 14, 2009 15:46

Ok, Michael, Thank you. It
 
Ok, Michael, Thank you.

It dawns me slowly - Ubbink (3.47) makes it some sense to me.
Since all the gradients in in natfv are explicitly interpolated the fvc::div(nHatf_) is just a tensor operation and does not need any approximation?
And the "outer" fvc::interpolation() is meant for the curvation at the cell center.

But why interpolate it from the neighbours, if can be calculated from -sigma*div(nHatf_) ?

It the force due to surface tension is treated like a source term, what happens to the cell volume?

sega January 21, 2009 15:31

Hello again. Thinking about
 
Hello again.

Thinking about the interface compression with phir I wondered about something. If the artificial velocity is directed into the direction of the interface normal, would this lead to a velocity directed to one side of the interface and compress it only on this side?

Is there an illustration of the behaviour of such a scheme to the interface?

S.

sega January 22, 2009 17:19

Well, the answers are getting
 
Well, the answers are getting thinner.

Lets try an easier one:
Using VOF means using blended density and viscosity in the Navier-Stokes-Equations:

rho = gamma*rho1+(1-gamma)*rho2

Is that applicable to the velocity, too?

u = gamma*u1+(1-gamma)*u2

liu January 22, 2009 17:58

Is that applicable to the velo
 
Is that applicable to the velocity, too?
u = gamma*u1+(1-gamma)*u2

That's a good question. But one should be clear about the two-fluid model concept. In two-fluid model (VOF here), only one velocity for both phases. NO u1 and u2.

sega January 22, 2009 18:43

Then, why is Jasak stating thi
 
Then, why is Jasak stating this in his slide

http://powerlab.fsb.hr/ped/kturbo/Op..._12Jan2005.pdf

on page 4?


All times are GMT -4. The time now is 06:36.