[Peter85]

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:
http://www.cfd-online.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 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:
http://www.cfd-online.com/Forums/ope...ess-model.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 SST-Model?

[Peter85]

Here more diagrams...
nut+.jpg

TurbulentShearStress.jpg

U+.jpg

Epsilon+_gnu.jpg

k+_gnu.jpg

[Peter85]

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]

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?

[FelixL]

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.cfd-online.com/Forums/ope...nce-model.html

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

Greetings,
Felix.

[Peter85]

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]

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

[FelixL]

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]

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:

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

    wall       
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

p:

    inlet           
    {
        type            zeroGradient;
    }

    wall       
    {
        type            zeroGradient;
    }

k:

    inlet
    {
        type            directMapped;
        value           uniform 0.015;
        setAverage      false;
        average        0.015;
    }

    wall
    {
        type            fixedValue;
        value           uniform 1e-11;
    }

omega:

    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.

[Peter85]

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

[FelixL]

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.:

    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.

[Peter85]

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