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

% This file is part of Code_Saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2018 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{cfqdmv}
%
\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Fonction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Pour les notations et l'algorithme dans son ensemble,
on se reportera \`a \fort{cfbase}.

Dans le premier pas fractionnaire (\fort{cfmsvl}), on a r\'esolu une
\'equation sur la masse volumique, obtenu une pr\'ediction de la pression
et un flux convectif "acoustique".
On consid\`ere ici un second pas fractionnaire au cours duquel seul varie
le vecteur flux de masse $\vect{Q}=\rho\vect{u}$
(seule varie la vitesse au centre des cellules).
On r\'esout l'\'equation de Navier-Stokes ind\'ependamment
pour chaque direction d'espace, et l'on utilise le flux de masse acoustique
calcul\'e pr\'ec\'edemment comme flux convecteur (on pourrait aussi utiliser
le vecteur quantit\'e de mouvement du pas de temps pr\'ec\'edent).
De plus, on r\'esout en variable $\vect{u}$ et non $\vect{Q}$.

Le syst\`eme \`a r\'esoudre entre $t^*$ et $t^{**}$ est (on exclut
la turbulence, dont le traitement n'a rien de particulier dans le
module compressible)~:

\begin{equation}\label{Cfbl_Cfqdmv_eq_qdm_cfqdmv}
\left\{\begin{array}{l}

\rho^{**}=\rho^{*}=\rho^{n+1}\\
\\
\displaystyle\frac{\partial \rho\vect{u}}{\partial t}+
\divv(\vect{u} \otimes \vect{Q}_{ac}) + \gradv{P}
= \rho \vect{f}_v + \divv(\tens{\Sigma}^v)\\
\\
e^{**}=e^{*}=e^n\\

\end{array}\right.
\end{equation}

La r\'esolution de cette \'etape est similaire \`a l'\'etape
de pr\'ediction des vitesses du sch\'ema de base de \CS.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discr\'etisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------
\subsection*{Discr\'etisation en temps}
%---------------------------------

On implicite le terme de convection, ventuellement
le gradient de pression (suivant la valeur de \var{IGRDPP}, en utilisant la
pression pr\'edite lors de l'\'etape acoustique) et le terme en gradient
du tenseur des contraintes visqueuses.
On explicite les autres termes du tenseur des contraintes visqueuses.
On implicite les forces
volumiques en utilisant $\rho^{n+1}$.

On obtient alors l'\'equation discr\`ete suivante~:
\begin{equation}\label{Cfbl_Cfqdmv_eq_vitesse_cfqdmv}
\begin{array}{l}
\displaystyle\frac{(\rho\vect{u})^{n+1}-(\rho\vect{u})^n}{\Delta t^n}
+ \divv(\vect{u}^{n+1} \otimes \vect{Q}_{ac}^{n+1})
- \divv(\mu^n \gradt{\vect{u}^{n+1}})\\
\\
\text{\ \ \ \ }= \rho^{n+1} \vect{f}_v - \gradv{\widetilde{P}}
+ \divv\left(\mu^n\ ^t\gradt{\vect{u}^n}
+ (\kappa^n-\frac{2}{3}\mu^n)\divs{\vect{u}^n}\ \tens{Id}\right)\\
\end{array}
\end{equation}
avec $\widetilde{P}=P^n\text{ ou }P^{Pred}$ suivant la valeur de \var{IGRDPP}
($P^n$ par dfaut).

En pratique, dans \CS, on r\'esout cette \'equation en faisant appara\^itre \`a
gauche l'\'ecart $\vect{u}^{n+1} - \vect{u}^n$. Pour cela, on \'ecrit la
d\'eriv\'ee en temps discr\`ete sous la forme suivante~:

\begin{equation}
\begin{array}{ll}
\displaystyle
\frac{(\rho \vect{u})^{n+1} - (\rho \vect{u})^n}{\Delta t^n}
& =
\displaystyle
\frac{\rho^{n+1}\, \vect{u}^{n+1} - \rho^n\, \vect{u}^n}{\Delta t^n}\\
& =
\displaystyle
\frac{\rho^{n}\, \vect{u}^{n+1} - \rho^n\, \vect{u}^n}{\Delta t^n}+
\frac{\rho^{n+1}\, \vect{u}^{n+1} - \rho^n\, \vect{u}^{n+1}}{\Delta t^n}\\
& =
\displaystyle
\frac{\rho^{n}}{\Delta t^n}\left(\vect{u}^{n+1} - \vect{u}^n\right)+
\vect{u}^{n+1}\frac{\rho^{n+1} - \rho^n}{\Delta t^n}
\end{array}
\end{equation}

et l'on utilise alors l'\'equation de la masse discr\`ete pour \'ecrire~:
\begin{equation}
\displaystyle
\frac{(\rho \vect{u})^{n+1} - (\rho \vect{u})^n}{\Delta t^n}
=
\frac{\rho^{n}}{\Delta t^n}\left(\vect{u}^{n+1} - \vect{u}^n\right)-
\vect{u}^{n+1}\dive\,\vect{Q}_{ac}^{n+1}
\end{equation}



%---------------------------------
\subsection*{Discr\'etisation en espace}
%---------------------------------

%.................................
\subsubsection*{Introduction}
%.................................

On int\`egre l'\'equation (\ref{Cfbl_Cfqdmv_eq_vitesse_cfqdmv})
sur la cellule $i$ de volume $\Omega_i$ et
on obtient l'\'equation discr\'etis\'ee en espace~:

\begin{equation}\label{Cfbl_Cfqdmv_eq_vitesse_discrete_cfqdmv}
\begin{array}{l}
\displaystyle\frac{\Omega_i}{\Delta t^n}
(\rho_i^{n+1}\vect{u}_i^{n+1}-\rho_i^n\vect{u}_i^n)
+ \displaystyle\sum\limits_{j\in V(i)}
(\vect{u}^{n+1} \otimes \vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij}
- \displaystyle\sum\limits_{j\in V(i)}
\left(\mu^n\gradt{\vect{u}^{n+1}}\right)_{ij} \cdot \vect{S}_{ij}\\
\\
= \Omega_i\rho_i^{n+1} {\vect{f}_v}_i
- \Omega_i(\gradv{\widetilde{P}})_i
+ \displaystyle\sum\limits_{j\in V(i)}
\left(\mu^n\ ^t\gradt{\vect{u}^n} + (\kappa^n-\frac{2}{3}\mu^n)
\divs{\vect{u}^n}\ \tens{Id}\right)_{ij}\vect{S}_{ij}\\
\end{array}
\end{equation}

%.................................
\subsubsection*{Discr\'etisation de la partie ``convective''}
%.................................

La valeur \`a la face s'\'ecrit~:
\begin{equation}
(\vect{u}^{n+1} \otimes \vect{Q}_{ac}^{n+1})_{ij}\cdot \vect{S}_{ij}
= \vect{u}_{ij}^{n+1}(\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij}
\end{equation}
avec un d\'ecentrement sur la valeur de $\vect{u}^{n+1}$ aux faces~:
\begin{equation}
\begin{array}{lllll}
\vect{u}_{ij}^{n+1}
& = & \vect{u}_i^{n+1}
& \text{si\ } & (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} \geqslant 0 \\
& = & \vect{u}_j^{n+1}
& \text{si\ } & (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} < 0 \\
\end{array}
\end{equation}
que l'on peut noter~:
\begin{equation}
\vect{u}_{ij}^{n+1}
 = \beta_{ij}\vect{u}_i^{n+1} + (1-\beta_{ij})\vect{u}_j^{n+1}
\end{equation}
avec
\begin{equation}
\left\{\begin{array}{lll}
\beta_{ij} = 1 & \text{si\ }
& (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} \geqslant 0 \\
\beta_{ij} = 0 & \text{si\ }
& (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} < 0 \\
\end{array}\right.
\end{equation}

%.................................
\subsubsection*{Discr\'etisation de la partie ``diffusive''}
%.................................

La valeur \`a la face s'\'ecrit~:
\begin{equation}
\left(\mu^n\gradt{\vect{u}^{n+1}}\right)_{ij}\vect{S}_{ij}
= \mu_{ij}^n
\displaystyle \left( \frac{\partial \vect{u}}{\partial n} \right)^{n+1}_{ij}
S_{ij}
\end{equation}
avec une interpolation lin\'eaire pour $\mu^n$ aux faces (en pratique avec
$\alpha_{ij}=\frac{1}{2}$)~:
\begin{equation}
\mu_{ij}^n
= \alpha_{ij}\mu_{i}^n+(1-\alpha_{ij})\mu_{j}^n
\end{equation}
et un sch\'ema centr\'e pour le gradient normal aux faces~:
\begin{equation}
\displaystyle \left( \frac{\partial \vect{u}}{\partial n} \right)^{n+1}_{ij}
= \displaystyle\frac{\vect{u}_{J'}^{n+1}-\vect{u}_{I'}^{n+1}}{\overline{I'J'}}
\end{equation}

%.................................
\subsubsection*{Discr\'etisation du gradient de pression}
%.................................

On utilise \fort{grdcel} standard. Suivant la valeur de \var{IMRGRA},
cela correspond  une reconstruction itrative ou par moindres carrs.

%.................................
\subsubsection*{Discr\'etisation du ``reste'' du tenseur des contraintes visqueuses}
%.................................

On calcule des gradients aux cellules et on utilise une
interpolation lin\'eaire aux
faces (avec, en pratique, $\alpha_{ij}=\frac{1}{2}$)~:
\begin{equation}
\begin{array}{r}
\left(\mu^n\ ^t\gradt{\vect{u}^n} + (\kappa^n-\frac{2}{3}\mu^n)
\divs{\vect{u}^n}\ \tens{Id}\right)_{ij}\cdot\vect{S}_{ij}
= \left\{\alpha_{ij} \left(\mu^n\ ^t\gradt{\vect{u}^n}
+ (\kappa^n-\frac{2}{3}\mu^n)\divs{\vect{u}^n}\ \tens{Id}\right)_i\right.\\
\\
\left.+ (1-\alpha_{ij}) \left(\mu^n\ ^t\gradt{\vect{u}^n}
+ (\kappa^n-\frac{2}{3}\mu^n)\divs{\vect{u}^n}\ \tens{Id}\right)_j
\right\} \cdot\vect{S}_{ij}\\
\end{array}
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
On r\'esout les trois directions d'espace du syst\`eme
(\ref{Cfbl_Cfqdmv_eq_vitesse_discrete_cfqdmv}) successivement et ind\'ependamment~:
\begin{equation}\label{Cfbl_Cfqdmv_eq_vitesse_finale_cfqdmv}
\left\{\begin{array}{l}
\displaystyle\frac{\Omega_i}{\Delta t^n}
(\rho_i^{n+1}{u_i}_{(\alpha)}^{n+1}-\rho_i^n{u_i}_{(\alpha)}^n)
+ \displaystyle\sum\limits_{j\in V(i)}
{u_{ij}}_{(\alpha)}^{n+1}(\vect{Q}_{ac}^{n+1})_{ij}\cdot\vect{S}_{ij}
- \displaystyle\sum\limits_{j\in V(i)}
\mu_{ij}^n\frac{{u_j}_{(\alpha)}^{n+1}-{u_i}_{(\alpha)}^{n+1}}{\overline{I'J'}}S_{ij}\\
\qquad\qquad\qquad\qquad= \Omega_i\rho_i^{n+1} {{f_v}_i}_{(\alpha)}
- \Omega_i{(\gradv{\widetilde{P}})_{i}}_{(\alpha)}\\
\qquad\qquad\qquad\qquad + \displaystyle\sum\limits_{j\in V(i)}
\left((\mu^n\ ^t\gradt{\vect{u}^n})_{ij}\cdot\vect{S}_{ij}\right)_{(\alpha)}
 + \displaystyle\sum\limits_{j\in V(i)} \left((\kappa^n-\frac{2}{3}\mu^n)
\divs{\vect{u}^n}\right)_{ij}{S_{ij}}_{(\alpha)}\\
i = 1\ldots N \qquad \text{et} \qquad (\alpha) = x,\ y,\ z\\
\end{array}\right.
\end{equation}

Chaque syst\`eme associ\'e \`a une direction est r\'esolu par une m\'ethode
d'incr\'ement et r\'esidu en utilisant une m\'ethode de Jacobi.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section*{Points \`a traiter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% propose en patch 1.2.1

%Compl\'eter le commentaire en ent\^ete de \fort{visecv} pour prendre en compte
%la viscosit\'e en volume.