CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

VOF method

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree9Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   December 24, 2008, 07:54
Default Hello Sebastian, The inlet
  #21
Senior Member
 
Mark Couwenberg
Join Date: Mar 2009
Location: Netherlands
Posts: 130
Rep Power: 8
markc is on a distinguished road
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,
markc is offline   Reply With Quote

Old   December 25, 2008, 02:02
Default Has anyone any comment on the
  #22
Member
 
Cem Albukrek
Join Date: Mar 2009
Posts: 50
Rep Power: 8
albcem is on a distinguished road
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
albcem is offline   Reply With Quote

Old   December 29, 2008, 06:35
Default Albcem, One thing which com
  #23
Senior Member
 
Mark Couwenberg
Join Date: Mar 2009
Location: Netherlands
Posts: 130
Rep Power: 8
markc is on a distinguished road
Albcem,

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

Brgds,

Mark
markc is offline   Reply With Quote

Old   December 29, 2008, 15:46
Default Hello Mark, Thanks for the
  #24
Member
 
Cem Albukrek
Join Date: Mar 2009
Posts: 50
Rep Power: 8
albcem is on a distinguished road
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
albcem is offline   Reply With Quote

Old   December 30, 2008, 08:00
Default Hello Albcem, I do not know
  #25
Senior Member
 
Mark Couwenberg
Join Date: Mar 2009
Location: Netherlands
Posts: 130
Rep Power: 8
markc is on a distinguished road
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 is offline   Reply With Quote

Old   December 30, 2008, 08:02
Default Hello Albcem, I do not know
  #26
Senior Member
 
Mark Couwenberg
Join Date: Mar 2009
Location: Netherlands
Posts: 130
Rep Power: 8
markc is on a distinguished road
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.

pd
U
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 is offline   Reply With Quote

Old   January 4, 2009, 08:56
Default Hello. Happy New Year to al
  #27
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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.
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 6, 2009, 04:56
Default Let me put my question in anot
  #28
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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?
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 12, 2009, 05:31
Default Hello Sebastian, I am also
  #29
New Member
 
HCC
Join Date: Mar 2009
Posts: 1
Rep Power: 0
amghsu is on a distinguished road
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
amghsu is offline   Reply With Quote

Old   January 12, 2009, 11:45
Default Thank you Jim. Then, how is
  #30
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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.
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 13, 2009, 08:07
Default Hallo Sebastian, As far as
  #31
New Member
 
Michael Wurst
Join Date: Mar 2009
Location: Munich, Germany
Posts: 6
Rep Power: 8
michi is on a distinguished road
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
michi is offline   Reply With Quote

Old   January 13, 2009, 15:11
Default Thank you Michael. Ok, sigm
  #32
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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?
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 13, 2009, 18:29
Default Hallo Sebastian, there are di
  #33
New Member
 
Michael Wurst
Join Date: Mar 2009
Location: Munich, Germany
Posts: 6
Rep Power: 8
michi is on a distinguished road
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
michi is offline   Reply With Quote

Old   January 14, 2009, 15:41
Default Ok, Michael, Thank you. It
  #34
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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_) ?
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 14, 2009, 15:43
Default Ok, Michael, Thank you. It
  #35
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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?
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 14, 2009, 15:46
Default Ok, Michael, Thank you. It
  #36
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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?
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 21, 2009, 15:31
Default Hello again. Thinking about
  #37
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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.
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 22, 2009, 17:19
Default Well, the answers are getting
  #38
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
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
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   January 22, 2009, 17:58
Default Is that applicable to the velo
  #39
liu
Senior Member
 
Xiaofeng Liu
Join Date: Mar 2009
Location: State College, PA, USA
Posts: 118
Rep Power: 8
liu is on a distinguished road
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.
__________________
Xiaofeng Liu, Ph.D., P.E.,
Assistant Professor
Department of Civil and Environmental Engineering
Penn State University
223B Sackett Building
University Park, PA 16802


Web: http://water.engr.psu.edu/liu/
liu is offline   Reply With Quote

Old   January 22, 2009, 18:43
Default Then, why is Jasak stating thi
  #40
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11
sega is on a distinguished road
Then, why is Jasak stating this in his slide

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

on page 4?
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Comparing between the Fractional step method and the SIMPLE method ghlee Main CFD Forum 1 April 10, 2012 16:59
Finite volume method vs finite difference method? superfool Main CFD Forum 4 October 21, 2006 14:37
Pressure Correction method, Finite Volume Method Seeker01 Main CFD Forum 2 January 13, 2003 03:49
hess-smith method and fvm method yangqing FLUENT 0 March 20, 2002 20:25
Projection method and Block-off method Leo Main CFD Forum 0 June 14, 2001 07:22


All times are GMT -4. The time now is 04:21.