[student666]

Hi all,

I'm running a steady state simulation, laminar with simpleFoam.
I have a rotating frame defined into fvOptions dict file.

I have the following dict for fvSolution:

/*--------------------------------*- 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          GAMG;
        tolerance       1e-08;
        relTol          0.001;//era 0.05
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }   
    
    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        tolerance       1e-07;
        relTol          0.01;// era 0.1
    }

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

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

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

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

SIMPLE
{
    nNonOrthogonalCorrectors 1;

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



relaxationFactors
{
    fields
    {
        p               0.1;//0.3;
    }
    equations
    {
        U               0.3;//0.7;
        k               0.7;
        epsilon         0.7;
        R               0.7;
        nuTilda         0.7;
    }
}

and for fvScheme:

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default        none; // Gauss linear;
     p     Gauss;
     U    Gauss;

}

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;
//    interpolate(U)  linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
}

I monitor the mass flow as well.


I cant't understand the strange behavior for pression, and I can't understand if this is a problem of meshing (4milion cells), numerical scheme or linear equation solution.


If I run checkMesh I have max non-orthogonality about 55 average 6.5; with 25 faces with high skew (9). Anyway I was able to obtain the same mesh with 1 skew cell, but the problem still there.


Any suggestion how to improve results? please help.


https://www.dropbox.com/s/vwg1vmb2gpy53ld/running1.png



Thanks a lot

[RodriguezFatz]

Can you give a short description (with picture) of your setup, that shows inlet and outlet / walls and all that stuff?
Is your case really laminar?

[student666]

Hi Mr. Rodriguez,

here's the picture of my domain:

https://www.dropbox.com/s/0hq9kndjzo9fkfl/domain.png

green is inlet and blue is outlet, while red is the rotating frame.

Indeed it is turbulence at all, but normally I start with 1°st order scheme (as suggested in many threads on cfd-on line) and no turbulence model activated.
Further I'm surely going to turn on turbulence and raise the order of the div scheme.

Further this is an improvement on a previous model (I modified the shape of yellow frame) and I started by reusing the same set up for fvSolution and for fvScheme. that's the reason of my confusion, because previously I had a good convergence: residual below 1e-4

[RodriguezFatz]

But if you refine the mesh and keep the turbulence model switched off, you are able to resolve more and more of the turbulence (DNS-like). Without the turbulence model you do not have a steady state solution, but the solver tries to get one.
I would rather not run a simulation without turbulence model. When it doesn't converge you can never be sure if it is due to the missing turbulence model or due to anything else.

[student666]

thanks for the hints!

if you refine the mesh and keep the turbulence model switched off, you are able to resolve more and more of the turbulence (DNS-like)

this is another question I always have in mind.
I'm not able to add layers at the walls due to complexity of the geometry (I'm using SHM, no commercial mesher provided by my company), I always think to refine the mesh near the wall reducing the size of the cells, in order to (in some manner) resolve the boundary layer with turbulence model switched off. Shall it works anyway if I have an y+ calculated lower than 10? SHM works hard to give hexahedral dominant domain, so maybe I'm able to have a sort of "prismatic layers" at walls anyway...

Without the turbulence model you do not have a steady state solution, but the solver tries to get one.

Do you mean that solver, even if the ddt scheme is steady state, runs a time marching solution anyway if turbulence is set to "laminar"?
I red something about iterative solvers that have to march in time anyway to resolve matrix (Peric) in order to reduce the error of the calculation. Why, if I turn the turbulence on, it becomes steady state? Can you clarify me? I'm bit confused...

[RodriguezFatz]

To the first part: I am not sure what you mean. You mean you set y+ = 1 in order to switch off the turbulence model? No, that's not the reason. You set y+ = 1 with turbulence model on (!) to be able to capture all kind different flows. The boundary sheath is pretty complex, also for RANS-model flows. It's just that most boundary sheaths actually look the same. That's why you can skip that part and use a wall function to emulate that layer. But if you have a detached / stall flow this wall function of course is not the correct choice.

To the second part: I think what you mean is a "pseudo-transient" solver. SimpleFoam is a steady-state solver that does not use any time at all. Changing the dt-scheme doesn't effect it, as far as I know. Also chaning the time step size will only effect the name of the directories where the data is stored, not the numerics. What I ment was, that since your solver is a steady-state solver it is designed to get a steady state solution. If you however solve a case that actually is time-dependent the solver will not be able to lower the residuals.

[student666]

first part:
what I mean is that I have big confusion, and know nothing about LES and DNS, and poorly of RANS so

... maybe a silly question...

My question is:
as RANS use to model the behavior of the flow by using correlation to model the reynold's stress and resolve the log-layer and buffer layer ecc.., I say, is it totaly wrong try to reduce cell size near walls (to the exterme to the order of the boundary layer) and let the computation go on to calculate the velocity gradient near the wall, even if laminar is set into turbulence dict?
probably yes (totally wrong) because I will never be able to model such small elements as calculation time blows up, and so I'm forced to use correlation anyway, but in this case I have to check the proper values for y+ for the turbulence model (RANS) I'm using?

told you I'm confused....

second part: ok, we meant the same thing.

[RodriguezFatz]

Originally Posted by student666

probably yes (totally wrong) because I will never be able to model such small elements as calculation time go blows up, and so I'm forced to use correlation anyway, but in this case I have to check the proper values for y+ for the turbulence model (RANS) I'm using?

Totally right, if you don't use any turbulence model, but just a small mesh that is called a DNS. If you model only a part of the turbulence, it would be an LES. But even LES (without wall model) needs ridiculous amount of cells in the boundary, that makes it pretty much unusable for any technical relevant simulation. Also, all of them (time-)resolve the flow, so you can not do that with a steady state solver.
And: RANS models the turbulence in the complete domain, not only in the boundary. It sounds like you mix this up...

[student666]

Thanks!

Regards.