[Dan1788]

Dear Foamers,

I have started doing an LES of pipe flow using pimpleFoam with cyclic boundary conditions in the streamwise direction to generate turbulent statistics and feed them into a premixed jet flame simulation. Some parameters of the pipe:
D = 15mm
L = 5*D
Ubar = 9.2 m/s
Re = 9000
Re_tau = 287
Turbulence model = oneEqEddy
delta = vanDriest

In the wall normal direction I have resolved the mesh with a y+ = 1. My axial mesh has a resolution of \Delta-z+ = 15. The mesh front view is attached herewith

https://www.dropbox.com/s/3ed02j8ojj..._mesh.png?dl=0

Here is my fvOptions file for the pressureGradientExplicitSource

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

momentumSource
{
    type            pressureGradientExplicitSource;
    active          on;            //on/off switch
    selectionMode   all;       //cellSet // points //cellZone

    pressureGradientExplicitSourceCoeffs
    {
        fieldNames  (U);
        Ubar        (0 0 9.2);
    }
}

The fvSchemes files is as follows:

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

ddtSchemes
{
    default         CrankNicolson 1;
}

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

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,B)      Gauss limitedLinear 1;
    div(B)          Gauss linear;
   // div(phi,nuTilda) Gauss limitedLinear 1;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         none;
    laplacian(nuEff,U) Gauss linear corrected;
    laplacian((1|A(U)),p) Gauss linear corrected;
    laplacian(DkEff,k) Gauss linear corrected;
    laplacian(DBEff,B) Gauss linear corrected;
   // laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
    interpolate(U)  linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
}


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

The fvSolution file is as follows:

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

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

    pFinal
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0;
    }

    "(U|k)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-05;
        relTol          0.1;
    }

    "(U|k)Final"
    {
        $U;
        tolerance       1e-05;
        relTol          0;
    }
}

PIMPLE
{
    nOuterCorrectors 2;
    nCorrectors     2;
    nNonOrthogonalCorrectors 1;
    pRefCell        1001;
    pRefValue       0;
}


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

I ran my code for 0.05 seconds (approximately 5 flow through times without averaging the results) and looked at the U-z (streamwise velocity profile) at the end time using ParaView. The figure is attached herewith.

My question is why does the instantaneous velocity have such a smooth laminar like structure? Why cant I see any eddies ??

Would be glad if anyone can help

[RodriguezFatz]

Hey Daniel,

1) Please post some of the log output of your run.
2) What is the x+ resolution of your mesh?
3) With y+=1, I guess, you don't use any wall functions?

[Dan1788]

Hi Philipp,

I just realized that the reason I was not getting any turbulence generated was because I wasnt specifying any manually generated turbulence at the inlet (pardon me. really stupid on my part )

I am now trying to generate turbulence using the boxTurb16 utility and then map it to my pipe for initial turbulence. A couple of questions I have are:

1. Why does the box grid in the boxTurb utility have to be in powers of 2^k?

2. In the transportProperties and the turbulenceProperties file are defined variables :

Ea              10;
k0              5;
UOsigma         0.090295;
UOalpha         0.81532;
UOKupper        10;
UOKlower        7;

How are these parameters chosen for a simulation ? Is it just trial and error?

[RodriguezFatz]

Daniel, you don't need to initialize the turbulence... it should be just much faster to get convergence. I didn't initialize it and get some nice LES of a periodic pipe.
I can recommend to track the kinetic energy in the log-output. If you plot it over time, you can see, if your simulation is converged.

I will post some pictures in a minute.

[RodriguezFatz]

This is a picture of a pipe flow just as you do, but with higher Re and wall modeled LES. It was initialized without any turbulence.
U.jpg

This is a picture of the kinetic energy (perpendicular to the main flow) to track the convergence of the case. You see, it is not fully converged after 40000 time steps. This means after 1.6s or about 100 flow through times. Of course, this is due to very bad initial guess, but you see, this can take some time.
kEn.jpeg

I track the energy by putting this in my controlDict:

functions
(
AverageResolvedKineticEnergyVxVy
    {
        type swakExpression;
        valueType cellZone;
        zoneName FLUID;
        accumulations (
            average
        );
        expression "U.x*U.x + U.y*U.y";
        verbose true;
    }
);

I only use Ux and Uy for tracking, because I only want to track the kinetic energy of the turbulent structures and Uz is my flow direction. I was too lazy to calculate the average Uz at each point. So this is a fast and good work around.

[Dan1788]

Thanks for your post Philipp. I am really curious to know how you achieve turbulence without initializing any at the beginning of the simulation.

What turbulence model do you use and what is the Reynolds number for your flow??

Is the tradeoff with your technique a longer run time? (I ask because the result I posted in my first post was for just 5 flow through times without initializing turbulence and I couldnt observe any turbulent structures in my flowfield )

[RodriguezFatz]

Daniel,
1) Turbulence is produced by the flow. So there is no need for any special initialization to get turbulence at all. The only exception that I know is, when you solve equations for turbulence and the production rate is proportional to some turbulence parameter (such as "k"). If you initialize "k" with zero everywhere, you are doomed, because the production rate will be zero for all time. I was using Smagorinsky model.
2) Two post above (#4) I wrote that it took more than 100 flow through times and it still is not converged.
3) Your problems can have many reasons. You need to answer the questions of post #2 if you want any help to solve them. Maybe the answer is just: "Wait even longer".

[Dan1788]

Hi Philipp,

Refering to your post#2 here is the log output of my run:

Courant Number mean: 0.318167 max: 0.49787
deltaT = 1.04167e-05
Time = 0.0999792

PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.0160958, Final residual = 0.000505205, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0150424, Final residual = 0.000417883, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 1.10966e-05, Final residual = 1.06779e-07, No Iterations 1
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9431
DICPCG:  Solving for p, Initial residual = 0.362879, Final residual = 0.0159966, No Iterations 16
DICPCG:  Solving for p, Initial residual = 0.016784, Final residual = 0.000816065, No Iterations 104
time step continuity errors : sum local = 9.65187e-11, global = 1.35337e-11, cumulative = -1.9175e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9424
DICPCG:  Solving for p, Initial residual = 0.0733988, Final residual = 0.00325182, No Iterations 6
DICPCG:  Solving for p, Initial residual = 0.00358734, Final residual = 0.000172535, No Iterations 223
time step continuity errors : sum local = 2.8346e-11, global = 1.35337e-11, cumulative = -1.91614e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9424
PIMPLE: iteration 2
DILUPBiCG:  Solving for Ux, Initial residual = 0.000610097, Final residual = 4.07385e-06, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.000570686, Final residual = 3.67736e-06, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 3.10073e-09, Final residual = 3.10073e-09, No Iterations 0
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9424
DICPCG:  Solving for p, Initial residual = 0.00980145, Final residual = 0.000418088, No Iterations 4
DICPCG:  Solving for p, Initial residual = 0.000437694, Final residual = 2.17607e-05, No Iterations 245
time step continuity errors : sum local = 1.53838e-11, global = 1.35337e-11, cumulative = -1.91479e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9425
DICPCG:  Solving for p, Initial residual = 0.00173541, Final residual = 8.43401e-05, No Iterations 3
DICPCG:  Solving for p, Initial residual = 8.7379e-05, Final residual = 9.85646e-07, No Iterations 271
time step continuity errors : sum local = 1.36178e-11, global = 1.35337e-11, cumulative = -1.91344e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9425
DILUPBiCG:  Solving for k, Initial residual = 2.96381e-05, Final residual = 3.65235e-07, No Iterations 1
bounding k, min: 0 max: 0.122827 average: 0.0593156
ExecutionTime = 8980.11 s  ClockTime = 9035 s

Courant Number mean: 0.318167 max: 0.497882
deltaT = 1.04167e-05
Time = 0.0999896

PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.015843, Final residual = 0.000523452, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0143468, Final residual = 0.000436435, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 1.10947e-05, Final residual = 1.06753e-07, No Iterations 1
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9421
DICPCG:  Solving for p, Initial residual = 0.340204, Final residual = 0.0165366, No Iterations 16
DICPCG:  Solving for p, Initial residual = 0.0178719, Final residual = 0.000870718, No Iterations 226
time step continuity errors : sum local = 8.40756e-11, global = -7.37073e-12, cumulative = -1.91417e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9416
DICPCG:  Solving for p, Initial residual = 0.0748232, Final residual = 0.00314183, No Iterations 7
DICPCG:  Solving for p, Initial residual = 0.00334815, Final residual = 0.000161511, No Iterations 89
time step continuity errors : sum local = 2.30996e-11, global = -7.37072e-12, cumulative = -1.91491e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9416
PIMPLE: iteration 2
DILUPBiCG:  Solving for Ux, Initial residual = 0.000584867, Final residual = 4.24894e-06, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.000545387, Final residual = 3.90244e-06, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 3.04409e-09, Final residual = 3.04409e-09, No Iterations 0
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9416
DICPCG:  Solving for p, Initial residual = 0.00976526, Final residual = 0.000399203, No Iterations 4
DICPCG:  Solving for p, Initial residual = 0.000416475, Final residual = 2.03138e-05, No Iterations 242
time step continuity errors : sum local = 9.11331e-12, global = -7.37073e-12, cumulative = -1.91565e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9416
DICPCG:  Solving for p, Initial residual = 0.00170644, Final residual = 8.18023e-05, No Iterations 3
DICPCG:  Solving for p, Initial residual = 8.49119e-05, Final residual = 9.48026e-07, No Iterations 312
time step continuity errors : sum local = 7.45187e-12, global = -7.37073e-12, cumulative = -1.91638e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9416
DILUPBiCG:  Solving for k, Initial residual = 2.96352e-05, Final residual = 3.65196e-07, No Iterations 1
bounding k, min: 0 max: 0.122824 average: 0.0593149
ExecutionTime = 8981.2 s  ClockTime = 9036 s

Courant Number mean: 0.318167 max: 0.497894
deltaT = 1.04167e-05
Time = 0.1

PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.0170933, Final residual = 0.000506882, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0157762, Final residual = 0.000423149, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 1.10933e-05, Final residual = 1.06765e-07, No Iterations 1
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9414
DICPCG:  Solving for p, Initial residual = 0.362403, Final residual = 0.0157429, No Iterations 16
DICPCG:  Solving for p, Initial residual = 0.0165154, Final residual = 0.000787541, No Iterations 230
time step continuity errors : sum local = 8.10639e-11, global = -6.68959e-12, cumulative = -1.91705e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9408
DICPCG:  Solving for p, Initial residual = 0.0734679, Final residual = 0.00324997, No Iterations 6
DICPCG:  Solving for p, Initial residual = 0.00359747, Final residual = 0.000178054, No Iterations 238
time step continuity errors : sum local = 2.2006e-11, global = -6.68959e-12, cumulative = -1.91772e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9408
PIMPLE: iteration 2
DILUPBiCG:  Solving for Ux, Initial residual = 0.000607425, Final residual = 4.14948e-06, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.00057021, Final residual = 3.84362e-06, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 3.0812e-09, Final residual = 3.0812e-09, No Iterations 0
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9408
DICPCG:  Solving for p, Initial residual = 0.00974758, Final residual = 0.00039325, No Iterations 4
DICPCG:  Solving for p, Initial residual = 0.000411935, Final residual = 1.98099e-05, No Iterations 111
time step continuity errors : sum local = 8.16657e-12, global = -6.68959e-12, cumulative = -1.91839e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9408
DICPCG:  Solving for p, Initial residual = 0.00173266, Final residual = 7.96663e-05, No Iterations 3
DICPCG:  Solving for p, Initial residual = 8.25938e-05, Final residual = 9.99876e-07, No Iterations 308
time step continuity errors : sum local = 6.77497e-12, global = -6.68959e-12, cumulative = -1.91906e-08
Pressure gradient source: uncorrected Ubar = 9.2, pressure gradient = 33.9408
DILUPBiCG:  Solving for k, Initial residual = 2.96323e-05, Final residual = 3.65171e-07, No Iterations 1
bounding k, min: 0 max: 0.122821 average: 0.0593141
ExecutionTime = 8982.74 s  ClockTime = 9038 s

End

2. If by x+ you mean streamwise resolution, then it is around 15.
3. My minimum y+ resolved at the wall is 1. I am not using any wall functions.

In fact I had made a mistake of assigning k to 0 using oneEqEdy model so that there wont be any production of turbulence. But later I changed the turbulence model to Smagorinsky that does not use a an k equation but is calculated. Here is the boundary condition file for k

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

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

internalField   uniform 0;

boundaryField
{
    Side
    {
        type            fixedValue;
        value           uniform 0;
    }
    Inlet
    {
        type            cyclic;
    }
    Outlet
    {
        type            cyclic;
    }
}


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

do you think despite not solving for k using Smagorinsky model the internalField of k should be nonzero?

[RodriguezFatz]

Originally Posted by Dan1788

In fact I had made a mistake of assigning k to 0 using oneEqEdy model so that there wont be any production of turbulence. But later I changed the turbulence model to Smagorinsky that does not use a an k equation but is calculated. Here is the boundary condition file for k

Daniel, when I have a look into the equations of oneEqEddy model, I see the source term is

volScalarField G(GName(), 2.0*nuSgs_*magSqr(symm(gradU)));

Thus, if you don't initialize your case with nuSgs = 0 everywhere, also the oneEqEddy model should work. "k" isn't part of the source term. Both models should work.
When I started my case with Smagorinsky, I had an internal nuSgs field of "0" everywhere. That worked! Remember if you compare the pictures: I had a Reynolds number of 820000. Not 9000.

You can also have a look at nuSgs in paraview, just as you plot the velocity. If there is anything to plot, the model started to work.

You need to change your solver settings. Firstly change pressure solver, the one you use currently takes too many iterations.

 
p
    {
        solver           GAMG;
        tolerance        1e-6;
        relTol           0.1;

    maxIter         100;    

        smoother         DICGaussSeidel;
    
        nPreSweeps       0;
        nPostSweeps      1;
    nFinestSweeps    2;

        cacheAgglomeration true;

        nCellsInCoarsestLevel 50;
        agglomerator     faceAreaPair;
        mergeLevels      1;
    };

pFinal
    {
        solver           GAMG;
        tolerance        1e-6;
        relTol            1.0e-2;

    maxIter         100;    

        smoother         DICGaussSeidel;

        nPreSweeps       0;
        nPostSweeps      1;
    nFinestSweeps    2;

        cacheAgglomeration true;

        nCellsInCoarsestLevel 50;
        agglomerator     faceAreaPair;
        mergeLevels      1;
    };

Then, you should switch off the non-orthogonal correction of PIMPLE, or is your mesh that bad? I guess it's good enough, but you can still post the "checkMesh" output if you are not sure.

Then, put the code I posted in your controlDict to track the turbulent energy. You can create a plotscript like this:

set title "Kinetic Energy - Pipe LES"
set ylabel 'Kinetic Energy'
set xlabel 'time step'
set key off
set grid

plot "< cat log | grep 'Expression AverageResolvedKineticEnergyVxVy :' | cut -d'=' -f2" title 'kEn' with dots

Just create a fresh file (in your case directory) and copy the code in it. It only works like that if your log-file is actually called "log". Otherwise you need to change the last line "log" to something else. In the terminal type "gnuplot logfile -persist" after some iterations. Tell me if it doesn't work.

Edit: One last thing: If you want to see turbulent structures in the picture, it is much better to plot a velocity that is perpendicular to the main flow. Because it's average is zero and you can see small structures much easier.

[Dan1788]

Thanks a ton for your insight and help Philipp. I input your code from post#3 in my controlDict and I got the following error:

[0] [1] [2] [3] 
[1] 
[1] --> FOAM FATAL ERROR: 
[1] Unknown function type swakExpression

Any library files that need to be included?

Btw here is my checkMesh output. The number of non-orthogonal faces were high which is why I have used a non-ortho correcter in fvSolution. Even then, please let me know your thoughts on this output (dont want to make a run on a bad mesh).

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                  
    Inlet               3036     3081     ok (non-closed singly connected)  
    Outlet              3036     3081     ok (non-closed singly connected)  
    Side                21912    22000    ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.0075 -0.0075 0) (0.0075 0.0075 0.075)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-2.93937e-17 -3.05889e-17 -3.25725e-16) OK.
    Max cell openness = 1.0304e-15 OK.
    Max aspect ratio = 21.0616 OK.
    Minimum face area = 3.01205e-09. Maximum face area = 1.64985e-07.  Face area magnitudes OK.
    Min volume = 1.7177e-12. Max volume = 4.96941e-11.  Total volume = 1.32423e-05.  Cell volumes OK.
    Mesh non-orthogonality Max: 180 average: 6.23511
 ***Number of non-orthogonality errors: 6072.
  <<Writing 6072 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
    Max skewness = 1.10262 OK.
    Coupled point location match (average 1.38778e-17) OK.

Failed 1 mesh checks.

End

[Dan1788]

I tried including the following library files in my control dict:

libs (
        "libOpenFOAM.so"
        "libsimpleSwakFunctionObjects.so"
        "libswakFunctionObjects.so"
        "libgroovyBC.so"
    );

Doing this I get the following error:

[1]     From function getCellZoneID(const fvMesh &mesh,const word &name)
[1]     in file CellZoneValueExpressionDriver.C at line 63.
[1] 
FOAM parallel run exiting
[1] 
[2] 
[2] 
[2] --> FOAM FATAL ERROR: 
[2] The cellZone FLUID was not found in

[alexeym]

Hi,

I guess this

functions
(
    AverageResolvedKineticEnergyVxVy
    {
        type swakExpression;
        valueType cellZone;
        zoneName FLUID;
        accumulations (
            average
        );
        expression "U.x*U.x + U.y*U.y";
        verbose true;
    }
);

should be changed to

functions
(
    AverageResolvedKineticEnergyVxVy
    {
        type swakExpression;
        valueType internalField;
        accumulations (
            average
        );
        expression "U.x*U.x + U.y*U.y";
        verbose true;
    }
);

Or you can create cellZone FLUID with topoSet (which will contain all cells of your mesh).

About your question

Why does the box grid in the boxTurb utility have to be in powers of 2^k?

Cause the utility uses FFT for initialization of the field.

[rolloblues]

Mesh non-orthogonality Max: 180 average: 6.23511
 ***Number of non-orthogonality errors: 6072.
  <<Writing 6072 non-orthogonal faces to set nonOrthoFaces

In my opinion a value of 180 for max non-orthogonality is way too much.
Try to improve it in order to keep the max orthogonality below 70-80

[RodriguezFatz]

Yes, indeed. For a simple pipe like you meshed it you should get something like:
Mesh non-orthogonality Max: 39.1121 average: 6.23583
Non-orthogonality check OK.

[rolloblues]

In the wall normal direction I have resolved the mesh with a y+ = 1. My axial mesh has a resolution of \Delta-z+ = 15. The mesh front view is attached herewith

I didn't noticed the picture before, I edited it a little to explain my point.
It seems like you used some O-Grid-based blocking techinque for your mesh ( is it ICEM by any chance? ):

Inside the red circles that's where I suspect your orthogonality goes.
Since you want a lot of nodes on the edge of the o-grid, which is good, you probably should "relax" the transition of the mesh along the diagonal by reducing the size of the internal square like the one I drew in pink.

Hope this helps!

[Dan1788]

Thus, if you don't initialize your case with nuSgs = 0 everywhere, also the oneEqEddy model should work. "k" isn't part of the source term. Both models should work.
When I started my case with Smagorinsky, I had an internal nuSgs field of "0" everywhere. That worked! Remember if you compare the pictures: I had a Reynolds number of 820000. Not 9000.

Hi Philipp,

So I tried my case with the Smagorinsky model. I initialized k and nuSgs with non-zero values in the internalField. Following is the BC file for k

internalField   uniform 1;

boundaryField
{
    Side
    {
        type            zeroGradient;
    }
    Inlet
    {
        type            cyclic;
    }
    Outlet
    {
        type            cyclic;
    }
}

I ran my simulation for 0.8 seconds (roughly 100 flow through times). However the problem is still the same. Turbulence is not reached.

I even tried your script for the kinetic energy. However unlike your case, in mine the KE is of the order of 1e-10 throughout the run, which means that that there is no production in the x and y direction.

I have attached my figures for k, nuSgs and U below. Dont know where the problem lies

I didn't noticed the picture before, I edited it a little to explain my point.
It seems like you used some O-Grid-based blocking techinque for your mesh ( is it ICEM by any chance? ):

Inside the red circles that's where I suspect your orthogonality goes.
Since you want a lot of nodes on the edge of the o-grid, which is good, you probably should "relax" the transition of the mesh along the diagonal by reducing the size of the internal square like the one I drew in pink.

Thanks for your suggestion Davide. I mesh my geometry in Pointwise. I did exercise the change you recommended but I still get the same error.

I even tried a test of running checkMesh for a 16 X 16 X 16 cube with cyclic patches and it gives me the same error of orthogonality. I wonder if there is an issue exporting periodic meshes from Pointwise to OF. I even posted it on the forum

http://www.cfd-online.com/Forums/ope...ity-error.html

please let me know if you have any information regarding the same. Thanks!

[RodriguezFatz]

As far as I know, mesh orthogonality shows the maximum angle between two cell-face and cell-center lines. Yours shows an angle of 180 degree, which means that something is really fucked up. Not just some regular bad cells.

Ok, you can try this: Export a complete regular mesh from your meshing tool. No special "periodic" boundary conditions. Just name the one side "inlet" and one side "outlet". You can change everything you need for periodic b.c. in the files in openFoam.
Your /constant/polyMesh/boundary file should be edited to something like this:

/*--------------------------------*- 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       polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

3
(
    PIPE_INLET
    {
        type            cyclic;
        nFaces          21285;
        startFace       19012158;
        neighbourPatch    PIPE_OUTLET;
    }
    PIPE_OUTLET
    {
        type            cyclic;
        nFaces          21285;
        startFace       19033443;
        neighbourPatch    PIPE_INLET;
    }
    PIPE_WALL
    {
        type            wall;
        nFaces          118404;
        startFace       19054728;
    }
)

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

So after importing the mesh, just change type of inlet and outlet to "cyclic" and also insert the "neighbourPatch XYZ".
This works for me (mesh comes from ICEM).

Edit: Can you upload the case?

[rolloblues]

Dan I'm no expert of Pointwise (I've read your other post) but I definitely agree with Philipp.
I suspect your problem may arise because of the mesh periodicity you defined.

For your case you don't need mesh periodicity but only b.c. Periodicity, which you can fully define in OF dictionaries

[Dan1788]

Ok, you can try this: Export a complete regular mesh from your meshing tool. No special "periodic" boundary conditions. Just name the one side "inlet" and one side "outlet". You can change everything you need for periodic b.c. in the files in openFoam.

The problem about the non-orthogonality was Pointwise specific. The way around the problem was having 'renumberMesh -overwrite' before checkMesh.

I have uploaded my case here. You can have a look.

[RodriguezFatz]

Daniel, at least for the one-Equation model: You try to do wall resolved LES, right? I am not an expert for that, but I don't think you can set k to zeroGradient at the walls, you need k=0.
For the Smagorinsky model I am not sure... is nuSgs boundary condition used at all? Anyway, is nuSgs zerogradient correct at low-Re walls?