Transcript Cours_1

Simulation en Dynamique
des Fluides
M2 SDFT, Université Paris-Sud
G. Kasperski, [email protected]
C.T. Pham, [email protected]
Organisation de l’UE
• Principes généraux de simulation (G.
Kasperski, 6h)
• Méthodes des différences finies et analyse
d’opérateurs discrets (C.T. Pham, 18h,
mercredi après-midi)
• Méthodes spectrales, résolution des
équations de Navier-Stokes (G. Kasperski,
18h, vendredi matin)
Ma partie…
• Généralités sur la simulation numérique
• Matlab pour les plus nuls que moi (s’il y en a)
• Méthodes des résidus pondérés (1ère
version…)
• Cours de méthodes spectrales
– bases théoriques
– nombreuses applications Matlab
– petit projet sous Matlab à réaliser
Pourquoi des simulations
numériques ?
- Investigation du comportement de systèmes complexes
(absence de solution analytique).
- Mise en œuvre et test de modèles théoriques pour
confrontation aux expériences.
- Cas inaccessibles aux méthodes expérimentales non
intrusives.
- Accès à des données expérimentales difficilement
mesurables.
-Souplesse phase R&D d’un projet industriel, coût moindre que
la réalisation de prototypes.
Navier-Stokes
2D
Convection
thermique
Écoulement réel.
Modélisation physique : approximation de Boussinesq, conditions aux limites
choisies …
Modèle mathématique : équations de Navier-Stokes (ici incompressibles).
Modèle numérique discret décrivant le problème comme de
grands systèmes linéaires à résoudre successivement.
Résolution des systèmes linéaires : Utilisation de Matlab.
Solution numérique.
Pourquoi s’en méfier ?
•
Sources d’erreurs inhérentes aux schémas, dispersion
et diffusion artificielles pour les méthodes aux noms se
terminant par « finies »…
•
Seule la physique modélisée est observée.
• Coût en temps de calcul augmentant très rapidement
avec la précision souhaitée pour la solution: des
écoulements turbulents 3D sont encore non
envisageables en DNS dans des géométries un peu
complexes.
Pourquoi des méthodes spectrales ?
Alors qu’elles sont
- un peu compliquées à appréhender
- souvent un peu lourdes à mettre en œuvre
- limitées à des géométries simples
Computational predictability of time-dependent natural convection flows in enclosures (including a benchmark
solution) Mark A. Christon, Philip M. Gresho and Steven B. Sutton, Int. J. Numer. Meth. Fluids 2002; 40:953–
980
Les méthodes de simulation numérique
ou, de façon (presque) équivalente,
Les méthodes des résidus pondérés
- Discrétisation : représentation des solutions sous forme discrète, comme
valeurs sur un maillage ou par décomposition sur une base de fonctions.
- Résidus : ils qualifient l’erreur commise en appliquant l’équation différentielle
à résoudre à une solution discrétisée… nous adopterons l’approche « base de
fonctions ».
(N.B. il existe beaucoup de définitions d’erreurs différentes)
- Pondérés : qualifie la démarche utilisée pour annuler ce résidu.
Problème différentiel sur un domaine Ω de frontière G:
Nous souhaitonsapprocherue , solutionexactedu problème
 L(ue )  g  0 dans 

 B(ue )  h  0 sur G
équationdifférentielle

conditionaux limites
 
G

n
exemple :
  

 v .    T  g
équation de la chaleur L(T )  
 t
avec des conditions aux limites : B (T )  T  
A noter :

T
h
n
 0
conditionde Dirichlet
 0
conditionde Neumann
 et   0 conditionde Robin
Discrétisation Spatiale
On cherche à approcher la solution exacte ue par une solution numérique
U:
N


U ( x , t )   ci (t )i ( x )
i 0
Les i sont appelées fonctions de forme, fonctions modales ou fonctions
test. Elles doivent être linéairement indépendantes les unes des autres
(définition unique de U).
La qualité de l’approximation dépendra du choix de ces fonctions : si on
souhaite représenter exactement une solution dans un espace de
fonctions donné, les i doivent en former une base.
La solution est représentée en fonction du temps par les N+1 valeurs ci(t)
(nous traiterons la discrétisation temporelle par la suite).
On définit le résidu du problèmedifférentiel par

R( x , t )  L(U )  g dans ,

R (x,t)  B(U )  h sur G.
Ces résidus sont non nuls a priori, car les équations ne sont exactement
vérifiées que si U=ue. Ils mesurent l’erreur locale commise sur la
satisfaction de l’équation différentielle.
Le but de la méthode numérique va être de chercher à rendre ce
résidu le plus petit possible.
Cela se fera en déterminant les N+1 coefficients ci de la discrétisation
pour un problème stationnaire.
Aparté matheux :
Espace fonctionnel : structure d’espace vectoriel.
Soient deux fonctions f et g définies sur un domaine D, l’intégrale


 
( f g )   f ( x ). g ( x ). w( x )dx
D
avec w une fonction strictement positive définit un produit scalaire des
fonctions f et g.
forme bilinéaire en f et g
symétrique en f et g
positive
définie
f g   f g   f g   f g
f g  g f
f f 0
f f 0 f 0
Nous l’utiliserons ici avec w=1.
On cherche à imposer l’annulation du résidu…
d’un certain point de vue :
en le rendant orthogonal à un ensemble
de fonctions choisies (P,P) .
G
 
n
P définie dans Ω
P définie sur Γ
N+1 inconnues pour la solution
N+1 équations à écrire
N+1 couples de fonctions (Pj,Pj) à choisir.




R
(
x
,
t
)
P
(
x
)
d


R
(
x
,
t
)
P
(
x
)dG  0, pour j  0, N
j
j



G
Pj et Pj sont appeléesfonctionsde pondération.
Le traitement des conditions aux limites peut faire l’objet de différentes
techniques via le choix des couples de fonctions de pondération.
Méthodes de collocation :
les fonctions de pondération sont des fonctions Dirac définies à partir
d’un jeu de points définis dans le domaine et sur le bord.

 

Pj ( x )   ( x  x j ) et Pj ( x )  0 pour les pointsdu domaine.


 
Pj ( x )  0 et Pj ( x )   ( x  x j ) pour les pointsde bord.
Les équations deviennent simplement :
G


R( x j )  0 si x j  


R ( x j )  0 si x j  G

Méthodes de sous-domaines
On définit un ensemble de N+1 sous-domaines Dj de .
Domaine intérieur



Pj ( x )  1 si x  Dj , Pj ( x )  0 sinon.
Implémentation spéciale des conditions aux limites : en général,
-Collocation pour les conditions de Dirichlet
-Intégration dans l’équation obtenue sur un sous-domaine de bord pour Neumann.
L’équation différentielle est intégrée analytiquement.
Le degré maximal de dérivation apparaissant
dans le problème est diminué d’une unité :
il s’agit d’une formulation faible du problème.
Dj
Intérêt : les fonctions j utilisées peuvent être
de classe de continuité inférieure à celle nécessaire
à la résolution de l’équation différentielle.
Une base de fonctions linéaires par morceaux peut être utilisée pour approcher
la solution d’un problème faisant intervenir des dérivées secondes (diffusion).
Méthode des moindres carrés
On utilise la base j de la décomposition modale.
R
Pj 
c j


 R
1 
2 
R
(
x
,
t
)
P
(
x
)
d


R
(
x
,
t
)
d


R
( x, t )d  0
j



c j
2 c j 
On cherche ainsi un minimum du carré de la norme L2 du résidu.
La méthode est en général complexe à mettre en œuvre et est
relativement peu utilisée.
1/ p
Aparté matheux :
p
NormeL d' une fonctionf : f
p


p

   f d 


Par passage à la limite, on définit la norme infinie de la fonction
comme le maximum de la valeur absolue de f sur le domaine.
Méthode de Galerkin
On utilise encore la base j de la décomposition modale.
Pj   j



 R( x, t ) Pj ( x )d   R( x, t ) j d  0


La méthode impose au résidu d’être orthogonal à l’espace de définition
de la solution.
Méthode de Galerkin en formulation faible
Si les fonctions j sont dérivables, il est possible d’intégrer par partie
(dans le cas 1D) ou d’utiliser une formule de Green (cas multidimensionnel)
pour diminuer le degré de dérivation maximal intervenant dans le problème
différentiel. On obtient alors une formulation faible du problème.
 

Formule de Green :  V .dG   div (V )d
G


dG est le vecteurd' élémentde contourdu domaine, normalau bord de 
et dirigé versl' extérieur.

Pour deux fonctionsu et v et un vecteur x unitairede directionarbitraire,
on en déduit :
 
 
 uv
 v 
 u 
d


uv
x
.
d
G

u
d


uv
x
.
d
G

v
 x
G
 x
G
 x d
On obtient la formulation faible en utilisant la fonction de pondération à la
place de u et le terme de plus haut degré de dérivation du résidu à la place de v.
Les grandes familles de méthodes numériques pour les EDP
Collocation
Sous-domaines
Galerkin
Moindres carrés
formulation
forte
formulation
faible
formulation
forte ou faible
(plutôt faible)
formulation
forte
- Idée des
différences
finies (mais
pas de
décomposition
modale)
-Collocation
spectrale
Chebyshev
-Volumes finis
- Eléments finis
-Volumes
spectraux
- Eléments
spectraux
Legendre
- vue une fois
dans une
publication.
…..
- Chebyshev
Galerkin
- Chebyshev Tau
…… On aura une précision spectrale ou…..
de précision finie selon le choix de la base j
Exercice illustratif
On veut approcher la solution d’une équation différentielle stationnaire
 2ue
 1 pour 0  x  1 avec ue( 0 )  0 et ue (1)  1.
2
x
N
U   c j j ( x )
j 0
- Poser les conditions de satisfaction des conditions aux limites (parti
pris)
- Ecrire le système obtenu en collocation utilisant des points équidistants
- Ecrire le système obtenu en méthode de sous-domaines sur les
domaines définis par les points équidistants précédents
- Ecrire le système obtenu par une méthode de Galerkin
- Ecrire le système obtenu par méthode des moindres carrés.
Un peu d’interpolation :
de l’intérêt d’une bonne base
de fonctions
Polynômes de Lagrange et théorème de l’unicité :
Un support d’interpolation est un ensemble de valeurs fp d’une fonction f définies
aux points particuliers x=xp.
Sur un support à n+1 points, il existe un seul polynôme de degré N passant par tous les
points du support: le polynôme de Lagrange.
N
Il s' écrit PN ( x)   a j x j et les coefficients sont déterminéspar :
j 0
Calcul: - résolution du système linéaire suivant:
P( x0 )  f 0
 1 x0

P ( x1 )  f1
1 x1


... ...
...

P( xN )  f N
 1 xN
- formule de Lagrange:
x 02
x12
...
xN2
... x 0N   a0   f 0 
   
... x1N   a1   f1 

... ...   ...   ... 
  
N 
... xN  a N   f N 
N
N
i 0
k 0
k i
PN ( x )   fi li ( x ), avec li ( x )  
x  xk
xi  xk
En effet : li ( x j )   ij
25
L’écriture polynomiale d’un polynôme de degré N est unique à partir de la
connaissance de ses valeurs en N+1 points distincts.
Si la fonction à interpoler est un polynôme de degré N, l’erreur d’approximation
est nulle mathématiquement, a priori très faible numériquement.
Vérifions-le avec un programme Matlab utilisant une fonction d’interpolation
de Lagrange :
function [inter] =lagrange(fxp,xp,x,n)
La fonction renvoie
la valeur de l’interpolé en
x
degré du polymôme
d’interpolation
point où on cherche
la valeur de l’interpolation
valeurs de la fonction
aux points du support
support d’interpolation
Dans le programme fourni : P1_test_Lagrange :
On définit une fonction fxp le support d’interpolation régulier de
x p (i)  1  (i  1)dx, dx  2 / N , i  1,...,N  1
x p (i) [1,1]
On calcule l’interpolation de f sur tous les points d’un maillage dix fois plus
raffiné.
x pextra (i)  1  (i  1).dx' , dx'  2 /(10 N ), i  0,...,10 N  1
On compare les valeurs de l’interpolé à la valeur exacte de f, pour différentes
valeurs de N et on trace l’erreur sur les points xpextra).
support régulier
L’écriture polynomiale d’un polynôme de degré N est unique à partir de la
connaissance de ses valeurs en N+1 points distincts.
Si la fonction à interpoler est un polynôme de degré N, l’erreur
d’approximation est nulle mathématiquement, a priori très faible
numériquement.
Vérifions-le avec la première une fonction
f ( x)  x n
- Voyons ce qui se passe pour une fonction non polynômiale f ( x ) 
- Changer de fonction, essayer avec f ( x)  cos( x)
- Revenir à f ( x ) 
1
mais changer le support d’interpolation,
1  25 x 2
 2i - 1 
x p (i)  -cos
, i  1,...,N  1
 2( N  1) 
1
1  25 x 2
Interpolation sur un support homogène : équivaut à une approximation
Fourier (ne fonctionne que sur des fonctions périodiques).
Interpolation sur le second support : équivaut à une approximation
Chebyshev (fonctionne toujours).
Intérêt : - approcher au mieux une solution
- annuler au mieux un résidu
La discrétisation temporelle
Une fois la méthode des résidus pondérés appliquée au problème, restent
à traiter les termes de dérivation temporelle. On ne traitera que les équations
d’évolution où ce terme apparaît seul et au premier ordre (cas des équations
de Navier-Stokes).
c
 f ( c, t ,...) où c représente le
Nous avons un problème du type
t
vecteur défini par les coefficients de la discrétisation.
Contrairement au domaine spatial, le domaine temporel à traiter n’est pas
défini à l’avance. Recherche de solution stationnaire, attente d’établissement
d’instabilité, …
Pour la discrétisation temporelle, on utilisera toujours une approximation
locale de la dérivée temporelle, faisant intervenir le vecteur solution au temps
du calcul et à quelques instants précédents.
La solutionsera supposéeconnueen des instantsdiscrets tn  nt,
n étantun entieret t le pas de tempsde la simulation.
On obtient alors la solution par approximation locale de la dérivée temporelle.
n
2 n

c n 1  c n 1  n
c n t 2  2c n

c

t

c
n

  c  t


...

c


 ...
2
2

t
t 
t
2 t
2 t
 t
La discrétisation équivaut à la dérivée temporelle moyennant une erreur
temporelle d’ordre 1 (puissance à laquelle intervient le pas de temps
dans l’erreur).
De même, cette dérivée temporelle peut être écrite
c n 1  c n 1  n 1 n 1
c n 1 t 2  2c n 1  c n 1 t  2c n 1
  c  c  t

... 

 ...
2
2
t
t 
t
2 t
t
2 t

Il suffit d’exprimer le second membre f au même temps que la dérivée
temporelle (ou approché au même temps avec le même ordre de précision)
pour obtenir un schéma complet.
33
Quelques schémas classiques :
Adams Bashforth :
Adams Moulton :
c n 1  c n
 fn
t
c n 1  c n 3 n 1 n 1
 f  f
t
2
2
c n 1  c n 23 n 16 n 1 5 n 2

f 
f 
f
t
12
12
12
c n 1  c n
 f n 1
t
c n 1  c n 1 n 1 n
 f  f 
t
2
c n 1  c n
5 n 1 8 n 1 n 1

f 
f 
f
t
12
12
12
Euler retardé d’ordre 2 :
3c
n 1
 4c  c
2t
n
n 1
 f n 1
Saute-moutons (leapfrog, ordre 2):
3c n 1  4c n  c n 1
 2 f n  f n 1
2t
c n 1  c n 1
 fn
2t
Différence entre les schémas de différentes couleurs ?
Pour un schéma explicite,
c n’intervient au temps n+1 que dans le terme de dérivée temporelle.
Le système linéaire résultant est plus simple à résoudre : éventuellement
diagonal.
Par contre le schéma devient conditionnellement stable (voir cours
C.T. Pham) : cela induit des restrictions sur les pas de temps utilisables.
Pour un schéma implicite,
c intervient au temps n+1 dans le terme de dérivée temporelle ainsi que
dans l’expression de f. Le système linéaire est plus complexe à résoudre.
Par contre, il est beaucoup plus stable (souvent inconditionnellement stable).
Exercice : vérifier que
3c n 1  4c n  c n 1
2t
fournit la dérivée temporelle de c à l’ordre 2 au temps n+1.
Traitement des termes non linéaires.
En règle générale, les termes non-linéaires d’une équation d’évolution
temporelle sont traités de façon explicite.
Ils sont donc simplement évalués aux instants précédents et n’entrent pas
dans l’opérateur à résoudre numériquement.
Des procédures de recherche de solution stationnaires de problèmes non
linéaires existent mais ne seront pas traités dans ce cours (voir les
méthodes de Newton-Raphson).
Schéma d’Euler
retardé explicite :
 f 
f n  f n 1  t. 
 t 
f
Approximation :
n 1
 f
 f
n 1
n 1

 f 
 2t. 
 t 
2f  f
Schéma
n 1
n
n 1
t   f 
n 1
2
2
n 1
 2 
2  t 
 O (t 3 )
 f
 2t  2
 t
n 1
2
2
 f
 t  2
 t
2
2






 O (t 3 )
n 1
 O (t 3 )
Erreur de l’approximation
Problème de diffusion instationnaire :
 ue  2ue
 t  x 2  S ( x, t )

ue ( x  0, t )  0
u ( x  1, t )  1
 e

Fonctions de base
1
maillage homogène
xi-1
x1=0
xi
xi+1
Première fonction
Fonction générique i
xN+1=1
x
Dernière fonction
La fonction sera linéaire par morceaux :
Ui
Ui+1
- Les coefficients de la
décomposition sont les valeurs
du champ aux points de
maillage.
-Méthode de sous-domaines
définis par les points zi à midistance des points de
maillage
Ui-1
zi 1

xi-1
zi
xi
x
zi+1
x
xi+1
zi
U
dx 
t
zi 1

zi
 2U
dx 
2
x
 i 1
U
Udx

t zi
x
z
zi 1
U

x
zi 1
 S ( x, t )dx
zi

zi
zi 1
 S ( x, t )dx
zi
Montrer que cela conduit à :
  Ui 1 3Ui Ui 1 
Ui 1  Ui Ui  Ui 1  Si 1 3Si Si 1 







x 
x
t  4
2
4 
x
x
2
4 
 4
  Ui 1 3Ui Ui 1  Ui 1  Ui Ui  Ui 1  Si 1 3Si Si 1 









2
2
t  4
2
4 
x
x
2
4 
 4
c
 f ( c, t ,...)
t
Adams-Bashforth
Application des deux schémas temporels
c n 1  c n
 fn
t
Uin11 3Uin1 Uin11 Uin1 3Uin Uin1 Uin1  Uin Uin  Uin1  Sin1 3Sin Sin1 








 


2
2
4
2
4
4
2
4
x
x
2
4 
 4
c n 1  c n
 f n 1
Adams-Moulton
t
Uin11 3Uin1 Uin11 Uin11  Uin1 Uin1  Uin11 Uin1 3Uin Uin1  Sin11 3Sin1 Sin11 








 


2
2
4
2
4
x
x
4
2
4  4
2
4 
… Ecrire les systèmes linéaires grâce à ces deux équations associées aux
conditions aux limites : AUn+1=b
Sous Matlab : Programmer les deux opérateurs.
Programme de résolution :
Initialisations
Définir un maillage homogène à N+1 points et le pas d’espace x.
Définir un pas de temps.
Définir deux vecteurs initiaux U0 et U1
Résolutions successives : pour les pas de temps successifs (boucle for)
Evaluer le temps du calcul ->définit éventuellement la valeur du terme source
Calculer b
Calculer U2=A-1b
Placer U1 dans U0 puis U2 dans U1 -> on ne stockera pas tous les champs en mémoire
Tracer la solution au fur et à mesure des calculs
Continuer les itérations
On utilisera les deux opérateurs dans deux programmes différents.
Démarche de validation du programme
On se définit une solution exacte (éventuellement satisfaisant les conditions aux limites).
ue  x  sin(x) sin(t )
On en déduit le terme source à imposer en fonction du temps pour l’obtenir :
ue  2ue
S ( x, t ) 
 2   sin(x ) cos(t )   2 sin(x ) sin(t )
t
x
On adapte les conditions aux limites à la solution exacte (ici pas nécessaire)
On initialise les champs initiaux avec les champs exacts :
t  0 :U 0  x
t  t : U 1  x  sin(x ) sin(t )
En cours de programme, on compare la solution calculée à la solution
exacte connue et sensée être obtenue.
Si on améliore (comme on l’attend) le résultat en diminuant les pas de
temps et d’espace, le programme est validé.