[pierellozzo]

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 sub-grid 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:

--> 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 run the simulation without averaging the k, everything runs fine, but I do not get any k file inside the time folders.

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

[kwardle]

Smagorinksy model does not use/write field k.

[kwardle]

Guess maybe I could be a little more helpful. The Smagorinsky is a zero-equation 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 sub-grid TKE. Take a look here:
http://www.cfd-online.com/Forums/cfx...tml#post334843

From that post (which is just a copy-paste 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

[pierellozzo]

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:

  //- 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()));
       }

Therefore, I am wondering if there is a way that k can be "routed" to be an output of the calculations and not only an internally referenced value.

Fabio

[pierellozzo]

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

[pierellozzo]

Any help with this issue would be greatly appreciated

Fabio

[hz283]

Hello,

About obtaining the k_sgs, why not directly use the following?

turbulence->k();

Originally Posted by kwardle

Guess maybe I could be a little more helpful. The Smagorinsky is a zero-equation 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 sub-grid TKE. Take a look here:
http://www.cfd-online.com/Forums/cfx...tml#post334843

From that post (which is just a copy-paste 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

[GiulioS]

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

Originally Posted by kwardle

Guess maybe I could be a little more helpful. The Smagorinsky is a zero-equation 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 sub-grid TKE. Take a look here:
http://www.cfd-online.com/Forums/cfx...tml#post334843

From that post (which is just a copy-paste 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

[iamacloud]

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 <K-total> = <k_res> by the "UprimeMean" in the on-the-fly-processing founction "fieldAverage" specified in controlDict file.

if any wrong , don't hesitate to correct me please

[HanSolo123]

Originally Posted by iamacloud

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 <K-total> = <k_res> by the "UprimeMean" in the on-the-fly-processing founction "fieldAverage" specified in controlDict file.

if any wrong , don't hesitate to correct me please

I dont agree with that, The SGS are indeed much smaller than the GS, but not that small that we can ignore it in a classic LES, where ~ 80% of the total kinetic energy should be resolved. That means, that 20% remains in the SGS part which cannot be ignored!

Especially around y+=11 the sgs energy rises, in my case to nearly 50% of the resolved kinetic energy.

Just my thoughts. Any suggestions?

[Uyan]

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.

[HanSolo123]

Originally Posted by Uyan

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??

You have to modify your solver to get k_sgs. k_sgs is just calculated temporarily inside the Smagorinsky model and then disappears after each timestep because the closure problem was solved.

To do so, create a new field in createFields.H, for example:

volScalarField ksgs_
(
    IOobject
    (
        "ksgs",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("ksgs", dimVelocity*dimVelocity, 0.0)
);

Then in your solver, for example pimpleFoam, you have to write out the modeled k (=k_sgs) before the runtime-command. So in pimpleFoam.C it should be:

while (runTime.run())
    {
        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "setDeltaT.H"

        runTime++;

        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity 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;
    }

This way, the solver creates a new field called ksgs each time step, which contains the subgridscale energy from your turbulence model. (I've just worked with pimpleFoam+Smagorinsky/oneEq so not sure about other combinations).
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

ksgs
            {
            mean      on;
            prime2Mean  off;
            base        time;
          }

This will average your subgridscale energy. For easy postprocessing and if you do the channelflow, I recommend you to modify the postChannel utility, so it will output your k_sgs over y (wallnormal distance).

I got very good results with this procedure. Dont forget to rename your utilities!

[wenxu]

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

[HanSolo123]

Hi Xu,

you can find the results of my LES Simulations with pimpleFoam and Smagorinsky-Model 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:

• coarse: xPlus = 40, zPlus = 26
• medium: xPlus = 31, zPlus = 19
• fine: xPlus = 26, zPlus = 15

In the figure you can see at yPlus = 11, that kPlus_sgs is about 2 and kPlus_total is about 4. So kPlus_sgs is about 50% what I mentioned in my post from June 16, 2014 03:28.

In total, the three different meshes resolved:

• coarse: 86%
• medium: 88%
• fine: 90%

of total kinetic energy (k_gs/k_tot). [k_gs = grid scale = resolved part].

I have done the average for 200 flow-throughs with Co_max = 0.3 after having a well-suited initial solution (took me 350 flow-throughs with Co = 0.5 with the initial solution from the tutorial. I was not using the perturbU-utility you can find on this forum - I didnt feel confident with that).

For post-processing I have used the postChannel-utility, that does an average over the x- and z-directions of the domain (if I remember correctly). Then you have the results over a line in wall-normal 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. ILR-0201. THTLAB, Dept. of Mech. Eng., University Tokyo, 2002

[wenxu]

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 x-y 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

[wenxu]

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

[HanSolo123]

Originally Posted by wenxu

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 x-y 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

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 time-averaging and area-averaging.

In this case, how can you ensure the results you got can be applied for all the computation domain in your simulation?

I think it is possible because of the boundary conditions of the channel flow. The sides of the computation domain are cyclic. If you have a good initial solution and if you time-average the flow in the whole domain for enough flow throughs, there is just a wallnormal (y) dependency of the flow profile.

As the following step, in my opinion it is allowed to do an area-average in x- and z-direction 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 z-direction, but I am not 100% sure.
In postChannel.C it is said:

        // Average fields over channel down to a line
#       include "collapse.H"

so the answer should be in collapse.H.


If you dont do the area-average and just get the results over a line in y-direction e.g. in paraFoam at an arbitrary x- and z- position, and plot this non-average energies over y, there should be just a very very small difference to the averaged-results (if you have a well suited flow development and enough flow-throughs for time-average)

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:

set term tikz standalone color solid size 14.8cm,15cm font ',12'
set output 'kPlus.tex'

The code for the plots is just standard .
To replace a certain axis-number, just say:

set ytics add ('$k_{gs}^+$' 5)

It will replace at y-axis at position 5 the digit with '$k_{gs}^+'.
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.

[wenxu]

First, thank you for your very kind and detailed reply. Really a nice guy. :>)

I do not have confusion with the time-averaging and area-averaging. :>)

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 x-y 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

[HanSolo123]

Ah ok, I think I got what you mean.

I have indeed used an uniform mesh in x- and z-direction. Just in wallnormal y-direction 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 :-)

[gu1]

Originally Posted by HanSolo123

You have to modify your solver to get k_sgs. k_sgs is just calculated temporarily inside the Smagorinsky model and then disappears after each timestep because the closure problem was solved.

To do so, create a new field in createFields.H, for example:

...

I got very good results with this procedure. Dont forget to rename your utilities!

Hans Barósz,

I can store only 'k'? How I can do this? I would like to be able to view the 'k' profile in an internal pipe flow ...