Serious bug in LES interface
 

Hello guys,

I want to report a serious bug in the LES interface. I had build my own SGS model in OpenFOAM and with the release of OpenFOAM 1.6 I decided to use the new OpenFOAM. After spending much work to the adaption of my SGS Model to the new interface of OpenFOAM 1.6, I was wondering about differences between the Solutions of OpenFOAM 1.5 and OpenFOAM 1.6. So I assumed a mistake of me while reprogramming my model and I spent many hours to find the mistake. At the end I decided to go back to a simpler model and check the differences between version 1.5 and 1.6 with the channel395 tutorial and a simple SGS Model in detail Smagorinsky. First I chose pisoFoam as base for my solver and added the driving pressure gradient. This gives me the green curve for reynolds stress u'v' in the figure shown at the end of this thread. Than I tested the channelFoam solver which produces another solution, in the figure the light blue curve. For comparison I solved the same case, identical input control files and mesh with OpenFOAM 1.5 and channelOodles which gives me the dark blue curve ... and what a surprise another solution. So I checked the differences in solver and sgs model in detail and I am able to reproduce the OpenFOAM 1.5 solution with OpenFOAM 1.6 (red curve with crosses) if I:

1) put the SGS model in my solver: turbulence->correct(); in front of my PISO loop instead at the end like in the new pisoFoam solver.

2) remove the call of updateSubGridScaleFields(fvc::grad(U)); in the constructor of the SGS model. This have to be done too, If not one will get another different solution :(

I also checked the solver without an SGS model (setting ck to zero for Smagorinsky) and I am getting the same results for both version. So in my opinion there is a serious bug in the LES interface somewhere. But I do not have the knowledge to understand everything whats going on in the SGS Model, wallModel, transport model and everywhere else. I use OpenFOAM on HPCs, where I am not able to make changes to the main source files. So please fix this bug as soon as possible or tell me my mistake. I really need a working LES interface.

kind regards,
Fabian

Attachment 1411

Edit: I am now checking the GIT Version OpenFOAM 1.6.x and there are some changes in the LES interface. A first test with channelFoam of the 1.6.x version leads me to the same results like the channelFoam of version 1.6. So there is still a gap between OpenFOAM 1.5 and 1.6 :(

I tested the same cases with OpenFOAM 1.6.x but the results are as bad as in OpenFOAM 1.6. I attached the graphs below. No body else experienced problems with the LES interface?

kind regards
Fabian

Attachment 1413

Sry guy its me again. I still try to find out whats wrong with OpenFoam. I noticed that the LESProperties file has changed between the channel395 case of OpenFOAM 1.5 and 1.6.x. So I thought: lets try this may be a parameter wasn't set correctly. And let me say, I am sry for that, but the results sucks. I have now 2 more curves or in other words 25 OpenFOAM channel cases and up to 10 different results. And I still have no idea why. Please help this OpenFOAM 1.6 makes me crazy. If anyone would try it himself, just take the channel395 tutorial and change LESModel from oneEqEddy to Smagorinsky and enjoy :-D

Attachment 1415

Changing the location of the LES model correction will affect the temporal truncation error and hence the solution. If you want to make sure the solution is independent of this choice you could add an outer loop as in pimpleFoam and check that the system of equations is completely converged at each time-step.

> 2) remove the call of updateSubGridScaleFields(fvc::grad(U)); in the constructor of the SGS model. This have to be done too, If not one will get another different solution

This doesn't make much sense because the initialization procedure for nuSgs should not affect the final averaged solution unless the averaging started from the beginning of the run and was not for long enough. How long are you running these cases before starting the averaging and how long are you averaging for?

The tutorial case is setup to run with valDriest delta, are the differences in the results sensitive to this? If you run with a simple geometric delta do you get the same level of difference?

H

I give you right from the numerical point of view. But if you look to the solver, e.g. pisoFoam the velocity and pressure field after the piso loop is the same as in front of the piso loop. So it have to be equal if one puts the call of the correct function in front of the piso correction or behind. This is not equal for initialization and in my optinion this is the reason for calling updateSubGridScaleFields(fvc::grad(U)); in the SGS Model.

You are right too but the main problem is the following: I use the channel395 tutorial which starts from given solution which hopefully is converged to a quasi steady state. And I start averaging from the given start point in the tutorial. In other words I take this tutorial for channelFoam, change the SGS Model and type into my command line: foamJob channelFoam . And it scares me that the initialization from a restart point affects the averaged solution. Because on our HPC system we have a batch system and there is no way to do an averaging over the whole necessary period with only one run. I have to restart during averaging. So I have to be sure that the initialization is correct. And the comparison of this simple tutorial between 1.5 and 1.6.x shows that there are differences, which may be affected by different initialization of the SGS Model.

Not not yet, but I will try this in a few minutes, after a coffee break :-D

The fields provided in the tutorial case are VERY old and should not be considered consistent with the current solvers. They are provided to avoid the need to run the case from a laminar flow through transition and to a fully-developed turbulent state. To run the case you should start from the fields provided, run for a while to develop the flow to be consistent with the solver before starting the averaging and then average for a sufficient time to obtain smooth statistics.

H

Allright, very interesting. So if I understand it right, I should run it a while starting from the provided fields in the tutorial and than start averaging and compare the results for both versions (1.5 and 1.6.x). Allright I will try this, hopefully the results will aproach :-D

thx
Fabian

So I have followed your suggestion henry. I took the channel395 tutorial, run 5000 time steps and than I started averaging for 5000 time steps. The results are shown in the attachment. Again channelOodles (green curve) looks very good and channelFoam (red curve) is to far away from the OpenFoam 1.5 solution. For my understanding there should be no difference, because in my opinion there is no difference between the solvers. If I set ck for Smagorinsky to zero I will get the same results for OpenFoam 1.5 and 1.6. So the difference must be located somewhere in the LES model.

kidn regards,
Fabian

Attachment 1435

I am not convinced that 5000 time steps is enough either for the start-up period or for the averaging. Could you look at the evolution of the averaged solution for both to check this? Could you also post the results you obtained setting ck to 0?

H

I perform now two calculations and double the time for start-up period and averaging. I will check the evolution for the avergaged solution as soon as possible. I think your right, because the graphs for both versions shows some waves, which could be a result of less averaging time. But still I am wondering why the solutions after identical time steps, with identical start-up are so different between 1.5 and 1.6. I will post my test with deaktivated SGS model. This is simple done by setting ck to zero and it leads to identical curves for the Renoldy stresses.

kind regards,
Fabian

Attachment 1436

It would be useful to know if the difference in behavior is related to the use of the van Driest damping function; could you try with a purely geometric delta?

H

How do I activate only geometric delta? In my LESProperties dictionary there is a line which indicates VanDriest wall damping:

delta vanDriest;

But in my Userguide there is only:
PrandtlDelta, cubeRootVolDelta, smoothDelta and no VanDriest delta. Do you know the keyword for geometric delta?

kind regards,
Fabian

cubeRootVolDelta is purely geometric as it is simply the cube-root of the cell volumes.

H

Note that the Smagorinsky coefficients set in the LESProperties file in the 1.5 tutorial case are not the same as the high-Re default values in the 1.6 code and to make the case run the same you will need to set the coefficients to be the same. We should set the coefficients in the LESProperties file in the 1.6 tutorial to override these defaults appropriate for the low-Re of this case as in the 1.5 tutorial case.

H

Yeah I allready noticed the difference and corrected it by supplying them in the LESProperties dictionery. Thx for the hint with the geometric delta. cubeRootVolDelta was the one I thought, but I wasn't sure. So I will try this an post the results.

thx
Fabian

Edit: Allright I am getting very different results if I increase the startup and/or the averaging time. So I will use the weekend to produce a solution with a high startup time (50000 time steps) for both versions, without vanDriest damping, to ensure that this is not the reason for the differences and on Monday I will do a long averaging run. Hopefully this will bring me one step forward :/

Alright, the tests with extended startup and averaging time are finished and the solution looks good in my opinion. There is a little gap, but this seems to be ok. The parameters of this run are: startup: 100000 time steps, averaging: 100000 time steps, Smagorinsky with adjusted constants on 1.6.x from 1.5, cubeRootVol for delta.
The next interesting thing would be: whats happens to the solution if I restart during averaging several times. So I will try this and see whats happens.

Fabian

Attachment 1476

It would also be interesting to re-check the results with van Driest damping and sufficient averaging. If those results are still dissimilar we will investigate in detail.

Thanks

H

What do you mean by "sufficient averaging"? I will start an additional case with OpenFoam 1.6.x and van Driest damping and we will see if the solution changes or not :D

Fabian

Edit: Ignore my question ... I know what you mean by "sufficient averaging" sry, I am not really awake today ;)

Test with the restart during averaging are finished and the results looks very good. See the brown graph below. So with pure geometric delta everything works fine. I will check whats happening when I use the vanDriest wall damping.

Fabian

Attachment 1489

The tests with VanDriest damping looks very good too. There are allways very small differences, but this could be due to some rounding errors or initialization, I dont know but they are very small and thus in my opinion neglectable. The more interessting thing to me, is that the results with vanDriest damping do not fit the DNS results of Moder, Kim and Masour for Re_tau = 395. But I think this is a problem of the mesh, timestep and so on. So at least: Please update the tutorials in the next OpenFoam versions. This is always the first try, if one detects some errors, and they have to work. Sry for any inconveniences :-D

kind regards,
Fabian

Attachment 1509

Thank you for performing these rigorous comparison tests; it is always useful to have our developments independently checked in the way.

Concerning the comparison with DNS we have also noticed that this case is not very accurate but the mesh is quite coarse. The other issue is that there is more than one way to normalize the results and the different approaches give different comparisons so you might want to investigate this.

H

I agree with you, the mesh seems to be to coarse for a validation. But I dont want to investigate this further. I will go back to my own model. Now I know OpenFoam 1.6 and my modified pisoFoam solver works correctly and there should be a problem somewhere in the implementation of my model. May be it is very simple, e.g. to short averaging time or something like that. I will have a look :-D Thx for your help.

kind regards,
Fabian