Simulation en Dynamique des Fluides

1. Simulation en Dynamique des Fluides M2 SDFT, Université Paris-Sud G. Kasperski, kasperski@fast.u-psud.fr C.T. Pham, Chi-Tuong.Pham@limsi.fr

2. 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)

3. 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

4. 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.

5. 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.

6. 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.

7. 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

8. 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

10. Les méthodes de simulation numérique ou, de façon (presque) équivalente, Les méthodes des résidus pondérés

11. - 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.

12. G W Problème différentiel sur un domaine Ω de frontière G: exemple : équation de la chaleur avec des conditions aux limites : A noter :

13. Discrétisation Spatiale On cherche à approcher la solution exacteue par une solution numérique U : Les fi 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 fi 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).

14. 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.

15. Aparté matheux : Espace fonctionnel : structure d’espace vectoriel. Soient deux fonctions f et g définies sur un domaine D, l’intégrale 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 Nous l’utiliserons ici avec w=1.

16. N+1 inconnues pour la solution N+1 équations à écrire N+1 couples de fonctions (Pj,Pj) à choisir. G W 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) . Le traitement des conditions aux limites peut faire l’objet de différentes techniques via le choix des couples de fonctions de pondération.

17. 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. Les équations deviennent simplement : G W

18. Dj Méthodes de sous-domaines On définit un ensemble de N+1 sous-domaines Djde W. Domaine intérieur • 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. Intérêt : les fonctions fj 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).

19. Méthode des moindres carrés On utilise la base fj de la décomposition modale. On utilise la base fj de la décomposition modale. 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. Aparté matheux : 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.

20. Méthode de Galerkin On utilise encore la base fj de la décomposition modale. La méthode impose au résidu d’être orthogonal à l’espace de définition de la solution.

21. Méthode de Galerkin en formulation faible Si les fonctions fjsont 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 : 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.

22. Les grandes familles de méthodes numériques pour les EDP Sous-domaines Moindres carrés Collocation Galerkin formulation forte formulation forte formulation faible formulation forte ou faible (plutôt faible) • Volumes finis • Volumes • spectraux • ….. - vue une fois dans une publication. • - Eléments finis • Eléments • spectraux Legendre • Chebyshev Galerkin • Chebyshev Tau • ….. • - Idée des différences • finies (mais pas de • décomposition modale) • Collocation • spectrale Chebyshev • …… On aura une précision spectrale ou de précision finie selon le choix de la base fj

23. Exercice illustratif On veut approcher la solution d’une équation différentielle stationnaire • 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.

24. Un peu d’interpolation :de l’intérêt d’une bonne basede fonctions

25. 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. Calcul: - résolution du système linéaire suivant: - formule de Lagrange: 25

26. 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

27. Dans le programme fourni : P1_test_Lagrange : On définit une fonction fxp le support d’interpolation régulier de On calcule l’interpolation de f sur tous les points d’un maillage dix fois plus raffiné. 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

28. 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

29. - Voyons ce qui se passe pour une fonction non polynômiale - Changer de fonction, essayer avec - Revenir à mais changer le support d’interpolation,

30. 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

31. La discrétisation temporelle

32. 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). où c représente le Nous avons un problème du type 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.

33. On obtient alors la solution par approximation locale de la dérivée temporelle. 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 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

34. Quelques schémas classiques : Adams Bashforth : Adams Moulton : Euler retardé d’ordre 2 : Saute-moutons (leapfrog, ordre 2): Différence entre les schémas de différentes couleurs ?

35. 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).

36. Exercice : vérifier que fournit la dérivée temporelle de c à l’ordre 2 au temps n+1.

37. 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 : Approximation : Schéma Erreur de l’approximation

38. Problème de diffusion instationnaire : Fonctions de base 1 maillage homogène xi+1 xi-1 xi xN+1=1 x x1=0 Dernière fonction Première fonction Fonction générique i

39. - 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 à mi-distance des points de maillage La fonction sera linéaire par morceaux : Ui Ui+1 Ui-1 xi-1 xi+1 xi zi zi+1 Montrer que cela conduit à :

40. Application des deux schémas temporels Adams-Bashforth Adams-Moulton … Ecrire les systèmes linéaires grâce à ces deux équations associées aux conditions aux limites : AUn+1=b

41. 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 dx. 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.

42. Démarche de validation du programme On se définit une solution exacte (éventuellement satisfaisant les conditions aux limites). On en déduit le terme source à imposer en fonction du temps pour l’obtenir : On adapte les conditions aux limites à la solution exacte (ici pas nécessaire) On initialise les champs initiaux avec les champs exacts : 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é.