
[Sponsors] 
March 26, 2010, 11:22 
potentialFoam and Tets

#1 
New Member
Scott Hogan
Join Date: Mar 2010
Location: Washington State, USA
Posts: 6
Rep Power: 7 
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); } 

March 26, 2010, 15:14 
more info

#2 
New Member
Scott Hogan
Join Date: Mar 2010
Location: Washington State, USA
Posts: 6
Rep Power: 7 
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?
Last edited by shogan50; March 26, 2010 at 15:15. Reason: More info 

March 26, 2010, 19:16 

#3 
Senior Member
Sebastian Gatzka
Join Date: Mar 2009
Location: Frankfurt, Germany
Posts: 729
Rep Power: 11 
I made bad experience with tetmeshes, but they are not uncommon.
If nothing is wrong with your boundary condition itself you could try to increase the number of nonorthogonalcorrectors. 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 tetelements. 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!" 

March 27, 2010, 23:58 
Interesting

#4 
New Member
Scott Hogan
Join Date: Mar 2010
Location: Washington State, USA
Posts: 6
Rep Power: 7 
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.82863e05, No Iterations 51 continuity error = 0.0420761 Interpolated U error = 1.70791e07 ExecutionTime = 1.62 s ClockTime = 1 s 

March 29, 2010, 12:41 

#5 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
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, 

April 3, 2010, 12:42 
skewness

#6 
New Member
Scott Hogan
Join Date: Mar 2010
Location: Washington State, USA
Posts: 6
Rep Power: 7 
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 1d2d3d in Salome to mesh a step output solidworks model.
Code:
scotth@scotthubuntu:~/OpenFOAM/scotth1.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.653b7f692aa41 Exec : checkMesh Date : Apr 03 2010 Time : 09:34:13 Host : scotthubuntu PID : 10279 Case : /home/scotth/OpenFOAM/scotth1.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 (nonclosed singly connected) inlet 624 350 ok (nonclosed singly connected) outlet 365 210 ok (nonclosed singly connected) exh 53 38 ok (nonclosed singly connected) symmetry 2007 1072 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0.109217 0 0) (0.109217 0.109151 0.3048) Mesh (nonempty, nonwedge) directions (1 1 1) Mesh (nonempty) directions (1 1 1) Boundary openness (5.4753e19 5.58808e20 4.65674e19) OK. Max cell openness = 1.3833e16 OK. Max aspect ratio = 4.69247 OK. Minumum face area = 1.40048e05. Maximum face area = 0.000151843. Face area magnitudes OK. Min volume = 2.60542e08. Max volume = 5.24061e07. Total volume = 0.00433954. Cell volumes OK. Mesh nonorthogonality Max: 51.3832 average: 16.7662 Nonorthogonality check OK. Face pyramids OK. Max skewness = 0.678317 OK. Mesh OK. 

April 3, 2010, 13:40 

#7 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
FYI, similar problems on meshes with tets were described here
Overestimated temperature values Best, 

April 27, 2010, 11:41 

#8 
Senior Member

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 FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

April 27, 2010, 11:56 

#9 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
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 nonorthogonal correctors (the number depends on how your specific convergence behaviour). Best, Alberto 

April 27, 2010, 12:00 

#10 
Senior Member

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 FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

April 27, 2010, 12:04 

#11 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
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, 

April 27, 2010, 12:07 

#12 
Senior Member

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 FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

April 27, 2010, 12:49 

#13 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
Lovely skewed cells... :D


April 28, 2010, 14:28 

#14 
Senior Member

Some one could please explain me following result?
I got this picture in icoFoam cavity case, tolerance for U and p is 1e30, 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 FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

April 28, 2010, 14:49 

#15 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
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.729e05 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 opensource implementation of quadraturebased moment methods Last edited by alberto; April 28, 2010 at 14:53. Reason: added note 

April 29, 2010, 01:29 

#16 
Senior Member

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 FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

April 29, 2010, 10:41 

#17 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,895
Rep Power: 26 
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 colocated 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 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 opensource implementation of quadraturebased moment methods Last edited by alberto; April 29, 2010 at 10:43. Reason: Added note 

April 29, 2010, 11:43 

#18 
Senior Member

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 FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at 

Tags 
potentialfoam tets 
Thread Tools  
Display Modes  

