CFD Online URL
[Sponsors]
Home > Forums > OpenFOAM

Non-physcial solution in the case of triangular mesh from interFoam

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

Reply
 
LinkBack Thread Tools Display Modes
Old   October 7, 2011, 01:30
Default Non-physcial solution in the case of triangular mesh from interFoam
  #1
ata
Senior Member
 
ata's Avatar
 
ata kamyabi
Join Date: Aug 2009
Location: Kerman
Posts: 322
Rep Power: 8
ata is on a distinguished road
Hi
I have a problem with interFoam. When I run my problem on a
quadrangular mesh every thing is OK but when I run my problem on a triangular mesh I have some non-physical big velocities. I know this problem is due to high density ratio (1000/1) and triangular mesh because when I change the densities to (1000/1000) or when I set the gravity to (0 0 0) the parasitic current eliminates.
Can any one tell me how I can get a good solution in the case of triangular mesh too?
Best regards
Ata
ata is offline   Reply With Quote

Old   October 8, 2011, 02:38
Default
  #2
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,890
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
You should post the schemes you are using. As a general suggestion, use bounded schemes, and limit gradients too.

Best,
__________________
Alberto

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.
GeekoCFD 32bit - The 32bit edition of GeekoCFD.
GeekoCFD text mode - A smaller version of GeekoCFD, text-mode only, with only OpenFOAM. Available in a variety of virtual formats.
alberto is offline   Reply With Quote

Old   October 8, 2011, 03:47
Default What scheme?
  #3
ata
Senior Member
 
ata's Avatar
 
ata kamyabi
Join Date: Aug 2009
Location: Kerman
Posts: 322
Rep Power: 8
ata is on a distinguished road
Hi Dear Alberto
These are schemes that I used:
ddtSchemes
{
default Euler;
}

gradSchemes
{
default Gauss linear;
}

divSchemes
{
div(rho*phi,U) Gauss upwind;
div(phi,alpha) Gauss vanLeer;
div(phirb,alpha) Gauss interfaceCompression;
}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

What is your suggestion?
Thank you very much
Best regards

Ata
ata is offline   Reply With Quote

Old   October 8, 2011, 03:57
Default
  #4
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,890
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
Try

Code:
gradSchemes
{
       default         cellLimited leastSquares 1; 
}
Also, check the mesh with checkMesh to see if it presents some problem.

Best,
__________________
Alberto

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.
GeekoCFD 32bit - The 32bit edition of GeekoCFD.
GeekoCFD text mode - A smaller version of GeekoCFD, text-mode only, with only OpenFOAM. Available in a variety of virtual formats.
alberto is offline   Reply With Quote

Old   October 8, 2011, 04:31
Default Non-physcial solution in the case of triangular mesh from interFoam
  #5
ata
Senior Member
 
ata's Avatar
 
ata kamyabi
Join Date: Aug 2009
Location: Kerman
Posts: 322
Rep Power: 8
ata is on a distinguished road
Hi Dear Alberto
Thank you very much for your quick reply.
I used your suggested scheme but the problem did not solve?
Do you have any other suggestion?
Best regards

Ata
ata is offline   Reply With Quote

Old   October 8, 2011, 15:06
Default
  #6
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,890
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
Could you show a) a picture of the mesh, b) contour plot of alpha, and c) contour plot of the velocity, showing where the problem happens?

Best,
__________________
Alberto

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.
GeekoCFD 32bit - The 32bit edition of GeekoCFD.
GeekoCFD text mode - A smaller version of GeekoCFD, text-mode only, with only OpenFOAM. Available in a variety of virtual formats.
alberto is offline   Reply With Quote

Old   October 9, 2011, 02:33
Default Non-physcial solution in the case of triangular mesh from interFoam
  #7
ata
Senior Member
 
ata's Avatar
 
ata kamyabi
Join Date: Aug 2009
Location: Kerman
Posts: 322
Rep Power: 8
ata is on a distinguished road
Hi Dear Alberto
Please see the attachments
Thank you very much
Best regards

Ata
Attached Images
File Type: jpg alpha1.jpg (97.5 KB, 101 views)
File Type: jpg U.jpg (97.3 KB, 91 views)
ata is offline   Reply With Quote

Old   October 9, 2011, 02:55
Default
  #8
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,890
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
It seems to me the problem is at the junction between the rectangular and the inclined part of the domain too, where the mesh is still hexahedral. The velocity peaks are all in hex cells too.

Your mesh seems fine to me, however I would play a bit with the schemes:

Code:
divSchemes
{
       div(rho*phi,U)  Gauss limitedLinearV 1;
       div(phi,alpha)  Gauss vanLeer;
       div(phirb,alpha) Gauss interfaceCompression;
}

laplacianSchemes
{
       default         Gauss linear limited 0.33;
}

snGradSchemes
{
   default limited 0.33;
}
Also, turn on the momentum predictor, and check if you have convergence difficulties for the equations.

Best,
__________________
Alberto

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.
GeekoCFD 32bit - The 32bit edition of GeekoCFD.
GeekoCFD text mode - A smaller version of GeekoCFD, text-mode only, with only OpenFOAM. Available in a variety of virtual formats.
alberto is offline   Reply With Quote

Old   October 10, 2011, 03:10
Default Non-physcial solution in the case of triangular mesh from interFoam
  #9
ata
Senior Member
 
ata's Avatar
 
ata kamyabi
Join Date: Aug 2009
Location: Kerman
Posts: 322
Rep Power: 8
ata is on a distinguished road
Hi Dear Alberto and thank you very much for your help.
I used your proposed schemes and it seems that alpha is good but the velocity is non-physical. Please see the attachments.
Do you have any other proposition.
Thank you very much.
Best regards

Ata
Attached Images
File Type: jpg alpha1.jpg (53.9 KB, 61 views)
File Type: jpg U.jpg (61.1 KB, 62 views)
ata is offline   Reply With Quote

Old   October 10, 2011, 11:49
Default
  #10
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,890
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
How is the solution behaving? Is it stable in time, or does the problem with alpha come back? Since the velocity is unphysical, I would expect alpha to be too.

What about the convergence history?

Best,
__________________
Alberto

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.
GeekoCFD 32bit - The 32bit edition of GeekoCFD.
GeekoCFD text mode - A smaller version of GeekoCFD, text-mode only, with only OpenFOAM. Available in a variety of virtual formats.
alberto is offline   Reply With Quote

Old   October 10, 2011, 12:07
Default
  #11
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Rotterdam, The Netherlands
Posts: 1,553
Rep Power: 23
ngj will become famous soon enoughngj will become famous soon enough
Hi Ata

I have been having a lot of "fun" with interFoam/MULES over the last couple of years, and what you report is a well known problem.

We, in our group, have tried a lot of different numerical schemes, etc, however, in vain. The problem lies in the evaluation of snGrad(rho), which is not evaluated in a manner, which is consistent with the underlying physics. This means, that you will not be able to run the "bucket of water" simulation with interFoam/MULES on a non-orthogonal, hexahedral mesh, or an orthogonal mesh with a gravity vector, which is not aligned with the mesh lines.

Do a search for "Wemmenhove" and "gravity-consistent" on Google, and you will find a proposition to solve this problem, however, the solution is not compatible with MULES, because it solve the hyperbolic differential equation for alpha, whereas the proposed solutions is meant for the PLIC approach for free surface flows.

Kind regards,

Niels

P.S. Are you looking at wave run-up on a breakwater?
ngj is offline   Reply With Quote

Old   October 10, 2011, 12:25
Default
  #12
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,890
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
Interesting. Thanks Niels :-)
__________________
Alberto

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.
GeekoCFD 32bit - The 32bit edition of GeekoCFD.
GeekoCFD text mode - A smaller version of GeekoCFD, text-mode only, with only OpenFOAM. Available in a variety of virtual formats.
alberto is offline   Reply With Quote

Old   October 10, 2011, 21:41
Default
  #13
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 355
Rep Power: 9
arjun is on a distinguished road
Quote:
Originally Posted by ngj View Post
Hi Ata

I have been having a lot of "fun" with interFoam/MULES over the last couple of years, and what you report is a well known problem.

We, in our group, have tried a lot of different numerical schemes, etc, however, in vain. The problem lies in the evaluation of snGrad(rho), which is not evaluated in a manner, which is consistent with the underlying physics. This means, that you will not be able to run the "bucket of water" simulation with interFoam/MULES on a non-orthogonal, hexahedral mesh, or an orthogonal mesh with a gravity vector, which is not aligned with the mesh lines.

Do a search for "Wemmenhove" and "gravity-consistent" on Google, and you will find a proposition to solve this problem, however, the solution is not compatible with MULES, because it solve the hyperbolic differential equation for alpha, whereas the proposed solutions is meant for the PLIC approach for free surface flows.

Kind regards,

Niels

P.S. Are you looking at wave run-up on a breakwater?

Hi thanks for pointing 'gravity consistent' part.
I have been testing my code (NOT OPENFOAM) with VOF on different types of meshes and I have yet to face this spurious velocity problem even on far worse meshes than what Ata showed. Off course I do lots of things different than openFOAM so it could be any of the thing.

One of the thing that I do differently is that I add gravity as (density * volume * g ) from the volume of cell and density of cell rather than calculation gravity source from integration from control volume faces (i think which is usual practice and assumed to be more stable). So in the way I do, gravity vector is not affected by type of mesh but by 'mass in cell'.
arjun is offline   Reply With Quote

Old   October 10, 2011, 21:53
Default
  #14
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,890
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
Why is the snGrad(rho) not consistent with the physics, and how should it be calculated?
__________________
Alberto

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.
GeekoCFD 32bit - The 32bit edition of GeekoCFD.
GeekoCFD text mode - A smaller version of GeekoCFD, text-mode only, with only OpenFOAM. Available in a variety of virtual formats.
alberto is offline   Reply With Quote

Old   October 10, 2011, 21:56
Default
  #15
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,890
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by arjun View Post
One of the thing that I do differently is that I add gravity as (density * volume * g ) from the volume of cell and density of cell rather than calculation gravity source from integration from control volume faces (i think which is usual practice and assumed to be more stable). So in the way I do, gravity vector is not affected by type of mesh but by 'mass in cell'.
Sounds easy to verify. :-)

However, how do you manage the Rhie-Chow interpolation in cases with strong density gradients then, since your body force does not seem to enter the interpolation formula? The body force should tend to destabilize the solution...
__________________
Alberto

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.
GeekoCFD 32bit - The 32bit edition of GeekoCFD.
GeekoCFD text mode - A smaller version of GeekoCFD, text-mode only, with only OpenFOAM. Available in a variety of virtual formats.
alberto is offline   Reply With Quote

Old   October 10, 2011, 22:05
Default
  #16
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 355
Rep Power: 9
arjun is on a distinguished road
Quote:
Originally Posted by alberto View Post
Why is the snGrad(rho) not consistent with the physics, and how should it be calculated?

Edited to remove, it seems something wrong with what I wrote. I will rewrite once I confirm it.

Sorry.
arjun is offline   Reply With Quote

Old   October 10, 2011, 22:09
Default
  #17
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 355
Rep Power: 9
arjun is on a distinguished road
Quote:
Originally Posted by alberto View Post
Sounds easy to verify. :-)

However, how do you manage the Rhie-Chow interpolation in cases with strong density gradients then, since your body force does not seem to enter the interpolation formula? The body force should tend to destabilize the solution...

This is true. However so far I have only tried water air system and did not notice problems.

NOTE: I do lots of other things to keep solver stable so it might be any one of the things that keeps solver stable and I did not notice it.
arjun is offline   Reply With Quote

Old   October 11, 2011, 03:16
Default
  #18
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Rotterdam, The Netherlands
Posts: 1,553
Rep Power: 23
ngj will become famous soon enoughngj will become famous soon enough
Hi Alberto

Their analysis with respect to the handling of the density on the faces is as follows:

Consider four computational cells

A B | 0.0 0.3333
C D | 0.6666 1.0000

where the letters are cell naming, and the numbers are the void ratio. The surface is a straight line, the mesh orthogonal in x-y, the gravity vector is (-1 -2), and U = 0. They perform an interpolation of rho onto the faces and uses these values to find the hydrostatic pressure in D based on either going from A->C->D or A->B->D.

Result: The two pressures in D differs!

The solution is to do a special interpolation routine, which depends on the unique distance from the cell centers to the reconstructed interface in a PLIC fashion. As MULES will smear the interface over more than one cell, such a distance will be difficult to define.

Their approach has removed the spurious currents at the interface. Another formulation of the same is that the face values of rho should fulfill

curl( rho g ) = 0

according to Wemmehove et al (do not have the full refernce here), however, I am not happy about this formulation, as it does not seem to be mathematically rigid, because they do a operator change over the interface, where "grad rho" is not continuous.

Hope this clarified the details.

Niels
ngj is offline   Reply With Quote

Old   October 11, 2011, 03:55
Default How we can solve the problem?
  #19
ata
Senior Member
 
ata's Avatar
 
ata kamyabi
Join Date: Aug 2009
Location: Kerman
Posts: 322
Rep Power: 8
ata is on a distinguished road
Hi Dear Friends
Very nice discussions.
I think that this problem can be solved in OpenFOAM by to ways:
1. Smoothed rho. I think we can use a normalized kernel function to smooth rho. I think this solves the problem. Because when I set gravity to zero or set the density ratio near one the problem has been solved. Therefore it seems that the problem is due to the discontinuity in the density near the surface and this can be solved using the smoother in order to grad(rho) calculate stably.
2. We can use pressure(p) instead of modified pressure(p_rgh) in order to grad(rho) eliminated from the momentum equation and in the wall boundaries we can use grad(p)=grad(rho*g&h) instead of grad(p_rgh)=0.

What is your opinion about these?
Best regards

Ata
ata is offline   Reply With Quote

Old   October 11, 2011, 06:35
Default
  #20
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Rotterdam, The Netherlands
Posts: 1,553
Rep Power: 23
ngj will become famous soon enoughngj will become famous soon enough
Hi Ata

With respect to "2", then this was tried in OF-1.6, however, they have returned to dynamic/excess pressure as of 1.7. The problem is simply that using total pressure you get very poor results.

Example:

Wave breaking on a coast with excess pressure gives the correct variation in the mean water level across the profile, i.e. set-down seaward of the breakpoint and subsequent set-up shoreward of the breakpoint.
If you use the formulation with total pressure, you will only get this behaviour under the severe constraint that the mesh lines are horizontal/vertical. If they are bending to become parallel with the bed, then you will experience a set-up from the beginning of the sloping beach profile even this point is seaward of the break point.

/ Niels
ngj 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
snappyHexMesh won't work - zeros everywhere! sc298 OpenFOAM Native Meshers: snappyHexMesh and Others 2 March 27, 2011 22:11
IcoFoam parallel woes msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 03:58
Discussion about Mesh independant solution Seb Main CFD Forum 13 May 22, 2001 14:37
CFL Condition Matt Umbel Main CFD Forum 14 January 12, 2001 15:34
unstructured vs. structured grids Frank Muldoon Main CFD Forum 1 January 5, 1999 11:09


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