CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Problem with SST-Model - strange behaviour

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

Like Tree1Likes
  • 1 Post By FelixL

Reply
 
LinkBack Thread Tools Display Modes
Old   October 29, 2010, 03:38
Default Problem with SST-Model - strange behaviour
  #1
New Member
 
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 5
Peter85 is on a distinguished road
Hi@all!

I did some testing on the SST-Model implemented in OpenFoam. I tried to simulate a simple rectangular channel in 2D with D_hyd = 30mm (equals 15mm height) at a Reynoldsnumber of 20000 (U=10 m/s). Therefore, I calculated the pressure loss using the correlations for a channel, so the flow was initiated by a pressure drop.
The simulation was periodic (using the "cyclic" bc with "fan" to get the pressure drop).
I did this simulation for several meshes with different distances of the closest cell to the wall (y+= 0.5 ; 1 ;2 ;15 and 30). I also tried all available nut-wall functions, referring to this topic:
Compressible kOmegaSST
For the other variables at the wall, I set K=0, omega=omegawallfunction, alphat = alphatwallfunction, p=zerogradient and U = 0.
I used schemes of second order accuracy for the calculations.

The results I got were quite amazing. I got the best results for y+2, but when the mesh was refined, the results got worse. The velocity profile looked more like a laminar profile, and the mean value was far below the expected 10 m/s.
Here the diagrams for y+1 and y+2:
u+_y+_y+1.jpg


u+_y+_y+2.jpg

as you can see, there is quite a difference between those profiles. As mentioned before, when the mesh was refined, the results got worse, the profile looked much more laminar and the mean value of the velocity got lower. In those diagrams you can see only one line, the reason is that all values of the different nut wall functions were equal. This is different when the mesh gets coarser, but this is not important here. The problem is that the SST-Modell simply delivers a completely wrong profile for high resolution meshes.

I did some further investigations on the model. I compared it to some low-re-models implemented in OpenFoam and to the EARSM-model posted here with my y+1 mesh:
New implemented algebraic Reynolds stress model

The results are here:
epsilon+.jpg
k+.jpg
U+_daten_eingelesen.jpg

Well...I just can attach 5 files, so I will post the other ones later.

Can someone explain the strange behaviour of the SST-Model?
Peter85 is offline   Reply With Quote

Old   October 29, 2010, 03:41
Default
  #2
New Member
 
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 5
Peter85 is on a distinguished road
Here more diagrams...
nut+.jpg

TurbulentShearStress.jpg

U+.jpg

Epsilon+_gnu.jpg

k+_gnu.jpg
Peter85 is offline   Reply With Quote

Old   October 29, 2010, 03:41
Default
  #3
New Member
 
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 5
Peter85 is on a distinguished road
and more...
TurbulentShearStress+_gnu.jpg

In some diagrams, there is a line calles "paper". This is valid not for the channel but for a flat plate and is just plotted for orientation
Peter85 is offline   Reply With Quote

Old   November 2, 2010, 13:29
Default
  #4
New Member
 
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 5
Peter85 is on a distinguished road
I did some more simulations on this channel. This time, the mesh was slightly changed, I used a mesh with y+2 distance of the closest cell to the wall.
The velocity profile is much better.
The main difference was observed in the turbulent viscosity:
nut+.jpg
For the low-Re models, it is unchanged, but for the SST-models, there is a huge difference.
The same SST-EARSM model was used with a y+1 circular channel and obtained good results, but both SST-models are not working for a rectangular channel with y+1.

Does anybody know the reason for this strange behaviour?
Has anyone used the SST-model with a rectangular channel and obtained good results for a y+1 mesh or had the same problems?
Peter85 is offline   Reply With Quote

Old   November 4, 2010, 17:14
Default
  #5
Senior Member
 
Felix L.
Join Date: Feb 2010
Location: Hamburg
Posts: 161
Rep Power: 8
FelixL is on a distinguished road
Hi, there,


I discovered that the current implementation of the kOmegaSST turbulence model in OpenFOAM leads to overpredicted values of the eddy viscosity when using fine meshes with near wall cell spacings of y+=1 and below. This is caused by a missing factor in the definition of nut.

It's difficult to distinguish the different cuves in your last attachement because some curves are of the same color, but it looks like your eddy viscosity profiles for the SST models seem to be overpredicted, too.

It would be great if you could try out the modification I proposed here: Wrong calculation of nut in the kOmegaSST turbulence model

Hopefully it works for you, too - I'd be glad to see your results!

Greetings,
Felix.
FelixL is offline   Reply With Quote

Old   November 5, 2010, 03:17
Default
  #6
New Member
 
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 5
Peter85 is on a distinguished road
Hi!
Thank you for your answer! The wrong calculation of nut seems to be one of the problems we observed. The other one is the strong grid dependence of nut. I am going to change the code as you suggested, I hope this works.
Sorry for the bad diagram. I have created another two, one with gnuplot and one with parafoam. The reason why I post both is the peak, that is different using both tools. It seems to be an interpolation problem in one of those tools.
By the way: nutquer is simply a dimensionless nut, means nut/nu
nutquer.jpg

nutquer_gnu.jpg
Peter85 is offline   Reply With Quote

Old   November 5, 2010, 04:35
Default
  #7
New Member
 
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 5
Peter85 is on a distinguished road
Hi again!
I have modified my SST-model as you suggested and simulated my case for a y+1 mesh. The results, compared to the unmodified model are here:
Epsilon+.jpg

k+.jpg

nutquer.jpg

U+.jpg

U+_complete.jpg

As you can see, the results are a little bit better. But the velocity profile looks still not like a turbulent velocity profile, more like a laminar one. But my Reynoldsnumber is 20000, so its turbulent.

Peter
Peter85 is offline   Reply With Quote

Old   November 15, 2010, 06:35
Default
  #8
Senior Member
 
Felix L.
Join Date: Feb 2010
Location: Hamburg
Posts: 161
Rep Power: 8
FelixL is on a distinguished road
Hi, Peter,


could you post your setup and boundary conditions plus some geometrical details, please, so I can reproduce that case and have a deeper look into it?

Thanks and Greetings,
Felix.
FelixL is offline   Reply With Quote

Old   November 15, 2010, 09:24
Default
  #9
Senior Member
 
Felix L.
Join Date: Feb 2010
Location: Hamburg
Posts: 161
Rep Power: 8
FelixL is on a distinguished road
Hey, Peter,


out of curiosity I just did a few calculations using the geometrical data und reynolds number you provided in the first post. I created two meshes with y+ at the wall of 2 and 1, respectively. The dimension of the domain in the axial direction was 0.15 meters, the front and back planes were set to empty to create a 2D solution.

I am not quite sure what you meant with using cyclic and fan BCs to get the pressure drop right. For the inlet BC I found it more suitable to use the directMapped boundary condition to map the velocity outlet field back to the inlet patch. With this BC you can specify an average value of U so you don't have to calculate and set a pressure drop over the domain.

Since you didn't provide any details on the inlet conditions for the turbulent quantities, I assumed a turbulence intensity of 1% and a turbulent length scale of 0.07 times the hydraulic diameter.

Inlet and wall boundary conditions for U, p, k and omega are listed here:

U:
Code:
    inlet           
    {
        type            directMapped;
        value           uniform (10 0 0);
        setAverage      true;
        average         (10 0 0);
    }

    wall       
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
p:
Code:
    inlet           
    {
        type            zeroGradient;
    }

    wall       
    {
        type            zeroGradient;
    }
k:
Code:
    inlet
    {
        type            directMapped;
        value           uniform 0.015;
        setAverage      false;
        average        0.015;
    }

    wall
    {
        type            fixedValue;
        value           uniform 1e-11;
    }
omega:
Code:
    inlet
    {
        type            directMapped;
        value           uniform 1.064794275e2;
        setAverage      false;
        average        1.064794275e2;
    }

     wall
    {
        type            omegaWallFunction;
        value        uniform 1.064794275e2;    
    }
The numerical schemes for all divergence terms were set to Gauss linear;.


Using these settings I wasn't able to reproduce the weird behavior you have observed. Actually the velocity profiles look very good for both meshes (y+ = 1 and y+ = 2) and resemble a turbulent velocity profile (unlike the "laminar" shape you observed at the finest mesh).
The profiles of the eddy viscosity don't match your results at all, especially in the central region. As you can see, when refining the mesh the peak value of nut increases, but the global shape of the profile doesn't change at all.

So my results are, that the kOmegaSST turbulence model is behaving as expected for meshes with grid spacings of the order of y+ = 1 near the walls.

I am pretty sure that the boundary conditions you specified were causing your trouble - it would be helpful if you could provide the boundary definitions of your case. Maybe we can find the error source leading to the strange behaviour you observed.

Greetings,
Felix.
Attached Images
File Type: png uPlus.png (20.9 KB, 41 views)
File Type: png nut.png (24.5 KB, 40 views)
FelixL is offline   Reply With Quote

Old   November 15, 2010, 15:04
Default
  #10
New Member
 
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 5
Peter85 is on a distinguished road
Hi Felix,

I can give you the case files of my simulation. The turbulence intensity was set to 5% in my case. I also did some simulations with your sqrt(2)-modification to the models, the results were a little bit better, but nu_t is still wrong. I tried to use "calculated" instead of the LowReWallFunction for nut, the results were the same.

I also did some simulations for a channel without the b.c. "cyclic", and using a mesh created with blockmesh and not with icem. Some results were good, some not. First, I did the simulation for 10 m/s, starting with internal field velocity = 0 m/s. The results were quite good, but I observed some oscillation in the residuals. To get better results, I changed the internal field at the beginning to 5 m/s. At the end of the simulation, I was surprised to see the same problem as in the periodic channel, a drop of nu_t towards the symmetry axis that I cannot explain.

This behavior is quite surprising and should not happen. Further investigation of the first channel showed a wrong behavior of the blend function switching between k-epsilon model and k-omega model, called "F1" in the code. It is 1 at the wall as it should be. Further away from the wall it drops, switching from k-omega to k-epsilon. But towards the middle of the channel, it raises again, switching back to the k-omega model, thats wrong. I don't know if that is just a result or the reason for the wrong behavior.

I tried to attach the complete cases, but unfortunately the board does not allow the size for this case (about 4mb).

Here are the "0" files:
0_periodic_channel.tar.gz

0_long_channel.tar.gz

I can also give you the system and constant folder, including the blockmeshdict for the second case, so you can reproduce my case:
system.tar.gzconstant.tar.gz

Finally, here is the nu_t diagram for the second case:
nutquer.jpg

as you can see, the results for the setting with 5 m/s for the internal field at the beginning are wrong. But both cases were completely converged.

Thank you !

Peter
Peter85 is offline   Reply With Quote

Old   November 16, 2010, 15:24
Default
  #11
Senior Member
 
Felix L.
Join Date: Feb 2010
Location: Hamburg
Posts: 161
Rep Power: 8
FelixL is on a distinguished road
Good Evening, Peter,


thanks for posting the cases. At first I want to say something about the weird behavior of the blending function F1 you observed for the first channel case: I am pretty confident, this isn't the source of the problem but another "side effect" of the erroneous behavior of your simulation. The function F1 usually works flawless for such simple cases - if you have a look at the implementation, you can see that it depends (among others) on the values of omega and k. So if your simulation (for whatever reasons) leads to too small values of k and too high values of omega near the symmetry plane of the channel, F1 switches back to zero. But this isn't a problem of the blending function itself, the source has to be somewhere else.


Anyways - I didn't have the time today to do full simulations of the second case you kindly provided to me, but I think this isn't really neccessary, as I am pretty sure I found the reason for all your headache.

I was able to reproduce the results you had with having a velocity of 5 m/s as a starting field for U: The eddy viscosity was way too low in the central region of the channel. (By the way: why 5 m/s as an initial guess and not 10 m/s? It's much closer to the final solution.) Furthermore the start of the simulation was pretty sensitive and unstable because you were fixing the pressure at the outlet to a value of 100000 mē/sē while the initial value of the internal field was set to zero - this huge pressure difference is giving the solver a really hard time and I assumed this was a mistake of yours so I felt free to change the outlet value to zero, improving the convergence behaviour of the simulation significantly.

However, this was not the source of your problems. When looking at the residual history for the first 1000 iterations, you see that with your tolerance settings for the linear solvers (tolerance 1e-06), the omega equation isn't being solved anymore after about 30 iterations:

residuals_tol_1e-06.png

This is crucial because the omega field depends on the velocity and turbulent kinetic energy fields and vice versa! This eventually leads to the wrong results for the eddy viscosity you observed and has of course a direct effect on the velocity distribution. Strictly speaking: your simulation didn't converge, although it seemed to have, because your tolerances for omega were too high.

I fixed this by lowering the tolerances for k and omega, viz.:
Code:
    k
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-10;
        relTol          0;
    }

    omega
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-10;
        relTol          0;
    }
Everthing else was left unchanged.

Now the residual history shows that the omega equation is being solved during the whole simulation:

residuals_tol_1e-10.png

The results obtained with these changes were significantly improved - the maximum value of nut was in the symmetry plane and the velocity profile resembled the typical turbulent velocity profile.

I am pretty sure this will also help with your first case using the cyclic BCs.

Greetings,
Felix.
pjohannes183 likes this.
FelixL is offline   Reply With Quote

Old   November 18, 2010, 01:32
Default
  #12
New Member
 
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 5
Peter85 is on a distinguished road
Hi Felix!

Thank you very much, you solved the problem! I just did the simulation for the cyclic channel with your changes, and the results were correct. I think this will work for all other cases, too.
I really appreciate your help...I tried to solve this problem for more than a month!
Well, I have to say, I've never done CFD-simulations before, so this problem was not easy to find for me.

Thanks again!

Peter
Peter85 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem on EDC model for coal combustion lei FLUENT 2 March 18, 2011 02:06
SST model swe704 Main CFD Forum 0 February 5, 2010 08:36
Understanding k-omega SST model source code tmhonka OpenFOAM Programming & Development 1 September 8, 2009 07:33
Superlinear speedup in OpenFOAM 13 msrinath80 OpenFOAM Running, Solving & CFD 17 August 22, 2009 03:59
strange behaviour of the one-equation SGS model cfdmarkus OpenFOAM Running, Solving & CFD 18 August 14, 2009 05:33


All times are GMT -4. The time now is 05:50.