CFD Online Discussion Forums - Favre or Reynolds Average in OpenFOAM?

Quote:

Originally Posted by patrex (Post 672782)

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.