[markc]

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]

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]

Albcem,

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

Brgds,

Mark

[albcem]

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]

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]

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

[sega]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

Then, why is Jasak stating this in his slide

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

on page 4?