%                      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}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L'étape de pr\'ediction de la vitesse est décrite dans \fort{preduv}. On adopte ici une notation d'opérateurs pour simplifier la compréhension des algorithmes implantés.
L'équation discr\`ete de quantité de mouvement s'écrit, à un instant de r\'esolution interm\'ediaire, dans chaque direction d'espace $\alpha$ ($\alpha \in \{1,2,3\}$) :

\begin{equation}
\displaystyle
\tens{M}_{\,\alpha}^n\ \tens{R}^n\ (\widetilde{\vect{V_\alpha}} - \vect{V_\alpha}^n) + 
\tens{A_\alpha}^n\ \widetilde{\vect{V_\alpha}} = 
-\ \tens{G_\alpha}\ \vect{P}^n + \vect{S_\alpha}^n + \tens{I}_{\,s,\alpha}\ \widetilde{\vect{V_\alpha}}
\label{Base_Resolp_QDM}
\end{equation}
\begin{itemize}
\item[$\star$] $\tens{M}_{\,\alpha}^n$, de  dimension $ \var{NCEL} \times \var{NCEL}$, est une matrice diagonale contenant le rapport du volume d'une cellule et du pas de temps local ($\ \displaystyle {\tens{M}_{\,\alpha}^n}(i,i) = \frac{|\Omega_i|}{{\Delta t}_{\,\alpha,I}^n} $\ ), ${\Delta t}_{\alpha,I}^n$ représente le pas de temps dans la direction d'espace $\alpha$ à la cellule $\Omega_i$ et au pas de temps $n$,
\item[$\star$] $ \tens{R}^n $, de dimension $\var{NCEL} \times \var{NCEL}$, est la matrice diagonale contenant la masse volumique (on la sépare de la matrice de masse pour l'étape de projection de $\rho \vect{u} $). Par d\'efinition, $ \tens{R}^n(i,i) = \rho^{\,n}_{\,I}$ et l'apparition du vide \'etant exclue, cette matrice est naturellement inversible,
\item[$\star$] $\widetilde{\vect{V_\alpha}}$, de dimension $\var{NCEL}$, est un tableau où est stockée la $\alpha^{\text{ème}}$ composante de la vitesse prédite $\widetilde{\vect{u}}$,
\item[$\star$] $\vect{V_\alpha}^n$, de dimension $\var{NCEL}$, est un tableau où est stockée la $\alpha^{\text{\`eme}}$ composante de la vitesse $\vect{u}^n$ à l'instant précédent $n$,
\item[$\star$] $\tens{A_\alpha}^n$ désigne l'opérateur de convection/diffusion (cet opérateur n'est pas forcément linéaire à cause des tests de pente éventuels dans le terme de convection, il peut d\'ependre de $\vect{V_\alpha}$),

\item[$\star$] $\tens{G_\alpha}$ est l'opérateur linéaire\footnote{en toute rigueur, à cause des limitations éventuelles de gradient qui peuvent être activée par l'utilisateur, cet opérateur peut ne plus être vraiment linéaire} gradient "cellules" dans la direction $\alpha$ (il est donc appliqu\'e \`a des vecteurs de dimension \var{NCEL}),
\item[$\star$] $\vect{P}^n$, de dimension $\var{NCEL}$, est un tableau où est
stockée la pression $p^n$ du pas de temps pr\'ec\'edent à chaque cellule,

\item[$\star$] $\vect{S_\alpha}^{\,n}$ est le tableau de dimension $\var{NCEL}$ qui contient tous les termes source explicites (voir \fort{preduv} pour plus de détails),

\item[$\star$] $\tens{I}_{\,s,\alpha}$ est le tenseur diagonal relatif aux termes source implicites des composantes de la vitesse.
\end{itemize}

L'\'etape de correction consiste \`a imposer :
\begin{equation}\label{Base_Resolp_Eq_Cont}
\dive(\,{\rho\,\vect{u}}\,) = \Gamma
\end{equation}
o\`u $\Gamma$ est un terme source de masse \'eventuel. 
%($\Gamma = 0$ pour un \'ecoulement incompressible et non dilatable).

Soit $\vect W$ le tableau de dimension $3 \times \var{NCEL}$ qui contient toutes les composantes de la quantité de mouvement ($\vect{V}$ d\'esigne $\vect{V}^n$, $\vect{V}^{n+1}$ ou $\widetilde{\vect{V}}$). 
$$ \vect{W} = \tens{R}^n\ \vect{V} = \left(
                    \begin{array} {c} 
                    \rho^n\ \vect{V_1} \\
                    \rho^n\ \vect{V_2} \\
                    \rho^n\ \vect{V_3}
                    \end{array}
              \right)
$$

Soit $\tens{D}$ l'opérateur de divergence.
L'équation de continuité (\ref{Base_Resolp_Eq_Cont}) s'écrit :
\begin{equation}\notag
\tens{D}\ \vect{W} = \vect{\Gamma}
\end{equation}
$\vect{\Gamma}$ contient la valeur de $\Gamma$ au centre des cellules.

D'après l'équation discr\`ete (\ref{Base_Resolp_QDM}), on écrit pour tout $\alpha \in \{1,2,3\}$ :
\begin{equation}
\displaystyle
(\tens{M}_{\,\alpha}^n + \tens{A_{\,\alpha}}^n\ {(\tens{R}^{n})}^{-1} -\ \tens{I}_{\,s,\alpha}\ {(\tens{R}^n)}^{-1})\ \tens{R}^n\ \widetilde{\vect{V_\alpha}} = -\ \tens{G_\alpha}\ \vect{P}^n + \vect{S_{\,\alpha}}^n + \tens{M_{\,\alpha}}\ \tens{R}^n\ \vect{V_\alpha}^n
\label{Base_Resolp_eqn2}
\end{equation}

On pose :  
$$ \tens{B_{\,\alpha}} =  \tens{M}_{\,\alpha}^n + \tens{A_{\,\alpha}}^n\ {(\tens{R}^n)}^{-1} - \tens{I}_{\,s,\alpha}\ {(\tens{R}^n)}^{-1}$$
$$ \tens{B} = \left(
                    \begin{array} {ccc} 
                    \tens{B_1} & 0 & 0 \\
                    0 & \tens{B_2} & 0 \\
                    0 & 0 & \tens{B_3}
                    \end{array}
              \right)
$$
$$ \tens{G} = \left(
                    \begin{array} {c} 
                    \tens{G_1} \\
                    \tens{G_2} \\
                    \tens{G_3}
                    \end{array}
              \right)
$$
$$ \vect{S}^{\,n} = \left(
                      \begin{array} {c} 
                      \vect{S_1}^n+\ \tens{M}_{\,1}\ \tens{R}^n\ \vect{V_1}^n\\
                      \vect{S_2}^n+\ \tens{M}_{\,2}\ \tens{R}^n\ \vect{V_2}^n\\
                      \vect{S_3}^n+\ \tens{M}_{\,3}\ \tens{R}^n\ \vect{V_3}^n
                      \end{array}
                \right)
$$

On peut donc écrire :
\begin{equation}
\tens{B}\ \widetilde{\vect{W}} = -\ \tens{G}\ \vect{P}^{\,n} + \vect{S}^{\,n}
\label{Base_Resolp_eqn3}
\end{equation}

La m\'ethode des pas fractionnaires se d\'ecompose en deux \'etapes :
\begin{enumerate}
\item r\'esolution de l'\'equation (\ref{Base_Resolp_QDM}) (cette \'equation nous donne l'\'equation(\ref{Base_Resolp_eqn3})), soit :
\begin{equation}
\displaystyle
\tens{M}_{\,\alpha}^n\ \tens{R}^n\ (\widetilde{\vect{V_\alpha}} - \vect{V_\alpha}^n) + 
\tens{A_{\,\alpha}}^n\ \widetilde{\vect{V_\alpha}} - \tens{I}_{\,s,\alpha}\ \widetilde{\vect{V_\alpha}}
 = 
-\ \tens{G_\alpha}\ \vect{P}^{\,n} + \vect{S_\alpha}^n
\label{Base_Resolp_QDM2}
\end{equation}

\item soustraction\footnote{On néglige la variation des termes sources explicites $\vect{S_\alpha}^n$ qui sont toujours estim\'es en $n$.} de l'\'equation de quantit\'e de mouvement prise \`a l'instant $(n+1)$ de (\ref{Base_Resolp_QDM2}) :
\begin{equation} 
\displaystyle
\tens{M}_{\,\alpha}^n\ \tens{R}^n\ (\vect{V_\alpha}^{n+1} - \widetilde{\vect{V_\alpha}}) + 
\tens{A_{\,\alpha}}^n\ (\vect{V_\alpha}^{n+1} - \widetilde{\vect{V_\alpha}}) - \tens{I}_{\,s,\alpha}\ (\vect{V_\alpha}^{n+1} - \widetilde{\vect{V_\alpha}})
 = 
-\ \tens{G_{\,\alpha}}\ (\vect{P}^{\,n+1}- \vect{P}^{\,n})
\label{Base_Resolp_frac2}
\end{equation} 
\end{enumerate}

L'\'equation (\ref{Base_Resolp_frac2}) nous donne :

\begin{equation}
\tens{B}\ (\vect{W}^{n+1}-\widetilde{\vect{W}}) = -\ \tens{G}\ (\vect{P}^{\,n+1}-\vect{P}^{\,n})
\label{Base_Resolp_eqn5}
\end{equation}

avec :
$$ \vect{W}^{\,n+1} = \tens{R}^n\ \vect{V}^{\,n+1}
$$

Or :
$$\tens{D}\ \vect{W}^{\,n+1} = \vect{\Gamma} $$

D'où, en posant $ \delta \vect{P} = \vect{P}^{\,n+1} - \vect{P}^{\,n}$,

\begin{equation}
 \tens{D}\ \tens{B}^{-1}\,\tens{G}\ \,\delta \vect{P} =\ \tens{D}\ \widetilde{W} -\ \vect{\Gamma}
\label{Base_Resolp_eqn5bis}
\end{equation}

Il nous reste donc à inverser le système (\ref{Base_Resolp_eqn5bis}) pour d\'eterminer
$\delta p$ (donc $p^{\,n+1}$) et corriger la vitesse pour obtenir
$\vect{u}^{\,n+1}$. La correction de la vitesse se fait dans \fort{navsto}, soit
en incrémentant la vitesse par le gradient de l'incrément de pression $\delta p$
calculé, soit en reconstruisant la vitesse à partir du flux de masse actualisé
par une m\'ethode de moindres carr\'es (en appelant le sous-programme \fort{recvmc}).

Le problème provient tout d'abord du calcul de $\tens{B}^{-1}$. En effet, il a été jugé trop coûteux de calculer l'inverse de $\tens{B}$. Les algorithmes "couplage faible vitesse-pression" et celui du "couplage vitesse-pression renforc\'e" correspondent à une approximation de cet opérateur.

Dans le cas de l'algorithme "couplage faible vitesse-pression", on suppose $\tens{B}_{\,\alpha}^{-1} = \tens{M}_{\,\alpha}^{-1}$ (on pourrait inclure aussi les  termes diagonaux issus de la convection, de la diffusion et des termes source implicites).

Pour le "couplage vitesse-pression renforc\'e", on inverse le système\footnote{$\vect{\Omega} = (|\Omega_1|,...,|\Omega_{\var{NCEL}}|,|\Omega_1|,...,|\Omega_{\var{NCEL}}|,|\Omega_1|,...,|\Omega_{\var{NCEL}}|)$.} $\tens{B}\ \vect{T} = \vect{\Omega} $ et on pose $\tens{B}_{\,\alpha}^{-1}~=~\text{diag}(T_\alpha) $. Cette étape a lieu dans le sous-programme \fort{preduv}.

Les écritures à l'aide des opérateurs posent un inconvénient majeur avec la
discrétisation colocalisée. En effet, l'opérateur\footnote{On insiste sur le
fait que l'op\'erateur $\tens{G}$ est l'op\'erateur "gradients cellules"
appliqu\'e \`a la pression explicite lors de la pr\'ediction de la vitesse.}
$\tens{D}\ \tens{B}^{-1}\ \tens{G}$  peut découpler les n{\oe}uds pairs et impairs sur
un maillage cartésien et r\'egulier\footnote{Si $u_i$ est la valeur d'une
variable aux centres des cellules sur un maillage cart\'esien 1D, le Laplacien
de $u$ calcul\'e par cet op\'erateur en $i$ s'\'ecrit : $\displaystyle
\frac{u_{i-2}+2u_i-u_{i+2}}{4h^2}$, o\`u $\displaystyle h$ est le pas
d'espace. D'o\`u le d\'ecouplage des cellules.}.
Pour \'eviter ce probl\`eme, on utilise l'opérateur $\tens{L}$ (issu
naturellement de la formulation volumes finis colocalis\'ee de l'op\'erateur
$\dive(\,\grad\,))$ défini en chaque cellule\footnote{On rappelle que $Vois(i)$ est l'ensemble des centres des cellules voisines de ${\Omega_i}$ et
$\gamma_b(i)$ l'ensemble des centres \'eventuels des faces de bord de ${\Omega_i}$.} $\Omega_i$ de centre $I$ par\footnote{Si $u_i$ est la valeur d'une variable aux centres des cellules sur un maillage cart\'esien 1D, le Laplacien de $u$ par ce dernier op\'erateur en $i$ s'\'ecrit : $\displaystyle \frac{u_{i-1}+2u_i-u_{i+1}}{4h^2}$, o\`u $\displaystyle h$ est le pas d'espace.}:
\begin{equation}\label{Base_Resolp_eqn7}
\begin{array}{ll}
&(\tens{L}\ \delta \vect{P})_I = \sum\limits_{j\in Vois(i)}[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}
+ \sum\limits_{k\in {\gamma_b(i)}}[\ \tens{T}^{\,n}_{\,b_{\,ik}}\ (\grad{\delta p})_{\,f_{\,b_{\,ik}}}]\,.\,\vect{S}_{\,b_{\,ik}}
\end{array} 
\end{equation}
$\tens{T}^n$ un tenseur diagonal d'ordre 2 contenant les pas de temps dans les trois directions d'espace et $\vect{S}_{\,ij}$ (respectivement $\vect{S}_{\,b_{\,ik}}$) le vecteur surface de la face purement interne $ij$ (resp. de la face de bord $ik$). Le gradient $(\grad{ \delta p})_{\,f_{  }}$  présent dans l'équation (\ref{Base_Resolp_eqn7}) est un gradient facette\footnote{Sur maillage orthogonal, $\displaystyle (\grad{p})_{\,f_{\,ij}} \ . \ \vect{S}_{\,ij} = \frac{ p_J - p_I }{\overline{IJ}} S_{\,ij}$. $I$ (resp. $J$) et $I'$ (resp. $J'$) sont en effet superpos\'es.}.

D'un point de vue continu, on peut écrire\footnote{Dans le cas de l'algorithme "couplage faible vitesse-pression", $\widetilde{T}^n_I = \Delta t^{\,n}_I$.}:
$$(\tens{L}\ \vect{P}^{n+1})_I = \int_{\Omega_i}\dive{(\ {\tens{T}}^n \ \grad{p^{n+1}})} \ d\Omega $$
%$$(\tens{L}\ \vect{P}^{n+1})_I = \int_{\Omega_i}\dive{(\ \widetilde{\tens{T}}^n \ \grad%{p^{n+1}})} \ d\Omega $$
L'opérateur $\tens{L}$ ne pose pas de problème de découplage pair/impair sur les maillages cart\'esiens et r\'eguliers.

Avec l'algorithme "couplage faible vitesse-pression" :
$ \tens{T}_{\,I}^{\,n} = \left(
                      \begin{array} {ccc} 
                      \Delta t_I^{\,n} & 0 & 0 \\
                      0 & \Delta t_I^{\,n} & 0 \\
                      0 & 0 & \Delta t_I^{\,n}
                      \end{array}
                \right) $\\
et avec l'algorithme "couplage vitesse-pression renforc\'e" :
$ \tens{T}_{\,I}^{\,n} = \left(
                      \begin{array} {ccc} 
                      T_{\,11,\,I}^{\,n} & 0 & 0 \\
                      0 &  T_{\,22,\,I}^{\,n}&0 \\
                      0 & 0 & T_{\,33,\,I}^{\,n}
                      \end{array}
                \right) $

La matrice du syst\`eme de pression, dans le cas de l'utilisation du "couplage vitesse-pression renforc\'e", n'est pas facilement calculable \`a cause du terme : 
$$[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta
p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}$$ pour une face interne $ij$ et du terme : $$[\ \tens{T}^{\,n}_{\,b_{\,ik}}\ (\grad{\delta p})_{\,f_{\,b_{\,ik}}}]\,.
\,\vect{S}_{\,b_{\,ik}}$$ pour une face de bord.\\
En effet, \`a cause du changement de la viscosit\'e suivant la direction de l'espace du fait de l'anisotropie de $\tens{T}^n$, il faudrait calculer les gradients directionnels en fonction du gradient normal.
%(\delta p) = \tens{T}^{\,n}\ (\grad{(\delta p)})_{\,f } = 
             
On pose pour tout scalaire $a$, sans pr\'eciser pour le moment la nature du gradient $\grad$ :
$$ \vect{\widetilde{G}}\ a = \tens{T}^{\,n}\ \grad{a} = 
              \left(
                    \begin{array} {c} 
                    \displaystyle T_{11}^{\,n}\ \frac{\partial a}{\partial x}\\
                    \displaystyle T_{22}^{\,n}\ \frac{\partial a}{\partial y}\\
                    \displaystyle T_{33}^{\,n}\ \frac{\partial a}{\partial z}
                    \end{array}
              \right)
$$
On peut \'egalement d\'efinir :\\
$$ \left[(\vect{\widetilde{G}}\ a )_{\,cell}\right]_{\,ij}\,.\,\vect{S}_{\,ij} \overset{\text{\it\small d\'ef}}{=} \left[\tens{T}^{\,n}\ (\grad{a})\right]_{\,ij}\,.\,\vect{S}_{\,ij} $$  o\`u 
$\grad{a}$ est le gradient cellule classique de $a$.\\

De m\^eme , on pose :
$$ \left[(\vect{\widetilde{G}}\ a )_{\,f}\right]_{\,ij}\,.\,\vect{S}_{\,ij} \overset{\text{\it\small d\'ef}}{=} \left[\tens{T}^{\,n}\ ((\grad{a}))_{\,f}\right]_{\,ij}\,.\,\vect{S}_{\,ij} $$
o\`u $(\grad{a})_{\,f}$ est le gradient facette de $a$.\\ 
  
On doit calculer $\vect{\widetilde{G}}\ (\delta p)\,.\,\vect{S}$ \`a la face.\\
S'il s'agit d'un gradient cellule, cela n'engendre aucune difficult\'e, les composantes du gradient cellule \'etant  toutes parfaitement calculables.\\
Par contre, s'il s'agit d'un gradient facette, seule la d\'ecomposition sur la normale \`a la face est exploitable\footnote{D'apr\`es la formule $ \displaystyle\frac { p_{\,J'} - p_{\,I'}}{\overline{I'J'}} S_{\,ij}$.}. On a  besoin, explicitement et s\'epar\'ement, des quantit\'es $\displaystyle \frac{\partial (\delta p)}{\partial x}$, $\displaystyle \frac{\partial (\delta p)}{\partial y}$ et $\displaystyle \frac{\partial (\delta p)}{\partial z}$ ce qui rend difficile l'utilisation du gradient normal de l'incr\'ement de la pression.\\
On fait donc une approximation sur le gradient de l'incr\'ement de la pression en supposant qu'il est \'egal \`a sa composante normale\footnote{{\it i.e.} on ne prend pas en compte la composante tangentielle $((\grad{(\delta p)})_{\,f}\,.\,\vect{\tau}).\vect{\tau}$, $\vect{\tau}$ \'etant un vecteur tangentiel unitaire.}, c'est \`a dire :
\begin{equation}
(\grad{(\delta p)})_{\,f} \approx ((\grad{(\delta p)})_{\,f}\,.\,\vect{n})\ \vect{n}
\end{equation}
o\`u $\vect{n}$ est la normale ext\'erieure unitaire.
On a donc :
$$ \vect{\widetilde{G}}\ (\delta p) \approx ((\grad{(\delta p)})_{\,f}\,.\,\vect{n})\ (\tens{T}^n\ \vect{n})$$
On peut donc se ramener au cas d'un pas de temps scalaire $\widetilde{T}^n_{\,ij}$ en posant :
\begin{equation}
\widetilde{T}^n_{\,ij} = (\tens{T}^n_{\,ij}\ \vect{n})\,.\,\vect{n}
\label{Base_Resolp_approx2}
\end{equation}
En effet :
\begin{equation}
\begin{array} {lll}
[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}& \approx (\tens{T}^n_{\,ij}\ \vect{n})\,.\,\vect{S}_{\,ij}\ (\grad{(\delta p)}_{f_{\,ij}}\,.\,\vect{n}) \\
 & =  (\tens{T}^n_{\,ij}\ \vect{n})\,.\,\vect{n} \ ({\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{n}) \ S_{\,ij} \\
 & =  \widetilde{T}^n_{\,ij} \ ({\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{n}) \ S_{\,ij}
 \end{array}
\label{Base_Resolp_approximation}
\end{equation}
or :
\begin{equation}
\begin{array} {lll}\label{Base_Resolp_Eq_Exacte}
\widetilde{T}^n_{\,ij} \ ({\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{n}) \ S_{\,ij} &= \widetilde{T}^n_{\,ij} \ {\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{S}_{\,ij}\\
&=\,\widetilde{T}^n_{\,ij}\ \left[ P_J - P_I
 + (\ \vect{JJ'} - \vect{II'}\ )\,\displaystyle \frac{1}{2}\ (\,{\grad P}_I + \,{\grad P}_J\,)\right]\, \displaystyle \frac{S_{\,ij}}{\overline{I'J'}}
\end{array}
\end{equation}
Ceci intervient lors de la reconstruction des gradients au second membre du syst\`eme final \`a r\'esoudre, par appel aux sous-programmes \fort{itrmas} et \fort{itrgrp}.\\
En fait, l'approximation ($\displaystyle \widetilde{T}^n_{\,ij}=(\tens{T}^n_{\,ij}\ \vect{n})\,.\,\vect{n}$) n'y est pas utilis\'ee. En effet, puisqu'on calcule un gradient directionnel de l'incr\'ement de pression avec le sous-programme \fort{grdcel}, on se permet d'utiliser le tenseur $\tens{T}$ pour corriger le terme $[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}$ .\\
En pratique, on le discr\'etise, pour une face interne $ij$, par\footnote{Le coefficient $\displaystyle \frac{1}{2}$ dans (\ref{Base_Resolp_grad1}) est mis pour des raisons de stabilit\'e. En effet, l'utilisation des coefficients de pond\'eration aux faces serait plus exact mais probablement moins stable.}
 :
\begin{equation}\label{Base_Resolp_grad1}
\begin{array} {lll}
&[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}
&=\left[\widetilde{T}^n_{\,ij}\ \displaystyle \frac{P_J - P_I}{\overline{I'J'}} + \displaystyle
 \frac{\vect{JJ'} -\vect{II'}}{\overline{I'J'}} \, \frac{1}{2}\ (\,\tens{T}^{\,n}_{\,I}\ {\grad P}_I + \tens{T}^{\,n}_{\,J}\ {\grad P}_J\,)\right]\,S_{\,ij}
\end{array}
\end{equation}
au lieu de (\ref{Base_Resolp_Eq_Exacte}).
On proc\`ede de m\^eme pour les termes de bord.\\
Par la suite, on notera quand m\^eme cette expression $\widetilde{T}^n_{\,ij} \ ({\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{S}_{\,ij})$.\\ 
Il reste le terme $\tens{D}\ \widetilde{\vect{W}}$, qui \textit{via} les gradients cellule de la pression présents dans l'équation de quantité de mouvement, découple les cellules paires et impaires sur une grille cartésienne.
On utilise pour éviter ce problème un filtre de type Rhie \& Chow qui permet de lisser la contribution de la pression de l'équation de quantité de mouvement.
On écrit :
\begin{equation}\label{Base_Resolp_eqn8}
\begin{array}{lll}
(\tens{D}\ \widetilde{\vect{W}})_I = &\sum\limits_{j\in Vois(i)}[\rho^{\,n} \widetilde{\vect{u}} + \alpha_{\,Arak}\ (\vect{\widetilde{G}}\ (p^n))_{\,cell}]_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij} - \alpha_{\,Arak}\sum\limits_{j\in Vois(i)}\ \widetilde{T}^n_{\,ij}\ (\grad{p^{n}})_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij} \\
&+ \sum\limits_{k\in {\gamma_b(i)}}[\rho^{\,n} \widetilde{\vect{u}} + \alpha_{\,Arak}\ (\vect{\widetilde{G}}\ (p^n))_{\,cell}]_{\,f_{\,b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}} - \alpha_{\,Arak}\sum\limits_{k\in {\gamma_b(i)}}\widetilde{T}^n_{\,b_{\,ik}}\ (\grad{p^{n}})_{\,f_{b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}}\\   
&= \sum\limits_{j\in Vois(i)} \widetilde{m}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{m}_{\,b_{\,ik}}
\end{array}
\end{equation}
$\widetilde{m}_{\,ij}$ (resp. $\widetilde{m}_{\,b_{\,ik}}$) est le flux de masse à la face purement interne $ij$ (resp.\`a la face de bord $ik$) modifié par le filtre Rhie \& Chow.\\
$(\grad{p^{n}})_{f_{\,ij}}$ est un gradient facette alors que $[\rho^{\,n} \widetilde{\vect{u}} + \alpha_{\,Arak}\ (\vect{\widetilde{G}}\ (p^n))_{\,cell}]_{\,f_{\,ij}}$ est une valeur \`a la face interpol\'ee \`a partir des valeurs estim\'ees aux cellules (le gradient de pression dans ce terme est un gradient cellule). Cette pr\'ecision s'applique \'egalement aux termes de bord.\\
Dans \CS, $ \alpha_{\,Arak} $ est appelé pour des raisons historiques "le coefficient d'Arakawa". Il est noté \var{ARAK}.
 
Il faut remarquer, encore une fois, que pour le filtre Rhie \& Chow, le tenseur $\tens{T}^n$ est utilis\'e dans le terme volumique $\widetilde{\vect{G}}$ alors que l'approximation (\ref{Base_Resolp_approx2}) est appliqu\'ee lors du calcul du gradient facette.

Finalement, d'apr\`es (\ref{Base_Resolp_eqn5bis}) et (\ref{Base_Resolp_eqn7}), on résout :
\begin{equation}
\sum\limits_{j\in Vois(i)}\ \widetilde{T}^n_{\,ij}\ (\grad{(\delta p)})_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{T}^n_{\,b_{\,ik}}\ (\grad{(\delta p)})_{\,f_{b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}} = \sum\limits_{j\in Vois(i)} \widetilde{m}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{m}_{\,b_{\,ik}} -\ \Gamma_I
\label{Base_Resolp_Poisson}
\end{equation}
Pour la prise en compte des non orthogonalités, on résout ce système linéaire par incrément, en ajoutant au second membre les termes de reconstruction. Si $\delta(\delta p)$ est l'incrément de l'incrément (l'incrément de la variable $\delta p$ à calculer) et $k$ l'indice d'itération sur le point fixe, on résout exactement :
\begin{equation}
\begin{array}{lcl}
&\sum\limits_{j\in Vois(i)}\ \widetilde{T}^n_{\,ij}\ \displaystyle \frac{(\delta(\delta p))_I^{\,k+1}-(\delta(\delta p))_J^{\,k+1}}{\overline{I'J'}}\ S_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{T}^n_{\,b_{\,ik}}\ \displaystyle\frac{(\delta(\delta p))_I^{\,k+1}-(\delta(\delta p))_{\,b_{ik}}^{\,k+1}}{\overline{I'F}}\ S_{\,b_{\,ik}} \\
&= \sum\limits_{j\in Vois(i)} \widetilde{m}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{m}_{\,b_{\,ik}}\\
& - \sum\limits_{j\in Vois(i)}\,\widetilde{T}^n_{\,ij}\ (\grad{(\delta{p})^k})_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij} - \sum\limits_{k\in {\gamma_b(i)}}\,\widetilde{T}^n_{\,b_{\,ik}}\ (\grad{(\delta{p})^k})_{\,f_{b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}} 
- \ \Gamma_I\\
&= \sum\limits_{j\in Vois(i)} m^{\,k}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}m^{\,k}_{\,b_{\,ik}} - \ \Gamma_I
\end{array}
\label{Base_Resolp_equation resolue}
\end{equation}
avec :
\begin{equation}\label{Base_Resolp_increment}
\left\{\begin{array}{lll}
(\delta(\delta p))^{0} &= 0\\ 
(\delta(\delta p))^{k+1} &= (\delta p)^{k+1}-(\delta p)^k &\text {$\forall k \in \grandN$}\\
\end{array}\right.
\end{equation}
Ce sont les gradients facette $(\grad{(\delta{p})^k})_{\,f_{\,ij}}$ et $(\grad{(\delta{p})^k})_{\,f_{b_{\,ik}}}$ qui seront reconstruits. 





