Mathématiques avec Python et Ruby/Joukovski et Ruby
L'écoulement d'un fluide autour d'un cylindre est aisé à calculer (par exemple pour étudier l'effet Magnus), et toute transformation conforme permet d'en déduire l'écoulement autour du transformé du cercle (la trace du cylindre). En particulier lorsque le transformé a la forme d'une aile d'avion, ce qui est le cas avec la transformée de Joukovski. Ce calcul sera effectué par Ruby grâce à ses nombres complexes, puis affiché en fabriquant une figure svg. En fait:
- On utilise l'inverse de la transformée de Joukovski pour calculer l'écoulement autour du cercle unité.
- On agrandit légèrement ce cercle pour qu'il ait une transformée de Joukovski plus intéressante;
- On applique au nouveau cercle la transformée de Joukovski: On a l'écoulement autour de l'aile.
C'est la méthode utilisée par Nikolaï Joukovski au début du vingtième siècle pour les premières études sur l'aéronautique.
Du cercle au segment
[modifier | modifier le wikicode]Pour calculer l'écoulement autour du cercle unité, il suffit de trouver une transformation conforme qui transforme un segment (dont l'écoulement est trivial) en le cercle unité. Or l'image d'un point quelconque eit du cercle unité par la transformée de Joukovski est eit+e-it=2 cos(t): La transformation de Joukovski J(z) transforme le cercle unité en le segment [-2;2].
Et l'écoulement autour du segment est le plus simple possible, puisque c'est comme si le segment n'était pas là: Des droites verticales représentent les équipotentielles (ou isobares), et des droites verticales représentent les lignes de courant (ou écoulement) qui montrent la vitesse du vent. La transformation conforme qui produit cette grille est la transformation identique id(z)=z.
Du segment au cercle
[modifier | modifier le wikicode]Transformée inverse
[modifier | modifier le wikicode]Il suffit donc d'appliquer à cet écoulement l'inverse J-1 de la transformée de Joukovski, qui envoie le segment [-2;2] sur le cercle unité. Or J-1 est multiforme. On utilisera donc la fonction fplus valant sur la moitié droite de l'image, et la fonction fmoins valant sur la moitié gauche. Ces fonctions sont complexes:
require 'cmath'
def fplus(z)
return (z+CMath.sqrt(z**2-4))/2
end
def fmoins(z)
return (z-CMath.sqrt(z**2-4))/2
end
R=100
Le paramètre R est un facteur d'échelle permettant de zoomer sans avoir à refaire tout le script.
La figure est enregistrée dans un fichier au format svg appelé JoukovskiRuby1.svg:
figure=File.open("JoukovskiRuby1.svg","w")
figure.puts('<?xml version="1.0" encoding="utf-8"?>')
figure.puts('<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"')
figure.puts('"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">')
figure.puts('<svg xmlns="http://www.w3.org/2000/svg" width="640" height="480">')
La figure occupera donc 640 pixels de large et 480 pixels de haut, et s'appellera figure tout simplement.
Tracé des isobares
[modifier | modifier le wikicode]Côté droit de la figure
[modifier | modifier le wikicode]Pour tracer les isobares, on prend un point d'affixe x-2i et on le transforme par J-1. Puis on commence un path de svg par ce point. Ensuite on transforme des points situés sur la même droite verticale que le premier point choisi, et on les ajoute au path (une chaîne de caractères appelée ligne). La courbe obtenue (le path) sera en jaune et est inscrite à la fin dans le fichier:
#moitié droite des isobares
(0..25).collect { |i|
z=fplus(Complex(i/10.0,-2.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne='<path d="M'
ligne+=xa.to_s+' '+ya.to_s
(-20..20).collect { |j|
z=fplus(Complex(i/10.0,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne+=' L'+xa.to_s+' '+ya.to_s
}
ligne+='" stroke="yellow" stroke-width="1" fill="none"/>'
figure.puts(ligne)
}
Côté gauche de la figure
[modifier | modifier le wikicode]Les isobares de gauche se dessinent de la même manière, en remplaçant fplus par fmoins:
#moitié gauche des isobares
(-25..0).collect { |i|
z=fmoins(Complex(i/10.0,-2.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne='<path d="M'
ligne+=xa.to_s+' '+ya.to_s
(-20..20).collect { |j|
z=fmoins(Complex(i/10.0,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne+=' L'+xa.to_s+' '+ya.to_s
}
ligne+='" stroke="yellow" stroke-width="1" fill="none"/>'
figure.puts(ligne)
}
Tracé des lignes de courant
[modifier | modifier le wikicode]Pour les mettre en exergue, on fera les lignes de courant en rouge.
Côté gauche de la figure
[modifier | modifier le wikicode]Il suffit d'intervertir les rôles de i (abscisse) et j (ordonnée):
#moitié gauche de l'écoulement
(-25..25).collect { |j|
z=fmoins(Complex(-3.0,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne='<path d="M'
ligne+=xa.to_s+' '+ya.to_s
(-30..0).collect { |i|
z=fmoins(Complex(i/10.0-0.000001,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne+=' L'+xa.to_s+' '+ya.to_s
}
ligne+='" stroke="red" stroke-width="1" fill="none"/>'
figure.puts(ligne)
}
Le décalage des abscisses permet d'éviter le passage par l'origine, où le changement de feuillet de la transformation donnerait des résultats graphiquement bizarres.
Côté droit de la figure
[modifier | modifier le wikicode]Même principe, en remplaçant fmoins par fplus:
#moitié droite de l'écoulement
(-25..25).collect { |j|
z=fplus(Complex(0,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne='<path d="M'
ligne+=xa.to_s+' '+ya.to_s
(0..30).collect { |i|
z=fplus(Complex(i/10.0,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne+=' L'+xa.to_s+' '+ya.to_s
}
ligne+='" stroke="red" stroke-width="1" fill="none"/>'
figure.puts(ligne)
}
Tracé du cercle
[modifier | modifier le wikicode]Pour ajouter le cercle lui-même sur la figure, on peut utiliser l'objet circle de svg; on en profite pour fermer la bannière svg et clore le fichier:
figure.puts('<circle cx="320" cy="240" r="'+R.to_s+'" fill="white" stroke="black" stroke-width="1" />')
figure.puts('</svg>')
figure.close
C'est tout ce qu'il faut pour obtenir l'écoulement que voici:
Le défaut des isobares en-dessous du cercle (la figure aurait dû être symétrique par rapport à l'axe des abscisses) est lié au changement de feuillet de J-1. Pour y remédier, on peut remplacer la moitié du bas par la symétrique de la moitié du haut, ce qui double la longueur du script (une fonction par quadrant). La correction est laissée en exercice.
Conception de l'aile d'avion
[modifier | modifier le wikicode]Si on applique à la figure ci-dessus (avec le cercle unité), la transformation J(z), on obtient l'écoulement autour de [-2;2] qui est trivial. Mais on peut modifier le cercle pour qu'il passe toujours par le point d'affixe 1, et la transformée de Joukovski du nouveau cercle aura toujours un point de rebroussement d'affixe 2, mais peut avoir une forme de lemniscate, ou de profil d'aile d'avion, selon le paramètre choisi.
Tracé du profil
[modifier | modifier le wikicode]Une fonction affine complexe laisse le point d'affixe 1 invariant, si elle est de la forme fa(z)=P(z-1)+1. Le meilleur moyen de choisir la valeur du coefficient directeur P est d'implémenter fa et J avec un logiciel de géométrie dynamique (par exemple, CaRMetal qui exporte au format svg), et de bouger le point d'affixe P avec la souris jusqu'à ce que l'image du cercle par J(fa(z)) ait la forme voulue:
Le choix de P=0,93+0,15i paraît satisfaisant.
Calcul de l'écoulement
[modifier | modifier le wikicode]Pour la fonction de Joukovski, on ajoute un très petit nombre à z pour que le nombre zéro passe entre les mailles du filet, ce qui évite le risque d'avoir une erreur de division par zéro. À part ça, il suffit d'appliquer J(fa(z)) au cercle unité, et de remplacer le cercle final par son image par J(fa(z)) (pour dessiner l'aile). Ce qui donne le script suivant:
require 'cmath'
def fa(z)
return Complex(0.93,-0.15)*(z-1)+1
end
def Joukovski(z)
return z+1/(z+0.0000001)
end
def fplus(z)
return Joukovski(fa((z+CMath.sqrt(z**2-4))/2))
end
def fmoins(z)
return Joukovski(fa((z-CMath.sqrt(z**2-4))/2))
end
R=100
figure=File.open("JoukovskiRuby2.svg","w")
figure.puts('<?xml version="1.0" encoding="utf-8"?>')
figure.puts('<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"')
figure.puts('"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">')
figure.puts('<svg xmlns="http://www.w3.org/2000/svg" width="640" height="480">')
#moitié gauche de l'écoulement
(-25..25).collect { |j|
z=fplus(Complex(-3.0,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne='<path d="M'
ligne+=xa.to_s+' '+ya.to_s
(-30..0).collect { |i|
z=fplus(Complex(i/10.0-0.000001,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne+=' L'+xa.to_s+' '+ya.to_s
}
ligne+='" stroke="red" stroke-width="1" fill="none"/>'
figure.puts(ligne)
}
#moitié droite de l'écoulement
(-25..25).collect { |j|
z=fmoins(Complex(0,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne='<path d="M'
ligne+=xa.to_s+' '+ya.to_s
(0..30).collect { |i|
z=fmoins(Complex(i/10.0,j/10.0))
xa=z.real*R+320
ya=240-z.imag*R
ligne+=' L'+xa.to_s+' '+ya.to_s
}
ligne+='" stroke="red" stroke-width="1" fill="none"/>'
figure.puts(ligne)
}
ligne='<path d="M'+(320+2*R).to_s+' 240'
(0..200).collect { |i|
t=Math::PI*i/100
z=Joukovski(fa(Complex(Math.cos(t),Math.sin(t))))
xa=z.real*R+320
ya=240-z.imag*R
ligne+=' L'+xa.to_s+' '+ya.to_s
}
ligne+='" stroke="black" stroke-width="1" fill="white"/>'
figure.puts(ligne)
figure.puts('</svg>')
figure.close
Figure obtenue
[modifier | modifier le wikicode]L'exécution du script produit la figure suivante:
Le fait que les lignes de courant sont plus courbées sur le dessus que sur le dessous entraîne par la force de Bernoulli, une dépression sur le dessus de l'aide, qui est à l'origine de la portance: L'avion vole!