CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Favre or Reynolds Average in OpenFOAM?

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree9Likes
  • 6 Post By patrex
  • 3 Post By LuckyTran

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 20, 2017, 05:52
Default Favre or Reynolds Average in OpenFOAM?
  #1
Member
 
NablaDyn's Avatar
 
Join Date: Oct 2015
Location: Germany
Posts: 84
Rep Power: 6
NablaDyn is on a distinguished road
Dear FOAMers,

I'm a bit confused. I cannot determine exactly whether OpenFOAM uses Reynolds (RANS) or Favre averaging (FANS) when performing compressible 'RAS' simulations. The OpenFOAM terminology clearly implies it's RANS, so the turbulent fluctuations are not density-normalised. This assumption proves reasonable for many compressible flows (and thus would fit my needs) but not unconditionally for supersonic flows or such including combustion etc. Instead, rationality suggests OpenFOAM solves FANS since (to my knowlegde) its 'RAS' models are applicable to incompressible as well as compressible cases with no need for any user interaction. Or does OpenFOAM choose the best suited approach by itself based on solver conditioning (fvSolution) or else?

I searched the forum and the source but did not find a clear answer.

Any help is appreciated.

Regards and thanks in advance
NablaDyn is offline   Reply With Quote

Old   November 24, 2017, 15:27
Default Favre or Reynolds Average in OpenFOAM?
  #2
New Member
 
Patrick H.
Join Date: Aug 2017
Location: Germany
Posts: 3
Rep Power: 4
patrex is on a distinguished road
Hello NablaDyn,

This is a question wich i ask myself a lot in the last time. So I am also really interested in the 'right' answer and i already think i did some progress in it.



Well I think the problem with the compressible Reynolds-averaged-Navier-Stokes-Equations is, if they are averaging is done with the time averaging:

\overline{\phi}(x)=\lim\limits_{\Delta t\rightarrow\infty}\frac{1}{\Delta t}\int_T^{T+\Delta t} \phi(x,t) \,dt

and the splitting into:

\phi=\overline{\phi}+\phi',

the resulting equations are very unhandy because of the products of fluctuations between density and othe variables like the velocity (c.f. [1], p.91). This is shown as example in [2]
for the continunity equation in index notation (c.f. [2], p.100):

\frac{\partial \overline{\varrho}}{\partial t} + \frac{\partial}{\partial x_i}(\overline{\varrho}\,\overline{u}_i+\overline{\varrho'u_i'})=0.

The equations for the momentum and e.g. the total energy will much more complex. Beacause of this a density-weighted averaged is introduced (c.f. e.g. [1], p.91):

\widetilde{\phi}=\frac{\overline{\varrho\phi}}{\overline{\varrho}}

and

\phi=\widetilde{\phi}+\phi''.


This leads to the following equations without sources and body forces in index notation (c.f. [3], Equations (12), (13) and (14)):

Continunity equation:

\frac{\partial \overline{\varrho}}{\partial t} +\frac{\partial}{\partial x_i}\left[ \overline{\varrho} \widetilde{u_i} \right] = 0

Momentum equation:

\frac{\partial}{\partial t}\left( \overline{\varrho} \widetilde{u_i} \right) +\frac{\partial}{\partial x_j}\left[\overline{\rho} \widetilde{u_i} \widetilde{u_j} + \overline{p} \delta_{ij} + \overline{\rho u''_i u''_j} - \overline{\tau_{ji}}\right]=0

Total energy:

\frac{\partial}{\partial t}\left( \overline{\varrho} \widetilde{e_0} \right) +\frac{\partial}{\partial x_j}\left[ \overline{\varrho} \widetilde{u_j} \widetilde{e_0} + \widetilde{u_j} \overline{p} + \overline{u''_j p} + \overline{\varrho u''_j e''_0} + \overline{q_j} - \overline{u_i \tau_{ij}}\right]= 0.

(You can also find these equations in [1], p.92 (without total energy equation) and in [4], p.176)

These equations are called Favre-averaged Navier-Stokes equations but are also termed as compressible Reynolds-averaged Navier-Stokes equations (c.f. [5]).

As shown in [3], after some assumptions the equations can be written as:

\frac{\partial \overline{\varrho}}{\partial t} +\frac{\partial}{\partial x_i}\left[ \overline{\varrho} \widetilde{u_i} \right] = 0

\frac{\partial}{\partial t}\left( \overline{\varrho} \widetilde{u_i} \right) +\frac{\partial}{\partial x_j}\left[\overline{\varrho} \widetilde{u_j} \widetilde{u_i}+ \overline{p} \delta_{ij}- \widetilde{\tau_{ji}^{tot}}\right]= 0

\frac{\partial}{\partial t}\left( \overline{\varrho} \widetilde{e_0} \right) +\frac{\partial}{\partial x_j}\left[ \overline{\varrho} \widetilde{u_j} \widetilde{e_0} + \widetilde{u_j} \overline{p} + \widetilde{q_j^{tot}} - \widetilde{u_i} \widetilde{\tau_{ij}^{tot}}\right] = 0,

where

\widetilde{\tau_{ij}^{tot}} \equiv \widetilde{\tau_{ij}^{lam}} + \widetilde{\tau_{ij}^{turb}}

\widetilde{\tau_{ij}^{lam}} \equiv \widetilde{\tau_{ij}} = \mu\left( \frac{\partial \widetilde{u_i} }{\partial x_j} + \frac{\partial \widetilde{u_j} }{\partial x_i} - \frac{2}{3} \frac{\partial \widetilde{u_k} }{\partial x_k} \delta_{ij}\right)

\widetilde{\tau_{ij}^{turb}} \equiv- \overline{\rho u''_i u''_j} \approx\mu_t\left( \frac{\partial \widetilde{u_i} }{\partial x_j} + \frac{\partial \widetilde{u_j} }{\partial x_i} - \frac{2}{3} \frac{\partial \widetilde{u_k} }{\partial x_k} \delta_{ij}\right) -\frac{2}{3} \overline{\rho} k \delta_{ij}.

In tensor notation it should be:

\frac{\partial\overline{\varrho}}{\partial t}+\nabla\bullet(\overline{\varrho}\widetilde{\mathbf{u}}) = 0
(1)

\frac{\partial}{\partial t}(\overline{\varrho}\widetilde{\mathbf{u}})+\nabla\bullet\left(\overline{\varrho}\widetilde{\mathbf{u}}\otimes\widetilde{\mathbf{u}}+\overline{p}\mathbf{I}-\widetilde{\mathbf{\tau^{tot}}}\right) = 0
(2)

\frac{\partial}{\partial t} (\overline{\varrho}\widetilde{e_0})+\nabla\bullet\left(\overline{\varrho}\widetilde{e_0}\widetilde{\mathbf{u}}+\widetilde{\mathbf{u}}\overline{p}+\widetilde{\mathbf{q}^{tot}}-\widetilde{\mathbf{\tau^{tot}}}\bullet\widetilde{\mathbf{u}}\right) = 0
(3)

\widetilde{\mathbf{\tau^{tot}}}=\widetilde{\mathbf{\tau^{lam}}}+\widetilde{\mathbf{\tau^{turb}}}
(4)

\widetilde{\mathbf{\tau^{lam}}}=2\mu\widetilde{\mathbf{D}}-\frac{2}{3}\mu(\nabla\bullet\widetilde{\mathbf{u}})\mathbf{I}
(5)

\widetilde{\mathbf{\tau^{turb}}}=2\mu_t\widetilde{\mathbf{D}}-\frac{2}{3}\mu_t(\nabla\bullet\widetilde{\mathbf{u}})\mathbf{I}-\frac{2}{3}\overline{\varrho}\mathbf{I}k.
(6)


Im looking on sonicFoam of OpenFOAM 4.1, which is a solver forr trans-sonic/supersonic, turbulent flow and there is first the continunity equation solved in rhoEqn.H (code-snippet of the equation):
Code:
     fvScalarMatrix rhoEqn
     (
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
         fvOptions(rho)
     );
which i think could be translated into:

\frac{\partial\varrho}{\partial t}+\nabla\bullet(\varrho \mathbf{u}) = 0
(7)

after that UEqn.H is called (code-snippet of the equation):

Code:
 fvVectorMatrix UEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
   + turbulence->divDevRhoReff(U)
  ==
     fvOptions(rho, U)
 );
.
.
.
 if (pimple.momentumPredictor())
 {
     solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);
 }

which i think could be translated into:

\frac{\partial}{\partial t}(\varrho\mathbf{u})+\nabla\bullet(\varrho\mathbf{u}\otimes\mathbf{u})=-\nabla p^*+\nabla\bullet\mathbf{\tau^{eff}}=-\nabla\bullet(p^*\mathbf{I})+\nabla\bullet\mathbf{\tau^{eff}}
(8)

(the +\nabla\bullet\mathbf{\tau^{eff}} is due to the negativ implementation of the terms in the function divDevRhoReff())

where (see [2], cap. 10.2):

\mathbf{\tau^{eff}}=2\mu_\text{eff}\left[\frac{1}{2}\left\{(\nabla\otimes\mathbf{u})+(\nabla\otimes\mathbf{u})^T\right\}\right]-\frac{2}{3}\mu_\text{eff}(\nabla\bullet\mathbf{u})\mathbf{I}
(9)
=2\mu\mathbf{D}-\frac{2}{3}\mu(\nabla\bullet\mathbf{u})\mathbf{I} +2\mu_t\mathbf{D}-\frac{2}{3}\mu_t(\nabla\bullet\mathbf{u})\mathbf{I}

and (c.f. [6]):

p^*=p+\frac{2}{3}\varrho k
(10)

if we use (10) with (8):

\frac{\partial}{\partial t}(\varrho\mathbf{u})+\nabla\bullet(\varrho\mathbf{u}\otimes\mathbf{u})=-\nabla\bullet((p+\frac{2}{3}\varrho k)\mathbf{I})+\nabla\bullet\mathbf{\tau^{eff}}
(11)

and sort it:

\frac{\partial}{\partial t}(\varrho\mathbf{u})+\nabla\bullet(\varrho\mathbf{u}\otimes\mathbf{u}+p\mathbf{I}+\frac{2}{3}\varrho k\mathbf{I}-\mathbf{\tau^{eff}})=0
(12)


and after that EEqn.H is called (code-snippet of the equation):

Code:
    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, e) + fvm::div(phi, e)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + fvc::div(fvc::absolute(phi/fvc::interpolate(rho), U), p, "div(phiv,p)")
      - fvm::laplacian(turbulence->alphaEff(), e)
     ==
        fvOptions(rho, e)
    );
which i think could be translated into:

\frac{\partial}{\partial t}\left(\varrho e+\frac{1}{2}\varrho|\mathbf{u}|^2\right)=-\nabla\bullet\left[\varrho\mathbf{u}\left(e+\frac{1}{2}|\mathbf{u}|^2\right)\right]+\nabla\bullet\left(\alpha_\text{eff}\nabla e\right)-\nabla\bullet(p^*\mathbf{u})
(13)

which is in coincidence with [7] if e_0= e+\frac{1}{2}|\mathbf{u}|^2 is used:

\frac{\partial}{\partial t}\left(\varrho e_0\right)=-\nabla\bullet\left(\varrho\mathbf{u} e_0\right)+\nabla\bullet\left(\alpha_\text{eff}\nabla e\right)-\nabla\bullet(p^*\mathbf{u}).
(14)

The term \left(\alpha_\text{eff}\nabla e\right) is a aproxximation for the heat flux -q^{eff} so we get after sorting and using (10) (c.f. [7]):

\frac{\partial}{\partial t}\left(\varrho e_0\right)+\nabla\bullet\left(\varrho\mathbf{u} e_0+p\mathbf{u}+q^{eff}+\frac{2}{3}\varrho k\mathbf{u}\right)=0
(15)

and with \mathbf{u}=\mathbf{I}\bullet \mathbf{u}

\frac{\partial}{\partial t}\left(\varrho e_0\right)+\nabla\bullet\left(\varrho\mathbf{u} e_0+p\mathbf{u}+q^{eff}+\frac{2}{3}\varrho\mathbf{I}k\bullet\mathbf{u}\right)=0.
(16)


So if I interpret the following variables as averaged and favre averaged in the equations (7), (12) and (15):

\mathbf{u} \rightarrow \widetilde{\mathbf{u}}

\varrho \rightarrow \overline{\varrho}

\mathbf{\tau^{eff}} \rightarrow \mathbf{\widetilde{\tau^{eff}}}

p \rightarrow \overline{p}

\mathbf{D} \rightarrow \widetilde{\mathbf{D}}

e_0 \rightarrow \widetilde{e_0}

q^{eff} \rightarrow \widetilde{q^{eff}}


I get out of (7):

\frac{\partial\overline{\varrho}}{\partial t}+\nabla\bullet(\overline{\varrho}\widetilde{\mathbf{u}}) = 0.

which is equivalent to (1).


I get out of (12):

\frac{\partial}{\partial t}(\overline{\varrho}\widetilde{\mathbf{u}})+\nabla\bullet(\overline{\varrho}\widetilde{\mathbf{u}}\otimes\widetilde{\mathbf{u}}+\overline{p}\mathbf{I}+\frac{2}{3}\overline{\varrho} k\mathbf{I}-\mathbf{\widetilde{\tau^{eff}}})=0

and if I use \mathbf{\tau^{tot}}=\mathbf{\widetilde{\tau^{eff}}}-\frac{2}{3}\overline{\varrho}\mathbf{I}k I get:

\frac{\partial}{\partial t}(\overline{\varrho}\widetilde{\mathbf{u}})+\nabla\bullet(\overline{\varrho}\widetilde{\mathbf{u}}\otimes\widetilde{\mathbf{u}}+\overline{p}\mathbf{I}-\mathbf{\widetilde{\tau^{tot}}})=0

which is equivalent to (2).

I get out of (15):

\frac{\partial}{\partial t}\left(\overline{\varrho} \widetilde{e_0}\right)+\nabla\bullet\left(\overline{\varrho}\widetilde{e_0}\widetilde{\mathbf{u}} +\widetilde{\mathbf{u}}\overline{p}+\widetilde{q^{eff}}+\frac{2}{3}\overline{\varrho} k\widetilde{\mathbf{u}}\right)=0.

This equation is not equivalent. The variable \widetilde{q^{eff}} is the approximation for \widetilde{\mathbf{q}^{tot}}. This leads me to the assumption, that the term \widetilde{\mathbf{\tau^{lam}}} is totaly neglected and also the part 2\mu_t\widetilde{\mathbf{D}}-\frac{2}{3}\mu_t(\nabla\bullet\widetilde{\mathbf{u}})\mathbf{I} of the tensor \widetilde{\mathbf{\tau^{turb}}} is neglected in the total energy equation.



Now this is my actual interpretation of the Code and the equations which were aviable to me until now. Does someone agree or disagree with the above text?


The questions which are arising in me are:

If the pressure, which is modeled in OpenFoam, is divided up like shown in equation (10), is this neglected in the Calculation of the equation of state?

In which cases would the neglection of the terms in the total energy equation lead to a discrepancy of the solution which is not negligible?


[1] Hirsch, Charles (2007): Fundamentals of computational fluid dynamics. 2. ed. Amsterdam: Elsevier/Butterworth-Heinemann (Numerical computation of internal and external flows, / Charles Hirsch ; Vol. 1).

[2] Tobias Holzmann (2016): Mathematics, Numerics, Derivations and OpenFOAMŪ: Holzmann CFD.

[3] https://www.cfd-online.com/Wiki/Favr...okes_equations

[4] Wilcox, David C. (1993): Turbulence modeling for CFD. La Caņada, Calif.: DCW Industries Inc.

[5] https://turbmodels.larc.nasa.gov/implementrans.html

[6] https://www.openfoam.com/documentati...lence-ras.html

[7] https://cfd.direct/openfoam/energy-equation/
fedvasu, Taataa, randolph and 3 others like this.
patrex is offline   Reply With Quote

Old   December 4, 2017, 12:44
Default
  #3
Senior Member
 
Lucky Tran
Join Date: Apr 2011
Location: Orlando, FL USA
Posts: 3,830
Rep Power: 46
LuckyTran has a spectacular aura aboutLuckyTran has a spectacular aura about
Quote:
Originally Posted by patrex View Post
If the pressure, which is modeled in OpenFoam, is divided up like shown in equation (10), is this neglected in the Calculation of the equation of state?
Yes there is the problem with pressure in which OF solves for p* and not p, which can be inferred from this line where you see only grad(p).
Code:
solve(UEqn == -fvc::grad(p));
So this is the only p that OF keeps track of, which it goes and uses for its thermo. It's also the p that gets written to file. So the pressure that you get out of OF is not the true mechanical/thermodynamic pressure but a modified pressure. In theory, one could subtract the contribution from the tke before calculating the thermo, but this is not done in general. I can't think of a solver which actually does this off the top-of-my-head. One would have to check the source code to make sure.

However, the difference goes like the I^2 so it is small except for very high turbulence intensity scenarios. Note the origin comes from the Boussinesq hypothesis. It's not a compressibility issue. You find the same modified pressure p* for both RANS and FANS. Star-CCM and Fluent do the same thing. There is also a mention of it in Pope's book. It seems the general approach is to assume that the k contribution to pressure is negligible.

This thread is useful, on calculating divdevreff but also contains some comments in the discussion on the modified pressure.
Calculating divDevReff

Quote:
Originally Posted by patrex View Post
In which cases would the neglection of the terms in the total energy equation lead to a discrepancy of the solution which is not negligible?
It affects the pressure and them subsequent thermo quantities derived from pressure.
First you need super high values of k. That should be apparent because k modifies the pressure. For thermo, you need an EOS which is pressure dependent. The other thing you need to check is properties. Very often, one assumes ideal gas where the specific heat and internal energy are only temperature dependent. So the modified pressure does not affect energy equation beyond pressure. So very often this modified pressure does not break anything. But you can imagine that you have a EOS that is highly dependent on pressure where this becomes a problem; and shouldn't forget to mention shockwaves.

Quote:
Originally Posted by patrex View Post
This leads me to the assumption, that the term \widetilde{\mathbf{\tau^{lam}}} is totaly neglected and also the part 2\mu_t\widetilde{\mathbf{D}}-\frac{2}{3}\mu_t(\nabla\bullet\widetilde{\mathbf{u}})\mathbf{I} of the tensor \widetilde{\mathbf{\tau^{turb}}} is neglected in the total energy equation.
Yes, the viscous terms are neglected in the energy equation in many solvers.

Btw to answer the original question:
Quote:
Originally Posted by NablaDyn View Post
Dear FOAMers,

I'm a bit confused. I cannot determine exactly whether OpenFOAM uses Reynolds (RANS) or Favre averaging (FANS) when performing compressible 'RAS' simulations. The OpenFOAM terminology clearly implies it's RANS, so the turbulent fluctuations are not density-normalised. This assumption proves reasonable for many compressible flows (and thus would fit my needs) but not unconditionally for supersonic flows or such including combustion etc. Instead, rationality suggests OpenFOAM solves FANS since (to my knowlegde) its 'RAS' models are applicable to incompressible as well as compressible cases with no need for any user interaction. Or does OpenFOAM choose the best suited approach by itself based on solver conditioning (fvSolution) or else?
The FANS equations should look identical to the RANS equations w/o density fluctuations. The other interpretation is that if you use the RANS equations, the equation is also valid for compressible flows as long as you remember that you should replace the Reynolds-average variables with their Favre-counterparts. This ability to completely interchange variables and have the same physics is precisely why Favre-averaging is useful. This is not 100% true because actually you do not always get the same equations w/ RANS as FANS. There are some assumptions you have to make to force to be able to recover RANS from FANS.

From a solver perspective, there is a fixed set of equations that gets solved (and they are the same equation). Hence the solver, (OF or any other solver) does not know whether these variables are the Reynolds ones or Favre ones because it's the same equation. Only the interpretation of the variables needs to be different. Remember in FANS you do not solve for the true velocity but a biased velocity, etc.
NablaDyn, cryabroad and patrex like this.
LuckyTran is offline   Reply With Quote

Old   January 13, 2018, 09:06
Default
  #4
Member
 
NablaDyn's Avatar
 
Join Date: Oct 2015
Location: Germany
Posts: 84
Rep Power: 6
NablaDyn is on a distinguished road
Dear Patrick and Lucky,

thank you for your very detailed and helpful contributions to this thread. It really cleared things up.

Quote:
Originally Posted by LuckyTran View Post
From a solver perspective, there is a fixed set of equations that gets solved (and they are the same equation). Hence the solver, (OF or any other solver) does not know whether these variables are the Reynolds ones or Favre ones because it's the same equation. Only the interpretation of the variables needs to be different. Remember in FANS you do not solve for the true velocity but a biased velocity, etc.
Well, this now seems ridiculously obvious to me... However, I think I should have a closer look at Pope's work to get a better, intuitive feel regarding FANS.

Thanks again!

Martin
NablaDyn is offline   Reply With Quote

Old   June 12, 2018, 17:31
Default
  #5
Member
 
Join Date: Oct 2013
Posts: 88
Rep Power: 8
fedvasu is on a distinguished road
Quote:
Originally Posted by patrex View Post
Hello NablaDyn,

This is a question wich i ask myself a lot in the last time. So I am also really interested in the 'right' answer and i already think i did some progress in it.



Well I think the problem with the compressible Reynolds-averaged-Navier-Stokes-Equations is, if they are averaging is done with the time averaging:

\overline{\phi}(x)=\lim\limits_{\Delta t\rightarrow\infty}\frac{1}{\Delta t}\int_T^{T+\Delta t} \phi(x,t) \,dt

and the splitting into:

\phi=\overline{\phi}+\phi',

the resulting equations are very unhandy because of the products of fluctuations between density and othe variables like the velocity (c.f. [1], p.91). This is shown as example in [2]
for the continunity equation in index notation (c.f. [2], p.100):

\frac{\partial \overline{\varrho}}{\partial t} + \frac{\partial}{\partial x_i}(\overline{\varrho}\,\overline{u}_i+\overline{\varrho'u_i'})=0.

The equations for the momentum and e.g. the total energy will much more complex. Beacause of this a density-weighted averaged is introduced (c.f. e.g. [1], p.91):

\widetilde{\phi}=\frac{\overline{\varrho\phi}}{\overline{\varrho}}

and

\phi=\widetilde{\phi}+\phi''.


This leads to the following equations without sources and body forces in index notation (c.f. [3], Equations (12), (13) and (14)):

Continunity equation:

\frac{\partial \overline{\varrho}}{\partial t} +\frac{\partial}{\partial x_i}\left[ \overline{\varrho} \widetilde{u_i} \right] = 0

Momentum equation:

\frac{\partial}{\partial t}\left( \overline{\varrho} \widetilde{u_i} \right) +\frac{\partial}{\partial x_j}\left[\overline{\rho} \widetilde{u_i} \widetilde{u_j} + \overline{p} \delta_{ij} + \overline{\rho u''_i u''_j} - \overline{\tau_{ji}}\right]=0

Total energy:

\frac{\partial}{\partial t}\left( \overline{\varrho} \widetilde{e_0} \right) +\frac{\partial}{\partial x_j}\left[ \overline{\varrho} \widetilde{u_j} \widetilde{e_0} + \widetilde{u_j} \overline{p} + \overline{u''_j p} + \overline{\varrho u''_j e''_0} + \overline{q_j} - \overline{u_i \tau_{ij}}\right]= 0.

(You can also find these equations in [1], p.92 (without total energy equation) and in [4], p.176)

These equations are called Favre-averaged Navier-Stokes equations but are also termed as compressible Reynolds-averaged Navier-Stokes equations (c.f. [5]).

As shown in [3], after some assumptions the equations can be written as:

\frac{\partial \overline{\varrho}}{\partial t} +\frac{\partial}{\partial x_i}\left[ \overline{\varrho} \widetilde{u_i} \right] = 0

\frac{\partial}{\partial t}\left( \overline{\varrho} \widetilde{u_i} \right) +\frac{\partial}{\partial x_j}\left[\overline{\varrho} \widetilde{u_j} \widetilde{u_i}+ \overline{p} \delta_{ij}- \widetilde{\tau_{ji}^{tot}}\right]= 0

\frac{\partial}{\partial t}\left( \overline{\varrho} \widetilde{e_0} \right) +\frac{\partial}{\partial x_j}\left[ \overline{\varrho} \widetilde{u_j} \widetilde{e_0} + \widetilde{u_j} \overline{p} + \widetilde{q_j^{tot}} - \widetilde{u_i} \widetilde{\tau_{ij}^{tot}}\right] = 0,

where

\widetilde{\tau_{ij}^{tot}} \equiv \widetilde{\tau_{ij}^{lam}} + \widetilde{\tau_{ij}^{turb}}

\widetilde{\tau_{ij}^{lam}} \equiv \widetilde{\tau_{ij}} = \mu\left( \frac{\partial \widetilde{u_i} }{\partial x_j} + \frac{\partial \widetilde{u_j} }{\partial x_i} - \frac{2}{3} \frac{\partial \widetilde{u_k} }{\partial x_k} \delta_{ij}\right)

\widetilde{\tau_{ij}^{turb}} \equiv- \overline{\rho u''_i u''_j} \approx\mu_t\left( \frac{\partial \widetilde{u_i} }{\partial x_j} + \frac{\partial \widetilde{u_j} }{\partial x_i} - \frac{2}{3} \frac{\partial \widetilde{u_k} }{\partial x_k} \delta_{ij}\right) -\frac{2}{3} \overline{\rho} k \delta_{ij}.

In tensor notation it should be:

\frac{\partial\overline{\varrho}}{\partial t}+\nabla\bullet(\overline{\varrho}\widetilde{\mathbf{u}}) = 0
(1)

\frac{\partial}{\partial t}(\overline{\varrho}\widetilde{\mathbf{u}})+\nabla\bullet\left(\overline{\varrho}\widetilde{\mathbf{u}}\otimes\widetilde{\mathbf{u}}+\overline{p}\mathbf{I}-\widetilde{\mathbf{\tau^{tot}}}\right) = 0
(2)

\frac{\partial}{\partial t} (\overline{\varrho}\widetilde{e_0})+\nabla\bullet\left(\overline{\varrho}\widetilde{e_0}\widetilde{\mathbf{u}}+\widetilde{\mathbf{u}}\overline{p}+\widetilde{\mathbf{q}^{tot}}-\widetilde{\mathbf{\tau^{tot}}}\bullet\widetilde{\mathbf{u}}\right) = 0
(3)

\widetilde{\mathbf{\tau^{tot}}}=\widetilde{\mathbf{\tau^{lam}}}+\widetilde{\mathbf{\tau^{turb}}}
(4)

\widetilde{\mathbf{\tau^{lam}}}=2\mu\widetilde{\mathbf{D}}-\frac{2}{3}\mu(\nabla\bullet\widetilde{\mathbf{u}})\mathbf{I}
(5)

\widetilde{\mathbf{\tau^{turb}}}=2\mu_t\widetilde{\mathbf{D}}-\frac{2}{3}\mu_t(\nabla\bullet\widetilde{\mathbf{u}})\mathbf{I}-\frac{2}{3}\overline{\varrho}\mathbf{I}k.
(6)


Im looking on sonicFoam of OpenFOAM 4.1, which is a solver forr trans-sonic/supersonic, turbulent flow and there is first the continunity equation solved in rhoEqn.H (code-snippet of the equation):
Code:
     fvScalarMatrix rhoEqn
     (
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
         fvOptions(rho)
     );
which i think could be translated into:

\frac{\partial\varrho}{\partial t}+\nabla\bullet(\varrho \mathbf{u}) = 0
(7)

after that UEqn.H is called (code-snippet of the equation):

Code:
 fvVectorMatrix UEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
   + turbulence->divDevRhoReff(U)
  ==
     fvOptions(rho, U)
 );
.
.
.
 if (pimple.momentumPredictor())
 {
     solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);
 }

which i think could be translated into:

\frac{\partial}{\partial t}(\varrho\mathbf{u})+\nabla\bullet(\varrho\mathbf{u}\otimes\mathbf{u})=-\nabla p^*+\nabla\bullet\mathbf{\tau^{eff}}=-\nabla\bullet(p^*\mathbf{I})+\nabla\bullet\mathbf{\tau^{eff}}
(8)

(the +\nabla\bullet\mathbf{\tau^{eff}} is due to the negativ implementation of the terms in the function divDevRhoReff())

where (see [2], cap. 10.2):

\mathbf{\tau^{eff}}=2\mu_\text{eff}\left[\frac{1}{2}\left\{(\nabla\otimes\mathbf{u})+(\nabla\otimes\mathbf{u})^T\right\}\right]-\frac{2}{3}\mu_\text{eff}(\nabla\bullet\mathbf{u})\mathbf{I}
(9)
=2\mu\mathbf{D}-\frac{2}{3}\mu(\nabla\bullet\mathbf{u})\mathbf{I} +2\mu_t\mathbf{D}-\frac{2}{3}\mu_t(\nabla\bullet\mathbf{u})\mathbf{I}

and (c.f. [6]):

p^*=p+\frac{2}{3}\varrho k
(10)

if we use (10) with (8):

\frac{\partial}{\partial t}(\varrho\mathbf{u})+\nabla\bullet(\varrho\mathbf{u}\otimes\mathbf{u})=-\nabla\bullet((p+\frac{2}{3}\varrho k)\mathbf{I})+\nabla\bullet\mathbf{\tau^{eff}}
(11)

and sort it:

\frac{\partial}{\partial t}(\varrho\mathbf{u})+\nabla\bullet(\varrho\mathbf{u}\otimes\mathbf{u}+p\mathbf{I}+\frac{2}{3}\varrho k\mathbf{I}-\mathbf{\tau^{eff}})=0
(12)


and after that EEqn.H is called (code-snippet of the equation):

Code:
    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, e) + fvm::div(phi, e)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + fvc::div(fvc::absolute(phi/fvc::interpolate(rho), U), p, "div(phiv,p)")
      - fvm::laplacian(turbulence->alphaEff(), e)
     ==
        fvOptions(rho, e)
    );
which i think could be translated into:

\frac{\partial}{\partial t}\left(\varrho e+\frac{1}{2}\varrho|\mathbf{u}|^2\right)=-\nabla\bullet\left[\varrho\mathbf{u}\left(e+\frac{1}{2}|\mathbf{u}|^2\right)\right]+\nabla\bullet\left(\alpha_\text{eff}\nabla e\right)-\nabla\bullet(p^*\mathbf{u})
(13)

which is in coincidence with [7] if e_0= e+\frac{1}{2}|\mathbf{u}|^2 is used:

\frac{\partial}{\partial t}\left(\varrho e_0\right)=-\nabla\bullet\left(\varrho\mathbf{u} e_0\right)+\nabla\bullet\left(\alpha_\text{eff}\nabla e\right)-\nabla\bullet(p^*\mathbf{u}).
(14)

The term \left(\alpha_\text{eff}\nabla e\right) is a aproxximation for the heat flux -q^{eff} so we get after sorting and using (10) (c.f. [7]):

\frac{\partial}{\partial t}\left(\varrho e_0\right)+\nabla\bullet\left(\varrho\mathbf{u} e_0+p\mathbf{u}+q^{eff}+\frac{2}{3}\varrho k\mathbf{u}\right)=0
(15)

and with \mathbf{u}=\mathbf{I}\bullet \mathbf{u}

\frac{\partial}{\partial t}\left(\varrho e_0\right)+\nabla\bullet\left(\varrho\mathbf{u} e_0+p\mathbf{u}+q^{eff}+\frac{2}{3}\varrho\mathbf{I}k\bullet\mathbf{u}\right)=0.
(16)


So if I interpret the following variables as averaged and favre averaged in the equations (7), (12) and (15):

\mathbf{u} \rightarrow \widetilde{\mathbf{u}}

\varrho \rightarrow \overline{\varrho}

\mathbf{\tau^{eff}} \rightarrow \mathbf{\widetilde{\tau^{eff}}}

p \rightarrow \overline{p}

\mathbf{D} \rightarrow \widetilde{\mathbf{D}}

e_0 \rightarrow \widetilde{e_0}

q^{eff} \rightarrow \widetilde{q^{eff}}


I get out of (7):

\frac{\partial\overline{\varrho}}{\partial t}+\nabla\bullet(\overline{\varrho}\widetilde{\mathbf{u}}) = 0.

which is equivalent to (1).


I get out of (12):

\frac{\partial}{\partial t}(\overline{\varrho}\widetilde{\mathbf{u}})+\nabla\bullet(\overline{\varrho}\widetilde{\mathbf{u}}\otimes\widetilde{\mathbf{u}}+\overline{p}\mathbf{I}+\frac{2}{3}\overline{\varrho} k\mathbf{I}-\mathbf{\widetilde{\tau^{eff}}})=0

and if I use \mathbf{\tau^{tot}}=\mathbf{\widetilde{\tau^{eff}}}-\frac{2}{3}\overline{\varrho}\mathbf{I}k I get:

\frac{\partial}{\partial t}(\overline{\varrho}\widetilde{\mathbf{u}})+\nabla\bullet(\overline{\varrho}\widetilde{\mathbf{u}}\otimes\widetilde{\mathbf{u}}+\overline{p}\mathbf{I}-\mathbf{\widetilde{\tau^{tot}}})=0

which is equivalent to (2).

I get out of (15):

\frac{\partial}{\partial t}\left(\overline{\varrho} \widetilde{e_0}\right)+\nabla\bullet\left(\overline{\varrho}\widetilde{e_0}\widetilde{\mathbf{u}} +\widetilde{\mathbf{u}}\overline{p}+\widetilde{q^{eff}}+\frac{2}{3}\overline{\varrho} k\widetilde{\mathbf{u}}\right)=0.

This equation is not equivalent. The variable \widetilde{q^{eff}} is the approximation for \widetilde{\mathbf{q}^{tot}}. This leads me to the assumption, that the term \widetilde{\mathbf{\tau^{lam}}} is totaly neglected and also the part 2\mu_t\widetilde{\mathbf{D}}-\frac{2}{3}\mu_t(\nabla\bullet\widetilde{\mathbf{u}})\mathbf{I} of the tensor \widetilde{\mathbf{\tau^{turb}}} is neglected in the total energy equation.



Now this is my actual interpretation of the Code and the equations which were aviable to me until now. Does someone agree or disagree with the above text?


The questions which are arising in me are:

If the pressure, which is modeled in OpenFoam, is divided up like shown in equation (10), is this neglected in the Calculation of the equation of state?

In which cases would the neglection of the terms in the total energy equation lead to a discrepancy of the solution which is not negligible?


[1] Hirsch, Charles (2007): Fundamentals of computational fluid dynamics. 2. ed. Amsterdam: Elsevier/Butterworth-Heinemann (Numerical computation of internal and external flows, / Charles Hirsch ; Vol. 1).

[2] Tobias Holzmann (2016): Mathematics, Numerics, Derivations and OpenFOAMŪ: Holzmann CFD.

[3] https://www.cfd-online.com/Wiki/Favr...okes_equations

[4] Wilcox, David C. (1993): Turbulence modeling for CFD. La Caņada, Calif.: DCW Industries Inc.

[5] https://turbmodels.larc.nasa.gov/implementrans.html

[6] https://www.openfoam.com/documentati...lence-ras.html

[7] https://cfd.direct/openfoam/energy-equation/
This is a clear and detailed explanation of how averaging is done in OpenFOAM and how to interpret the notation.


Thanks.

Last edited by fedvasu; June 12, 2018 at 23:38.
fedvasu is offline   Reply With Quote

Old   October 13, 2018, 13:25
Default
  #6
Senior Member
 
Join Date: Mar 2011
Posts: 158
Rep Power: 11
tH3f0rC3 is on a distinguished road
Thank you for the detailed and illustrative contribution. It really helped me to better understand how OpenFoam deals with Reynolds and Favre averaging.


But what is still unclear for me is the following:
How do I get to know which solver uses RANS or FANS equations?


Do all solvers using the compressible NS equations solve for the FANS (=compressible RANS) equations?


When I assume, as you already indicated, that the same equations are solved, where is the massweighted time averaging done?



I do not find a hint in the code, where the velocity and specific enthalpy (which are the variables which are averaged with the density) are treated separately.


What also confuses me a little bit: No word about Favre averaging is spent here:

https://www.openfoam.com/documentati...lence-ras.html
So how are the compressible NS equations time averaged when using RANS equations?



Best Regards
tH3f0rC3 is offline   Reply With Quote

Old   October 13, 2018, 13:34
Default
  #7
Senior Member
 
Lucky Tran
Join Date: Apr 2011
Location: Orlando, FL USA
Posts: 3,830
Rep Power: 46
LuckyTran has a spectacular aura aboutLuckyTran has a spectacular aura about
Quote:
Originally Posted by tH3f0rC3 View Post
How do I get to know which solver uses RANS or FANS equations?

You never know. The solver solves an equation like Ax=b and x can be anything. You have to infer based on the models whether x is a RANS of FANS or neither.

In the case of RANS and FANS, the equations are identical (the A and b are the same). So you can say all solvers always solve the FANS. No one ever talks about it, because the equations are the same. Once you are given a solver you don't get to decide what x is. The solver solves for x.
LuckyTran is offline   Reply With Quote

Old   October 13, 2018, 14:18
Default
  #8
Senior Member
 
Join Date: Mar 2011
Posts: 158
Rep Power: 11
tH3f0rC3 is on a distinguished road
OK, yes, that's right.

If the equations are the same, x is the time averaged velocity when using the RANS equations and x is the massweighted time averaged velocity, for example.

Thus, somewhere in the code there must be hidden the corresponding averaging equation, which provides the x to the solver.
And I can't find such a hint in the code. Do you get my point?

Last edited by tH3f0rC3; October 13, 2018 at 14:20. Reason: Misspelling
tH3f0rC3 is offline   Reply With Quote

Old   October 14, 2018, 14:52
Default
  #9
New Member
 
Patrick H.
Join Date: Aug 2017
Location: Germany
Posts: 3
Rep Power: 4
patrex is on a distinguished road
It is the following way:

If you had all values continuous in time and spatial direction for a given problem. You could calculate the how ever shaped averaged value.

If there is a flow, for which the simplifications, which are made for the incompressible Reynolds-averaging, are valid, you could now calculate the averaged values from the continuous values. The averaged values which a RANS-Solver would calculate directly.

The same applies for a flow, for which the simplifications, made for the FANS-equation, are valid. As shown in different Sources and my previous post the Equations for incompressible RANS and the compressible FANS are under somehow shaped simplifications the same. So now it depends on your calculated flow how you have to interpret the of your solver calculated values. So if you only look at the equations (FANS and RANS) there could be two different flows, one compressible and one incompressible, which are different shaped in their actual values, but with proper averaging (FANS- and RANS-averaging consecutively) are the same in their averaged values.

Means: OpenFOAM does no averaging. The averaged values are calculated directly with the solver.
patrex is offline   Reply With Quote

Old   October 15, 2018, 07:01
Default
  #10
Senior Member
 
Join Date: Mar 2011
Posts: 158
Rep Power: 11
tH3f0rC3 is on a distinguished road
Thanks patrex for your explanation.
This point is also clear to me.


I will try it this way to explain my problem:

Let's assume we have a flow, for which the simplifications, which are made for the incompressible Reynolds-averaging, are valid.


Which solver do I use in contrast to a case, where we assume a flow, for which the simplifications, made for the FANS-equation, are valid.


After having solved the compressible Navier Stokes equations using for example the solver buoyantSimpleFoam I get values for u, p, T.


After solving the incompressible Nasvier-Stokes equations using for example the solver buoyantBoussineswSimpleFoam I also get values for u, p, T.


As far as I understand, the Reynolds-averaging is used to treat turbulence for incompressible problems and Favre-averaging to treat turbulence for compressible problems (to secure the validity of the continuity equation for the averaged flow).




Comming back to my first question:
-> How do I get to know which solver uses RANS or FANS equations?
tH3f0rC3 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
OpenFOAM 5.0 Released CFDFoundation OpenFOAM Announcements from OpenFOAM Foundation 11 June 5, 2018 23:48
Reynolds Stresses Tensor from OpenFoam Irstea OpenFOAM Post-Processing 0 August 29, 2017 09:04
Map of the OpenFOAM Forum - Understanding where to post your questions! wyldckat OpenFOAM 9 March 30, 2017 05:19
OpenFOAM v3.0+ ?? SBusch OpenFOAM 22 December 26, 2016 14:24
URGENT - Simulation setup using Reynolds Tensor Transport (R) - OpenFOAM beluiz93 OpenFOAM 1 December 2, 2016 06:51


All times are GMT -4. The time now is 01:09.