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

Artificial high velocities at the interface using interFoam

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

Reply
 
LinkBack Thread Tools Display Modes
Old   February 24, 2011, 15:10
Default Artificial high velocities at the interface using interFoam
  #1
Senior Member
 
Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 8
Arnoldinho is on a distinguished road
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


Attached Files
File Type: txt fvSchemes.txt (1.5 KB, 97 views)
File Type: txt fvSolution.txt (1.8 KB, 64 views)
Arnoldinho is offline   Reply With Quote

Old   February 24, 2011, 17:34
Default
  #2
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
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,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   February 25, 2011, 05:21
Default
  #3
Senior Member
 
Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 8
Arnoldinho is on a distinguished road
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
Arnoldinho is offline   Reply With Quote

Old   February 25, 2011, 23:33
Default
  #4
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by Arnoldinho View Post
The Gauss limitedLinearV 1 seems more reliable and accurate to me.
More accurate, yes. More stable, I am not so sure

Quote:
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   February 28, 2011, 10:52
Default
  #5
Senior Member
 
Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 8
Arnoldinho is on a distinguished road
Hi Alberto,

Quote:
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
Arnoldinho is offline   Reply With Quote

Old   February 28, 2011, 14:47
Default
  #6
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by Arnoldinho View Post
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).

Quote:
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
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   June 6, 2011, 14:23
Default
  #7
Senior Member
 
Pablo
Join Date: Mar 2009
Posts: 102
Rep Power: 8
pablodecastillo is on a distinguished road
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
pablodecastillo is offline   Reply With Quote

Old   June 6, 2011, 14:49
Default
  #8
Senior Member
 
Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 8
Arnoldinho is on a distinguished road
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
Arnoldinho is offline   Reply With Quote

Old   June 6, 2011, 16:54
Default
  #9
Senior Member
 
Pablo
Join Date: Mar 2009
Posts: 102
Rep Power: 8
pablodecastillo is on a distinguished road
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
pablodecastillo is offline   Reply With Quote

Old   June 7, 2011, 02:36
Default
  #10
Senior Member
 
Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 8
Arnoldinho is on a distinguished road
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.
Arnoldinho is offline   Reply With Quote

Old   June 7, 2011, 05:09
Default
  #11
Senior Member
 
Pablo
Join Date: Mar 2009
Posts: 102
Rep Power: 8
pablodecastillo is on a distinguished road
Now it is very clear and easy. Thank you very much.

Pablo
pablodecastillo is offline   Reply With Quote

Old   June 9, 2011, 15:39
Default
  #12
Senior Member
 
Pablo
Join Date: Mar 2009
Posts: 102
Rep Power: 8
pablodecastillo is on a distinguished road
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
pablodecastillo is offline   Reply With Quote

Old   June 13, 2011, 14:28
Default
  #13
New Member
 
Join Date: Jun 2011
Posts: 5
Rep Power: 6
cfd_user2011 is on a distinguished road
Quote:
Originally Posted by Arnoldinho View Post
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
cfd_user2011 is offline   Reply With Quote

Old   June 13, 2011, 17:04
Default
  #14
Senior Member
 
Pablo
Join Date: Mar 2009
Posts: 102
Rep Power: 8
pablodecastillo is on a distinguished road
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
pablodecastillo is offline   Reply With Quote

Old   June 14, 2011, 01:29
Default
  #15
Senior Member
 
Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 8
Arnoldinho is on a distinguished road
Quote:
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?
Arnoldinho is offline   Reply With Quote

Old   June 14, 2011, 06:45
Default
  #16
Senior Member
 
Pablo
Join Date: Mar 2009
Posts: 102
Rep Power: 8
pablodecastillo is on a distinguished road
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
pablodecastillo is offline   Reply With Quote

Old   June 14, 2011, 10:41
Default
  #17
Senior Member
 
Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 8
Arnoldinho is on a distinguished road
Quote:
Originally Posted by pablodecastillo View Post
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
Arnoldinho is offline   Reply With Quote

Old   December 7, 2011, 09:33
Default
  #18
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
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)
Code:
- ghf*fvc::snGrad(rho)
By
Code:
- 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.
Bernhard is offline   Reply With Quote

Old   December 20, 2011, 08:14
Default
  #19
Senior Member
 
Pablo
Join Date: Mar 2009
Posts: 102
Rep Power: 8
pablodecastillo is on a distinguished road
Thanks to share Bernhard, i will try to test ASAP.

Pablo
pablodecastillo is offline   Reply With Quote

Old   February 15, 2012, 12:25
Default
  #20
New Member
 
Join Date: May 2011
Posts: 15
Rep Power: 6
joris.hey is on a distinguished road
Hi, the solution
Quote:
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 :

Quote:
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
joris.hey 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
Wind turbine simulation Saturn CFX 34 October 16, 2014 05:27
InterFoam Artificially High Velocities andersson.j OpenFOAM Running, Solving & CFD 0 February 8, 2011 11:43
RPM in Wind Turbine Pankaj CFX 9 November 23, 2009 05:05
Multicomponent fluid Andrea CFX 2 October 11, 2004 05:12
Replace periodic by inlet-outlet pair lego CFX 3 November 5, 2002 21:09


All times are GMT -4. The time now is 12:41.