Calcul d’Écoulements

, par MiKaël Navarro

Le but de ce tp est de calculer un écoulement potentiel par différences finies autour d’un cockpit d’avion.

Construction du Maillage

Solveur elliptique

On souhaite obtenir un maillage quadrangulaire sur un domaine topologiquement semblable à un carré. Pour cela on associe une "carte" aussi régulière que possible entre le carré unité et le domaine.

On associe donc un couple de coordonnées curvilignes (’epsilon’ et ’nu’ du domaine) à un couple de coordonnées cartésiennes (x et y), avec des conditions de régularité et de monotonie.

On peut remplir ces deux solutions en imposant :

En exprimant, cette relation par le biais des fonctions inverses :

on arrive à : Px = Py = 0

où P est l’opérateur différentiel suivant :

Ces relations sont du second ordre, non-linéaires et elliptiques.

Discrétisation

On dicrétise le problème sur le carré unité : epsilon=j/jM , eta=k/kM , j = 0,2,...,jM et k = 0,2,...,kM.

L’équation Pu=0 peut s’approcher de la manière suivante (S) :

k=1,2,...,kM-1 et j=1,2,...,jM-1
et où les coefficients alphaj,k, betaj,k, gammaj,k sont :

On en déduit donc, l’algorithme de résolution :

  • Initialisation de xj,k et yj,k ;
  • Itération (sur les points intérieurs)
    • Linéarisation : calculer les coefficients alphaj,k, betaj,k, gammaj,k par différences finies
    • Gauss Seidel : appliquer un certain nombre d’itérations du système (S) pour xj,k et yj,k.

Résultats

Aspect du maillage :

On s’apercoit que la convergence est vraiment très rapide : visuellement, on ne voit plus de différence sur le maillage après la 10ième linéarisation avec seulement 5 itérations du système (S) par linéarisation.

ceci est le maillage initial, taille = 64*24.

Après passage par solveur elliptique :

Cette image présente le résultat du solveur elliptique sur le maillage initial.
Ce maillage a été obtenu après seulement 9 linéarisations, pour 5 itérations de Gauss Seidel par linéarisation.
Visuellement, le maillage semble déjà suffisement régulier.

Cette technique permet donc d’obtenir un maillage régulier avec un temps de calcul très court !

Etudes des résidus

Pour tracer les courbes de résidus, nous avons effectué 25 étapes de linéarisations, avec 200 itérations de Gauss Seidel du système (S) par linéarisation. Nous présentons d’une part le résidu "non linéaire", c’est à dire le résidu calculé après l’étape de linéarisation, et d’autre part le résidu lors d’un cycle d’itération (le 10ième) de Gauss Seidel (résidu "linéaire").

Cela nous permet donc de distinquer la qualité de la linéarisation, de celle de la résolution par itérations de Gauss Seidel.

Résidu non linéaire :

log10 (Résidu / Résiduinitial) = f(nombre linéarisation)
Cette image montre comment évolue le résidu après chaque linéarisation
Globalement, le résidu se comporte de manière linéaire.

Résidu linéaire :

log10 (Résidu / Résiduinitial) = f(nombre itération)

Cette image montre comment évolue le résidu de l’équation (S), lors du 10ième cycle d’itération de Gauss Seidel.

Le residu se comporte de facon quasi-lineaire apres 40 iterations.

Pour profiter de la pente plus importante lors des premieres iterations il serait peut-etre judicieux de faire une linearisation des que le residu ne decroit plus suffisamment vite.

Calcul de l’écoulement

Ecoulement potentiel

On simule un écoulement plan et irrotationnel. Dans ce cas, on trouve que les lignes de courants sont les lignes de niveau d’une fonction Psi(x,y) qui vérifie Laplacien(Psi)=0, on a donc encore : P Psi = 0

Les conditions aux limites sont les suivantes :

 Sur le cockpit : Psi = constante = 0 ; condition de type Dirichlet ;
 sur le bord : d(Psi)/dn = - Vinf ny ; condition de type Neumann (n est la normale extérieure unitaire).

Discrétisation

On discrétise P Psi = 0 à l’intérieur du maillage par le même systeme (S)

La condition de Neumann se discrétise en :

PsiA = betaA PsiB + (1-betaA) PsiC - gammaA.

où betaA et gammaA s’écrivent en fonction des cordonnées des points A,B,C,D considérés.

Résultats

Aspect de l’écoulement :

La convergence est beaucoup moins rapide que pour le maillage.
Voici quelques résultats que nous avons visualisé avec Maple.

5 linéarisations , 2000 itérations par linéarisation.

15 linéarisations , 2000 itérations par linéarisation.

25 linéarisations , 2000 itérations par linéarisation.

Etudes des résidus

Le calcul de l’ecoulement evolue de facon lineaire, les perturbations en fin de courbe sont peut-etre dus a des erreurs d’arondi.

Téléchargement

Pour pouvoir faire tourner ce code, vous devez :

 Avoir le programme GnuPlot (gratuit), pour la visualisation des maillages et des résidus ;
 Avoir le programme Maple, pour la visualisation de l’écoulement ;
 Télécharger et décompresser l’archive qui contient les sources ;
 Compiler les sources (un makefile est fourni) ; (nb : sous windows, remplacer l’appel système à rm par un del).

Note :
GnuPlot doit être dans le même répertoire ou dans le PATH.
Il existe une version de Gnuplot sous windows.

Auteurs

ESSI3/CSI
 Cauvin Florian
 Navarro Mikael