CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)
-   -   SST SimpleFoam Convergence Problems (https://www.cfd-online.com/Forums/openfoam/92484-sst-simplefoam-convergence-problems.html)

camoesas September 15, 2011 11:10
SST SimpleFoam Convergence Problems
 
2 Attachment(s)
Hey Everybody,

I´am trying to simulate the flow over a backwardfacing step. (2D, SST, steadystate) But I don´t get any convergence. (I am looking at residuals and flow solution) I´ve imported the mesh from icem and as the geometry is quite simple I presume theres no mesh error. (6450 cells)
I am running lowRe with the yPlus values:
Patch 3 named BOTTOM y+ : min: 0.303561 max: 6.06216 average: 2.17811
Patch 4 named HOTWALL y+ : min: 0.0991544 max: 6.86093 average: 2.33885
Patch 5 named STEP y+ : min: 0.458733 max: 5.17977 average: 2.02599

Here are the fvSchemes and fvSolution. I have also attached the whole setup.
Thanks for Help!

Camoesas


fvSchemes:
Code:
ddtSchemes
{
    default        steadyState;
}

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

divSchemes
{
    default        none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,k)      Gauss upwind;
    div(phi,omega)  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:
Code:
solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance      1e-06;
        relTol          0.1;
    }

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

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

SIMPLE
{
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    p              0.3;
    U              0.7;
    k              0.7;
    epsilon        0.7;
    R              0.7;
    nuTilda        0.7;
}

FelixL September 15, 2011 12:29
Good evening, Camoesas,


I presume your mesh is fully orthogonal and the wall is appropriately resolved cell-size-ratio-wise?

There are two things I'd comment to:

1. You say youre running the case on LowRe (Whats the Reynolds number of your case anyway?) and the y+ values at your wall are quite okay for that. So the mesh seems to be ready for a simulation with fully resolved wall boundary layers.
But your BCs for nut and k still impose wall functions! (for omega this is okay as omegaWallFunction works for all y+ values) This is incorrect since the wall functions are only suitable for y+ > 30!

2. Tighten the tolerance for omega in fvSolution a bit. Set it to something like 1e-20, as you can see in your residuals plot the omega equation converges rather quickly and isn't being solved anymore after that. It's better for the solution process to keep the omega tolerance quite tight.


Greetings,
Felix.

blacksquirrel September 15, 2011 12:41
Hi Camoesas,

I simulate a backwardFacing Step, too and had problems with the convergence. Are you doing the Driver and Seegmiller backwardsFacing Step?

I coudn't open your attached file, but it's maybe a problem with my system. Therefore I don't know what Boundary Conditions you are using.

I have an acceptable solution for incompressible,steadyState and kOmegaSST. I'm using the solver simpleFOAM.

I used following BC:
Wall: U = fixed Value 0; p= zeroGradient; k= fixedValue 0; omega = fixedValue 1e-8, nu_t=calculated
Inlet: U = fixedValue U_in; p=zeroGradient
Outlet: U=zeroGradient, p=fixedValue;

Your fv_solution seems OK. I used in the fv_schemes everywhere Gausslinear corrected, but I think your schemes are working as well.

I hope this helps you a little bit.
Greetings
blacksquirrel

camoesas September 19, 2011 04:35
HI Felix,
HI blacksquirrel,

1) I Have an Reynolds Number of 28000 (calculated with step high)
2) I am using the Paper of Vogel: "Combined Heat Transfer and Fluid Dynamic Measurements Downstream of a Backward-Facing Step" (1985)
The aim is to add heat transfer some day...

3) I have decreased tolerance a for omega a little bit ( 15 orders of magnitude:-) from 1e-5 to 1e-20

4) Now I am using the following boundary conditions:

Inlet: U = fixed 11.3; p = zero G; k = fixed 7.6e-4; omega = fixed 0.55; nut calculated 0;
Outlet: U = zero Gr; p = 0; k = fixedValue 7.6e-4; omega = fixed 0.55; nut calculated 0;
Walls: U = 0; p = zero Gr; k = fixedValue 7.6e-4; omega = fixed 0.55; nut calculated 0;

5) What boundary conditions do I have to set for k omega and nut for a free slip boundary condition? Or where can I find a Answer to this question?

I still get weird solutions for my problem. I think the mesh cant be a Problem, because its simple and checkMesh says its ok. Is there any other source of stupid errors I could have made?

thanks Camoesas

blacksquirrel September 19, 2011 05:19
Hello Camoesas,

Sorry, I can't answer your question no. 5.
But at the wall the boundary Condition of k should be 0. Because k is calculated from U', which is 0 at the Wall.

greetings
Black Squirrel

FelixL September 19, 2011 07:35
Good afternoon, Camoesas,


your wall boundary conditions are still wrong. Like the black squirrel said, k has to be set to 0 at the wall (for boundary layer resolving meshes!). Also, omega needs to be set to a value according to your first cell layer height - see http://turbmodels.larc.nasa.gov/sst.html for the formula. I recommend using omegaWallFunction (like I said in the second post!), because it does the calculation automatically.

Free slip condition means no slip? Either use the noslip boundary condition or symmetry. These are basically the same.


Greetings,
Felix

camoesas September 19, 2011 12:02
HI Everybody!

I have choosen my boundary inlet condition computed by formulas from this page: http://www.cfd-online.com/Wiki/Turbu...ary_conditions
At the wall U, k, omega have to be zero. But why is nut zero everywhere?
Am I right that there is a difference between nut and nutilda?

Now I have choosen my boundary conditions as follows: Now they must be right?! :eek:
omega:
Code:
boundaryField
{
    INLET
    {
        type            fixedValue;
        value          uniform 0.55;
    }
    OUTLET
    {
        type            zeroGradient;
    }
    TOP
    {
        type            symmetryPlane;
    }
    BOTTOM
    {
        type            omegaWallFunction;
        Cmu            0.09;
        kappa          0.41;
        E              9.8;
        beta1          0.075;
        value          uniform 0.55;
    }
    HOTWALL
    {
        type            omegaWallFunction;
        Cmu            0.09;
        kappa          0.41;
        E              9.8;
        beta1          0.075;
        value          uniform 0.55;
    }
    STEP
    {
        type            omegaWallFunction;
        Cmu            0.09;
        kappa          0.41;
        E              9.8;
        beta1          0.075;
        value          uniform 0.55;
    }
    EMPTY
    {
        type            empty;
    }
}
nut:
Code:
boundaryField
{
    INLET
    {
        type            calculated;
        value          uniform 0;
    }
    OUTLET
    {
        type            calculated;
        value          uniform 0;
    }
    TOP
    {
        type            symmetryPlane;
    }
    BOTTOM
    {
        type            calculated;
        value          uniform 0;
    }
    HOTWALL
    {
        type            calculated;
        value          uniform 0;
    }
    STEP
    {
        type            calculated;
        value          uniform 0;
    }
    EMPTY
    {
        type            empty;
    }
}
k:
Code:
boundaryField
{
    INLET
    {
        type            fixedValue;
        value          uniform 0.000766;
    }
    OUTLET
    {
        type            fixedValue;
        value          uniform 0;
    }
    TOP
    {
        type            symmetryPlane;
    }
    BOTTOM
    {
        type            fixedValue;
        value          uniform 0;
    }
    HOTWALL
    {
        type            fixedValue;
        value          uniform 0;
    }
    STEP
    {
        type            fixedValue;
        value          uniform 0;
    }
    EMPTY
    {
        type            empty;
    }
}
U:
Code:
boundaryField
{
    INLET
    {
        type            fixedValue;
        value          uniform (11.3 0 0);
    }

    OUTLET
    {
        type            zeroGradient;
    }

    TOP
    {
        type            symmetryPlane;
    }

    BOTTOM
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }
   
    HOTWALL
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }
   
    STEP
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }

    EMPTY
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }
}
p:
Code:
boundaryField
{
    INLET
    {
        type            zeroGradient;
    }

    OUTLET
    {
        type            fixedValue;
        value          uniform 0;
    }

    TOP
    {
        type            symmetryPlane;
    }

    BOTTOM
    {
        type            zeroGradient;
    }
   
    HOTWALL
    {
    type        zeroGradient;
    }
   
    STEP
    {
    type        zeroGradient;
    }

    EMPTY
    {
        type            empty;
    }
}
thanks camoesas

FelixL September 19, 2011 13:38
At the OUTLET, k should be set to zeroGradient. The rest looks good, though!

Greetings,
Felix

camoesas September 20, 2011 03:12
Hey Felix,

thanks for revising once again my boundary conditions! I will test them and tell results.

Camoesas

camoesas September 20, 2011 08:35
I still get nonsense solutions :confused:
But I have tried my mesh with the setup of the pitzDaily tutorial and it works fine. So I can be sure the mesh´s ok. Now trying to change the turbulence modell from keps to SST ...

camoesas September 20, 2011 10:12
5 Attachment(s)
HI Everybody,

is it right to set:
nut: calculated 0 everywhere?

How about this explanation: The solution is not converging because the problem is transient. Im getting this impression by looking on my solution values (pictures)

Can I tell from examining the residuals? next post

Camoesas

camoesas September 20, 2011 10:14
1 Attachment(s)
here the residuals

blacksquirrel September 20, 2011 10:31
Hi Camoesas,

I set nut on a fixed value at the Inlet and on zero Gradient at the outlet. And I estimated the value of nut via nu/nut= 10~100.

Here are my Boundary conditions:
nut
Code:
dimensions      [0 2 -1 0 0 0 0];

internalField  uniform 1.511e-4;

boundaryField
{
    Inlet
    {
        type            fixedValue;
        value          uniform 1.511e-4 ;
    }
    Outlet
    {
        type            zeroGradient;
    }
    Wall
    {
        type            calculated;
    }
    empty
    {
        type            empty;
    }
}
k
Code:
dimensions      [0 2 -2 0 0 0 0];

internalField  uniform 9.526;

boundaryField
{
    Inlet
    {
        type            fixedValue;
        value          uniform 9.526;
    }
    Outlet
    {
        type            zeroGradient;
    }
    Wall
    {
      type              fixedValue;
        value          uniform 0;
    }
    empty
    {
        type            empty;
    }
}
omega
Code:
dimensions      [0 0 -1 0 0 0 0];

internalField  uniform 63044.616;

boundaryField
{
    Inlet
    {
        type            fixedValue;
        value          uniform 63044.616;
    }
    Outlet
    {
        type            zeroGradient;
    }
    Wall
    {
      type              fixedValue;
      value        uniform 1e-8;
    }
    empty
    {
        type            empty;
    }
}
Greetings
Black Squirrel

camoesas September 20, 2011 11:55
Hey Blacksquirrel,

I ´ve tested your BC, first I copied the types and then i copied them 1:1 But no converged solution.

Why are youre omega values so high? I have computed my values with formulas from this page: http://www.cfd-online.com/Wiki/Turbu...ary_conditions
But mine are five orders of magnitude smaller. Whats your value for turbulence length scale? Why are you setting omega 1e-8 at the wall?
I thought nut is supposed to be greater than nu?! Do you have a reference for your estimation?

Thanks Camoesas

FelixL September 20, 2011 13:49
Hello, guys,


black squirrel:

Setting omega to 1e-8 at the wall is clearly wrong. Omega approaches infinity when y+ -> 0, so it has to be some finite, but pretty high value. See e.g. the book by WILCOX for more detailled information.

camoesas:

Please tell us something about the inlet turbulence characteristics of your case. I assume the inflow has a low turbulence intensity? From your U and k BCs I can estimate a TI of 0.2%, which is pretty low. This results in a quite low value for omega using the equations you referenced to.
My experience is that the SST model has difficulties with low values for omega and e.g. the paper "Effective Inflow Conditions for Turbulence Models in Aerodynamic Calculations" of SPALART et.al. seems to confirm this.
Usually my reasoning is: very low inflow turbulence intensity requires a high dissipation (omega), because then the turbulent viscosity is of the order of the molecular viscosity. Low k and Low omega might lead to very, very high turbulent viscosity values which clearly shouldn't be true for a low-turbulence inflow.

Setting nut to calculated everywhere is fine, since nut = f(k,omega). k and omega are known everywhere after solving the transport equations, thus nut can easily be calculated.

EDIT: I just had a look at your contour plots. Are you sure your domain is extended sufficiently in the streamwise direction? I'd double or even quadruple the distance, just to be sure. Shouldn't be too costly to do that.

Greetings,
Felix.

camoesas September 21, 2011 05:00
1 Attachment(s)
HI Felix,

Yea youre right! Turbulent intensity at the Inlet is very low: 0.2% In my opinion really low, but the paper tells so, but remains silent about turbulent length scale. So I had to guess: 1.38mm (calculated from the pitzdaily tutorial)
Now I doubled the outlet.
According to the motorbike tutorial (which is the only known tutorial for simpleFoam with SST) I changed in fvSolution to linear solvers:

p: PCG --> GAMG
U,k,omega: PBiCG --> smooth Solver

And now i get this fine converged solution! Does this mean the other settings are wrong? Whats linear solver? Do I have loss in accuracy like first order?

camoesas September 21, 2011 05:04
1 Attachment(s)
and here the velocity for the fully resolved wall boundary layer

camoesas September 21, 2011 05:23
1 Attachment(s)
HI,

I´ve also made a simulation with wallFunctions. And this simulation is just converging fine. I am using the following boundary conditions:

k:
inlet: value 7.66e-4
outlet: zeroGradient
walls: kqRWallFunction

omega:
inlet: value 0.55;
outlet: zeroGradient;
walls: omegaWallFunction

nut:
inlet: calculated 0;
outlet: calculated 0;
walls: nutkWallFunction;

The fvSchemes and fvSolution are still the same:
Code:
ddtSchemes
{
    default        steadyState;
}

gradSchemes
{
    default        Gauss linear;
}

divSchemes
{
    default        none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,k)      Gauss upwind;
    div(phi,omega)  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;
}
Code:
solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance      1e-06;
        relTol          0.01;
    }

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

    omega
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance      1e-20;
        relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;

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

relaxationFactors
{
    p              0.3;
    U              0.7;
    k              0.7;
    omega          0.7;
}
The solution is converging fine but then it starts oszillating. Can I accept this solution?
As I am using the same fvSchemes and fvSolution settings for the lowRe and for the highRe simulation, I conclude that for the lowRe I still have errors in the boundary conditions!? Because one is converging and one not.


Edit:
Yplus Values:
Patch 3 named HOTWALL y+ : min: 6.18212 max: 56.0136 average: 39.387
Patch 5 named BOTTOM y+ : min: 1.70152 max: 30.41 average: 12.862
Patch 6 named STEP y+ : min: 6.34481 max: 24.2188 average: 16.7388

Here are my y+ Values they are still to small, aren't they? But I am hesitating to coarse the grid any more in my opinion its already to coarse

camoesas September 21, 2011 05:38
2 Attachment(s)
Felix,

youre experience is also right here. I have increased omega by 5 and 7 orders of magnitude and convergence is getting better but its still not satisfying. See pictures.

But I am hesitating to fake the boundary conditions.

FelixL September 21, 2011 12:50
Good evening, camoesas,


I don't understand your post #16. Did you extend the grid in flow direction AND change the linear solvers?
The linear solvers are only responsible for solving the discretized system of equations (see CFD text books for that!). They shouldn't affect the convergence behaviour of the solution. But extending the domain seems to be the right way to me. It should be possible to use PBiCG/PCG for the linear system solving, this doesn't affect the order of accuracy of your solution (despite of the convergence error, of course, but this is easily eliminated), though, so you can also keep GAMG and smoothSolver.

Wall functions shall only be used when (at least!) the average y+ is > 30. But your results in #16 seem to be quite reasonable so why use wall functions anyway?
My recommendation: use the inlet values for k and omega you initially specified and keep on extending the domain until there's no reasonable effect on the solution anymore.
Oh, and if you don't want your solution to degrade to 1st order, don't forget to use something more accurate than upwind for the divergence interpolation of k and omega!


Greetings,
Felix

camoesas September 23, 2011 10:26
4 Attachment(s)
Hey Felix,

Rereading my post you´re right its not clear enough. I extended the domain with no difference in the solution and convergence. And then I have changed PCG > GAMG and PCiCG > smooth solver. And I have got a better solution. But I could not reproduce this result. Must have been something else.

Now I ´ve build Up my case completely new (I have lost some files rm -rf * :eek:)
And now its working... Unfortunately I dont know what I have changed.

As you suggested I changed from 1. to second Order. In an first step just for U and then for k and omega. Of course convergence is getting worse and it takes longer to converge. The solution stops telling:
Code:
SIMPLE solution converged in 1176 iterations
Do you think I should run it a little bit longer. The residuals are still quite big.
Do I need second order for k and omega? The tutorial just uses Upwind.

I have tried a setup with wallFunctions because I wanted gather some experience. In my final geometry I wont be able to fully resolve boundary Layer.

I have also added my latest setup with that I achieve something like convergence.

Thanks for your help

Camoesas

FelixL September 24, 2011 03:09
Good morning, camoesas,


well, well - the mysteries of computational fluid dynamics :D I'm glad it works now.

Yeah it's quite common that the residuals convergence rate decreases when increasing the accuracy of the interpolation schemes. Sometimes they even seem to stall. This is not that big an issue, though, because the residuals are just a hint for convergence - they don't neccessarily have to vanish, it's a question of their definition.
I'd monitor other integral parameters to check for convergence - in your case e.g. the pressure drop along the computational domain. If this is constant you have your convergence.

Regarding upwinding k and omega: "The chain is only as strong as its weakest link". So if you have a first order scheme somewhere, your overall accuracy degrades. Nevertheless, if you are unsure you can do a grid convergence check and see how it influences your solution.


Greetings
Felix

camoesas September 26, 2011 08:10
HI Felix,

How could I monitor pressure drop along the domain? Do I have to write an expression or is there already a usable function.

Grid convergence is one of my next topics. Thanks for the reminder ;)

camoesas

FelixL September 27, 2011 16:27
Using the drag coefficient along the channel's walls should do the trick ;) This should be constant in the end, too (for a steady state problem!).

Greetings,
Felix


All times are GMT -4. The time now is 21:15.