1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
|
%-------------------------------------------------------------------------------
% 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{visecv}
%
\hypertarget{visecv}{}
\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Fonction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dans ce sous-programme sont calcul\'es les termes de gradient transpos et de viscosit
secondaire (fluide newtonien). Ils seront traits en explicite
pur. En effet, si on traitait ces termes en implicite, cela reviendrait
coupler les quations des diffrentes composantes de la vitesse, ce qui n'est
pas compatible avec le sch\'ema de r\'esolution actuel ({\it cf. } sous-programme \fort{navstv}).
See the \doxygenfile{visecv_8f90.html}{programmers reference of the dedicated subroutine} for further details.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discr\'etisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L'intgration des termes de gradient transpos
$\dive\,(\mu_{\,tot}\,^t\,\ggrad(\underline{v}))$ et de viscosit
secondaire \\
$-\displaystyle \frac{2}{3}\,\grad (\mu_{\,tot}\,\dive(\underline{v}))$ est la suivante\footnote{la viscosit de volume
$\kappa$ est suppose nulle, cf. \fort{navstv}} :
\begin{equation}\notag
\begin{array}{llll}
&\displaystyle \int_{\Omega_i}\dive (\mu_{\,tot}
\,^t\,\ggrad(\underline{v})\,)\,d\Omega \\
&= \sum\limits_{l=x,y,z}\left[ \sum\limits_{j \in
Vois(i)} \mu_{\,tot,ij} ((\displaystyle \frac {\partial v_x}{\partial
l})_{\,moy,ij}\,n_{\,ij}^{\,x} +(\frac {\partial v_y}
{\partial l})_{\,moy,ij}\,n_{\,ij}^{\,y}+ \displaystyle(\frac {\partial v_z}
{\partial l})_{\,moy,ij}\,n_{\,ij}^{\,z}) S_{\,ij}\right.\\
&+\left.\displaystyle \sum\limits_{k \in \gamma_b(i)} \mu_{\,tot,\,b_{ik}}
((\frac {\partial v_x}{\partial l})_{\,moy,\,b_{ik}}\,n_{\,b_{ik}}^{\,x} +(\frac
{\partial v_y}{\partial l})_{\,moy,\,b_{ik}}\,n_{\,b_{ik}}^{\,y}+
\displaystyle(\frac {\partial v_z}{\partial
l})_{\,moy,\,b_{ik}}\,n_{\,b_{ik}}^{\,z}) S_{\,b_{ik}}\right] \underline
{e}_{\,l} \\\\
&-\displaystyle \frac{2}{3}\,\int_{\Omega_i}\grad (\mu_{\,tot} \dive\,(\underline{v}))
\,d\Omega \\
&= -\displaystyle \frac{2}{3}\,\sum\limits_{l=x,y,z}\left[\sum\limits_{j \in Vois(i)}(\mu_{\,tot}\,\dive(\underline{v}))_{\,ij} S_{\,ij}^{\,l} + \displaystyle \sum_{k \in \gamma_b(i)}(\mu_{\,tot}\,\dive(\underline{v}))_{\,b_{ik}} S_{\,b_{ik}}^{\,l}\right]\,\underline {e}_{\,l}
\end{array}
\end{equation}
Le terme de viscosit $\mu_{\,tot}$ est calcul en fonction du modle de turbulence utilis :
\begin{equation}\notag
\mu_{\,tot}=
\begin{cases}
\mu + \mu_t & \text{pour les modles viscosit turbulente ou en LES}, \\
\mu& \text{pour les modles au second ordre ou en laminaire}.
\end{cases}
\end{equation}
o $\mu$ et $\mu_t$ repr\'esentent respectivement la viscosit\'e dynamique molculaire et la viscosit dynamique turbulente.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection*{Terme de gradient transpos}
Avant de calculer l'intgration des deux termes, on calcule dans un premier
temps la viscosit totale en fonction du modle de turbulence considr.\\
Les termes calculs dans ce sous-programme, appel par \fort{predvv}, interviennent dans le second membre
de l'quation de quantit de mouvement, et sont donc directement stocks dans le
tableau correspondant \var{TRAV}.\\
Le terme $\dive\,(\mu_{\,tot}\,^t\ggrad(\underline{v}))$ est calcul en
oprant comme suit.\\
On effectue une boucle sur les composantes $v_\alpha$ o $\alpha=x,y,z$ de la
vitesse ($\alpha$ correspond \var{ISOU} dans le code) :
\begin{itemize}
\item on calcule le gradient cellule de $v_{\,\alpha}$ par un appel au
sous-programme \fort{grdcel}.
\item on initialise un tableau nomm $\var{W}_6$ la valeur $1$ pour les cellules
internes, et la valeur $0$ pour les cellules de bord. Ce tableau sert par la suite ne
pas considrer la contribution du terme de gradient transpos sur les cellules
de bord. En effet, on ne sait pas \'ecrire de conditions aux limites correctes
pour le gradient transpos\'e. On pr\'ef\`ere donc tout simplement annuler son
effet sur les cellules de bord (\emph{cf.} paragraphe \ref{Base_Visecv_paragraphe4}).
\item pour chaque direction $l$ ($l=x,y,z$), $l$ correspondant \var{IDIM}
dans le code, on calcule pour chaque cellule $\Omega_i$ dont le centre
correspond la variable \var{II} ou
\var{JJ} (pour les centres voisins) dans le code :\\
\begin{equation}\notag
\begin{array}{ll}
&\var{TRAV}(i,l)=\var{TRAV}(i,l)\\
&+\var{W}_6(i)\ \left[\displaystyle \sum\limits_{j \in
Vois(i)} \mu_{\,tot,\,ij}\,(\frac {\partial v_{\,\alpha}}{\partial l}
)_{\,moy,ij}\,S_{\,ij}^{\,\alpha}\ + \displaystyle \sum_{k \in
\gamma_b(i)} \mu_{\,tot,\, {b_{\,ik}}}\,(\frac {\partial v_{\,\alpha}}{\partial
l})_{\,moy,{b_{\,ik}}}\,S_{\,{b_{\,ik}}}^{\,\alpha}\ \right]\\
&\displaystyle\text{avec }
(\frac{\partial v_{\,\alpha}}{\partial l})_{\,moy,ij}=\frac{1}{2}\left[
(\frac{\partial v_{\,\alpha}}{\partial l})_i
+(\frac{\partial v_{\,\alpha}}{\partial l})_j\right]
\end{array}
\end{equation}
\end{itemize}
Fin de la boucle sur les composantes de la vitesse.\\
\\
\subsection*{Terme de viscosit secondaire}
Le terme de seconde viscosit
$-\displaystyle \frac{2}{3}\grad (\mu_{\,tot} \dive\,(\underline{v}))$
est calcul de la faon suivante :
\begin{itemize}
\item on calcule la valeur de la vitesse sur la face $ij$ en divisant le flux
masse connu la face par la densit moyenne $\rho_{\,moy,\,ij}$ de la face
($\rho_{\,moy,\,ij}=\displaystyle \frac{\rho_{\,i}+\rho_{\,j}}{2}$).
\item on calcule ensuite l'intgrale volumique de la divergence de la vitesse sur chaque cellule en appelant le sous-programme
\fort{divmas}.
\item on calcule alors pour chaque cellule $\Omega_i$ le terme $-\displaystyle
\frac{2}{3}(\mu_{\,tot}\,\dive(\underline{v}))_{\,i}$ que l'on met dans le
tableau de travail $\var{W}_4$. La valeur de ce terme sur la face interne $ij$
est obtenue en prenant la moyenne arithmtique des valeurs des deux cellules
avoisinantes (tableau \var{VISCF}) et celle sur la face de bord est prise gale
la valeur de la cellule avoisinante (tableau \var{VISCB}).
\item on calcule alors pour chaque direction $l$ le terme final,
\emph{i.e.} :\\
\begin{equation}\notag
\var{TRAV}(i,l)= \var{TRAV}(i,l)-\displaystyle \frac{2}{3}\left[\sum\limits_{j \in Vois(i)} (\mu_{\,tot}\,\dive(\underline{v}))_{\,moy,\,ij}\ S_{\,ij}^{\,l}
+\sum\limits_{k \in \gamma_b(i)} (\mu_{\,tot}
\dive\,(\underline{v}))_{\,moy,{b_{\,ik}}}\ S_{\,b_{ik}}^{\,l}\right]
\end{equation}
\end{itemize}
\vspace{0.5cm}
Le traitement est similaire pour le terme de viscosit de volume dans le module compressible.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Points \`a traiter} \label{Base_Visecv_paragraphe4}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L'intgration du terme de gradient transpos pose un problme de compatibilit.
En effet, le gradient transpos provient de l'criture de la divergence
du tenseur des contraintes (\emph{cf.} \fort{predvv}), soit :
\begin{equation}\notag
\dive\,(\tens{\sigma}) = \dive\,(- p \tens{Id} + \tens{\tau })
\end{equation}
o :
\begin{equation}\notag
\tens{\tau} = 2 \mu\ [\underbrace{\displaystyle \frac{1}{2}\,(\ggrad \underline
{v}\ +\ ^t\ggrad \underline {v})}_{\text{partie 1}} -\underbrace{\displaystyle \frac{2}{3}\,tr(\,\frac{1}{2}\,(\ggrad \underline
{v}\ +\ ^t\ggrad \underline {v}))\ \tens {Id}\ }_{\text{partie 2}}]
\end{equation}
Or, lorsque l'on intgre la premire partie de la divergence de
$\tens{\tau}$, on implicite le terme $\dive\,(\mu\,\ggrad \underline{v})$ et on
explicite le gradient transpos $\dive\,(\mu ^t\ggrad \underline{v})$.
Ce traitement fait intervenir la vitesse au centre des cellules. Elle ne
vrifie pas exactement la condition $\dive\,(\rho \underline{v})\,=\,0$. En
effet, au cours de l'tape de correction, on utilise un filtre Rhie et Chow
(\emph{cf.} \fort{resopv}) et la vitesse n'est mise jour qu' la fin de
l'tape. Par contre, lorsque l'on intgre la deuxime partie de
la divergence de $\tens{\tau}$ de faon explicite, on utilise
la vitesse issue du flux masse aux faces qui vrifie la condition $\dive\,(\rho
\underline{v})\,=\,0$ (du moins $\rho$ constant, l'interpolation de $\rho$
la face tant galement un point considrer). Ainsi, la discrtisation de ces
deux parties n'est pas totalement cohrente. Il serait utile de baser la
discrtisation de ces termes sur une vitesse vrifiant la contrainte
$\dive\,(\rho \underline{v})\,=\,0$.\\
Pour la mme raison, il est difficile de connatre les conditions aux limites du
terme en gradient transpos. Sur les cellules de bord, on sait uniquement que la
contrainte totale normale doit quilibrer le frottement et toutes les autres
forces. Or, le tenseur des contraintes est scind en une partie explicite et une partie implicite, donc c'est un peu difficile d'utiliser cette condition physique.\\
Actuellement, la contribution aux cellules de bord du terme de gradient transpos est
annule, ce qui limine l'influence des conditions aux limites mais n'est naturellement pas satisfaisant. Quelques essais d'intgration des conditions aux limites pour ce terme n'ont pas t concluants jusqu' prsent. Cependant, des essais supplmentaires
sont envisageables.
|
