Convergence problems with simpleFoam on human airway
 

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 steady-state solution. However, the results are very bad while the case converges perfectly in Fluent.

First, what do I want to do:

• Obtain a laminar steady-state solution
• The model has three patches: inlet, outlet, airway
  • inlet: pressure inlet of 0 Pa (atmosphere)
  • outlet: pressure outlet of -20 Pa (under-pressure in the lungs)
  • airway: wall
• Unstructured volume mesh created in TGrid and converted to OpenFoam with the command: fluent3DMeshToFoam -scale 0.001 meshfile

Using checkMesh on the converted mesh gives me following output:

As one can see, the quality of the mesh is ok. However, due to the complex geometry, there are some non-orthogonal faces (> 60°). Therefore, I have used limited schemes in my fvSchemes file:

My fvSolution file is a rather standard one with nCellsInCoarsestLevel set to a bit more than sqrt(meshsize). The value for nNonOrthogonalCorrectors is set to 0 as this is recommended for a steady-state SIMPLE solution. When setting this to a higher value, the residuals for p even become worse.

As I have mentioned before, we want to simulate an under-pressure of 20 Pa at the outlet, while the inlet should stay 0 Pa (atmosphere). Therefore my p file looks like this:

And the corresponding U file like this:

Now, when running the case, I am expecting total pressure drop of about -14 Pa and a mass flow of about 0.00023 kg/s. The solution converges to these values until the 200th iteration, but after that it goes terribly wrong and the system diverges (see iteration 300 and further). One can see that the time step continuity errors become very high. In the attached logfiles.zip file, you can see the logfile (logfile_pvi_io.txt) and the residual plot (residuals_pvi_io.eps). The reported mass flows are in kg/s and pressures in Pa (everything is multiplied by rho).

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:

Using these boundary conditions, the case does not diverge any more and the values are rather good. However, the values of mass flow and total pressure drop still fluctuate and the residual plot is still not very good. Also, I do not think these boundary conditions are physically correct for my case. My first U file looks much more like the commonly used options for pressure inlet and pressure outlet. In the attached logfiles.zip file, you can see the logfile (logfile_freestream.txt) and the residual plot (residuals_freestream.eps).

Has anybody an idea to make my case converge?

Thanks in advance

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:

My time step continuity errors are just way to high to get a converged solution.

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?

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 CT-scan 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

Just a few additions to what said above. You can limit the gradient (meaning what computed with fvc::grad(..)) with

and you can limit also the gradient used in the linearUpwind scheme as said by Matej, with

replacing linearUpwindV with linearUpwind if instead if a vector field you have a scalar.

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,

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 patient-specific basis, so my geometries are segmented from CT-scans. 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 under-pressures. 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:

• I have decreased relTol on p to 0.001 and 0.01 on U
• I have set the tolerance for p to 10e-8, and for U to 10e-7
• My fvSchemes file was already as all you suggested so I did not change it a lot, I have just set all the gradSchemes to "cellLimited Gauss linear 1" instead of only grad(U). Using cellMDLimited instead of cellLimited did not change the results. Setting div(phi,U) to Gauss upwind also did not improve the stability.

First, I have tried to set a velocity inlet (U = flowRateInletVelocity and p = zeroGradient) and a pressure outlet as suggested by matejfor. With these boundary conditions, the solutions converges very well and all residuals < 1e-3! In the logfiles2.zip file, you can see the logfile (logfile_Uinlet.txt) and the residual plot (residuals_Uinlet.eps). Using a velocity inlet solves all my problems. However, I really want to use a pressure inlet as all my results in Fluent are calculated that way.

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

From what you said after, the problem is not in the mesh. However "automated meshing" doesn't mean "I can't control the quality" :D

You fix both the pressures, so the pressure drop is fixed. Do you mean it's profile is wrong?

Best,

I meant that I have already created an optimal meshing strategy with tetrahedrals:) The skewness is alright and the solutions are mesh independent in Fluent.

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.

That's great.

Hint: try with a simple 2D case with a rectangular channel to figure out how to fix the boundaries :D

Best,