
[Sponsors] 
August 21, 2013, 09:16 
Smagorinsky LES: output and average k value

#1 
New Member

Dear Foamers,
I am simulating a wind turbine in a wind tunnel, and I am performing a LES simulation using the Smagorinsky model and the pisoFoam solver. I am also averaging the velocity, pressure and subgrid scale kinetic energy ( p U k) by the use of the fieldAverage library. While everything works fine for p and U, I cannot average the k and I get the following error: Code:
> FOAM FATAL ERROR: Requested field k does not exist in the database From function Foam::fieldAverage::initialize() in file fieldAverage/fieldAverage/fieldAverage.C at line 104. FOAM exiting If I switch LES model to oneEqEddy f.ex., then everything runs smooth, I get all the averages including k and the k is written in the output time folders. My question is: is there a way to force OF to write and average the k value when using the Smagorinsky model? I have tried to use the swak4foam utility but did not get any luck. Thanks in advance for your help! Fabio 

August 21, 2013, 10:39 

#2 
Senior Member
Kent Wardle
Join Date: Mar 2009
Location: Illinois, USA
Posts: 218
Rep Power: 14 
Smagorinksy model does not use/write field k.


August 21, 2013, 11:20 

#3 
Senior Member
Kent Wardle
Join Date: Mar 2009
Location: Illinois, USA
Posts: 218
Rep Power: 14 
Guess maybe I could be a little more helpful. The Smagorinsky is a zeroequation model so there is no additional transport eqn for k. You can compute this but it depends if you are looking for the resolved TKE or just the subgrid TKE. Take a look here:
http://www.cfdonline.com/Forums/cfx...tml#post334843 From that post (which is just a copypaste from Pope, 2000) the total TKE is the sum of the resolved part: k_res = 1/2 * avg(U) * avg(U) and the SGS part k_sgs = 1/2 * tau_ii_R = 1/2 * avg(U^2) * [avg(U) * avg(U)] I think this is correct anyway. Hope this helps. Kent 

August 21, 2013, 19:14 

#4 
New Member

Hello Kent,
thanks for your help, your solution is good enough for my case. I can use swak4foam to output k at runtime according to the equations. Anyway I'd like to know, out of curiosity/future need, one more thing. From Smagorinsky.C,line 115, looks to me that k is actually calculated as: Code:
// Return SGS kinetic energy // calculated from the given velocity gradient tmp<volScalarField> k(const tmp<volTensorField>& gradU) const { return (2.0*ck_/ce_)*sqr(delta())*magSqr(dev(symm(gradU))); } // Return SGS kinetic energy virtual tmp<volScalarField> k() const { return k(fvc::grad(U())); } Fabio 

August 22, 2013, 06:33 

#5 
New Member

I have been thinking over your answer, and I think that I cannot calculate the subgrid scale energy in that way, since it is impossible to calculate the average of the squared velocity field:
0.5*avg(U^2) which is also one of the reasons why we need subgrid scale modeling, if I am not wrong? Fabio 

August 26, 2013, 16:35 

#7  
Senior Member
Join Date: Nov 2012
Posts: 168
Rep Power: 7 
Hello,
About obtaining the k_sgs, why not directly use the following? turbulence>k(); Quote:


November 28, 2013, 05:38 
mistake

#8  
New Member
Giulio S
Join Date: Oct 2013
Location: Italy
Posts: 4
Rep Power: 6 
I would like to report a mistake:
k_sgs== 1/2 * tau_ii_R = 1/2 * filtered(U U)  [filtered(U) * filtered(U)] But you don't know the filtered(U U) so perhaps you have to use: return (2.0*ck_/ce_)*sqr(delta())*magSqr(dev(symm(gradU))) which is very similar to the rate of transfer energy yo the residual motion: Pr=2*nuSgs* filtered(Sij) * filtered (Sij) But I attend any suggestion Best Regard Giulio Quote:


March 8, 2014, 09:41 

#9 
New Member
Zhong, Hong Min
Join Date: Oct 2012
Location: UESTC, Chengdu in China
Posts: 8
Rep Power: 7 
Dear OFers, for kinetic energy K in LES, my understanding is (maybe not correct),
1. Note that if u = Filtered(u) + uprime, then resolved velocity vector Filtered(u) is not usually equal to Filter(Filter(u)), while Filtered(uprime) isn't equal to zero, where u is turbulence instantaneous velocity field, or u = u(x, t). However, <uprime> maybe is equal to zero under the isotropic assumptions, <.> denotes Reynolds averaging or time averaging. 2. According to Giulio and Pope, kinetic energy K_total = k_res + k_sgs, where k_res is the resolved kinetic energy and k_sgs the subgrid kinetic energy. Firstly, for k_sgs, if Smagorinsky model, then k_sgs = (2.0*ck_/ce_)*sqr(delta())*magSqr(dev(symm(gradU))), note that it is obtained by the resolved velocity field Filtered(u) (=U). So k_sgs can be got easily. Secondly, k_res = 0.5*Filtered(u)_i*Filtered(u)_i = 0.5U_i*U_i, note also that it is instantaneous or is at current timestep, obviously，we can get it. However, most of the cases, we hope to get <K_total> because we simply concern a equilibrium flow rather than a evolutive flow. If really so, we can ignore <k_sgs> and the reason is that the scale of the subgrid energy vortex (k_sgs) is much smaller than the resolved vortex scale (k_res) based on the Kolmogorov energy cascade. In brief, for OFers we can get <Ktotal> = <k_res> by the "UprimeMean" in the ontheflyprocessing founction "fieldAverage" specified in controlDict file. if any wrong , don't hesitate to correct me please Last edited by iamacloud; March 19, 2014 at 11:20. 

June 16, 2014, 04:28 

#10  
New Member
Hans Barósz
Join Date: May 2014
Posts: 22
Rep Power: 5 
Quote:
Especially around y+=11 the sgs energy rises, in my case to nearly 50% of the resolved kinetic energy. Just my thoughts. Any suggestions? 

December 23, 2014, 19:57 

#11 
Member
Join Date: Feb 2014
Posts: 55
Rep Power: 5 
Hi All,
I am also in this problem of calculating Ksgs and Kres. But i am unable to get K field as an output from Smagorinsky Model. Could anyone help me out there?? I agree with HanSolo. I think that instantaneous Ktotal = 1/2(u^2 v^2 w^2) , and Ksgs is calculated from the LES model. Kres = Ktotal  Ksgs. Last edited by Uyan; January 13, 2015 at 16:01. 

December 25, 2014, 10:12 

#12  
New Member
Hans Barósz
Join Date: May 2014
Posts: 22
Rep Power: 5 
Quote:
To do so, create a new field in createFields.H, for example: Code:
volScalarField ksgs_ ( IOobject ( "ksgs", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh, dimensionedScalar("ksgs", dimVelocity*dimVelocity, 0.0) ); Code:
while (runTime.run()) { #include "readTimeControls.H" #include "CourantNo.H" #include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; //  Pressurevelocity PIMPLE corrector loop while (pimple.loop()) { #include "UEqn.H" //  Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" } if (pimple.turbCorr()) { turbulence>correct(); } } ksgs_ = turbulence>k(); runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } With that new field, you can and of course have to average to get proper results. So in your /systems/controlDict inside the functions>fieldAverage use Code:
ksgs { mean on; prime2Mean off; base time; } I got very good results with this procedure. Dont forget to rename your utilities! 

April 2, 2016, 00:13 

#13 
Senior Member
Freedom
Join Date: May 2014
Posts: 205
Rep Power: 6 
Thank you for clarify this issue with the code manner. Thank you, Hans!
Now I want to know that in LES, the resolved kinetic energy should be a larger part in total kinetic energy (say 80%). That means, Kres/(Kres+ksgs)>80% Then, this criteria should be satisfied in all of the computational domain, right? Did you have some examples of the result you have plotted about the ratio of the resolved and SGS part of the turbulent kinetic energy? Could you past some figures or do you have the related papers? Do you get the result along a line of domain, or over the entire domain? Thanks! Best, Xu 

April 2, 2016, 13:26 

#14 
New Member
Hans Barósz
Join Date: May 2014
Posts: 22
Rep Power: 5 
Hi Xu,
you can find the results of my LES Simulations with pimpleFoam and SmagorinskyModel attached in this post. The black line shows the DNS Results from [1]. Red line belongs to a coarse mesh, green medium and blue finest mesh I have used. All three meshes have a yPlus of 1 at first cell. Other mesh statistics:
In total, the three different meshes resolved:
I have done the average for 200 flowthroughs with Co_max = 0.3 after having a wellsuited initial solution (took me 350 flowthroughs with Co = 0.5 with the initial solution from the tutorial. I was not using the perturbUutility you can find on this forum  I didnt feel confident with that). For postprocessing I have used the postChannelutility, that does an average over the x and zdirections of the domain (if I remember correctly). Then you have the results over a line in wallnormal direction y. These are parts of my master thesis I have finished 2 years ago. Greetings.  [1] Iwamoto, K.: Database of Fully Developed Channel Flow. THTLAB Internal Report, No. ILR0201. THTLAB, Dept. of Mech. Eng., University Tokyo, 2002 

April 2, 2016, 21:55 

#15 
Senior Member
Freedom
Join Date: May 2014
Posts: 205
Rep Power: 6 
Many thanks for your detailed reply, Hans.
I do not know if I understand it correctly. Let me rephrase your results. For resolved kinetic energy Kgs, we calculate it as follows, Kgs = 0.5*(<U>^2 + <V>^2 + <W>^2) For unresolved kinetic energy (SGS) part, we calculated it with the LES SGS model, e.g. Smagorinsky used here, Ksgs = (2.0*ck_/ce_)*sqr(delta())*magSqr(dev(symm(gradU))); I do not know whether it is an error in this attached figure. The DNS should calculate ALL of the energy, so, ksgs_DNS should be equal to ZERO. Rather than equal to Kg or Ktotal. Then, if we want to know how much the kinetic energy is solved with our specific discretization scheme and mesh resolution, we can calculate the turbulence resolution as follows, M = Kgs/(Kgs+ Ksgs) For LES, M should be larger than 80% (Pope, NJP). ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~ But the confusing thing is, the above criteria is applied for EVERY cell, that means in EVERY cell, the resolved kinetic energy should be larger than 80%. But in your case, if I understand it correctly, you sample the data in the y+ direction, when in the other directions (x+, z+), you do average on them. In this case, how can you ensure the results you got can be applied for all the computation domain in your simulation? I mean, if you use the stretched mesh in the flow direction, then in the downstream, the mesh is coarse, even if you get the 90% resolved TKE in the upstream in the spanwise direction, then the results should not be convincing. Because in the downstream of the spanwise direction, there may be only 60% TKE is resolved where the mesh is coarse. Thus, my suggestion is to extract the data from the xy plane in the flow direction, in this case, we can get the result with all levels of mesh resolution. My result is attached. Many thanks for your kind reply. I am looking forward to hear from you for your suggestions. Best regards, Xu 

April 2, 2016, 22:00 

#16 
Senior Member
Freedom
Join Date: May 2014
Posts: 205
Rep Power: 6 
Also do you plot the figure with gnuplot? What the size of these figures (0.3, 1) and how can you replace the number 5 with the character "Kg+". Could you share your plot code? Or could you describe it how to get it? Thank you!
Now I am not familiar which size is good to write the research paper. The figure should be located in one column. Xu 

April 8, 2016, 13:45 

#17  
New Member
Hans Barósz
Join Date: May 2014
Posts: 22
Rep Power: 5 
Quote:
Hi Xu, sorry for late reply. I fear that I cant answer your questions to your satisfaction  too much time has past since I have put all my effort in this topic. First of all you are right. The DNS plots in the figues with k_gs and k_sgs are not correct if you just look at this picture. But I have mentioned in the text in my thesis that I plot the DNS data there to give the reader the chance to see the rate of k_gs and k_sgs compared to the DNS data on the first sight. The sencond part of your question is difficult for me to answer. We have to seperate timeaveraging and areaaveraging. Quote:
As the following step, in my opinion it is allowed to do an areaaverage in x and zdirection with the postChannel utility to get the data that is shown in the picture i attached because there should be no dependency in those directions of the flow. I think the postChannel utility does exactly this, an average in x and zdirection, but I am not 100% sure. In postChannel.C it is said: Code:
// Average fields over channel down to a line # include "collapse.H" If you dont do the areaaverage and just get the results over a line in ydirection e.g. in paraFoam at an arbitrary x and z position, and plot this nonaverage energies over y, there should be just a very very small difference to the averagedresults (if you have a well suited flow development and enough flowthroughs for timeaverage) If it is the case  in my opinion  it can be said that the results can be applied to the whole domain. Im sorry that I cant say anything to your attaced figure. I havent seen such a plot in my study. 3rd part: Yes, I have used gnuplot to make my figures. I have made the size of all figures in a way that they fit well on a DIN A4 page. I am not familiar with scientifi papers either. I configured it with: Code:
set term tikz standalone color solid size 14.8cm,15cm font ',12' set output 'kPlus.tex' To replace a certain axisnumber, just say: Code:
set ytics add ('$k_{gs}^+$' 5) Note that it is LaTex code. gnuplot will generate a file 'kPlus.tex' as mentioned above. Then I needed to translate this with my latex editor and it gave me a pdf with the figure and the text. 

April 8, 2016, 22:12 

#18 
Senior Member
Freedom
Join Date: May 2014
Posts: 205
Rep Power: 6 
First, thank you for your very kind and detailed reply. Really a nice guy. :>)
I do not have confusion with the timeaveraging and areaaveraging. :>) I mean in the x direction (the flow direction), if you use the stretched mesh, then for a give x point, then you sample the data along the spanwise direction (y direction), you may get an unreal result. In other words, if the point x is in the upstream where you have refined your mesh, then you may get a very small SGS kinetic energy. However, the possible (also bad) thing is that in the downstream where you stretch your mesh (much coarse than that in the upstream), then the SGS kinetic energy can be very large, say below 80%. Then how can you justify your result if you only based your data in the upstream flow? So, to resolve the potential problem, my suggestion is to sample the data over the domain OR only in the xy plane. In this way, we have covered all possible mesh. The result is more believable. But if the uniform mesh is used in you case, the problem I list above will not happen. Since the turbulence length scale in the upstream will be certainly much small than that in the downstream. This paper is suggested to read: LES OF THE SYDNEY SWIRL FLAME SERIES: AN INITIAL INVESTIGATION OF THE FLUID DYNAMICS. Many thanks for your kind reply. I really know it is hard to answer something if that has past a long time. Have a nice day! My best wishes, Xu 

April 10, 2016, 03:42 

#19 
New Member
Hans Barósz
Join Date: May 2014
Posts: 22
Rep Power: 5 
Ah ok, I think I got what you mean.
I have indeed used an uniform mesh in x and zdirection. Just in wallnormal ydirection I have had configured a ratio. The xPlus and zPlus values I have posted at April 2, 2016, 12:26 are constant for the whole mesh. For yPlus I had 1 at the first wall cell with increasing yPlus to the middle of the domain. A very nice discussion :) 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Utilities: post average turbulence fields and create turbulence fields for LES  Hanzo  OpenFOAM Running, Solving & CFD  10  August 18, 2017 18:33 
y+ and u+ values with lowRe RANS turbulence models: utility + testcase  florian_krause  OpenFOAM  110  April 21, 2016 11:54 
[swak4Foam] outputTime in Swak function  immortality  OpenFOAM Community Contributions  16  February 14, 2016 17:02 
RNG diverged during the analysis the flow over a multi element airfoil, why?  s.m  OpenFOAM Running, Solving & CFD  0  August 5, 2013 08:39 
average weighted area of volume fraction at output  umar959  Main CFD Forum  0  January 20, 2012 03:47 