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:

  1. On utilise l'inverse de la transformée de Joukovski pour calculer l'écoulement autour du cercle unité.
  2. On agrandit légèrement ce cercle pour qu'il ait une transformée de Joukovski plus intéressante;
  3. 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

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

Transformée inverse

modifier

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

Côté droit de la figure

modifier

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

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

Pour les mettre en exergue, on fera les lignes de courant en rouge.

Côté gauche de la figure

modifier

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

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

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

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

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

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

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!