[ottodrissen]

Hi Everybody,

I’d like to ask how I can get SRFSimpleFoam (openFoam 3.0) to converge on a wind turbine blade in a periodic sector. Even in pure laminar flow it doesn’t work out. The first 1000 or even 1500 iterations it goes fine, but after that the residuals increase and the solver either breaks down or keeps oscillating around high values.

I tried to play with the relaxation factors, but that didn’t help, if the pressure residual goes down, the velocity residual goes up and vice versa (running laminar only for now).

• What can I do to improve convergence?
• Are my boundary conditions correctly stated?
• Is there something I can change to the schemes?

I’m trying to run a 120° sector (size 5·Drotor in all distances) with periodicity on the side walls, a fixed velocity inlet and pressure outlets on the back face and the outer shroud (the turbine’s rotor brakes the flow, so the air should be escaping through the outer shroud). The fluid domain is sketched in the picture Domain.png that I attached.

The mesh is of a more than reasonable quality.

checkMesh results:

Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           2056944
    faces:            6072059
    internal faces:   5971939
    cells:            2007640
    faces per cell:   5.9990825
    boundary patches: 7
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     2005798
    prisms:        1812
    wedges:        0
    pyramids:      30
    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                  
    AXIS                6908     7064     ok (non-closed singly connected)  
    OUTERSHROUD         6980     7086     ok (non-closed singly connected)  
    INLET               12760    13095    ok (non-closed singly connected)  
    BLADE               16632    16654    ok (non-closed singly connected)  
    OUTLET              12760    13095    ok (non-closed singly connected)  
    PERIODICPORT        22040    22407    ok (non-closed singly connected)  
    PERIODICSTAR        22040    22407    ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-165 -95.2628021 0.299999982) (165 95.2628021 110)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (1.14666433e-016 4.19347538e-015 -6.4542397e-015) OK.
    Max cell openness = 2.26783601e-015 OK.
    Max aspect ratio = 912.558483 OK.
    Minimum face area = 5.59392527e-007. Maximum face area = 160.838438.  Face area magnitudes OK.
    Min volume = 1.86467905e-009. Max volume = 330.931674.  Total volume = 3711995.15.  Cell volumes OK.
    Mesh non-orthogonality Max: 80.3005844 average: 20.0762812
   *Number of severely non-orthogonal (> 70 degrees) faces: 2101.
    Non-orthogonality check OK.
  <<Writing 2101 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
    Max skewness = 2.16006169 OK.
    Coupled point location match (average 9.20047212e-007) OK.

Mesh OK.

End

To check for convergence I monitor the residuals, the thrust and the torque on the rotor.

I used the following fvSchemes:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.x                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default                          none;
    div(phi,Urel)                    bounded Gauss linearUpwind limited;
    div(phi,k)                       bounded Gauss upwind;
    div(phi,epsilon)                 bounded Gauss upwind;
    div(phi,omega)                   bounded Gauss upwind;
    div(phi,R)                       bounded Gauss upwind;
    div(R)                           Gauss linear;
    div(phi,nuTilda)                 bounded Gauss upwind;
    div((nuEff*dev2(T(grad(Urel))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

// ************************************************************************* //

And these solvers (fvSolution):

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.x                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-07;
        relTol          0.00;
    }

    Urel
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0.0;
    }

    k
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0.0;
    }

    epsilon
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0.0;
    }

    omega
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0.0;
    }

    R
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0.0;
    }

    nuTilda
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0.0;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    pRefCell                 0;
    pRefValue                0.0;

    residualControl
    {
        p               1e-8;
        Urel            1e-8;
        k               1e-8;
        epsilon         1e-8;
        omega           1e-8;
        R               1e-8;
        nuTilda         1e-8;
    }
}

relaxationFactors
{
    p               0.3;
    Urel            0.7;
    k               0.5;
    epsilon         0.5;
    omega           0.5;
    R               0.5;
    nuTilda         0.5;
}


// ************************************************************************* //

I tried to play with the relaxation factors, but that didn’t help, if the pressure residual goes down, the velocity residual goes up and vice versa (running laminar only for now).

• What can I do to improve convergence?
• Are my boundary conditions correctly stated?
• Is there something I can change to the schemes?

Thanks,

Otto

[ottodrissen]

48 views by now, but no one seems to have any advice.

Being a beginner with openFoam I thought my question to be quite simple.

I wonder whether the question is more difficult than I thought or that I'm just among other beginners like me.

[akidess]

Your mesh is quite poor, so reduce the non-orthogonality of your mesh. If you *really* cannot do so, you need to limit your gradients. Have a look at the following presentation, especially the slides starting at slide 48:
http://www.dicat.unige.it/guerrero/o...sandtricks.pdf

[ottodrissen]

Thank you very much Anton!! Finally someone who answers.

To be honest, I thought my mesh to be very good, but it seems I’m wrong.
The mesh is elliptical in a bar around the blade (see BladeAirfoil.png) and the rest of the elements to fill the domain are as squared as a cube as a matter of speaking.

There are only 2101 severely non orthogonal faces of 6072059 faces in total, that’s 0.03% of the mesh. From these faces, the maximum seems to be 80.3°, but the mesh average of 20° is quite good (at least I thought so).

The non orthogonal faces are located at a certain distance on top of the blade (see NonOrthogonalElements.png).

Ok, it’s not clean laminar parallel flow overthere, but neither is it a zone where I expect gradients to be very high.

Anyway, thanks again for your advice, I will let you know if it works better with the tips from the PDF as soon as I can start the computation again.

Otto

[arjun]

cant help you with openfoam but your mesh is quite acceptable.

Originally Posted by ottodrissen

Thank you very much Anton!! Finally someone who answers.

To be honest, I thought my mesh to be very good, but it seems I’m wrong.
The mesh is elliptical in a bar around the blade (see BladeAirfoil.png) and the rest of the elements to fill the domain are as squared as a cube as a matter of speaking.

There are only 2101 severely non orthogonal faces of 6072059 faces in total, that’s 0.03% of the mesh. From these faces, the maximum seems to be 80.3°, but the mesh average of 20° is quite good (at least I thought so).

The non orthogonal faces are located at a certain distance on top of the blade (see NonOrthogonalElements.png).

Ok, it’s not clean laminar parallel flow overthere, but neither is it a zone where I expect gradients to be very high.

Anyway, thanks again for your advice, I will let you know if it works better with the tips from the PDF as soon as I can start the computation again.

Otto

[akidess]

Originally Posted by arjun

cant help you with openfoam but your mesh is quite acceptable.

So are you suggesting him to use a different solver? Because you'll get unbounded results using corrected schemes on that mesh with OpenFOAM in my experience.

[arjun]

Originally Posted by akidess

So are you suggesting him to use a different solver? Because you'll get unbounded results using corrected schemes on that mesh with OpenFOAM in my experience.

This is why I said, I cant help with openFOAM part.

suggesting to improve the mesh is usually good advice because in real life applications making perfect mesh is not possible so there is always room for improvement.

Given this thing that making perfect mesh in general is not possible, i think the quality of the mesh that I expect solver to be able to run. This mesh in my opinion is good enough that a CFD solver shall be able to perform simulation on. (This is what I meant).
Now whether openFOAM could do it is another matter. I dont know the answer to this, but with my experience with cfd solver I would expect that openfoam shall be able to run it.

I know that the same mesh fluent for example would have no issues running simulation with.

[arjun]

Originally Posted by ottodrissen

There are only 2101 severely non orthogonal faces of 6072059 faces in total, that’s 0.03% of the mesh. From these faces, the maximum seems to be 80.3°, but the mesh average of 20° is quite good (at least I thought so).

I think, first try first order upwind schemes. Then start from here.

[akidess]

OpenFOAM can do it too, if you adapt the discretization (see my post above) and live with the reduced accuracy. I'm confident Fluent does this as well, except you might perhaps not realize it because it's done behind the scenes.

[arjun]

Originally Posted by akidess

OpenFOAM can do it too, if you adapt the discretization (see my post above) and live with the reduced accuracy. I'm confident Fluent does this as well, except you might perhaps not realize it because it's done behind the scenes.

I worked at Cd adapco and at least for CCM I know what is done behind the scenes. (Expect fluent to do the same as ccm is developed by people who were fluent developers before).

There is nothing special done behind the scenes in ccm. There exist gradient limiter that only limits if bounds are broken. The same I guess shall be available to openfoam too.

Stablity comes from combination of things, this is why each solver behaves differently.

PS: Edited to add. Cant speak for fluent but for ccm, the amount of stress that is given to accuracy of solver is very high in ccm. Stablity is second concern there so I expect ccm to be less stable actually.

[ottodrissen]

Back again,

So far, changing the schemes didn't work, unfortunately. Nevertheless thanks for the PDF, it's full with good advice.

I agree with arjun that it's hard to get OpenFoam to converge, at least in comparison with Fluent. Yes, maybe Fluent does the same tricks under the bonnet that one should apply in OpenFoam (but automatically), but the info that's available for OpenFoam's schemes, their physical significance, and mathematical effects can at best be called extremely scarce. The link

http://cfd.direct/openfoam/user-guide/fvschemes/

for example doesn't even cover a quarter of all options that I see passing by in tutorial cases and on top of that it's only descriptive, not explanatory.

I will try to improve the mesh (that's gonna take time, as I only have some spare hours for this).

Otto

[arjun]

Originally Posted by ottodrissen

Back again,

So far, changing the schemes didn't work, unfortunately. Nevertheless thanks for the PDF, it's full with good advice.

I agree with arjun that it's hard to get OpenFoam to converge, at least in comparison with Fluent. Yes, maybe Fluent does the same tricks under the bonnet that one should apply in OpenFoam (but automatically), but the info that's available for OpenFoam's schemes, their physical significance, and mathematical effects can at best be called extremely scarce. The link

http://cfd.direct/openfoam/user-guide/fvschemes/

for example doesn't even cover a quarter of all options that I see passing by in tutorial cases and on top of that it's only descriptive, not explanatory.

I will try to improve the mesh (that's gonna take time, as I only have some spare hours for this).

Otto

The main problem here is (already pointed out) gradients in these cells where this skew problem happen are wrong. Trying to apply second or higher order interpolations could create problems.

To remedy this solvers usually use gradient limiters, how they apply this makes most difference.

What you actually need is a solver that automatically switches off higher order interpolations on these faces where skew is bad. These are handful of cells and you should not be making new mesh because of them. You accept the trade off and move on.

This is exactly how FVUS/wildkatze solver that i created is designed. If the skew is bad, interpolations are moving to first order scheme and the degree of blending is dependent on skew.

Off course on the same mesh, in second order schemes my solver shows more convergence and stablity compared to fluent and ccm. (But as far as I know mine is only solver that does it, ie skew based blending.)

[arjun]

Here is a sample torture test for my solver.

Second order upwind and run 500 iterations.

Delibrately created extremely bad mesh to test solver robustness.

Ensight gold files here along with residuals in .csv

http://www.dravvya.co.in/outgoing/skDir.tar

Have fun.