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

Le terme convectif en $\dive(\underline{u} \otimes \rho\,\underline{u})$
introduit une non lin\'earit\'e et un couplage des composantes de la vitesse
$\vect{u}$ dans l'équation (\ref{Base_Preduv_eqqdm}). Une lin\'earisation et un d\'ecouplage
des trois composantes de la 
vitesse sont r\'ealis\'es lors de la discr\'etisation de cette \'etape de
pr\'ediction.\\
En effet, soit :
\begin{equation}
\vect{\widetilde{u}}= \vect{u}^n + \delta \vect{u} 
\end{equation}
La contribution exacte du terme convectif \`a prendre en compte dans cette
\'etape de pr\'ediction serait :\\
\begin{equation}\label{Base_Preduv_Conv_exact}
\begin{array}{ll}
\dive(\vect{\widetilde{u}} \otimes \rho\,\vect{\widetilde{u}}) =
\dive(\vect{u}^{n} \otimes \rho\,\vect{u}^{n}) + \dive(\delta \vect{u} \otimes
\rho\,\vect{u}^{n}) +  \underbrace { \dive(\vect{u}^{n} \otimes
\rho\,\delta \vect{u})}_{\text {terme couplant lin\'eaire}} +  \underbrace { \dive(\delta \vect{u} \otimes
\rho\,\delta \vect{u})}_{\text {terme couplant et non lin\'eaire}}\\
\end{array} 
\end{equation}
Les deux derniers termes de l'expression (\ref{Base_Preduv_Conv_exact}) sont {\em a priori} négligés
de manière à obtenir un syst\`eme en vitesse qui soit d\'ecoupl\'e et donc,
éviter l'inversion d'une matrice pouvant \^etre de tr\`es grande taille. Ces
deux termes peuvent néanmoins être pris en compte de manière plus ou moins
approchée par extrapolation explicite du flux de masse en $n+\theta_F$ (pour le
terme couplant linéaire provenant de la convection de $\vect{u}^{n}$ par $\delta
\vect{u}$) et utilisation d'un point-fixe par sous itération sur le sous
programme \fort{navsto} (pour le terme non-linéaire, en spécifiant $\var{NTERUP}>1$).

L'équation (\ref{Base_Preduv_eqqdm}) est discrétisée au temps $n+\theta$ à l'aide d'un
$\theta$-schéma, et le tenseur des pertes de charges décomposé en une partie
diagonale $\tens{K}_{d}$ et une extradiagonale $\tens{K}_{e}$ (soit
 $\tens{K}_{pdc}=\tens{K}_{d}+\tens{K}_{e}$).\\
$\bullet$ La pression est supposée connue en $n-1+\theta$ (décalage temporel
pression-vitesse) et le gradient naturellement calculé à cet instant.\\ 
$\bullet$ Les termes sources de viscosité secondaire, de gradient transpos\'e,
ceux provenant du modèle de turbulence\footnote{excepté $\dive (\mu_t\ (\ggrad
\underline {u}))$}, $\rho\,\tens{K}_{\,e}\ \underline{u}$, $(\rho -\rho_0)
\underline {g}$ ainsi que $\underline{T}_{s}^{\,exp}$ et
$\Gamma\,\underline{u}_{\,i}$ sont pris de manière explicite au temps $n$, ou
extrapolés suivant le schéma en temps choisi pour les propriétés physique et les
termes sources.\\ 
$\bullet$ Les termes sources $\underline{u}\,\,\dive (\rho\,\underline {u})$,
$\Gamma\,\,\underline{u}$, $T_{s}^{\,imp}\,\,\underline{u}$ et
$-\rho\,\tens{K}_{\,d}\,\,\underline{u}$ sont implicités est calculés à
l'instant $n+\theta$.\\ 
$\bullet$ Le terme de diffusion $\dive (\mu_{\,tot}\,\ggrad \underline{u})$ est
implicité : la vitesse est prise à l'instant $n+\theta$ et la viscosité
explicitée ou extrapolée.\\ 
$\bullet$ Enfin, le terme de convection en $\dive(\,\underline{u} \otimes
(\rho\underline{u})\,)$ est implicité : la composante résolue de la vitesse est
prise en $n+\theta$, et le flux de masse, explicité, ou extrapolé en
$n+\theta_F$. 

Par souci de clarté, on suppose, en l'absence d'indication, que les propriétes
physiques $\Phi$ ($\rho,\,\mu_{tot},\,...$) et le flux de masse
$(\rho\underline{u})$ sont pris respectivement aux instants $n+\theta_\Phi$ et
$n+\theta_F$, où $\theta_\Phi$ et $\theta_F$ dépendent des schémas en temps
spécifiquement utilisés pour ces grandeurs\footnote{cf. \fort{introd}}. 

La discrétisation temporelle de l'équation (\ref{Base_Preduv_eqqdm}) s'écrit alors comme suit : 

\begin{equation}\label{Base_Preduv_eq_di1}
 \begin{array}{c}
\displaystyle \rho\,\ \frac{ \underline {\widetilde{u}}^{n+1} -\underline {u}^{n} }
{\Delta t} + \dive(\,\underline{\widetilde{u}}^{n+\theta} \otimes (\rho\underline{u})\,) -\dive
(\mu_{\,tot}\,\ggrad \underline{\widetilde{u}}^{n+\theta}) =
\\
\displaystyle
 - \grad p^{n-1+\theta} + \dive (\rho\,\underline {u})\,\underline{\widetilde{u}}^{n+\theta} +(\Gamma\,\underline{u}_{\,i})^{n+\theta_S}-\Gamma^n\,\,\underline{\widetilde{u}}^{n+\theta}
\\
\begin{array}{c}
\displaystyle
- \rho\,\tens{K}_{\,d}^{n}\,\,\underline{\widetilde{u}}^{n+\theta} - (\rho\,\tens{K}_{\,e}\ \underline{u})^{n+\theta_S} + (\underline{T}_{s}^{\,exp})^{\,n+\theta_S} + T_{s}^{\,imp}\,\,\underline{\widetilde{u}}^{n+\theta}
\\
\displaystyle
+[\dive (\mu_{\,tot}\,^t\ggrad \underline {u})]^{n+\theta_S}-\frac {2} {3}[\,\grad (\mu_{\,tot}\,\dive \underline {u})]^{n+\theta_S} + (\rho -\rho_0) \underline {g}
 - (\underline{turb})^{n+\theta_S}
\end{array}
\end{array}
\end{equation}
o\`u, par souci de simplification, on a pos\'e :
\begin{equation}
\mu_{\,tot}=
\begin{cases}
\mu+\mu_t & \text{pour les modèles à viscosité turbulente ou en LES}, \\
\mu & \text{pour les modèles au second ordre ou en laminaire}
\end{cases} \ 
\end{equation}
\\
et :
\begin{equation}
\underline{turb}^{n}=
\begin{cases}
\displaystyle\frac {2}{3}\grad (\rho^{n}\,k^{n}) & \text{pour les modèles à viscosité turbulente}, \\
\dive(\rho^{n}\,\tens{R}^n) & \text{pour les modèles au second ordre},\\
0 & \text{en laminaire ou en LES}\\
\end{cases}
\end{equation}
Par analogie avec l'écriture du $\theta$-schéma pour une variable scalaire, $\,
\underline {\widetilde{u}}^{n+\theta}$ est interpolée à partir de la vitesse
prédite $\underline {\widetilde{u}}^{n+1}$ de la mani\`ere suivante\footnote{si
$\theta=1/2$, ou qu'une extrapolation est utilisée, l'ordre 2 n'est obtenu que si
le pas de temps $\Delta t$ est uniforme en temps et en espace.}~: 
\begin{equation}
\underline {\widetilde{u}}^{n+\theta}=\theta\, \underline
{\widetilde{u}}^{n+1}+(1-\theta)\, \underline {u}^{n}\\ 
\end{equation}
Avec :
\begin{equation}
\left\{
\begin{array}{ll}
\theta = 1   & \text{Pour un sch\'ema de type Euler implicite d'ordre 1.}\\
\theta = 1/2 & \text{Pour un sch\'ema de type Cranck-Nicolson d'ordre 2.}\\
\end{array}
\right.
\end{equation}

L'équation (\ref{Base_Preduv_eq_di1}) est alors réécrite sous la forme :

\begin{equation}\label{Base_Preduv_eq_di2}
\begin{array}{c}
\displaystyle \underbrace{\left(\frac{\rho}{\Delta t} -\theta \,\dive (\rho\,\underline {u}) +\theta \,\, \Gamma^n +
\theta \,\, \rho\,\tens{K}_{\,d}^n-\theta \,T_s^{\,imp} \right)}_{\displaystyle f_s^{imp}}\, (\underline {\,\widetilde{u}}^{n+1} -\underline {u}^{n})
\\
 +\, \theta\, \dive(\underline {\widetilde{u}}^{n+1} \otimes (\rho\underline{u}))-\, \theta\,\dive (\mu_{\,tot}\,\ggrad \underline {\widetilde{u}}^{n+1}) =
\\
-\,(1-\theta)\, \dive(\underline {u}^{n} \otimes (\rho\underline{u})) +\,(1-\theta)\,\dive (\mu_{\,tot}\,\ggrad \underline {u}^{n})
\\
f_s^{exp}\left\{
\begin{array}{c}
\displaystyle 
- \grad p^{n-1+\theta} + \dive (\rho\,\underline {u})\,\underline{u}^{n} +\,(\,\Gamma^{n}\,\underline{u}_{\,i}\,)^{n+\theta_S}- \Gamma^n\,\,\underline{u}^{n}
\\
\displaystyle
-(\,\rho\,\tens{K}_{\,e}\ \underline{u}\,)^{n+\theta_S} -\rho\,\tens{K}_{\,d}^n\ \underline{u}^{n}+ (\underline{T}_{s}^{\,exp})^{\,n+\theta_S} + T_s^{\,imp}\,\,\underline {u}^{n} 
\\
\displaystyle
+[\dive (\mu_{\,tot}\,^t\ggrad \underline {u}\,)]^{n+\theta_S}-\frac {2} {3}[\,\grad (\mu_{\,tot}\,\dive \underline {u}\,)]^{n+\theta_S} + (\rho -\rho_0) \underline {g}-(\underline{turb})^{n+\theta_S}
\end{array}
\right.
\end{array}
\end{equation}

d'où l'équation résolue par le sous-programme \fort{codits} :
\begin{equation}\begin{array}{c}
\displaystyle
f_s^{\,imp}(\underline {\widetilde{u}}^{n+1}-\underline {u}^{n}) + \theta\, \dive(\underline{\widetilde{u}}^{n+1} \otimes (\rho
\underline{u})) - \theta\,\dive (\,\mu_{\,tot}\,\ggrad \underline{\widetilde{u}}^{n+1}) = 
\\\\
\displaystyle
-(1-\theta)\,\dive(\underline{u}^{n} \otimes (\rho \underline{u}))+(1-\theta)\,\dive (\,\mu_{\,tot}\,\ggrad \underline{u}^{n})
+ \underline{f}_{\,s}^{\,exp}
\end{array}
\end{equation}
La m\'ethode de discr\'etisation spatiale est d\'evelopp\'ee dans le sous-programme \fort{codits}.\\



\minititre{Remarques :}
{\tiny$\blacksquare$} Dans le cas standard sans extrapolation, le terme
$-\,T_s^{\,imp}$ n'est ajouté à $f_s^{\,imp}$ que s'il est positif afin de ne
pas affaiblir la dominance de la diagonale de la matrice à inverser.\\ 
{\tiny$\blacksquare$} Si une extrapolation est utilisée, par contre,
$\,T_s^{\,imp}$ est ajouté à $f_s^{\,imp}$ quel que soit son signe. En effet, l'idée intuitive qui
consiste à prendre~: 
\begin{equation}
\begin{cases}
\displaystyle
(\underline{T}_{s}^{\,exp} + T_{s}^{\,imp}\,\underline {u})^{\,n+\theta_S} &
\text{si } T_{s}^{\,imp} > 0\\ 
\displaystyle
(\underline{T}_{s}^{\,exp})^{\,n+\theta_S} + T_{s}^{\,imp}\,\underline{u}^{n+\theta} &\text{sinon}\\
\end{cases}
\end{equation} 
aboutit à une incohérence dans le traitement si $T_s^{imp}$ change de signe
entre deux pas de temps.\\ 
{\tiny$\blacksquare$} la partie diagonale $\tens{K}_{\,d}$ du terme
de perte de charge est utilisée dans $f_s^{\,imp}$. En fait, pour \^etre rigoureux,
il faudrait ne retenir que les contributions positives (point signal\'e dans le
sous-programme utilisateur associ\'e \fort{uskpdc}). Cette prise en compte sera \`a am\'eliorer.\\
{\tiny$\blacksquare$} Le terme $\theta\,\Gamma^{n}-\theta\,\dive
(\rho\,\underline {u})$ ne pose pas de problème pour la 
dominance de la diagonale de la matrice car il est exactement compensé par le
terme de convection (cf. \fort{covofi}). 

