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

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

\hypertarget{vortex}{}

\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Fonction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ce sous-programme est ddi  la gnration des conditions d'entre
turbulente utilises en LES.


La mthode des vortex est base sur une approche de tourbillons
ponctuels. L'ide de la mthode consiste  injecter des tourbillons 2D dans le
plan d'entre du calcul, puis  calculer le champ de vitesse induit par ces
tourbillons au centre des faces d'entre.

See the \doxygenfile{vortex_8f90.html}{programmers reference of the dedicated subroutine} for further details.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discr\'etisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Pour utiliser la mthode, on se place tout d'abord dans un repre local dfini
de manire  ce que le plan $(0yz)$, o sont injects les vortex, soit confondu
avec le plan d'entre du calcul (voir figure \ref{Base_Vortex_entree}).

\begin{figure}[h]
\centerline{\includegraphics[height=6cm]{entree}}
\caption{\label{Base_Vortex_entree} Dfiniton des diffrentes grandeurs dans le repre local
correspondant  l'entre d'une conduite de section carre.}
\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'entre. $\overline{U}(y,z)$ reprsente ici la vitesse
principale moyenne impose par l'utilisateur dans le plan d'entre.

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

Pour cela, on suppose que la vorticit gnre par un vortex $p$ au point $M$ de
coordonne $(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 sparant le point $M$ du point $P$.

Dans la mthode implante, la fonction de forme est de type gaussienne modifie :
\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 dduites 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 gnral, ce systme peut tre intgr  l'aide de la formule de Biot et Savart.

Dans le cas d'une distribution de vorticit de type gaussienne modifie, les
composantes de vitesse vrifient :
\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 faon immdiate 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*{Paramtres physiques}
%================================

%-------------------------------
\subsubsection*{Marche en temps}
%-------------------------------
La position initiale de chaque vortex est tire de manire alatoire. On calcul
les dplacements successifs de chacun des vortex dans le plan d'entre par
intgration 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 assimils  des particules ponctuelles qui sont convectes
par le champ $(V,W)$. Ce champ peut tre impos par des tirages alatoires ou
bien dduit de la vitesse induite par les vortex dans le plan d'entre. 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 donnes.

%---------------------------------------------------
\subsubsection*{Intensit et dure de vie des vortex}
%---------------------------------------------------
Il serait possible,  partir de l'quation de transport de la vorticit,
d'obtenir un modle d'volution pour l'intensit du vecteur tourbillon
$\omega_p$ associ  chacun des vortex. En pratique, on prfre utiliser un
modle simplifi dans lequel la circulation des tourbillons ne dpend que de la
postion de ces derniers  l'instant considr. 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'entre, $N$ le nombre de vortex, et $k$
l'nergie cintique turbulente au point o se trouve le vortex  l'instant
considr. Le signe de $\Gamma_p$ correspond au sens de rotation du vortex et
est tir alatoirement.

Ce paramtre est celui qui contrle l'intensit des fluctuations. Sa dpendance
en $k$ exprime que, plus l'coulement est turbulent, plus les vortex sont
intenses. La valeur de $k$ est spcifie par
l'utilisateur. Elle peut tre constante ou impose  partir de profils d'nergie
cintique turbulente en entre.

Pour viter que des structures trop allonges ne se dveloppent au niveau de
l'entre, l'utilisateur doit galement spcifier un temps limites $\tau_p$ au
bout duquel le vortex $p$ va tre dtruit. 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$ reprsentent respectivement la vitesse moyenne
principale et la dissipation turbulente au point o est initialement gnr le
vortex ($C_{\mu}=0,09$).
\\
Lorsque le temps coul depuis la cration du vortex $p$ est suprieur 
$\tau_p$, le vortex est dtruit et un nouveau vortex gnr (sa position et le
signe de $\Gamma_p$ sont tirs de faon alatoire).

%--------------------------------
\subsubsection*{Taille des vortex}
%--------------------------------
La taille des vortex peut tre prise constante, ou calcule  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 cintique
turbulente et la dissipation turbulente au point o se trouve le vortex.

Dans tous les cas, la taille du vortex doit tre suprieure  la taille des
mailles en entre afin que le vortex soit effectivement simul.

%==================================
\subsection*{Conditions aux limites}
%==================================
Le champ de vitesse gnr  l'aide de la relation \ref{Base_Vortex_compvit} ne tient pas
compte {\em a priori} des conditions aux limites appliques sur les bords du plan
d'entre. Pour obtenir des valeurs de la vitesse qui soient cohrentes sur les
frontires du domaine d'entre, des ``vortex images'', pseudo-vortex situs en
dehors du domaine d'entre, sont gnrs  des positions particulires et leur
champ de vitesse associ est superpos au champ prcdemment calcul.\\
Seuls les cas d'une conduite rectangulaire et d'une conduite circulaire
permettent la gnration de ces pseudo-vortex.
On distingue pour cela trois types de conditions aux limites.

\begin{figure}[h]
\centerline{\includegraphics[height=6cm]{condlimite}}
\caption{\label{Base_Vortex_condli} Principe de gnration des ``vortex images'' suivant le
type de conditions aux limites dans une conduite carre.}
\end{figure}

%----------------------------------
\subsubsection*{Condition de paroi}
%----------------------------------
On cre, pour chaque vortex $P$ contenu dans le plan d'entre, un vortex image
$P'$ identique  $P$ (\textit{i.e} de mme caractristiques) et symtrique 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
carre. Dans ce cas les coordonnes du vortex situ en $P'$ vrifient
$(y_{p'}+y_{p})/2 = y_{J}$ et $(z_{p'}+ z_{p})/2 = z_{J}$. Le champ de vitesse
peru depuis le point $M$ au niveau du point $J$ est nul, ce qui est bien
l'effet recherch.

%------------------------------------
\subsubsection*{Condition de symtrie}
%-------------------------------------
La technique est identique  celle utilise pour les conditions de paroi, mais
seule la composante pour la vitesse normale au bord est modifie dans ce cas.

%---------------------------------------
\subsubsection*{Condition de priodicit}
%---------------------------------------
On cre pour chaque vortex, un vortex images $P'$ identique  $P$ mais translat
d'une quantit $L$ correspondant  la longueur qui spare les deux plans de la
section d'entre o sont appliques les conditions de priodicit. Dans le cas
o il y a deux directions de priodicit, on cre deux vortex image.

%=============================================
\subsection*{Composante de vitesse principale}
%=============================================
La mthode des vortex ne gnrant pas de fluctuation $u$ de la vitesse dans la
direction principale, la dernire composante est calcule  partir d'une
quation de Langevin. Les coefficients de cette quation sont dtermins 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 altoire Gaussienne de variance
$\sqrt{dt}$.

En pratique, l'quation de Langevin n'amliore pas vraiment les rsultats. Elle
n'est utilise que dans le cas de conduites circulaires.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{itemize}
\item[$\star$] Aprs une tape de prparation de la mmoire (\fort{memvor}), on
repre dans \fort{usvort} les faces d'entre pour lesquelles la mthode va tre
utilise.
\item[$\star$] Vrification des dimensions rentres (\fort{vervor}).
\\
\item[$\star$] Le sous-programme \fort{vorpre} se charge ensuite de prparer le
calcul (transmission de la gomtrie des entres  tous les processeurs en cas
de paralllisme, et construction d'un tableau de connectivit). Le
sous-programme procde ainsi :
\\
\begin{itemize}

\item[$\bullet$] On compte, pour chaque entre \var{IENT}, le nombre de faces o
est applique la mthode. Celui-ci est stock dans le tableau
\var{ICVOR(IENT)}. Un passage dans la sous-routine \fort{memvor} (avec
\var{IAPPEL = 2}) permet d'allouer la mmoire ncessaire  cette phase de
prparation.

\item[$\bullet$] Pour chaque processeur, on stocke les coordonnes des faces
d'entre repres prcdemment dans les tableaux de travail
\var{RA(IW1X),RA(IW1Y),RA(IW1Z),...}

\item[$\bullet$]  On regarde ensuite pour chaque processeur (boucle
\var{IPROC=1, NRANGP-1}), si le processeur \var{IPROC} a des donnes  envoyer
aux autres processeurs (afin que tous disposent des coordonnes).
\begin{itemize}
\item Si c'est le cas : \var{ICVOR(IENT)>0}, et on place les donnes  envoyer
dans les tableaux de travail \var{RA(IW2X),RA(IW2Y),RA(IW2Z),...}. La valeur
\var{NCOMV = ICVOR(IENT)} correspond alors  la longueur des tableaux  envoyer.
\item Sinon, on ne fait rien et \var{NCOM=0}.
\end{itemize}
\item[$\bullet$] Le processeur numro \var{IPROC} distribue  tous les autres
processeurs la valeur \var{NCOM}. Si \var{NCOM > 0}, il envoie galement les
donnes contenues dans les tableaux de travails \var{RA(IW2X),...}. Ces donnes
sont ensuite stockes par tous les processeurs dans les tableaux
\var{RA(IXYZV+III),...} afin de librer les tableaux de travail pour la
communication suivante, et l'indice \var{III = III + NCOM} est incrment de
manire  ranger les valeurs de faon chronologique.
\\\\
$\rightarrow$ Au final de la boucle sur \var{IPROC}, chaque processeur dispose
des coordonnes des faces d'entre pour lesquelles la mthode va tre utilise,
et il est donc simple de construire la connectivit.
\\
\item[$\bullet$] Construction de la connectivit. Au final, la vitesse au centre
de la \var{II} me face d'entre utilisant la mthode est contenue  la
\var{IA(IIFAGL+II)} me ligne du tableau \var{RA(IUVORT)}.

\item[$\bullet$] La routine se termine par un appel au sous-programme
\fort{memvor} ( avec \var{IAPPEL = 3}) afin de rserver la mmoire utile  la
mthode des vortex.
\end{itemize}
\end{itemize}
\bigskip

Cette phase d'initialisation est ralise une seule fois au dbut du
calcul. C'est aprs cette phase seulement que commence la mthode des vortex
proprement dite.
\\
\begin{itemize}
\item[$\star$] Initialisation des variables avant intervention utilisateur (\fort{inivor}).
\item[$\star$] Appel au sous-programme utilisateur \fort{usvort} (\var{IAPPEL = 2}).
\item[$\star$] Vrification des paramtres rentrs (\fort{vervor}).
\item[$\star$] Calcul de la vitesse par la mthode des vortex (\fort{vortex})
\begin{itemize}
\item[$\bullet$] Initialisation du calcul gnration du champ initial par appel
au sous-programme \fort{vorini} :

\begin{itemize}
\item Construction du repre local (et calcul de l'quation du plan d'entre
suivant les cas), localisation du centre de l'entre, et transformation des
coordonnes de l'entre dans le repre local. Les tableaux \var{YZCEL(II,1)} et
\var{YZCEL(II,2)} contiennent les coordonnes des faces du plan d'entre une
fois ramenes dans le repre $(0yz)$ (\var{II} est compris entre 1 et
\var{NCEVOR} o \var{NCEVOR}=\var{ICVOR} reprsente le nombre de faces pour
lesquelles la mthode va tre utilise a cette entre).
\item Lecture du fichier de donnes, et initialisation des tableaux \var{XDAT},
\var{YDAT}, \var{UDAT}, \var{VDAT}, \var{WDAT}, \var{DUYDAT}, \var{KDAT},
\var{EPSDAT}, ...
\item Si on ne fait pas de suite (\var{ISUIVO=0}) ou que l'on rinitialise le
calcul (\var{INITVO=1}), tirage alatoire de la position des vortex et de leur
sens de rotation, ainsi que calcul de leur dure de vie. Les positions sont
stockes dans les tableaux \var{YZVOR(IVOR,1)} et \var{YZVOR(IVOR,2)}
(\var{IVOR} dsignant le numro du vortex).
\item Stockage de la vitesse principale moyenne au centre de la cellule dans le
tableau \var{XU}, et recherche pour chaque vortex, de la face d'entre qui lui
est la plus proche.
\end{itemize}

\item[$\bullet$] Dplacement des vortex par appel au sous-programme \fort{vordep} :
\begin{itemize}
\item Convection des vortex.
\item Traitement des conditions aux limites. Les vortex qui sortent du domaine
de calcul sont replacs  leur position d'origine.
\item Rgnration des vortex ``morts''. Si le temps de vie cumul
\var{TEMPS(II)} du vortex \var{II} est suprieur  sont temps de vie limite
\var{TPSLIM(II)}, alors le vortex est dtruit, et un nouveau vortex est gnr.
\item Recherche pour chaque vortex de la face d'entre qui lui est la plus
proche aprs dplacement (mise  jour du tableau \var{IVORCE}).
\end{itemize}

\item[$\bullet$] Calcul du champ de vitesse induit par appel au sous-programme \fort{vorvit} :
\begin{itemize}
\item Calcul de l'intensit du vortex.
\item Calcul de la taille du vortex.
\item Calcul du champ de vitesse induit par l'ensemble des vortex au centre des
faces d'entre.
\item Traitement suivant les cas, des conditions de priodicit de symtrie et
des conditions de paroi par gnration de vortex images.
\item Ajout de la vitesse moyenne dans les directions transverse aux tableaux
\var{XV} et \var{XW}.
\end{itemize}

\item[$\bullet$] Gnration des fluctuations de vitesse dans la direction
principale par appel au sous-programme \fort{vorlgv}.
\end{itemize}

\item[$\star$] appel au sous-programme \fort{vor2cl} :
\item[$\bullet$] Communication en cas de paralllisme de la vitesse calcule en
entre par le processeur 0 aux autres processeurs.
\item[$\bullet$] Application des conditions aux limites aprs utilisation d'un
changement de repre ventuel.
\end{itemize}


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

Il serait possible de gagner de la mmoire en liberant l'espace alou aux
tableaux \var{IW1X},...,\var{IW2V} aprs le passage dans \fort{vorpre}.