# Artificial high velocities at the interface using interFoam

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

February 24, 2011, 15:10
Artificial high velocities at the interface using interFoam
#1
Senior Member

Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 10
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
 fvSchemes.txt (1.5 KB, 172 views) fvSolution.txt (1.8 KB, 123 views)

 February 24, 2011, 17:34 #2 Senior Member   Alberto Passalacqua Join Date: Mar 2009 Location: Ames, Iowa, United States Posts: 1,910 Rep Power: 28 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, vsammartano and mizzou like this. __________________ Alberto Passalacqua GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM - An open-source implementation of quadrature-based moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using.

 February 25, 2011, 05:21 #3 Senior Member   Arne Stahlmann Join Date: Nov 2009 Location: Hanover, Germany Posts: 209 Rep Power: 10 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

February 25, 2011, 23:33
#4
Senior Member

Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,910
Rep Power: 28
Quote:
 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

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 in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.

February 28, 2011, 10:52
#5
Senior Member

Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 10
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

February 28, 2011, 14:47
#6
Senior Member

Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,910
Rep Power: 28
Quote:
 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).

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 in both physical and virtual formats.
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.

 June 6, 2011, 14:23 #7 Senior Member   Pablo Join Date: Mar 2009 Posts: 102 Rep Power: 10 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

 June 6, 2011, 14:49 #8 Senior Member   Arne Stahlmann Join Date: Nov 2009 Location: Hanover, Germany Posts: 209 Rep Power: 10 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

 June 6, 2011, 16:54 #9 Senior Member   Pablo Join Date: Mar 2009 Posts: 102 Rep Power: 10 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

 June 7, 2011, 02:36 #10 Senior Member   Arne Stahlmann Join Date: Nov 2009 Location: Hanover, Germany Posts: 209 Rep Power: 10 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.

 June 7, 2011, 05:09 #11 Senior Member   Pablo Join Date: Mar 2009 Posts: 102 Rep Power: 10 Now it is very clear and easy. Thank you very much. Pablo

 June 9, 2011, 15:39 #12 Senior Member   Pablo Join Date: Mar 2009 Posts: 102 Rep Power: 10 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

June 13, 2011, 14:28
#13
New Member

Join Date: Jun 2011
Posts: 5
Rep Power: 7
Quote:
 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

 June 13, 2011, 17:04 #14 Senior Member   Pablo Join Date: Mar 2009 Posts: 102 Rep Power: 10 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

June 14, 2011, 01:29
#15
Senior Member

Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 10
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?

 June 14, 2011, 06:45 #16 Senior Member   Pablo Join Date: Mar 2009 Posts: 102 Rep Power: 10 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 vsammartano likes this.

June 14, 2011, 10:41
#17
Senior Member

Arne Stahlmann
Join Date: Nov 2009
Location: Hanover, Germany
Posts: 209
Rep Power: 10
Quote:
 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

 December 7, 2011, 09:33 #18 Senior Member   Bernhard Join Date: Sep 2009 Location: Delft Posts: 790 Rep Power: 14 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. vsammartano and mo_na like this.

 December 20, 2011, 08:14 #19 Senior Member   Pablo Join Date: Mar 2009 Posts: 102 Rep Power: 10 Thanks to share Bernhard, i will try to test ASAP. Pablo

February 15, 2012, 12:25
#20
New Member

Join Date: May 2011
Posts: 15
Rep Power: 7
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::fvPatchField, Foam::volMesh>::.Foam:imensionedField , Foam::volMesh>::.Foam::Field >::.Foam::List >::.Foam::UList:perator[] [with T = Foam::Vector](celli) = 0.0’ /opt/openfoam210/src/OpenFOAM/lnInclude/Vector.H:61: note: candidates are: Foam::Vector& Foam::Vector:perator=(const Foam::Vector&) /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 ?

Joris

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Saturn CFX 45 February 8, 2016 05:42 andersson.j OpenFOAM Running, Solving & CFD 0 February 8, 2011 11:43 Pankaj CFX 9 November 23, 2009 05:05 Andrea CFX 2 October 11, 2004 05:12 lego CFX 3 November 5, 2002 21:09

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