Transcript Cours_9_Equations de Navier-Stokes3

Les algorithmes de découplage vitesse-pression
du problème de Stokes
Formulation du problème.


 
v
 v  p  f , dans   t  0
t
 
  v  0, dans     

 
v
v  V (ou
 ...),sur 
n
équivalent à :
2
 
 p   f

      
    v      f ,   v  0
 t

Discrétisation temporelle : formulation semi-discrète
P our une discrétisation temporelle
 n 1 dgt 1  n q
 0v    q v
q 0
 n 1
 n 1 
 v  p  f , dans   t  0
t
  n 1
  v  0, dans     


 n 1
v
n 1
v  V (ou
 ...),sur 
n
 n 1
f ne regroupeque des termessources connusau tempsn  1
ou traitésde façonexplicite.
dgt 1
 nq
 qv
 n 1 q 0
 n 1 

   v  p  f 
, dans 
0

t

H
S
On obtientle système:

 n 1 
H v  p  S , dans 
  n 1
  v  0, dans     


 n 1
v
n 1
v  V (ou
 ...),sur 
n
Nous allons chercher à découpler la vitesse de la pression en imposant
au mieux une vitesse à divergence nulle.
Les différentes familles de solveurs :
Découplage
Consistance
Coût de calcul
Div(v)=0
Uzawa
Oui
Cher
Exact
Green/Matrice d’influence Oui
Cher
Exact
Time-splitting
Non
Raisonnable
Approché
Projection/Diffusion
Oui
Raisonnable
Approché
Coût raisonnable : résolution tensorisée des problèmes de Poisson et
de Helmholtz .
Algorithme d’Uzawa
Le système d’équations est discrétisé en espace :
 



H v int  Dp int  Sint ,
 
D  v  0,
 
 
v  V (ou Dn  v  ...),sur les noeuds de frontière
Les composantes de vitesse et la pression sont exprimés en vecteurs colonne.
La première équation est écrite pour les nœuds intérieurs mais fait intervenir
les valeurs de vitesse et de pression en frontière.
H est l’opérateur de Helmholtz discret, écrit aux points intérieurs au domaine.
D est l’opérateur gradient (ou divergence) discret.
L’élimination des conditions de frontière conduit à un problème modifié :
"
"
"
"



Sx 
 H u uint   Dx 
 "

 " "
"
 H v vint    Dy  p   S x ,
 " " 
 H w   " D " 
S
 w int   z 
 x 
 
D  v  0.
Les opérateurs H sont les opérateurs de Helmholtz réduits.
Les matrices « D » sont les matrices rectangulaires Nint lignes et N colonnes
modifiées par l’élimination des conditions aux limites.
Les vecteurs « S » sont les second membres modifiés par l’inclusion des
conditions aux limites.
N: nombre de points total, Nint : nombre de points intérieurs
 Hu  Sx 
 uint   H u  Dx 




1 "
"
1 "
"
 vint    H v  D y  p   H v  S y ,
 1 " " 
 w   H 1" D " 
H

S
z 
w
z 
 int   w

1 "
"
1 "
"
 
 
 
D  v  Sdiv " D"  vint  0 " D"  vint   Sdiv
Les conditions aux limites sont utilisées pour déterminer Sdiv.


" D"  " Dx " " D y " " Dz "
"

Dx" ," Dy " ," Dz " , matricesN lignes N int colonnes
On en déduit l’équation à satisfaire pour la pression en prenant la divergence
de l’équation précédente :
 H u 1" Dx" 
 1

"
"
D
D
D
H

Dy  p   " Dx" " Dy " " Dz"

" x" "
y" " z"
v
 1 " " 
 H w  Dz 



 H u 1" S x" 
 1 " 
 H v " S y   Sdiv ,
 1 " " 
 H w  Sz 

Opérateur d’Uzawa, N lignes, N colonnes.
Opérateur à noyau : les modes parasites de pression sont présents.
En 3D, 100 points par direction, N=106…
Résolution itérative par préconditionnement.
Le calcul de la vitesse revient à résoudre des problèmes de Helmholtz
classiques pour chaque composante de vitesse.
Méthode de Green, Matrice d’influence.


 
 v
  v  p  f , dans   t  0
 t
  v  0,

2
 
par application de la divergence,  p    f

 
et donc      v  0.
 t

La divergence de la vitesse obéit à une équation de diffusion instationnaire.
Si elle est nulle à t=0 partout et si elle est nulle en frontière pour tout t>0,
alors elle reste nulle.
On se définit :
N b  N  N int est le nombrede pointsde bord
xi  ( xi , yi , zi ), i  1,...,Nb ces pointsde bord.
On se définit de plus :

 
N b champsde pressionqi(x) et champsde vitessevi(x) vérifiant:
2 
 qi ( x )  0, qi ( x j )   ij
 2   


 vi ( x )  qi ( x )  0, vi ( x j )  0.
Ces champs de vitesse n’étant pas à divergence nulle, on définit :

 
DIVji    vi

xj
Ces champs sont calculés une fois pour toutes et stockés.
A chaque pas de temps :
2
 

 p0 ( x )    f , p0 ( x j )  0,
 

   2   
   v0 ( x )  p0 ( x )  f , v0 ( x j )  0
 t

On calcule :

 
Q j    v0

xj
,    DIV 1Q
La solution à divergence nulle aux points de bord est alors :
Nb
  Nb 
p  p0  i qi , v  v0  i vi .
i 1
i 1
Cette méthode est coûteuse :
Il faut calculer un grand nombre de solutions homogènes des équations (autant
que de points de bords).
Le calcul de l’inverse de la matrice DIV peut devenir coûteux en 3D
Méthode de splitting temporel
 n 1 dgt 1  n q
 0v    q v
q 0
 n 1
 n 1 
 v  p  f , dans   t  0
t
  n 1
  v  0, dans     
est remplacépar :
dgt 1
 n q
vˆ   q v
q 0
t
 n 1 

 p  f ,   vˆ  0,
 n 1
 n 1
 n 1
 0v  vˆ
 v , associé aux conditionsaux limitessur v .
t
Le calcul de la pression est découplé du calcul de la vitesse au temps n+1.
Il est effectué en prenant la divergence de la première équation, et en imposant
que le champ intermédiaire soit à divergence nulle (étape dite de projection).
Cela amène donc à une équation de Poisson pour la pression (nous discuterons
des conditions aux limites à imposer par la suite).
Le champ intermédiaire est donc imposé à divergence nulle, quelle que soit
la solénoïdalité des champs précédents. Aucune condition aux limites n’est
imposée au champ intermédiaire qui est simplement évalué dans l’intérieur du
domaine.
Il intervient alors comme terme source de la vitesse au temps n+1, lui aussi
sinusoïdal, via la seconde équation. La vitesse au temps n+1 est solution d’un
problème de Helmholtz.
L’étape de projection :
 dgt 1  n q

  q v

 n 1 
  q 0
p   
f
, dans 


t




 
  n 1  n 1   n 1 
n.p  n. v  f
 vb .n , sur .
t


On utilise l’équation de Navier-Stokes projetée en frontière pour assurer
les conditions aux limites, de Neumann partout.
La valeur de la pression doit donc être imposée en frontière. Le mode constant
est le seul mode parasite de pression non filtré par la méthode. Div(V)=0 n’est
pas explicitement imposé en frontière et l’information n’est pas redondante
avec l’imposition des conditions aux limites.
Dans la condition aux limites, la vitesse au temps n+1 doit être extrapolée :


   n 1
 n 1    n 1    n 1
v   v
     v      v
dgt 1
   nq
   q     v
 o t dgt 


q 0
Le coût de la méthode est minimal : un opérateur de poisson pour la
pression et un opérateur de Helmholtz pour la vitesse.
Consistance du schéma :
en posant f=0, et en laissant de côté les conditions aux limites, on peut
considérer le découplage du système comme un étalement sur deux
équations du terme de dérivation temporelle.


 

vˆ
v
v vˆ

 p,   vˆ  0, (1   )   v .
t
t
t t
Découpler les équations introduit une équation d’évolution temporelle
sur la pression :



(
1


)



p  0
t


qui viole l’hypothèse d’incompressibilité.
Méthode de projection-Diffusion
Contrairement au schéma de découplage précédent, celui-ci s’introduit sans
nécessiter de discrétisation temporelle.



 a
 p  f
   
(1)   a 
Projection
 0
 
  V
 


 v   n 
a  n  
CL extrapolée en temps

t




 

 
(2)   v  a , associé aux conditionsaux limit essur la vitesse.
 t

Diffusion
On sépare les termes de l’équation de divergence nulle
 
 


.a  0  int .aint  bord .abord  0

int est la restriction à l' intérieurdu domainede l' opérateurdivergence.

bord est le restede l' opérateurs' appliquantaux pointsde bord.

a bord est imposécommeconditionaux limites.



aint   p  f


int


 




int .aint  int .  p  f int  bord .abord
 

 

int . p int  bord .abord  int . f int
 
On obtient un opérateur « quasi Poisson » pour la pression en remplaçant
les composantes de l’accélération à l’intérieur du domaine par les
composantes du gradient de pression.
   
" " P  B  aB    f
Aucune condition aux limites sur la pression n’est imposée !
L’étape de diffusion est classique, mais il faut noter que par troncature,
le second membre ne prendra en compte qu’en partie de l’accélération
calculée. La divergence de la vitesse ne sera donc pas nulle mais
convergera exponentiellement vers 0 avec la discrétisation spatiale.
Problème modèle 1D de l’étape de projection :
dp

a  dz  f ( z ), z   1,1

a ( z  1)  a

 da
 dz  g ( z ), z   1,1
Problème modèle : g(z)=0 impliquerait une solution triviale.
1
Il faut respecter

1
da
dz  a  a
1 dz
1
g ( z )dz  
Ecrire le système d’équations discret et en déduire le système à
résoudre pour déterminer la pression.
 a1 
DZ a   g   DZ  f  p1   g


 aN 1 
0
 a1 
DZ  p1   DZ  f   g  DZ
 0 
aN 1 


 p1 
 a1 
DZ  p   DZ  f   g




 pN 1 
aN 1 

Ligne ou colonne manquante conduisant à des opérateurs rectangulaires
L’opérateur possède des modes parasites de pression : les filtrer en
calculant son inverse.
Pour déterminer un cas test :
- Choisir une fonction aexacte
- En déduire g par dérivation de aexacte
-Choisir une fonction pour la pression pexacte
- En déduire f par dérivation de pexacte
Rédiger le programme et imposer f et g.
Comparer les champs a et p calculés à leur valeurs exactes.
Programme E_9_toy.m
Ce qui n’a pas été abordé
En Chebyshev :
- la méthode tau
- les procédures de décomposition de domaine (domaines en L par exemple)
- la résolution par préconditionnement d’opérateurs non tensorisables (passage
en géométries non orthogonale, physique complexe …)
En Legendre :
L’approche éléments spectraux (éléments finis d’ordre élevé), donc
en formulation Galerkin
Intérêt :
décomposition de domaines assurée directement par la méthode
gestion naturelle de géométries déformées