
[Sponsors] 
October 19, 2010, 12:45 

#21  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
Quote:
Best, 

October 19, 2010, 12:52 

#22  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:


October 19, 2010, 12:54 

#23 
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 

October 21, 2010, 12:20 

#24  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:
as I said I've tried to follow your suggestions, but unfortunately the results were not satisfactory...In particular, I've tried the following new settings:  URF factors for k and epsilon reduced to 0.4  Tolerances of the solver reduced to 10^10 for p and 10^08 for the other quantities  Relative tolerance for p reduced to 10^03 (from 10^02)  divSchemes set to "Gauss linearUpwindV cellMDLimited Gauss linear" for the velocities and the same without "V" for k and epsilon Additionally, before switching to linearUpwind I let the simulation go on with all the divSchemes set to the basic upwind scheme for 250 iterations, hoping this would improve the convergence of the run...But, as I said before, the results were not good. In particular, after a few iterations the residuals of k fall suddenly to a very low value, so in practice for the rest of the run the kequation is no more solved. The other equations instead, keep to be solved and the relative residuals reduce slowly to lower values (with the exception of epsilon, whose residuals remain fairly constant). Apart from this, the solution after 1500 iterations is definetly not in accordance with what one should expect from the Ahmed body case: in particular, the tke distribution in the symmetry plane is completely unphysical (but this should be expected from a practically not resolved field), and the drag coefficient associated with the nose part of the body is heavily underestimated comparing with both experimental and other numerical data (including my previous allupwind run). For more clearness, I will attache the images concerning the convergence patterns, the tke distribution in the symmetry plane and the nose drag coefficients calculated in two different runs (the linearUpwind one and a previous allupwind one). Of course any help would be really appreciated... Regards V. 

October 21, 2010, 12:53 

#25 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hello,
the solution is not converged. You are doing a steady state simulation: residuals should be significantly lower. Set the convergence criterion at 10^10 and wait until that value is reached or the residuals become completely flat for a number of iterations. I would not start with upwind and then switch scheme. It is of no use. If the problem is set up correctly, it will converge also with the secondorder scheme. Since k is not being solved anymore, the underrelaxation of 0.4 is probably too strong (assuming the physical setup is correct). You might want to increase the URF again. Best, 

October 21, 2010, 13:24 

#26  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:
ok, I'll try with these new settings but let me add some more observations:  The same behaviour of tke's residuals can be viewed with the URFs set to 0.7: does it mean that even 0.7 is a too strong URF for k and epsilon?  Let's say that lowering further the tolerance I will come to some satisfactory solution: it seems to me (but maybe I'm wrong about it) that this is a quite restrictive condition... should I amplify many times the total runtime in any situation that would require a lower numerical diffusion? Or maybe, as you said, this is simply a typical behaviour of the steadystate solvers (in this last case, It seems like it would be much more convenient to run a transient simulation and wait until it will converge to a steadystate configuration)? However, thanks once again for your replies regards V. 

October 21, 2010, 14:39 

#27  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
Quote:
Running an unsteady case would require significantly longer in general, since obtaining those residual values is fast with steady solvers if the mesh is good and the setup is correct. Best, 

October 22, 2010, 09:21 

#28  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:
Apart from this, I lowered the tolerance of the solvers to 10^11 for p and to 10^10 for the other quantities and (as it could be seen from the residuals pattern atached) there are no substantial differences, because tke's residual still falls down below the imposed tolerance after a few iterations, and after that all the other residuals does not converge at all...I made a review of all the settings of the case, but honestly I cannot find a reason for which the kequation behaves like this (apart from the numerical scheme chosen for the convection term). I can post again my fvSchemes and fvSolution dictionary, maybe they could add some more information which probably I'm missing... Regards V. fvSchemes: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) Gauss linear; } divSchemes { default none; div(phi,U) Gauss linearUpwindV cellMDLimited Gauss linear 1; div(phi,k) Gauss linearUpwind cellMDLimited Gauss linear 1; div(phi,epsilon) Gauss linearUpwind cellMDLimited Gauss linear 1; div(R) Gauss linear; div((nuEff*dev(grad(U).T()))) Gauss linear; } laplacianSchemes { default none; laplacian(nuEff,U) Gauss linear corrected; laplacian((1A(U)),p) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; } // ************************************************** *********************** // fvSolution: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { p { solver PCG; preconditioner DIC; tolerance 1e11; relTol 0.001; } U { solver PBiCG; preconditioner DILU; tolerance 1e10; relTol 0.1; } k { solver PBiCG; preconditioner DILU; tolerance 1e10; relTol 0.1; } epsilon { solver PBiCG; preconditioner DILU; tolerance 1e10; relTol 0.1; } } SIMPLE { nNonOrthogonalCorrectors 3; } relaxationFactors { p 0.3; U 0.7; k 0.7; epsilon 0.7; } // ************************************************** *********************** // 

October 22, 2010, 14:51 

#29  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Quote:
Quote:
Best, 

October 23, 2010, 04:56 

#30 
Senior Member
Markus Rehm
Join Date: Mar 2009
Location: Erlangen (Germany)
Posts: 180
Rep Power: 8 
Hello Vesselin,
have you tried to set k and epsilon tolerance to real machine precision (1e16)? Does make a difference? Regards, Markus. 

October 25, 2010, 04:41 

#31  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:
Hi, Markus if you're talking about the absolute tolerances I didn't try such a setting (at most I have lowered the tolerances to 10^10, 10^11), but i don't think that it could make some difference, because the run starts to have some problems with the kequation after a few iterations and after this point the residuals stop to converge at all (you can check this behaviour in the image attached to my previous post), remaining at values of about 10^03/10^04. Regards V. 

October 25, 2010, 11:17 

#32 
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
However, my case is a 3D external aerodynamics case, with a symmetric domain. The domain itself is composed by inlet and outlet sections, the symmetry plane, a floor and top and side patches. I've tried the following settings:
setting 1: top and side treated as generic patches, with a zeroGradient condition for all the quantities involved in the calculation Of course the inlet and outlet are generic patches and the symmetry plane is symmetryPlane. Here there is the content of my 0/ folder for this setting. U: internalField uniform (40 0 0); boundaryField { inlet { type fixedValue; value uniform (40 0 0); } outlet { type zeroGradient; } body { type fixedValue; value uniform (0 0 0); } nose { type fixedValue; value uniform (0 0 0); } slant { type fixedValue; value uniform (0 0 0); } back { type fixedValue; value uniform (0 0 0); } floor { type fixedValue; value uniform (0 0 0); } side { type zeroGradient; } top { type zeroGradient; } symmetry { type symmetryPlane; } } p: internalField uniform 0; boundaryField { inlet { type zeroGradient; } outlet { type fixedValue; value uniform 0; } body { type zeroGradient; } nose { type zeroGradient; } slant { type zeroGradient; } back { type zeroGradient; } floor { type zeroGradient; } symmetry { type symmetryPlane; } side { type zeroGradient; } top { type zeroGradient; } } k: internalField uniform 0.0096; boundaryField { inlet { type fixedValue; value uniform 0.0096; } outlet { type zeroGradient; } body { type kqRWallFunction; value uniform 0.0096; } nose { type kqRWallFunction; value uniform 0.0096; } slant { type kqRWallFunction; value uniform 0.0096; } back { type kqRWallFunction; value uniform 0.0096; } floor { type kqRWallFunction; value uniform 0.0096; } side { type zeroGradient; } top { type zeroGradient; } symmetry { type symmetryPlane; } } epsilon: internalField uniform 0.155; boundaryField { inlet { type fixedValue; value uniform 0.00155; } outlet { type zeroGradient; } body { type epsilonWallFunction; value uniform 0.00155; } nose { type epsilonWallFunction; value uniform 0.00155; } slant { type epsilonWallFunction; value uniform 0.00155; } back { type epsilonWallFunction; value uniform 0.00155; } floor { type epsilonWallFunction; value uniform 0.00155; } side { type zeroGradient; } top { type zeroGradient; } symmetry { type symmetryPlane; } } nut: internalField uniform 0; boundaryField { inlet { type calculated; value uniform 0; } outlet { type calculated; value uniform 0; } body { type nutWallFunction; value uniform 0; } nose { type nutWallFunction; value uniform 0; } slant { type nutWallFunction; value uniform 0; } back { type nutWallFunction; value uniform 0; } floor { type nutWallFunction; value uniform 0; } side { type zeroGradient; value uniform 0; } top { type zeroGradient; value uniform 0; } symmetry { type symmetryPlane; } } k and epsilon values are calculated imposing a turbulence inlet intensity of 0.2% (reffered to the inlet bulk velocity). setting 2: is very close to setting 1, but the top and side patches are imposed to be symmetry planes. All the entries in the 0/ folder are changed consequently As you already know, both of the settigns does not work with higher order convection schemes, even with a mesh with quite good quality parameters (you can check all of them in my previous posts). Any additional ideas? V. 

October 27, 2010, 06:04 

#33 
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Does anybody think that changing the linear solver of some of the quantities (for instance, switching to GAMG for p) would help?
V. 

November 2, 2010, 09:34 

#34  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:
Thank you in advance Best Regards V. fvSchemes: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default steadyState; } gradSchemes { default leastSquares; grad(p) leastSquares; grad(U) leastSquares; } divSchemes { default none; div(phi,U) Gauss linearUpwindV faceMDLimited leastSquares 1; div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind; div(R) Gauss linear; div((nuEff*dev(grad(U).T()))) Gauss linear; } laplacianSchemes { default Gauss linear limited 0.5; // laplacian(nuEff,U) Gauss linear corrected; // laplacian((1A(U)),p) Gauss linear corrected; // laplacian(DkEff,k) Gauss linear corrected; // laplacian(DepsilonEff,epsilon) Gauss linear corrected; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default limited 0.5; } fluxRequired { default no; p ; } // ************************************************** *********************** // fvSolution: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { // p // { // solver PCG; // preconditioner DIC; // tolerance 1e11; // relTol 0.001; // } p { solver GAMG; tolerance 1e11; relTol 0.0001; smoother GaussSeidel; cacheAgglomeration true; nCellsInCoarsestLevel 1000; agglomerator faceAreaPair; mergeLevels 1; } U { solver PBiCG; preconditioner DILU; tolerance 1e10; relTol 0.1; } k { solver PBiCG; preconditioner DILU; tolerance 1e10; relTol 0.1; } epsilon { solver PBiCG; preconditioner DILU; tolerance 1e10; relTol 0.1; } } SIMPLE { nNonOrthogonalCorrectors 3; } relaxationFactors { p 0.3; U 0.7; k 0.5; epsilon 0.5; } // ************************************************** *********************** // checkMesh report: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 357845 faces: 3218365 internal faces: 3105296 cells: 1527609 boundary patches: 10 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 0 prisms: 213225 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 1314384 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 symmetry 22600 12584 ok (nonclosed singly connected) floor 14660 7609 ok (nonclosed singly connected) inlet 426 242 ok (nonclosed singly connected) outlet 432 245 ok (nonclosed singly connected) side 2228 1204 ok (nonclosed singly connected) top 1648 910 ok (nonclosed singly connected) body 54899 27768 ok (nonclosed singly connected) nose 8704 4461 ok (nonclosed singly connected) slant 3948 2058 ok (nonclosed singly connected) back 3524 1841 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (1.26 0 0.1945) (6.27 1.44 1.1945) Mesh (nonempty, nonwedge) directions (1 1 1) Mesh (nonempty) directions (1 1 1) Boundary openness (8.32191e20 1.52667e19 7.68048e20) OK. Max cell openness = 2.7637e16 OK. Max aspect ratio = 5.97752 OK. Minumum face area = 1.71072e06. Maximum face area = 0.0159908. Face area magnitudes OK. Min volume = 1.10751e09. Max volume = 0.000575106. Total volume = 10.7878. Cell volumes OK. Mesh nonorthogonality Max: 54.8744 average: 14.5203 Nonorthogonality check OK. Face pyramids OK. Max skewness = 1.98187 OK. Mesh OK. End 

November 2, 2010, 09:52 

#35  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:


November 2, 2010, 11:41 

#36 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hi,
I have used the setup I told you successfully in almost all my runs, so I would doublecheck the boundary conditions. Your checkMesh does not have anything wrong. Out of curiosity, what happens to the residuals if you let the case run more than 2000 iterations? Do they keep increasing or do they oscillate? Best, 

November 2, 2010, 12:07 

#37  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:
Symmetrycal domain composed by a half Ahmedd model and the surrounding wind tunnel The symmetry plane is (obviously) set as symmetryPlane bc The inlet is set to fixed values for U, k and epsilon and to zeroGradient for p The outlet is set to zeroGradient for all the quantities except p, which is fixed to 0. The floor and the body are set as walls, so fixedValue (0 0 0) for U, zeroGradient for p and wallFunctions for k, epsilon and nut (I'm using the highRe realizable kepsilon model). For the side and top patches I decided to use the zeroGradient condition for all the quantities. I also tried to switch to symmetryPlane, without (seemingly) any benefits. To my knowledge, and having a look around for the literature on the Ahmed body (incompressible) simulations, such a BC setting seems to my quite standard. About the second issue, the same run is now going on further: I could tell you more about it maybe this evening or tomorrow morning Best Regards V. 

November 3, 2010, 06:29 

#38  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:
V. Last edited by vkrastev; November 3, 2010 at 12:05. 

November 4, 2010, 04:46 

#39 
Senior Member
Matthias Voß
Join Date: Mar 2009
Location: Berlin, Germany
Posts: 446
Rep Power: 11 
hi,
since i´m not very familiar with all of the impacts of the schemes and everything this is maybe smth. totally wrong: if i look at the flow and keep in mind that my residuals aren't converging tight enough... isn't that supposing a transient effect (back flow etc.) in my flow field leading to the "poor" convergence? > check for the high residuals within the domain.. if they are at the back of the body this could be an hint on transient effects. neewbie 

November 4, 2010, 06:21 

#40  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 361
Rep Power: 10 
Quote:
about the backflow I don't see any trace of it during the runs (the outlet patch is at 5 body lenghts downstream of the body itself: at this distance the flow is again quite unidirectional, so having a backflow there would be a really uncommon situation...). About to check where the residuals are higher in the domain, can you give me a hint on how it could be done? Thank you V. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
4th order schemes in channelOodles  maka  OpenFOAM Bugs  9  January 19, 2009 12:58 
Help for fourth order accurate convective schemes  Z.C.Wang  Main CFD Forum  0  January 15, 2009 07:53 
2nd order conservative schemes  taw  Main CFD Forum  1  September 16, 2008 07:05 
CFL condition for higher order schemes  Shyam  Main CFD Forum  2  February 14, 2008 15:24 
High order compact finite difference schemes  Mikhail  Main CFD Forum  6  August 5, 2003 10:36 