[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

[CedricVH]

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:

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

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

[josp]

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?

[matejfor]

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

[alberto]

Originally Posted by matejfor

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

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

gradSchemes
{
    default         cellLimited Gauss linear 1;
}

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

div(phi,U)      Gauss linearUpwindV cellLimited Gauss linear 1;

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,

[CedricVH]

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

[alberto]

Originally Posted by CedricVH

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.

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

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.

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

Best,

[CedricVH]

Originally Posted by alberto

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



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.

[alberto]

Originally Posted by CedricVH

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.

That's great.

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.

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

Best,