[pete_87]

Hey all

I'm new to OpenFoam and my understanding of fluid dynamics is somewhat basic since I never got further but the basic courses.
Still, I now got the task to run simulations of an experiment (turbulent flow over backward facing step) and am trying to get the model output close to the data of the experiments. I've set up all the files to the best of my knowledge and it does run smoothly. But no matter how I tweak around with boundary conditions, timing and the inlet values that are not given (center streamline velocity is given), the simulation simply won't show turbulence, it always converges to a steady state. The flow is supposed to have Reynolds numbers of up to 10'000 and higher... It is given that I must use the RNGkEpsilon model and I been solving it with pisoFoam, pimpleFoam and simpleFoam as well (since it converged anyway...).

Below you find my files, I'd be extremely thankful for any advice where I might have gone wrong or what I missed.

Thank you all very much!

Steve

blockMeshDict:

convertToMeters 0.001; // metrics in mm

vertices // all edge points of blocks in the model
(
    (2000 0 20) // point 0
    (2500 0 20) // point 1
    (2500 0 -20) // point 2
    (2000 0 -20) // point 3
    (2000 41 20) // point 4
    (2500 41 20) // point 5
    (2500 41 -20) // point 6
    (2000 41 -20) // point 7
    (0 41 -20) // point 8
    (0 41 20) // point 9
    (0 122 20) // point 10
    (0 122 -20) // point 11
    (2000 122 -20) // point 12
    (2000 122 20) // point 13
    (2500 122 20) // point 14
    (2500 122 -20) // point 15
);

blocks
(
    hex (3 2 6 7 0 1 5 4) (150 10 1) simpleGrading (10 0.1 1) // block 0: 
    hex (7 6 15 12 4 5 14 13) (150 20 1) simpleGrading (10 10 1) // block 1
    hex (8 7 12 11 9 4 13 10) (100 20 1) simpleGrading (0.01 10 1) // block 2
);

edges
(
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (9 10 11 8)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (1 2 6 5)
        (5 6 15 14)
        );
    }
    upperWall
    {
        type wall;
        faces
        (
            (10 13 12 11)
        (13 14 15 12)
        );
    }
    lowerWall
    {
        type wall;
        faces
        (
            (0 3 2 1)
        (9 8 7 4)
        );
    }
    sideWall
    {
    type wall;
    faces
    (
        (0 4 7 3)
    );
    }
    frontAndBack
    {
        type empty;
        faces
        (
            (0 1 5 4)
        (4 5 14 13)
        (9 4 13 10)
        (3 7 6 2)
        (7 12 15 6)
        (8 11 12 7)
        );
    }
);

mergePatchPairs
(
);

controlDict:

/*--------------------------------*- 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;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     pimpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         1000;

deltaT          0.0001;

writeControl    adjustableRunTime;

writeInterval   0.001; //only last time will be written...

adjustTimeStep    yes;

maxCo        1.15; // max Courant number

maxDeltaT    1;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

functions
{
  interestpoints
  {
    type probes;
    functionObjectLibs ("libsampling.so");
    outputControl timeStep;
    outputInterval 1;
    probeLocations ((2.041 0.07954 0)(2.246 0.04059 0)); // center streamline measure point: (1.9 0.0815 0)
    fields ( U );
  }
};

fvSchemes:

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

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

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
}

fvSolution:

solvers
{
    "(p|pFinal)"
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0;//0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        cacheAgglomeration on;
        agglomerator    faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels     1;
    }


    "(U|k|epsilon|R|nuTilda)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;//0.1;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    "(U|k|epsilon)Final"
    {
        $U;
        relTol          0;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;

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

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
} 

relaxationFactors
{
    fields
    {
        p               0.3; // 0.3 originally
    }
    equations // all 0.7 originally - 0.9 crashes for standard values...
    {
        U               0.8;
        k               0.8;
        epsilon         0.8;
        R               0.8;
        nuTilda         0.8;
    }
}

PIMPLE
{
    nNonOrthogonalCorrectors 0;
    nCorrectors         2;
}

transportProperties:

transportModel  Newtonian;

nu              nu [ 0 2 -1 0 0 0 0 ] 1e-06;// water!
rho        rho [ 1 -3 0 0 0 0 0 ] 1e03;

CrossPowerLawCoeffs
{
    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
    m               m [ 0 0 1 0 0 0 0 ] 1;
    n               n [ 0 0 0 0 0 0 0 ] 1;
}

BirdCarreauCoeffs
{
    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
    k               k [ 0 0 1 0 0 0 0 ] 0;
    n               n [ 0 0 0 0 0 0 0 ] 1;
}

turbulenceProperties:

simulationType    RASModel;

RASProperties:

RASModel        RNGkEpsilon;
turbulence      on;
printCoeffs     on;
RNGkEpsilonCoeffs
{
    Cmu             0.0845;
    C1              0.7194;
    C2              0.7194;
    sigmak          1.42;
    sigmaEps        1.68;
    eta0            4.38;
    beta            0.012;
}

Initial values: ***********
U:

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.0997 0 0);//uniform (0.086466666 0 0); //using average velocity as input  0.0997 - leads to wanted centerline velocity
    }

    outlet
    {
        type            zeroGradient;
    }

    upperWall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    lowerWall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    sideWall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    frontAndBack
    {
        type            empty;
    }
}

p:

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
    type            fixedValue;
        value           uniform 0;
    }

    upperWall
    {
        type            zeroGradient;
    }

    lowerWall
    {
        type            zeroGradient;
    }
    
    sideWall
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}

epsilon: (calculated from velocity and hydraulic diameter of pipe)

dimensions      [0 2 -3 0 0 0 0];

internalField   uniform 5.735756440858130e-06;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 5.735756440858130e-06;
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            epsilonWallFunction;
        value           uniform 5.735756440858130e-06;
    }
    lowerWall
    {
        type            epsilonWallFunction;
        value           uniform 5.735756440858130e-06;
    }
    sideWall
    {
        type            epsilonWallFunction;
        value           uniform 5.735756440858130e-06;
    }
    frontAndBack
    {
        type            empty;
    }
}

k:

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 3.471953550469051e-05;

boundaryField
{
    inlet
    {
    type            fixedValue;
        value           uniform 3.471953550469051e-05;
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            kqRWallFunction;
        value           uniform 3.471953550469051e-05;
    }
    lowerWall
    {
        type            kqRWallFunction;
        value           uniform 3.471953550469051e-05;
    }
    sideWall
    {
        type            kqRWallFunction;
        value           uniform 3.471953550469051e-05;
    }
    frontAndBack
    {
        type            empty;
    }
}

nut:

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            calculated;
        value           uniform 0;
    }
    outlet
    {
        type            calculated;
        value           uniform 0;
    }
    upperWall
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    lowerWall
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    sideWall
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    frontAndBack
    {
        type            empty;
    }
}

nuTilda:

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    outlet
    {
        type            zeroGradient;
    }

    upperWall
    {
        type            zeroGradient;
    }

    lowerWall
    {
        type            zeroGradient;
    }

    sideWall
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}

[Lieven]

Hi Steve,

It is very logical that you obtain a steady state result when using RANS turbulence modeling. I recommend you to do some reading in the differences between RANS and LES. A quick wikipedia search on RANS would already explain you that:

The RANS equations are primarily used to describe turbulent flows. These equations can be used with approximations based on knowledge of the properties of flow turbulence to give approximate time-averaged solutions to the Navier-Stokes equations.

So you see that RANS-models result in a time-averaged solution, hence you won't (or at least should not) see fluctuations...

Cheers,

Lieven

[pete_87]

Hm, I almost suspected something like that, but it seemed contradictory to the task. But maybe I'm wrong again: if this basically yields me only the averaged velocities, then I have no way to calculate reynolds shear stresses, do I? The way I understood it those are caused by the fluctuations and I need the latter to calculate them, am I wrong?

Thank you for your answer!

[pete_87]

So upon that answer I went on another extensive google search trail and discovered the openFoam utilities for Reynolds stress and stress tensors. How does openFoam calculate those? Can it actually be derived from the kappa and epsilon fields?

[pete_87]

Follow up question: I need the Reynolds stress at specific locations but apparently the R utility can't do that. Combining the R field with "probes" doesn't work either. Can I interpolate it from the R field output or something?