
[Sponsors] 
May 21, 2010, 05:48 
Convergence problems with simpleFoam on human airway

#1 
Member
Cedric Van Holsbeke
Join Date: Dec 2009
Location: Belgium
Posts: 81
Rep Power: 10 
I'm trying to simulate the flow trough a part of the human airways. I'm using simpleFoam (with an added volScalarField ptot in createFields.H in order to write out the total pressure) to obtain a steadystate solution. However, the results are very bad while the case converges perfectly in Fluent.
First, what do I want to do:
Code:
Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 122864 faces: 1203339 internal faces: 1098113 cells: 575363 boundary patches: 3 point zones: 0 face zones: 1 cell zones: 1 Overall number of cells of each type: hexahedra: 0 prisms: 0 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 575363 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 airway 96829 48666 ok (nonclosed singly connected) inlet 4629 2466 ok (nonclosed singly connected) outlet 3768 1990 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0.106216 0.0977228 0.103362) (0.144808 0.119839 0.0257335) Mesh (nonempty, nonwedge) directions (1 1 1) Mesh (nonempty) directions (1 1 1) Boundary openness (1.38743e17 5.89666e17 2.40551e16) OK. Max cell openness = 2.57245e16 OK. Max aspect ratio = 12.5695 OK. Minumum face area = 6.96372e10. Maximum face area = 4.44277e07. Face area magnitudes OK. Min volume = 1.84058e14. Max volume = 9.0144e11. Total volume = 1.02185e05. Cell volumes OK. Mesh nonorthogonality Max: 75.7532 average: 19.7554 *Number of severely nonorthogonal faces: 8. Nonorthogonality check OK. <<Writing 8 nonorthogonal faces to set nonOrthoFaces Face pyramids OK. Max skewness = 1.25992 OK. Mesh OK. End Code:
ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) cellLimited Gauss linear 1; } divSchemes { default none; div(phi,U) Gauss linearUpwind cellLimited Gauss linear 1; div((nuEff*dev(grad(U).T()))) Gauss linear; } laplacianSchemes { default Gauss linear limited 0.333; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default limited 0.333; } fluxRequired { default no; p; } Code:
solvers { p { solver GAMG; tolerance 1e05; relTol 0.01; smoother GaussSeidel; nCellsInCoarsestLevel 900; agglomerator faceAreaPair; mergeLevels 1; cacheAgglomeration true; }; U { solver smoothSolver; tolerance 1e05; relTol 0.1; smoother GaussSeidel; }; } SIMPLE { nNonOrthogonalCorrectors 0; } relaxationFactors { p 0.3; U 0.7; } Code:
dimensions [0 2 2 0 0 0 0]; internalField uniform 0; boundaryField { inlet { type fixedValue; value uniform 0; } outlet { type fixedValue; // Pa = 1 kg/(m·s^2). Dimensions for incompressible solver: m^2/s^2. // 20 Pa = 16.3265306122449 m^2/s^2 (rho = 1.225 kg/s) value uniform 16.3265306122449; } airway { type zeroGradient; } } Code:
dimensions [0 1 1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { inlet { type pressureInletVelocity; value uniform (0 0 0); } outlet { type inletOutlet; inletValue uniform (0 0 0); value uniform (0 0 0); } airway { type fixedValue; value uniform (0 0 0); } } As all my files look good, I thought that it had something to do with the boundary conditions. For that reason, I altered the U file: Code:
dimensions [0 1 1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { inlet { type freestream; freestreamValue uniform (0 0 0); } outlet { type freestream; freestreamValue uniform (0 0 0); } airway { type fixedValue; value uniform (0 0 0); } } Has anybody an idea to make my case converge? Thanks in advance Last edited by CedricVH; June 2, 2010 at 08:32. 

May 21, 2010, 07:44 

#2 
Member
Cedric Van Holsbeke
Join Date: Dec 2009
Location: Belgium
Posts: 81
Rep Power: 10 
In order to reduce mesh scaling errors and rounding errors, I have set my writeFormat to binary. I have also tried to run potentialFoam first, but this makes things even worse as the initial conditions are assumed wrong (even with high (50) nNonOrthogonalCorrectors). As you can see it then goes wrong from the first time step:
Code:
Time = 1 smoothSolver: Solving for Ux, Initial residual = 0.33810348994305, Final residual = 0.021350422410913, No Iterations 3 smoothSolver: Solving for Uy, Initial residual = 0.36170510392225, Final residual = 0.021805047282949, No Iterations 3 smoothSolver: Solving for Uz, Initial residual = 0.15638801269015, Final residual = 0.010144109442471, No Iterations 3 GAMG: Solving for p, Initial residual = 1, Final residual = 0.0088716566191601, No Iterations 4 time step continuity errors : sum local = 222.61631013276, global = 10.24764395806, cumulative = 10.24764395806 ExecutionTime = 0.97 s ClockTime = 1 s MassFlows: inlet = 0.023941635447993 outlet = 0.024069911643539 Averages of ptot : inlet = 3618.8845244285 outlet = 6304.3978782402 

May 21, 2010, 08:48 

#3 
Member
Johan Spång
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 35
Rep Power: 10 
Try lowering relTol on p to 0.001 and 0.01 on U. Did you try zeroGradient on U inlet?
Can you post the case? 

May 22, 2010, 03:54 

#4 
Senior Member
matej forman
Join Date: Mar 2009
Location: Brno, Czech Republic
Posts: 133
Rep Power: 10 
Hi,
have you tried to run some flow with a given inlet velocity? just to see where are the problems. I'm doing human airways simulations myself but except real geometry from CTscan I'm pretty sure the mesh could always be of better quality. I'd run the case with given velocity to see where the mesh is wrong, and then I'd try two things: A/ U at inlet set to simple fixedGradient and then I'd try more limiting as: cellLimitedGrad  The scalar limiter based on limiting the extrapolated face values between the maximum and minumum cell and cell neighbour values and is applied to all components of the gradient. $FOAM/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C in fvSchemes: div(phi,U) Gauss linearUpwind cellLimited Gauss linear 1; cellMDLimitedGrad  the same as cellLimitedGrad with explicit limiter btw. I'd like to know how did you find out 20Pa pressure drop in airways. Do you have papers or reports on your work? good luck matej 

May 23, 2010, 02:58 

#5  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
Quote:
Code:
gradSchemes { default cellLimited Gauss linear 1; } Code:
div(phi,U) Gauss linearUpwindV cellLimited Gauss linear 1; You can additionally use cellMDLimited instead than cellLimited, if you want the limiter to be applied independently in the three directions. Additionally, lower the tolerance required for the pressure to at least one order less than the one used for the velocity. I'd set the pressure to 10^8, and U to 10^7, with a convergence criteria for the the steady solution ~ 10^6 (or the minimum you can reach, when residuals stop changing). Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

June 2, 2010, 08:27 

#6 
Member
Cedric Van Holsbeke
Join Date: Dec 2009
Location: Belgium
Posts: 81
Rep Power: 10 
Thanks for all the replies! I'm sorry that I did not answer faster, but I had other projects going on. I'm doing studies on patientspecific basis, so my geometries are segmented from CTscans. As I have a lot of patients to deal with, the meshing procedure has to be almost fully automatic (tri/tet in TGrid 5.0.6) so my mesh quality can not improve.
The 20Pa is just arbitrary in order to obtain a laminar solution, it has nothing to do with the real underpressures. I have just started my PhD in October and do not have any A1 publications. However, I have just send in my first paper and I will post it here when it is published. I have changed several things according to your suggestions:
When using my original boundary conditions for the inlet (U = pressureInletVelocity and p = fixedValue 0), the system still diverges and gives very unphysical values. In the logfiles2.zip file, you can see the logfile (logfile_pvi_io2.txt) and the residual plot (residuals_pvi_io2.eps). Like josp suggested, I have changed my boundary conditions for my pressure inlet to (U = zeroGradient and p = fixedValue 0). I have always learned that zeroGradient for U is not recommended as the system can become instable due to possible inward flow (and that therefore pressureInletVelocity was created). The results seem better and the case does not diverge as fast. However, at 1000 iterations, the mass flow rate is rather good, but the total pressure drop is completely wrong. Also the residuals start to increase. In the logfiles2.zip file, you can see the logfile (logfile_Uzerogradient.txt) and the residual plot (residuals_Uzerogradient.eps). Using freestream for U and fixedValue for p (for both inlet and outlet) still gives the best results if one wants to have a pressure inlet and pressure outlet. Using your suggestions, the residual plot looks a bit better than before and also the values for mass flow and pressure drop are physical and rather steady. It seems that these boundary condition are the best for my cases although they are not so obvious. In the logfiles2.zip file, you can see the logfile (logfile_freestream2.txt) and the residual plot (residuals_freestream2.eps). If somebody has another suggestion to improve my results, just bring it on! The logfiles2.zip file can be found at: http://www.2shared.com/file/jnlWk1Us/logfiles2.html Last edited by CedricVH; June 2, 2010 at 09:06. 

June 2, 2010, 11:58 

#7  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
Quote:
Quote:
Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

June 3, 2010, 03:28 

#8  
Member
Cedric Van Holsbeke
Join Date: Dec 2009
Location: Belgium
Posts: 81
Rep Power: 10 
Quote:
I fix the static pressures at inlet and outlet, however, I am interested in the total pressure at inlet and outlet in order to measure the pressure drop. These values are completely wrong. 

June 3, 2010, 09:05 

#9  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
Quote:
Quote:
Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
SimpleFoam convergence problems  brahim  OpenFOAM Running, Solving & CFD  20  June 9, 2015 09:09 
Getting faster convergence in simpleFoam  basneb  OpenFOAM  8  February 9, 2010 05:20 
Convergence of CFX field in FSI analysis  nasdak  CFX  2  June 29, 2009 01:17 
Convergence problems  Simone  Siemens  5  June 29, 2005 10:48 
Convergence problems in CFX5  Soren  CFX  18  March 23, 2002 20:32 