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

Pour utiliser la méthode, on se place tout d'abord dans un repère local défini
de manière à ce que le plan $(0yz)$, où sont injectés les vortex, soit confondu
avec le plan d'entrée du calcul (voir figure \ref{Base_Vortex_entree}). 

\begin{figure}[h]
\centerline{\includegraphics[height=6cm]{../Base/Vortex/Images/entree.pdf}}
\caption{\label{Base_Vortex_entree} Définiton des différentes grandeurs dans le repère local
correspondant à l'entrée d'une conduite de section carrée.} 
\end{figure}

$u$, $v$ et $w$  sont les composantes de la vitesse fluctuante (principale et
transverse) dans ce plan, et
$\displaystyle \omega(y,z) = \frac{\partial w}{\partial y}-\frac{\partial v}{\partial z}$
la vorticité dans la direction
normale au plan d'entrée. $\overline{U}(y,z)$ représente ici la vitesse
principale moyenne imposée par l'utilisateur dans le plan d'entrée. 

Chaque vortex $p$ va être caractérisé par sa fonction de forme $\xi$ (identique
pour tous les vortex), sa
circulation $\Gamma_p$, son rayon $\sigma_p$ et les coordonnées $(y_p,z_p)$ du
point $P$ où est situé le vortex dans le plan $(0yz)$. 

Pour cela, on suppose que la vorticité générée par un vortex $p$ au point $M$ de
coordonnée $(y,z)$ s'écrit : 
\begin{equation}\notag
\omega_p(y,z)= \Gamma_p \, \xi_{\sigma_p}(r)
\end{equation}
où $r = \sqrt{(y-y_p)^2+(z-z_p)^2}$ est la distance séparant le point $M$ du point $P$.

Dans la méthode implantée, la fonction de forme est de type gaussienne modifiée :
\begin{equation}\notag
\displaystyle
\xi_\sigma (r) = \frac{1}{2\pi \sigma^2} 
\left(2 e^{-\frac{r^2}{2\sigma^2}}-1\right) e^{-\frac{r^2}{2\sigma^2}}
\end{equation}

Le champ de vitesse induit par cette distribution de vorticité s'obtient par
inversion des deux équations de poisson suivantes qui sont déduites de la
condition d'incompressibilité dans la plan\footnote{\textit{i.e}
$\displaystyle \frac{\partial v}{\partial y}+\frac{\partial w}{\partial w} = 0$} :
\begin{equation}\notag
\begin{array}{lcr}
\displaystyle
\frac{\partial \omega}{\partial y} = \Delta w
&
\text{    et    }
&
\displaystyle
\frac{\partial \omega}{\partial y} = -\Delta v
\\
\end{array}
\end{equation}

Dans le cas général, ce système peut être intégré à l'aide de la formule de Biot et Savart.

Dans le cas d'une distribution de vorticité de type gaussienne modifiée, les
composantes de vitesse vérifient : 
\begin{equation}\notag
\left\{
\begin{array}{c}
\displaystyle
v_p(y,x) = - \frac{1}{2\pi} \frac{(z-z_p)}{r^2}\left(1 -
e^{-\frac{r^2}{2\sigma^2}}\right)\,e^{-\frac{r^2}{2\sigma^2}} 
\\
\displaystyle
w_p(y,z) = \frac{1}{2\pi} \frac{(y-y_p)}{r^2}\left(1 -e^{-\frac{r^2}{2\sigma^2}}
\right)\,e^{-\frac{r^2}{2\sigma^2}} 
\end{array}
\right.
\end{equation}

Ces relations s'étendent de façon immédiate au cas de $N$ vortex distincts.
Le champ de vitesse induit par la distribution de vorticité 
\begin{equation}
\omega(y,z) = \sum_{p=1}^N \Gamma_p \, \xi_{\sigma_p}(r)
\end{equation}
vaut au point $M$ :
\begin{equation}\notag
\begin{array}{lcr}
v(x,y) = \sum_{p=1}^N \Gamma_p\, v_p(y,z) 
&
\text{    et    }
&
w(y,z) = \sum_{p=1}^N \Gamma_p\, w_p(y,z)
\\
\label{Base_Vortex_compvit}
\end{array}
\end{equation}
%================================
\subsection{Paramètres physiques}
%================================

%-------------------------------
\subsubsection{Marche en temps}
%-------------------------------
La position initiale de chaque vortex est tirée de manière aléatoire. On calcul
les déplacements successifs de chacun des vortex dans le plan d'entrée par
intégration explicite du champ de vitesse lagrangien : 
\begin{equation}\notag
\begin{array}{lcr}
\displaystyle
\frac{dy_p}{dt} = V(y,z)
&
\text{    et    }
&
\displaystyle
\frac{dz_p}{dt} = W(y,z)
\\
\end{array}
\end{equation}
Les vortex sont alors assimilés à des particules ponctuelles qui sont convectées
par le champ $(V,W)$. Ce champ peut être imposé par des tirages aléatoires ou
bien déduit de la vitesse induite par les vortex dans le plan d'entrée. Dans ce
cas $V(x,y) = \overline{V}(y,z) + v (y,z)$ et $W(y,z)= \overline{W}(y,z) +
w(y,z)$ où $\overline{V}$ et $\overline{W}$ sont les composantes de la vitesse
transverse moyenne qu'impose l'utilisateur à l'aide des fichiers de données. 

%---------------------------------------------------
\subsubsection{Intensité et durée de vie des vortex}
%---------------------------------------------------
Il serait possible, à partir de l'équation de transport de la vorticité,
d'obtenir un modèle d'évolution pour l'intensité du vecteur tourbillon
$\omega_p$ associé à chacun des vortex. En pratique, on préfère utiliser un
modèle simplifié dans lequel la circulation des tourbillons ne dépend que de la
postion de ces derniers à l'instant considéré. La circulation initiale de chaque
vortex est alors obtenue à partir de la relation : 
\begin{equation}\notag
|\Gamma_p| = 4 \sqrt{\frac{\pi\,S\,k}{3N\,[2ln(3)-3ln(2)]}}
\end{equation}
où $S$ est la surface du plan d'entrée, $N$ le nombre de vortex, et $k$
l'énergie cinétique turbulente au point où se trouve le vortex à l'instant
considéré. Le signe de $\Gamma_p$ correspond au sens de rotation du vortex et
est tiré aléatoirement. 

Ce paramètre est celui qui contrôle l'intensité des fluctuations. Sa dépendance
en $k$ exprime que, plus l'écoulement est turbulent, plus les vortex sont
intenses. La valeur de $k$ est spécifiée par
l'utilisateur. Elle peut être constante ou imposée à partir de profils d'énergie
cinétique turbulente en entrée. 

Pour éviter que des structures trop allongées ne se développent au niveau de
l'entrée, l'utilisateur doit également spécifier un temps limites $\tau_p$ au
bout duquel le vortex $p$ va être détruit. Ce temps $\tau_p$ peut être pris
constant ou estimé au moyen de la relation : 
\begin{equation}\notag
\tau_p = \frac{5 C_{\mu}k^{\frac{3}{2}}}{\varepsilon\,\overline{U}}
\end{equation}

$\overline{U}$ et $\varepsilon$ représentent respectivement la vitesse moyenne
principale et la dissipation turbulente au point où est initialement généré le
vortex ($C_{\mu}=0,09$). 
\\
Lorsque le temps écoulé depuis la création du vortex $p$ est supérieur à
$\tau_p$, le vortex est détruit et un nouveau vortex généré (sa position et le
signe de $\Gamma_p$ sont tirés de façon aléatoire). 

%-------------------------------- 
\subsubsection{Taille des vortex}
%--------------------------------
La taille des vortex peut être prise constante, ou calculée à partir des
relations :
\begin{equation}\notag
\begin{array}{ccc}
\displaystyle
\sigma = \frac{C_{\mu}^{\frac{3}{4}}k^{\frac{3}{2}}}{\varepsilon} 
& \text{    ou    } &
\sigma = max[L_t,L_k]
\\
\end{array}
\end{equation}
avec:
\begin{equation}\notag
\begin{array}{ccc}
\displaystyle
L_t = \sqrt{\left( \frac{5 \nu k}{\varepsilon} \right)} 
& \text{    et    } & 
\displaystyle
L_k = 200\, \left(\frac{\nu^3}{\varepsilon}\right)^{\frac{1}{4}}
\end{array}
\end{equation}
où $\nu$, $k$ et $\varepsilon$ sont la viscosité dynamique, l'énergie cinétique
turbulente et la dissipation turbulente au point où se trouve le vortex. 

Dans tous les cas, la taille du vortex doit être supérieure à la taille des
mailles en entrée afin que le vortex soit effectivement simulé. 

%==================================
\subsection{Conditions aux limites}
%==================================
Le champ de vitesse généré à l'aide de la relation \ref{Base_Vortex_compvit} ne tient pas
compte {\em a priori} des conditions aux limites appliquées sur les bords du plan
d'entrée. Pour obtenir des valeurs de la vitesse qui soient cohérentes sur les
frontières du domaine d'entrée, des ``vortex images'', pseudo-vortex situés en
dehors du domaine d'entrée, sont générés à des positions particulières et leur
champ de vitesse associé est superposé au champ précédemment calculé.\\
Seuls les cas d'une conduite rectangulaire et d'une conduite circulaire
permettent la génération de ces pseudo-vortex.
On distingue pour cela trois types de conditions aux limites. 

\begin{figure}[h]
\centerline{\includegraphics[height=6cm]{../Base/Vortex/Images/condlimite.pdf}}
\caption{\label{Base_Vortex_condli} Principe de génération des ``vortex images'' suivant le
type de conditions aux limites dans une conduite carrée.} 
\end{figure}

%----------------------------------
\subsubsection{Condition de paroi}
%----------------------------------
On crée, pour chaque vortex $P$ contenu dans le plan d'entrée, un vortex image
$P'$ identique à $P$ (\textit{i.e} de même caractéristiques) et symétrique de $P$
par rapport au
point $J$ ($J$ étant la projection orthogonalement à la paroi du point $M$
correspondant au centre de la face où l'on cherche à calculer la vitesse). La
figure \ref{Base_Vortex_condli} illustre la technique dans le cas d'une conduite
carrée. Dans ce cas les coordonnées du vortex situé en $P'$ vérifient
$(y_{p'}+y_{p})/2 = y_{J}$ et $(z_{p'}+ z_{p})/2 = z_{J}$. Le champ de vitesse
perçu depuis le point $M$ au niveau du point $J$ est nul, ce qui est bien
l'effet recherché. 

%------------------------------------
\subsubsection{Condition de symétrie}
%-------------------------------------
La technique est identique à celle utilisée pour les conditions de paroi, mais
seule la composante pour la vitesse normale au bord est modifiée dans ce cas. 

%---------------------------------------
\subsubsection{Condition de périodicité}
%---------------------------------------
On crée pour chaque vortex, un vortex images $P'$ identique à $P$ mais translaté
d'une quantité $L$ correspondant à la longueur qui sépare les deux plans de la
section d'entrée où sont appliquées les conditions de périodicité. Dans le cas
où il y a deux directions de périodicité, on crée deux vortex image.

%=============================================
\subsection{Composante de vitesse principale}
%=============================================
La méthode des vortex ne générant pas de fluctuation $u$ de la vitesse dans la
direction principale, la dernière composante est calculée à partir d'une
équation de Langevin. Les coefficients de cette équation sont déterminés par
identification des expressions obtenues pour les contraintes de Reynolds en
$R_{ij}-\varepsilon$. Dans le cas d'un écoulement en canal plan, cette technique
conduit à l'équation : 
\begin{equation}\notag
\displaystyle
\frac{du}{dt} = - \frac{C_1}{2T} u + \left(\frac{2}{3}C_2-1\right)\frac{\partial
U}{\partial y} v + \sqrt{C_0\varepsilon} dW_i 
\end{equation}

avec $\displaystyle T=\frac{k}{\varepsilon}$, $C_1 = 1,8$, $C_2 = 0,6$,
$C_0=\frac{14}{15}$, et $dW_i$ une variable alétoire Gaussienne de variance
$\sqrt{dt}$. 

En pratique, l'équation de Langevin n'améliore pas vraiment les résultats. Elle
n'est utilisée que dans le cas de conduites circulaires. 
