[CedricVH]

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:

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 (non-closed singly connected)  
    inlet               4629     2466     ok (non-closed singly connected)  
    outlet              3768     1990     ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (0.106216 0.0977228 -0.103362) (0.144808 0.119839 -0.0257335)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-1.38743e-17 5.89666e-17 2.40551e-16) OK.
    Max cell openness = 2.57245e-16 OK.
    Max aspect ratio = 12.5695 OK.
    Minumum face area = 6.96372e-10. Maximum face area = 4.44277e-07.  Face area magnitudes OK.
    Min volume = 1.84058e-14. Max volume = 9.0144e-11.  Total volume = 1.02185e-05.  Cell volumes OK.
    Mesh non-orthogonality Max: 75.7532 average: 19.7554
   *Number of severely non-orthogonal faces: 8.
    Non-orthogonality check OK.
  <<Writing 8 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
    Max skewness = 1.25992 OK.

Mesh OK.

End

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:

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;
}

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.

solvers
{
    p
    {
        solver           GAMG;
        tolerance        1e-05;
        relTol           0.01;
        smoother         GaussSeidel;
        nCellsInCoarsestLevel 900;
        agglomerator     faceAreaPair;
        mergeLevels      1;
        cacheAgglomeration true;
    };

    U
    {
        solver           smoothSolver;
        tolerance        1e-05;
        relTol           0.1;
        smoother         GaussSeidel;
    };
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    p               0.3;
    U               0.7;
}

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:

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;
    }
}

And the corresponding U file like this:

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);
    }
}

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:

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);
    }
}

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