%                      Code_Saturne version 1.3
%                      ------------------------
%
%     This file is part of the Code_Saturne Kernel, element of the
%     Code_Saturne CFD tool.
% 
%     Copyright (C) 1998-2007 EDF S.A., France
%
%     contact: saturne-support@edf.fr
% 
%     The Code_Saturne Kernel is free software; you can redistribute it
%     and/or modify it under the terms of the GNU General Public License
%     as published by the Free Software Foundation; either version 2 of
%     the License, or (at your option) any later version.
% 
%     The Code_Saturne Kernel is distributed in the hope that it will be
%     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
%     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with the Code_Saturne Kernel; if not, write to the
%     Free Software Foundation, Inc.,
%     51 Franklin St, Fifth Floor,
%     Boston, MA  02110-1301  USA
%
%-----------------------------------------------------------------------
%
\programme{turrij}

\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Fonction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Le but de ce sous-programme est de r\'esoudre le syst\`eme des \'equations des
tensions de Reynolds et de la dissipation $\varepsilon$ de mani\`ere totalement d\'ecoupl\'ee dans le cadre de l'utilisation du mod\`ele $R_{ij}-\varepsilon$  LRR\footnote{la description du SSG est prévue pour une version ultérieure de la documentation} (option $\var{ITURB}=30$ dans \fort{usini1}).\\
Le tenseur sym\'etrique des tensions de Reynolds est not\'e $\tens{R}$. Les composantes de ce tenseur repr\'esentent le moment d'ordre deux de la vitesse : $R_{ij} = \overline{u_iu_j}$.

Pour chaque composante $R_{ij}$, on r\'esout :

\begin{equation}
\begin{array}{ll}
\displaystyle
\rho\frac{\partial R_{ij}}{\partial t} + 
\dive(\rho \vect{u}\,R_{ij} - \mu\,\grad{R_{ij}}) = &
\mathcal{P}_{ij} + \mathcal{G}_{ij}+\Phi_{ij} + \it{d}_{ij} - \varepsilon_{ij} +
R_{ij}\,\dive{(\rho \vect{u})} \\
& \displaystyle + \Gamma(R^{\,in}_{ij}-R_{ij}) + \alpha_{R_{ij}} R_{ij} + \beta_{R_{ij}}
\end{array}
\end{equation}

$\tens{\mathcal{P}}$ est le tenseur de production par cisaillement moyen : 

\begin{equation}
\displaystyle \mathcal{P}_{ij} = \displaystyle -\rho \left[ R_{ik} \frac{\partial u_j}{\partial x_k} + R_{jk} \frac{\partial u_i}{\partial x_k} \right]
\end{equation}


$\tens{\mathcal{G}}$ est le tenseur de production par gravit\'e :

\begin{equation}
\displaystyle
\mathcal{G}_{ij}= \left[ G_{ij} - C_3 (G_{ij}-\frac{1}{3} \delta_{ij} G_{kk}) \right]
\end{equation}

avec

\begin{equation}
\left\{
\begin{array} {c}
\displaystyle G_{ij} = - \frac{3}{2} \frac{C_{\mu}}{\sigma_{t}} \frac{k}{\varepsilon} (r_i g_j + r_j g_i) \\
\displaystyle k = \frac{1}{2} R_{ll} \\
\displaystyle r_i = R_{ik} \frac{\partial \rho}{\partial x_k} 
\end{array}\right.
\end{equation}

Dans ce qui pr\'ec\`ede, $k$ repr\'esente l'\'energie turbulente\footnote{Les
sommations se font sur l'indice $l$ et on applique plus
g\'en\'eralement la convention de sommation d'Einstein.}, $g_i$ la composante de
la gravit\'e dans la direction $i$, $\sigma_{t}$ le nombre de Prandlt turbulent  et $C_{\mu}$, $C_3$ des constantes d\'efinies dans Tab.~\ref{Base_Turrij_table_Cstes}.


$\tens{\Phi}$ est le terme de corr\'elations pression-d\'eformation. Il est mod\'elis\'e avec le terme de dissipation $\tens{\varepsilon}$ de la mani\`ere suivante :

\begin{equation}
\displaystyle
\Phi_{ij} - [\varepsilon_{ij}- \frac{2}{3} \rho \ \delta_{ij} \varepsilon] = \phi_{ij,1} + \phi_{ij,2} + \phi_{ij,w}
\end{equation}

Il en r\'esulte :

\begin{equation}
\displaystyle
\Phi_{ij} - \varepsilon_{ij} = \phi_{ij,1} + \phi_{ij,2} + \phi_{ij,w}  -\frac{2}{3} \rho \ \delta_{ij} \varepsilon
\end{equation}

Le terme $\phi_{ij,1}$ est un terme "lent" de retour \`a l'isotropie. Il est donn\'e par : 

\begin{equation}
\displaystyle
\phi_{ij,1} = -\rho\,C_1 \frac{\varepsilon}{k} (R_{ij} - \frac{2}{3} k \delta_{ij})
\end{equation}

Le terme $\phi_{ij,2}$ est un terme "rapide" d'isotropisation de la production. Il est donn\'e par : 
\begin{equation}
\displaystyle
\phi_{ij,2} = -\rho\,C_2 (\mathcal{P}_{ij} - \frac{2}{3} \mathcal{P} \delta_{ij}) 
\end{equation}

avec,

$$\displaystyle \mathcal{P} = \frac{1}{2} \mathcal{P}_{kk}$$

Le terme $\phi_{ij,w}$ est appel\'e "terme d'echo de paroi". Il n'est pas
utilis\'e par d\'efaut dans \CS, mais peut \^etre activ\'e avec $\var{IRIJEC} = 1$. Si $y$ repr\'esente la distance \`a la paroi :

\begin{equation}
\begin{array} {ll}
\displaystyle
\phi_{ij,w}  = &
\displaystyle \rho\,C'_1 \frac{k}{\varepsilon} \left[ R_{km} n_k n_m \delta_{ij} - 
\frac{3}{2} R_{ki} n_k n_j - 
\frac{3}{2} R_{kj} n_k n_i \right] f(\frac{l}{y})  \\
&
+\displaystyle \rho\,C'_2 \left[ \phi_{km,2} n_k n_m \delta_{ij} - 
\frac{3}{2} \phi_{ki,2} n_k n_j - 
\frac{3}{2} \phi_{kj,2} n_k n_i \right] f(\frac{l}{y})
\end{array}
\end{equation}

$f$ est une fonction d'amortissement construite pour valoir 1 en paroi et tendre
vers 0 en s'\'eloignant des parois.\\
La longueur $l$ repr\'esente 
$\displaystyle\frac{k^{\,\frac{3}{2}}}{\varepsilon}$, caract\'eristique de la turbulence. On prend :

\begin{equation}
f(\frac{l}{y}) = min(1, \ C^{\,0,75}_{\mu} \
\frac{k^{\,\frac{3}{2}}}{\varepsilon\ \kappa y}) 
\end{equation}


$\it{d}_{ij}$ est un terme de diffusion turbulente\footnote{Dans la litt\'erature, ce terme contient en g\'en\'eral la dissipation par viscosit\'e mol\'eculaire.} qui vaut :

\begin{equation}
\it{d}_{ij} = C_{S} \frac{\partial}{\partial x_k} (\rho \frac{k}{\varepsilon} R_{km} \frac{\partial R_{ij}}{\partial x_m})
\end{equation}

On notera par la suite $\displaystyle \tens{A} = C_S\,\rho\,\frac{k}{\varepsilon}\,\tens{R}$. Ainsi, $\displaystyle d_{ij} = \dive(\,\tens{A}\,\grad(R_{ij}))$ est une diffusion avec un coefficient tensoriel.

Le terme de dissipation turbulente $\tens{\varepsilon}$ est trait\'e dans ce qui pr\'ec\`ede avec le terme $\tens{\Phi}$.

$\Gamma$ est le terme source de masse\footnote{Dans ce cas, l'\'equation de continuit\'e s'\'ecrit : $\displaystyle \frac{\partial \rho}{\partial t} + \dive{(\rho \vect{u})} = \Gamma$.}, $R^{\,in}_{ij}$ est la valeur de $R_{ij}$ associ\'ee \`a la masse inject\'ee ou retir\'ee.

($\alpha_{R_{ij}}\,R_{ij} + \beta_{R_{ij}}$) repr\'esente le terme source
utilisateur (sous-programme \fort{ustsri}) \'eventuel avec une d\'ecomposition
permettant d'impliciter la partie $\alpha_{R_{ij}}\,R_{ij}$ si $\alpha_{R_{ij}} \geqslant 0$.

De m\^eme, on r\'esout une \'equation de convection/diffusion/termes sources pour la dissipation $\varepsilon$. Cette \'equation est tr\`es semblable \`a celle du mod\`ele $k-\varepsilon$ (voir \fort{turbke}), seuls les termes de viscosit\'e turbulente et de gravit\'e changent. On r\'esout :

\begin{equation}
\begin{array} {ll}
\displaystyle \rho\frac{\partial \varepsilon}{\partial t} +
\dive\left[\rho \vect{u}\,\varepsilon-
(\mu \grad{\varepsilon})\right] = & 
\displaystyle \it{d}_{\,\varepsilon}
+ C_{\varepsilon_1}\frac{\varepsilon}{k}\left[\mathcal{P}
+\mathcal{G}_{\varepsilon}\right]
-\rho C_{\varepsilon_2}\frac{\varepsilon^2}{k}
+\varepsilon\dive(\rho\vect{u})\\
&
\displaystyle
+\Gamma(\varepsilon^{\,in}-\varepsilon)
+\alpha_\varepsilon \varepsilon +\beta_\varepsilon
\end{array}
\end{equation}


$\it{d}_{\,\varepsilon}$ est le terme de diffusion turbulente : 
\begin{equation}
\displaystyle
\it{d}_{\,\varepsilon} = C_{\varepsilon} \displaystyle \frac{\partial}{\partial x_k} \left( \rho \frac{k}{\varepsilon} R_{km} \frac{\partial \varepsilon}{\partial x_m} \right)
\end{equation}
On notera par la suite $\tens{A'} = \displaystyle \rho \,C_{\varepsilon} \frac{k}{\varepsilon} \tens{R}$.
Le terme de diffusion turbulente est donc mod\'elis\'e par : $$\it{d}_{\,\varepsilon} =
\dive(\tens{A'}\,\grad(\varepsilon))$$
La viscosit\'e turbulente usuelle ($\nu_t$) en mod\`ele $k-\varepsilon$ est remplac\'ee par le tenseur visqueux~$\tens{A'}$.

$\mathcal{P}$ est le terme de production par cisaillement moyen :
$\mathcal{P} =\displaystyle \frac{1}{2} \mathcal{P}_{kk}$. Ce terme est
mod\'elis\'e avec la notion de viscosit\'e turbulente dans le cadre du mod\`ele
$k-\varepsilon$. Dans le cas pr\'esent, il est calcul\'e \`a l'aide des tensions
de Reynolds (\`a partir de $\mathcal{P}_{ij}$).

$\mathcal{G}_{\varepsilon}$ est le terme de production des effets de gravit\'e pour la variable $\varepsilon$.
\begin{equation}
\mathcal{G}_{\varepsilon} = max(0,\frac{1}{2}G_{kk})
\end{equation}
\begin{table}
{\scriptsize
\begin{center}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|}
\hline
$C_\mu$  & $C_{\varepsilon}$  & $C_{\varepsilon_1}$ &
$C_{\varepsilon_2}$  & $C_1$ & $C_2$ & $C_3$ & $C_S$
& $C'_1$ & $C'_2$ \\
\hline
$0,09$ & $ 0,18$ & $1,44$ & $1,92$ & $1,8$ & $0,6$ & $0,55$ & $0,22$ & $0,5$ & 
$0,3$ \\
\hline
\end{tabular}
\end{center}
}   
\caption{D\'efinition des constantes utilis\'ees.}\label{Base_Turrij_table_Cstes}
\end{table}

