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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
La num\'ero de la phase trait\'ee fait partie des arguments de \fort{turrij}. On
omettra volontairement de le pr\'eciser dans ce qui suit, on indiquera par $(\ )$ la
notion de tableau s'y rattachant.

\etape{Calcul des termes de production $\tens{\mathcal{P}}$}
\begin{itemize}
\item [$\star$] Initialisation \`a z\'ero du tableau \var{PRODUC} dimensionn\'e \`a $\var{NCEL}\times 6$.
\item [$\star$] On appelle trois fois \fort{grdcel} pour calculer les gradients des composantes de la vitesse $u$, $v$ et
$w$ prises au temps $n$. 

Au final, on a :\\
$\displaystyle
\begin{array} {ll}
\var{PRODUC(1,IEL)} = & \displaystyle - 2 \left[ R_{11}^{\,n} \frac{\partial u^{\,n}} {\partial x} +R_{12}^{\,n} \frac{\partial u^{\,n}} {\partial y}+R_{13}^{\,n} \frac{\partial u^{\,n}} {\partial z} \right] \text{	(production de $R_{11}^{\,n}$)}\\ 
\var{PRODUC(2,IEL)} = & \displaystyle - 2 \left[ R_{12}^{\,n} \frac{\partial v^{\,n}} {\partial x} +R_{22}^{\,n} \frac{\partial v^{\,n}} {\partial y}+R_{23}^{\,n} \frac{\partial v^{\,n}} {\partial z} \right] \text{	(production de $R_{22}^{\,n}$)}\\
\var{PRODUC(3,IEL)} = & \displaystyle - 2 \left[ R_{13}^{\,n} \frac{\partial w^{\,n}} {\partial x} +R_{23}^{\,n} \frac{\partial w^{\,n}} {\partial y}+R_{33}^{\,n} \frac{\partial w^{\,n}} {\partial z} \right] \text{	(production de $R_{33}^{\,n}$)}\\
\var{PRODUC(4,IEL)} = & \displaystyle - \left[ R_{12}^{\,n} \frac{\partial u^{\,n}} {\partial x} +R_{22}^{\,n} \frac{\partial u^{\,n}} {\partial y}+R_{23}^{\,n} \frac{\partial u^{\,n}} {\partial z} \right] \\
& \displaystyle - \left[ R_{11}^{\,n} \frac{\partial v^{\,n}} {\partial x} +R_{12}^{\,n} \frac{\partial v^{\,n}} {\partial y}+R_{13}^{\,n} \frac{\partial v^{\,n}} {\partial z} \right] \text{	(production de $R_{12}^{\,n}$)} \\
\var{PRODUC(5,IEL)} = & \displaystyle - \left[ R_{13}^{\,n} \frac{\partial u^{\,n}} {\partial x} +R_{23}^{\,n} \frac{\partial u^{\,n}} {\partial y}+R_{33}^{\,n} \frac{\partial u^{\,n}} {\partial z} \right] \\ 
& \displaystyle - \left[ R_{11}^{\,n} \frac{\partial w^{\,n}} {\partial x} +R_{12}^{\,n} \frac{\partial w^{\,n}} {\partial y}+R_{13}^{\,n} \frac{\partial w^{\,n}} {\partial z} \right] \text{	(production de $R_{13}^{\,n}$)} \\
\var{PRODUC(6,IEL)} = & \displaystyle - \left[ R_{13}^{\,n} \frac{\partial v^{\,n}} {\partial x} +R_{23}^{\,n} \frac{\partial v^{\,n}} {\partial y}+R_{33}^{\,n} \frac{\partial v^{\,n}} {\partial z} \right] \\
& \displaystyle - \left[ R_{12}^{\,n} \frac{\partial w^{\,n}} {\partial x} +R_{22}^{\,n} \frac{\partial w^{\,n}} {\partial y}+R_{23}^{\,n} \frac{\partial w^{\,n}} {\partial z} \right]  \text{	(production de $R_{23}^{\,n}$)}
\end{array}
$
\end{itemize}

\etape{Calcul du gradient de la masse volumique $\rho^n$ prise au d\'ebut du pas
de temps courant\footnote{{\it i.e.} calcul\'ee \`a partir des
variables du pas de temps pr\'ec\'edent $n$ si n\'ecessaire.} $(n+1)$}
Ce calcul n'a lieu que si les termes de gravit\'e doivent \^etre pris en compte
($\var{IGRARI()} =1$). 
\begin{itemize}
\item [$\star$] Appel de \fort{grdcel}  pour calculer le gradient de $\rho^n$
dans les trois directions de l'espace. Les conditions aux limites sur $\rho^n$
sont des conditions de Dirichlet puisque la valeur de $\rho^n$ aux faces de bord
$ik$ (variable \var{IFAC}) est connue et vaut $\rho_{\,b_{\,ik}}$. Pour \'ecrire les conditions aux limites
sous la forme habituelle, $$\rho_{\,b_{\,ik}} = \var{COEFA} + \var{COEFB}
\,\rho^n_{\,I'}$$ on pose alors $\var{COEFA}=
\var{PROPCE(IFAC,IPPROB(IROM(IPHAS)))}$ et $\var{COEFB} = \var{VISCB} = 0$.\\
$\var{PROPCE(1,IPPROB(IROM(IPHAS)))}$ (resp.$\var{VISCB}$) est utilis\'e en lieu
et place de l'habituel \var{COEFA} ($\var{COEFB}$), lors de l'appel \`a \fort{grdcel}.\\
On a donc :\\
$\displaystyle \var{GRAROX}= \frac{\partial \rho^n}{\partial x}\ $,$\displaystyle \ \var{GRAROY}= \frac{\partial
\rho^n}{\partial y}$ et $
\displaystyle \ \var{GRAROZ}= \frac{\partial \rho^n}{\partial z}\ $.

\end{itemize}

Le gradient de $\rho^n$ servira \`a calculer les termes de production par effets de gravit\'e si ces derniers sont pris en compte.

\etape{Boucle \var{ISOU} de $1$ \`a $6$ sur les tensions de Reynolds}
Pour $\var{ISOU} = 1,2,3,4,5,6$, on r\'esout respectivement et dans
l'ordre  les
\'equations de $R_{11}$, $R_{22}$, $R_{33}$, $R_{12}$, $R_{13}$ et $R_{23}$ par
l'appel au sous-programme \fort{resrij}.\\
La r\'esolution se fait par incr\'ement $\delta {R}_{ij}^{\,n+1,k+1}$ , en utilisant la m\^eme m\'ethode que
celle d\'ecrite dans le sous-programme \fort{codits}. On adopte ici les m\^emes notations.
\var{SMBR} est le second membre du syst\`eme \`a inverser, syst\`eme portant sur
les incr\'ements de la variable. \var{ROVSDT} repr\'esente la diagonale de la
matrice, hors convection/diffusion.\\
On va r\'esoudre l'\'equation (\ref{Base_Turrij_Eq_Temp_Rij}) sous forme incr\'ementale en
utilisant \fort{codits}, soit :
\begin{equation}\label{Base_Turrij_Eq_Temp_deltaRij}
\begin{array}{ll}
&\displaystyle \underbrace{\left(\frac {\rho^n_L}{\Delta t^n} 
+ \rho^n_L \,C_1\,\frac{\varepsilon^n_L}{k^n_L}(1-\frac{\delta_{ij}}{3})
 - m^n_{\,lm} + \Gamma_L\,+ max(-\alpha^n_{R_{ij}},0)\right)\,|\Omega_l|}
_{\text {$\var{ROVSDT}$ contribuant
\`a la diagonale de la matrice simplifi\'ee de \fort{matrix}}}\,(\delta{R}_{ij}^{\,n+1,p+1})_{\,L}\\\\
&  \underbrace{+\sum\limits_{m\in Vois(l)}\displaystyle \left[
 m^n_{\,lm} \delta R_{ij,\,f_{\,lm}}^{\,n+1,p+1}
- (\mu^n_{\,lm} + \gamma^n_{\,lm})\
\frac{({\delta R}_{ij}^{\,n+1,p+1})_{M}-({\delta R}_{ij}^{\,n+1,p+1})_{L})}{\overline{L'M'}}\,
S_{\,lm} \right]}_{\text { convection upwind pur et diffusion non reconstruite
relatives \`a la matrice simplifi\'ee de \fort{matrix}\footnotemark}} \\ 
% voir le texte de la footmark plus bas
&= - \displaystyle\frac {\rho^n_L}{\Delta t^n}\,\left(\,(R^{\,n+1,p}_{ij})_L - (R^{\,n}_{ij})_L\,\right)\\  
&-\,\underbrace{\displaystyle\int_{\Omega_l} \left(
\dive\,[\,(\rho\,\vect{u})^n\,R^{\,n+1,p}_{ij} - (\mu^n\,+ \gamma^n\,)\,
\grad{R^{\,n+1,p}_{ij}}\,]\right)\,d\Omega}_{\text {convection et diffusion
trait\'ees par \fort{bilsc2}}}\\
&+\displaystyle \int_{\Omega_l} \left[\,\mathcal{P}^{\,n+1,p}_{ij} + \mathcal{G}^{\,n+1,p}_{ij}
- \displaystyle\rho^n \,C_1\,\frac{\varepsilon^n}{k^n}\left[R^{\,n+1,p}_{ij}- 
\frac{2}{3}\,k^n\,\delta_{ij}\right] + \phi^{\,n+1,p}_{ij,2} +
\phi^{\,n+1,p}_{ij,w}\,\right]\, d\Omega \\
& + \displaystyle\int_{\Omega_l} \left[- \frac{2}{3} \rho^n \varepsilon^n \delta_{ij} 
 + \Gamma\,(\,R^{\,in}_{ij} - R^{\,n+1,p}_{ij}\,) +
\alpha^n_{R_{ij}}\,R^{\,n+1,p}_{ij}+ \beta^n_{R_{ij}}\right]\, d\Omega\\
&+ \sum\limits_{m\in
Vois(l)}\displaystyle \left[\ \tens{E}^n\,\grad{R}^{\,n+1,p}_{ij} \right]_{\,lm}\,.\,\vect{n}_{\,lm}S_{\,lm}\\
&+ \sum\limits_{m\in Vois(l)}\displaystyle \left[\
\tens{D}^n\,\grad{R}^{\,n+1,p}_{ij} \right]_{\,lm}\,.\,\vect{n}_{\,lm}S_{\,lm}\\
&- \sum\limits_{m\in Vois(l)} \gamma^n_{\,lm} \left( \grad{R}^{\,n+1,p}_{ij}\,.\,\vect{n}_{\,lm} \right)  S_{\,lm}\\
&+ \sum\limits_{m\in Vois(l)}  m^n_{\,lm}\,(R^{\,n+1,p}_{ij})_L\\
\end{array}
\end{equation}
% si on ne fait pas comme ca, il n'apparait pas
\footnotetext[\thefootnote]{Si $\var{IRIJNU} = 1$, on remplace  $\mu^n_{\,lm}$ par $(\mu +
\mu_t)^n_{\,lm}$ dans l'expression de la diffusion non reconstruite
associ\'ee \`a la matrice simplifi\'ee de \fort{matrix} ($\mu_t$ d\'esigne la
viscosit\'e turbulente calcul\'ee comme en $k-\varepsilon$).}

o\`u on rappelle :\\
pour $n$ donn\'e entier positif, on d\'efinit la suite 
 $({R}_{ij}^{\,n+1,p})_{p \in \grandN}$
 par :
\begin{equation}\notag
\left\{\begin{array}{l}
{R}_{ij}^{\,n+1,0} = {R}_{ij}^{\,n}\\
{R}_{ij}^{\,n+1,p+1} = {R}_{ij}^{\,n+1,p} + \delta{R}_{ij}^{\,n+1,p+1} \\
\end{array}\right.
\end{equation}
$(\delta{R}_{ij}^{\,n+1,p+1})_{\,L}$ d\'esigne la valeur sur l'\'el\'ement
$\Omega_l$ du $\text{$(\,p+1\,)$-i\`eme}$ incr\'ement de ${R}_{ij}^{\,n+1}$, 
$ m^n_{\,lm}$ le flux de masse \`a l'instant $n$ \`a travers la face $lm$,  
$\delta R_{ij,\,f_{\,lm}}^{\,n+1,p+1}$ vaut $({\delta
R}_{ij}^{\,n+1,p+1})_{L}$  si $ m^n_{\,lm} \geqslant 0$, $({\delta
R}_{ij}^{\,n+1,p+1})_{M}$ sinon,  
$\mathcal{P}^{\,n+1,p}_{ij}$, $\phi^{\,n+1,p}_{ij,2}$, $\phi^{\,n+1,p}_{ij,w}$ les valeurs
des quantit\'es associ\'ees correspondant \`a l'incr\'ement
$(\delta{R}_{ij}^{\,n+1,p})$.\\



Tous ces termes sont calcul\'es comme suit :
\begin{itemize}
\item Terme de gauche de l'\'equation (\ref{Base_Turrij_Eq_Temp_deltaRij})\\
Dans \fort{resrij} est calcul\'ee la variable \var{ROVSDT}. Les autres
termes sont compl\'et\'es par \fort{codits}, lors de la construction de la matrice simplifi\'ee , {\it via} un
appel au sous-programme \fort{matrix}. La quantit\'e
 $(\mu^n_{\,lm} + \gamma^n_{\,lm})$ \`a la face $lm$ est calcul\'ee lors de l'appel \`a
\fort{visort}.\\ 
\item Second membre de l'\'equation (\ref{Base_Turrij_Eq_Temp_deltaRij})\\
Le premier terme non d\'etaill\'e est calcul\'e par le sous-programme
\fort{bilsc2}, qui applique le sch\'ema convectif choisi par l'utilisateur, qui
reconstruit ou non selon le souhait de l'utilisateur les gradients intervenants
dans la convection-diffusion.\\ 
Les termes sans accolade sont, eux, compl\`etement explicites et ajout\'es au fur et
\`a mesure dans \var{SMBR} pour former 
l'expression $f^{\,exp}_s$ de \fort{codits}.
\end{itemize}
On d\'ecrit ci-dessous les \'etapes de \fort{resrij} :
\begin{itemize}

\item DELTIJ = 1, pour $\var{ISOU} \leqslant 3$ et DELTIJ = 0  Si $\var{ISOU} >
3$. Cette valeur repr\'esente le symbole de Kroeneker $\delta_{ij}$.

\item Initialisation \`a z\'ero de \var{SMBR} (tableau contenant le second
membre) et \var{ROVSDT} (tableau contenant la diagonale de la matrice sauf celle
relative \`a la contribution de la
diagonale des op\'erateurs de convection et de diffusion lin\'earis\'es
\footnote{qui correspondent aux sch\'emas convectif upwind pur et diffusif sans
reconstruction.}), tous deux de dimension $\var{NCEL}$.

\item Lecture et prise en compte des termes sources utilisateur pour la variable $R_{ij}$ 

Appel \`a \fort{ustsri} pour charger les termes sources utilisateurs. Ils sont
stock\'es comme suit. Pour la cellule $\Omega_l$ de centre $L$, repr\'esent\'ee par $\var{IEL}$, on a :\\
\begin{equation}\notag
\left\{\begin{array}{lll}
&\var{ROVSDT(IEL)}&= |\Omega_l| \ \alpha_{R_{ij}}\\
&\var{SMBR(IEL)}&=|\Omega_l| \ \beta_{R_{ij}}\\
\end{array}\right.
\end{equation}
On affecte alors les valeurs ad\'equates au second membre \var{SMBR} et \`a la
diagonale \var{ROVSDT} comme suit :
\begin{equation}\notag
\left\{\begin{array}{lll}
&\var{SMBR(IEL)} &= \var{SMBR(IEL)} +\ |\Omega_l| \ \alpha_{R_{ij}} \ (R^n_{ij})_L \\
&\var{ROVSDT(IEL)}&= \text{max }(-\ |\Omega_l| \ \alpha_{R_{ij}},0)\\
\end{array}\right.
\end{equation}
La valeur de $ \var{ROVSDT}$ est ainsi calcul\'ee pour des raisons de stabilit\'e
num\'erique. En effet, on ne rajoute sur la diagonale que les valeurs positives,
ce qui correspond physiquement \`a impliciter les termes de rappel uniquement. 
\item{Calcul du terme source de masse  si $\Gamma_L > 0$}

Appel de \fort{catsma} et incr\'ementation si n\'ecessaire de \var{SMBR} et
\var{ROVSDT} {\it via} :\\
\begin{equation}\notag
\left\{\begin{array}{lll}
\displaystyle \var{SMBR(IEL)} = \var{SMBR(IEL)} + |\Omega_l| \ \Gamma_L \
\left[(R^{\,in}_{ij})_L - (R^{\,n}_{ij})_L \right] \\
\displaystyle \var{ROVSDT(IEL)}=\var{ROVSDT(IEL)} + |\Omega_l| \ \Gamma_L
\end{array}\right.
\end{equation}
\item Calcul du terme d'accumulation de masse et du terme instationnaire  

On stocke $\displaystyle \var{W1}= \int_{\Omega_l}\dive\,(\rho\,\vect{u})\,d\Omega$
calcul\'e par \fort{divmas} \`a l'aide des flux de masse aux faces internes 
$ m^n_{\,lm}=\sum\limits_{m\in Vois(l)}{(\rho \vect{u})_{\,lm}^n} \text{.}\,
\vect{S}_{\,lm} $ (tableau \var{FLUMAS}) et des flux de masse aux bords  $ m^n_{\,b_{lk}} = \sum\limits_{k\in{\gamma_b(l)}}{(\rho \vect{u})_{\,{b}_{lk}}^n} \text{.}\,
\vect{S}_{\,{b}_{lk}} $ (tableau \var{FLUMAB}).
On incr\'emente ensuite \var{SMBR} et \var{ROVSDT}.
\begin{equation}\notag
\left\{\begin{array}{lll}
&\var{SMBR(IEL)} &= \var{SMBR(IEL)} + \var{ICONV}\  (R^n_{ij})_L\,(\displaystyle
\int_{\Omega_l}\dive\,(\rho\,\vect{u})\ d\Omega) \\
&\var{ROVSDT(IEL)}& = \var{ROVSDT(IEL)} +  \var{ISTAT}\,\displaystyle
\frac{\rho^n_L \ |\Omega_l|}{\Delta t^n} -  \var{ICONV}\ (\displaystyle
\int_{\Omega_l}\dive\,(\rho\,\vect{u})\ d\Omega) \\
\end{array}\right.
\end{equation}
\item Calcul des termes sources de production, des termes $\displaystyle
\phi_{\,ij,1}+\phi_{\,ij,2}$ et de la dissipation~$\displaystyle-\frac{2}{3} \varepsilon\,\delta_{\,ij}$ :

On effectue une boucle d'indice \var{IEL} sur les cellules $\Omega_l$ de centre $L$ :
\begin{itemize}
\item [$\Rightarrow$] $\displaystyle \var{TRPROD}= \frac{1}{2} (\mathcal{P}^n_{ii})_L = \frac{1}{2} \left[ \var{PRODUC(1,IEL)} +  \var{PRODUC(2,IEL)} +  \var{PRODUC(3,IEL)} \right] $
\item [$\Rightarrow$] $\displaystyle \var{TRRIJ }= \frac{1}{2} (R^n_{ii})_L $ 
\item [$\Rightarrow$] $\displaystyle \var{SMBR(IEL)} =\ \var{SMBR(IEL)}\ +$\\
$\ \displaystyle\rho^n_L |\Omega_l| \left[ \displaystyle
\frac{2}{3}\,\delta_{\,ij} \left( \ \displaystyle \frac{ C_2}{2}\,(\mathcal{P}^n_{ii})_L\ +
(C_1-1)\ \varepsilon^n_L\, \right)\right.$\\
$ + \left.\ (1-C_2) \ \var{PRODUC(ISOU,IEL)} -
\displaystyle C_1\ \frac{2\,\varepsilon^n_L}{(R^n_{ii})_L}\ (R^n_{ij})_L \right]$
\item [$\Rightarrow$] $\displaystyle \var{ROVSDT(IEL)} = \var{ROVSDT(IEL)} +
\rho^n_L \ |\Omega_l| \ (- \displaystyle \frac{1}{3} \ \,\delta_{\,ij} + 1) \ C_1
\ \frac{2\ \varepsilon^n_L}{(R^n_{ii})_L}$
\end{itemize}
\item Appel de \fort{rijech} pour le calcul des termes d'\'echo de paroi
 $\phi^n_{ij,w}$ si $\var{IRIJEC()}=1$ et ajout dans \var{SMBR}.\\ 
$\var{SMBR} = \var{SMBR} + \phi^n_{ij,w}$\\
Suivant son mode de calcul (\var{ICDPAR}), la distance à la paroi est directement accessible 
par \var{RA(IDIPAR+IEL-1)} (\var{|ICDPAR|} = 1) ou bien 
est calcul\'ee \`a partir de $\var{IA(IIFAPA(IPHAS)+IEL - 1)}$, 
qui donne pour l'\'el\'ement $\var{IEL}$ le num\'ero de la face de bord
paroi la plus  proche (\var{|ICDPAR|} = 2). Ces tableaux ont \'et\'e renseign\'e une fois pour toutes au
d\'ebut de calcul.

\item  Appel de \fort{rijthe} pour le calcul des termes de gravit\'e $\mathcal{G}^n_{ij}$ et ajout dans \var{SMBR}. 

Ce calcul n'a lieu que si $\var{IGRARI()} = 1$.
$ \var{SMBR} = \var{SMBR} + \mathcal{G}^n_{ij}$ 
\item Calcul de la partie explicite du terme de diffusion
 $\dive{\,\left[\tens{A}\,\grad{R}^{\,n}_{ij}\right]}$, plus pr\'ecis\'ement
des contributions du terme extradiagonal pris aux faces purement internes
(remplissage du tableau \var{VISCF}), puis aux faces de bord (remplissage du
tableau \var{VISCB}).
\begin{itemize}
\item [$\star$] Appel de \fort{grdcel} pour le calcul du gradient de
$R^{\,n}_{ij}$ dans chaque direction. Ces gradients sont respectivement
stock\'es dans les tableaux de travail \var{W1}, \var{W2} et \var{W3}.

\item [$\star$] boucle d'indice \var{IEL} sur les cellules $\Omega_l$ de centre
$L$ pour le
calcul de $\tens{E}^n\,\grad{R}^{\,n}_{ij}$ aux cellules dans un premier temps :\\
\begin{itemize}
\item [$\Rightarrow$] $\displaystyle \var{TRRIJ}= \frac{1}{2} (R^{\,n}_{ii})_L $ 
\item [$\Rightarrow$] $\displaystyle \var{CSTRIJ} = \rho^n_L\ C_S \ \displaystyle\frac{(R^n_{ii})_L}{2\,\varepsilon^n_L}$
\item [$\Rightarrow$] $\displaystyle \var{W4(IEL)} = \rho^n_L\ C_S\
\displaystyle\frac{(R^n_{ii})_L}{2\,\varepsilon^n_L} \left[\,(R^{\,n}_{12})_L \ \var{W2(IEL)} +
(R^{\,n}_{13})_L \ \var{W3(IEL)}\,\right]$
\item [$\Rightarrow$] $\displaystyle \var{W5(IEL)} = \rho^n_L\ C_S\
\displaystyle\frac{(R^n_{ii})_L}{2\,\varepsilon^n_L} \left[\,(R^{\,n}_{12})_L \ \var{W1(IEL)} +
(R^{\,n}_{23})_L \ \var{W3(IEL)}\,\right]$
\item [$\Rightarrow$] $\displaystyle \var{W6(IEL)} = \rho^n_L\ C_S\
\displaystyle\frac{(R^n_{ii})_L}{2\,\varepsilon^n_L} \left[\,(R^{\,n}_{13})_L \ \var{W1(IEL)} + (R^{\,n}_{23})_L \ \var{W2(IEL)}\,\right]$
\end{itemize}



\item [$\star$] Appel de \fort{vectds}\footnote{Le r\'esultat est stock\'e dans
\var{VISCF} et \var{VISCB}. Dans \fort{vectds}, les valeurs aux faces internes
sont interpol\'ees lin\'eairement sans reconstruction et \var{VISCB} est mis \`a
z\'ero.} pour assembler $\displaystyle\left[ \tens{E}^n\,\grad{R}^{\,n}_{ij}
\right]\,.\,\vect{n}_{\,lm}S_{\,lm}$ aux faces $lm$.
\item [$\star$] Appel de \fort{divmas} pour calculer la divergence du flux d\'efini par \var{VISCF} et \var{VISCB}.
Le r\'esultat est stock\'e dans \var{W4}.\\
Ajout au second membre \var{SMBR}.\\
\var{SMBR} = \var{SMBR} + \var{W4}
\end{itemize}

A l'issue de cette \'etape, seule la partie extradiagonale de la diffusion prise
enti\`erement explicite~:
 $$\sum\limits_{m\in
Vois(l)}\left[\ \tens{E}^n\,\grad{R}^{\,n}_{ij} \right]_{\,lm}\,.\,\vect{n}_{\,lm}S_{\,lm}$$ a \'et\'e calcul\'ee.\\

\item Calcul de la partie diagonale du terme de diffusion\footnote{Seule la
composante normale  du  gradient de $R_{ij}$ aux faces sera implicite.} :\\
Comme on l'a d\'eja vu, une partie de cette contribution sera trait\'ee en
implicite (celle relative \`a la composante normale du gradient) et donc
ajout\'ee au second membre par \fort{bilsc2} ; l'autre 
partie sera explicite et prise en compte dans $f_s^{\,exp}$.
\begin{itemize}
\item [$\star$] On effectue une boucle d'indice \var{IEL} sur les cellules
$\Omega_l$ de centre $L$ :
\begin{itemize}
\item [$\Rightarrow$] $\displaystyle \var{TRRIJ }= \frac{1}{2} (R^{\,n}_{ii})_L $ 
\item [$\Rightarrow$] $\displaystyle \var{CSTRIJ} = \rho^n_L \ C_S \ \frac{(R^{\,n}_{ii})_L}{2\,\varepsilon^n_L}$
\item [$\Rightarrow$] $\displaystyle \var{W4(IEL)} = \rho^n_L \ C_S \
\frac{(R^{\,n}_{ii})_L}{2\,\varepsilon^n_L} \ (R^{\,n}_{11})_L$
\item [$\Rightarrow$] $\displaystyle \var{W5(IEL)} = \rho^n_L \ C_S \ \frac{(R^{\,n}_{ii})_L}{2\,\varepsilon^n_L}\ (R^n_{22})_L$
\item [$\Rightarrow$] $\displaystyle \var{W6(IEL)} = \rho^n_L \ C_S \ \frac{(R^{\,n}_{ii})_L}{2\,\varepsilon^n_L} \ (R^n_{33})_L$
\end{itemize}

%\item Traitement du parall\'elisme et de la p\'eriodicit\'e.

\item [$\star$] On effectue une boucle d'indice \var{IFAC} sur les faces
purement internes $lm$ pour remplir le tableau \var{VISCF} :
\begin{itemize}
\item [$\Rightarrow$] Identification des cellules $\Omega_l$ et $\Omega_m$ de
centre respectif $L$ (variable \var{II}) et $M$ (variable \var{JJ}), se trouvant de chaque c\^ot\'e de la face
$lm$\footnote{La normale \`a la face est orient\'ee de L vers M.}.
\item [$\Rightarrow$] Calcul du carr\'e de la surface de la face. La valeur est
stock\'ee dans le tableau \var{SURFN2}.
\item [$\Rightarrow$] Interpolation du gradient de $R^{\,n}_{ij}$ \`a la face
$lm$ (gradient facette $\left[\grad{R}^{\,n}_{ij}\right]_{\,lm}$) :
\begin{equation}\notag
\left\{\begin{array}{ll}
\var{GRDPX} &= \displaystyle \frac{1}{2} \left(\var{W1(II)} + \var{W1(JJ)}
\right) \\
&\\
\var{GRDPY} &= \displaystyle \frac{1}{2} \left(\var{W2(II)} + \var{W2(JJ)} \right) \\
&\\
\var{GRDPZ} &= \displaystyle \frac{1}{2} \left(\var{W3(II)} + \var{W3(JJ)} \right)
\end{array}\right.
\end{equation}
\item [$\Rightarrow$] Calcul du gradient de $R^{\,n}_{ij}$ normal \`a la face
$lm$, $\left[\grad{R}^{\,n}_{ij}\right]_{\,lm}.\vect{n}_{\,lm}\,S_{\,lm}$.\\

$\displaystyle \var{GRDSN} =  \var{GRDPX} \ \var{SURFAC(1,IFAC)} + \var{GRDPY} \ \var{SURFAC(2,IFAC)} +  \var{GRDPZ} \ \var{SURFAC(3,IFAC)}$
$\var{SURFAC}$ est le vecteur surface de la face \var{IFAC}.


\item [$\Rightarrow$] calcul de
 $\left[\grad{R^{\,n}_{ij}} - (\grad
R^{\,n}_{ij}\,.\,\vect{n}_{\,lm})\vect{n}_{\,lm}\right]$, les vecteurs \'etant
calcul\'es \`a la face $lm$ :
\begin{equation}\notag
\left\{\begin{array}{lll} 
&\displaystyle \var{GRDPX} &= \var{GRDPX} - \displaystyle\frac{\var{GRDSN}}{\var{SURFN2}} \ \var{SURFAC(1,IFAC)}\\
&&\\
&\displaystyle \var{GRDPY} &= \var{GRDPY} - \displaystyle\frac{\var{GRDSN}}{\var{SURFN2}} \ \var{SURFAC(2,IFAC)} \\
&&\\
&\displaystyle \var{GRDPZ} &= \var{GRDPZ} - \displaystyle\frac{\var{GRDSN}}{\var{SURFN2}} \ \var{SURFAC(3,IFAC)}
\end{array}\right.
\end{equation}
\item [$\Rightarrow$] finalisation du calcul de l'expression totalement
explicite 
 $$\left[ \tens{D}^n\,\left( \grad{R^{\,n}_{ij}} - (\grad R^{\,n}_{ij}\,.\,\vect{n}_{\,lm})\,\vect{n}_{\,lm}\right) \right]\,.\,\vect{n}_{\,lm}$$
\begin{equation}\notag
\begin{array} {ll}
\displaystyle \var{VISCF} = &
 \displaystyle\frac{1}{2} (\ \var{W4(II)} +\ \var{W4(JJ)}) \ \var{GRDPX} \
\var{SURFAC(1,IFAC)})\ + \\
&\\
&  \displaystyle\frac{1}{2} (\ \var{W5(II)} +\ \var{W5(JJ)}) \ \var{GRDPY} \
\var{SURFAC(2,IFAC)})\ + \\
&\\
&  \displaystyle\frac{1}{2} (\ \var{W6(II)} +\ \var{W6(JJ)}) \ \var{GRDPZ} \ \var{SURFAC(3,IFAC)})
\end{array}
\end{equation}
\end{itemize}

\item [$\star$] Mise \`a z\'ero du tableau \var{VISCB}.

\item [$\star$] Appel de \fort{divmas} pour calculer la divergence de~:
 $$\tens{D}^{\,n}\,\left( \grad{R^{\,n}_{ij}} - (\grad R^{\,n}_{ij}\,.\,\vect{n}_{\,lm})\vect{n}_{\,lm}\right)$$ d\'efini aux faces dans \var{VISCF} et \var{VISCB}.

Le r\'esultat est stock\'e dans le tableau \var{W1}.\\
Ajout au second membre \var{SMBR}.\\
$\var{SMBR} = \var{SMBR} + \var{W1}$
\end{itemize}
\item Calcul de la viscosit\'e orthotrope $\gamma^n_{\,lm}$ \`a la face $lm$ de la variable principale
$R^{\,n}_{ij}$\\ 
Ce calcul permet au sous-programme \fort{codits} de compl\'eter le second membre
\var{SMBR} par :
\begin{equation}
\begin{array} {ll}
& \sum\limits_{m\in Vois(l)}
\mu^n_{\,lm}\,\left(\grad{R}^{\,n}_{ij}\,.\,\vect{n}_{\,lm}\right)S_{\,lm}
 + \sum\limits_{m\in Vois(l)} \left(\grad{R}^{\,n}_{ij}
\,.\,\vect{n}_{\,lm}\right)\left[\tens{D}^{\,n}\,\vect{n}_{\,lm}\right]_{\,lm}\,.\,\vect{n}_{\,lm}\
S_{\,lm}\\
& = \sum\limits_{m\in Vois(l)}(\,\mu^n_{\,lm}\, + \,\gamma^n_{\,lm}\,)\,\left(\grad{R}^{\,n}_{ij}\,.\,\vect{n}_{\,lm}\right)S_{\,lm} 
\end{array}
\end{equation}
sans pr\'eciser la nature de la face $lm$, {\it via} l'appel \`a \fort{bilsc2}  et de disposer de la quantit\'e
$(\mu^n_{\,lm}\, + \gamma^n_{\,lm})$ pour construire sa
matrice simplifi\'ee.\\
\begin{itemize}
\item [$\star$] On effectue une boucle d'indice \var{IEL} sur les cellules
$\Omega_l$ :
\begin{itemize}
\item [$\Rightarrow$] $\displaystyle \var{TRRIJ }= \frac{1}{2} (R^{\,n}_{ii})_L $ 
\item [$\Rightarrow$] $\displaystyle \var{RCSTE} = \rho^n_L \ C_S \ \frac{ (R^{\,n}_{ii})_L}{2\,\varepsilon^n_L} $ 
\item [$\Rightarrow$] $\displaystyle \var{W1(IEL)} = \mu^n +\rho^n_L \ C_S \ \frac{
(R^{\,n}_{ii})_L}{2\,\varepsilon^n_L}\ (R^n_{11})_L$
\item [$\Rightarrow$] $\displaystyle \var{W2(IEL)} = \mu^n + \rho^n_L \ C_S \ \frac{ (R^{\,n}_{ii})_L}{2\,\varepsilon^n_L}\ (R^n_{22})_L$
\item [$\Rightarrow$] $\displaystyle \var{W3(IEL)} = \mu^n + \rho^n_L \ C_S \ \frac{ (R^{\,n}_{ii})_L}{2\,\varepsilon^n_L}\ (R^n_{33})_L$
\end{itemize}

\item [$\star$] Appel de \fort{visort} pour calculer la viscosit\'e orthotrope
\footnote{Comme dans le sous-programme \fort{viscfa}, on multiplie la viscosit\'e par
$\displaystyle \frac{S_{\,lm}}{\overline{L'M'}}$, o\`u $S_{\,lm}$ et
$\overline{L'M'}$ repr\'esentent respectivement la surface de la face $lm$ et la
mesure alg\'ebrique du segment reliant les projections des centres des cellules
voisines sur la normale \`a la face. On garde dans ce sous-programme  la possibilit\'e d'interpoler la viscosit\'e aux cellules lin\'eairement ou d'utiliser une moyenne harmonique. La viscosit\'e au bord est celle de la cellule de bord correspondante.}
$ \gamma^n_{\,lm} = (\tens{D}^{\,n}\,\vect{n}_{\,lm}).\vect{n}_{\,lm}$ aux faces $lm$

Le r\'esultat est stock\'e dans les tableaux \var{VISCF} et \var{VISCB}.
\end{itemize}

\item appel de \fort{codits} pour la r\'esolution de l'\'equation de
convection/diffusion/termes sources de la variable $R_{ij}$. Le terme source est
\var{SMBR}, la viscosit\'e \var{VISCF} aux faces purement internes (
resp. \var{VISCB} aux faces de bord) et \var{FLUMAS} le flux de masse interne
 ( resp. \var{FLUMAB} flux de masse au bord). Le r\'esultat est la variable $R_{ij}$ au temps
$n+1$, donc $R^{\,n+1}_{ij}$. Elle est stock\'ee dans le tableau \var{RTP} des
variables mises \`a jour.

\end{itemize}

\etape{Appel de \fort{reseps} pour la r\'esolution de la variable $\varepsilon$}

On d\'ecrit ci-dessous le sous-programme \fort{reseps}, les commentaires du sous-programme \fort{resrij} ne seront pas r\'ep\'et\'es puisque les deux sous-programmes ne diff\`erent que par leurs termes sources.

\begin{itemize}
\item Initialisation \`a z\'ero de \var{SMBR} et \var{ROVSDT}.

\item{Lecture et prise en compte des termes sources utilisateur pour la variable $\varepsilon$ :}

Appel de \fort{ustsri} pour charger les termes sources utilisateurs. Ils sont
stock\'es dans les tableaux suivants :\\
pour la cellule $\Omega_l$ repr\'esent\'ee par $\var{IEL}$ de centre $L$, on a :
\begin{equation}\notag
\left\{\begin{array}{lll}
&\var{ROVSDT(IEL)}&= |\Omega_l| \ \alpha_{\varepsilon}\\
&\var{SMBR(IEL)}&=|\Omega_l| \ \beta_{\varepsilon}\\
\end{array}\right.
\end{equation}
On affecte alors les valeurs ad\'equates au second membre \var{SMBR} et \`a la
diagonale \var{ROVSDT} comme suit :
\begin{equation}\notag
\left\{\begin{array}{lll}
&\var{SMBR(IEL)} &= \var{SMBR(IEL)} +\ |\Omega_l| \ \alpha_{\,\varepsilon} \
\varepsilon^n_L \\
&\var{ROVSDT(IEL)}&= \text{max }(-\ |\Omega_l| \ \alpha_{\,\varepsilon},0)\\
\end{array}\right.
\end{equation}

\item{Calcul du terme source de masse si $\Gamma_L > 0$ :
\begin{equation}\notag
\left\{\begin{array}{lll}
&\displaystyle \var{SMBR(IEL)} = \var{SMBR(IEL)} + |\Omega_l| \ \Gamma_L \
(\varepsilon^{\,in}_L -\varepsilon^n_L) \\
&\displaystyle \var{ROVSDT(IEL)}= \var{ROVSDT(IEL)} + |\Omega_l| \ \Gamma_L
\end{array}\right.
\end{equation}
\item Calcul du terme d'accumulation de masse et du terme instationnaire \\
On stocke $\displaystyle \var{W1}= \int_{\Omega_l}\dive\,(\rho\,\vect{u})\,d\Omega$
calcul\'e par \fort{divmas} \`a l'aide des flux de masse internes et aux bords.\\
On incr\'emente ensuite \var{SMBR} et \var{ROVSDT}.
\begin{equation}\notag
\left\{\begin{array}{lll}
&\var{SMBR(IEL)} &= \var{SMBR(IEL)} + \var{ICONV}\ \varepsilon^n_L\,(\displaystyle
\int_{\Omega_l}\dive\,(\rho\,\vect{u})\ d\Omega) \\
&\var{ROVSDT(IEL)}& = \var{ROVSDT(IEL)} +  \var{ISTAT}\,\displaystyle
\frac{\rho^n_L \ |\Omega_l|}{\Delta t^n} -  \var{ICONV}\ (\displaystyle
\int_{\Omega_l}\dive\,(\rho\,\vect{u})\ d\Omega) \\
\end{array}\right.
\end{equation}

\item Traitement du terme de production
 $\displaystyle \rho\,C_{\varepsilon_1}\,\frac{\varepsilon}{k}\,\mathcal{P}$
 et du terme de dissipation $-\,\displaystyle \rho\,C_{\varepsilon_2}\,\frac{\varepsilon}{k}\,\varepsilon$ \\
pour cela, on effectue une boucle d'indice \var{IEL} sur les cellules $\Omega_l$
de centre $L$ :
\begin{itemize}
\item [$\Rightarrow$] $\displaystyle \var{TRPROD}= \frac{1}{2} (\mathcal{P}^n_{ii})_L = \frac{1}{2} \left[ \var{PRODUC(1,IEL)} +  \var{PRODUC(2,IEL)} +  \var{PRODUC(3,IEL)} \right] $
\item [$\Rightarrow$] $\displaystyle \var{TRRIJ }= \frac{1}{2} (R^n_{ii})_L $ 
\item [$\Rightarrow$] $\displaystyle \var{SMBR(IEL)} = \var{SMBR(IEL)} + \rho^n_L
|\Omega_l| \left[ -C_{\varepsilon_2} \ \frac{2\,(\varepsilon^n_L)^2}{(R^n_{ii})_L} + C_{\varepsilon_1} \ \frac{\varepsilon^n_L}{(R^n_{ii})_L}\ (\mathcal{P}^n_{ii})_L \right] $
\item [$\Rightarrow$] $\displaystyle \var{ROVSDT(IEL)} = \var{ROVSDT(IEL)} + C_{\varepsilon_2} \ \rho^n_L \ |\Omega_l| \ \frac{2\,\varepsilon^n_L}{(R^n_{ii})_L}$
\end{itemize}

\item Appel de \fort{rijthe} pour le calcul des termes de gravit\'e $\mathcal{G}^n_{\varepsilon}$ et ajout dans \var{SMBR}.

$ \var{SMBR} = \var{SMBR} + \mathcal{G}^n_{\varepsilon}$\\
Ce calcul n'a lieu que si $\var{IGRARI()} = 1$.

\item Calcul de la diffusion de $\varepsilon$ \\
 Le terme $\dive \left[\mu\, \grad(\varepsilon) + \tens{A'}\,\grad(\varepsilon)
\right]$ est calcul\'e exactement de la m\^eme mani\`ere que pour les tensions
de Reynolds $R_{ij}$ en rempla\c cant $\tens{A}$ par $\tens{A'}$.

\item Appel de \fort{codits} pour la r\'esolution de l'\'equation de
convection/diffusion/termes sources de la variable principale $\varepsilon$. Le
r\'esultat $\varepsilon^{\,n+1}$ est stock\'e dans le tableau \var{RTP} des
variables mises \`a jour.
}
\end{itemize}

\etape{clippings finaux}
On passe enfin dans le sous-programme  \fort{clprij} pour faire un clipping \'eventuel
des variables $R^{\,n+1}_{ij}$ et $\varepsilon^{\,n+1}$. Le sous-programme
\fort{clprij} est appel\'e\footnote{L'option
$\var{ICLIP} = 1$ consiste en un clipping minimal des variables $R_{ii}$ et
$\varepsilon$ en prenant la valeur absolue de ces variables puisqu'elles ne
peuvent \^etre que positives.} avec $\var{ICLIP} = 2$ . Cette option
\footnote{Quand la valeur des grandeurs $R_{ii}$ ou $\varepsilon$ est
n\'egative, on la remplace par le minimum entre sa valeur absolue et (1,1)
fois la valeur obtenue au pas de temps pr\'ec\'edent.} contient l'option $\var{ICLIP} = 1$  et permet de v\'erifier l'in\'egalit\'e de Cauchy-Schwarz sur les grandeurs extra-diagonales du tenseur $\tens{R}$ (pour $i \neq j$, $|R_{ij}|^2 \le R_{ii} R_{jj}$).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Points \`a traiter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Sauf mention explicite, $\phi$ repr\'esentera une tension de Reynolds ou la dissipation turbulente ($\phi = R_{ij} \ \text{ou} \ \varepsilon$).

\begin{itemize}
\item {La vitesse utilis\'ee pour le calcul de la production est explicite. Est-ce qu'une implicitation peut am\'eliorer la pr\'ecision temporelle de $\phi$ \footnote{Cette remarque peut \^etre g\'en\'eralis\'ee. En effet, peut-on envisager d'actualiser les variables d\'ej\`a r\'esolues (sans r\'eactualiser les variables turbulentes apr\`es leur r\'esolution)? Ceci obligerait \`a modifier les sous-programmes tels que \fort{condli} qui sont appel\'es au d\'ebut de la boucle en temps.} ?}
\item {Dans quelle mesure le terme d'\'echo de paroi est-il valide ? En effet, ce terme est remis en question par certains auteurs.}
\item {On peut envisager la r\'esolution d'un syst\`eme hyperbolique pour les
tensions de Reynolds afin d'introduire un couplage avec le champ de vitesse.}
\item {Le flux au bord \var{VISCB} est annul\'e dans le sous-programme
\fort{vectds}. Peut-on envisager de mettre au bord la valeur de la variable
concern\'ee \`a la cellule de bord correspondant? De m\^eme, il faudrait se
pencher sur les hypoth\`eses sous-jacentes \`a l'annulation des contributions
aux bords de \var{VISCB} lors du calcul de : $$\left[ \tens{D}^n\,\left( \grad{R^{\,n}_{ij}} - (\grad R^{\,n}_{ij}\,.\,\vect{n}_{\,lm})\,\vect{n}_{\,lm}\right) \right]\,.\,\vect{n}_{\,lm}.$$}
\item {Un probl\`eme de pond\'eration appara\^\i t plus g\'en\'eralement. Si on prend la partie explicite de $\tens{D}\,\grad(\phi)$, la pond\'eration aux faces internes utilise le coefficient $\displaystyle\frac{1}{2}$ avec pond\'eration s\'epar\'ee de $\tens{D}$ et $\grad(\phi)$, alors que pour $\tens{E}\,\grad(\phi)$, on calcule d'abord ce terme aux cellules pour ensuite l'interpoler lin\'eairement aux faces \footnote{Cette interpolation se fait dans \fort{vectds} avec des coefficients de pond\'eration aux faces.}. Ceci donne donc deux types d'interpolations pour des termes de m\^eme nature.}
\item {On laisse la possibilit\'e dans \fort{visort} d'utiliser une moyenne
harmonique aux faces. Est-ce que ceci est valable puisque les interpolations
utilis\'ees lors du calcul de la partie explicite de $\tens{A}\,\grad{\phi}$
sont des moyennes arithm\'etiques ?}
\item {Les techniques adopt\'ees lors du clipping sont \`a revoir.}
\item {On utilise dans le cadre du mod\`ele $\displaystyle R_{ij}-\varepsilon $ une semi-implicitation de termes comme $\displaystyle \phi_{ij,1}$ ou $\displaystyle -\rho\,C_{\varepsilon_2}\,\frac{\varepsilon}{k}\,\varepsilon$. On peut envisager le m\^eme type d'implicitation dans \fort{turbke} m\^eme en pr\'esence du couplage $\displaystyle k-\varepsilon$.}
\item L'adoption d'une r\'esolution d\'ecoupl\'ee fait perdre l'invariance par rotation.
\item La formulation et l'implantation des conditions aux limites de paroi
devront \^etre v\'erifi\'ees. En effet, il semblerait que, dans certains cas, des ph\'enom\`enes
``oscillatoires'' apparaissent, sans qu'il soit ais\'e d'en d\'eterminer la cause.
\item L'implicitation partielle (du fait de la r\'esolution d\'ecoupl\'ee) des
conditions aux limites conduit souvent \`a des calculs instables. Il
conviendrait d'en conna\^\i tre la raison. L'implicitation partielle avait
\'et\'e mise en \oe uvre afin de tenter d'utiliser un pas de temps plus grand
dans le cas de jets axisym\'etriques en particulier.

\end{itemize}
