
[Sponsors] 
October 29, 2010, 03:38 
Problem with SSTModel  strange behaviour

#1 
New Member
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 13 
Hi@all!
I did some testing on the SSTModel 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 nutwall functions, referring to this topic: http://www.cfdonline.com/Forums/ope...komegasst.html 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 SSTModell simply delivers a completely wrong profile for high resolution meshes. I did some further investigations on the model. I compared it to some lowremodels implemented in OpenFoam and to the EARSMmodel posted here with my y+1 mesh: http://www.cfdonline.com/Forums/ope...essmodel.html 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 SSTModel? 

October 29, 2010, 03:41 

#2 
New Member
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 13 


October 29, 2010, 03:41 

#3 
New Member
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 13 
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 

November 2, 2010, 12:29 

#4 
New Member
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 13 
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 lowRe models, it is unchanged, but for the SSTmodels, there is a huge difference. The same SSTEARSM model was used with a y+1 circular channel and obtained good results, but both SSTmodels are not working for a rectangular channel with y+1. Does anybody know the reason for this strange behaviour? Has anyone used the SSTmodel with a rectangular channel and obtained good results for a y+1 mesh or had the same problems? 

November 4, 2010, 16:14 

#5 
Senior Member
Felix L.
Join Date: Feb 2010
Location: Hamburg
Posts: 165
Rep Power: 15 
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: http://www.cfdonline.com/Forums/ope...ncemodel.html Hopefully it works for you, too  I'd be glad to see your results! Greetings, Felix. 

November 5, 2010, 02:17 

#6 
New Member
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 13 
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 

November 5, 2010, 03:35 

#7 
New Member
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 13 
Hi again!
I have modified my SSTmodel 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 

November 15, 2010, 06:35 

#8 
Senior Member
Felix L.
Join Date: Feb 2010
Location: Hamburg
Posts: 165
Rep Power: 15 
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. 

November 15, 2010, 09:24 

#9 
Senior Member
Felix L.
Join Date: Feb 2010
Location: Hamburg
Posts: 165
Rep Power: 15 
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); } Code:
inlet { type zeroGradient; } wall { type zeroGradient; } Code:
inlet { type directMapped; value uniform 0.015; setAverage false; average 0.015; } wall { type fixedValue; value uniform 1e11; } Code:
inlet { type directMapped; value uniform 1.064794275e2; setAverage false; average 1.064794275e2; } wall { type omegaWallFunction; value uniform 1.064794275e2; } 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. 

November 15, 2010, 15:04 

#10 
New Member
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 13 
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 kepsilon model and komega 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 komega to kepsilon. But towards the middle of the channel, it raises again, switching back to the komega 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 

November 16, 2010, 15:24 

#11 
Senior Member
Felix L.
Join Date: Feb 2010
Location: Hamburg
Posts: 165
Rep Power: 15 
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 1e06), the omega equation isn't being solved anymore after about 30 iterations: residuals_tol_1e06.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 1e10; relTol 0; } omega { solver PBiCG; preconditioner DILU; tolerance 1e10; relTol 0; } Now the residual history shows that the omega equation is being solved during the whole simulation: residuals_tol_1e10.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. 

November 18, 2010, 01:32 

#12 
New Member
Peter
Join Date: Aug 2010
Posts: 16
Rep Power: 13 
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 CFDsimulations before, so this problem was not easy to find for me. Thanks again! Peter 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Problem on EDC model for coal combustion  lei  FLUENT  4  September 3, 2015 09:39 
Superlinear speedup in OpenFOAM 13  msrinath80  OpenFOAM Running, Solving & CFD  18  March 3, 2015 05:36 
SST model  swe704  Main CFD Forum  0  February 5, 2010 08:36 
Understanding komega SST model source code  tmhonka  OpenFOAM Programming & Development  1  September 8, 2009 07:33 
strange behaviour of the oneequation SGS model  cfdmarkus  OpenFOAM Running, Solving & CFD  18  August 14, 2009 05:33 