# SST SimpleFoam Convergence Problems

 Register Blogs Members List Search Today's Posts Mark Forums Read

September 15, 2011, 11:10
SST SimpleFoam Convergence Problems
#1
Senior Member

Join Date: Mar 2009
Posts: 138
Rep Power: 9
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         Gauss linear;
}

divSchemes
{
default         none;
div(phi,k)      Gauss upwind;
div(phi,omega)  Gauss upwind;
}

laplacianSchemes
{
default         Gauss linear corrected;
}

interpolationSchemes
{
default         linear;
}

{
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;
}```
Attached Images
 Residuals.jpg (37.0 KB, 152 views)
Attached Files
 step.zip (30.0 KB, 36 views)

 September 15, 2011, 12:29 #2 Senior Member   Felix L. Join Date: Feb 2010 Location: Hamburg Posts: 165 Rep Power: 10 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. Tauseef likes this.

 September 15, 2011, 12:41 #3 Member   Join Date: Jun 2011 Posts: 53 Rep Power: 7 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 Last edited by blacksquirrel; September 15, 2011 at 13:08.

 September 19, 2011, 04:35 #4 Senior Member   Join Date: Mar 2009 Posts: 138 Rep Power: 9 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

 September 19, 2011, 05:19 #5 Member   Join Date: Jun 2011 Posts: 53 Rep Power: 7 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

 September 19, 2011, 07:35 #6 Senior Member   Felix L. Join Date: Feb 2010 Location: Hamburg Posts: 165 Rep Power: 10 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

 September 19, 2011, 12:02 #7 Senior Member   Join Date: Mar 2009 Posts: 138 Rep Power: 9 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?! 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

 September 19, 2011, 13:38 #8 Senior Member   Felix L. Join Date: Feb 2010 Location: Hamburg Posts: 165 Rep Power: 10 At the OUTLET, k should be set to zeroGradient. The rest looks good, though! Greetings, Felix

 September 20, 2011, 03:12 #9 Senior Member   Join Date: Mar 2009 Posts: 138 Rep Power: 9 Hey Felix, thanks for revising once again my boundary conditions! I will test them and tell results. Camoesas

 September 20, 2011, 08:35 #10 Senior Member   Join Date: Mar 2009 Posts: 138 Rep Power: 9 I still get nonsense solutions 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 ...

September 20, 2011, 10:12
#11
Senior Member

Join Date: Mar 2009
Posts: 138
Rep Power: 9
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
Attached Images
 1e3.jpg (11.8 KB, 54 views) 2e3.jpg (11.4 KB, 40 views) 3e3.jpg (11.8 KB, 40 views) 4e3.jpg (11.7 KB, 34 views) 5e3.jpg (11.9 KB, 44 views)

September 20, 2011, 10:14
#12
Senior Member

Join Date: Mar 2009
Posts: 138
Rep Power: 9
here the residuals
Attached Images
 Residuals1.png (66.2 KB, 76 views)

 September 20, 2011, 10:31 #13 Member   Join Date: Jun 2011 Posts: 53 Rep Power: 7 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

 September 20, 2011, 11:55 #14 Senior Member   Join Date: Mar 2009 Posts: 138 Rep Power: 9 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

 September 20, 2011, 13:49 #15 Senior Member   Felix L. Join Date: Feb 2010 Location: Hamburg Posts: 165 Rep Power: 10 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.

September 21, 2011, 05:00
#16
Senior Member

Join Date: Mar 2009
Posts: 138
Rep Power: 9
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?
Attached Images
 Residuals_linearSolver.png (38.6 KB, 101 views)

September 21, 2011, 05:04
#17
Senior Member

Join Date: Mar 2009
Posts: 138
Rep Power: 9
and here the velocity for the fully resolved wall boundary layer
Attached Images
 longoutlet.jpg (13.2 KB, 68 views)

September 21, 2011, 05:23
#18
Senior Member

Join Date: Mar 2009
Posts: 138
Rep Power: 9
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
walls: kqRWallFunction

omega:
inlet: value 0.55;
walls: omegaWallFunction

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

The fvSchemes and fvSolution are still the same:
Code:
```ddtSchemes
{
}

{
default         Gauss linear;
}

divSchemes
{
default         none;
div(phi,k)      Gauss upwind;
div(phi,omega)  Gauss upwind;
}

laplacianSchemes
{
default         Gauss linear corrected;
}

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

{
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
Attached Images
 Residuals2.png (51.3 KB, 74 views)

September 21, 2011, 05:38
#19
Senior Member

Join Date: Mar 2009
Posts: 138
Rep Power: 9
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.
Attached Images
 Residuals_omega0.55e5.png (40.9 KB, 72 views) Residuals_omega0.55e7.png (40.1 KB, 63 views)

 September 21, 2011, 12:50 #20 Senior Member   Felix L. Join Date: Feb 2010 Location: Hamburg Posts: 165 Rep Power: 10 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

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post ziyan7 FLUENT 0 March 8, 2011 07:13 boddyouareboy OpenFOAM 0 January 25, 2011 04:49 A.Devesa OpenFOAM Running, Solving & CFD 0 November 9, 2010 05:43 CedricVH OpenFOAM Running, Solving & CFD 8 June 3, 2010 09:05 david aquilina CD-adapco 4 October 27, 2000 09:18

All times are GMT -4. The time now is 18:35.