二次元翼理論

2022 AdventCalender 数学・物理 の10日目です。

二次元完全流体

流体の基礎方程式は解きづらいことで有名(?)ですが、完全流体の仮定を置いてよいときにはかなり簡単化されます。完全流体とは、

• 非圧縮：密度が一定
• 粘性なし：応力が圧力による面法線方向のみ

を満たす流体のことです。非圧縮なので密度$\rho$は考えなくて良くなります。

流体の速度場$\mathbf{v}$の回転と発散、$\nabla \times \mathbf{v},\nabla \cdot \mathbf{v}$をそれぞれ渦、湧き出しと呼びます。湧き出しはそのまま流体が発生したり消滅したりすることなので、特別そのような装置を考えたりしなければ$0$を仮定してよいでしょう。加えて、しばし渦がない流れを考えます。実際には渦は現実の流体では度々発生して重要な役割を果たしますが、完全流体では渦は局所的にも保存することが知られています。ということで、ある時点で渦なしの状態が実現していたとして、渦なしの流れを考えます。(渦なしが期待できる状況がどういうものなのかちょっと筆者もわかっていません)。

さらに速度場が時間的に定常であるような状況を考えます。

ここまで仮定するとベルヌーイの定理から、

$$p + \frac{1}{2}\rho \|\mathbf{v}\|^2 + \rho f = \mathrm{const.}$$

が成り立つことが知られています。ここで$f$は外力ポテンシャルとします。つまり、圧力は速度場から明示的に与えられるため、流体の運動は速度場$\mathbf{v}=(v_x,v_y,v_z)^\top$だけを考えれば良くなります。空気など重力の影響が小さいなどで外力を無視する場合は$f$も省略します。

回転と発散がないということは、速度場は調和関数ということになります。このため、境界条件を色々固定するだけで速度場の解が決まってしまい、流体の運動方程式を実質的に考える必要がなくなります。

さらに3次元方向を無視します。つまり、$v_z = 0$固定とし、$v_x,v_y$も$z$方向に均一だとします。ここまで状況を理想化すると、速度場を今度は複素関数を用いて解析できるようになり、具体的にはある特徴を持つ断面をもった翼の性能を一般的に計算することができる ようになります。この二次元翼理論について解説を行います。予備知識として、ベクトル解析(発散回転と各種積分を知っている)・複素解析(積分定理とローラン展開)の初歩を仮定します。

以下適当に数学的な話もしますが、もし数学的な話に飽きたら最後の節に飛んで目の保養をしましょう。

二次元速度場と複素関数

二次元速度場$\mathbf{v}=(v_x,v_y)^\top$を微分形式で

$$v = v_x dx + v_y dy$$

と表現することにします。このとき、回転(渦)は

$$dv = (\partial_x v_y - \partial_y v_x) dx\wedge dy = (\nabla \times \mathbf{v})dx\wedge dy$$

発散(湧き出し)は

$$\delta v = \partial_x v_x + \partial_y v_y$$

に相当します。ただし二次元であるため、ベクトル場の回転はその$z$成分スカラーとします。
なおホッジ作用素は

$$\star dx = dy\\ \star dy = -dx$$

と振る舞い$\star^{-1} = (-1)\star$です。

ベクトル解析的には$(x,y)$平面単純閉軌道$C=\partial D$による２種類の積分とそのストークス定理がありますが、これは以下のように対応します。まず軌道線積分について

$$\int_C \mathbf{v} \cdot d\mathbf{r} = \int_D (\nabla\times \mathbf{v}) dS = \int_C v = \int_D dv$$

軌道法線積分について、$\bm{n}$を、軌道$dl$の右側法線として

$$\int_C \mathbf{v} \cdot dl \mathbf{n} = \int_D (\nabla\cdot \mathbf{v}) dS = \int_C \star v = \int_D \star\delta v$$

速度場に対するこれらをそれぞれ循環と流量と呼ぶことにします。

渦なし$dv = 0$であれば、単連結領域においては$v = d\phi$を満たすポテンシャル(0形式スカラー)$\phi$があります。また湧き出しなし$\delta v = 0$であれば、単連結領域においては、$v = -\delta (\psi dx\wedge dy)$を満たすスカラー$\psi$があります。これを流れ関数と言います。

ポテンシャル、流れ関数がそれぞれ存在する場合は

$$v_x = \partial_x \phi = \partial_y\psi\\ v_y = \partial_y \phi = - \partial_x \psi$$

となります。渦なし湧き出しなしの完全流体を考えますが、二次元翼理論では、平面上に翼断面が穴を開けているようなシチュエーションを考えます。このため、速度場が渦なし、つまり環形式であっても、境界形式とは限りません。これは翼の周りの線積分が0でない=循環が存在し得ることを意味しますが、このことが翼理論では重要になってきます。

以下、平面上の区分的滑らかな単純閉曲線$W$を左回りにとって、これを翼断面の境界線とします。$W$の表面上の境界条件を考慮しておきましょう。まず、$W$法線方向の速度はゼロです。ゼロでなければ流体が翼にめり込んでしまいます。$(dx,dy)^\top$をこの翼境界の線素とすると

$$v_x dy - v_y dx = 0$$

という条件にあたります。つまり法線ベクトルが

$$dl\mathbf{n} = (dy,-dx)^\top$$

で、それと直交するということです。従って翼型境界での速度場の軌道法線積分は常にゼロになります。

$$\int_W \mathbf{v} \cdot dl \mathbf{n} = 0$$

ところで($\mathbf{v}\cdot dl \mathbf{n}$)の積分は一般に

$$\int \nabla \psi \cdot d\mathbf{r}$$

と書き直すことができます。ここで、$C$が任意の閉軌道だとしましょう。それが翼型を含まないときは、湧き出し無し条件から常にゼロです。もし翼型を含むなら、積分値に影響することなく翼型境界の軌道$W$を追加できます。そうすると軌道$C-W$はある領域の境界になるので、これを$D^\prime$としましょう。つまり$C-W = \partial D^\prime$です。この積分はふたたび湧き出し無し条件から結果はゼロです。

$$\int_C \nabla \psi \cdot d\mathbf{r} = \int_{C-W} \nabla \psi \cdot d\mathbf{r} = \int_{C-W} \star v =\int_{D^\prime} \star\delta v = 0$$

$C$は任意だったので、これは$\psi$が大域的に一価なことを意味します。

したがって、翼型境界で流体がめり込まない境界条件を課しておけば、流れ関数$\psi$は常に存在するものとして考えることができます。流れ関数からは速度場が得られるので、速度場の代わりに流れ関数を考えてもよいことになります。

二次元では、速度場を複素関数として表現できます。$z= x+iy$とし

$$w(z) = v_x(x,y) - i v_y(x,y)$$

を複素速度といいます。ここで残念なお知らせがあり、$v_y$側の係数が$-i$であるのはタイポではありません。ここには分野の慣習の衝突が起きています。

よいお知らせは、渦なし湧き出し無し条件が実は$w(z)$のコーシーリーマン関係式そのものだと言うことです。したがって$w(z)$は正則関数です。これは速度場が調和であったことから期待されたことです。

翼境界上の境界条件を複素関数の観点から書き直しておくとあとで便利です。$dz = dx + i dy$とすると、

$$v_x dy - v_y dx = 0$$

とは

$$\mathrm{Im} w^*(z) dz^* = 0$$

と書き直すことができます。これは

$$w(z) dz = w^*(z) dz^*$$

と同じことです。

複素速度については正則なので複素積分をしたくなります。これを実行すると

$$w(z) dz = (v_x - iv_y)(dx + idy)\\ = (v_x dx + v_y dy) + i(v_x dy - v_y dx) \\ = \bm{v}\cdot d\bm{r} + i \bm{v}\cdot dl\bm{n}$$

となるので循環は

$$\mathrm{Re} \int_C w(z) dz = \int_C \mathbf{v}\cdot d\mathbf{r}$$

であり、流量は

$$\mathrm{Im} \int_C w(z) dz = \int_C \mathbf{v}\cdot dl\mathbf{n}$$

と複素積分を用いて表現できます。

ブラウジスの定理

今、二次元平面/複素平面があって、そこに翼断面に相当する領域とその境界である区分的連続な単純閉曲線$W =\partial D$があります。速度場$\bm{v}$はこの外側で定義されており、これを複素関数で表現した場合の$w(z)$は正則です。

このような状況ではローラン展開が可能です。ただし、ここで追加で次のように境界条件を課してしておきましょう。$z\rightarrow \infty$で$w(z)$が収束するものとします。つまり、無限遠で有界な値をとるものとします。また、原点は翼型の内側にとるものとします。

そうなると、$w(z)$のローラン展開は負のべきのみとなります。

$$w(z) = \sum_{n=0}^\infty b_n z^{-n}$$

$b_0$は、物理的には、無限遠では一様な速度場$(\mathrm{Re}b_0,-\mathrm{Im}b_0)$であるという境界条件を表現しています。

翼型を左に回る軌道で複素速度を積分した場合、積分定理から

$$\int_C w(z)dz = 2\pi i b_1$$

です。ところで、翼型の表面での流量は0であり湧き出しもないので、この値の虚部は0であり、実部は循環です。したがって$\Gamma$を循環(実数値)として、

$$b_1 = -\frac{i\Gamma}{2\pi}$$

とおくことができます。

翼は剛体だとします。そうすると力学的に重要なのは、かかる力の総和と、モーメントです。

流体への外力を無視してよいなら、ベルヌーイの定理から、圧力が

$$p(z) = -\frac{1}{2}\rho \|\mathbf{v}\|^2 + \mathrm{const.}$$

と分かっているので、この速度場から力とモーメントを計算することができます。$W$を再び翼型の境界とします。法線は外向きなので、翼にかかる力はこれに圧力の負をかけたものの総和です。

$$(F_x,F_y)^\top = - \int_W p(x,y)dl\mathbf{n}\\ = \rho\int_W \left(\frac{1}{2} \|\mathbf{v}\|^2 + c\right) dl\mathbf{n} \\ = \rho\int_W \frac{1}{2} \|\mathbf{v}\|^2 dl\mathbf{n}$$

定数項を単純閉曲線で線積分すると相殺して消えるため、これは無視します。これを複素積分で書き直します。

$$F_x - i F_y = \rho\int_W \frac{1}{2} \|w(z)\|^2 (dy+idx)\\ = \frac{i\rho}{2} \int_W w(z)w^*(z) dz^*\\ = \frac{i\rho}{2} \int_W (w(z))^2 dz$$

ここで、翼型境界では

$$w(z)dz = w^* (z)dz^*$$

が成り立つことを用いています。
同様に、モーメントについても

$$M = -\int_W p \mathbf{r}\times dl\mathbf{n}\\ = \rho\int_W \left(\frac{1}{2}\|\mathbf{v}\|^2+c\right) \mathbf{r}\times dl\mathbf{n}$$

ですが、

$$\int_W \mathbf{r}\times dl \mathbf{n}\\ = \int_W \mathbf{r} \cdot d\mathbf{r}\\ = \int_D 0 dS = 0$$

なのでこの項を再び無視すると、

$$M = -\rho\int_W \frac{1}{2}\|\mathbf{v}\|^2 (xdx + ydy)\\ = -\rho\int_W \frac{1}{4}w(z)w^*(z) d(zz^*)\\ = -\rho\int_W \frac{1}{4}w(z)w^*(z) (zdz^* + z^* dz)\\ = -\rho\int_W \frac{1}{4}((w(z))^2zdz + (w^*(z))^2z^*dz^*)\\ = -\mathrm{Re}\frac{\rho}{2}\int_W (w(z))^2zdz$$

これらの公式はブラウジスの定理と呼ばれています(クラウジウスではありません)。これは複素速度の複素積分になっており、正則関数の積分なので、積分軌道を翼型境界から、翼を含む任意の左回り閉曲線に取り直してもよいです。ローラン展開を考慮すると、負べきの複素関数の閉軌道積分では、$-1$乗以外の項は残りません。したがって、翼からみれば、複素速度のローラン展開のうち、$0,-1,-2$次の項($b_0,b_1,b_2$)だけが、力とモーメントに関する情報を担っていることになります。

$$w(z)^2 = b_0^2 + 2b_0b_1z^{-1} + (2b_0b_2+b_1^2)z^{-2} + o(z^{-2})$$

例えば、翼の性能として揚力というものがありますが、揚力については循環$\Gamma$が重要です。揚力とは翼にかかる上向きの力ですが、ここでの"上"は対気速度に垂直な方向のことです。ブラウジスの定理から

$$F_x - iF_y = \frac{i\rho}{2}2 b_0 b_1 2\pi i\\ = - \rho b_0b_1 2\pi\\ = i \rho\Gamma b_0$$

$b_0$は無限遠での速度場ですが、これは対気速度を表すと考えられるので、水平方向に流れがあるとして$b_0 = U (\in \mathbb{R})$としましょう。そうすると

$$F_y = - \rho U\Gamma$$

となり、揚力が気体密度、翼まわりの循環、対気速度で決まることがわかります。これはクッタ・ジューコフスキーの定理と呼ばれています。(後に説明する同名の仮定とは別)。速度を上げれば揚力が得られるのは経験上も頷けます。その係数が翼まわりの循環$\Gamma$であるというのがポイントとなります。高い揚力を生み出す翼型は、翼まわりに効率的に循環を発生させるものということです。(符号がマイナスなのは、循環の正の向きが左回りであるためです。気流が右向き($U > 0$)のとき、右回り循環が起きれば上向きの揚力となります)

なおこのとき$F_x=0$になっており、抗力を一切受けないという直観に反する結果が出ます。これはダランベールのパラドックスで知られ、完全流体として粘性を無視したことに起因し、現実の翼は多少なり抗力を受けます。

リーマンの写像定理と円筒翼

リーマンの写像定理によれば、単連結(円盤と同相)領域の間の正則な同相写像がつねに存在します。今考えている速度場の定義域は、円盤領域というよりはその外側ですが、考えている複素速度が無限遠で収束するという条件を課しているので、無限遠点を追加して反転することで同様の主張を得ておきましょう。つまり、任意の単連結な翼型について、その外側領域同士を結ぶ正則な同相写像があることになります。

これは、二次元翼理論にある種の万能性を保証します。$z=x+iy$平面での翼型$W=\partial D$による複素速度が$w(z)$だとします。$\zeta = \xi + i \eta$平面上、翼型$X = \partial E$があるとします。このとき、翼領域$D,E$の外側を正則に写像する同相$f$があるはずです。そこで

$$x(\zeta) = w(f^{-1}(\zeta))\frac{df^{-1}}{d\zeta}(\zeta)$$

を考えると、これは$E$の外側で定義されていて正則であり、無限遠で有界、かつ、翼型に流体がめり込まないという境界条件を継承しています。

$$x(\zeta)d\zeta = w(f^{-1}(\zeta))dz = w^*(f^{-1}(\zeta))dz^* = x^*(\zeta)d\zeta^*$$

したがって、翼型$W= \partial D$において、速度場を完全に求めておき、別の翼型$X= \partial E$については写像定理で速度場を求める、という戦略が可能になります。$W= \partial D$は単純な形にして解析的に解を求めておき、実際に計算したい複雑な形$X= \partial E$については、$D^c$を$E^c$に写像する正則関数を見つけることで速度場を構成するということです。

円筒翼は対称性が高いため、速度場を簡単に書き下せます。天下りですが、原点の半径$a$の円筒翼まわりの循環$\Gamma$の流れは

$$w(z) = b_0 \left( 1 - \frac{a^2}{e^{2i\arg b_0}z^2} \right) - \frac{i\Gamma}{2\pi} \frac{1}{z}$$

で与えることができます。実際に$|z| = a$での速度場を計算すると

$$w(ae^{i\theta}) = |b_0|e^{i\arg b_0}(1-e^{-2i(\theta+\arg b_0)}) - \frac{i\Gamma}{2a\pi}e^{-i\theta}\\ v_x = \left(|b_0|\sin(\theta+\arg b_0) - \frac{\Gamma}{2a\pi} \right) \sin\theta\\ v_y = -\left(|b_0|\sin(\theta+\arg b_0) - \frac{\Gamma}{2a\pi} \right) \cos\theta$$

となり、円筒翼の接ベクトルに並行になっていて境界条件を満たしています。

これより高次の負冪項は、境界条件の観点からは線形独立になるため他の項をここに加えることはできません。したがって、円筒翼周りの複素速度場はこれで尽きることになります。(円筒翼境界条件を満たさせるという観点からみると、複素速度の$z$のベキは、次数が1ずれたフーリエsin級数のようなものです。したがって、$-n-2$次と$n$次が相殺することで境界条件を満足し、$-1$次が単独で境界条件を満たします。しかし今無限遠での境界条件から正のベキは採用できず、$0$次と$-2$次の組み合わせと、$-1$次だけが許容できる自由度になります。これが$b_0,\Gamma$に相当することになります)

クッタ・ジューコフスキー条件

前節で円筒翼の速度場を与えました。しかし 円筒翼の速度場には、その循環パラメータ$\Gamma$が不定のまま残っています。実は循環$\Gamma$は正則関数で写像しても変わりません。

$$\int_C w(z)dz = \int_{f(C)} w(f^{-1}(\zeta))\frac{df^{-1}}{d\zeta}(\zeta)d\zeta = \int_{f(C)}x(\zeta) d\zeta$$

つまり、ここから任意の翼型に写像してその翼型について調べようとしても、もっとも重要な数値が不明なままです。 翼型境界方向の境界条件を決める必要があります。

ところが 完全流体の仮定のもとでは、翼表面付近の流体の振る舞いを記述することができません。 実際には粘性などによって翼の表面付近に薄い渦の層が発生し、これによって翼表面でいくらか流体が「滑る」ことができ....といった感じで翼境界方向の境界条件が決まるはずですが、完全流体は粘性を考慮に入れていない上、今渦なしの流れを考えているため、そうした議論ができません。

しかし実は、完全流体のまま、こうした境界条件を固定する絶妙な経験則があります。これがクッタ・ジューコフスキーの条件です。このポイントは、翼型を「尖らせる」ことにあります。

翼の後端を非常に薄く尖らせた場合、翼の上面と下面からきた流れは「なめらかに合流する」と期待できます(これこそが経験則であり、ここでは証明されません)。もとの円筒翼の速度場の循環項$\Gamma$を色々に変えると、実は翼の上面と下面からきた流れは、一般の値ではなめらかに合流せず、特定の値でだけ滑らかに合流します。 (この様子をあとで可視化します)

この滑らかに合流するときの$\Gamma$を採用することで、速度場が決定され、翼に掛かる力やモーメントも算出できるようになります。

逆に言えば、翼型後端が鋭く尖っていない場合、クッタ・ジューコフスキーの条件を運用できません。クッタ・ジューコフスキーの条件は、本来完全流体の観点からは議論できない翼境界付近での挙動を、翼境界を「摘んで押しつぶす」ことで、「翼後端付近の流れが滑らかに合流しているか」という単純な規準に落とし込んでしまうという(やや強引な)テクニックです。

この代償として、クッタ・ジューコフスキーの仮定で速度場を決定可能な翼型は、常に鋭く薄い後端をもつことになります。薄いということは強度を稼げないという工学上のデメリットがあります。

ジューコフスキー翼の可視化

クッタ・ジューコフスキーの条件は前節のように言葉で説明されてもピンと来ないと思われます。そこで、具体的にこの条件を運用するとどうなるのか可視化しましょう。

円/翼型の外側領域の間の正則写像をもち、翼型の一点が特異的に尖った形をもつ翼型としてジューコフスキー翼があります。そしてこれを実現する関数もジューコフスキー変換と呼ばれています。これを精査することで、ここまでの話を実際に実現してみましょう。つまり、写像定理で得られる写像がジューコフスキー変換だった場合に、実際にクッタ・ジューコフスキー条件で循環を決定してみます(二次元翼理論としては、一般に一点で尖った翼型とそれを円筒からの写像として与える正則関数があれば、同様の議論ができますが、具体例としてジューコフスキー変換を扱うということです)。ここで、ジューコフスキー変換とは以下です。

$$J(z) = z + \frac{a^2}{z}$$

これの逆関数は二価です。

$$J^{-1}_\pm (\zeta) = \frac{\zeta\pm \sqrt{\zeta^2 -4a^2}}{2}$$

解と係数の関係から

$$J^{-1}_+(-) J^{-1}_-(-) = a^2$$

が成り立ちます。

半径$a$の円は、この変換では$x$軸線分に潰れます。

$$J(\{z:|z| = a\}) = \{x+i0: x \in [-2a,2a]\}$$

ジューコフスキー翼とは、複素平面上$(a,0)$と交差するような円の、$J$による像のこととします。 よってその自由度は二次元です。ただしこの中心$c$は$\mathrm{Re} c < \frac{a}{2}$に制限します。あとで可視化しますが、この制限の外側だと、ジューコフスキー翼が空間全体に広がってしまい、有界な翼にならないからです。この円を便宜的にパラメータ円と呼ぶことにしましょう。そしてパラメータ円を$z \mapsto \frac{a^2}{z}$で反転した円を、便宜的に反転円とよび、写像前の円筒翼として用います。

これも図示することでわかるのですが、$c$の実部が負のときでは、反転した円が元の円の内側に来て、翼型としては同じものが手に入ります。さらに$c$の複素共役を取るとこれは単に上下を反転しただけなので翼型も相似になります。したがって、ジューコフスキー翼の対応するパラメータ円について、$c$をその中心として、

$$\mathrm{Re}c \in \left[ 0,\left.\frac{a}{2}\right.\right)\\ \mathrm{Im}c \in \left[ 0,\left.\infty\right.\right)$$

に制限しても一般性を失いません。

この条件を満たしているとき、円と反転円はともに原点を含みます。ジューコフスキー変換は明らかに原点で正則でないので、原点を円筒翼で隠す必要がありますが、反転円を円筒翼とみなすことでこれが満足されます。

このとき、パラメータ円の半径は$|c-a|$ですが、幾何学的考察により、反転円の半径は$\frac{|c-a|a}{a - 2\mathrm{Re}c}$で、中心は$a+ \frac{(c-a)a}{a-2\mathrm{Re}c}$となります。

ジューコフスキー変換は逆関数が二価なので、ここまでに説明した方法で翼周りの速度場を求めるには、実際に考慮する翼型を考慮に入れた上で、逆像を適切に選ぶ必要があります。

これ以上は数式をこねくり回して理解するのは直観が追いつかず大変です。ということで、計算機の力で可視化します。

このあたりの可視化の良いツールがあるのかもしれませんが、平面上の各点で計算を高速にできるのでGLSLを用います。

VScodeのExtensionでShaderToyExtensionを用いると、以下のコードを適当な新規ファイルにコピペ+Preview(cmd or ctl + shift + P, Show GLSL Preview)で動作確認ができます。

GLSLでジューコフスキー変換とその逆変換を実装します。まずはパラメータ円と、その反転円と、ジューコフスキー翼がどのような形になるか見てみましょう。

const vec4 c1 = vec4(1.0,0.0,0.0,1.0);
const vec4 c2 = vec4(0.0,1.0,0.0,1.0);
const vec4 c3 = vec4(0.0,1.0,1.0,1.0);
const vec4 c4 = vec4(1.0,0.0,1.0,1.0);
const vec4 c5 = vec4(0.5,0.5,1.0,1.0);
const vec4 c6 = vec4(1.0,1.0,0.0,1.0);
const float a = 2.5;
const float eps = 0.01; 


vec2 anotherPair(vec2 pt) {
    return vec2(pt.x,-pt.y)*a*a/pow(length(pt),2.0);
}

vec2 reScale(vec2 pt,float scale){
    float l = min(iResolution.x,iResolution.y);
    return (pt.xy - (iResolution.xy / 2.0)) / (l*scale);
}

bool onCircle(vec2 pt,float r,vec2 center) {
    return abs(length(pt-center) - r) < eps;
}
bool inCircle(vec2 pt,float r,vec2 center) {
    return length(pt-center) < r;
}

vec2 joukowskiForward(vec2 pt) {
    float r = length(pt);
    float x = pt.x + ((a*a)*pt.x)/(r*r);
    float y = pt.y - ((a*a)*pt.y)/(r*r);
    return vec2(x,y);
}

const float PI = 3.14159265359;

vec2 poler(float r, float theta) {
    return vec2(r*cos(theta),r*sin(theta));
}

vec2 sqrtZZminusAA4(vec2 pt,float a) {
    float r = length(pt);
    float d = pow(
        pow(r,4.0)
        -8.0*pow(a,2.0)*(pt.x*pt.x - pt.y*pt.y)
        +16.0*pow(a,4.0)
        ,0.25);
    float zzminusaa4_y = (2.0*pt.x*pt.y);
    float zzminusaa4_x = ((pt.x*pt.x - pt.y*pt.y) - (4.0 * pow(a,2.0)));
    float angle = atan(zzminusaa4_y/zzminusaa4_x);
    bool cos2thetaover = zzminusaa4_x > 0.0;
    float angleModified;
    if (pt.x > 0.0 && cos2thetaover) {
        angleModified = 0.5 * angle;
    }
    if (pt.y > 0.0 && !cos2thetaover) {
        angleModified = 0.5 * (PI+angle);
    }
    if (pt.y < 0.0 && !cos2thetaover) {
        angleModified = 0.5 * (-PI+angle);
    }
    if (pt.x < 0.0 && pt.y > 0.0 && cos2thetaover) {
        angleModified = 0.5 * (2.0*PI+angle);
    }
    if (pt.x < 0.0 && pt.y < 0.0 && cos2thetaover) {
        angleModified = 0.5 * ((-2.0*PI)+angle);
    }
    return poler(d,angleModified);
}

vec2 joukowskiInv1(vec2 pt) {
    return (pt + sqrtZZminusAA4(pt,a))*0.5;
}

vec2 joukowskiInv2(vec2 pt) {
    return (pt - sqrtZZminusAA4(pt,a))*0.5;
}

void main(){
    float scale = 0.1;
    vec2 pt = reScale(gl_FragCoord.xy,scale);
    vec2 mouse = reScale(iMouse.xy,scale);
    vec2 anotherMouse = anotherPair(mouse);
    gl_FragColor = vec4(1.0,1.0,1.0,1.0);

    // パラメタ円の反転円
    float outer_rate = a / (a - 2.0 * mouse.x);
    float outer_radias = length(mouse-vec2(a,0)) * abs(outer_rate);
    vec2 outer_center =vec2(a,0.0) + ((mouse-vec2(a,0))* outer_rate);

    // 翼の内側の逆写像：ただし、パラメータの実部が正のときだけ動作する。
    if (!inCircle(pt,length(mouse-vec2(a,0)),mouse) && (outer_rate<0.0) != inCircle(pt,outer_radias,outer_center)) {
        gl_FragColor = c6;
    }
    // 翼の内側：ただし、パラメータの実部が正のときだけ動作する。
    vec2 j1 = joukowskiInv1(pt);
    vec2 j2 = joukowskiInv2(pt);
    if (length(j1-mouse) > length(mouse-vec2(a,0)) && length(j2-mouse) > length(mouse-vec2(a,0)) ) {
        gl_FragColor = c5;
    }

    // 座標軸
    if (abs(pt.x) < eps || abs(pt.y) < eps) {
        gl_FragColor = vec4(0.0,0.0,0.0,1.0);
    }

    // 原点のa円
    if (onCircle(pt,a,vec2(0,0))) {
        gl_FragColor = vec4(0.0,0.0,0.0,1.0);
    }
    // パラメータ円
    if (onCircle(pt,length(mouse-vec2(a,0)),mouse)) {
        gl_FragColor = c1;
    }
    // の中心 
    if (length(pt-mouse)< eps+0.03) {
        gl_FragColor = c1;
    }
    // 反転円
    if (onCircle(pt,outer_radias,outer_center)) {
        gl_FragColor = c2;
    }
    // 翼境界1
    if (onCircle(joukowskiInv1(pt),length(mouse-vec2(a,0)),mouse)) {
        gl_FragColor = c3;
    }
    // 翼境界2
    if (onCircle(joukowskiInv2(pt),length(mouse-vec2(a,0)),mouse)) {
        gl_FragColor = c4;
    }
}

逆変換の実装がちょっと面倒です。根気よく場合分けをします。ジューコフスキー変換はあくまで翼型を与える関数の一例だといいましたが、逆変換をする必要があることを考えると正則といえあまり複雑な関数は使いたくないですね...ともかく表示しましょう。

コードではExtensionの機能で、マウス座標をパラメータ円の中心$c$と取っているので、このコードをコピペしてShaderToyExtension Previewでクリックするといろいろなジューコフスキー翼を表示できます。

ちょっと太い。

$c$が$y$軸に近づくと薄くなります。

これはイケメン。

$\frac{a}{2}$を超えると翼が非有界になります。

ここで、赤い色が、パラメータ円とその中心です。緑は反転円かつ参照する円筒翼であり、ジューコフスキー変換で赤い円と同じ場所に写像されます。黒の円は半径$a$の円、ピンクとスカイブルーで囲まれた青のところがジューコフスキー翼に対応します。

パラメータ円(赤)を$(a,0)$で交差させたために、ジューコフスキー円は鋭い後端を持っています。

ジューコフスキー逆変換を調べましょう。ジューコフスキー翼の輪郭はわかりましたが、この翼の内側と外側は写像前のどこに対応するでしょうか？

写像前後の領域を追ってみると、実はジューコフスキー翼の内側の逆像とは、パラメータ円(赤)の反転円(緑)から、元のパラメータ円(赤)を引いた領域 (上図の黄色部分)だとわかります。　上図でいいう青領域と黄領域が1:2で対応しているということです(黄色の三日月(?)領域はジューコフスキー変換によって、半径$a$の円を折り目として自己交差して翼型へ折りたたまれます)

$J^{-1}_\pm(\zeta)$は、実は、「ともに反転円に含まれてパラメータ円に含まれない」か「一方がパラメータ円に含まれ、もう一方は反転円の外側に居る」かのどちらかの場合しかありません。また、この時逆変換の２つの像は、半径$a$の内側と外側に一つずつ現れます。

このことから、ある点がジューコフスキー翼の内側であることは 逆変換$J^{-1}_\pm(\zeta)$の2つの値がともにパラメータ円に含まれないという条件で表すことができます。 反転円を円筒翼として参照するので、翼の外側のジューコフスキー逆変換については、$J^{-1}_\pm(\zeta)$のうち、反転円に入らない方を常に採用すればよいです。

円筒翼の速度場を使って、一般翼形の流れを求めるために、円筒翼の速度場も可視化しましょう。ところで、速度場の可視化をするには、流れ場そのものよりも、流れ関数$\psi$のほうを求めておくほうが便利です。流れ関数の読み方としては、その等高線は流線で、流線の密度が流れの速さとなります。
円筒翼の複素速度は

$$w(z) = b_0 \left( 1 - \frac{a^2}{e^{2i\arg b_0}z^2} \right) - \frac{i\Gamma}{2\pi} \frac{1}{z}$$

でしたが、これに対応する流れ関数は

$$\psi = |b_0|\left(r-\frac{a^2}{r}\right)\sin(\theta +\arg b_0)-\frac{\Gamma}{2\pi} \log r$$

と与えることができます。ここで、$b_0$は$(\mathrm{Re}b,-\mathrm{Im}b_0)$の無限遠での流れであることに注意します。つまり、$\arg b_0$は、無限遠での流れ角度ですが、これは通常の極座標の角度とは逆方向を正の向きとしています。

これもGLSLで描いてみましょう。GLSLは関数の等高線を表示するのが得意なのでこちらは楽勝です。

const float PI = 3.14159265359;


const float r = 1.0;
const vec2 center = vec2(0.0,0.0); 
const float gamma_coeff = 10.0;

float arg(vec2 pt) {
    if (pt.x>0.0) {
        return atan(pt.y/pt.x);
    } else {
        return PI - atan(-pt.y/pt.x);
    }
}

// 円柱まわりの循環を持つ流れ関数
float flowFunctionAroundCircle(vec2 pt,vec2 center,vec2 b_0,float a,float g) {
    float r = length(pt-center);
    float u = length(b_0);
    float angle = -arg(b_0);
    return u*(r-((a*a)/r))*sin(arg(pt-center)-angle)-(g/(PI*2.0))*log(r);
}


bool contour(float div, float width,float v) {
    return abs(fract(iTime + v/div) - width*0.5) < (width*0.5);
}

bool inCircle(vec2 pt,vec2 center,float a) {
    return length(pt-center) < a;
}

vec2 reScale(vec2 pt,float scale){
    float l = min(iResolution.x,iResolution.y);
    return (pt.xy - (iResolution.xy / 2.0)) / (l*scale);
}

void main(){
    float scale = 0.1;
    vec2 mouse = reScale(iMouse.xy,scale)-center;
    vec2 mouse_yinv = vec2(mouse.x,-mouse.y);
    vec2 center = vec2(0.0,0.0);
    vec2 pt = reScale(gl_FragCoord.xy,scale);
    // 円の内側
    float a = 1.0;
    if (inCircle(pt,center,a)) {
        gl_FragColor = vec4(0.0,0.0,1.0,1.0);
        return;
    }
    float f = flowFunctionAroundCircle(pt,center,mouse_yinv/length(mouse_yinv),a,-length(mouse)*gamma_coeff);
    if (contour(1.0,0.03,f)) {
        gl_FragColor = vec4(1.0,0.0,0.0,1.0);
    } else {
        gl_FragColor = vec4(1.0,1.0,1.0,1.0);
    }
}

コード上はShaderToyExtensionのuniform機能でアニメーションしていますが、このアニメーションの方向は流れの方向ではありません。赤い線が流れの方向で、その密度が速さです。循環を設定することで、円筒翼の上下の流線密度が変わっています。 密な方が速度が早く、圧力が低いため、圧力差によって揚力が発生します。循環が大きい場合は円筒翼の外側に淀み点が発生します。

さて、循環$\Gamma$はここまでの議論では決定できないのでした。

ジューコフスキー変換/逆変換と円筒翼の流れ関数が手に入ったので、円筒翼の流れ関数をジューコフスキー翼に流し込みましょう。

各座標について、ジューコフスキー逆変換を実施し、先の規準にしたがって、翼の内側になった場合はそれを描画します。外側になった場合は、反転円の外側を常に選択し、反転円を円筒翼とした流れ関数を返し、等高線を描画します。つまり、円筒翼流れ関数をジューコフスキー逆変換で引き戻します。

// ジューコフスキー翼
// 循環:
const float gamma =-2.0;
// 円筒翼のローランb_0値：迎角と速度
const vec2 b_0 = vec2(1.0,-0.4);
// ジューコフスキー変換の元円盤の中心(xは正で２倍がaを超えないこと)
const vec2 inner_center = vec2(0.3,1.0);
// ジューコフスキー変換と単位円の大きさ
const float a = 3.0;
// 描画スケール
const float scale = 0.05; // 描画サイズを変える。
const vec2 offset = vec2(0.0,0.0); // 描画位置をずらす。
const float countour_div = 1.0; // 大きくすると粗くなる。
// 適当な色
const vec4 c1 = vec4(1.0,0.0,0.0,1.0);
const vec4 c2 = vec4(0.0,1.0,0.0,1.0);
const vec4 c3 = vec4(0.0,1.0,1.0,1.0);
const vec4 c4 = vec4(1.0,0.0,1.0,1.0);
const vec4 c5 = vec4(0.5,0.5,0.5,1.0);

// 等高線を描画する場合の厚み
const float eps = 0.01; 

// a^2/z
vec2 anotherPair(vec2 pt) {
    return vec2(pt.x,-pt.y)*a*a/pow(length(pt),2.0);
}

vec2 reScale(vec2 pt,float scale){
    float l = min(iResolution.x,iResolution.y);
    return (pt.xy - (iResolution.xy / 2.0)) / (l*scale);
}

// 円周に乗ってる
bool onCircle(vec2 pt,float r,vec2 center) {
    return abs(length(pt-center) - r) < eps;
}

vec2 joukowskiForward(vec2 pt) {
    float r = length(pt);
    float x = pt.x + ((a*a)*pt.x)/(r*r);
    float y = pt.y - ((a*a)*pt.y)/(r*r);
    return vec2(x,y);
}

const float PI = 3.14159265359;

vec2 poler(float r, float theta) {
    return vec2(r*cos(theta),r*sin(theta));
}

// 複素関数(z^2-4a^2)^0.5
vec2 sqrtZZminusAA4(vec2 pt,float a) {
    float r = length(pt);
    float d = pow(
        pow(r,4.0)
        -8.0*pow(a,2.0)*(pt.x*pt.x - pt.y*pt.y)
        +16.0*pow(a,4.0)
        ,0.25);
    float zzminusaa4_y = (2.0*pt.x*pt.y);
    float zzminusaa4_x = ((pt.x*pt.x - pt.y*pt.y) - (4.0 * pow(a,2.0)));
    float angle = atan(zzminusaa4_y/zzminusaa4_x);
    bool cos2thetaover = zzminusaa4_x > 0.0;
    float angleModified;
    if (pt.x > 0.0 && cos2thetaover) {
        angleModified = 0.5 * angle;
    }
    if (pt.y > 0.0 && !cos2thetaover) {
        angleModified = 0.5 * (PI+angle);
    }
    if (pt.y < 0.0 && !cos2thetaover) {
        angleModified = 0.5 * (-PI+angle);
    }
    if (pt.x < 0.0 && pt.y > 0.0 && cos2thetaover) {
        angleModified = 0.5 * (2.0*PI+angle);
    }
    if (pt.x < 0.0 && pt.y < 0.0 && cos2thetaover) {
        angleModified = 0.5 * ((-2.0*PI)+angle);
    }
    return poler(d,angleModified);
}

// a円盤の外側にある方の解。
vec2 joukowskiInv1(vec2 pt) {
    return (pt + sqrtZZminusAA4(pt,a))*0.5;
}

// a円盤の内側にある方の解。
vec2 joukowskiInv2(vec2 pt) {
    return (pt - sqrtZZminusAA4(pt,a))*0.5;
}

// 偏角
float arg(vec2 pt) {
    if (pt.x>0.0) {
        return atan(pt.y/pt.x);
    } else {
        return PI - atan(-pt.y/pt.x);
    }
}

// 円柱まわりの循環を持つ流れ関数
float flowFunctionAroundCircle(vec2 pt,vec2 center,vec2 b_0,float circlerad,float g) {
    float r = length(pt-center);
    float u = length(b_0);
    float angle = -arg(b_0);
    return u*(r-((circlerad*circlerad)/r))*sin(arg(pt-center)-angle)-(g/(PI*2.0))*log(r);
}

// 等高線判定
bool contour(float div, float width,float v) {
    return abs(fract(iTime + v/div) - width*0.5) < (width*0.5);
}

// 円に属するか
bool inCircle(vec2 pt,vec2 center,float a) {
    return length(pt-center) < a;
}

void main(){
    // vec2 mouse = reScale(iMouse.xy,scale);
    float inner_radias = length(inner_center - vec2(a,0.0));
    vec2 pt = reScale(gl_FragCoord.xy,scale)-offset;
    // 背景と座標軸
    gl_FragColor = vec4(1.0,1.0,1.0,1.0);
    if (abs(pt.x) < eps*1.5 || abs(pt.y) < eps*1.5) {
        gl_FragColor = vec4(0.0,0.0,0.0,1.0);
    }
    // 逆変換
    vec2 j1 = joukowskiInv1(pt);
    vec2 j2 = joukowskiInv2(pt);

    // パラメタ円の反転円
    float outer_rate = a / (a - 2.0 * inner_center.x);
    float outer_radias = inner_radias * outer_rate;
    vec2 outer_center =vec2(a,0.0) + ((inner_center-vec2(a,0))* outer_rate);

    // 翼の内側判定と翼型塗りつぶし
    bool insideFoil = length(j1-outer_center) < outer_radias && length(j2-outer_center) < outer_radias;
    if (insideFoil) {
        gl_FragColor = c5;
    }
    // 翼の外側判定と円筒翼流れ関数の参照
    float f;
    if (!insideFoil) {
        // 反転大円を円筒翼として外側にあるものを参照
        if (length(j2-outer_center) > outer_radias) {
            f = flowFunctionAroundCircle(j2,outer_center,b_0,outer_radias,gamma);
        } 
        if (length(j1-outer_center) > outer_radias) {
            f = flowFunctionAroundCircle(j1,outer_center,b_0,outer_radias,gamma);
        } 
        if (contour(countour_div,0.03,f)) {
            gl_FragColor = vec4(1.0,0.0,0.0,1.0);
        }
    }
}

ジューコフスキー翼まわりに流線が描かれます。

コード上の

// 循環:
const float gamma =-2.0;
// 円筒翼のローランb_0値：迎角と速度
const vec2 b_0 = vec2(1.0,-0.4);
// ジューコフスキー変換の元円盤の中心(xは正で２倍がaを超えないこと)
const vec2 inner_center = vec2(0.3,1.0);

を変えることで、循環と無限遠の一様流の向き、パラメータ円の中心を調整します。ここで適当に描画すると、翼の後端で流れが折れ曲がっていることがあります。

そしてここがクッタ・ジューコフスキー条件の出番となります。鋭く尖った翼後端であるにもかかわらず、その近傍で流れがなめらかに合流せず折れ曲がっているのは奇妙ですね？ そこで、参照している円筒翼周りの循環である$\Gamma$を色々な値に変えて、翼後端の流れがなめらかに見えるように調整します。

こんなところでしょうか。
上記のパラメータだと、$\Gamma =-32$位で筆者の目にはなめらかになっているように見えました。$b_0 = 1.0 - 0.4i$ということは迎え角26度くらいでまぁまぁ傾いてますね。ということで、26度くらいの迎え角で、速度に対して32(?)のレートで揚力を得られる翼型だとわかりました。この方法は、ジューコフスキー変換でなくても、クッタ・ジューコフスキー条件を運用可能な鋭い翼後端を与える正則関数であれば適用できます(逆関数を計算したくないけど)。

$c = 0.3 + 1.0i$の虚部$1.0$については、大きくすると翼が大きく湾曲します。実部は$\frac{a}{2}$に近づけると太くなり、空間を覆ってしまうようになります。逆にゼロに近づけると薄い翼型になります。

パラメータを色々に変えると、しばし極端な翼型とその流れが得られたりします。

// 循環:
const float gamma =-80.9;
// 円筒翼のローランb_0値：迎角と速度
const vec2 b_0 = vec2(1.0,-0.4);
// ジューコフスキー変換の元円盤の中心(xは正で２倍がaを超えないこと)
const vec2 inner_center = vec2(0.3,4.0);

こうした形については、クッタ・ジューコフスキー条件を満足しようとすると、翼後端が右に回り込んでいるため、循環がかなり大きくなります(上記は$\Gamma =-80.9$)。ということは高揚力ということですが、しかしこうした極端な翼型は実際には動作しないでしょう。おそらく、あまりに湾曲した翼型やキツすぎる迎え角では、流れをピッタリ翼型に張り付かせることは現実にはできず、そこは完全流体渦なし近似が破綻する領域になってくる だろうと思われます。そういう意味でこの二次元翼理論はかなり理想化されてはいますが、翼の性能において循環が重要な役割を果たすことを示すトイモデル といえましょう。

追記

上の可視化ではクッタ・ジューコフスキー条件を満たす$\Gamma$を目視で探しましたが、これを明示的に導出することも可能でした(@subarusatosi氏からの指摘・議論によります)。これを考えてみます。

ジューコフスキー翼の場合、翼後端は常に$\zeta = 2a$で、この逆像は重解であって$z = a$です。円筒翼の流れを引き戻したので、ジューコフスキー翼周りの複素速度は

$$x(\zeta) = w(J^{-1}(\zeta))\frac{dJ^{-1}}{d\zeta}(\zeta)$$

です。ところがジューコフスキー変換の定義から$\frac{dJ}{dz}(a) = 0$であるため、この近傍で同程度のオーダーで$w(a)=0$とならなければ速度場は特異的になります。クッタ・ジューコフスキー条件を満たさない流れの流線が滑らかでないのはこれが原因だと考えられます。この特異性を回避できている状況とは$w(a)=0$すなわち、 「翼後端の逆像が円筒翼表面の淀み点である」 ということです。この条件を考えましょう。

参照している円筒翼は、今原点にはなく、その大きさは反転円のそれであり、円筒翼中心は反転円の中心です。したがって、ジューコフスキー翼が参照している円筒翼の複素速度は

$$w(z) = b_0\left(1-\frac{\tilde{r}^2}{e^{2i\arg b_0}(z-\tilde{c})^2} \right)-\frac{i\Gamma}{2\pi}\frac{1}{z-\tilde{c}}\\ \tilde{r} = \frac{|c-a|a}{a-2\mathrm{Re}c}\\ \tilde{c} = a + \frac{(c-a)a}{a-2\mathrm{Re}c}$$

です。$z = a$にてこれがゼロである、つまり淀み点であるという条件を化すと、

$$\Gamma = -\frac{4a\pi}{a-2\mathrm{Re}c}\mathrm{Im}(b_0(c-a))$$

となります。

そこで、先のコードで

...
// 循環の理論値
float t_gamma() {
    float num = 2.0 * PI *a;
    float den = a - 2.0 * inner_center.x;
    float r = b_0.x * inner_center.y + b_0.y * (inner_center.x-a);
    return -(num / den) * 2.0 *r;
}
...

void main() {
	...
	float gamma = t_gamma();
	...
}

とすることで、常にクッタ・ジューコフスキー条件を満足した流線を描画できます。