[alberto]

Originally Posted by foamWang

Hi panda,

I updated the fvSchemes with QUICK phi. So far, the simulation doesn't easily blow up. I will try other method to further improve.

Roro

QUICK is not a suitable scheme for LES, and the upwind scheme is even worse. They are too dissipative.

You should use linear, or higher order central schemes. Generally, avoid limited schemes, they're not a good choice if you want accurate LES.

If the solution "blows up" using the linear scheme, it is probably because your grid is not fine enough to have the local Pé < 2, as required by the scheme.

Best,

[xiao]

Thank you, Alberto, for sharing your experience.

It seems that mostly commonly used scheme is still central differencing (at least in the tutorials). Could anyone explain why the cubic schemes are not used? Such as "Gauss cubic" or "Gauss cubicCorrection", or cubic Corrected (ref. User guide 114 to 116)? These schemes are used in the DNS tutorial. Is it because of stability considerations?

Also, why is upwinded schemes not suitable for LES? Intuitively, I would expect a 2nd order upwinded scheme to be better than a 2nd order central scheme. If both are 2nd order, the diffusion should be of the same order. Maybe I am missing something? In User guide pp116, Table 4.10, it says: linearUpwind, QUICK, TVD, NVD are all first/sescond order. Does that mean it is 2nd order for smooth solution and 1st order for shocks? That's what I remember from CFD classes. Then, for incompressible flows, it should be 2nd order, and should perform better than central scheme. Could anyone shed some light on this issue? Alberto, could you explain a bit? Thanks very much!

Best,
Heng

Originally Posted by alberto

QUICK is not a suitable scheme for LES, and the upwind scheme is even worse. They are too dissipative.

You should use linear, or higher order central schemes. Generally, avoid limited schemes, they're not a good choice if you want accurate LES.

If the solution "blows up" using the linear scheme, it is probably because your grid is not fine enough to have the local Pé < 2, as required by the scheme.

Best,

[alberto]

Originally Posted by xiao

Thank you, Alberto, for sharing your experience.

It seems that mostly commonly used scheme is still central differencing (at least in the tutorials). Could anyone explain why the cubic schemes are not used? Such as "Gauss cubic" or "Gauss cubicCorrection", or cubic Corrected (ref. User guide 114 to 116)? These schemes are used in the DNS tutorial. Is it because of stability considerations?

I think it is a tradeoff between stability and accuracy. In a finite volume code, your global accuracy is second order anyway. If you use higher order schemes for the convection term you can have some benefit, but you have also to evaluate how important this benefit compared to the standard linear scheme.

Also, why is upwinded schemes not suitable for LES? Intuitively, I would expect a 2nd order upwinded scheme to be better than a 2nd order central scheme.
If both are 2nd order, the diffusion should be of the same order. Maybe I am missing something?
In User guide pp116, Table 4.10, it says: linearUpwind, QUICK, TVD, NVD are all first/sescond order. Does that mean it is 2nd order for smooth solution and 1st order for shocks? That's what I remember from CFD classes. Then, for incompressible flows, it should be 2nd order, and should perform better than central scheme. Could anyone shed some light on this issue? Alberto, could you explain a bit? Thanks very much!

I'll try to give you some information of why central schemes are generally preferred in LES and then point you to some reference that is surely more authoritative than me

The order of accuracy gives you information on the dependency of the error on the local size of the discretization (See Ferziger and Peric book). So you might be led to think you simply need to increase the order of the scheme to be safe and do not dissipate too much, to avoid losing information on your turbulence structures due to the dissipation. Unfortunately the story is not so short.

R. Mittal, P. Moin (AIAA Journal, 0001-1452 vol.35 no.8 pp. 1415-1417, 2007, doi: 10.2514/2.253) showed that the numerical dissipation of high-order upwind schemes removes a significant amount of energy from the resolved range of wave numbers, affecting especially the high wavenumbers part, which becomes significantly contaminated by diffusion and dispersion errrors.
This did not happen with energy-conserving central schemes, which do not introduce dissipation (diffusion error).

Why this difference? You find an explanation for example in A. Aprovitola, F. M. Denaro, J. Comp. Phys, 194(1), 329-343 ( http://dx.doi.org/10.1016/j.jcp.2003.09.027 ) , who also describes some upwinded-biased scheme that should be suitable for LES. What the say is essentially that if you discretize the derivatives of your equation with schemes that have non-symmetric stencils you obtain modified wavenumbers with imaginary part that does not tend to disappear. This imaginary part is the responsible of the diffusion error that leads to energy loss.

This problem does not appear in central schemes, which are characterized by real modified wave-numbers. The real part is responsible of the dispersion error, leading to energy pile-up, but the diffusion error is not present.

Best,

[panda60]

Dear Heng Xiao,
I agree with you. "dynamic SGS modeling not needing wall function is only theoretical. Using a wall model (like wall damping, or wall stress model) may sometimes be desired". But I also have a lot of questions.

1.smoothDelta also can be called "wall damping" ?
in smoothDelta.H , it says like this:

Class
Foam::smoothDelta

Description
Smoothed delta which takes a given simple geometric delta and applies
smoothing to it such that the ratio of deltas between two cells is no
larger than a specified amount, typically 1.15.
2. Eugene recommended that we use backward for ddtSchemes, and Gauss
filteredLinear for div(phi,U). So I use these for my LES simulation. But my professor
always ask me why I use these schemes, I can't say, Because my C++ is poor. In
OpenFOAM's Program Guide, doesn't have explanation of these schemes. Who can
tell me the exact formula of backward and Gauss filteredLinear ?
3. I also want to talk about inflow condition for LES. Because my LES result is not so good, both my professor and I think it is due to the lack of inflow turbulence. OpenFOAM has directMapped method, directly maps the internal data back to inlet, but this method doesn't take into account boundary layer roughness. Lund developed a method in 1998, but this method is a little complex. And then a Japanese researcher called Kataoka simplefy Lund's method in 2003. Kataoka's formula is very simple, and easy to use, now is widely used in wind engineering. Eugene said maybe G.R.Tabor has add this method to OpenFOAM. Yes , indeed, recently I found G.R.Tabor's paper.

M.H.Baba-Ahmadi, G.Tabor. Inlet conditions for LES using mapping and feedback control. Computer & Fluids. 2009.
G.R.Tabor, M.H.Baba-Ahmadi. Inlet conditions for large eddy simulation: A review. Computer & Fluids. 2010.

I also want to add this to OpenFOAM in the reference of directMapped utility. But my C++ is poor, so I hope somebody could give me some guide.

[alberto]

Originally Posted by panda60

Dear Heng Xiao,
I agree with you. "dynamic SGS modeling not needing wall function is only theoretical. Using a wall model (like wall damping, or wall stress model) may sometimes be desired".

Pay attention to one detail. In OpenFOAM the dynamic model computes the cofficient dynamically but this is averaged on the whole domain. This is not how the original dynamic Smagorinsky model is supposed to work, because the coefficient has to be a local, and not a single value for the whole domain, to make the model work as intended.

The dynamic Smagorinsky model is capable to work correctly in the low-Re regions, because of its local coefficient. If you average on the whole domain, you simply reduce it to a normal Smagorinsky model with a coefficient that is adapted on the base of the global flow conditions at each time step!

You can see the details of the implementation in OpenFOAM taking a look at .../src/turbulenceModels/incompressible/LES/dynSmagorinsky.C, where you read

dimensionedScalar dynSmagorinsky::cD(const volSymmTensorField& D) const
{
    volSymmTensorField LL = dev(filter_(sqr(U())) - (sqr(filter_(U()))));

    volSymmTensorField MM =
        sqr(delta())*(filter_(mag(D)*(D)) - 4*mag(filter_(D))*filter_(D));

    dimensionedScalar MMMM = average(magSqr(MM));

    if (MMMM.value() > VSMALL)
    {
        return average(LL && MM)/MMMM;
    }
    else
    {
        return 0.0;
    }
}

dimensionedScalar dynSmagorinsky::cI(const volSymmTensorField& D) const
{
    volScalarField KK = 0.5*(filter_(magSqr(U())) - magSqr(filter_(U())));

    volScalarField mm =
        sqr(delta())*(4*sqr(mag(filter_(D))) - filter_(sqr(mag(D))));

    dimensionedScalar mmmm = average(magSqr(mm));

    if (mmmm.value() > VSMALL)
    {
        return average(KK*mm)/mmmm;
    }
    else
    {
        return 0.0;
    }
}

Best,

[panda60]

Dear alberto ,
Thank you very much. Your explanation is very good. So if it means dynamic Smagorinsky will be no useful ? And the dynamic Smagorinsky result will not be better than normal Smagorinsky ? Because the coefficient changing with space is more better than changing with time, so
normal Smagorinsky+VanDriest will be more better than dynamic Smagorinsky+smooth.
Am I right ?
If we can use dynamic Smagorinsky+VanDriest ?


Are you one of OpenFOAM's developer ? I can't see your pHD paper in http://foamcfd.org/resources/theses.html.

[panda60]

This is Kataoka's formula.

[alberto]

Originally Posted by panda60

Dear alberto ,
Thank you very much. Your explanation is very good. So if it means dynamic Smagorinsky will be no useful ? And the dynamic Smagorinsky result will not be better than normal Smagorinsky ?

It is still better than the normal Smagorinsky, because it can adapt to the average flow conditions, meaning that if the flow globally has low-Re, the coefficient will adapt.

The averaging is introduced to prevent the values of the coefficient to become negative, which would lead to negative viscosities. This is usually in simple cases (channel flow, for example) averaging along the homogeneous direction. Clearly this is not possible in general, because you do not necessarily have a homogeneous direction, and you have to find some other solution.

Because the coefficient changing with space is more better than changing with time, so
normal Smagorinsky+VanDriest will be more better than dynamic Smagorinsky+smooth.
Am I right ?
If we can use dynamic Smagorinsky+VanDriest ?

From a theoretical point of view the dynamic procedure removes the need of damping functions, but this assumes your coefficient can change at least when approaching walls.

If you want to use the dynamic models, I think you should modify the code so that the coefficient can actaully be local and not averaged on the whole computational domain.

Are you one of OpenFOAM's developer ? I can't see your pHD paper in http://foamcfd.org/resources/theses.html.

No, I'm not a developer of OpenFOAM(r). I'm an old-time user and I do research in CFD.

Best,

[xiao]

Hi Alberto,

Thank you for your explanations and for the reference. I have read the paper in detail, and things start to becoming much clearer to me now.

I would recommend your explanations and your summary of the relevant literature to anyone who has doubts on the proper choice of numerical scheme for LES simulations. Hope they can find your post here.

Your discussion about the dynamic LES model is also very help. I will look at the code and reference in more detail, and get back to the thread later.

BTW, do you have any writings or publications related to OpenFoam? E.g. your thesis (master/phd), report, etc? About modifying the implementation of dynamic LES model as you mentioned in this thread, did you give it a try during your work, or did you know of anyone who may have an faithful implementation of the original idea? I guess they (the developers) must have had a reason to chose to do the spatial averaging (e.g. robustness considerations?).

Best,
Heng

Originally Posted by alberto

I think it is a tradeoff between stability and accuracy. In a finite volume code, your global accuracy is second order anyway. If you use higher order schemes for the convection term you can have some benefit, but you have also to evaluate how important this benefit compared to the standard linear scheme.

....

This problem does not appear in central schemes, which are characterized by real modified wave-numbers. The real part is responsible of the dispersion error, leading to energy pile-up, but the diffusion error is not present.

Best,

[alberto]

Originally Posted by xiao

BTW, do you have any writings or publications related to OpenFoam? E.g. your thesis (master/phd), report, etc?

My PhD work was done in part with OpenFOAM.

About modifying the implementation of dynamic LES model as you mentioned in this thread, did you give it a try during your work, or did you know of anyone who may have an faithful implementation of the original idea?

I'm not aware of that. To try it, you can create your own version of dynSmagorinsky, remove the global averaging and implement some strategy to avoid negative values of the coefficient itself.

I guess they (the developers) must have had a reason to chose to do the spatial averaging (e.g. robustness considerations?).

I think this is the reason. Probably a weaker limitation on the coefficient would be enough to make the code reasonably stable. You could try looking for how the model is implemented in general codes.

Best,

[xiao]

Hi Alberto,

Thanks very much for your responses! It is very help for me and for future FOAMers as well. A few follow-up questions:

1. If you have finished and published your PhD work, could you give me a pointer where we could find them?

2. If I understand correctly, the negative coefficient should be allowed, in order to account for "back-scatter" of the energy, however, it should be not allowed to be "too negative" so that the whole viscosity is negative, which is surely the recipe for disaster. Is that correct? Or the negative Cs should be forbidden overall?

3. Could you be more specific about in which code one may be able to find such an implementation? I guess one wouldn't be able to see the source code of commercial codes (fluent etc.).

Best,
Heng

Originally Posted by alberto

My PhD work was done in part with OpenFOAM.

I'm not aware of that. To try it, you can create your own version of dynSmagorinsky, remove the global averaging and implement some strategy to avoid negative values of the coefficient itself.

I think this is the reason. Probably a weaker limitation on the coefficient would be enough to make the code reasonably stable. You could try looking for how the model is implemented in general codes.

Best,

[alberto]

Hi, sorry for the late reply.

Originally Posted by xiao

1. If you have finished and published your PhD work, could you give me a pointer where we could find them?

Send me an email if you want my thesis. There is nothing in it about the implementation in OpenFOAM of the models however.

2. If I understand correctly, the negative coefficient should be allowed, in order to account for "back-scatter" of the energy, however, it should be not allowed to be "too negative" so that the whole viscosity is negative, which is surely the recipe for disaster. Is that correct? Or the negative Cs should be forbidden overall?

Yes, that's correct. You could limit the coefficient so that, if the coefficient is negative, the absolute value of the SGS viscosity cannot become the same or higher than the molecular viscosity, so that the global viscosity cannot be negative.
I guess the idea of keeping the coefficient positive is to assume the backscatter negligible and to just take the adaptive coefficient, so that it goes to zero when the flow becomes laminar.

3. Could you be more specific about in which code one may be able to find such an implementation? I guess one wouldn't be able to see the source code of commercial codes (fluent etc.).

I meant you could check their documentation and look for references, but probably they do not put anything specific. You might also want to search for other open-source projects that might implement them.

Best,

[foamWang]

Hi alberto,

Thank you for your explaination and I learnt from your posts.

In terms of your claim, I am not fully convinced.

In literature, I found some guy is using QUICK for LES, e.g.

1. Y. Zang, “On the Development of Tools for the Simulation of Geophysical Flows.” Ph.D Thesis, Stanford University (1995)
2. M. V. Salvetti, Y. Zang, R. L. Street, and S. Banerjee, “Large-eddy simulation of free-surface decaying turbulence with dynamic subgrid-scale models”, Phys. Fluids 9, 2405-19 (1997).
3. A. Cui and R. L. Street, “Large-Eddy Simulation of Coastal Upwelling Flow”, Environ. Fluid Mech. 4, 197-223 (2004).
4. Y. Zang and R. L. Street, “Numerical simulation of coastal upwelling and interfacial instability of a rotating and stratified fluid”, J. Fluid Mech. 305, 47-75 (1995).
5. L. L. Yuan, R. L. Street, and J. H. Ferziger, “Large-eddy simulations of a round jet in crossflow”, J. Fluid Mech. 379, 71-104 (1999).


They are using "Dynamic Mixed Model (DMM)".


I don't know whether the case is suitable to mention here. But please give comment.
BTW, I changed all the schemes back to Gauss linear with an improved grid. Now I get rid of the blow up problem using dynSmagorinsky and dynMixedSmagorinsky. My experience is that avoiding wedge cells and make the grid space increase as smooth as possible. Any sudden jump in grid space may lead to blow up.


I am comparing different LES model now and hope to get back with some posts.

Thanks.


roro

Originally Posted by alberto

QUICK is not a suitable scheme for LES, and the upwind scheme is even worse. They are too dissipative.

You should use linear, or higher order central schemes. Generally, avoid limited schemes, they're not a good choice if you want accurate LES.

If the solution "blows up" using the linear scheme, it is probably because your grid is not fine enough to have the local Pé < 2, as required by the scheme.

Best,

[alberto]

Originally Posted by foamWang

Hi alberto,

Thank you for your explaination and I learnt from your posts.

In terms of your claim, I am not fully convinced.

Those claims are extensively discussed in the literature, and well known to who does accurate turbulent simulations.

In literature, I found some guy is using QUICK for LES, e.g.

<cut>

They are using "Dynamic Mixed Model (DMM)".


I don't know whether the case is suitable to mention here. But please give comment.

Do they report a comparison with other schemes? I'll check the references.

BTW, I changed all the schemes back to Gauss linear with an improved grid. Now I get rid of the blow up problem using dynSmagorinsky and dynMixedSmagorinsky. My experience is that avoiding wedge cells and make the grid space increase as smooth as possible. Any sudden jump in grid space may lead to blow up.

Yes, that's kind of expected with central schemes

Best,

[alberto]

I think they used QUICK because they had oscillatory behaviour from central schemes, and wanted to remove them to obtain a smooth solution.

There are other strategies propose to do this (see, for example, the paper cited above), where, with the due care, upwinded schemes are developed for LES to find a solution to the important dissipative error.

Best,

[panda60]

Dear alberto,
I want to know where the coefficient is easy to become negative. If this usually happens near the wall, I think dynamic model+VanDriest damping can solve this problem.

I saw Baba-Ahmadi and G.Tabor's paper already use like this.

M.H.Baba-Ahmadi, G.Tabor. Inlet conditions for LES using mapping and feedback control. Computers&Fluids, 2009.

in page1303 of their paper, it wrote like this:
SGS modelling is provided by the dynamic one-equation model.
The mesh is generated from two blocks in the y direction, allowing mesh grading towards the walls coupled with VanDriest damping to deal with the near-wall flow.

Please correct me if I am wrong.

[xiao]

The best way to prove (or disprove) is to take a channel chase, and do a few tests. An implementation of the "correct" dynamic model (with Cs varying both spatially and temporally, as opposed to only temporally in OpenFOAM, pointed out by Alberto) should not be challenging, however, a comprehensive testing is. As I said in a previous post, the developers must have had some reason to implement the dynamic models this way. If it were easy, they should have implemented according to the original formulation.

I will have a student to investigate this issue at a later time as his "semester project" (after the Easter break). Certainly we will share the findings here. If anyone has any input, we would be happy to know.

Best,
Heng

Originally Posted by panda60

Dear alberto,
I want to know where the coefficient is easy to become negative. If this usually happens near the wall, I think dynamic model+VanDriest damping can solve this problem.

I saw Baba-Ahmadi and G.Tabor's paper already use like this.

M.H.Baba-Ahmadi, G.Tabor. Inlet conditions for LES using mapping and feedback control. Computers&Fluids, 2009.

in page1303 of their paper, it wrote like this:
SGS modelling is provided by the dynamic one-equation model.
The mesh is generated from two blocks in the y direction, allowing mesh grading towards the walls coupled with VanDriest damping to deal with the near-wall flow.

Please correct me if I am wrong.

[alberto]

Originally Posted by panda60

Dear alberto,
I saw Baba-Ahmadi and G.Tabor's paper already use like this.

M.H.Baba-Ahmadi, G.Tabor. Inlet conditions for LES using mapping and feedback control. Computers&Fluids, 2009.

in page1303 of their paper, it wrote like this:
SGS modelling is provided by the dynamic one-equation model.
The mesh is generated from two blocks in the y direction, allowing mesh grading towards the walls coupled with VanDriest damping to deal with the near-wall flow.

Please correct me if I am wrong.

The dynamic one-equation model is not the dynamic Smagorinsky model. In the one-equation model you have a transport equation for the SGS turbulent kinetic energy, and the coefficients are determined according to the dynamic procedure.

I would suggest you to refer to the original paper about the dynamic procedure, and then look in the literature (not necessarily OpenFOAM-related, for the points discussed above) to see comparisons of standard Smagorinsky and dynamic Smagorinsky.

Best,
Alberto

[alberto]

Originally Posted by xiao

The best way to prove (or disprove) is to take a channel chase, and do a few tests. An implementation of the "correct" dynamic model (with Cs varying both spatially and temporally, as opposed to only temporally in OpenFOAM, pointed out by Alberto) should not be challenging, however, a comprehensive testing is. As I said in a previous post, the developers must have had some reason to implement the dynamic models this way. If it were easy, they should have implemented according to the original formulation.

As a suggestion to save time: search the literature, the dynamic model has been used many times and compared to others and DNS.
Stanford CTR is surely a good source of information. For example

http://ctr.stanford.edu/Summer02/jeanmart.pdf

explains what they did to find the coefficient as a function of the distance from the wall.

Best,
Alberto

[jiejie]

Originally Posted by panda60

Dear foamers,

dynSmagorinskyCoeffs
{
filter simple;
ce 1.05;
}

Can I ask what is the "ce 1.05" in the brackets? Do I have to adjust it for different cases?

Thank you very much.