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

%=================================
\subsection{Introduction}
%=================================

%---------------------------------
\subsubsection{Objectif}
%---------------------------------

On r\'esume ici les diff\'erentes conditions aux limites utilis\'ees pour l'algorithme 
compressible afin de fournir une vue d'ensemble. Pour atteindre cet objectif, 
il est n\'ecessaire de faire r\'ef\'erence \`a  des \'el\'ements relatifs 
\`a la discr\'etisation et au mode d'implantation des conditions aux limites. 

Lors de l'implantation, on a cherch\'e \`a pr\'eserver la coh\'erence avec l'approche 
utilis\'ee dans le cadre standard de l'algorithme incompressible de \CS.  
Il est donc conseill\'e d'avoir pris connaissance du mode de traitement des conditions 
aux limites incompressibles avant d'aborder les d\'etails de l'algorithme compressible. 

Comme pour l'algorithme incompressible, les conditions aux limites sont impos\'ees 
par le biais d'une valeur de bord associ\'ee \`a chaque variable. De plus, 
pour certaines 
fronti\`eres (parois \`a temp\'erature impos\'ee ou \`a flux thermique impos\'e), 
on dispose de deux valeurs de bord pour la m\^eme variable, l'une d'elles \'etant d\'edi\'ee au calcul du 
flux diffusif. 
Enfin, sur certains types d'entr\'ee et de sortie, on d\'efinit \'egalement  
une valeur du flux convectif au bord.  

Comme pour l'algorithme incompressible, l'utilisateur peut d\'efinir, 
pour chaque face 
de bord, des conditions aux limites pour chaque variable, mais on conseille cependant 
d'utiliser uniquement les types pr\'ed\'efinis 
d\'ecrits ci-apr\`es (entr\'ee, 
sortie, paroi, sym\'etrie) qui ont l'avantage d'assurer la coh\'erence entre 
les diff\'erentes variables et les diff\'erentes \'etapes de calcul.  


%---------------------------------
\subsubsection{Parois}
%---------------------------------

{\bf Pression} : on doit disposer d'une condition pour le calcul du gradient 
qui intervient dans l'\'etape de quantit\'e de mouvement. 
On dispose de deux types de condition, au choix de l'utilisateur~: 
\begin{itemize}
\item par d\'efaut, la pression impos\'ee au bord est proportionnelle 
\`a la valeur interne (la pression au bord est obtenue comme solution 
d'un probl\`eme de Riemann sur les \'equations d'Euler  
avec un \'etat miroir~; on distingue les cas de choc et de 
d\'etente et, dans le cas d'une d\'etente trop forte, une condition de 
Dirichlet homog\`ene est utilis\'ee pour \'eviter de voir appara\^itre une 
pression n\'egative), 
\item si l'utilisateur le souhaite (\var{ICFGRP(IPHAS)=1}), le gradient de 
pression est impos\'e \`a partir du profil de pression hydrostatique. 
\end{itemize}
\bigskip

{\bf Vitesse et turbulence}~: traitement standard (voir la documentation des
sous-programmes \fort{condli} 
et \fort{clptur}). 

{\bf Scalaires passifs}~: traitement standard    
(flux nul par d\'efaut impos\'e dans \fort{typecl}). 

{\bf Masse volumique}~: traitement standard des scalaires 
(flux nul par d\'efaut impos\'e dans \fort{typecl}). 

{\bf \'Energie et temp\'erature\footnote{Le gradient de temp\'erature est 
{\it a priori} inutile, mais peut \^etre requis par l'utilisateur.}}~: 
traitement standard des scalaires 
(flux nul par d\'efaut impos\'e dans \fort{typecl}), hormis pour le 
calcul du flux diffusif dans le cas de parois \`a temp\'erature impos\'ee 
ou \`a flux thermique impos\'e. 

{\bf Flux diffusif pour l'\'energie en paroi}~:  
l'utilisateur peut choisir (dans \fort{uscfcl}) entre une temp\'erature de paroi impos\'ee
et un flux thermique diffusif (ou "conductif") impos\'e. 
S'il ne pr\'ecise rien, on consid\`ere que la paroi est adiabatique 
(flux thermique diffusif impos\'e et de valeur nulle).  
Dans tous les cas, il faut donc disposer d'un moyen d'imposer le flux diffusif
souhait\'e. Pour cela, on d\'etermine une valeur de bord pour l'\'energie
qui, introduite dans la formule donnant le flux discret, permettra 
d'obtenir la contribution attendue 
(voir le paragraphe~\ref{Cfbl_Cfxtcl_section_cl_flux_diffusif_energie_cfener}). 
Conform\'ement \`a l'approche classique de \CS, cette valeur est 
stock\'ee sous la forme d'un couple de coefficients 
(de type \var{COEFAF}, \var{COEFBF}). 
Il est important de souligner que cette valeur de bord 
ne doit \^etre utilis\'ee que pour le calcul 
du flux diffusif~: dans les autres situations pour lesquelles 
une valeur de bord de l'\'energie ou de la temp\'erature est requise 
(calcul de gradient par exemple), on utilise une condition de flux nul 
(traitement standard des scalaires). Pour cela, on dispose d'une  
seconde valeur de bord qui est stock\'ee au moyen d'un 
couple de coefficients (\var{COEFA}, \var{COEFB}) distinct du pr\'ec\'edent.
	
{\bf Flux convectifs}~: le flux de masse dans la direction normale \`a la paroi est 
pris nul. De ce fait, les flux convectifs seront nuls quelle que soit les valeurs 
de bord impos\'ees pour les diff\'erentes variables transport\'ees.

%---------------------------------
\subsubsection{Sym\'etrie}
%--------------------------------- 

Les conditions appliqu\'ees sont les conditions classiques de l'algorithme 
incompressible (vitesse normale nulle, flux nul pour les autres variables). 

Elles sont impos\'ees dans le sous-programme \fort{typecl} essentiellement. Pour 
la pression, la condition de flux nul est impos\'ee dans \fort{cfxtcl} 
(au d\'ebut des d\'eveloppements, on appliquait le m\^eme traitement qu'en paroi, 
mais une condition de flux nul a \'et\'e pr\'ef\'er\'ee afin de s'affranchir des 
probl\`emes potentiels dans les configurations 2D).  


%---------------------------------
\subsubsection{Entr\'ees et sorties}
%--------------------------------- 

On obtient, par r\'esolution d'un probl\`eme de Riemann au bord, compl\'et\'e par 
des relations de thermodynamique (\fort{uscfth}), des valeurs de bord pour toutes 
les variables (on suppose qu'en entr\'ee, toutes les composantes de la vitesse 
sont fournies~; elles sont suppos\'ees nulles par d\'efaut, hormis pour les
entrées à $(\rho,\vect{u})$ imposés, \var{IERUCF}, pour lesquelles il faut
fournir la vitesse explicitement).

Ces valeurs de bord sont utilis\'ees de deux fa\c cons~: 
\begin{itemize}
\item elles sont utilis\'ees pour calculer les flux convectifs, en faisant 
appel au sch\'ema de Rusanov (sauf en sortie supersonique)~; ces flux sont 
directement int\'egr\'es au second membre des \'equations \`a r\'esoudre. 
\item elles servent de valeur de Dirichlet dans toutes les autres configurations 
pour lesquelles une valeur de bord est requise (calcul de flux diffusif, 
calcul de gradient...) 
\end{itemize}

Deux cas particuliers~: 
\begin{itemize}
\item aux entr\'ees ou sorties pour lesquelles toutes les variables sont impos\'ees
(\var{IESICF}), on utilise une condition de Neumann homog\`ene pour la pression 
(hormis pour le calcul du gradient intervenant dans l'\'equation de la quantit\'e 
de mouvement, qui est pris en compte par le flux convectif d\'etermin\'e par le 
sch\'ema de Rusanov). Ce choix est arbitraire (on n'a pas test\'e le comportement
de l'algorithme si l'on conserve une condition de Dirichlet sur la pression), mais 
a \'et\'e fait en supposant qu'une condition de Neumann homog\`ene serait {\it a 
priori} moins d\'estabilisante, dans la mesure o\`u, pour ce type de fronti\`ere, 
l'utilisateur peut imposer une valeur de pression tr\`es diff\'erente de 
celle r\'egnant \`a l'int\'erieur du domaine (la valeur imposée est utilisée
pour le flux convectif). 
\item pour les grandeurs turbulentes et les scalaires utilisateur, si 
le flux de masse est entrant et que l'on a fourni  
une valeur de Dirichlet (\var{RCODCL(*,*,1)} dans \fort{uscfcl}), 
on l'utilise, pour le calcul du flux convectif et du flux diffusif~;
sinon, on utilise une condition de Neumann homog\`ene (le concept
de sortie de type 9 ou 10 est couvert par cette approche). 
\end{itemize}
			      


%=================================
\subsection{Probl\`eme de Riemann au bord}
\label{Cfbl_Cfxtcl_section_pb_riemann_cfener}
%=================================

%---------------------------------
\subsubsection{Introduction}
%---------------------------------

On cherche \`a obtenir un \'etat au bord, 
pour les entr\'ees, les sorties et les parois. 

Pour cela, on fait abstraction des flux diffusifs et des sources.
Le syst\`eme r\'esultant est alors appel\'e syst\`eme
d'\'equations d'Euler. On se place de plus dans un
rep\`ere orient\'e suivant la normale au bord consid\'er\'e
$(\vect{\tau}_1,\vect{\tau}_2,\vect{n})$
et l'on ne consid\`ere que les variations suivant cette normale.
Le syst\`eme devient donc~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_euler_cfxtcl}
\begin{array}{lllll}
\displaystyle\frac{\partial \vect{W}}{\partial t}
+ \frac{\partial}{\partial n}\vect{F}_n(\vect{W})
= 0
&\text{avec}
& \vect{F}_n(\vect{W})
 = \displaystyle\sum\limits_{i=1}^3 n_i \vect{F}_i(\vect{W})
& \text{et}
& \displaystyle\frac{\partial}{\partial n}
= \displaystyle\sum\limits_{i=1}^3 n_i \frac{\partial}{\partial x_i}
\end{array}
\end{equation}

Pour d\'eterminer les valeurs des variables au bord, on recherche
l'\'evolution du probl\`eme instationnaire suivant,
appel\'e probl\`eme de Riemann~:

\unitlength=1cm
\begin{picture}(20,2.6)
\put(1.5,0){\framebox(12,2.5){}}
\put(7.5,0){\line(0,1){2.5}}
\put(7.5,2.2){bord}
\put(7.5,1.7){\vector(1,0){0.7}}
\put(8.25,1.65){$\vect{n}$}
\multiput(7.5,0)(0,0.5){4}{\line(2,3){.4}}
\put(2,1.2){$\begin{array}{c}
\text{int\'erieur}\\
\vect{W}_{\,i}\\
\text{\'etat constant dans la cellule $i$}
\end{array}$}
\put(9.5,1.2){$\begin{array}{c}
\text{ext\'erieur}\\
\vect{W}_{\,\infty}\\
\text{\'etat constant}
\end{array}$}
\end{picture}

avec $\vect{W}_{\,\infty}$ d\'ependant du type de bord et diff\'erent
de $\vect{W}_{\,i}$ {\it a priori}.

\vspace{0.3cm}

Pour r\'esoudre ce probl\`eme de Riemann, on utilisera les variables
non-conservatives $\widetilde{\vect{W}}=\ ^t(\rho, \vect{u}, P)$
et l'on retrouvera l'\'energie gr\^ace \`a l'\'equation d'\'etat.

Pour all\'eger l'\'ecriture, dans le pr\'esent 
paragraphe~\ref{Cfbl_Cfxtcl_section_pb_riemann_cfener},  
on notera aussi $\vect{W}$ le vecteur 
$^t(\rho, \vect{u}, P)$ et $\vect{u} = \vect{u}_\tau + u\,\vect{n}$
(en posant $u=\vect{u} \cdot \vect{n}$
et $\vect{u}_\tau = \vect{u} - (\vect{u} \cdot \vect{n})\vect{n}$).

La solution est une suite d'\'etats constants, dont les valeurs
d\'ependent de $\vect{W}_{\,i}$ et $\vect{W}_{\,\infty}$,
s\'epar\'es par des ondes se d\'epla\c cant \`a des vitesses donn\'ees
par les valeurs propres du syst\`eme $(\lambda_i)_{i=1\ldots 5}$.
On repr\'esente les caract\'eristiques du syst\`eme sur le sch\'ema suivant~: 

\unitlength=1cm
\begin{picture}(20,3.5)
\put(3,0){\vector(1,0){8}}
\put(7,0){\vector(0,1){3}}
\put(11,0.1){$x$}
\put(6.8,2.9){$t$}
\put(7,3.1){bord}
\put(7,2.6){\vector(1,0){0.5}}
\put(7.55,2.65){$\vect{n}$}
\multiput(7,0)(0,0.5){5}{\line(2,3){.4}}
\put(7,0){\line(-1,2){1.4}}
\put(4,2.9){$\lambda_1=u-c$}
\put(7,0){\qbezier[15](0,0)(0.7,1.4)(1.4,2.8)}
\put(8.3,2.9){$\lambda_{2,3,4}=u$}
\put(7,0){\line(2,1){3}}
\put(9.5,1.6){$\lambda_5=u+c$}
\put(5,1.5){$\vect{W}_{\,i}$}
\put(6.3,2){$\vect{W}_{\,1}$}
\put(8.4,1.6){$\vect{W}_{\,2}$}
\put(9.5,0.5){$\vect{W}_{\,\infty}$}
\end{picture}

Comme valeurs des variables au bord, on prendra les valeurs correspondant \`a
l'\'etat constant qui contient le bord ($\vect{W}_1$ dans l'exemple
pr\'ec\'edent).

Il faut remarquer que la solution du probl\`eme de Riemann d\'epend de la 
thermodynamique et devra donc \^etre calcul\'ee et cod\'ee 
par l'utilisateur si la thermodynamique n'a pas \'et\'e pr\'evue 
(en version 1.2, la seule 
thermodynamique pr\'evue est celle des gaz parfaits). 



%---------------------------------
\subsubsection{En paroi, pour la condition de pression (sans effet de gravit\'e)}
%---------------------------------

Pour les faces de paroi, on d\'efinit \`a l'ext\'erieur du domaine
un \'etat miroir $\vect{W}_{\,\infty}$ par~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_paroi_cfxtcl}
\begin{array}{lllll}
\vect{W}_{\,i} &=&
\left(\begin{array}{l}
\rho_i\\ {\vect{u}_\tau}_i\\ u_i\\ P_i
\end{array}\right)
\qquad
\vect{W}_{\,\infty} &=&
\left(\begin{array}{lll}
\rho_\infty &=& \rho_i\\
{\vect{u}_\tau}_\infty &=& {\vect{u}_\tau}_i\\
u_\infty &=& -u_i\\
P_\infty &=& P_i
\end{array}\right)
\end{array}
\end{equation}

\vspace{0.5cm}

Les solutions d\'ependent de l'orientation de la vitesse dans la cellule
de bord~:

\vspace{0.5cm}

\begin{enumerate}

\item Si $u_i \leqslant 0$,
la solution est une double d\'etente sym\'etrique.

\unitlength=1cm
\begin{picture}(20,4)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){paroi}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-1,1){2.5}}
\put(4,0.5){\line(-5,4){2.8}}
\put(4,0.5){\line(-4,5){2}}
\put(0,3.1){$\lambda_1=u-c\ (<0)$}
\put(4.2,2.7){$\lambda_{2,3,4}=0$}
\put(4,0.5){\line(1,1){2.5}}
\put(4,0.5){\line(5,4){2.8}}
\put(4,0.5){\line(4,5){2}}
\put(6.5,3.1){$\lambda_5=u+c\ (>0)$}
\put(1.5,1.5){$\vect{W}_{\,i}$}
\put(3.1,2){$\vect{W}_{\,1}$}
\put(4.5,2){$\vect{W}_{\,2}$}
\put(6,1.5){$\vect{W}_{\,\infty}$}
\put(8.5,2){$\vect{W}_{\,paroi} = \vect{W}_{\,1} = \vect{W}_{\,2}$}
\put(12,2)
{$\left\{\begin{array}{l}
\rho_p = \rho_1 = \rho_2\\
u_p = u_1 = u_2\\
P_p = P_1 = P_2
\end{array}\right.$}
\end{picture}

La conservation de la vitesse tangentielle \`a travers la 1-onde donne
${\vect{u}_\tau}_p = {\vect{u}_\tau}_i$.
Par des consid\'erations de sym\'etrie on trouve $u_p = 0$.
Puis on obtient $\rho_p$ et $P_p$ en \'ecrivant la conservation
des invariants de Riemann \`a travers la 1-d\'etente~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_invariants_detente_cfxtcl}
\begin{array}{lll}
\left\{\begin{array}{l}
u_1 + \displaystyle\int_0^{\rho_1} \frac{c}{\rho} d\rho
= u_i + \displaystyle\int_0^{\rho_i} \frac{c}{\rho} d\rho\\
\\
s_1 = s_i
\end{array}\right.
&
\Rightarrow
\left\{\begin{array}{ll}
\displaystyle\int_{\rho_i}^{\rho_1} \frac{c}{\rho} d\rho = u_i
& \Rightarrow \rho_p=\rho_1\\
\\
s(P_1,\rho_1) = s(P_i,\rho_i)
& \Rightarrow P_p=P_1
\end{array}\right.
\end{array}
\end{equation}

\vspace{0.5cm}

\item  Si $u_i > 0$,
la solution est un double choc sym\'etrique.

\unitlength=1cm
\begin{picture}(20,4)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){paroi}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-1,1){2.5}}
\put(0,3.1){$\lambda_1=u-c\ (<0)$}
\put(4.4,2.7){$\lambda_{2,3,4}=0$}
\put(4,0.5){\line(1,1){2.5}}
\put(6.5,3.1){$\lambda_5=u+c\ (>0)$}
\put(1.5,1.5){$\vect{W}_{\,i}$}
\put(3,2){$\vect{W}_{\,1}$}
\put(4.7,2){$\vect{W}_{\,2}$}
\put(6,1.5){$\vect{W}_{\,\infty}$}
\put(8.5,2){$\vect{W}_{\,paroi} = \vect{W}_{\,1} = \vect{W}_{\,2}$}
\put(12,2)
{$\left\{\begin{array}{l}
\rho_p = \rho_1 = \rho_2\\
u_p = u_1 = u_2\\
P_p = P_1 = P_2
\end{array}\right.$}
\end{picture}

De m\^eme que pr\'ec\'edemment,
on trouve ${\vect{u}_\tau}_p = {\vect{u}_\tau}_i$ et $u_p = 0$,
puis $\rho_p$ et $P_p$ en \'ecrivant les relations de saut
\`a travers le 1-choc~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_saut_choc_cfxtcl}
\begin{array}{lll}
\left\{\begin{array}{l}
\rho_1 \rho_i (u_1 - u_i)^2
= (P_1 - P_i)(\rho_1 - \rho_i)\\
\\
2\rho_1 \rho_i (\varepsilon_1 - \varepsilon_i)
= (P_1 + P_i)(\rho_1 - \rho_i)
\end{array}\right.
& \text{avec } \varepsilon = \varepsilon(P,\rho)
&
\Rightarrow
\left\{\begin{array}{l}
\rho_p=\rho_1\\
\\
P_p=P_1
\end{array}\right.
\end{array}
\end{equation}

\bigskip
Pour les gaz parfaits, 
avec $M_i = \displaystyle\frac{\vect{u}_i \cdot \vect{n}}{c_i}$
(Nombre de Mach de paroi), on a~: 

\begin{itemize}

\item Cas d\'etente ($M_i \leqslant 0$)~:\\
$$
\begin{array}{l}
\left\{\begin{array}{lll}
P_p=0 & \text{si} & 1 + \displaystyle\frac{\gamma-1}{2}M_i<0\\
P_p = P_i \left(1 + \displaystyle\frac{\gamma-1}{2}M_i\right)
^{\frac{2\gamma}{\gamma-1}} & \text{sinon}\\
\end{array}\right.\\
\\
\rho_p=\rho_i \left(\displaystyle\frac{P_p}{P_i}\right)^{\frac{1}{\gamma}}\\
\end{array}
$$

\item Cas choc ($M_i > 0$)~:\\
$$
\begin{array}{l}
P_p = P_i \left(1 + \displaystyle\frac{\gamma(\gamma+1)}{4}M_i^2
+\gamma M_i \displaystyle\sqrt{1+\displaystyle\frac{(\gamma+1)^2}{16}M_i^2}\right)
\\
\\
\rho_p=\rho_i \left(\displaystyle\frac{P_p-P_i}
{P_p-P_i-\rho_i (\vect{u}_i \cdot \vect{n})^2}\right)\\
\end{array}
$$

\end{itemize}

\end{enumerate}

En pratique, le flux convectif normal \`a la paroi est nul et seule 
la condition de pression d\'etermin\'ee ci-dessus est effectivement 
utilis\'ee (pour le calcul du gradient sans effet de gravit\'e). 

%---------------------------------
\subsubsection{En sortie}
%---------------------------------

Il existe deux cas de traitement des conditions en sortie,
selon le nombre de Mach normal \`a la face de bord
($c_i$ est la vitesse du son dans la cellule de bord)~:
$$M_i = \displaystyle\frac{u_i}{c_i}
= \displaystyle\frac{\vect{u}_i \cdot \vect{n}}{c_i}$$

\paragraph{Sortie supersonique (condition ISSPCF de 
\fort{uscfcl})~:}
$M_i \geqslant 1 \Rightarrow u_i - c_i \geqslant 0$
\nopagebreak
\linebreak
\unitlength=1cm
\begin{picture}(20,4.5)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){bord}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(1,2){1.4}}
\put(5.3,3.4){$\lambda_1=u-c\ (\geqslant 0)$}
\put(4,0.5){\qbezier[20](0,0)(1.2,1.2)(2.4,2.4)}
\put(6.5,2.9){$\lambda_{2,3,4}=u$}
\put(4,0.5){\line(2,1){3}}
\put(7.1,2){$\lambda_5=u+c$}
\put(2,2){$\vect{W}_{\,i}$}
\put(5.2,2.5){$\vect{W}_{\,1}$}
\put(6,2){$\vect{W}_{\,2}$}
\put(6.5,1){$\vect{W}_{\,\infty}$}
\put(10,2){$\vect{W}_{\,bord} = \vect{W}_{\,i}$}
\end{picture}

Toutes les caract\'eristiques sont sortantes,
on conna\^it donc toutes les conditions au bord~:

\begin{equation}
\left\{\begin{array}{l}
\rho_b = \rho_i\\
{\vect{u}_\tau}_b = {\vect{u}_\tau}_i\\
u_b = u_i\\
P_b = P_i
\end{array}\right.
\end{equation}

\paragraph{Sortie subsonique (condition ISOPCF de 
\fort{uscfcl})~:}
$0 \leqslant M_i < 1 \Rightarrow (u_i \geqslant 0 \text{ et } u_i - c_i < 0)$

\unitlength=1cm
\begin{picture}(20,4.5)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){bord}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-1,2){1.4}}
\put(1,3.4){$\lambda_1=u-c\ (<0)$}
\put(4,0.5){\qbezier[15](0,0)(0.7,1.4)(1.4,2.8)}
\put(5.3,3.4){$\lambda_{2,3,4}=u\ (\geqslant 0)$}
\put(4,0.5){\line(2,1){3}}
\put(6.5,2.1){$\lambda_5=u+c$}
\put(2,2){$\vect{W}_{\,i}$}
\put(3.3,2.5){$\vect{W}_{\,1}$}
\put(5.4,2.1){$\vect{W}_{\,2}$}
\put(6.5,1){$\vect{W}_{\,\infty}$}
\put(10,2){$\vect{W}_{\,bord} = \vect{W}_{\,1}$}
\put(12.5,2)
{$\left\{\begin{array}{l}
\rho_b = \rho_1\\
u_b = u_1\\
P_b = P_1
\end{array}\right.$}
\end{picture}

On a une caract\'eristique entrante,
on doit donc imposer une seule condition  au bord
(en g\'en\'eral la pression de sortie $P_{ext}$).

On conna\^it alors $P_b = P_{ext}$ et ${\vect{u}_\tau}_b = {\vect{u}_\tau}_i$
(par conservation de la vitesse tangentielle \`a travers la 1-onde).
Pour trouver les inconnues manquantes ($\rho_b$ et $u_b$)
on doit r\'esoudre le passage de la 1-onde~:

\begin{enumerate}

\item Si $P_{ext} \leqslant P_i$,
on a une 1-d\'etente.

On \'ecrit la conservation
des invariants de Riemann \`a travers la 1-d\'etente~:
\begin{equation}
\begin{array}{lll}
\left\{\begin{array}{l}
s_1 = s_i\\
\\
u_1 + \displaystyle\int_0^{\rho_1} \frac{c}{\rho} d\rho
= u_i + \displaystyle\int_0^{\rho_i} \frac{c}{\rho} d\rho
\end{array}\right.
&
\Rightarrow
\left\{\begin{array}{ll}
s(P_{ext},\rho_1) = s(P_i,\rho_i)
& \Rightarrow \rho_b=\rho_1\\
\\
u_1 = u_i - \displaystyle\int_{\rho_i}^{\rho_1} \frac{c}{\rho} d\rho
& \Rightarrow u_b = u_1
\end{array}\right.
\end{array}
\end{equation}

\item Si $P_{ext} > P_i$,
on a un 1-choc.

On \'ecrit les relations de saut \`a travers le 1-choc~:
\begin{equation}
\begin{array}{lll}
\left\{\begin{array}{l}
\rho_1 \rho_i (u_1 - u_i)^2
= (P_{ext} - P_i)(\rho_1 - \rho_i)\\
\\
2\rho_1 \rho_i (\varepsilon(P_{ext},\rho_1) - \varepsilon(P_i,\rho_i))
= (P_{ext} + P_i)(\rho_1 - \rho_i)
\end{array}\right.
&
\Rightarrow
\left\{\begin{array}{l}
\rho_b=\rho_1\\
\\
u_b = u_1
\end{array}\right.
\end{array}
\end{equation}

\bigskip
Pour les gaz parfaits, on a~:
\begin{itemize}

\item Cas d\'etente ($P_{ext} \leqslant P_i$)~:\\
$$
\begin{array}{l}
P_b=P_{ext}\\
\\
\rho_b=\rho_i \left(\displaystyle\frac{P_{ext}}{P_i}\right)^{\frac{1}{\gamma}}
\end{array}
$$

\item Cas choc ($P_{ext} > P_i$)~:\\
$$
\begin{array}{l}
P_b=P_{ext}\\
\\
\rho_b=\rho_i \left(\displaystyle\frac{P_{ext}-P_i}{P_{ext}-P_i-\rho_i
(\vect{u}_i \cdot \vect{n} - \vect{u}_b \cdot \vect{n})^2}\right)
= \rho_i \left(\displaystyle\frac{(\gamma+1)P_{ext}+(\gamma-1)P_i}
{(\gamma-1)P_{ext}+(\gamma+1)P_i}\right)\\
\end{array}
$$

\end{itemize}

\end{enumerate}

La valeur de la masse volumique au bord intervient en particulier  
dans le flux de masse.


%---------------------------------
\subsubsection{En entr\'ee}
%---------------------------------

L'utilisateur impose les valeurs qu'il souhaite pour les variables
en entr\'ee~:
$$
\begin{array}{lllll}
\vect{W}_{\,ext} &=&
\left(\begin{array}{l}
\rho_{ext}\\ {\vect{u}_\tau}_{ext}\\ u_{ext}\\ P_{ext}
\end{array}\right)
\end{array}
$$

De m\^eme que pr\'ec\'edemment, il existe deux cas de traitement
des conditions en entr\'ee,
pilot\'es par le nombre de Mach entrant, normalement \`a la face de bord
(avec $c_{ext}$ la vitesse du son en entr\'ee)~:
$$M_{ext} = \displaystyle\frac{u_{ext}}{c_{ext}}
= \displaystyle\frac{\vect{u}_{ext} \cdot \vect{n}}{c_{ext}}$$

\paragraph{Entr\'ee supersonique (condition IESICF de 
\fort{uscfcl})~:}
$M_{ext} \leqslant -1 \Rightarrow u_{ext} + c_{ext} \leqslant 0$

\unitlength=1cm
\begin{picture}(20,4.5)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){bord}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-2,1){3}}
\put(0,2.1){$\lambda_1=u-c$}
\put(4,0.5){\qbezier[20](0,0)(-1.2,1.2)(-2.4,2.4)}
\put(0,2.9){$\lambda_{2,3,4}=u$}
\put(4,0.5){\line(-1,2){1.4}}
\put(1,3.4){$\lambda_5=u+c\ (\leqslant 0)$}
\put(1.5,1){$\vect{W}_{\,i}$}
\put(2,1.7){$\vect{W}_{\,1}$}
\put(2.2,2.5){$\vect{W}_{\,2}$}
\put(6,2){$\vect{W}_{\,\infty}$}
\put(10,2){$\vect{W}_{\,bord} = \vect{W}_{\,\infty} = \vect{W}_{\,ext}$}
\end{picture}

Toutes les caract\'eristiques sont entrantes,
toutes les conditions au bord sont donc impos\'ees par l'utilisateur.

\begin{equation}
\left\{\begin{array}{l}
\rho_b = \rho_{ext}\\
{\vect{u}_\tau}_b = {\vect{u}_\tau}_{ext}\\
u_b = u_{ext}\\
P_b = P_{ext}
\end{array}\right.
\end{equation}


\paragraph{Entr\'ee subsonique (condition IERUCF de 
\fort{uscfcl})~: }
$$-1 < M_{ext} \leqslant 0
\Rightarrow (u_{ext} \leqslant 0 \text{ et } u_{ext} + c_{ext} > 0)$$ 


\unitlength=1cm
\begin{picture}(20,4.5)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){bord}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-2,1){3}}
\put(0,2.1){$\lambda_1=u-c$}
\put(4,0.5){\qbezier[15](0,0)(-0.7,1.4)(-1.4,2.8)}
\put(1.1,3.4){$\lambda_{2,3,4}=u\ (\leqslant 0)$}
\put(4,0.5){\line(1,2){1.4}}
\put(5.3,3.4){$\lambda_5=u+c\ (>0)$}
\put(1.5,1){$\vect{W}_{\,i}$}
\put(2,2.1){$\vect{W}_{\,1}$}
\put(3.3,2.5){$\vect{W}_{\,2}$}
\put(6,2){$\vect{W}_{\,\infty}$}
\put(10,2){$\vect{W}_{\,bord} = \vect{W}_{\,2}$}
\put(12.5,2)
{$\left\{\begin{array}{l}
\rho_b = \rho_2\\
u_b = u_2\\
P_b = P_2
\end{array}\right.$}
\end{picture}


On a une caract\'eristique sortante.
L'utilisateur doit donc laisser un degr\'e de libert\'e.

En g\'en\'eral, on impose le flux de masse en entr\'ee, donc $\rho_{ext}$
et $u_{ext}$, et l'on calcule la pression au bord en r\'esolvant
le passage des 1$\sim$4-ondes.
On conna\^it aussi ${\vect{u}_\tau}_b = {\vect{u}_\tau}_{ext}$,
par conservation de la vitesse tangentielle \`a travers la 5-onde.

\begin{enumerate}

\item Si $u_{ext} \geqslant u_i$,
on a une 1-d\'etente.

On \'ecrit la conservation
des invariants de Riemann \`a travers la 1-d\'etente
et la conservation de la vitesse et de la pression \`a travers le contact~:
\begin{equation}
\begin{array}{lll}
\begin{array}{l}
\left\{\begin{array}{l}
u_1 + \displaystyle\int_0^{\rho_1} \frac{c}{\rho} d\rho
= u_i + \displaystyle\int_0^{\rho_i} \frac{c}{\rho} d\rho\\
\\
s_1 = s_i
\end{array}\right.\\
\\
\left\{\begin{array}{l}
u_1 = u_2 = u_{ext}\\
\\
P_1 = P_2
\end{array}\right.
\end{array}
&
\Rightarrow
\left\{\begin{array}{ll}
\displaystyle\int_{\rho_i}^{\rho_1} \frac{c}{\rho} d\rho
= u_i - u_{ext}
& \Rightarrow \rho_1\\
\\
s(P_2,\rho_1) = s(P_i,\rho_i)
& \Rightarrow P_b = P_2
\end{array}\right.
\end{array}
\end{equation}

\item Si $u_{ext} < u_i$,
on a un 1-choc.

On \'ecrit les relations de saut \`a travers le 1-choc
et la conservation de la vitesse et de la pression \`a travers le contact~:
\begin{equation}
\begin{array}{lll}
\begin{array}{l}
\left\{\begin{array}{l}
\rho_1 \rho_i (u_1 - u_i)^2
= (P_1 - P_i)(\rho_1 - \rho_i)\\
\\
2\rho_1 \rho_i (\varepsilon_1 - \varepsilon_i)
= (P_1 + P_i)(\rho_1 - \rho_i)\\
\\
\varepsilon = \varepsilon(P,\rho)
\end{array}\right.\\
\\
\left\{\begin{array}{l}
u_1 = u_2 = u_{ext}\\
\\
P_1 = P_2
\end{array}\right.
\end{array}
&
\Rightarrow
\left\{\begin{array}{l}
\rho_1\\
\\
P_b = P_2
\end{array}\right.
\end{array}
\end{equation}

\bigskip
Pour les gaz parfaits, on a~: 

\begin{itemize}

\item Cas d\'etente ($\delta M \leqslant 0$)~:\\
$$
\begin{array}{l}
\left\{\begin{array}{lll}
P_b=0 & \text{si} & 1 + \displaystyle\frac{\gamma-1}{2}\delta M<0\\
P_b = P_i \left(1 + \displaystyle\frac{\gamma-1}{2}\delta M\right)
^{\frac{2\gamma}{\gamma-1}} & \text{sinon}\\
\end{array}\right.\\
\\
\rho_b= \rho_{ext}\\
\end{array}
$$

\item Cas choc ($\delta M > 0$)~:\\
$$
\begin{array}{l}
P_b = P_i \left(1 + \displaystyle\frac{\gamma(\gamma+1)}{4}\delta M^2
+\gamma \delta M \displaystyle\sqrt{1+\displaystyle\frac{(\gamma+1)^2}{16}\delta M^2}\right)
\\
\\
\rho_b=\rho_{ext}\\
\end{array}
$$

\end{itemize}

\end{enumerate}


%=================================
\subsection{Condition de pression en paroi avec effets de gravit\'e}
%=================================

Le probl\`eme de Riemann consid\'er\'e pr\'ec\'edemment
ne prend pas en compte les effets de la gravit\'e. 
Or, dans certains cas, si l'on ne prend pas en compte le gradient de 
pression ``hydrostatique'',  on peut obtenir une solution erron\'ee  
(en particulier, par exemple, 
on peut cr\'eer une source de quantit\'e de mouvement
non physique dans un milieu initialement au repos).

\'Ecrivons l'\'equilibre local dans la maille de bord~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_equilibre_local_cfxtcl}
\gradv{P} = \rho \vect{g}
\end{equation}

Pour simplifier la r\'esolution, on peut utiliser la formulation 
de (\ref{Cfbl_Cfxtcl_eq_equilibre_local_cfxtcl}) en incompressible 
(c'est cette approche qui a \'et\'e adopt\'ee dans \CS)~: 
\begin{equation}\label{Cfbl_Cfxtcl_eq_equilibre_incompressible_cfxtcl}
\begin{array}{lll}
\left(\gradv{P}\right)_i = \rho_i \vect{g}
& \text{ce qui donne}
& P_{paroi} = P_i + \rho_i \vect{g} \cdot (\vect{x}_{paroi} - \vect{x}_i)
\end{array}
\end{equation}


Une autre approche (d\'ependante de l'\'equation d'\'etat)
consiste \`a r\'esoudre l'\'equilibre local avec la formulation
compressible (\ref{Cfbl_Cfxtcl_eq_equilibre_local_cfxtcl}), en supposant de plus que
la maille est isentropique~:
\begin{equation}
\left\{\begin{array}{lll}
\gradv{P} = \rho \vect{g}\\
\\
P = P(\rho,s_i)
\end{array}\right.
\end{equation}
Ce qui donne, pour un gaz parfait~:
\begin{equation}
\label{Cfbl_Cfxtcl_eq_equilibre_compressible_cfxtcl}
 P_{paroi} = P_i \left( 1+ \displaystyle\frac{\gamma -1}{\gamma}
\displaystyle\frac{\rho_i}{P_i} \vect{g} \cdot (\vect{x}_{paroi} - \vect{x}_i)
\right)^{\frac{\gamma}{\gamma -1}}
\end{equation}

\paragraph{Remarque~:} 
la formule issue de l'incompressible (\ref{Cfbl_Cfxtcl_eq_equilibre_incompressible_cfxtcl})
est une lin\'earisation de la formule (\ref{Cfbl_Cfxtcl_eq_equilibre_compressible_cfxtcl}).
Dans les cas courants elle s'\'eloigne tr\`es peu de la formule exacte.
Dans des conditions extr\^emes, 
si l'on consid\`ere par exemple 
de l'air \`a $1000K$ et $10bar$, avec une acc\'el\'eration
de la pesanteur $g=1000m/s^2$ et une diff\'erence de hauteur entre
le centre de la cellule et le centre de la face de bord de $10m$,
l'expression (\ref{Cfbl_Cfxtcl_eq_equilibre_compressible_cfxtcl}) donne $P_{paroi} = 1034640,4Pa$
et l'expression (\ref{Cfbl_Cfxtcl_eq_equilibre_incompressible_cfxtcl}) donne $P_{paroi} = 1034644,7Pa$,
soit une diff\'erence relative de moins de $0,001\%$.
On voit aussi que la diff\'erence entre la pression calcul\'ee au centre
de la cellule et celle calcul\'ee au bord est de l'ordre de~$3\%$.

%=================================
\subsection{Sch\'ema de Rusanov pour le calcul de flux convectifs au bord}
%=================================


%---------------------------------
\subsubsection{Introduction}
%---------------------------------

Le sch\'ema de Rusanov est utilis\'e pour certains types de conditions aux
limites afin de passer du vecteur d'\'etat calcul\'e au bord comme indiqu\'e
pr\'ec\'edemment (solution du probl\`eme de Riemann) \`a un flux convectif de
bord (pour la masse, la quantit\'e de
mouvement et l'\'energie). L'utilisation de ce sch\'ema (d\'ecentr\'e amont) 
permet de gagner en  stabilit\'e.

Le sch\'ema de Rusanov est appliqu\'e aux fronti\`eres auxquelles on consid\`ere
qu'il est le plus probable de rencontrer des conditions en accord imparfait  
avec l'\'etat r\'egnant dans le domaine, conditions qui sont donc susceptibles de
d\'estabiliser le calcul~: il s'agit des entr\'ees et des sorties (fronti\`eres
de type IESICF, ISOPCF, IERUCF, IEQHCF). En sortie
supersonique (ISSPCF) cependant, le sch\'ema de Rusanov est inutile et 
n'est donc pas appliqu\'e~:  
en effet, pour ce type de fronti\`ere, l'\'etat impos\'e au bord est exactement
l'\'etat amont et le d\'ecentrement du sch\'ema de Rusanov n'apporterait donc
rien.   

%---------------------------------
\subsubsection{Principe}
%---------------------------------

Pour le calcul du flux d\'ecentr\'e de Rusanov, on consid\`ere 
le syst\`eme hyperbolique 
constitu\'e des seuls termes convectifs issus 
des \'equations de masse, quantit\'e de mouvement et \'energie. Ce
syst\`eme est \'ecrit, par changement de variable, en non conservatif 
(on utilise la relation 
$\displaystyle P=\frac{\rho\varepsilon}{\gamma-1}$ et  
on note $u_\xi$ les composantes de $\vect{u}$)~:

\begin{equation}
\left\{\begin{array}{lllll}
\displaystyle\frac{\partial\rho}{\partial t} 
&+&\rho\divv{\,\vect{u}} + \vect{u}\,\grad{\,\rho}&=& 0 \\
\displaystyle\frac{\partial u_\xi}{\partial t}
&+& \vect{u}\,\grad{u_\xi}+\displaystyle\frac{1}{\rho}\,\frac{\partial
P}{\partial \xi} &=& 0 \\
\displaystyle\frac{\partial P}{\partial t}
&+&\gamma\,P\,\dive{\vect{u}}+\vect{u}\,\grad{P}&=& 0 
\end{array}\right.
\end{equation}

En notant le vecteur d'\'etat $\vect{W}= (\rho,\vect{u},P)^t$, 
ce syst\`eme est not\'e~: 
\begin{equation}
\displaystyle\frac{\partial \vect{W}}{\partial t} +\dive{\,\vect{F}(\vect{W})} = 0
\end{equation}

Avec $\delta\,\vect{W}$ l'incr\'ement temporel du vecteur d'\'etat, $\vect{n}$ la
normale \`a une face, $ij$ la face interne partag\'ee par les cellules $i$ et
$j$ et $ik$ la face de bord $k$ associ\'ee \`a la cellule $i$,  
la discr\'etisation spatiale conduit \`a~: 
\begin{equation}
\displaystyle\frac{|\Omega_i|}{\Delta t}\delta\,\vect{W}_i 
+\sum\limits_{j\in\,Vois(i)}\int_{S_{ij}} \vect{F}(\vect{W})\,\vect{n}\,dS 
+\sum\limits_{k\in {\gamma_b(i)}}\int_{S_{\,{b}_{ik}}} \vect{F}(\vect{W})\,\vect{n}\,dS
=0
\end{equation}

Sur une face de bord donn\'ee, 
on applique le sch\'ema de Rusanov pour calculer le flux
comme suit~: 
\begin{equation}
\frac{1}{|S_{\,{b}_{ik}}|}\int_{S_{\,{b}_{ik}}} \vect{F}(\vect{W})\,\vect{n}\,dS
=\frac{1}{2}\left(\vect{F}(\vect{W}_i)+\vect{F}(\vect{W}_{\,{b}_{ik}})\right)\cdot\vect{n}_{\,{b}_{ik}}
-\frac{1}{2}\rho_{rus\,{b}_{ik}}\left(\vect{W}_{\,{b}_{ik}}-\vect{W}_i\right)=\vect{F}_{rus\,{b}_{ik}}(\vect{W})
\end{equation}

Dans cette relation, $\vect{W}_{\,{b}_{ik}}$ est  le vecteur d'\'etat 
$\vect{W}_{\infty}$, connu au bord (tel
qu'il r\'esulte de la r\'esolution du probl\`eme de Riemann au bord
pr\'esent\'ee plus haut pour chaque type de fronti\`ere consid\'er\'e). 

%---------------------------------
\subsubsection{Param\`etre de d\'ecentrement $\rho_{rus\,{b}_{ik}}$}
%---------------------------------

Pour chaque face de bord, le scalaire $\rho_{rus\,{b}_{ik}}$ est la 
plus grande valeur du rayon spectral de la matrice jacobienne
$\displaystyle\frac{\partial\,\vect{F}_n(\vect{W})}{\partial \vect{W}}$ 
obtenu pour les vecteurs d'\'etat $\vect{W}_i$ et $\vect{W}_{\,{b}_{ik}}$. 

$\vect{F}_n$ est la composante du
flux $\vect{F}$ dans la direction de la normale \`a la face de bord,
$\vect{n}_{\,{b}_{ik}}$. Utiliser $\vect{F}_n$ 
pour la d\'etermination du
param\`etre de d\'ecentrement $\rho_{rus\,{b}_{ik}}$
rel\`eve d'une approche classique qui consiste 
\`a remplacer le syst\`eme tridimensionnel
initial par le syst\`eme unidimensionnel projet\'e dans la direction 
normale \`a la face, en n\'egligeant les variations du vecteur d'\'etat 
$\vect{W}$ dans la direction tangeante \`a la face~: 
\begin{equation}
\displaystyle\frac{\partial \vect{W}}{\partial t} +\frac{\partial\,\vect{F}_n(\vect{W})}{\partial
\vect{W}}\,\frac{\partial \vect{W}}{\partial n} = 0
\end{equation}

De mani\`ere plus explicite, si l'on se place dans un rep\`ere de calcul ayant
$\vect{n}_{\,{b}_{ik}}$ comme  vecteur de base, et si l'on note $u$ la
composante de vitesse associ\'ee, le syst\`eme est le suivant (les \'equations
portant sur les composantes transverses de la vitesse sont d\'ecoupl\'ees,
associ\'ees \`a la valeur propre $u$, comme le serait un scalaire simplement
convect\'e et ne sont pas \'ecrites ci-apr\`es)~: 
\begin{equation}
\left\{\begin{array}{lllll}
\displaystyle\frac{\partial\rho}{\partial t} 
&+&\displaystyle\rho\frac{\partial\,u}{\partial\,n} + u\,\frac{\partial\,\rho}{\partial\,n}&=& 0 \\
\displaystyle\frac{\partial u}{\partial t}
&+&\displaystyle u\,\frac{\partial\,u}{\partial\,n}+\frac{1}{\rho}\,\frac{\partial
P}{\partial n} &=& 0 \\
\displaystyle\frac{\partial P}{\partial t}
&+&\displaystyle\gamma\,P\,\frac{\partial\,u}{\partial\,n}+u\,\frac{\partial\,P}{\partial\,n}&=& 0 
\end{array}\right.
\end{equation}

La matrice jacobienne associ\'ee est donc~:
\begin{equation}
\left(\begin{array}{lll}
\displaystyle u & \rho 		& 0				\\
\displaystyle 0 & u    		& \displaystyle\frac{1}{\rho}	\\
\displaystyle 0 & \gamma\, P	& 0 				\\
\end{array}\right)
\end{equation}

Les valeurs propres sont $u$ et $\displaystyle\,u\pm c$ (avec 
$c=\sqrt\frac{\gamma\,P}{\rho}$). Le rayon spectral est donc 
$|u|+c$ et le param\`etre de d\'ecentrement s'en d\'eduit~: 
\begin{equation}
\rho_{rus\,{b}_{ik}} = max\left(|u_i|+c_i,|u_{{b}_{ik}}|+c_{{b}_{ik}}\right)
\end{equation}


%---------------------------------
\subsubsection{Expression des flux convectifs}
%---------------------------------

Les flux convectifs calcul\'es par le sch\'ema de Rusanov 
pour les variables masse, quantit\'e de mouvement
et \'energie repr\'esentent donc la discr\'etisation des termes suivants~: 
\begin{equation}
\left\{\begin{array}{l}
\displaystyle\dive(\vect{Q})\\
\displaystyle\divv(\vect{u}\otimes\vect{Q})+\grad\,P\\
\displaystyle\dive\left(\vect{Q}\,(e+\frac{P}{\rho})\right)
\end{array}\right.
\end{equation}

Pour une face de bord $ik$ adjacente \`a la cellule $i$ et 
avec la valeur pr\'ec\'edente de $\rho_{rus\,{b}_{ik}}$, on a~:
\begin{equation}
\left\{\begin{array}{lll}
\displaystyle\int_{S_{\,{b}_{ik}}}\vect{Q}\cdot\vect{n}\,dS
&=&
\displaystyle\frac{1}{2}\left(
(\vect{Q}_i+\vect{Q}_{\,{b}_{ik}})\cdot\vect{n}_{\,{b}_{ik}}\right)\,S_{\,{b}_{ik}}\\
&&\displaystyle \qquad\qquad\qquad\qquad\qquad\qquad
-\frac{1}{2}\,\rho_{rus\,{b}_{ik}}
\left(\rho_{\,{b}_{ik}}-\rho_i \right)\,S_{\,{b}_{ik}}\\
%
\displaystyle\int_{S_{\,{b}_{ik}}}(\vect{u}\otimes\vect{Q}+\grad\,P)\cdot\vect{n}\,dS
&=&
\displaystyle\frac{1}{2}\left(
 \vect{u}_i(\vect{Q}_i\cdot\vect{n}_{\,{b}_{ik}})
+P_i\,\vect{n}_{\,{b}_{ik}}
+\vect{u}_{\,{b}_{ik}}(\vect{Q}_{\,{b}_{ik}}\cdot\vect{n}_{\,{b}_{ik}})
+P_{\,{b}_{ik}}\,\vect{n}_{\,{b}_{ik}}\right)\,S_{\,{b}_{ik}}\\
&&\displaystyle \qquad\qquad\qquad\qquad\qquad\qquad
-\frac{1}{2}\,\rho_{rus\,{b}_{ik}}
\left(\vect{Q}_{\,{b}_{ik}}-\vect{Q}_i \right)\,S_{\,{b}_{ik}}\\
%
\displaystyle\int_{S_{\,{b}_{ik}}}(e+\frac{P}{\rho})\,\vect{Q}\cdot\vect{n}\,dS
&=&
\displaystyle\frac{1}{2}\left(
(e_i+\frac{P_i}{\rho_i})\,(\vect{Q}_i\cdot\vect{n}_{\,{b}_{ik}})
+(e_{\,{b}_{ik}}+\frac{P_{\,{b}_{ik}}}{\rho_{\,{b}_{ik}}})(\vect{Q}_{\,{b}_{ik}}\cdot\vect{n}_{\,{b}_{ik}})
\right)\,S_{\,{b}_{ik}}\\
&&\displaystyle \qquad\qquad\qquad\qquad\qquad\qquad
-\frac{1}{2}\,\rho_{rus\,{b}_{ik}}
\left(\rho_{\,{b}_{ik}}\,e_{\,{b}_{ik}}-\rho_i\,e_i \right)\,S_{\,{b}_{ik}}\\
\end{array}\right.
\end{equation}



%=================================
\subsection{Conditions aux limites pour le flux diffusif d'\'energie}
\label{Cfbl_Cfxtcl_section_cl_flux_diffusif_energie_cfener}
%=================================

%---------------------------------
\subsubsection{Rappel}
%---------------------------------

Pour le flux de diffusion d'\'energie, les conditions aux limites sont
impos\'ees de mani\`ere similaire \`a ce qui est d\'ecrit dans 
la documentation de \fort{clptur} et de 
\fort{condli}. La figure~(\ref{Cfbl_Cfxtcl_fig_flux_cfxtcl}) rappelle quelques notations
usuelles et l'\'equation~(\ref{Cfbl_Cfxtcl_eq_flux_cfxtcl}) traduit la conservation du flux
normal au bord pour la variable $f$. 

\begin{figure}[htp]
\centerline{\includegraphics[height=7cm]{../Cfbl/Cfxtcl/Images/fluxbord.pdf}}
\caption{\label{Cfbl_Cfxtcl_fig_flux_cfxtcl}Cellule de bord.}
\end{figure}

\begin{equation}\label{Cfbl_Cfxtcl_eq_flux_cfxtcl}
\begin{array}{l}
    \underbrace{h_{int}(f_{b,int}-f_{I'})}_{\phi_{int}}
  = \underbrace{h_{b}(f_{b,ext}-f_{I'})}_{\phi_{b}}
  = \left\{\begin{array}{ll}
    \underbrace{h_{imp,ext}(f_{imp,ext}-f_{b,ext})}_{\phi_{\text{\it r\'eel
impos\'e}}} &\text{(condition de Dirichlet)}\\
    \underbrace{\phi_{\text{\it imp,ext}}}_{\phi_{\text{\it r\'eel impos\'e}}}
            &\text{(condition de Neumann)}
           \end{array}\right.
\end{array}
\end{equation}


L'\'equation~(\ref{Cfbl_Cfxtcl_eq_fbint_cfxtcl}) rappelle la formulation des
conditions aux limites pour une variable $f$. 
\begin{equation}\label{Cfbl_Cfxtcl_eq_fbint_cfxtcl}
f_{b,int}
  = \left\{\begin{array}{cccccl}
    \displaystyle\frac{h_{imp,ext}}{h_{int}+h_r h_{imp,ext} }&f_{imp,ext}&+&
    \displaystyle\frac{h_{int}+h_{imp,ext}(h_r-1)}{h_{int}+h_r h_{imp,ext} }&f_{I'} 
                         &\text{(condition de Dirichlet)}\\
    \displaystyle\frac{1}{h_{int}}&\phi_{\text{\it imp,ext}}&+&
    \ &f_{I'} 
            &\text{(condition de Neumann)}
           \end{array}\right.
\end{equation}

Les coefficients d'\'echange sont d\'efinis comme suit\footnote{On rappelle que, comme
dans \fort{condli}, $\alpha$ d\'esigne $\lambda+C_p\,\frac{\mu_t}{\sigma_t}$ 
si $f$ est la temp\'erature,
$\frac{\lambda}{C_p}+\frac{\mu_t}{\sigma_t}$ si $f$ repr\'esente l'enthalpie. 
Le coefficient $C$ repr\'esente $C_p$ pour la temp\'erature et vaut
$1$ pour l'enthalpie. La grandeur adimensionnelle $f^+$ est obtenue par
application d'un principe de similitude en paroi~: pour la temp\'erature, 
elle d\'epend du nombre de
Prandlt mol\'eculaire, du nombre de Prandtl turbulent et de la distance adimensionnelle \`a la paroi $y^+$ dans la cellule de bord.}~: 
\begin{equation}
\left\{\begin{array}{lll}
h_{int}&=&\displaystyle\frac{\alpha}{\overline{I'F}}\\
h_r&=&\displaystyle\frac{h_{int}}{h_{b}} \\
h_b&=&\displaystyle\frac{\phi_b}{f_{b,ext}-f_{I'}}=\frac{\rho\,C\,u_k}{f^+_{I'}}
\end{array}\right.
\end{equation}

Dans \CS, on note les conditions aux limites de mani\`ere g\'en\'erale sous
la forme suivante~: 
\begin{equation}
f_{b,int}=A_b + B_b\,f_{I'}
\end{equation}
avec $A_b$ et $B_b$ d\'efinis selon le type des conditions~:
\begin{equation}
\begin{array}{c}
\text{Dirichlet}\left\{\begin{array}{ll}
    A_b = &\displaystyle\frac{h_{imp,ext}}{h_{int}+h_r h_{imp,ext} } f_{imp,ext}\\
    B_b = &\displaystyle\frac{h_{int}+h_{imp,ext}(h_r-1)}{h_{int}+h_r h_{imp,ext} }
                  \end{array}\right.                                         
\text{\ \  Neumann}\left\{\begin{array}{ll}
    A_b = &\displaystyle\frac{1}{h_{int}}\phi_{\text{\it imp,ext}}\\
    B_b = &1
                  \end{array}\right.
\end{array}
\end{equation}

%---------------------------------
\subsubsection{Flux diffusif d'\'energie}
%---------------------------------

Dans le module compressible, on r\'esout une \'equation sur l'\'energie, qui s'\'ecrit, si
l'on excepte tous les termes hormis le flux de diffusion et le terme
instationnaire, pour faciliter la pr\'esentation~: 

\begin{equation}
\begin{array}{lll}
\displaystyle\frac{\partial \rho e}{\partial t} &=& - \dive{\,\vect{\Phi}_s}\\
&=& \displaystyle\dive{(K\,\grad{T})} \text{\ \ avec \ \ } K=\lambda+C_p\,\frac{\mu_t}{\sigma_t}\\
&=& \displaystyle\dive{\left(K\,\grad{\frac{e-\frac{1}{2}\,u^2-\varepsilon_{sup}}{C_v}}\right)} \\
&=& \displaystyle\dive{\left(\frac{K}{C_v}\,\grad{(e-\frac{1}{2}\,u^2-\varepsilon_{sup})}\right)} \text{\ \
si \ \ } C_v \text{\ est constant}\\
&=& \displaystyle\dive{\left(\frac{K}{C_v}\,\grad\,e\right)}
-\dive{\left(\frac{K}{C_v}\,\grad{(\frac{1}{2}\,u^2+\varepsilon_{sup})}\right)} \\

\end{array}
\end{equation}

La d\'ecomposition en $e$ et $\frac{1}{2}\,u^2+\varepsilon_{sup}$ est purement
math\'ematique (elle r\'esulte du fait que l'on r\'esout en \'energie alors que
le flux thermique s'exprime en fonction de la temp\'erature). Aussi,  pour imposer un
flux de bord ou une temp\'erature de bord (ce qui revient au m\^eme puisque l'on
impose toujours finalement la conservation du flux normal), on {\it choisit}
de reporter la totalit\'e de la condition \`a la limite sur le terme
$\displaystyle\frac{K}{C_v}\,\grad\,e$
et donc d'annuler le flux associ\'e au terme 
$\displaystyle\frac{K}{C_v}\,\grad{(\frac{1}{2}\,u^2+\varepsilon_{sup})}$ 
(en pratique, pour l'annuler, on se contente de ne pas l'ajouter  
au second membre de l'\'equation). Conform\'ement \`a l'approche retenue dans \CS\ et
rappel\'ee pr\'ec\'edemment, on d\'eterminera donc une valeur de bord {\it
fictive} de l'\'energie qui permette de reconstruire le flux diffusif total 
attendu \`a partir
de la discr\'etisation du seul terme $\displaystyle\frac{K}{C_v}\,\grad\,e$. 

Remarque : dans la version 1.2.0, 
on utilise $\displaystyle 
\frac{K}{C_v}=\left(\frac{\lambda}{C_v}+\frac{\mu_t}{\sigma_t}\right)$, \`a 
partir de 1.2.1, on utilise la valeur  $\displaystyle 
\frac{K}{C_v}=\left(\frac{\lambda}{C_v}+\frac{C_p}{C_v}\frac{\mu_t}{\sigma_t}\right)$.
On notera que le nombre de Prandtl turbulent $\sigma_t$ est associ\'e \`a la variable
r\'esolue et peut \^etre fix\'e par l'utilisateur. 


%---------------------------------
\subsubsection{Condition de Neumann}
%---------------------------------

La conservation du flux s'\'ecrit~: 

\begin{equation}
    \underbrace{h_{int}(e_{b,int}-e_{I'})}_{\phi_{int}}
    =\underbrace{\phi_{\text{\it imp,ext}}}_{\phi_{\text{\it r\'eel impos\'e}}}
\end{equation}

On a donc dans ce cas~:
\begin{equation}
\left\{\begin{array}{lll}
  A_b &= &\displaystyle\frac{1}{h_{int}}\phi_{\text{\it imp,ext}}\\
  B_b &= &1
\end{array}\right.
\end{equation}


%---------------------------------
\subsubsection{Condition de Dirichlet}
%---------------------------------

On suppose que la condition de Dirichlet porte sur la temp\'erature $T_{b,ext}$.


La conservation du flux s'\'ecrit~: 
\begin{equation}\label{Cfbl_Cfxtcl_eq_conservation_flux_cfxtcl}
    \underbrace{h_{int}(e_{b,int}-e_{I'})}_{\phi_{int}\text{\ (forme num\'erique
du flux)}}
  = \underbrace{h_{b}(T_{b,ext}-T_{I'})}_{\phi_{b}\text{ qui int\`egre l'effet
de couche limite}}
  = 
    \underbrace{h'_{imp,ext}(T_{imp,ext}-T_{b,ext})}_{\phi_{\text{\it r\'eel
impos\'e}}}
\end{equation}

Avec pour les coefficients d'\'echange~:
\begin{equation}
\left\{\begin{array}{lll}
h_{int}&=&\displaystyle\frac{K}{C_v\,\overline{I'F}}\\
h_b&=&\displaystyle\frac{\phi_b}{T_{b,ext}-T_{I'}}=\frac{\rho\,C_p\,u_k}{T^+_{I'}}
\end{array}\right.
\end{equation}

On tire $T_{b,ext}$
de la seconde partie de l'\'egalit\'e~(\ref{Cfbl_Cfxtcl_eq_conservation_flux_cfxtcl}) 
traduisant la conservation du flux~:
\begin{equation}
\displaystyle T_{b,ext} = \frac{h'_{imp,ext}\,T_{imp,ext}+h_b\,T_{I'}}{h_b+h'_{imp,ext}}
\end{equation}

En utilisant cette valeur et la premi\`ere partie de l'\'equation de conservation
du flux~(\ref{Cfbl_Cfxtcl_eq_conservation_flux_cfxtcl}), on obtient~:
\begin{equation}
e_{b,int} = \frac{h_b\,h'_{imp,ext}}{h_{int}\,(h_b+h'_{imp,ext})}\,(T_{imp,ext}-T_{I'})+e_{I'}
\end{equation}

On utilise alors
$\displaystyle T_{I'}=\frac{1}{C_v}\left(e_{I'}-\frac{1}{2}u^2_{i}-\varepsilon_{sup,i}\right)$ pour
\'ecrire (sans reconstruction pour la vitesse et $\varepsilon_{sup}$)~:
\begin{equation}
\displaystyle e_{b,int} = 
\frac{ -\frac{h_b\,h'_{imp,ext}}{C_v}+h_{int}\,(h_b+h'_{imp,ext}) }
     { h_{int}\,(h_b+h'_{imp,ext}) } \,e_{I'}
+\frac{h_b\,h'_{imp,ext}}{h_{int}\,(h_b+h'_{imp,ext})}\,
  \left(T_{imp,ext}+\frac{\frac{1}{2}u^2_{i}+\varepsilon_{sup,i}}{C_v}\right)
\end{equation}


Et on a donc, avec $\displaystyle h'_r=\frac{h_{int}}{\frac{h_b}{C_v}}$~:
\begin{equation}
\displaystyle e_{b,int} = 
\underbrace{\frac{ h'_{imp,ext} }{ C_v\,h_{int}+h'_r\,h'_{imp,ext} }\,
  \left(C_v\,T_{imp,ext}+\frac{1}{2}u^2_{i}+\varepsilon_{sup,i}\right)}_{A_b}
+\underbrace{\frac{ C_v\,h_{int}+h'_{imp,ext}(h'_r-1) }{ C_v\,h_{int}+h'_r\,h'_{imp,ext} }}_{B_b}\,e_{I'}
\end{equation}

Avec ces notations, $h_b$ est le coefficient habituellement calcul\'e pour la 
temp\'erature. 

Le coefficient $h'_{imp,ext}$ est le coefficient d'\'echange externe qui est 
impos\'e pour la temp\'erature\footnote{Le coefficient $h'_{imp,ext}$  
est utile pour les cas o\`u l'on 
souhaite relaxer la condition \`a la limite~: 
pour la temp\'erature, cela correspond \`a imposer une valeur sur la face 
externe d'une paroi unidimensionnelle id\'eale, sans inertie, 
caract\'eris\'ee par un simple coefficient d'\'echange.}. 
Pour obtenir l'\'equivalent dimensionnel de $h'_{imp,ext}$ pour l'\'energie, 
il faut diviser sa valeur par $C_v$ (ce qui ne fait pas de diff\'erence dans 
la majorit\'e des cas, car il est habituellement pris infini).   



 
 


 


