
[Sponsors] 
September 23, 2016, 09:01 
icoFoam Courant number growth problem

#1 
Member
Join Date: Jun 2016
Posts: 31
Rep Power: 9 
Hello,
I'm new to OpenFOAM 4.0 and only beginning to explore the solver. I want to run a laminar pipe flow case with either icoFoam or pisoFoam, but I've run into some issues. The pipe is 0.1m in diameter and 0.5m in length in the positive ydirection, the fluid is an oil with eta=0.1 Pa*s and rho=900kg/m³. The smallest tetraelement is about 1mm and with an expected u_max of 0.2m/s (inflow u_mean=0.1m/s), the smallest time step shouldn't exceed 0.005s for a Conumber of 1. I know that this issue is a rather common one while getting to know the solver and CFD in general, but I can't seem to find the root of it in my specific case. I've attached the necessary files. I tried a wide range of time steps, but the Conumber just keeps growing, so I suspect the fvsettings (taken from the icoFoam cavity tutorial) or rather my model setup to be the cause. I also tried to get it to run with pisoFoam, but I didn't know what to write in the turbulenceProperties file after setting the turbulenceModel to laminar. The boundary conditions for inflow and outflow patches don't seem to be overconstrained, so I don't know what to change there apart from setting the fixed inflow velocity to a mean mass flux, but the keyword for the mass flux doesn't seem to be described in the documentation, although I found some mentionings in older versions of OF via Google but not for recent ones. Also, how do I know which options to set in fvSchemes ? Especially which types of keywords I need in ddt, d2dt2, grad, div, laplacianSchemes  for examples, when do I need div(phi,U), when do I need more settings and how to choose them? I'd appreciate any help. 

September 23, 2016, 13:01 

#2 
Senior Member
Mahdi Hosseinali
Join Date: Apr 2009
Location: NB, Canada
Posts: 273
Rep Power: 18 
Hi,
1. I couldn't find any attached file. 2. If your turbulence model is laminar then you don't need any options there. 3. When using standard solvers which have been used extensively, usually the reason for divergence is the boundary conditions, can you share them with us? 4. If you are trying to use periodic boundary condition with pressure mass flux, try checking $FOAM_TUTORIALS/incompressible/pimpleFoam/LES/channel395 5. For those options I'd suggest to read user manual which gives you a good understanding what every choice will do. To make sure if you need it or not, you can comment it and if OF complains about not finding it you know it's needed, also you can write banana for the options, OF tells you he doesn't know banana, but here's your available keywords. 

September 23, 2016, 13:43 

#3 
Member
Join Date: Jun 2016
Posts: 31
Rep Power: 9 
Oh sorry, I'm new to this forum. I uploaded the files as .txt before posting though, I don't know what went wrong
Good to know, there was an error message complaining about the missing file. So I recon it must exist, but can be empty? See the attachments, I hope you can see them now No periodicities, but I might have a look nonetheless The keyword trick is a good one, I'll remember that 

September 27, 2016, 16:51 

#4 
Member
Join Date: Jun 2016
Posts: 31
Rep Power: 9 
Alright, I revised my BCs and found that the inflow direction was off. My workstation ran into some trouble and I had lost track of my models. The inflow direction is in negative y direction, but it somehow did not change anything. I got pisoFoam to run, but the Courantnumber just keeps growing, even with adjustTimeStep activated! And no matter the step size using icoFoam, it's the same. I checked a few setup videos like this one: https://www.youtube.com/watch?v=zDFk0GCZG68 and didn't find any obvious differences concerning the BC setup. Here: http://www.cfdonline.com/Forums/ope...flowpipe.html is some more input, but no improvement either. I noticed that Hypermesh did not create physicalType keywords in the boundary file, but adding them did'nt change anything. I could really use some advice as I don't know what to look for anymore.
p dimensions [0 2 2 0 0 0 0]; internalField uniform 0; boundaryField { inflow { type zeroGradient; } wall { type zeroGradient; } outflow { type fixedValue; value uniform 0; } } U dimensions [0 1 1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { inflow { type fixedValue; value uniform (0 0.1 0); } wall { type fixedValue; value uniform (0 0 0); } outflow { type zeroGradient; } } transportProperties transportModel Newtonian; nu [0 2 1 0 0 0 0] 0.1; controlDict application pisoFoam; startFrom startTime; startTime 0.000; stopAt endTime; endTime 0.1; deltaT 0.00025; writeControl timeStep; writeInterval 80; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression uncompressed; timeFormat general; timePrecision 6; runTimeModifiable yes; adjustTimeStep yes; maxCo 0.8; fvSchemes ddtSchemes { default Euler; } gradSchemes { default Gauss linear; grad(p) Gauss linear; } divSchemes { default none; div(phi,U) Gauss linear; div((nuEff*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear orthogonal; } interpolationSchemes { default linear; } snGradSchemes { default orthogonal; } fvSolution solvers { p { solver PCG; preconditioner DIC; tolerance 1e4; relTol 0; } pFinal { $p; relTol 0; } U { solver PBiCG; preconditioner DILU; tolerance 1e5; relTol 0; } } PISO { nCorrectors 2; nNonOrthogonalCorrectors 0; } turbulenceProperties simulationType laminar; 

September 28, 2016, 08:27 

#5 
Member
Join Date: Jun 2016
Posts: 31
Rep Power: 9 
I succeeded by taking the pipe flow tutorial from FOAMHouse and replacing the mesh with my own one (see http://thefoamhouse5.webnode.es/pr...idefoamcases/). The mesh consists now of 25mm tetras without boundary layer, so it's very rough. deltaT is set to 0.0005s and it runs smoothly, maxCo is 0.021 at the end of the simulation. However, if I set deltaT to 0.005s, Co explodes instantly. How come!? The Courant number should be around 10 times the previous Co, so what is happening there?


September 28, 2016, 12:41 

#6 
Senior Member
Mahdi Hosseinali
Join Date: Apr 2009
Location: NB, Canada
Posts: 273
Rep Power: 18 
Hi,
The best way to share your case is to remove all the time directories and constant/polyMesh if you have blockMeshDict so others can download it and run it. I'm glad you got somewhere, try changing your gradSchemes from "linear" to "limitedLinear 1" 

September 28, 2016, 14:17 

#7 
Member
Join Date: Jun 2016
Posts: 31
Rep Power: 9 
Hi,
I use Hypermesh, so no blockMeshDict available. I got some advice in the meantime, my choices in fvSchemes and fvSolution weren't all that great in general. I was just so focused on the BCs that I did not look there more than I thought to be necessary. 

June 20, 2017, 19:06 
fvSchemes and fvSolution

#8  
Member
Sugajen
Join Date: Jan 2012
Location: Tempe, USA
Posts: 52
Rep Power: 14 
Quote:
I am able to run the same flow with 0 bar at outlet. But when I include the pressure of 60 bar (which is 6010 sq.m/sq.s in OF units) the Co number increases. Also, I checked my case for a cuboid domain instead of the one with curved edges as seen below. For the cuboid domain with 60 bar, OF solves it perfectly without Co increase. So either it is the mesh or solver. The checkMesh and checkMesh allGeometry allTopology yielded no error. So I am pretty sure it is the fvScheme or fvSolution or both. All ideas and suggestions are most welcome! Do you mind sharing the fvSchemes and fvSolution that fixed your issue ? It would really be of great help. Thanks in advance! Currently I am using the following fvSchemes and fvSolution fvSchemes Code:
ddtSchemes { default Euler; // trusted answer } gradSchemes { default leastSquares; } divSchemes { default Gauss upwind; } laplacianSchemes { default Gauss linear uncorrected; } interpolationSchemes { default linear; } snGradSchemes { default orthogonal; } fvSolution Code:
solvers { p { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0.05; } pFinal { $p; relTol 0; } U { solver smoothSolver; smoother symGaussSeidel; tolerance 1e05; relTol 0; } } PISO { nCorrectors 1; nNonOrthogonalCorrectors 2; pRefCell 0; pRefValue 0; } 

June 21, 2017, 00:20 

#9 
Senior Member
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 744
Rep Power: 15 
High pressure gradient lead to high local velocity. It is quite normal that the max Co number increases then. I recommend refining the mesh in the areas with the pressure drop or increase.
__________________
Uwe Pilz  Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950) 

June 22, 2017, 14:28 
Solved!

#10  
Member
Sugajen
Join Date: Jan 2012
Location: Tempe, USA
Posts: 52
Rep Power: 14 
Quote:
Code:
ddtSchemes { default Euler; } gradSchemes { default cellLimited Gauss linear 1; } divSchemes { default Gauss upwind; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } 

June 23, 2017, 00:21 

#11 
Senior Member
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 744
Rep Power: 15 
> As you suggested I refined the mesh but it didn't make much of a difference
If you want to control the Co number you need to refine the mesh only in the areas of fast changes and high velocity (areas where the high Co occurs). This way you minimize the required wall time for a simulation. Or, in other words , you may use a coarser mesh in areas with low CO number.
__________________
Uwe Pilz  Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950) 

June 27, 2017, 07:19 

#12 
Member
Join Date: Jun 2016
Posts: 31
Rep Power: 9 
Hi Sugajen, I'm using the settings set by HelyxOS (I tried it once, unfortunately no longer supported by newer versions of OF). They're working quite well and I've used them ever since. However, I don't know if there are possibilities to improve the settings manually since I didn't change much except the usage of GAMG:
fvSchemes: Code:
/** C++ **\  o    o o  HELYXOS   o O o  Version: v2.3.1   o o  Web: http://www.engys.com   o   \**/ FoamFile { version 2.0; format ascii; class dictionary; location system; object fvSchemes; } ddtSchemes { default Euler; } gradSchemes { default Gauss linear; grad(nuTilda) cellLimited Gauss linear 1; grad(k) cellLimited Gauss linear 1; grad(kl) cellLimited Gauss linear 1; grad(omega) cellLimited Gauss linear 1; grad(epsilon) cellLimited Gauss linear 1; grad(q) cellLimited Gauss linear 1; grad(zeta) cellLimited Gauss linear 1; grad(v2) cellLimited Gauss linear 1; grad(f) cellLimited Gauss linear 1; grad(sqrt(kt)) cellLimited Gauss linear 1; grad(kt) cellLimited Gauss linear 1; grad(sqrt(kl)) cellLimited Gauss linear 1; } divSchemes { //default Gauss linear; div(phi,U) bounded Gauss linearUpwindV grad(U); div(phi,U) Gauss linear grad(U); div(phi,k) bounded Gauss linearUpwind grad(k); div(phi,epsilon) bounded Gauss linearUpwind grad(epsilon); div(phi,zeta) bounded Gauss linearUpwind grad(zeta); div(phi,q) bounded Gauss linearUpwind grad(q); div(phi,omega) bounded Gauss linearUpwind grad(omega); div(phi,nuTilda) bounded Gauss linearUpwind grad(nuTilda); div(phi,T) bounded Gauss limitedLinear 1; div(phi,kl) Gauss limitedLinear 1; div(phi,kt) Gauss limitedLinear 1; div(phi,R) Gauss upwind; div(R) Gauss linear; div((nuEff*dev(grad(U).T()))) Gauss linear; div((nuEff*dev2(T(grad(U))))) Gauss linear; div(phi,v2) bounded Gauss linearUpwind grad(v2); div(phi,f) bounded Gauss linearUpwind grad(f); } interpolationSchemes { default linear; interpolate(HbyA) linear; } laplacianSchemes { default Gauss linear limited 0.333; } snGradSchemes { default limited 0.333; } fluxRequired { default no; p ; } Code:
/** C++ **\  o    o o  HELYXOS   o O o  Version: v2.3.1   o o  Web: http://www.engys.com   o   \**/ FoamFile { version 2.0; format ascii; class dictionary; location system; object fvSolution; } PISO { nCorrectors 1; nOuterCorrectors 1; nNonOrthogonalCorrectors 1; nOuterCorrectors 1; } solvers { p { solver GAMG; agglomerator faceAreaPair; mergeLevels 1; cacheAgglomeration true; nCellsInCoarsestLevel 200; tolerance 1e6; relTol 0.001; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; minIter 1; } U { solver GAMG; agglomerator faceAreaPair; mergeLevels 1; cacheAgglomeration true; nCellsInCoarsestLevel 200; tolerance 1e6; relTol 0.001; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; minIter 1; } UFinal { solver GAMG; agglomerator faceAreaPair; mergeLevels 1; cacheAgglomeration true; nCellsInCoarsestLevel 200; tolerance 1e6; relTol 0.001; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; minIter 1; } pFinal { solver GAMG; smoother GaussSeidel; agglomerator faceAreaPair; mergeLevels 1; cacheAgglomeration true; nCellsInCoarsestLevel 200; tolerance 1e6; relTol 0.001; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; minIter 1; } p_rgh { solver GAMG; agglomerator faceAreaPair; mergeLevels 1; cacheAgglomeration true; nCellsInCoarsestLevel 200; tolerance 1e10; relTol 0.01; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; minIter 1; } f { solver GAMG; agglomerator faceAreaPair; mergeLevels 1; cacheAgglomeration true; nCellsInCoarsestLevel 200; tolerance 1e10; relTol 0.01; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; minIter 1; } p_rghFinal { solver GAMG; smoother GaussSeidel; agglomerator faceAreaPair; mergeLevels 1; cacheAgglomeration true; nCellsInCoarsestLevel 200; tolerance 1e6; relTol 0.0001; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; minIter 1; } fFinal { solver GAMG; smoother GaussSeidel; agglomerator faceAreaPair; mergeLevels 1; cacheAgglomeration true; nCellsInCoarsestLevel 200; tolerance 1e6; relTol 0.0001; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; minIter 1; } k { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } kl { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } kt { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } q { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } zeta { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } epsilon { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } R { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } nuTilda { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } omega { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } h { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } T { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } v2 { solver smoothSolver; smoother GaussSeidel; tolerance 1e6; relTol 0; minIter 1; } rho { solver PCG; preconditioner DIC; tolerance 0; relTol 0; minIter 1; } rhoFinal { solver PCG; preconditioner DIC; tolerance 0; relTol 0; minIter 1; } e { solver PBiCG; preconditioner DILU; tolerance 1e06; relTol 0; minIter 1; } } relaxationFactors { fields { p 1; p_rgh 1; rho 1; } equations { epsilon 1; f 1; h 1; k 1; kl 1; kt 1; nuTilda 1; omega 1; q 1; R 1; T 1; U 0.3; v2 1; zeta 1; "alpha.*" 1; epsilonFinal 1; fFinal 1; hFinal 1; kFinal 1; klFinal 1; ktFinal 1; nuTildaFinal 1; omegaFinal 1; qFinal 1; RFinal 1; TFinal 1; UFinal 1; v2Final 1; zetaFinal 1; "alpha.*Final" 1; } } 

May 24, 2023, 14:00 
Courant Number and Residual Exponential Growth icoFoam

#13 
New Member
Andre Nel
Join Date: Sep 2022
Location: Skippack PA
Posts: 2
Rep Power: 0 
I had the same problem. I finally figured out that it might be the numerical scheme might be the problem.
So I outlined the problem to chatGPT4 and asked it to write a fvScheme. It did and I copied the text into the fvScheme. Worked. I am using icoFoam. Below is the fvShemes entry that the AI provided: /** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: v2012   \\ / A nd  Website: www.openfoam.com   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } ddtSchemes { default Euler; } gradSchemes { default Gauss linear; limited cellLimited Gauss linear 1; } divSchemes { default none; div(phi,U) Gauss linearUpwind grad(U); div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind; div(phi,R) Gauss upwind; div(R) Gauss linear; div(phi,nuTilda) Gauss upwind; div(phi,omega) Gauss linearUpwind grad(omega); div((nuEff*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; } // ************************************************** *********************** // 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
[Other] Can't Shake Erros: patch type 'patch' not constraint type 'empty'  BrendaEM  OpenFOAM Meshing & Mesh Conversion  12  April 3, 2022 18:32 
SigFpe when running ANY application in parallel  Pj.  OpenFOAM Running, Solving & CFD  3  April 23, 2015 14:53 
decomposePar pointfield  flying  OpenFOAM Running, Solving & CFD  28  December 30, 2013 15:05 
same geometry,structured and unstructured mesh,different behaviour.  sharonyue  OpenFOAM Running, Solving & CFD  13  January 2, 2013 22:40 
IcoFoam parallel woes  msrinath80  OpenFOAM Running, Solving & CFD  9  July 22, 2007 02:58 