LES of channel flow: data, case files, technical report.
 

Dear foamers,

I would like to share the results of a simulation campaign that was performed by us here at Uppsala University.
I was supposed to present this at the ongoing OF workshop, but I couldn't make it because of visa issues :(.

We have performed LES of channel flow using OpenFOAM. In our study we have computed a large amount of statistical quantities, including high-order moments, two-point correlations etc. Also, we have considered three different grids, so one clearly sees how the grid affects the results, and how the results converge towards DNS data as the grid gets better.

We have collected all the acquired profiles in this repository:
https://bitbucket.org/lesituu/channel_flow_data/

A ready-to-run case file can also be found there, so the data should be relatively easily reproduced by others.

A techincal report that documents our study is available for download here:
http://www.it.uu.se/research/publica...orts/2015-014/

This project was basically a way for me to start my Phd and learn about LES and OpenFOAM at the same time, so the starting point was the existing channel flow tutorial. The majority of the setting were left unchanged.
In order to compute most of the statistical quantities we used swak4foam.
Since I think it is quite common for people to start out with channel flow, I hope that the repository we provide will make life a bit easier for future new-commers :). In some sense it can also be seen as a swak4foam tutorial.

Any kind of feedback is of course appreciated.

P.S. I was not sure if I should start a new thread or just post this into the existing thread for LES of channel flow. I thought that this thread can serve for discussing the particular resuts/setup/whatever found in our work and the general thread can continue just being general. I leave it to the moderators to decide whether that was a good call :).

Hi Timofey,

First of all thanks for this handy report it was very useful for me to read it.
I am trying to validate the k-ε model in openfoam based on the same DNS data.

So far I have one comment regarding your report I am not sure if I misunderstood your point or it is really a mistake. In table 3.2 you put the number of cells along axis and corresponding Y-plus. I have calculated Y-plus based on this information and the equation for the mesh grading used in blockMesh and I found that the the corresponding Y-plus values for this grid should be the double. So I think that these numbers correspond to number of cells in δ not in H, I am I right?

Hi Mohamed!

I did a quick test using the UMean profile for mesh M1 and got a a value of 1.74.
We have
y^+ = y_p*sqrt((nu+nuSgs)*U_p/y_p)/nu,
where the _p values are in the cell center closest to the wall.
I ignore the nuSgs, because I don't want to search for the value right now.
Then we get
y^+ = 0.0048*sqrt(2e-5*0.0127449/0.0048)/2e-5 = 1.7489

In the table, the value is 1.9 which is due to the nuSgs I suppose.

In any case, the numbers in the table correspond to h, not delta.
Since you got double values, have you perhaps taken y_p to be the location of the first *point* of the mesh (i.e. the location of the first off-the wall face), not the cell center? Just a guess. Anyhow, please share the values you use and we'll get to the bottom of this :).

Best,
Timofey

Hi Timofey,

You are right I am using the height of first cell to calculate Y+, if the value in the table are for the first cell center, then that is perfect.

I have one more question regarding initial and boundary conditions. The value you have used as initial conditions (for U, p, nuSgs and k) seems to be from other simulation (i.e you run for specific time steps to initiate these fields) not calculated based on analytical formula (velocity defect law for example), am I right? if yes can you give some elaboration on this point.

Regarding boundary condition, you are using cyclic B.C., I am not familiar with this B.C. and really I can't understand how it works. for example does it require additional information to define the velocity profile at time step 0?

for my case I am using nonuniform values for U, k and ε at the inlet which I've calculated from the analytical formulas and fix that for all time steps, do you think that this approach is physically correct?

Hi!

Regarding the ICs, I used the channel395 tutorial case that is availabnle in openfoam as a starting point. Some initial field is given there, and I did not bother looking deeper into that since, well, it worked for me.
However recently I had to do work initialization and the perturbU utility which you can google for seems to be a great tool. I have successfully used it to initialize an LES.

The cyclic bc is just a periodic boundary condition, so what flows out of the outlet comes back into the inlet. To make this work you have to make sure the domain is sufficiently large.

BUT. You say that you want to verify the k-eps model, so you are doing RANS. Most of this becomes irrelevant to you then. For RANS channel flow is a 1-d problem: you solve for the mean, and the mean only depends on y! You just need a line in the wall-normal direction and that is your mesh :p.

I think you will benifit from checking out this report:
http://uu.diva-portal.org/smash/reco...89&dswid=-9746
It deals with exactly your problem: testing how different RANS models perform for the channel flow!

Best,
Timofey

Thanks Timofey, that is extremely helpful :)

Dear Timofey,

thank you for your contribution regarding this specific topic.

I have two questions.

1.) You have used oneEqEddy model together with wall damping? Or I missed something?
2.) I tried to perform LES channel flow computations and I cannot get a proper averaged pressure gradient (i.e. wall shear stress or uTau) values. When I compare computed uTau values with DNS, errors are around 20%. These errors seems to be to large as other researchers reported errors below 2% (for same subgrid models and similar mesh resolutions). Did you manage to calculate a proper uTau value?

Best regards!

1. Yes, see https://bitbucket.org/lesituu/channe.../LESProperties


I'd recommend using WALE though. Or no SGS model at all.



2. You can see the Re_tau that I got in the report. For a fine mesh, you should get good results. I'd recommend delta z+ = 10 and delta x+= 25 as a good start.

Dear Timofey,

thank you.

However, I'm still curious regarding the pressure gradient value in channel395 computations. For instance, when I tried to run channel395 tutorial in of6, the average pressure gradient is still significantly lower than theoretical (i.e. computational is equal to cca 4e-05, whereas theoretical is equal to cca 6.2e-5). These results produce significant error when computational friction coefficient is compared with theoretical value. I've read Phd thesis written by Eugene de Villiers and his results were much better for same setup.

Can you please comment this situation?

Thank you in advance! I appreciate your help.

Best regards!

How do you measure the pressure gradient? Do you mean the value reported by the fvOption?

Dear Timofey,

I tried to run tutorial case for larger amount of time (i.e. for 7000 s). After that I exported the pressure gradient values from log file and time average it using libreOffice calc.

P.S. I use foam-extend 3.2 so I use channelFoam solver as foam-extend distribution does not have implemented fvOptions.

OK, I see. I have to say that I've never looked at the pressure gradient values that much, but rather at u_tau, tau_w and Re_tau. How does the error look for those quantities, is it different from what you get for the pressure gradient?

Well, I supposed that uTau can be calculated as sqrt(gradP), where gradP basically represents Tau_w (I've read that in Eugenes PhD). So I used this methodology for calculating uTau, Tau_w and Re_tau. Did I miss something?

How did you calculate these values in your technical report?

Look at how I compute y+ in a post in this very thread. Basically, you can just compute tau_w by definition as (nu + nuSgs) du/dy, where the gradient is computed using the first off-wall value from the mean velocity profile. Note that nuSgs should be 0 on the wall, but I added it so that you can double check that :).

Thank you! I'll check it.

However, I am still curious regarding that pressure gradient. Eugene noted in his PhD (page 163):

To ensure consistency, the streamwise bulk velocity, Ub, through the channel is adjusted to be equal to the average DNS value by varying an imposed streamwise constant pressure gradient, ∂P/∂x. - this basically means Eugene used something similar (or same) to the channelFoam solver.
The time averaged value of this quantity is equivalent to the mean wall shear stress and can be compared to the corresponding DNS value. - Regarding to this statement, the tau_w should be calculated as sqrt(gradP), where gradP is value given in simulation log file.

I mean, yes, the gradP and tau_w are trivially connected in the case of channel flow, so in principle you should get similar results using the two methods. I prefer computing it by definition because


- It works for all flows and not only channel flow.
- It uses a field that is readily available without the need to parse the log.

I agree with your statements.

I'll try to investigate this topic a little bit more and give you some feedback as soon as I find some solution.

One more question. What procedure did you use to calculate the vorticity fluctuations that you compared with available DNS results (i.e. https://turbulence.oden.utexas.edu/d...han395.vortvar )?

I use swak4foam for these type of things, see

https://bitbucket.org/lesituu/channe...em/controlDict


Basically, I create a field defined as curl(U) and then add it to the fieldAverage. I am pretty sure there is a native function object for vorticity as well.

Great. Thank you for all comments. :)

P.S. I've seen your paper "A Library for Wall-Modelled LES Based on OpenFOAM Technology". Really nice paper. :) Will this library be available (or perhaps it is already available?) to OpenFoam community soon?

https://bitbucket.org/lesituu/libwallmodelledles


Enjoy :)

Nice! Thank you! :)

Respected Sir, it was so nice of you to share this topic here at CFD Forum. I am a newbie in OpenFoam and trying my best to learn most out of it but most of the time I suffer in finding the proper answer to my doubt in google.
Anyway I was doing some research with Channel 395 case with some modification like the boundary condition. In an unchanged channel flow the U and p folder look like the attachment shared along with this post. From my obserservation what I depicted is as follows:-
1.Bottom and top wall is having a no-slip boundary condition.
2. The solution has been provided with some initial condition i.e for 60000 cells they have provided some initial values.
3. The flow that is coming out of outlet is again fed to inlet and in this way the cycle goes on.
What I didnt understand is what is the significance of the side using it as a cyclic boundary condition?
Once I m clear with this part then I ll aim for my next target that is feeding some surface value to the inlet patch using mapped field and run the simulation. So in what way can I modify the bc in order to feed the inlet patch with some predefined value instead of the flow coming out from the outlet?

The idea is that the domain is infinite in the wall-parallel directions. We simulate that with cyclic boundaries. If you have walls on the side, that is going to lead to different results, since you'd have associated velocity gradients etc.



Regarding changing the boundary to some sort of inlet, it depends on what you intend to "feed". In any case, you put some sort of Dirichlet bc for U and zeroGradient for the pressure at the inlet, and the outlet set pressure to 0 and zeroGradient for velocity.

I also realized that once I remove the initial condition for all the cells from non-uniform to uniform the simulation is not running properly. But, for my use of timeVaryingMappedFixedValue at inlet I ve to make the initial condition as uniform otherwise how is it possible? Basically the cyclic bc runs when I put initial condition to all the cells but how can I run the same by feeding face-center value to the inlet patch?
Also attached with this post is the boundary field which I m trying to feed to the inlet, u can have a look.

If you'll have an inlet and an outlet then initial conditions don't matter that much at all, they will be flushed out of the domain pretty fast anyway.

Do you have any idea how can I make it working?

I'm sorry but I don't really understand what the problem is.

Hi @Tiam,

Thanks for your paper and the codes that you shared, it's a very old post, so not sure, if it's still active. I am also trying to perform LES of Re_tau = 180 for channel flow in a little bit bigger domain, by imposing the pressure gradient using fvoptions and calculating the pressure gradient as

,

and calculating u_\tau using the Re_\tau value, since h is known, viscosity I predefined for my simulations. Unfortunately I am unable to see any turbulence in my simulations, could you please share your 0 folder? I assume that when you started your simulations it could have been started from 0 time. Any comment or suggestion would be welcomed, many thanks in advance.

180 is difficult to make turbulent. I recommend that you change viscosity and run at a higher Re for a bit, just to get turbulence. Than u can switch tje viscosity back. Also, be ready that when u set the p gradient instead of Ub, it takes a lobger time to reach a statistically steady state.

Hello

Good morning, many thanks for your prompt reply, but is there anything I can do to trigger the turbulence by boundary condition or initial condition?
To save time for getting the statistically stationery state, I am using simpleFoam to get converged results and then mapping the fields to LES solver, it this fine ?
Any further help would be very much appreciated, thanks in advance.

Dear Chandra,

you should artificially introduce turbulence to the flow. There is an utility written by Eugene de Villiers that modifies the initial velocity flow in order to trigger turbulent motion in the channel. Basically, the utility will modify the initial velocity field near the wall boundary to mimic the coherent vortex-like structures based on the paper mentioned in the utility code. I tested it and it works just fine. But have in mind that you should probably modify code depending on the distribution of OpenFoam you use in order to compile it.

Find utility attached to this post. I hope this will help you with your work.

Best regards.

Thanks a lot Camel for your wonderful and helpful comment. I was trying to install it, but after unziping your sent zip file, i could see that it contain Make directory, perturbU.c and perturbUDict, could you please elaborate how to install and use it. I am using OpenFOAM v1906 and/or v2112. Any further comments/suggestions would be very much appreciated, many thanks in advance.

Hi again Camel,

Many thanks for your help in this forum, could you please help me to use the utility perturbU in the OF v1906/v2112. How to install/use this utility, I am still stuglling to get turbulence for the channel flow Re_\tau = 180, with specified pressure gradient condition. I managed to get the steady state solution from RANS, using that I am trying to run LES case. I used WALE scheme for LES but unfortunately seems like that is killing turbulence, thus now switching towards the KEqn model. If this utility could help, then it would save enormous amount of time, hope you understand and will help me, thanks in advance.

channel395 tutorial case
 

Hi Timofey,
I have been trying to run channel395 tutorial case to reproduce the moser’s results at Re_tau=180 initially. I used to run it till t=10s. I also tried to run this case till t=2000s but results were similar as for t=10s. In this tutorial, there are folders like 0 and 0.orig, so I used 0.orig for my simulations. However, I could not reproduce moser’s results yet. I have visited https://bitbucket.org/lesituu/channel_flow_data/ and there is a folder named initial_conditions which does not readable, however, other folders can be readable. Would you (or others) please share the initial_conditions folder again so that I can read it.
Regards,

Thank you for sharing your work ! Useful !