[ilpaso]

Hi,
I'm trying to set up a case with simpleFoam solver without turbolence.

The geometry is a pipe with a diameter of 3mm and a length of 117mm.


U:
The volume flow rate at inflow is 5 ml/s, so the velocity field is

inflow
{
   type fixedValue;
   value uniform (0.707355 0 0)
}
outflow
{
   type zeroGradient;
}
wall
{
   type fixedValue;
   value uniform (0 0 0);
}

nu:
The density rho is 1065kg/m^3 (blood)
The dynamic viscosity is 3.4e-3 Pa*s (blood)
So my cinematic viscosity is:

nu    nu [0  2 -1 0 0 0 0 ] 3.2e-6

p:
pressure boundary condictions are

inflow
{
   type zeroGradient;
}
outflow
{
   type fixedValue;
   value uniform 0;
}
wall
{
   type zeroGradient;
}

turbolenceProperties:
I'd like to use a laminar model, so I set:

simulationType laminar;
turbolenceModel laminar;
tubolence off;

The solution doesn't converge but if I change the inflow velocity to

inflow
{
   type fixedValue;
   value uniform (0.707355e-3 0 0)
}

and the nu to:

nu    nu [0  2 -1 0 0 0 0 ] 3.2e-3

the solution converges after 43 iterations


What is wrong in my case?

[alexeym]

Hi,

you've reduced inlet velocity and increased viscosity, so Re went from 663.15 to 0.00066315, that's why it converged quickly.

What's the number of time steps in your original simulation (with higher Re), what's in your SIMPLE dictionary?

[ilpaso]

thank you for the reply

after 1000 timesteps the residuals are horizontal but there is no convergence

this is the SIMPLE dictionary

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    residualControl
    {
        p               1e-2;
        U               1e-3;
        "(k|epsilon)"   1e-3;
    }

[alexeym]

Well...

1. what are the values of residuals for pressure and velocity (flat ones)?
2. I suppose there's no error in the direction of the velocity (i.e. it really goes along the pipe not in normal direction)
3. is it cylindrical pipe? can you show checkMesh output?

[ilpaso]

1) residuals:

smoothSolver:  Solving for Ux, Initial residual = 0.000598572, Final residual = 4.88836e-05, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 0.00176527, Final residual = 0.000138654, No Iterations 3
smoothSolver:  Solving for Uz, Initial residual = 0.00174763, Final residual = 0.000137999, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.00861638, Final residual = 0.000187439, No Iterations 2
time step continuity errors : sum local = 2.47772e-08, global = 1.5928e-10, cumulative = -8.88454e-05

2) yes, the direction is ok

3) checkMesh output:

Mesh stats
    points:           2387
    faces:            9722
    internal faces:   8578
    cells:            3660
    faces per cell:   5
    boundary patches: 3
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     0
    prisms:        3660
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing 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                  
    wall                900      930      ok (non-closed singly connected)  
    inflow              122      77       ok (non-closed singly connected)  
    outflow             122      77       ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (0 -1.4918 -1.5) (50 1.4918 1.5)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (6.36006e-18 3.39489e-18 -1.83353e-18) OK.
    Max cell openness = 1.57376e-16 OK.
    Max aspect ratio = 16.779 OK.
    Minimum face area = 0.0297778. Maximum face area = 1.11068.  Face area magnitudes OK.
    Min volume = 0.0496099. Max volume = 0.217021.  Total volume = 350.847.  Cell volumes OK.
    Mesh non-orthogonality Max: 24.1209 average: 6.15059
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.333255 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

[alexeym]

What's in your fvSolution? What if you switch GAMG/smoothSolver to PCG/PBiCG?

[ilpaso]

What means change GAMG/smoothSolver to PCG/PBiCG is obscure for me (I've to study the solver before run the simulation!! I know).
Please can you help me to set up the dictionary?

The second question is: I've another geometry with a variable diameter. The lower value is 1.2mm and the reynolds at that section is about 1600. I think I've to use a turbolence model. How change the dictionaries in order to activate this feature?

Thank you very much

this is my fvSolution.

solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;

        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        solver          GAMG;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;

        tolerance       1e-06;
        relTol          0;
    }

    "(U|k|epsilon)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }

    "(U|k|epsilon)Final"
    {
        solver          PBiCG;
        preconditioner  DILU;

        tolerance       1e-05;
        relTol          0;
    }
}

PIMPLE
{
    nOuterCorrectors 4;
    nCorrectors     1;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    residualControl
    {
        p               1e-2;
        U               1e-3;
        "(k|epsilon)"   1e-3;
    }
}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        k               0.7;
        "epsilon.*"     0.7;
    }
}

cache
{
    grad(U);
}

[alexeym]

Well,

you can just take fvSchemes/fvSolution dictionaries from $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily.

In the same tutorial you'll find all necessary modifications to run case with turbulence.

by switching from GAMG/smoothSolver to PCG/PBiCGI meant changing:

solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;

        tolerance       1e-06;
        relTol          0.05;
    }

    "(U|k|epsilon)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }
...
}

to

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.01;
    }

    U
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-05;
        relTol          0.1;
    }
...
}

[ilpaso]

It doesn't converge.
But the residuals are low! or not?
Do you think the problem is in the mesh?

these are the residuals:

DILUPBiCG:  Solving for Ux, Initial residual = 0.000595124, Final residual = 1.49812e-07, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0018561, Final residual = 3.85916e-07, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.00181595, Final residual = 4.51476e-07, No Iterations 1
DICPCG:  Solving for p, Initial residual = 0.00817729, Final residual = 7.80714e-05, No Iterations 27
time step continuity errors : sum local = 1.06637e-08, global = -2.70428e-10, cumulative = 9.34916e-06

[alexeym]

Are you sure your mesh is 3 mm in diameter? checkMesh thinks it is 3 m in diameter

Overall domain bounding box (0 -1.4918 -1.5) (50 1.4918 1.5)

(AFAIK these number are in meters)

in this case, you've got Re = UD/Nu = 663145.3125 and it's obviously non-laminar case and therefore you've got no convergence as you're running without turbulence.

[ilpaso]

now the mesh is scaled

Overall domain bounding box (0 -0.0014918 -0.0015) (0.05 0.0014918 0.0015)

but there is no convergence. I'm desperate!

DILUPBiCG:  Solving for Ux, Initial residual = 0.000565348, Final residual = 2.30843e-05, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.00543364, Final residual = 0.000228817, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.0052694, Final residual = 0.000218604, No Iterations 1
DICPCG:  Solving for p, Initial residual = 0.00364692, Final residual = 3.63526e-05, No Iterations 55
time step continuity errors : sum local = 0.000746084, global = 5.79239e-06, cumulative = -0.000574256

[ilpaso]

attached here there is the plot of the residuals

[alexeym]

I don't know what I'm doing wrong but attached case converged in 143 iterations. I've taken velocity and viscosity from your post. Two points that are different from your case:

1. Fully hexagonal mesh.
2. van Leer scheme for velocity discretisation.

(if you'd like to run attached case, you'll need Gmsh to create mesh)

[ilpaso]

thank you Alexey!
I'll try with gmsh.

Attached here you can find my polymesh directory.

[alexeym]

With attached polyMesh directory case converged in 42 iterations.

[ilpaso]

and here the case

[alexeym]

There's not convergence with "limitedLinear 1.0" scheme but if I change it to upwind it converges in 42 iteration (64 iterations with GammaV, 68 iteration with vanLeerV, 65 with linear).

[ilpaso]

with your fvSolution and fvScheme it converges in 67 iterations.
Where can I change the scheme from "limitedLinear 1.0" scheme to upwind in order to reduce the number o iterations?
A lower number of iterations is very important because I need to run the case recursively because I'll couple it with a 0D model.

[alexeym]

Your current fvSchemes:

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss limitedLinearV 1;
    div(phi,k)      bounded Gauss limitedLinear 1;
    div(phi,epsilon) bounded Gauss limitedLinear 1;
    div(phi,R)      bounded Gauss limitedLinear 1;
    div(R)          Gauss linear;
    div(phi,nuTilda) bounded Gauss limitedLinear 1;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

change it to

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss upwind;
    div(phi,k)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div(phi,R)      bounded Gauss upwind;
    div(R)          Gauss linear;
    div(phi,nuTilda) bounded Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

[ilpaso]

Ok. I'll try this evening with the different divSchemes.
Now It converges.

Thank you very much Alexey for your help!! I've lost a lot of hours with this problem!