CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Calculating divDevReff (https://www.cfd-online.com/Forums/openfoam-solving/58214-calculating-divdevreff.html)

jposunz December 11, 2008 22:57

Hi All, I'm using simpleFoa
 
Hi All,

I'm using simpleFoam and turbFoam with k-epsilon and will be trying different turbulence models in the future. Before that though, I'm trying to understand the implementation exactly.

I've searched a lot and I don't understand why the calculation for turbulence->divDevReff(U) is:

divDevReff(U) =
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))

instead of:

divDevReff(U) =
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*fvc::grad(U)().T())

ie. why do we take the deviatoric?

I've read several papers including Hrvoje's "A tensorial approach to computational continuum mechanics using object-oriented techniques" and still haven't found an answer.

Is there another paper that explains the implementation?

Thanks very much
John

jposunz December 14, 2008 18:48

Anyone? Or should I have been
 
Anyone? Or should I have been able to find this in the literature somewhere already?

louisgag July 16, 2009 17:59

Hi John,

did you find the answer to that? If so I'd be curious to know what it is. Thanks!

-Louis

jet November 14, 2009 14:09

Hi,

have you already found an answer to your question? I tried to understand this equation, too. But I failed!

louisgag November 30, 2009 23:04

I think it has to do with the fact that

dev(fvc::grad(U)().T()))

is, at convergence, the same as

(fvc::grad(U)().T()))

for some fluids

-Louis

4xF December 2, 2009 12:39

fvc::div(nuEff()*dev(fvc::grad(U)().T())) and fvc::div(nuEff()*fvc::grad(U)().T()) are the same. The reason is that the missing term is the diagonal of the tensor fvc::grad(U)().T() which is the trace of fvc::grad(U)().T() which is div(phi). For an incompressible fluid, this value is strictly zero by virtue of mass conservation.

mchurchf December 22, 2009 17:45

John,

I don't know if anyone has answered your question yet, but I've come across the same thing, and I have questions about the incompressible implementation of divDevReff also. As you said, in the incompressible models, turbulence->divDevReff(U) is:

- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))

Shouldn't dev2 be used instead of dev? In the compressible models, that is the case as turbulence->divDevRhoReff is:

- fvm::laplacian(muEff(), U)
- fvc::div(muEff()*dev2(fvc::grad(U)().T()))

In src\OpenFOAM\primitives\Tensor\TensorI.H, dev and dev2 are defined as dev(A) = A-1/3*I*trace(A) and dev2(A) = A-2/3*I*trace(A).

The stress term in the incompressible momentum equation is as follows:
d/dx_j {nu_eff*[(du_i/dx_j + du_j/dx_i)-2/3*du_k/dx_k delta_ij]},

which is the same as:
d/dx_j {nu_eff*(du_i/dx_j+du_j/dx_i)} - 2/3*d/dx_i {nu_eff(du_k/dx_k)}

which is in vector notation:
laplacian(nu_eff*U)
+ div{nu_eff*transpose[grad(U)]}
- div{2/3*nu_eff*I*trace[transpose(grad(U))]}

and reduces to:
laplacian(nu_eff*U) + div{nu_eff*dev2(transpose[grad(U)])}

Notice that dev2 is here instead of dev, which agrees with the compressible coding. However, in the incompressible coding, dev is used, which seem incorrect.

Now, the du_k/dx_k = trace[transpose(grad(U))] should be equal to zero in incompressible flow due to the continuity equation, so it seems that it should not matter if dev or dev2 is used. However, in the PISO algorithm, velocity may be predicted from a field that does not satisfy continuity due to initial conditions or insufficient pressure correction, so it seems like it would be wise to keep du_k/dx_k.

I am wondering if dev2 is an addition to OpenFOAM 1.6 and the incompressible divDevReff (and divDevBeff for LES) didn't get updated to reflect the change. Can someone with the correct knowledge please help?

Thank you,

Matt






Quote:

Originally Posted by jposunz (Post 183644)
Hi All,

I'm using simpleFoam and turbFoam with k-epsilon and will be trying different turbulence models in the future. Before that though, I'm trying to understand the implementation exactly.

I've searched a lot and I don't understand why the calculation for turbulence->divDevReff(U) is:

divDevReff(U) =
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))

instead of:

divDevReff(U) =
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*fvc::grad(U)().T())

ie. why do we take the deviatoric?

I've read several papers including Hrvoje's "A tensorial approach to computational continuum mechanics using object-oriented techniques" and still haven't found an answer.

Is there another paper that explains the implementation?

Thanks very much
John


christju January 2, 2010 12:04

Hello,
In fact, I've come across the same thing and I found the same conclusion than you Matthew...
The reynolds stress tensor is defined as :
-rho u'_i u'_j = 2 mu_t S_ij - 2/3 rho k delta_ij

where k = (u'_i u'_i) /2 is the kinetic energy, and S_ij is the deviatoric part of the mean strain rate :
S_ij = 1/2 (du_i/dx_j + du_j/dx_i) - 1/3 d_uk/dx_k delta_ij

As you said Matthew, the stress term in the incompressible momentum equation is as follows:
d/dx_j {nu_eff*[(du_i/dx_j + du_j/dx_i)-2/3*du_k/dx_k delta_ij]},

And in fact, from my point of view, a 2 is missing in the implementation of turbulence->divDevReff(U) of the incompressible model....

My other question is :
Where the term "- 2/3 rho k delta_ij" of the Reynold stress is taken into account?
Is it through a modified pressure defined as P = p/rho+2/3k and solved in the pressure equation?
And then, what is writen in the output files? this modified pressure?

If you have any idea....
Thank you in advance for your help,
Julien

poplar August 25, 2010 01:51

Hi,mchurchf.

The complete form of div(StressTensor)is:
div(muEff()*dev(fvc::grad(U)())) + div(muEff()*dev(fvc::grad(U)().T()))
So -1/3*I*trace(A) is lncluded in both dev(grad(U)) and dev(grad(U).T()),

in fact,
laplacian(nuEff(), U) = div(muEff()*grad(U)) =div(muEff()*dev(fvc::grad(U))) + 1/3*I*trace(A)

as a result, one 1/3*I*trace(A) is addded in laplacian(nuEff(), U), so a minus one must be subtracted in the later term.

then,dev2(A) = A-2/3*I*trace(A).

du_k/dx_k in incompressible fluids may be removed. If one retains it for correction, 1/3 or 2/3 is only a scale factor.

Quote:

Originally Posted by mchurchf (Post 240765)
John,

I don't know if anyone has answered your question yet, but I've come across the same thing, and I have questions about the incompressible implementation of divDevReff also. As you said, in the incompressible models, turbulence->divDevReff(U) is:

- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))

Shouldn't dev2 be used instead of dev? In the compressible models, that is the case as turbulence->divDevRhoReff is:

- fvm::laplacian(muEff(), U)
- fvc::div(muEff()*dev2(fvc::grad(U)().T()))

In src\OpenFOAM\primitives\Tensor\TensorI.H, dev and dev2 are defined as dev(A) = A-1/3*I*trace(A) and dev2(A) = A-2/3*I*trace(A).

The stress term in the incompressible momentum equation is as follows:
d/dx_j {nu_eff*[(du_i/dx_j + du_j/dx_i)-2/3*du_k/dx_k delta_ij]},

which is the same as:
d/dx_j {nu_eff*(du_i/dx_j+du_j/dx_i)} - 2/3*d/dx_i {nu_eff(du_k/dx_k)}

which is in vector notation:
laplacian(nu_eff*U)
+ div{nu_eff*transpose[grad(U)]}
- div{2/3*nu_eff*I*trace[transpose(grad(U))]}

and reduces to:
laplacian(nu_eff*U) + div{nu_eff*dev2(transpose[grad(U)])}

Notice that dev2 is here instead of dev, which agrees with the compressible coding. However, in the incompressible coding, dev is used, which seem incorrect.

Now, the du_k/dx_k = trace[transpose(grad(U))] should be equal to zero in incompressible flow due to the continuity equation, so it seems that it should not matter if dev or dev2 is used. However, in the PISO algorithm, velocity may be predicted from a field that does not satisfy continuity due to initial conditions or insufficient pressure correction, so it seems like it would be wise to keep du_k/dx_k.

I am wondering if dev2 is an addition to OpenFOAM 1.6 and the incompressible divDevReff (and divDevBeff for LES) didn't get updated to reflect the change. Can someone with the correct knowledge please help?

Thank you,

Matt


T.D. October 13, 2010 13:54

new nuEff()
 
Hi guys
i need to add (implement) and use a new nuEff2(), how to do that?
where to add the new nuEff() ?

i tried but it's hard, since it was only one original "nu" defined in all the transportModels and viscosityModels, and is hard to change in all the files.

help please

haze_1986 October 26, 2012 10:14

I am curious about this too

haze_1986 January 25, 2013 05:27

Quote:

Originally Posted by poplar (Post 272612)
Hi,mchurchf.

The complete form of div(StressTensor)is:
div(muEff()*dev(fvc::grad(U)())) + div(muEff()*dev(fvc::grad(U)().T()))
So -1/3*I*trace(A) is lncluded in both dev(grad(U)) and dev(grad(U).T()),

in fact,
laplacian(nuEff(), U) = div(muEff()*grad(U)) =div(muEff()*dev(fvc::grad(U))) + 1/3*I*trace(A)

as a result, one 1/3*I*trace(A) is addded in laplacian(nuEff(), U), so a minus one must be subtracted in the later term.

then,dev2(A) = A-2/3*I*trace(A).

du_k/dx_k in incompressible fluids may be removed. If one retains it for correction, 1/3 or 2/3 is only a scale factor.

Hi, that sounds correct, but I do not understand why for the compressible case divDevRhoReff(U), dev2 is used instead?
http://www.cfd-online.com/Forums/ope...sonicfoam.html

MatzeS March 31, 2014 09:22

Hello,

although it's some time ago and several replies were posted, I still don't understand, if the implementation is correct...

Don't we have to use dev2() instead of dev() for the incompressible solvers (at least until convergence is reached)?
Or did I just miss the correct reply?

Regards,
Matthias

galap November 6, 2014 04:31

Quote:

Originally Posted by christju (Post 241306)
Hello,
In fact, I've come across the same thing and I found the same conclusion than you Matthew...
The reynolds stress tensor is defined as :
-rho u'_i u'_j = 2 mu_t S_ij - 2/3 rho k delta_ij

where k = (u'_i u'_i) /2 is the kinetic energy, and S_ij is the deviatoric part of the mean strain rate :
S_ij = 1/2 (du_i/dx_j + du_j/dx_i) - 1/3 d_uk/dx_k delta_ij

As you said Matthew, the stress term in the incompressible momentum equation is as follows:
d/dx_j {nu_eff*[(du_i/dx_j + du_j/dx_i)-2/3*du_k/dx_k delta_ij]},

And in fact, from my point of view, a 2 is missing in the implementation of turbulence->divDevReff(U) of the incompressible model....

My other question is :
Where the term "- 2/3 rho k delta_ij" of the Reynold stress is taken into account?
Is it through a modified pressure defined as P = p/rho+2/3k and solved in the pressure equation?
And then, what is writen in the output files? this modified pressure?

If you have any idea....
Thank you in advance for your help,
Julien

Exactly reprint of my doubts. The issue in the incompressible is not the problem for me - since I am dealing with compressible flows. There the implementation seems to be correct. But what I don't see:

Where the term "- 2/3 rho k delta_ij" of the Reynold stress is taken into account? Has anyone find the answer?

Tobi July 15, 2015 16:07

divDevRhoReff and divDevReff
 
Hi all,

as I read before, there are still some people who are confused about the equations and the behavior of the implemented code. Therefore, I want to make the derivations and want to show you that the implementation is correct for both cases. We consider compressible and incompressible fluids. Further more I only take into account laminar flow and explain the way to get to divDevRhoReff and divDevReff.

Note: A complete derivation can be found in my public book with more hints and comparison with the c++ code.


Momentum equation

Lets start with the obvious momentum equation:

\frac{\partial \rho \textbf{U}}{\partial t} + \nabla \bullet (\rho  \textbf{U} \textbf{U}) = \nabla \bullet \boldsymbol \sigma + \textbf{f}

Due to the fact that everyone are more or less familiar with the equations, I do not mention the variables meaning. For the further proceeding I do not consider body forces. If we have newtonian fluid the stress tensor \boldsymbol \sigma can be expressed as:

\boldsymbol \sigma = - \left(p + \frac{2}{3} \mu \nabla \bullet \textbf{U}\right) \textbf{I} + \mu \left[ \nabla \textbf{U} + (\nabla \textbf{U})^T\right]

Some modifications and we end up with:

\boldsymbol \sigma = - p\textbf{I} - \frac{2}{3} \mu \nabla \bullet  \textbf{U} \textbf{I} + \mu \nabla \textbf{U} + \mu  (\nabla  \textbf{U})^T

Combining both, and neglect all other body forces, we end up with:

\frac{\partial \rho \textbf{U}}{\partial t} + \nabla \bullet (\rho   \textbf{U} \textbf{U}) = \nabla \bullet \left[ - p\textbf{I} - \frac{2}{3} \mu \nabla \bullet  \textbf{U} \textbf{I} + \mu \nabla \textbf{U} + \mu  (\nabla  \textbf{U})^T\right]

Lets focus only only to the RHS (step by step):

\nabla \bullet  (- p\textbf{I}) + \nabla \bullet \left[-\frac{2}{3} \mu \nabla \bullet   \textbf{U} \textbf{I} + \mu \nabla \textbf{U} + \mu  (\nabla   \textbf{U})^T\right]

Now it can be shown that the first term is the gradient of the pressure:

\nabla \bullet  (- p\textbf{I}) = -\nabla p

The analytic proof:

\nabla \bullet  (- p\textbf{I}) = \nabla \bullet \left[ -p 
\begin{pmatrix}
1 & 0 &0 \\ 
0 & 1 &0 \\ 
 0 & 0 & 1
\end{pmatrix}
\right]
 =
\nabla \bullet
\begin{pmatrix}
-p & 0 &0 \\ 
0 & -p &0 \\ 
 0 & 0 & -p
\end{pmatrix}
= \begin{pmatrix}
 (-\frac{\partial p}{\partial x}) + 0 + 0\\ 
0 + (-\frac{\partial p}{\partial y}) + 0\\
0 + 0 + (- \frac{\partial p}{\partial x})
\end{pmatrix} =

=
\begin{pmatrix}
 (-\frac{\partial p}{\partial x})\\ 
(-\frac{\partial p}{\partial y})\\
(-\frac{\partial p}{\partial x})
\end{pmatrix}
=- \nabla p

The second term at the RHS is often symboled with \boldsymbol \tau and stand for the viscose part of the stress tensor (normally without the divergence symbol):

\nabla \bullet \boldsymbol \tau = \nabla \bullet \left[-\frac{2}{3} \mu \nabla \bullet   \textbf{U}  \textbf{I} + \mu \nabla \textbf{U} + \mu  (\nabla   \textbf{U})^T\right]

This guy is calculated by calling the function divDevRhoReff(U) or divDevReff(U). Now I will show that the implemented functions are correct.


Compressible

In compressible flow we have the full viscose stress tensor. Now I demonstrate that this term is exactly the divDevRhoReff(U) function. Lets start to modify the equation (again step by step). But first I want to introduce the necessary math operation. \textbf{a} denotes an arbitary vector:

\nabla \bullet
\textbf{a} = \mathrm{tr}(\nabla\textbf{a}) = \mathrm{tr}(\nabla \textbf{a})^T

using this formula and sorting the equation, we get the following:

\nabla \bullet \left[\mu \nabla \textbf{U} + \mu  \left((\nabla   \textbf{U})^T -\frac{2}{3} \mathrm{tr}(\nabla \textbf{U})^T   \textbf{I} \right)\right]

\nabla \bullet (\mu \nabla \textbf{U}) 
+ \nabla \bullet 
\left( \mu  \left[(\nabla   \textbf{U})^T - \frac{2}{3}  \mathrm{tr}(\nabla  \textbf{U})^T   \textbf{I} \right]\right)

Thats it. If we call turbulence->divDevRhoReff() we will exactly get this equation. Using a substituion for the gradient of the transponed velocity field, the name of the function dev2 gets clearer:

(\nabla   \textbf{U})^T = \textbf{A}

\nabla \bullet (\mu \nabla \textbf{U}) 
+ \nabla \bullet 
\left( \mu  \left[(\textbf{A} - \frac{2}{3}\mathrm{tr}(\textbf{A})   \textbf{I} \right]\right)

The deviatoric part of a matrix is defined as:

\mathrm{dev}(\textbf{A}) = \textbf{A} - \frac{1}{3}\mathrm{tr} (\textbf{A})\textbf{I}

Now we see, that the second term include the deviatoric part of the matrix (gradU). In compressible we call the dev2 due to the fact that the factor 2 is in the trace. To proove the equation, we can check out the code in the compressible tubulence file (laminar.C): http://foam.sourceforge.net/docs/cpp...ce.html#l00219

Code:

  219 tmp<fvVectorMatrix> laminar::divDevRhoReff(volVectorField& U) const   220 {
  221    return
  222    (
  223      - fvm::laplacian(muEff(), U)
  224      - fvc::div(muEff()*dev2(T(fvc::grad(U))))
  225    );
  226 }



\nabla \bullet (\mu \nabla \textbf{U})

is equal to (muEff = mu for laminar):
Code:

fvm::laplacian(muEff, U)
\nabla \bullet 
\left( \mu  \left[(\nabla   \textbf{U})^T - \frac{2}{3}  \mathrm{tr}(\nabla  \textbf{U})^T   \textbf{I} \right]\right)

is equal to
Code:

  224      - fvc::div(muEff()*dev2(T(fvc::grad(U))))
At least dev2 calculates exactly the deviatoric part. The code is here: http://foam.sourceforge.net/docs/cpp/a08631_source.html

Code:

  303 //- Return the deviatoric part of a symmetric tensor  304 template<class Cmpt>
   305 inline SymmTensor<Cmpt> dev2(const SymmTensor<Cmpt>& st)
  306 {
  307    return st - SphericalTensor<Cmpt>::twoThirdsI*tr(st);
  308 }

This is equal to:

\textbf{A} - \frac{2}{3}\mathrm{tr} (\textbf{A})\textbf{I} = \left[(\nabla   \textbf{U})^T - \frac{2}{3}  \mathrm{tr}(\nabla  \textbf{U})^T   \textbf{I} \right]

Check :) or no check :confused:, I think, its obvious.


Incompressible
For incompressible fluids (devide by \rho) we again start with the viscous stress tensor \boldsymbol \tau:

\nabla \bullet \boldsymbol \tau = \nabla \bullet \left[-\frac{2}{3} \nu  \underline{\nabla \bullet   \textbf{U}}  \textbf{I} + \nu \nabla \textbf{U} + \nu   (\nabla   \textbf{U})^T\right]

Hence, the density is constant the first term (underlined) at the RHS is zero due to the mass conservation equation:

\nabla\bullet \textbf{U} = 0

\nabla \bullet \boldsymbol \tau = \nabla \bullet \left[\nu \nabla \textbf{U} + \nu   (\nabla   \textbf{U})^T\right]

After a huge amount of mathematics we get :D

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet (\nu   (\nabla   \textbf{U})^T)

So there is no deviatoric part till now. But in FOAM we calculate the viscouse part of the stress tensor using divDevReff. Now we add the hydrostatic part multiply by the viscosity and substract it again:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu   (\nabla   \textbf{U})^T + \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I} -  \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right)

it follows:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[  (\nabla    \textbf{U})^T + \frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right] -   \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right)

The last term is zero due to mass conservation and can be removed:

-   \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I} =  -   \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})\textbf{I} =  -   \nu \nabla \bullet \textbf{U}

We end up with:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[  (\nabla    \textbf{U})^T + \frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right]\right)

If we substitute the transponed gradient U matrix with A:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[ A + \frac{1}{3}\mathrm{tr}(A)\textbf{I}\right]\right)

Again we have the laplacian term and again the deviatoric part (here calcated with oneThird, so the correct definition of dev).

And this is exactly what we get here:

Code:

  tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(nuEff(), U)
      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
    );
}

The code for the deviatoric part (dev()) is equal to dev2. The only difference is the factor 2/3 which is introduced by the divergence term as shown before.

I hope everything is clear now and that there are no mistakes.
Additionally I hope it makes stuff clearer and help anybody to get a better understanding how things work.

Kind regards


Thanks to Alexander Vakhrushev for all support!

ssherman July 16, 2015 10:10

Tobi,

Thanks for the clear and concise explanation. However, in the incompressible case, there are a few places where mass conservation can be invoked to set terms to zero, but are not. Why do we set the terms that we do to zero?

I ask because nonNewtonianIcoFoam does not compute the viscous terms through divDevReff. Instead it goes through a couple more steps (that appear valid to me) to get to an alternate expression:
Code:

- fvm::laplacian(fluid.nu(), U) - (fvc::grad(U) & fvc::grad(fluid.nu())
Is nonNewtonianIcoFoam 'wrong'? What is the motivation for the different terms?

In the thread below, I show that the two forms yield noticeably different results for a model problem on a highly non-Newtonian fluid.
http://www.cfd-online.com/Forums/ope...tml#post555712

Thanks,
Steve

Tobi July 24, 2015 03:17

Hi,

why we set some term (that could be set to zero) not to zero is the way I derivated the equations (:
If we will go straight forward there is no need to do this. See below.

Tobi August 6, 2015 04:23

Hi all,

the derivation is correct but not straight forward from the beginning. If I have time I will make a PDF where everything is straight forward. Then you can see why once we have 1/3 and the other time 2/3 and its not due to the fact that we extend a equation (like I did in the incompressible case). Of course, like I showed it is possible to derivate the 2/3 and 1/3 like I did but as I mentioned it is not straight forward.

galap August 6, 2015 04:48

I would highly appreciate it! Thank you - good work

syavash October 23, 2015 09:35

Quote:

Originally Posted by Tobi (Post 555624)
Hi all,

as I read before, there are still some people who are confused about the equations and the behavior of the implemented code. Therefore, I want to make the derivations and want to show you that the implementation is correct for both cases. We consider compressible and incompressible fluids. Further more I only take into account laminar flow and explain the way to get to divDevRhoReff and divDevReff.

Momentum equation

Lets start with the obvious momentum equation:

\frac{\partial \rho \textbf{U}}{\partial t} + \nabla \bullet (\rho  \textbf{U} \textbf{U}) = \nabla \bullet \boldsymbol \sigma + \textbf{f}

Due to the fact that everyone are more or less familiar with the equations, I do not mention the variables meaning. For the further proceeding I do not consider body forces. If we have newtonian fluid the stress tensor \boldsymbol \sigma can be expressed as:

\boldsymbol \sigma = - \left(p + \frac{2}{3} \mu \nabla \bullet \textbf{U}\right) \textbf{I} + \mu \left[ \nabla \textbf{U} + (\nabla \textbf{U})^T\right]

Some modifications and we end up with:

\boldsymbol \sigma = - p\textbf{I} - \frac{2}{3} \mu \nabla \bullet  \textbf{U} \textbf{I} + \mu \nabla \textbf{U} + \mu  (\nabla  \textbf{U})^T

Combining both, and neglect all other body forces, we end up with:

\frac{\partial \rho \textbf{U}}{\partial t} + \nabla \bullet (\rho   \textbf{U} \textbf{U}) = \nabla \bullet \left[ - p\textbf{I} - \frac{2}{3} \mu \nabla \bullet  \textbf{U} \textbf{I} + \mu \nabla \textbf{U} + \mu  (\nabla  \textbf{U})^T\right]

Lets focus only only to the RHS (step by step):

\nabla \bullet  (- p\textbf{I}) + \nabla \bullet \left[-\frac{2}{3} \mu \nabla \bullet   \textbf{U} \textbf{I} + \mu \nabla \textbf{U} + \mu  (\nabla   \textbf{U})^T\right]

Now it can be shown that the first term is the gradient of the pressure:

\nabla \bullet  (- p\textbf{I}) = -\nabla p

The analytic proof:

\nabla \bullet  (- p\textbf{I}) = \nabla \bullet \left[ -p 
\begin{pmatrix}
1 & 0 &0 \\ 
0 & 1 &0 \\ 
 0 & 0 & 1
\end{pmatrix}
\right]
 =
\nabla \bullet
\begin{pmatrix}
-p & 0 &0 \\ 
0 & -p &0 \\ 
 0 & 0 & -p
\end{pmatrix}
= \begin{pmatrix}
 (-\frac{\partial p}{\partial x}) + 0 + 0\\ 
0 + (-\frac{\partial p}{\partial y}) + 0\\
0 + 0 + (- \frac{\partial p}{\partial x})
\end{pmatrix} =

=
\begin{pmatrix}
 (-\frac{\partial p}{\partial x})\\ 
(-\frac{\partial p}{\partial y})\\
(-\frac{\partial p}{\partial x})
\end{pmatrix}
=- \nabla p

The second term at the RHS is often symboled with \boldsymbol \tau and stand for the viscose part of the stress tensor (normally without the divergence symbol):

\nabla \bullet \boldsymbol \tau = \nabla \bullet \left[-\frac{2}{3} \mu \nabla \bullet   \textbf{U}  \textbf{I} + \mu \nabla \textbf{U} + \mu  (\nabla   \textbf{U})^T\right]

This guy is calculated by calling the function divDevRhoReff(U) or divDevReff(U). Now I will show that the implemented functions are correct.


Compressible

In compressible flow we have the full viscose stress tensor. Now I demonstrate that this term is exactly the divDevRhoReff(U) function. Lets start to modify the equation (again step by step). But first I want to introduce the necessary math operation. \textbf{a} denotes an arbitary vector:

\nabla \bullet
\textbf{a} = \mathrm{tr}(\nabla\textbf{a}) = \mathrm{tr}(\nabla \textbf{a})^T

using this formula and sorting the equation, we get the following:

\nabla \bullet \left[\mu \nabla \textbf{U} + \mu  \left((\nabla   \textbf{U})^T -\frac{2}{3} \mathrm{tr}(\nabla \textbf{U})^T   \textbf{I} \right)\right]

\nabla \bullet (\mu \nabla \textbf{U}) 
+ \nabla \bullet 
\left( \mu  \left[(\nabla   \textbf{U})^T - \frac{2}{3}  \mathrm{tr}(\nabla  \textbf{U})^T   \textbf{I} \right]\right)

Thats it. If we call turbulence->divDevRhoReff() we will exactly get this equation. Using a substituion for the gradient of the transponed velocity field, the name of the function dev2 gets clearer:

(\nabla   \textbf{U})^T = \textbf{A}

\nabla \bullet (\mu \nabla \textbf{U}) 
+ \nabla \bullet 
\left( \mu  \left[(\textbf{A} - \frac{2}{3}\mathrm{tr}(\textbf{A})   \textbf{I} \right]\right)

The deviatoric part of a matrix is defined as:

\mathrm{dev}(\textbf{A}) = \textbf{A} - \frac{1}{3}\mathrm{tr} (\textbf{A})\textbf{I}

Now we see, that the second term include the deviatoric part of the matrix (gradU). In compressible we call the dev2 due to the fact that the factor 2 is in the trace. To proove the equation, we can check out the code in the compressible tubulence file (laminar.C): http://foam.sourceforge.net/docs/cpp...ce.html#l00219

Code:

  219 tmp<fvVectorMatrix> laminar::divDevRhoReff(volVectorField& U) const   220 {
  221    return
  222    (
  223      - fvm::laplacian(muEff(), U)
  224      - fvc::div(muEff()*dev2(T(fvc::grad(U))))
  225    );
  226 }



\nabla \bullet (\mu \nabla \textbf{U})

is equal to (muEff = mu for laminar):
Code:

fvm::laplacian(muEff, U)
\nabla \bullet 
\left( \mu  \left[(\nabla   \textbf{U})^T - \frac{2}{3}  \mathrm{tr}(\nabla  \textbf{U})^T   \textbf{I} \right]\right)

is equal to
Code:

  224      - fvc::div(muEff()*dev2(T(fvc::grad(U))))
At least dev2 calculates exactly the deviatoric part. The code is here: http://foam.sourceforge.net/docs/cpp/a08631_source.html

Code:

  303 //- Return the deviatoric part of a symmetric tensor  304 template<class Cmpt>
   305 inline SymmTensor<Cmpt> dev2(const SymmTensor<Cmpt>& st)
  306 {
  307    return st - SphericalTensor<Cmpt>::twoThirdsI*tr(st);
  308 }

This is equal to:

\textbf{A} - \frac{2}{3}\mathrm{tr} (\textbf{A})\textbf{I} = \left[(\nabla   \textbf{U})^T - \frac{2}{3}  \mathrm{tr}(\nabla  \textbf{U})^T   \textbf{I} \right]

Check :) or no check :confused:, I think, its obvious.


Incompressible
For incompressible fluids (devide by \rho) we again start with the viscous stress tensor \boldsymbol \tau:

\nabla \bullet \boldsymbol \tau = \nabla \bullet \left[-\frac{2}{3} \nu  \underline{\nabla \bullet   \textbf{U}}  \textbf{I} + \nu \nabla \textbf{U} + \nu   (\nabla   \textbf{U})^T\right]

Hence, the density is constant the first term (underlined) at the RHS is zero due to the mass conservation equation:

\nabla\bullet \textbf{U} = 0

\nabla \bullet \boldsymbol \tau = \nabla \bullet \left[\nu \nabla \textbf{U} + \nu   (\nabla   \textbf{U})^T\right]

After a huge amount of mathematics we get :D

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet (\nu   (\nabla   \textbf{U})^T)

So there is no deviatoric part till now. But in FOAM we calculate the viscouse part of the stress tensor using divDevReff. Now we add the hydrostatic part multiply by the viscosity and substract it again:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu   (\nabla   \textbf{U})^T + \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I} -  \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right)

it follows:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[  (\nabla    \textbf{U})^T + \frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right] -   \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right)

The last term is zero due to mass conservation and can be removed:

-   \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I} =  -   \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})\textbf{I} =  -   \nu \nabla \bullet \textbf{U}

We end up with:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[  (\nabla    \textbf{U})^T + \frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right]\right)

If we substitute the transponed gradient U matrix with A:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[ A + \frac{1}{3}\mathrm{tr}(A)\textbf{I}\right]\right)

Again we have the laplacian term and again the deviatoric part (here calcated with oneThird, so the correct definition of dev).

And this is exactly what we get here:

Code:

  tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(nuEff(), U)
      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
    );
}

The code for the deviatoric part (dev()) is equal to dev2. The only difference is the factor 2/3 which is introduced by the divergence term as shown before.

I hope everything is clear now and that there are no mistakes.
Additionally I hope it makes stuff clearer and help anybody to get a better understanding how things work.

Kind regards


Thanks to Alexander Vakhrushev for all support!

Dear Tobias,

I suppose that a deviatoric tensor should be trace-free, while you have added the trace to the original tensor in the following part:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[ A + \frac{1}{3}\mathrm{tr}(A)\textbf{I}\right]\right)


Please correct me if I am wrong!!

Thanks,
Syavash

Tobi November 24, 2015 11:37

Quote:

Originally Posted by syavash (Post 569921)
Dear Tobias,

I suppose that a deviatoric tensor should be trace-free, while you have added the trace to the original tensor in the following part:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[ A + \frac{1}{3}\mathrm{tr}(A)\textbf{I}\right]\right)


Please correct me if I am wrong!!

Thanks,
Syavash

You are correct. The deviatoric tensor is traceless but as you see out of the derivation this term will remain. I did more derivations and put all together. I think I am going to share this with you. Also divDevRhoReff and divDevReff are included.

agustinvo December 2, 2015 06:14

Quote:

Originally Posted by Tobi (Post 555624)
Incompressible
For incompressible fluids (devide by \rho) we again start with the viscous stress tensor \boldsymbol \tau:

\nabla \bullet \boldsymbol \tau = \nabla \bullet \left[-\frac{2}{3} \nu  \underline{\nabla \bullet   \textbf{U}}  \textbf{I} + \nu \nabla \textbf{U} + \nu   (\nabla   \textbf{U})^T\right]

Hence, the density is constant the first term (underlined) at the RHS is zero due to the mass conservation equation:

\nabla\bullet \textbf{U} = 0

\nabla \bullet \boldsymbol \tau = \nabla \bullet \left[\nu \nabla \textbf{U} + \nu   (\nabla   \textbf{U})^T\right]

After a huge amount of mathematics we get :D

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet (\nu   (\nabla   \textbf{U})^T)

So there is no deviatoric part till now. But in FOAM we calculate the viscouse part of the stress tensor using divDevReff. Now we add the hydrostatic part multiply by the viscosity and substract it again:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu   (\nabla   \textbf{U})^T + \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I} -  \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right)

it follows:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[  (\nabla    \textbf{U})^T + \frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right] -   \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right)

The last term is zero due to mass conservation and can be removed:

-   \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I} =  -   \nu\frac{1}{3}\mathrm{tr}(\nabla \textbf{U})\textbf{I} =  -   \nu \nabla \bullet \textbf{U}


We end up with:

\nabla \bullet (\nu \nabla \textbf{U}) + \nabla \bullet \left(\nu \left[  (\nabla    \textbf{U})^T + \frac{1}{3}\mathrm{tr}(\nabla \textbf{U})^T\textbf{I}\right]\right)

Hello

In the OpenFOAM programmer's guide, it is said that:

dev (A) = A - 1/3 * trace(A) * I

So I think the bold part I quoted is not correct at all: the other term should be removed instead the one you used, if I am not wrong. At the end, since it is an incompressible case, it should not be a problem during your derivation.

Tobi December 2, 2015 10:03

Again,

you are right. The deviatoric part is tracefree but I also added the trace because it is a term of numerical stabilisation. This will happen due to the fact of my derivation. Also you are right that the deviatoric part of a matrix is defined as:

dev\textbf{A} = \textbf{A} - \frac{1}{3} tr(\textbf{A}) \textbf{I}

But you mix up the way I derived the equation. Starting from the incompressible term; that means that the conti is zero and also the secondary viscosity \lambda = - \frac{2}{3}\mu is zero it follows for the shear-rate tensor:

\boldsymbol \tau = [\nu \nabla \textbf{U} + \nu(\nabla \textbf{U})^T]

to be more consistent it would be actually:

\boldsymbol \tau = [\nu \nabla \otimes \textbf{U} + \nu(\nabla \otimes \textbf{U})^T]

What I was doing now was adding the trace of gradient U and removed it again:

\boldsymbol \tau = [\nu \nabla \otimes \textbf{U} + \nu(\nabla \otimes \textbf{U})^T + \nu \frac{1}{3} tr (\nabla \textbf{U}) \textbf{I} - \nu \frac{1}{3} tr (\nabla \textbf{U}) \textbf{I} ]

It also can be shown that:

tr (\nabla \textbf{U}) = \nabla \bullet \textbf{U}

which is actually the continuity equation which is zero for incompressible fluids. I think so far we are fine. Due to the fact that the terms of tr (\nabla \textbf{U}) = 0 we remove the last term and keep the first one. That results in the euqation I showed you above.

Now you can also do it vice versa. Then you will end up with the negative trace which will actually lead to the defintion of the deviatoric part.

So it's only a sign question we are talking about at the moment. Due to the fact that the trace term has no influence if you are converged, we can neglect it. We can feel free to add or remove it.

Of course if we consider the deviatoric definition we will at least end up with the negative trace; this derivation I also did in the paper/summary I will publish. AND of course it has a more physical way to use the negative trace because it is at least the definition of the deviatoric part (:

Also I am not sure if the solution will change if we would add the trace part. At least to clean everything up. I should have used the negative trace in my derivation, because then we get the same stuff as in the FOAM code and we correspond with the definition of the deviatoric part.

I hope it is clear now!

The only thing that is wrong at the above equations is that the sign of the trace if you start from the deviatoric definition is different. Then you get a minus sign at the trace. BUT again this part can also be neglected due to the continuity equation and one may implement it without the trace. This term is only kept for numerical stabilization.

niko0807 December 4, 2015 05:40

Hi,

I have been trying to figure out the exact implementation of incompressible RANS equation based on this thread and this article

http://openfoamwiki.net/index.php/Bu...sinesqPisoFoam

Nevertheless, there are a few things I do not understand, and i think there are a few mistakes in the article.

Anyway, I will start by showing you my derivation from RANS to OpenFOAM:
__________________________________________________ ____________

The incompressible RANS equation is given by,

u_j\frac{\partial u_i}{\partial	x_j}=-\frac{1}{\rho}\frac{\partial p}{\partial x_i}+\frac{\partial}{\partial x_j}\left[\nu\left(\frac{\partial u_i}{\partial x_j	}+\frac{\partial u_j}{\partial x_i}\right)-\overline{u'_iu'_j} \right]

where the Boussinesq approximation states:

\overline{u'_iu'_j}= -\nu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right)  + \frac{2}{3}\overline{k} \delta_{ij}

Expressing the Reynolds stresses with the Boussinesq approximation gives cause for,

\bar{u}_j\frac{\partial\bar{u}_i}{\partial\bar{u}_j} =-\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} \left[(\nu+\nu_{t}) \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i}\right)\delta_{ij} -\frac{2}{3}k \delta_{ij}\right]

where the kinetic energy term, \frac{2}{3}k \delta_{ij}\, drops out.

What would seem to pose a paradox in connection with the incompressibility assumption is introduced. The isotropic part of the stress tensor is added and subtracted again. This entails:

\bar{u}_j\frac{\partial\bar{u}_i}{\partial\bar{u}_j} =-\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j}\nu_{Eff} \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i}+ \frac{1}{3}\frac{\partial \bar{u}_k}{\partial x_k}\delta_{ij} -\frac{1}{3}\frac{\partial \bar{u}_k}{\partial x_k}\delta_{ij} \right)

Due to incompressibility the positive isotropic part of the stress tensor can be canceled and an expansion of the equation gives rise to:

\bar{u}_j\frac{\partial\bar{u}_i}{\partial\bar{u}_j} =-\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j}\nu_{Eff} \frac{\partial \bar{u}_i}{\partial x_j}+ \frac{\partial}{\partial x_j} \nu_{Eff}  \left( \frac{\partial \bar{u}_j}{\partial x_i} - \frac{1}{3}\frac{\partial \bar{u}_k}{\partial x_k}\delta_{ij}\right)

This is the incompressible turbulent RANS equation for fluids where \nu and thus \nu_{Eff} cannot be taken outside the spatial derivative due to it being a variable as a consequence of large temperature variations (or other reasons e.g. non-newtonian).

In OpenFOAM the terms can be identified as,

\underbrace{u_j\frac{\partial u_i}{\partial u_j}}_{\text{div(phi, U)}} =-\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \underbrace{\underbrace{\frac{\partial}{\partial x_j}\nu_{Eff} \frac{\partial u_i}{\partial x_j}}_{\text{laplacian(nuEff(), U)}}+ \underbrace{\frac{\partial}{\partial x_j} \nu_{Eff}  \left( \frac{\partial u_j}{\partial x_i} - \frac{1}{3}\frac{\partial u_k}{\partial x_k}\delta_{ij}\right)}_{\text{div(nuEff()*dev(T(grad(U))))}}}_{\text{divDevReff(U)}}

which exactly what is given in the source code of the momentum equation,

Code:

    tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
    );

and the kOmegaSST turbulence model:

Code:

tmp<fvVectorMatrix> kOmegaSST::divDevReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(nuEff(), U)
      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
    );
}

The deviatoric is given by:

Code:

inline Tensor<Cmpt> dev(const Tensor<Cmpt>& t)
{
    return t - SphericalTensor<Cmpt>::oneThirdI*tr(t);
}

__________________________________________________ _____________

My questions:

- Why does the kinematic energy term drop out? Is it because the gradient of it is very small?

- Why is the isotropic part of the stress tensor added and subtracted again when it is zero due to incompressibility? So it is for numerical reasons, but which numerical reasons?

Tobi December 4, 2015 07:06

Quote:

Originally Posted by niko0807 (Post 576183)
Hi,

My questions:

- Why does the kinematic energy term drop out? Is it because the gradient of it is very small?

- Why is the isotropic part of the stress tensor added and subtracted again when it is zero due to incompressibility? So it is for numerical reasons, but which numerical reasons?


The reason is the following:

- I don't know it exactly without having a look into the code but it is easy to understand

- You have an estimation of the velocity field and hence you know the fluxes

- Then you calculate the pressure field due to the fluxes. Then you get new fluxes due to a new velocity field.

- NOW: for incompressible fluids \nabla \bullet \textbf{U} = 0 but will never be exactly zero due to numerics (imagine abort error)

- Therefore this term is not zero especially at the beginning of a calculation

- If we use this term we introduce somehow "additional flux" which then will lead to faster convergence


Till now I did not check your derivation but your kinetmatic energy term seems to be the secondary viscosity or first lames coefficient that is zero for incompressible fluids. You should start really from the beginning with the Cauchy stress tensor and start to derive the Navier-Stokes equations without any averaging and turbulence models. If you do so you get:

\frac{\mathrm{D}\rho\textbf{U}}{\mathrm{D}t} = \nabla \bullet \boldsymbol  \sigma + S_\mathrm{f}

where the last term represents other sources like buoyancy, gravitational acceleration and so on (neglecting now).

Then we know that we can split the Cauchy stress tensor to a hydrostatic and deviatoric part. At least the hydrostatic part is the negative pressure and the deviatoric part is for Newtonian fluids the shear rate tensor \boldsymbol\tau. We get:

\frac{\mathrm{D}\rho\textbf{U}}{\mathrm{D}t} = - \nabla p + \nabla \bullet \boldsymbol  \tau

It can be derived that the deviatoric part is actually:

\boldsymbol \tau = 2 \mu \underbrace{\frac{1}{2} [ \nabla \otimes \textbf{U} + (\nabla \otimes \textbf{U})^T}_{\mathrm{tensor~} \textbf{D}}] + \left\{(\lambda + \kappa) \nabla \bullet \textbf{U}\right\} \textbf{I}

D is also known as deformation rate tensor. Here we have the gradients of Velocity. If we are in solid mechanics, we get the gradient of displacements. The last term is due to compression and expansion effects. \kappa is named bulk viscosity that is zero for low gases and almost neglect able for dense gases and liquids. At least we end up with that equation:

\boldsymbol \tau = 2 \mu \textbf{D} + \left\{(\lambda) \nabla \bullet  \textbf{U}\right\} \textbf{I}

\lambda is the first Lamés coefficient, also known as secondary viscosity which is actually for Newtonian fluids \lambda = -\frac{2}{3} \mu. If we are in solid mechanics we have to replace \lambda with other equations. We end up:

\boldsymbol \tau = 2 \mu \textbf{D} - \left\{\frac{2}{3} \mu \nabla \bullet  \textbf{U}\right\}\textbf{I}

that is actually the same like:

\boldsymbol \tau = 2 \mu \textbf{D} - \left\{\frac{2}{3} \mu \mathrm{tr} \textbf{D}\right\} \textbf{I}

Now for incompressible fluids, the secondary viscosity is zero (no expansion and compression), also the divergence of velocity is zero (continuity). Therefore the last term is removed. If we also introduce the assumption of constant viscosity we will end up with a simple formulation for the shear-rate tensor and at least (laplace equation), if the shear-rate tensor is zero we get the famous Euler equation.

Now if we introduce some turbulence models, nothing will be changed in the shear rate tensor. Only the viscosity will increase due to the fact that that the assumption of the common turbulence models are based on the increase of viscosity.

Therefore I don't get the point why you have -\frac{2}{3}k in your equation.

niko0807 December 4, 2015 08:31

Hi Tobias,

Thank you for your thorough answer.

As I understand, you suggest that I have mixed up turbulent kinematic energy and bulk viscosity. I do not think that is the case. According to this reference

http://www.cfd-online.com/Wiki/Bouss...ity_assumption

The Boussinesq eddy assumption includes mean turbulent kinetic energy, and it has nothing to do with bulk viscosity.

I still do not understand why the mean turbulence kinetic energy drops out.

For consistency, I have now included the basics of the derivation from Cauchy's equation of motion to the implementation of the RANS kOmegaSST turbulence model in OpenFOAM:

----------------------------------------------------------------------------------------
Cauchy's equation of motion states:

\rho\bigg(\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial u_j}\bigg)=\rho g_i+\frac{\partial}{\partial x_j}\tau_{ij}

It relates the fluid-particle acceleration to the net body force (here just from gravity), \rho q_i, and the surface force, \partial\tau_{ij}/\partial x_j. The stress tensor, \tau_{ij}, is composed of a deviatoric part, \sigma_{ij}, for viscous shear stresses and an isotropic part, -p\delta_{ij}, for fluid static normal stresses:

\tau_{ij}=\sigma_{ij}-p\delta_{ij}

The viscous shear stress tensor is given by,

\sigma_{ij}=\left[2\mu S_{ij}+\bigg(\mu_v-\frac{2}{3}\mu\bigg)\frac{\partial u_k}{\partial x_k}\delta_{ij}\right]

and the strain-rate is:

S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)

Combining everything produces the NS in its full form:

\rho\bigg(\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial u_j}\bigg)=-\frac{\partial p}{\partial x_i}+\rho g_i+\frac{\partial}{\partial x_j}\bigg[\mu\bigg(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\bigg)+\bigg(\mu_v-\frac{2}{3}\mu\bigg)\frac{\partial u_k}{\partial x_k}\delta_{ij}\bigg]

Assuming incompressibility, steady state and neglecting gravity gives rise to:

u_j\frac{\partial u_i}{\partial	x_j}=-\frac{1}{\rho}\frac{\partial p}{\partial x_i}+\frac{\partial}{\partial x_j}\nu\left(\frac{\partial u_i}{\partial x_j	}+\frac{\partial u_j}{\partial x_i}\right)

Reynolds decomposition,

u_i=\overline{u_i}+u'_i \quad \quad	p=\overline{p}+p'

produces the RANS equation which is given by:

u_j\frac{\partial u_i}{\partial	x_j}=-\frac{1}{\rho}\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial}{\partial x_j}\left[\nu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\overline{u'_iu'_j} \right]

According to http://www.cfd-online.com/Wiki/Bouss...ity_assumption ,with the notation from the reference, the Boussinesq approximation for an incompressible flow states:

-\rho \overline{u'_i u'_j} = \mu_t \, \left( \frac{\partial U_i}{\partial x_j} + \frac{\partial U_j}{\partial x_i} - \frac{2}{3} \frac{\partial U_k}{\partial x_k} \delta_{ij} \right) - \frac{2}{3} \rho k \delta_{ij}

For an incompressible flow it is thus:

\overline{u'_iu'_j}= -\nu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right)  + \frac{2}{3}\overline{k} \delta_{ij}

Expressing the Reynolds stresses with the Boussinesq approximation gives cause for,

\bar{u}_j\frac{\partial\bar{u}_i}{\partial\bar{u}_j} =-\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} \left[(\nu+\nu_{t}) \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i}\right) -\frac{2}{3}k \delta_{ij}\right]

where the kinetic energy term, \frac{2}{3}k \delta_{ij}\, drops out.

What would seem to pose a paradox in connection with the incompressibility assumption is introduced. The isotropic part of the stress tensor is added and subtracted again. This entails:

\bar{u}_j\frac{\partial\bar{u}_i}{\partial\bar{u}_j} =-\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j}\nu_{Eff} \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i}+ \frac{1}{3}\frac{\partial \bar{u}_k}{\partial x_k}\delta_{ij} -\frac{1}{3}\frac{\partial \bar{u}_k}{\partial x_k}\delta_{ij} \right)

Due to incompressibility the negative isotropic part of the stress tensor can be canceled and an expansion of the equation gives rise to:

\bar{u}_j\frac{\partial\bar{u}_i}{\partial\bar{u}_j} =-\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j}\nu_{Eff} \frac{\partial \bar{u}_i}{\partial x_j}+ \frac{\partial}{\partial x_j} \nu_{Eff}  \left( \frac{\partial \bar{u}_j}{\partial x_i} + \frac{1}{3}\frac{\partial \bar{u}_k}{\partial x_k}\delta_{ij}\right)

This is the incompressible turbulent RANS equation for fluids where \nu and thus \nu_{Eff} cannot be taken outside the spatial derivative due to it being a variable as a consequence of large temperature variations (or other reasons e.g. non-newtonian).

In OpenFOAM the terms can be identified as,

\underbrace{u_j\frac{\partial u_i}{\partial u_j}}_{\text{div(phi, U)}} =-\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \underbrace{\underbrace{\frac{\partial}{\partial x_j}\nu_{Eff} \frac{\partial u_i}{\partial x_j}}_{\text{laplacian(nuEff(), U)}}+ \underbrace{\frac{\partial}{\partial x_j} \nu_{Eff}  \left( \frac{\partial u_j}{\partial x_i} + \frac{1}{3}\frac{\partial u_k}{\partial x_k}\delta_{ij}\right)}_{\text{div(nuEff()*dev(T(grad(U))))}}}_{\text{divDevReff(U)}}

which exactly what is given in the source code of the momentum equation,

Code:

    tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
    );

and the kOmegaSST turbulence model:

Code:

tmp<fvVectorMatrix> kOmegaSST::divDevReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(nuEff(), U)
      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
    );
}

The deviatoric is given by:

Code:

inline Tensor<Cmpt> dev(const Tensor<Cmpt>& t)
{
    return t - SphericalTensor<Cmpt>::oneThirdI*tr(t);
}


Tobi December 4, 2015 10:13

Hi,

okay as far as I read your equations your are right but it is interesting that everybody denotes everything with other symbols. In all literature that I read, the Cauchy stress tensor is denoted by \boldsymbol \sigma and the deviatoric part of this tensor (the shear-rate tensor) as \boldsymbol \tau. In your case it is exactly vice versa. Anyway ...

I just wanted to mention this for people who are not so familiar with the equations (if someone will read all the stuff) :D

Okay. Your derivation seems to be well done and I also found in my literature also the definition / approximation:

-\rho \overline{u_i' u_j'} = \mu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} +  \frac{\partial \bar{u}_j}{\partial x_x}\right) - \frac{2}{3}\rho\delta_{ij}k

But as far as I know the equations we get the following: The Navier-Stokes-Equation with the Reynolds-Average method is given as:

\frac{\partial \bar {u_i}}{\partial t} + \frac{\partial}{\partial x_j} \left(\rho \bar u_i \bar u_j + \rho \overline{u_i' u_j'} \right)
=
- \frac{\partial \bar p}{\partial x_i} + \frac{\partial  \bar \tau_{ij}}{\partial x_j}

Up to know I think we are at the same level. The stress-tensor \boldsymbol \tau is calculated as before... the only difference is that we take the mean value of the velocity and the effective viscosity. Therefore the calculation of the U field is okay for the shear-rate tensor.

Now if we check out the equation above and compare it with yours, we get the same but in your equation the fluctuation term is put to the right hand side and is combined with the stress tensor.

I think so far we are fine. The question now is, what happen to the term including the kinetic energy. Till now I dont know but I will check it out.

Maybe I am not correct but if I insert everything into the reynolds averaged equation we get:

LHS:

\frac{\partial \bar {u_i}}{\partial t} + \frac{\partial}{\partial x_j}  \left(\rho \bar u_i \bar u_j - \mu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} +  \frac{\partial  \bar{u}_j}{\partial x_x}\right) + \frac{2}{3}\rho\delta_{ij}k \right)

RHS:

=
- \frac{\partial \bar p}{\partial x_i} + \frac{\partial}{\partial x_j} \underbrace{\left[\mu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} +  \frac{\partial  \bar{u}_j}{\partial x_x}\right) - \frac{2}{3}\mu_t  \frac{\partial u_k}{\partial x_k} \delta_{ij} \right]}_{\mathrm{shear-rate~tensor}}

Then we can make some simple math to get at least your eqaution that is side sth like (RHS):

\frac{\partial}{\partial x_i} \left[2 \mu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} +  \frac{\partial   \bar{u}_j}{\partial x_x}\right)\right] + \frac{\partial}{\partial x_i} \left(-\frac{2}{3}\mu_t  \frac{\partial  u_k}{\partial x_k} \delta_{ij} -  \frac{2}{3}k \delta_{ij}\right)


I think you agree with me. If the second term is zero we actually get the implementation of icoFoam with the factor 2. Of course it is clear what will happen if we have incompressible flow. But what will happen with the term that include k, I don't know. But as we can see, in the calculation of the deviatoric part of the stress tensor there is no twoSymm(grad(U)) included and I think the terms that appear (its the first one in this post) is treated somewhere else. Maybe in the pressure equation somehow?

But I hope I could show you the way that the shear-rate tensor calculation (laplacian->turbulence...) has no influence of the reynolds averaging (RHS)


At least I would mention that you should use the turbulence viscosity for the shear-rate tensor (you used the laminar viscosity).

Tobi December 4, 2015 10:25

Hi,

I talked to my college and he told me that this term with k is not included (buoyancy turbulent) because it is so small and can be neglected.

niko0807 December 4, 2015 10:56

Hi,

Thanks! I expected k to be small and negligible.

Correct me if I am wrong, but I do not think you are right in connection with the viscosity.

The RANS equation expressed with the Boussinesq approximation should include the sum of the molecular (laminar) viscosity \nu and the turbulent eddy viscosity \nu_t i.e. \nu_{Eff}=\nu+\nu_t. Like this:

\frac{\partial}{\partial x_j} \left[(\nu+\nu_{t}) \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i}\right)\right]

In your derivation (where I have allowed myself to correct what I suppose is simply a misspelling):

\frac{\partial}{\partial x_i} \left[2 \mu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} +  \frac{\partial   \bar{u}_j}{\partial x_i}\right)\right]

you only have the turbulent eddy viscosity. That does not make sense to me, as diffusion would not be present in laminar regimes of a flow problem. Also

2\left( \frac{\partial \bar{u}_i}{\partial x_j} +  \frac{\partial   \bar{u}_j}{\partial x_x}\right)=4 S_{ij}

does not make sense, as it equals four times the strain rate instead of only two. Check the Cauchy's equation again.

Tobi December 4, 2015 12:22

Hi,

Of course you are right. In terms of buoyancy turbulence you should use the effective viscosity.

How you get 4 x S?

In my previous derivatives it would be better to use the effective viscosity instead of the turbulent. That is obvious. Unfortunately I am not available till Wednesday.

But I will replay after I read my books

niko0807 December 4, 2015 12:47

Hi,

The viscous shear stress tensor (the one I call \sigma_{ij} and you call \tau_{ij}) is given by,

\sigma_{ij}=\left[2\mu S_{ij}+\bigg(\mu_v-\frac{2}{3}\mu\bigg)\frac{\partial u_k}{\partial x_k}\delta_{ij}\right]

and the strain rate:

S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)

Combining,

\sigma_{ij}=\left[2\mu \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)+\bigg(\mu_v-\frac{2}{3}\mu\bigg)\frac{\partial u_k}{\partial x_k}\delta_{ij}\right]

produces:

\sigma_{ij}=\left[\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)+\bigg(\mu_v-\frac{2}{3}\mu\bigg)\frac{\partial u_k}{\partial x_k}\delta_{ij}\right]

You thus need to omit the 2 in this:

\frac{\partial}{\partial x_i} \left[2 \mu \left( \frac{\partial \bar{u}_i}{\partial x_j} +  \frac{\partial   \bar{u}_j}{\partial x_i}\right)\right]

Just to be clear, for the above everything is laminar.

About the turbulent viscosity:

I do not see why using the effective velocity is only necessary for buoyancy turbulence. The effective viscosity should always be used, right?

Tobi December 7, 2015 04:33

Hi Niko,

thanks for the replay but I know how the equations look like (:
You are right that the number "2" is shortened with the 1/2 from the strain-rate tensor. For the 2 in my derivation I have to check.

Maybe I will do it today.

For the viscosity. The best way is to use the effective viscosity but sometimes you will find that people only use the turbulent viscosity. For example, I used the flamelet model build by Alberto Cuoci et. al. in my masterthesis. This model is based on high turbulence and therefore laminar effects are really neglectible and they only used turbulence->mu(). However I prefer to use the effective viscosity.

Tobi December 7, 2015 07:07

1 Attachment(s)
Dear Nico,

you are right. I made a mistake with the viscosities ;).
Here are my derivation that I did now... I have to check it again these days but it seems that everything is correct.

Tobi December 9, 2015 02:22

Hi all,

I checked the derivation yesterday again and there should be a mistake with the sign. I will check it again today.

Regards,
Tobi

charitonas December 9, 2015 14:08

Hi All,

I am trying to implement a new stress tensor for compressible flows. instead of divDevRhoReff(U) i am trying to do:

volVectorField UK("UK", U*10);

and change the stress tensor to

divDevRhoReff(UK)

I just want to test how the simulation will react for higher velocity. I also tried the
volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(UK)))); and my simulation crash as well. Does anyone knows why it happens for high values of the velocity?

Thank you in advance

Tobi December 17, 2015 06:18

Niko I have a question to you:

\left(\mu_v - \frac{2}{3}\mu\right)

The first term in your equation is the volumetric viscosity also known as bulk viscosity that is neglected in most/all equations (Bird et. al). I ask because in solid mechanics we have the deviatoric part of the stress tensor given by:

\mathrm{dev}(\boldsymbol \sigma) = \boldsymbol \tau = 2\mu \textbf{D} + \lambda \mathrm{tr}(\textbf{D})\textbf{I}

If we talk about a newtonian fluid we have:

\mathrm{dev}(\boldsymbol \sigma) = \boldsymbol \tau = 2\mu \textbf{D} + \left(-\frac{2}{3} + \mu_v\right) \mathrm{tr}(\textbf{D})\textbf{I}

Now I am not sure if I can say that:

\lambda = \left(-\frac{2}{3}\mu + \mu_v\right)

or maybe:

\lambda = -\frac{2}{3}\mu

And the bulk viscosity is zero for solids. I read that \lambda is also known as dilitational viscosity or in solid mechanics first Lámes coefficent. If its so, I think the second one is valid. Hope you get the point.
Again... I am asking if we can generalize the shear-rate tensor (for solids and fluids) as:

\mathrm{dev}(\boldsymbol \sigma) = \boldsymbol \tau = 2\mu \textbf{D} +  \left(\lambda + \mu_v\right) \mathrm{tr}(\textbf{D})\textbf{I}

where \mu_v is zero for solids and neglectable for fluids and \lambda = -\frac{2}{3}\mu if we talk about netwonian fluids or is it even better or are we allowed to define it like:

\mathrm{dev}(\boldsymbol \sigma) = \boldsymbol \tau = 2\mu \textbf{D} + \lambda \mathrm{tr}(\textbf{D})\textbf{I}

where \lambda is defined for newtonian liquids as:

\lambda = \left(-\frac{2}{3} + \mu_v\right)

Hope it is clear (:


@charitonas:

If you think about conservation of equations you make a big mistake!

Tobi January 12, 2016 04:13

Hi all,

now the derivation are done clearly and smooth:

http://www.cfd-online.com/Forums/ope...tml#post580452

GerhardHolzinger January 12, 2016 04:32

Quote:

Originally Posted by niko0807 (Post 576201)
Hi Tobias,

For an incompressible flow it is thus:

\overline{u'_iu'_j}= -\nu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right)  + \frac{2}{3}\overline{k} \delta_{ij}

Expressing the Reynolds stresses with the Boussinesq approximation gives cause for,

\bar{u}_j\frac{\partial\bar{u}_i}{\partial\bar{u}_j} =-\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} \left[(\nu+\nu_{t}) \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i}\right) -\frac{2}{3}k \delta_{ij}\right]

where the kinetic energy term, \frac{2}{3}k \delta_{ij}\, drops out.

The last term ( \frac{2}{3}k \delta_{ij}\ ) is not neglected.

The pressure and the isotropic part of the Reynolds stress tensor can be absorbed into a modified mean pressure p^\prime [208].

As Pope states on page 88:

Quote:

[...] showing that the isotropic component (\frac{2}{3}k) can be absorbed in a modified mean pressure.
Thus:

\bar{u}_j\frac{\partial\bar{u}_i}{\partial\bar{u}_j}  =
-\frac{1}{\rho}\frac{\partial}{\partial x_i} \left[p + \frac{2}{3}k \delta_{ij} \right]
+  \frac{\partial}{\partial x_j} \left[(\nu+\nu_{t}) \left( \frac{\partial  \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial  x_i}\right) \right]

with

\delta_{ij} \neq 0 \quad \text{ for } i=j

we finally get

\bar{u}_j\frac{\partial\bar{u}_i}{\partial\bar{u}_j}  =
-\frac{1}{\rho}\frac{\partial}{\partial x_i} \underbrace{\left[p + \frac{2}{3}k \right]}_{p^\prime}
+  \frac{\partial}{\partial x_j} \left[(\nu+\nu_{t}) \left(  \frac{\partial  \bar{u}_i}{\partial x_j} + \frac{\partial  \bar{u}_j}{\partial  x_i}\right) \right]

References

[208] S. P. Pope. Turbulent Flows. Cambridge University Press, 2000.

Tobi January 12, 2016 04:42

Thanks for the interesting hint. That means that we calculate p' in OpenFOAM? Well I have to check it out in the literature.

Good to know. So at least the term exist in FOAM?

GerhardHolzinger January 12, 2016 04:55

We are solving for all kinds of stuff.
If you apply Reynolds-averaging to your governing equations and apply the Boussinesq hypothesis, then you solve for the Reynolds averaged flow quantities.
The pressure is the sum of the actual, Reynolds averaged pressure and the isotropic part of the Reynolds stress tensor.
The actual pressure is also the isotropic part of a general fluid stress tensor. However, for probably historical reasons, we separate the pressure and the (shear) stresses. In solid mechanics there is no such distinction, at least up to the point I got in contact with solid mechanics and FE simulation.

If you apply spatial filtering to derive the governing equations with LES turbulence modelling, then you end up with equations for the grid-scale flow quantities. The sub-grid stress tensor in LES modelling is analogous to the Reynolds stress tensor in RANS modelling.
Also in LES modelling, the pressure you solve for is the sum of the actual pressure and the isotropic part of the sub-grid stress tensor.

The beauty of this derivations is, although, we get equations for different things (temporal or spatially filtered flow quantities), the equations have the same mathematical structure. This allows the makers of CFD simulation software to implement a generic handling of turbulence.


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