[Arnoldinho]

Hi all,

I'm simulation flow of water around an object using the interFoam solver and k-Omega SST turbulence model. At the moment, I'm struggling with some artificial high velocities which occur at the interface in the upper air zone. At the beginning of the simulation, some air interface cells have velocities of up to 20 times the normal maximum velocities in my water zone, which really reduce my time steps or let the simulation crash.

For the meshing, I used blockMesh and sHM including layers. checkMesh gives

...
Minumum face area = 3.3007033e-07. Maximum face area = 0.0026526524. Face area magnitudes OK.
Min volume = 4.735839e-10. Max volume = 5.2955194e-05. Total volume = 1.6946638. Cell volumes OK.
Mesh non-orthogonality Max: 74.576967 average: 6.9868146
*Number of severely non-orthogonal faces: 70.
Non-orthogonality check OK.
<<Writing 70 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 2.4098523 OK.

Mesh OK.


As I'm not interested in what is happening on the air part, I already changed the interFoam solver to ignore the convective term on the air side. This improved the simulation, but it's still not running very well.

Can someone give me a hint on how to modify my solver settings or improve my mesh?
A picture and the settings are attached/given below.


BTW: Is it somehow possible to 'sharpen' the interface by changing the settings? Or is this only possible by refining the mesh at the interface?

Arne

[alberto]

Hi,

try using a limited scheme for div(rho*phi,U). You are currently using central differences (linear), which is not necessarily stable (is your cell Re < 2)?

Maybe try

div(rho*phi,U) Gauss limitedLinearV 1;

which should be quite accurate, if it works. If you need something more robust:

div(rho*phi,U) Gauss linearUpwindV cellLimited Gauss linear 1;

Best,

[Arnoldinho]

Thanks Alberto,

I guess you saved a lot of my time! I tried both schemes, and both simulations are working. The Gauss limitedLinearV 1 seems more reliable and accurate to me.

Nevertheless, artificial velocities at the interface still occure, but are much lower. They are "generated" directly at the simulation start, and never disappear (the latter maybe because the convective term on the air side is ignored?).
Do you have an idea where they are coming from and, what is more important, how I can completely avoid them? Any hints on which direction I could have a closer look?

Initial and BC: For the initial state, the whole water domain has a unique velocity given by setFields. Air has zero velocity. At the inlet, velocity is set as constant at the water part using groovyBC and zero for the air part during the simulation.
For the initial state, I also changed this to a smoother transition for alpha and U giving the interface cells half values (0.5 and half water velocity). This did not solve it as well...

Arne

[alberto]

Originally Posted by Arnoldinho

The Gauss limitedLinearV 1 seems more reliable and accurate to me.

More accurate, yes. More stable, I am not so sure

Nevertheless, artificial velocities at the interface still occure, but are much lower. They are "generated" directly at the simulation start, and never disappear (the latter maybe because the convective term on the air side is ignored?).

Try with the unmodified code to see if this still happens.

Best,
Alberto

[Arnoldinho]

Hi Alberto,

Try with the unmodified code to see if this still happens.

Although I changed and primarily tried different fvSolution settings for a comparison of PCG and GAMG solver, I also did calculations with modified and unmodified interFoam code. For the 'original' one, these artificial velocities occur as well, in the near-field of the mediums air, water and solid structure (with layers).
I noticed that some small cells in the boundary layer around the structure at the water/air transition have a very high value for k, in an order of 10 times the normal value for cells around the structure.
But I don't exactly know what to do with them...

Another question: Do you also have a suggestion for the 'right' fvSolution settings (PCG, GAMG) in terms of speed and accuracy for interFoam and around 2 million cells?

Arne

[alberto]

Originally Posted by Arnoldinho

Hi Alberto,

Although I changed and primarily tried different fvSolution settings for a comparison of PCG and GAMG solver, I also did calculations with modified and unmodified interFoam code. For the 'original' one, these artificial velocities occur as well, in the near-field of the mediums air, water and solid structure (with layers).
I noticed that some small cells in the boundary layer around the structure at the water/air transition have a very high value for k, in an order of 10 times the normal value for cells around the structure.

You might want to see if it still happens without turbulence model (laminar).

Another question: Do you also have a suggestion for the 'right' fvSolution settings (PCG, GAMG) in terms of speed and accuracy for interFoam and around 2 million cells?

It is more a question of speed than accuracy. I tent to use GAMG for pressure, and CG methods for the rest, which is probably the default setting in the tutorials. If GAMG gives you troubles, use PCG for pressure too: it will take more iterations, but it seems to be a little more robust.

Best,
Alberto

[pablodecastillo]

Hello Arne,

I am interested in ignore the convective term on the air side, can you let me known how did you do?, is it improving your simulation? (i mean if you are avoiding very small time steps?

Advanced thanks

Pablo

[Arnoldinho]

Pablo, it did not really improve the simulation with regard to smaller time steps and high velocities in the air phase.

Depending on what you want to do and if the air phase is not of great importance for you, you could try modifying the solver (e.g. interFoam, copied and compiled to my_interFoam) and set the velocities in alpha1 (all cels smaller than a value of lets say 0.05) to zero every timestep. For my case, this saved computational time. Nevertheless, in case of surface waves, you have to be careful with regard to wave damping.

Ignoring the convective term would also be done in the solver within the U equation.

Arne

[pablodecastillo]

Hello Arne,

Set velocities 0 when alpha1<0.05 is clear for me, it is like a BC in every timestep, but modify the convertive term, i can not see too clear to implement.

We have
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
So we can modify U or rhoPhi, but how??, can you be more explicit?

Advanced thanks

[Arnoldinho]

Please have a look at the slides from Eric Paterson: http://www.google.de/url?sa=t&source...bWpkzA&cad=rja

Here, gamma/alpha1 is explicitly included in the equation and therefore the term becomes zero in air phase.

[pablodecastillo]

Now it is very clear and easy. Thank you very much.

Pablo

[pablodecastillo]

Hi Arne,

My version from interFoam is rotational and inviscid, and i modify with ;

fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ alpha1*fvm::div(rhoPhi, U)

but surprised that i got unstable solution. Any idea? if i need to modify another piece of code?

Pablo

[cfd_user2011]

Originally Posted by Arnoldinho

Pablo, it did not really improve the simulation with regard to smaller time steps and high velocities in the air phase.

Depending on what you want to do and if the air phase is not of great importance for you, you could try modifying the solver (e.g. interFoam, copied and compiled to my_interFoam) and set the velocities in alpha1 (all cels smaller than a value of lets say 0.05) to zero every timestep. For my case, this saved computational time. Nevertheless, in case of surface waves, you have to be careful with regard to wave damping.

Ignoring the convective term would also be done in the solver within the U equation.

Arne

Hello,
I am trying to set the velocities to 0 in the air phase as you have mentioned. could you please share how did you do that? I am new to OpenFoam.
Thanks

[pablodecastillo]

You can try something like this:

forAll(U, celli)
{
if (alpha1[celli] < 0.01) {
U[celli] = 0.0;
}

}

Info <<"/////////////////////// update U Air ///////////////////////////////// " << nl;

You can write in a file like mycorrectAir.H, and call from UEqn.H, after convective term.

Pablo

[Arnoldinho]

fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ alpha1*fvm::div(rhoPhi, U)

but surprised that i got unstable solution. Any idea? if i need to modify another piece of code?

That's what I tried as well, but also got unstable results, depending on the mesh. Therefore I did not use it.

@ cfd_user2011: Concerning setting the velocity in air to zero: You have to modify the (interFoam) solver, like Pablo said. Please have a look at this tutorial http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam to see how a copy of the solver can be modified and compiled.

@ Pablo: I'm encountering that the deltaT significantly drops (by a factor 2) durin the whole simulation (using modified solver), compared to normal interFoam. Do you get the same?

[pablodecastillo]

Hi Arne,

I am not using exactly 0 for velocities at the air, i think that is too agressive. I am relaxing only the air phase and i am getting nice results, i mean i am running mesh that before was stopped because unstabilities at the air.

About alpha1*fvm::div(rhoPhi, U) the idea looks right but maybe the mesh , allways i got unstable.

Pablo

[Arnoldinho]

Originally Posted by pablodecastillo

I am not using exactly 0 for velocities at the air, i think that is too agressive.

Yes, that might be true. Setting the velocity to 0.7 times computed velocity results in a speedup of around 15 percent in my case...

Arne

[Bernhard]

I'd like to add something to this thread, since I encoutered similar problems for my case. This problems disappeared with reducing the density ratio in my case from 1000 to 10.

Since I don't like the clipping in the air velocity I was able to get better results by applying the argument of Brackbill ( http://www.sciencedirect.com/science...2199919290240Y ) for surface tension, on the additional source term in the Navier-Stokes equations resutling for the reformulation of the pressure.

Basically what I did, was replacing this (in pEqn.H and UEqn.H of interFoam)

- ghf*fvc::snGrad(rho)

By

- ghf*fvc::snGrad(rho)*2.0*fvc::interpolate(alpha1)

The argument by Brackbill is summarized as follows: as the interface thickness goes to zero, you can multiply the body force by a function that is 1 at the interface. Since \grad\rho is zero except at the interface and 2alpha1 is 1 there, I think this is justified (please correct me if I'm wrong)

For me it solved the issues with the high air velocities, and this implementation is at least more justified than just clipping U.

[pablodecastillo]

Thanks to share Bernhard, i will try to test ASAP.

Pablo

[joris.hey]

Hi, the solution

forAll(U, celli)
{
if (alpha1[celli] < 0.01) {
U[celli] = 0.0;
}

}

Info <<"/////////////////////// update U Air ///////////////////////////////// " << nl;

is not compiling for me... I got this answer :

Uair.H:5: error: no match for ‘operator=’ in ‘U.Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>::<anonymous>.Foam:imensionedField <Foam::Vector<double>, Foam::volMesh>::<anonymous>.Foam::Field<Foam::Vect or<double> >::<anonymous>.Foam::List<Foam::Vector<double> >::<anonymous>.Foam::UList<T>:perator[] [with T = Foam::Vector<double>](celli) = 0.0’
/opt/openfoam210/src/OpenFOAM/lnInclude/Vector.H:61: note: candidates are: Foam::Vector<double>& Foam::Vector<double>:perator=(const Foam::Vector<double>&)
/opt/openfoam210/src/finiteVolume/lnInclude/readTimeControls.H:38: warning: unused variable ‘maxDeltaT’
make: *** [Make/linuxGccDPOpt/OpenChannelFoam.o] Error 1

It seems that the =0.0 is not working for a vector class... How should I write it better ?

thanks in advance !
Joris