%-------------------------------------------------------------------------------

% This file is part of Code_Saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2021 EDF S.A.
%
% This program 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.
%
% This program 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
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.

%-------------------------------------------------------------------------------

\programme{turbke}\label{ap:turbke}

\hypertarget{turbke}{}

\vspace{1cm}
%-------------------------------------------------------------------------------
\section*{Fonction}
%-------------------------------------------------------------------------------
Le but de ce sous-programme est de r\'esoudre le syst\`eme des \'equations de
$k$ et $\varepsilon$ de mani\`ere semi-coupl\'ee.\\
Le syst\`eme d'\'equations r\'esolu est le suivant :

\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho\frac{\partial k}{\partial t} +
\dive\left[\rho \vect{u}\,k-(\mu+\frac{\mu_t}{\sigma_k})\grad{k}\right] =
\mathcal{P}+\mathcal{G}-\rho\varepsilon+k\dive(\rho\vect{u})
+\Gamma(k_i-k)\\
\multicolumn{1}{c}{+\alpha_k k +\beta_k}\\
\displaystyle
\rho\frac{\partial \varepsilon}{\partial t} +
\dive\left[\rho \vect{u}\,\varepsilon-
(\mu+\frac{\mu_t}{\sigma_\varepsilon})\grad{\varepsilon}\right] =
C_{\varepsilon_1}\frac{\varepsilon}{k}\left[\mathcal{P}
+(1-C_{\varepsilon_3})\mathcal{G}\right]
-\rho C_{\varepsilon_2}\frac{\varepsilon^2}{k}
+\varepsilon\dive(\rho\vect{u})\\
\multicolumn{1}{c}{+\Gamma(\varepsilon_i-\varepsilon)
+\alpha_\varepsilon \varepsilon +\beta_\varepsilon}
\end{array}\right.
\end{equation}

$\mathcal{P}$ est le terme de production par cisaillement moyen :
\begin{displaymath}
\begin{array}{rcl}
\mathcal{P} & = & \displaystyle -\rho R_{ij}\frac{\partial u_i}{\partial x_j}
= -\left[-\mu_t \left(\frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i}\right)
+\frac{2}{3}\mu_t\frac{\partial u_k}{\partial x_k}\delta_{ij}
+\frac{2}{3}\rho k\delta_{ij}\right]
\frac{\partial u_i}{\partial x_j}\\
& = & \displaystyle \mu_t \left(\frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}
-\frac{2}{3}\mu_t(\dive\vect{u})^2-\frac{2}{3}\rho k \dive(\vect{u})\\
& = & \displaystyle \mu_t \left[
2\left(\frac{\partial u}{\partial x}\right)^2+
2\left(\frac{\partial v}{\partial y}\right)^2+
2\left(\frac{\partial w}{\partial z}\right)^2+
\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2+
\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)^2+
\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2
\right]\\
\multicolumn{3}{r}%
{\displaystyle -\frac{2}{3}\mu_t(\dive\vect{u})^2-\frac{2}{3}\rho k \dive(\vect{u})}
\end{array}
\end{displaymath}

$\mathcal{G}$ est le terme de production par gravit\'e :
$\displaystyle
\mathcal{G}=-\frac{1}{\rho}\frac{\mu_t}{\sigma_t}
\frac{\partial\rho}{\partial x_i}g_i$

La viscosit\'e turbulente est
$\displaystyle \mu_t=\rho C_\mu\frac{k^2}{\varepsilon}$.

Les constantes sont :\\
$C_\mu=0,09$ ;
$C_{\varepsilon_2}=1,92$ ; $\sigma_k=1$ ; $\sigma_\varepsilon=1,3$\\
$C_{\varepsilon_3}=0$ si $\mathcal{G}\geqslant0$ (stratification instable) et
$C_{\varepsilon_3}=1$ si $\mathcal{G}\leqslant0$ (stratification stable).

$\Gamma$ est un \'eventuel terme source de masse (tel que l'\'equation de
conservation de masse devienne
$\displaystyle \frac{\partial \rho}{\partial t}+\dive(\rho\vect{u})=\Gamma$).
$\varphi_i$ ($\varphi=k$ ou $\varepsilon$) est la valeur de $\varphi$
associ\'ee \`a la masse inject\'ee ou retir\'ee. Dans le cas o\`u on retire de
la masse ($\Gamma<0$), on a forc\'ement $\varphi_i=\varphi$. De m\^eme, quand on
injecte de la masse, on sp\'ecifie souvent aussi $\varphi_i=\varphi$. Dans ces
deux cas, le terme dispara\^\i t de l'\'equation. Dans la suite du document, on
qualifiera d'{\em injection forc\'ee} les cas o\`u on a $\Gamma>0$ et
$\varphi_i\ne\varphi$.

$\alpha_k$, $\beta_k$, $\alpha_\varepsilon$, $\beta_\varepsilon$ sont des termes
sources utilisateur \'eventuels, conduisant \`a une implicitation partielle, impos\'es le cas
\'ech\'eant par le sous-programme \fort{ustske}.

See the \doxygenfile{turbke_8f90.html}{programmers reference of the dedicated subroutine} for further details.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discr\'etisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
La r\'esolution se fait en trois \'etapes, afin de coupler partiellement les
deux variables $k$ et $\varepsilon$. Pour simplifier, r\'e\'ecrivons le
syst\`eme de la fa\c con suivante :

\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho\frac{\partial k}{\partial t} =
D(k) + S_k(k,\varepsilon)+k\dive(\rho\vect{u})+\Gamma(k_i-k)+\alpha_k k +\beta_k\\
\displaystyle
\rho\frac{\partial \varepsilon}{\partial t}  =
D(\varepsilon) + S_\varepsilon(k,\varepsilon)
+\varepsilon\dive(\rho\vect{u})
+\Gamma(\varepsilon_i-\varepsilon)+\alpha_\varepsilon \varepsilon +\beta_\varepsilon
\end{array}\right.
\end{equation}

$D$ est l'op\'erateur de convection/diffusion.
$S_k$ (resp. $S_\varepsilon$) est le terme source de $k$ (resp. $\varepsilon$).

\minititre{Premi\`ere phase : bilan explicite}

On r\'esout le bilan explicite :
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k_e-k^{(n)}}{\Delta t} =
D(k^{(n)}) + S_k(k^{(n)},\varepsilon^{(n)})
+k^{(n)}\dive(\rho\vect{u})+\Gamma(k_i-k^{(n)})+\alpha_k k^{(n)} +\beta_k\\
\displaystyle
\rho^{(n)}\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}  =
D(\varepsilon^{(n)}) + S_\varepsilon(k^{(n)},\varepsilon^{(n)})
+\varepsilon^{(n)}\dive(\rho\vect{u})
+\Gamma(\varepsilon_i-\varepsilon^{(n)})
+\alpha_\varepsilon \varepsilon^{(n)} +\beta_\varepsilon
\end{array}\right.
\end{equation}

(le terme en $\Gamma$ n'est pris en compte que dans le cas de l'injection forc\'ee)

\minititre{Deuxi\`eme phase : couplage des termes sources}

On implicite les termes sources de mani\`ere coupl\'ee :
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t} =
D(k^{(n)}) + S_k(k_{ts},\varepsilon_{ts})
+k^{(n)}\dive(\rho\vect{u})+\Gamma(k_i-k^{(n)})+\alpha_k k^{(n)} +\beta_k\\
\displaystyle
\rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}  =
D(\varepsilon^{(n)}) + S_\varepsilon(k_{ts},\varepsilon_{ts})
+\varepsilon^{(n)}\dive(\rho\vect{u})
+\Gamma(\varepsilon_i-\varepsilon^{(n)})
+\alpha_\varepsilon \varepsilon^{(n)} +\beta_\varepsilon
\end{array}\right.
\end{equation}
soit
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t} =
\rho^{(n)}\frac{k_e-k^{(n)}}{\Delta t}
+S_k(k_{ts},\varepsilon_{ts})-S_k(k^{(n)},\varepsilon^{(n)})\\
\displaystyle
\rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}  =
\rho^{(n)}\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}
+S_\varepsilon(k_{ts},\varepsilon_{ts})-S_\varepsilon(k^{(n)},\varepsilon^{(n)})
\end{array}\right.
\end{equation}

Le terme en $\dive(\rho\vect{u})$ n'est pas implicit\'e car il est li\'e au
terme $D$ pour assurer que la matrice d'implicitation sera \`a diagonale
dominante. Le terme en $\Gamma$ et les termes sources utilisateur ne sont
pas implicit\'es non plus, mais ils le seront dans la troisi\`eme phase.

Et on \'ecrit (pour $\varphi=k$ ou $\varepsilon$)
\begin{equation}
S_\varphi(k_{ts},\varepsilon_{ts})-S_\varphi(k^{(n)},\varepsilon^{(n)})
=(k_{ts}-k^{(n)})
\left.\frac{\partial S_\varphi}{\partial k}\right|_{k^{(n)},\varepsilon^{(n)}}
+(\varepsilon_{ts}-\varepsilon^{(n)})
\left.\frac{\partial S_\varphi}{\partial \varepsilon}\right|_{k^{(n)},\varepsilon^{(n)}}
\end{equation}

On r\'esout donc finalement le syst\`eme $2\times 2$ :
\begin{equation}
\left(\begin{array}{cc}
\displaystyle \frac{\rho^{(n)}}{\Delta t}
-\left.\frac{\partial S_k}{\partial k}\right|_{k^{(n)},\varepsilon^{(n)}}
&\displaystyle
-\left.\frac{\partial S_k}{\partial \varepsilon}\right|_{k^{(n)},\varepsilon^{(n)}}\\
\displaystyle
-\left.\frac{\partial S_\varepsilon}{\partial k}\right|_{k^{(n)},\varepsilon^{(n)}}
&\displaystyle
\displaystyle \frac{\rho^{(n)}}{\Delta t}
-\left.\frac{\partial S_\varepsilon}{\partial \varepsilon}\right|_{k^{(n)},\varepsilon^{(n)}}
\end{array}\right)
\left(\begin{array}{c}
(k_{ts}-k^{(n)})\\(\varepsilon_{ts}-\varepsilon^{(n)})
\end{array}\right)
=\left(\begin{array}{c}
\displaystyle\rho^{(n)}\frac{k_e-k^{(n)}}{\Delta t}\\
\displaystyle\rho^{(n)}\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}
\end{array}\right)
\end{equation}

\vspace*{0.2cm}

\minititre{Troisi\`eme phase : implicitation de la convection/diffusion}

On r\'esout le syst\`eme :
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k^{(n+1)}-k^{(n)}}{\Delta t} =
D(k^{(n+1)}) + S_k(k_{ts},\varepsilon_{ts})
+k^{(n+1)}\dive(\rho\vect{u})+\Gamma(k_i-k^{(n+1)})\\
\multicolumn{1}{r}{+\alpha_k k^{(n+1)} +\beta_k}\\
\displaystyle
\rho^{(n)}\frac{\varepsilon^{(n+1)}-\varepsilon^{(n)}}{\Delta t}  =
D(\varepsilon^{(n+1)}) + S_\varepsilon(k_{ts},\varepsilon_{ts})
+\varepsilon^{(n+1)}\dive(\rho\vect{u})
+\Gamma(\varepsilon_i-\varepsilon^{(n+1)})\\
\multicolumn{1}{r}{+\alpha_\varepsilon \varepsilon^{(n+1)} +\beta_\varepsilon}
\end{array}\right.
\end{equation}
soit
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k^{(n+1)}-k^{(n)}}{\Delta t} =
D(k^{(n+1)})-D(k^{(n)})+\rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t}
+(k^{(n+1)}-k^{(n)})\dive(\rho\vect{u})\\
\multicolumn{1}{r}{-\Gamma(k^{(n+1)}-k^{(n)})+\alpha_k(k^{(n+1)}-k^{(n)})}\\
\displaystyle
\rho^{(n)}\frac{\varepsilon^{(n+1)}-\varepsilon^{(n)}}{\Delta t}  =
D(\varepsilon^{(n+1)})-D(\varepsilon^{(n)})
+\rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}
+(\varepsilon^{(n+1)}-\varepsilon^{(n)})\dive(\rho\vect{u})\\
\multicolumn{1}{r}{-\Gamma(\varepsilon^{(n+1)}-\varepsilon^{(n)})
+\alpha_\varepsilon(\varepsilon^{(n+1)}-\varepsilon^{(n)})}
\end{array}\right.
\end{equation}

Le terme en $\Gamma$ n'est l\`a encore pris en compte que dans le cas de
l'injection forc\'ee. Le terme en $\alpha$ n'est pris en compte que si $\alpha$ est
n\'egatif, pour \'eviter d'affaiblir la diagonale de la matrice qu'on va
inverser.


\minititre{Pr\'ecisions sur le couplage}
Lors de la phase de couplage, afin de privil\'egier la stabilit\'e et la
r\'ealisabilit\'e du r\'esultat, tous les termes ne sont pas pris en
compte. Plus pr\'ecis\'ement, on peut \'ecrire :

\begin{equation}
\left\{\begin{array}{l}
\displaystyle
S_k =
\rho C_\mu\frac{k^2}{\varepsilon}\left(\tilde{\mathcal{P}}+\tilde{\mathcal{G}}\right)
-\frac{2}{3}\rho k \dive(\vect{u})
-\rho\varepsilon\\
\displaystyle
S_\varepsilon =
\rho C_{\varepsilon_1} C_\mu k\left(\tilde{\mathcal{P}}
+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}\right)
-\frac{2}{3}C_{\varepsilon_1}\rho \varepsilon \dive(\vect{u})
-\rho C_{\varepsilon_2}\frac{\varepsilon^2}{k}
\end{array}\right.
\end{equation}

en notant
$\displaystyle\tilde{\mathcal{P}}
= \left(\frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}
-\frac{2}{3}(\dive\vect{u})^2$\\
et
$\displaystyle\tilde{\mathcal{G}}
= -\frac{1}{\rho\sigma_t}
\frac{\partial\rho}{\partial x_i}g_i$

On a donc en th\'eorie
\begin{equation}
\left\{\begin{array}{l}
\displaystyle \frac{\partial S_k}{\partial k}=
2\rho C_\mu\frac{k}{\varepsilon}\left(\tilde{\mathcal{P}}+\tilde{\mathcal{G}}\right)
-\frac{2}{3}\rho \dive(\vect{u})\\
\displaystyle \frac{\partial S_k}{\partial \varepsilon}= -\rho\\
\displaystyle \frac{\partial S_\varepsilon}{\partial k}=
\rho C_{\varepsilon_1} C_\mu \left(\tilde{\mathcal{P}}
+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}\right)
+\rho C_{\varepsilon_2}\frac{\varepsilon^2}{k^2}\\
\displaystyle \frac{\partial S_\varepsilon}{\partial \varepsilon}=
-\frac{2}{3}C_{\varepsilon_1}\rho \dive(\vect{u})
-2\rho C_{\varepsilon_2}\frac{\varepsilon}{k}
\end{array}\right.
\end{equation}

En pratique, on va chercher \`a assurer $k_{ts}>0$ et $\varepsilon_{ts}>0$. En se
basant sur un calcul simplifi\'e, ainsi que sur le retour d'exp\'erience
d'ESTET, on montre qu'il est pr\'ef\'erable de ne pas prendre en compte
certains termes. Au final, on r\'ealise le couplage suivant :

\begin{equation}
\left(\begin{array}{cc}
A_{11}&A_{12}\\
A_{21}&A_{22}
\end{array}\right)
\left(\begin{array}{c}
(k_{ts}-k^{(n)})\\(\varepsilon_{ts}-\varepsilon^{(n)})
\end{array}\right)
=\left(\begin{array}{c}
\displaystyle\frac{k_e-k^{(n)}}{\Delta t}\\
\displaystyle\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}
\end{array}\right)
\end{equation}
avec
\begin{equation}
\left\{\begin{array}{l}
\displaystyle A_{11}=\frac{1}{\Delta t}
-2 C_\mu\frac{k^{(n)}}{\varepsilon^{(n)}}
\Min\left[\left(\tilde{\mathcal{P}}+\tilde{\mathcal{G}}\right),0\right]
+\frac{2}{3}\Max\left[\dive(\vect{u}),0\right]\\
\displaystyle A_{12}= 1\\
\displaystyle A_{21}=
- C_{\varepsilon_1} C_\mu \left(\tilde{\mathcal{P}}
+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}\right)
- C_{\varepsilon_2}\left(\frac{\varepsilon^{(n)}}{k^{(n)}}\right)^2\\
\displaystyle A_{22}=\frac{1}{\Delta t}
+\frac{2}{3}C_{\varepsilon_1}\Max\left[\dive(\vect{u}),0\right]
+2 C_{\varepsilon_2}\frac{\varepsilon^{(n)}}{k^{(n)}}
\end{array}\right.
\end{equation}

(par d\'efinition de $C_{\varepsilon_3}$,
$\tilde{\mathcal{P}}+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}$
est toujours positif)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\etape{Calcul du terme de production}
On appelle trois fois \fort{grdcel} pour calculer les gradients de $u$, $v$ et
$w$. Au final, on a \\
$\displaystyle \var{TINSTK}=
2\left(\frac{\partial u}{\partial x}\right)^2+
2\left(\frac{\partial v}{\partial y}\right)^2+
2\left(\frac{\partial w}{\partial z}\right)^2+
\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2+
\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)^2+
\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2$\\
et\\
$\displaystyle \var{DIVU}=
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}
+\frac{\partial w}{\partial z}$

(le terme $div(\vect{u})$ n'est pas calcul\'e par \fort{divmas}, pour
correspondre exactement \`a la trace du tenseur des d\'eformations qui est
calcul\'e pour la production)


\etape{Lecture des termes sources utilisateur}
Appel de \fort{cs\_user\_turbulence\_source\_terms} pour charger les termes sources utilisateurs. Ils sont
stock\'es dans les tableaux suivants :\\
$\var{W7}=\Omega\beta_k$\\
$\var{W8}=\Omega\beta_\varepsilon$\\
$\var{DAM}=\Omega\alpha_k$\\
$\var{W9}=\Omega\alpha_\varepsilon$

Puis on ajoute le terme en $(div\vect{u})^2$ \`a \var{TINSTK}. On a donc \\
$\var{TINSTK}=\tilde{\mathcal{P}}$

\etape{Calcul du terme de gravit\'e}
Calcul uniquement si $\var{IGRAKE}=1$.\\
On appelle \fort{grdcel} pour \var{ROM}, avec comme conditions aux limites
$\var{COEFA}=\var{ROMB}$ et \mbox{$\var{COEFB}=\var{VISCB}=0$}.\\
$\var{PRDTUR}=\sigma_t$ est mis \`a 1 si on n'a pas de scalaire temp\'erature.

$\tilde{\mathcal{G}}$ est calcul\'e et les termes sources sont mis \`a jour :\\
$\var{TINSTK}=\tilde{\mathcal{P}}+\tilde{\mathcal{G}}$\\
$\var{TINSTE}=\tilde{\mathcal{P}}+\Max\left[\tilde{\mathcal{G}},0\right]
=\tilde{\mathcal{P}}+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}$

Si $\var{IGRAKE}=0$, on a simplement\\
$\var{TINSTK}=\var{TINSTE}=\tilde{\mathcal{P}}$

\etape{Calcul du terme d'accumulation de masse}
On stocke
$\displaystyle \var{W1}=\Omega\dive(\rho\vect{u})$
calcul\'e par \fort{divmas} (pour correspondre aux termes de convection de la
matrice).

\etape{Calcul des termes sources explicites}
On affecte les termes sources explicites de $k$ et $\varepsilon$ pour la
premi\`ere \'etape.\\
$\displaystyle\var{SMBRK}=\Omega\left(\mu_t(\tilde{\mathcal{P}}+\tilde{\mathcal{G}})
-\frac{2}{3}\rho^{(n)} k^{(n)}\dive{\vect{u}}
-\rho^{(n)} \varepsilon^{(n)}\right)+\Omega k^{(n)}\dive(\rho\vect{u})$\\
$\displaystyle\var{SMBRE}=\Omega\frac{\varepsilon^{(n)}}{k^{(n)}}
\left(C_{\varepsilon_1}\left(
\mu_t(\tilde{\mathcal{P}}+(1-C_{\varepsilon_3})\tilde{\mathcal{G}})
-\frac{2}{3}\rho^{(n)} k^{(n)}\dive{\vect{u}}\right)
-C_{\varepsilon_2}\rho^{(n)}\varepsilon^{(n)}\right)
+\Omega\varepsilon^{(n)}\dive(\rho\vect{u})$

soit $\var{SMBRK}=\Omega S_k^{(n)}+\Omega k^{(n)}\dive(\rho\vect{u})$
et $\var{SMBRE}=\Omega S_\varepsilon^{(n)}+\Omega\varepsilon^{(n)}\dive(\rho\vect{u})$.


\etape{Calcul des termes sources utilisateur}
On ajoute les termes sources utilisateur explicites \`a \var{SMBRK} et
\var{SMBRE}, soit :\\
$\var{SMBRK}=\Omega S_k^{(n)}+\Omega k^{(n)}\dive(\rho\vect{u})+\Omega\alpha_k k^{(n)} +\Omega\beta_k$\\
$\var{SMBRE}=\Omega S_\varepsilon^{(n)}+\Omega\varepsilon^{(n)}\dive(\rho\vect{u})
+\Omega\alpha_\varepsilon \varepsilon^{(n)} +\Omega\beta_\varepsilon$

Les tableaux \var{W7} et \var{W8} sont lib\'er\'es, \var{DAM} et \var{W9} sont
conserv\'es pour \^etre utilis\'es dans la troisi\`eme phase de r\'esolution.

\etape{Calcul des termes de convection/diffusion explicites}
\fort{bilsc2} est appel\'e deux fois, pour $k$ et pour $\varepsilon$, afin
d'ajouter \`a \var{SMBRK} et \var{SMBRE} les termes de convection/diffusion
explicites $D(k^{(n)})$ et $D(\varepsilon^{(n)})$. Ces termes sont d'abord
stock\'es dans \var{W7} et \var{W8}, pour \^etre conserv\'es et r\'eutilis\'es
dans la troisi\`eme phase de r\'esolution.


\etape{Termes source de masse}
Dans le cas d'une injection forc\'ee de mati\`ere, on passe deux fois dans
\fort{catsma} pour ajouter les termes en
$\Omega \Gamma (k_i-k^{(n)})$ et
$\Omega \Gamma (\varepsilon_i-\varepsilon^{(n)})$ \`a \var{SMBRK} et
\var{SMBRE}. La partie implicite ($\Omega\Gamma$) est stock\'ee dans les
variables \var{W2} et \var{W3}, qui seront utilis\'ees lors de la troisi\`eme
phase (les deux variables sont bien n\'ecessaires, au cas o\`u on aurait une
injection forc\'ee sur $k$ et pas sur $\varepsilon$, par exemple).

\etape{Fin de la premi\`ere phase}
Ceci termine la premi\`ere phase. On a \\
$\displaystyle \var{SMBRK}=\Omega \rho^{(n)}\frac{k_e-k^{(n)}}{\Delta t}$\\
$\displaystyle \var{SMBRE}=\Omega \rho^{(n)}\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}$

\etape{Phase de couplage}
(uniquement si $\var{IKECOU}=1$)

On renormalise \var{SMBRK} et \var{SMBRE} qui deviennent les seconds membres du
syst\`eme de couplage.\\
$\displaystyle \var{SMBRK}=\frac{1}{\Omega\rho^{(n)}}\var{SMBRK}
=\frac{k_e-k^{(n)}}{\Delta t}$\\
$\displaystyle \var{SMBRE}=\frac{1}{\Omega\rho^{(n)}}\var{SMBRE}
=\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}$\\
et $\displaystyle \var{DIVP23}=\frac{2}{3}\Max\left[\dive(\vect{u}),0\right]$.

On remplit la matrice de couplage\\
$\displaystyle \var{A11}=\frac{1}{\Delta t}
-2 C_\mu\frac{k^{(n)}}{\varepsilon^{(n)}}
\Min\left[\left(\tilde{\mathcal{P}}+\tilde{\mathcal{G}}\right),0\right]
+\frac{2}{3}\Max\left[\dive(\vect{u}),0\right]$\\
$\displaystyle \var{A12}= 1$\\
$\displaystyle \var{A21}=
- C_{\varepsilon_1} C_\mu \left(\tilde{\mathcal{P}}
+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}\right)
- C_{\varepsilon_2}\left(\frac{\varepsilon^{(n)}}{k^{(n)}}\right)^2$\\
$\displaystyle \var{A22}=\frac{1}{\Delta t}
+\frac{2}{3}C_{\varepsilon_1}\Max\left[\dive(\vect{u}),0\right]
+2 C_{\varepsilon_2}\frac{\varepsilon^{(n)}}{k^{(n)}}$

On inverse le syst\`eme $2\times 2$, et on r\'ecup\`ere :\\
$\displaystyle \var{DELTK}=k_{ts}-k^{(n)}$\\
$\displaystyle \var{DELTE}=\varepsilon_{ts}-\varepsilon^{(n)}$

\etape{Fin de la deuxi\`eme phase}
On met \`a jour les variables \var{SMBRK} et \var{SMBRE}.\\
$\displaystyle \var{SMBRK}=\Omega \rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t}$\\
$\displaystyle \var{SMBRE}=
\Omega \rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}$

Dans le cas o\`u on ne couple pas ($\var{IKECOU}=0$), ces deux variables gardent
la m\^eme valeur qu'\`a la fin de la premi\`ere \'etape.

\etape{Calcul des termes implicites}
On retire \`a \var{SMBRK} et \var{SMBRE} la partie en convection diffusion au
temps $n$, qui \'etait stock\'ee dans \var{W7} et \var{W8}.\\
$\displaystyle \var{SMBRK}=\Omega \rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t}
-\Omega D(k^{(n)})$\\
$\displaystyle \var{SMBRE}=
\Omega \rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}
-\Omega D(\varepsilon^{(n)})$


On calcule les termes implicites, hors convection/diffusion, qui correspondent
\`a la diagonale de la matrice.\\
$\displaystyle \var{TINSTK}=\frac{\Omega \rho^{(n)}}{\Delta t}
-\Omega\dive(\rho\vect{u})+\Omega\Gamma+\Omega\Max[-\alpha_k,0]$\\
$\displaystyle \var{TINSTE}=\frac{\Omega \rho^{(n)}}{\Delta t}
-\Omega\dive(\rho\vect{u})+\Omega\Gamma+\Omega\Max[-\alpha_\varepsilon,0]$\\
($\Gamma$ n'est pris en compte qu'en injection forc\'ee, c'est-\`a-dire qu'il
est forc\'ement positif et ne risque pas d'affaiblir la diagonale de la matrice).

\etape{R\'esolution finale}
On passe alors deux fois dans \fort{codits}, pour $k$ et $\varepsilon$,
pour r\'esoudre les \'equations du type :

$\var{TINST}\times(\varphi^{(n+1)}-\varphi^{(n)}) = D(\varphi^{(n+1)})+\var{SMBR}$.

\etape{clipping final}
On passe enfin dans la routine \fort{clipke} pour faire un clipping \'eventuel
de $k^{(n+1)}$ et $\varepsilon^{(n+1)}$.