[jck]

Good morning,

i created a 2D NACA airfoil (angle of attack = 14°) with blockMesh and already checked it with checkMesh. Everything seems to be ok.

The calculation (kOmegaSST) is runnning without any mistakes and the results look right at the first sight but the pressure at the suction side, especially at the tip is to low.

Physically, the minimum pressure coefficient should be around -6, but i get a minimum of around -1 with a fine grid (4800 points on the airfoil). The more i refine the grid, the pressure value scales down in the wrong direction.

I changed the grid many times but I'm not sure if it's a mesh mistake or a definition mistake in one of the system-/, constant/- or 0/-files.

Does anyone have an idea or had the same problems?

Thanks in advance

Claus

[teuk]

Hi jck,

I think people will need a little more information. What are your BCs on the airfoil? What is your y+? etc...


regards,
teuk

[jck]

Hey, that's true. I'm working with wallfunctions.

k:
boundaryField
{
wing 
    {
        type            kqRWallFunction;
        value           uniform 0.0096;
    }
}

nut:
boundaryField
{
    wing 
    {
        type            nutkWallFunction;
        Cmu                0.09;
        kappa            0.41;
        E                9.8;
        value           uniform 0;
    }
}

omega:
boundaryField
{
    wing 
    {
        type            omegaWallFunction;
        Cmu                0.09;
        kappa            0.41;
        E                9.8;
        beta1            0.075;
        value           uniform 9.8;
    }
}

p:
boundaryField
{
    wing 
    {
        type            zeroGradient;
    }

}

U:
boundaryField
{
    wing 
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

fvSchemes:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
    grad(U)         Gauss linear;
}

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

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p;
}

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

solvers
{
    p
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;
        minIter             5;
        maxIter             100;
        smoother         GaussSeidel;
        nPreSweeps       0;
        nPostSweeps      2;
        nFinestSweeps     2;
        cacheAgglomeration on;
        agglomerator     faceAreaPair;
        nCellsInCoarsestLevel 0;
        mergeLevels      1;
    }

    pFinal
    {
        $p;
        tolerance       1e-7;
        relTol          0.01;
    }

    U
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-8;
        relTol           0.01;
        nSweeps          1;
    }

    k
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-8;
        relTol           0.01;
        nSweeps          1;
    }

    omega
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-8;
        relTol           0.01;
        nSweeps          1;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

RASProperties:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      RASProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

RASModel            kOmegaSST;

turbulence          on;

printCoeffs         on;

[teuk]

Hi Claus,

there is a bunch of possible reasons for this.

(4800 points on the airfoil)

What is your chordlength in relation to this number of points?

There is a unlikely possibilty to have a to fine mesh. -> What is your cellsize close to the airfoil (and y+)?

Is this a low-Re simulation?

Is your simulation in stall condition? -> Do you expect a seperation bubble at this point?

By "tip" you mean the Leading Edge?


Since your results are getting better with a coarser mesh probably this is the solution. Unfortunately its easier to refine a mesh then coarsen it :-P

I'm not sure about GAMG for p. Could "blur" your field a little...



regards,
teuk

[jck]

Hey teuk,

-chordlength= 1m

-Re=1.85e+6

-yes the airfoil should be in stall condition

-with "tip" i mean "leading edge"

Is there a function for the cell size in y-direction or for the y+ value?

I will try to calculate it with a coarser mesh on monday and tell you about the results because i only have access to OF during the week.

Regards
Claus

[jck]

Hey,

i took the function yPlusRAS to calculate the y+ value on the wing.
I get values from 0.7 in low-speed area to ~30 in higher speed area (in the finer mesh calculation).

Regards
Claus

[Kire]

Hi Jck

If it is in stall condition, you probably should consider a 3D model. You can use periodic boundary conditions in the span-wise direction. Two-dimensional models are not very good at computing forces in massive flow separation.

Sometimes 2D unsteady SST or SA return reasonably good force values for flow past a cylinder at a few Reynolds numbers only. But looking at a wide range of Re = 3900 to 10^6, forces computed by 2D RANS are usually off on this geometry.

By the way, with SST you can apply a very small wall distance so that y+~1.

[jck]

Hey Kire,

the NASA made a similar calculation with a 2D grid and got good results.

How can i apply a very small wall distance without waiting 6 days for the results of a very fine mesh?
I read something about the function refineWallLayer.
Could it work with it?

[Kire]

Hi Jck

1) Can you link me to the relevant NASA page?

2) Since you are running a transient solver, consider faster linear system solver settings:

solvers
{
    p
    {
        solver           GAMG;
        tolerance        1e-6;
        relTol           1e-3;
        smoother         GaussSeidel;
        nPreSweeps       0;
        nPostSweeps      2;
        cacheAgglomeration on;
        agglomerator     faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels      1;
    }

    pFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-06;
            relTol           0;
            smoother        DICGaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nFinestSweeps   2;
            cacheAgglomeration true;
            nCellsInCoarsestLevel 10;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }
        tolerance       1e-06;
        relTol           0;
    }

    U
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-5;
        relTol            0;
        nSweeps          1;
    }

    k
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-6;
        relTol             0;
        nSweeps          1;
    }

    omega
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-6;
        relTol             0;
        nSweeps          1;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

3) I have never used the function refineWallLayer. What mesh generator do you use? It can be difficult to create low-Re cells in snappyHex.

[jck]

Hey,

1) http://turbmodels.larc.nasa.gov/naca4412sep_val.html

Here is link to the refineWallLayer function:
https://openfoamwiki.net/index.php/RefineWallLayer

I use blockMesh to generate my mesh in this high Re calculation

[Kire]

Hi Jck

1) I am certain you can control grid clustering easily in blockMesh. You can alter simpleGrading to control the wall distance. It is not necessary to use refineWallLayer.

2) It seems to me that the NASA validation was done in a compressible steady state solver at Mach = 0.09, where compressible effects are believed to not have much influence on the results. You can use incompressible/compressible SIMPLEC in this case.

[jck]

Hi Kire,

1) yeah it's possible to control the wall distance with simpleGrading. I tried both ways. Sometimes the aspect ratio or the skewness were increasing.

2) Do you thing that the compressible solver could solve the problem with the pressure coefficient? I would be surprised because of the small Mach number.

3) Is it right, that there are two ways to deal with y+:

3.1) no wallfunctions -> 1-10 cells in viscous sublayer, y+ < 1

3.2) wallfunctions -> coarse grid, 30 < y+ < 300 (log law layer)

[Kire]

Hi Jck

1) I don't think aspect ratio is your biggest concern. A mesh can and will solve even when aspect ratio in the first layer is larger than 1000, the value at which OpenFOAM starts showing warnings. If it is a steady state solver, you can use leastSquares or its extended versions when having large skewness to try and preserve result accuracy.

2) I believe a pressure-based compressible solver will give you reasonably good results at such Mach number. How about trying it out?

3)
a. No wallfunction for nut when y+～1. Simply use nut = 0 or nutLowReWallFunction ( nut = 0 ). Ten layers might be too many, try 5. Btw, I use 20 or more prism layers to cover the entire boundary layer.

b. Yes, wall function as in wall function for nut, but still need 5~8 or more prism layers to cover the entire boundary layer.

[jck]

Hey Kire,

1)I will tried it with the compressible solver and without wallfunctions.

2)Do i use the right definitions for k and omega:

k=(3/2)*(U*Tu), Tu=turbulence intensity

omega=(rho*k)/(mut), mut=turbulent viscosity, rho=density ?

With these definitions i get: k=0.00081656 and omega=544,3729

[jck]

Good morning,

there is no big difference between the compressible and the incompressible calculation. The pressure coefficient is wrong in both cases.

1) What's about the boundary conditions? I don't believe that it's a mesh problem.

2) What is the best way to calculate y+? I read that yPlusRAS just calculates y* and not y+.

Regards