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

potentialFoam and Tets

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

Reply
 
LinkBack Thread Tools Display Modes
Old   March 26, 2010, 12:22
Default potentialFoam and Tets
  #1
New Member
 
Scott Hogan
Join Date: Mar 2010
Location: Washington State, USA
Posts: 6
Rep Power: 6
shogan50 is on a distinguished road
I'm playing with a wedge geometry which is meshed with tets. (netgen 1D..3D) in Salome. At the sharp corner of the wedge at the outlet, potential foam generates very large magnitude velocities. I'm beginning to suspect it is due to the tet mesh. I've played with various BC's. The phenomenon seems to occur at both inlet and outlet if the velocity isn't completely constrained.
Suggestions?
My BC's:

inlet:
p:zeroGradient
U0 0 -100);
p: 1e5
U:
outlet
{
type inletOutlet;
value uniform (0 0 0);
inletValue uniform ( 0 0 0);

}
shogan50 is offline   Reply With Quote

Old   March 26, 2010, 16:14
Default more info
  #2
New Member
 
Scott Hogan
Join Date: Mar 2010
Location: Washington State, USA
Posts: 6
Rep Power: 6
shogan50 is on a distinguished road
I modeled this as a half pipe with symmetryPlane. The outlet has the same phenomenon - U values on the order of 70000 in all directions. This partially rules out my previous theory - that it has to do with the point in the wedge. Suggestions?
Attached Images
File Type: png PotentialFoamBCIssue.png (49.7 KB, 77 views)

Last edited by shogan50; March 26, 2010 at 16:15. Reason: More info
shogan50 is offline   Reply With Quote

Old   March 26, 2010, 20:16
Default
  #3
Senior Member
 
sega's Avatar
 
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 10
sega is on a distinguished road
I made bad experience with tet-meshes, but they are not uncommon.

If nothing is wrong with your boundary condition itself you could try to increase the number of non-orthogonal-correctors.
This is the entry in fvSolution called nNonOrthogonalCorrectors.

Increasing this entry to much may result in a severe increase of computational time!
But with the right value you will be able to overcome the problems concerning correct flux calculations when dealing with tet-elements.

Unfortunately there is no guidline how to choose the right value. I suppose you will have to try ...
__________________
Schrödingers wife: "What did you do to the cat? It's half dead!"
sega is offline   Reply With Quote

Old   March 28, 2010, 00:58
Default Interesting
  #4
New Member
 
Scott Hogan
Join Date: Mar 2010
Location: Washington State, USA
Posts: 6
Rep Power: 6
shogan50 is on a distinguished road
Thanks very much for the help. Apparently this falls under the SIMPLE heading in the fvSolution file. It was set to zero and increasing it decreases the phenomenon. At 3 it starts to look good. 10 seems to approach an error asymptote. Below is an output snippet.

I was also confused by functionality in paraview. I had tried this previously. Paraview wasn't automatically changing the scaling, so I didn't notice the change.




Calculating potential flow
DICPCG: Solving for p, Initial residual = 1, Final residual = 0.000951584, No Iterations 61
DICPCG: Solving for p, Initial residual = 0.909937, Final residual = 0.000886357, No Iterations 43
DICPCG: Solving for p, Initial residual = 0.775395, Final residual = 0.000722846, No Iterations 54
DICPCG: Solving for p, Initial residual = 0.847058, Final residual = 0.000837037, No Iterations 39
DICPCG: Solving for p, Initial residual = 0.758512, Final residual = 0.000746753, No Iterations 53
DICPCG: Solving for p, Initial residual = 0.897108, Final residual = 0.000881784, No Iterations 34
DICPCG: Solving for p, Initial residual = 0.89887, Final residual = 0.000770348, No Iterations 55
DICPCG: Solving for p, Initial residual = 0.870719, Final residual = 0.000765984, No Iterations 32
DICPCG: Solving for p, Initial residual = 0.638767, Final residual = 0.000635883, No Iterations 53
DICPCG: Solving for p, Initial residual = 0.26727, Final residual = 0.000256407, No Iterations 32
DICPCG: Solving for p, Initial residual = 0.080843, Final residual = 7.82863e-05, No Iterations 51
continuity error = 0.0420761
Interpolated U error = 1.70791e-07
ExecutionTime = 1.62 s ClockTime = 1 s
shogan50 is offline   Reply With Quote

Old   March 29, 2010, 13:41
Default
  #5
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,888
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
Hello,

did you check the skewness of the mesh?

Your pressure equation has a hard time to converge, and after 11 correctors the initial residual is still around 0.1, which means "not converged" ;-)

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   April 3, 2010, 13:42
Default skewness
  #6
New Member
 
Scott Hogan
Join Date: Mar 2010
Location: Washington State, USA
Posts: 6
Rep Power: 6
shogan50 is on a distinguished road
I'm not sure if this is the exact mesh (may have made some changes), but potential foam shows similar results. Here is the checkMesh output. Is there a good description of skewness somewhere? I'm using netgen 1d-2d-3d in Salome to mesh a step output solidworks model.

Code:
scotth@scotth-ubuntu:~/OpenFOAM/scotth-1.6/run/trumpetLES$ checkMesh
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.6                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 1.6-53b7f692aa41
Exec   : checkMesh
Date   : Apr 03 2010
Time   : 09:34:13
Host   : scotth-ubuntu
PID    : 10279
Case   : /home/scotth/OpenFOAM/scotth-1.6/run/trumpetLES
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           6907
    faces:            65092
    internal faces:   58628
    cells:            30930
    boundary patches: 5
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     0
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    30930
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch               Faces    Points   Surface topology                  
    wall                3415     1804     ok (non-closed singly connected)  
    inlet               624      350      ok (non-closed singly connected)  
    outlet              365      210      ok (non-closed singly connected)  
    exh                 53       38       ok (non-closed singly connected)  
    symmetry            2007     1072     ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.109217 0 0) (0.109217 0.109151 0.3048)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (5.4753e-19 -5.58808e-20 -4.65674e-19) OK.
    Max cell openness = 1.3833e-16 OK.
    Max aspect ratio = 4.69247 OK.
    Minumum face area = 1.40048e-05. Maximum face area = 0.000151843.  Face area magnitudes OK.
    Min volume = 2.60542e-08. Max volume = 5.24061e-07.  Total volume = 0.00433954.  Cell volumes OK.
    Mesh non-orthogonality Max: 51.3832 average: 16.7662
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.678317 OK.

Mesh OK.
shogan50 is offline   Reply With Quote

Old   April 3, 2010, 14:40
Default
  #7
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,888
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
FYI, similar problems on meshes with tets were described here

Overestimated temperature values

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   April 27, 2010, 12:41
Default
  #8
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 210
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Regarding my ivestigations at Overestimated temperature values:

I made calculations with foamCalc:

1. foamCalc div phi maximal error of order 10^-4
2. foamCalc div U maximal errors of order 10^2!!!

For U I set "div(U) Gauss linear;" in fvSchemes
Does that a reason for early described errors?
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   April 27, 2010, 12:56
Default
  #9
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,888
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
Hello,

you should look at div(phi) since phi is the conservative flux, and yes, that error is telling you the solution is not converging well.

Try reducing the tolerances on the pressure linear solver to 10^-15 or 10^-20 (it will slow down the solution a lot, usually with these strict requirements conjugate gradients methods are faster than GAMG in my experience). In addition, consider more non-orthogonal correctors (the number depends on how your specific convergence behaviour).

Best,
Alberto
__________________
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   April 27, 2010, 13:00
Default
  #10
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 210
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Thanks, I'm trying your guides right now...
Should div(U) also go to zero or is that something apart from flux conservation?
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   April 27, 2010, 13:04
Default
  #11
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,888
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
It should be a lot smaller than what you have. The difference between div(U) and div(phi) is due to the interpolation. The flux phi is the variable responsible of ensuring the conservation, while you can look at U as a secondary variable.

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   April 27, 2010, 13:07
Default
  #12
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 210
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Ok! Trying to treat it more, because I found that most errors for div(U) are at the same cells, which appear as an numerical heat sources in my case...
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   April 27, 2010, 13:49
Default
  #13
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,888
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
Lovely skewed cells... :-D
__________________
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   April 28, 2010, 15:28
Default
  #14
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 210
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Some one could please explain me following result?
I got this picture in icoFoam cavity case, tolerance for U and p is 1e-30, steady state achieved.

divU.jpg

Picture represents logarithmic plot of the interpolated velocities divergence.
Utility, which I used to achieve such result is also attached.

divCheck.tar.gz

You can compile utility with wmake command and use then:

divCheck -time YourTime


Waiting for your comments!
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   April 28, 2010, 15:49
Default
  #15
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,888
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
If you compute fvc::div(phi), which is what you have to check, since phi is the conservative flux, you obtain (cavity tutorial, time = 0.5, probably not at perfect steady state):

Divergens absolute max/min/avg values: 0.0004 0 4.729e-05

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.

Last edited by alberto; April 28, 2010 at 15:53. Reason: added note
alberto is offline   Reply With Quote

Old   April 29, 2010, 02:29
Default
  #16
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 210
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Alberto,

thank you for your comment! I also have checked it before with "foamCalc div phi" and calculated the flow for 5 seconds, sorry I didn't mention that.
But my question is - why flux, interpolated from U cell values, does not correspond to continuity equation

div U = 0?

Is this a PISO / SIMPLE feature? Or perhaps some other scheme apart from "linear" should be used for interpolate(U)? And do you mean, that if "div phi = 0" I should not worry about mass conservation?
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   April 29, 2010, 11:41
Default
  #17
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,888
Rep Power: 25
alberto will become famous soon enoughalberto will become famous soon enough
Hi,

in OpenFOAM you enforce the mass conservation by imposing

div(phi) = 0

where phi is the mass or volumetric flux (meaning the velocity interpolated on faces, dotted with the surface + terms due to co-located grids). This is what you actually do also on codes based on the staggered grid arrangement, but there you do not risk any confusion, because the velocity does not need to be interpolated, since you have it already on the face.

From the previous equation, you compute a pressure that ensures the the flux phi is conservative after the correction is applied, and then, you use phi in the linearization of the momentum equation and in the evaluation of the other convective terms (for example in scalar transport equations).

So, strictly speaking, when you check if the mass conservation is enforced correctly, you should check div(phi). See for example continuityErrs.H

I would say that, if you want to calculate the difference between div(phi), with phi coming from the pEqn, and the interpolated U, you should compute

div(fvc::interpolate(U) & mesh.Sf())

after the velocity has been corrected with the updated pressure field. At this point you can calculate the difference of the correct flux and the interpolated one.

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.

Last edited by alberto; April 29, 2010 at 11:43. Reason: Added note
alberto is offline   Reply With Quote

Old   April 29, 2010, 12:43
Default
  #18
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 210
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Thank you, Alberto! Now it is more clear for me.
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Reply

Tags
potentialfoam tets

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



All times are GMT -4. The time now is 07:20.