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

