CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Convectional term -fvm::Sp(fvc::div(phi_), epsilon_) in RNGkEpsilon OF 1.7.1 missing (https://www.cfd-online.com/Forums/openfoam-solving/88770-convectional-term-fvm-sp-fvc-div-phi_-epsilon_-rngkepsilon-1-7-1-missing.html)

 makaveli_lcf May 25, 2011 11:56

Convectional term -fvm::Sp(fvc::div(phi_), epsilon_) in RNGkEpsilon OF 1.7.1 missing

Hi!

I have a quistion: why this terms are missing in RNGkEpsilon model in epsilon and k equations?

Code:

```        fvm::ddt(epsilon_)       + fvm::div(phi_, epsilon_)       - fvm::Sp(fvc::div(phi_), epsilon_)       - fvm::laplacian(DepsilonEff(), epsilon_) ...         fvm::ddt(k_)       + fvm::div(phi_, k_)       - fvm::Sp(fvc::div(phi_), k_)       - fvm::laplacian(DkEff(), k_)```
... and special thanks to Prof. Jasak, I'am using his updated version with
the

Code:

```      + fvm::SuSp(-fvc::div(phi_), epsilon_) ...       + fvm::SuSp(-fvc::div(phi_), k_)```
which helps to reduce error due to the non-orthgonality.
Only one remark: why there is no such treatment for kEpsilon model in OF-1.6-ext (in kEpsilon.C)?

BTW the residuals normalization used now in OF is not appropriated, because it is connected with the source term and for the turbulent transport equations produces under-estimated residuals level.

Regards,
Alexander

 alberto May 26, 2011 02:35

Quote:
 Originally Posted by makaveli_lcf (Post 309221) Hi! I have a quistion: why this terms are missing in RNGkEpsilon model in epsilon and k equations? Code: ```        fvm::ddt(epsilon_)       + fvm::div(phi_, epsilon_)       - fvm::Sp(fvc::div(phi_), epsilon_)       - fvm::laplacian(DepsilonEff(), epsilon_) ...         fvm::ddt(k_)       + fvm::div(phi_, k_)       - fvm::Sp(fvc::div(phi_), k_)       - fvm::laplacian(DkEff(), k_)```
The code currently in 1.7.x looks consistent with the standard form of the incompressible RNG k-epsilon model. Why should that term be there?

It seems you are trying to treat the convective term in non-conservative form (why?). If you think to:

div(kU) = U . grad(k) + k div(U),

you are discretizing U . grad(k), since you have the full convective term in conservative form at the line

fvm::div(phi_, epsilon_)

Best,

 makaveli_lcf May 26, 2011 05:14

... I ask because IT IS, for example, in realizableKE (incompressible) in OF 1.7.1:

Code:

```    // Dissipation equation     tmp<fvScalarMatrix> epsEqn     (         fvm::ddt(epsilon_)       + fvm::div(phi_, epsilon_)       - fvm::Sp(fvc::div(phi_), epsilon_)       - fvm::laplacian(DepsilonEff(), epsilon_)     ==         C1*magS*epsilon_       - fvm::Sp         (             C2_*epsilon_/(k_ + sqrt(nu()*epsilon_)),             epsilon_         )     ); ...     // Turbulent kinetic energy equation     tmp<fvScalarMatrix> kEqn     (         fvm::ddt(k_)       + fvm::div(phi_, k_)       - fvm::Sp(fvc::div(phi_), k_)       - fvm::laplacian(DkEff(), k_)     ==         G - fvm::Sp(epsilon_/k_, k_)     );```
as well as in standard k-epsilon, and also in kOmegaSST...

Heh..., it is funny, the flux phi is divergence free... So why is it included?

for incompressible liquid. That is well known transition between conservative and non-conservatives forms of the convectional term.

Alberto, I am sorry, but would you be kind to clarify: is convectional term in RNG somehow different from those in standard and realizable k-epsilon?
And are you tallking about 1.7.x, which is not 1.7.1? Do you mean some development version?

 makaveli_lcf May 26, 2011 05:36

So I agree with Alberto,

Hmm..., why? Is it some error?

For example, standard k-epsilon, transport of kinetic energy is:

dk/dt + div(U K) = div(D_k grad(k)) + G - epsilon

There is no terms, from which we can get U*grad(k).

To avoid miss-understanding:

I am not trying to put "div(U*k) - k*div(U) = U*grad(k)" in the OF source code, it is all already there)))

I only noticed the difference between approaches how convectional term is discretized in RNGkEpsilon.C and in other incompressible RAS models. I recognized it while modifying RAS models for my solidification part.

 alberto May 26, 2011 10:15

Hi,

Quote:
 Originally Posted by makaveli_lcf (Post 309314) So I agree with Alberto, that U*grad(k) or U*grad(epsilon) or U*grad(omega) are discretized. Hmm..., why? Is it some error? For example, standard k-epsilon, transport of kinetic energy is: dk/dt + div(U K) = div(D_k grad(k)) + G - epsilon There is no terms, from which we can get U*grad(k). Any comments?
well, as you noticed, div(U), which in OpenFOAM is div(phi), is zero.
As a consequence, they are implementing U . grad(k) = div(U, k) - k div(U), to avoid using U & grad(k), and use the conservative flux.

What is not clear to me is why not implementing directly div(phi, k), as in the original equation from the literature?

 makaveli_lcf May 26, 2011 11:04

Sometimes for the discretization (e.g. for the finite differences method) it is easier to use nonconservative form of the convectional term

But here we get the wrong term. So I will try direct conservative form and compare the results.

 alberto May 26, 2011 11:25

Quote:
 Originally Posted by makaveli_lcf (Post 309381) Sometimes for the discretization (e.g. for the finite differences method) it is easier to use nonconservative form of the convectional term div(U*k) = k*div(U)+U & grad(k)
Yes, and I do not like using the non-conservative form, if not strictly necessary. :D

Quote:
 But here we get the wrong term. So I will try direct conservative form and compare the results.
If you consider div(U) = 0, the terms they have are identical, since

if div(U) = 0. Now, since numerically it is not exactly zero, they probably enforce the "incompressible form" by using the term

U & grad(k) = div(U*k) - k*div(U)

This in OpenFOAM becomes

fmv::div(phi, k) - fvm::Sp(fvc::div(phi), k),

since using phi and not U does not introduce errors, being phi corrected to enforce conservation after solving for p.

 makaveli_lcf May 26, 2011 16:14

Alberto, thanks! I got it:

Code:

```U = U* + delta U* = U - delta U* - exact solution, delta - error vector div(U*) = 0                                              (1) div(U) = div(U*) + div(delta) = div(delta)        (2)```
So if we would get the exact solution, it should be:
Code:

`div(U*.k) = U*&grad(k)`
Now we estimate an calculation error:
Code:

```Err1 = div(U.k) - U*&grad(k)       = k.div(U*) + U*&grad(k) + k.div(delta) + delta&grad(k) - U*&grad(k)```
From (1) and (2)
Code:

`Err1 = k.div(delta) + delta&grad(k) = k.div(U) + delta&grad(k)`
So, performing
div(U.k) - k.div(U) we reduce the calculation error to the delta&grad(k)

)))

I am curious, does it worth to use such approach for other scalar equations? E.g. for heat transfer... It would be really nice to here some comment from developers, why they decided to use this method for the turbulence modeling?

 lakeat July 19, 2012 11:53

I have the same questions, I feel very confused:

1. Is using "fvm::SuSp(-fvc::div(phi_), epsilon_)" better using "- fvm::Sp(fvc::div(phi_), epsilon_)"?
2. For other transport equations, should we use this approach? Why it is commented out from RNG k-epsilon?

 makaveli_lcf July 20, 2012 03:24

Hi Daniel!

1. Using SuSp operator is to increase matrix diagonal influence. With Sp operator it goes always to the diagonal and can reduce diagonal elements which makes its solution more problematic for the linear solver.

2. If you check OpenFOAM-ext such approach is used in all RANS model (*reffering to the source code.... yap! it's there!).

I am just using it for RANS models and for energy equation as well. But I did not make any test if it is better for the scalar equation in general. Let's do it!
What are your suggestions for the benchmark? It should be some simple case...

 lakeat July 20, 2012 10:42

pitzDaily maybe or other BFS flow?

 lakeat July 20, 2012 11:06

Hmm, Two more questions,

1. When to use source term explicitly and when shall one use it implicitly using Sp or Susp?

2. In ke formulation, why epsilon is solved first instead of k?

Thanks

 All times are GMT -4. The time now is 15:25.