Calculating divDevReff

jposunz · December 11, 2008, 22:57

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

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

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

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

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

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

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

, 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

$\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

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

, 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

$\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

December 11, 2008, 22:57	Hi All, I'm using simpleFoa	#1
jposunz New Member John O\'Sullivan Join Date: Mar 2009 Location: Auckland, New Zealand Posts: 7 Rep Power: 17	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

December 14, 2008, 18:48	Anyone? Or should I have been	#2
jposunz New Member John O\'Sullivan Join Date: Mar 2009 Location: Auckland, New Zealand Posts: 7 Rep Power: 17	Anyone? Or should I have been able to find this in the literature somewhere already?

December 2, 2009, 12:39		#6
4xF New Member Frank Albina Join Date: Mar 2009 Location: Switzerland Posts: 14 Rep Power: 17	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. Tushar@cfd, sharonyue, shang and 5 others like this.

January 2, 2010, 12:04		#8
christju New Member Julien Christophe Join Date: Nov 2009 Posts: 4 Rep Power: 17	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/3du_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 galap, kiddmax, dongchao yang and 2 others like this.

October 13, 2010, 13:54	new nuEff()	#10
T.D. Senior Member Join Date: Sep 2010 Posts: 226 Rep Power: 17	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

July 16, 2009, 17:59		#3
louisgag Senior Member Louis Gagnon Join Date: Mar 2009 Location: Stuttgart, Germany Posts: 338 Rep Power: 18	Hi John, did you find the answer to that? If so I'd be curious to know what it is. Thanks! -Louis

November 14, 2009, 14:09		#4
jet New Member Join Date: Nov 2009 Posts: 12 Rep Power: 17	Hi, have you already found an answer to your question? I tried to understand this equation, too. But I failed!

November 30, 2009, 23:04		#5
louisgag Senior Member Louis Gagnon Join Date: Mar 2009 Location: Stuttgart, Germany Posts: 338 Rep Power: 18	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

October 26, 2012, 10:14		#11
haze_1986 Senior Member Join Date: Jul 2011 Posts: 120 Rep Power: 15	I am curious about this too Last edited by haze_1986; October 27, 2012 at 08:35.

March 31, 2014, 09:22		#13
MatzeS New Member Matthias Stammen Join Date: Oct 2010 Posts: 8 Rep Power: 16	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

July 15, 2015, 16:07	divDevRhoReff and divDevReff	#15
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	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>::twoThirdsItr(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 , 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 $\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! makaveli_lcf, olivierG, salehi144 and 80 others like this. __________________ Keep foaming, Tobias Holzmann Last edited by Tobi; April 27, 2017 at 05:37. Reason: mu = nu in incompressible; Sign mistake and renamed phi to a because phi = scalar normally*

July 16, 2015, 10:10		#16
ssherman New Member Join Date: Feb 2015 Posts: 2 Rep Power: 0	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

July 24, 2015, 03:17		#17
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	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. __________________ Keep foaming, Tobias Holzmann Last edited by Tobi; August 6, 2015 at 05:05. Reason: Reject the formulation...

August 6, 2015, 04:23		#18
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	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. __________________ Keep foaming, Tobias Holzmann

August 6, 2015, 04:48		#19
galap Member Fabian E. Join Date: Nov 2009 Posts: 38 Rep Power: 17	I would highly appreciate it! Thank you - good work Tobi likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Calculating Vorticity	David	CFX	28	October 29, 2017 19:02
Calculating Drag	Ajay Rao	FLUENT	8	February 15, 2010 10:15
what is Fluent calculating?	tomek	FLUENT	1	July 24, 2006 18:52
errors in calculating	ustcer	FLUENT	1	April 4, 2004 14:08
Calculating Coefficients	shabah	CFX	2	June 19, 2001 00:59