自由流線理論による揚力計算

何となしに目にした1951年の古い論文に、面白い計算が書かれていた。

死水域を伴う任意翼型の特性
https://doi.org/10.2322/jjsass1948.4.30

一様流中に(任意の迎角を持って)置かれた平板翼で、後端と上面の途中に剥離点がある場合、揚力がどうなるか計算している。非粘性流体力学のバリエーションなので、剥離位置自体は、理論的に決めるわけではなく、外部パラメータとして与えられる。特に、自由流線理論に基づく計算に於いて、剥離点が限りなく後端に近い極限を取れば、後流が無限に薄い流れをモデル化できそうな気がする(これは、正しくないけど)

論文の計算方法は、今井功流体力学(前編)』(裳華房)の7章"不連続流"に説明されているのと同じ。この本も、初版は1973年とあるので、古い本ではある。不連続流の理論や、死水理論と呼ばれることもあるけど、自由流線理論(free streamline theory)という名前も使われるようなので、ここでは、この名前を使う(自由流線は、圧力が流線に沿って一定であることを指す;今井功流体力学(前編)』§51)

迎角を大きくしていく時、後端から前端へ剥離点が移動するパターンが見られるのは、後縁失速を起こす場合で、NACA4421翼型などで見られるらしい。自由流線理論によるモデル化が一番適切っぽいのは、このケースかもしれない。

平板翼では、前端で剥離した流れが、上面途中で再付着する剥離泡の形成が、迎角がかなり小さい時でも起きるらしいけど、自由流線理論で計算される流れでは、剥離泡は出ない(標準的な揚力計算でも出ないけど)(今井功流体力学(前編)』§34,Q3参照)。平板翼などで迎角を大きくしていくと、再付着点が後端に移動していき、やがて、前端で剥離した流れが再付着しなくなって揚力を失うらしい(薄翼失速と呼ばれる)

上論文中で扱われているNACA4412は、前縁失速と呼ばれるパターンを示し、迎角15度くらいで、突然、前縁で剥離するようになって揚力を失うっぽい。多分、現代のCFDなら、迎角と翼型(とレイノルズ数)を決めた時、剥離パターンや剥離位置(や再付着がある場合は再付着点)を決めることはできるかもしれないけど、解析的に計算できるケースは、ほぼないと思われる。

歴史的には、標準的な揚力計算の方は、平板翼や円弧翼だけでなく、実際の翼型で数値計算して、揚力係数のみでなく、翼面上の圧力分布を、定量的に、かなり再現することが確認されている。とはいえ、自由流線理論による計算でも、実際の翼型に対して、定量的によい結果を与える可能性は否定できない。

【翼型】翼型のデータは、検索すると拾って来れる。以下は、いくつか描画してみたもの。少なくとも、私には、これを見ても何もわからないが
f:id:m-a-o:20190921152626p:plain:w300f:id:m-a-o:20190921152637p:plain:w300f:id:m-a-o:20190921152646p:plain:w300

【補足:揚力の数値計算(CFD以前)】平板翼や円弧翼、Joukowski翼型などで揚力が解析的に決定された後、1920年代に、薄翼理論(thin airfoil theory,今井功流体力学(前編)』§49の内容)が作られた。真に任意翼型で適用できる計算法は、1930年代に開発された(今では、教科書で見ることもない。1960年代に現れたパネル法は、航空工学の教科書で説明されているのを見かける)。最初のものは、アメリカのTheodore Theodorsenが考案した方法で、1931年に、NACA technical report 411として、出版された。ここで、Clark Y翼型に対して、実験値との比較も行われている。また、NACA technical report 452では、NACA M6翼型に対して、同様の比較が行われている
[NACA-TR-411] Theory of wing sections of arbitrary shape
https://ntrs.nasa.gov/search.jsp?R=19930091485

[NACA-TR-452] General Potential Theory of Arbitrary Wing Section
https://ntrs.nasa.gov/search.jsp?R=19930091527

1933年に、Theodorsenの方法を、実際の翼型20個と、いくつかの異なる迎角に対して適用して、圧力分布を計算した報告が出ている。
[NACA-TR-465] Determination of the theoretical pressure distribution for twenty airfoils
https://ntrs.nasa.gov/search.jsp?R=19930091540

任意形状翼型に適用できる別の方法として、1938年には、守屋富次郎という人が、別の方法を提案している。こっちは、Fourier級数の形で与えた翼型形状に対して、速度ポテンシャルの形をFourier係数を使って具体的に書き、循環を計算している。このような手法は、層流翼の開発で役立ったらしい
任意の翼型の特性を求める一つの方法
https://doi.org/10.14822/jjsass1934.5.33_7


【補足:"Kuttaの条件"の呼称】Kuttaの条件は、後端に於ける流速が有限になる(後縁角が正の場合は0になる)という境界条件の一種で、今井功流体力学(前編)』§33では、"Kuttaの条件あるいはJoukowskiの仮定"と書かれている。ランダウ・リフシッツの教科書では、Zhukovskii-Chaplyginの条件と呼ばれている。誰が何を考えたのか、原論文を参照できなかったけど、以下の論説に、事情が書かれている。

クッタとジュコフスキーの翼理論
https://www.jstage.jst.go.jp/article/nagare1970/5/4/5_4_1/_article/-char/ja/

まず、1902年に、Kuttaが(迎角0度の)円弧翼に働く揚力を決定した(リリエンタールの測定結果に影響を受けたと言われるが、リリエンタールの測定とは有限翼幅補正を入れても定量的には合わないので、この不一致をどう考えたのかは分からない。とにかく、迎角0度でも揚力が働き、流速の二乗に比例するという結果は得られた)。従って、Kuttaは"Kuttaの条件"そのものか、それと同値な仮定をしたのだろう。この仕事は、特に知られることもなく、1906年に、Joukowskiが、いわゆるKutta-Joukowskiの定理を一般的な状況下で示した。けれど、Joukowskiは、循環の値を決める方法は、与えなかったようである。上論説中には『これに対して後縁の条件は,Kuttaひとりに帰すべきものである。少なくとも1910年以前のZhukovskiの論文には,この条件に言及したものがないのである。したがって,"Kuttaの条件"と呼ぶのが妥当であって,Zhukovskiの名を加える必要はない』と書いている。

少し後の方には、『ここで付け加えておきたいのは,同じ1910年に発表されたロシアのChaplyginの論文である。ChaplyginはKuttaの発表(1910)を知らないで,まず迎角のある円弧翼の揚力を計算し,さらに二つの円弧から成る三日月翼断面その他の場合を取扱っている』ともあり、Chaplyginは、循環の値を決定しているようなので、Kuttaの条件に相当する条件を置いたのだろう。なので、Kutta-Chaplyginの条件と呼ぶのが妥当かもしれない。そう呼ぶ人もいないわけではないっぽい。

Chaplyginの1910年の論文というのは、以下で読めるっぽいけど、(ロシア語が)読めない
О давлении плоскопараллельного потока на преграждающие тела (к теории аэроплана)
On the pressure of plane-parallel flow on the blocking bodies (to the theory of an airplane)
http://gpntb.dlibrary.org/ru/nodes/5240-chaplygin-s-a-o-davlenii-ploskoparallelnogo-potoka-na-pregrazhdayuschie-tela-k-teorii-aeroplana-m-1910

20世紀初頭の物理学者にとっては、複素解析や等角写像の理論は、まだ標準的なレパートリーになってなかったかもしれない。リーマンの写像定理が1851年で、Christoffelが、Schwartz-Christoffel変換の論文を書いたのが1867年らしい。流体力学自体は古い分野だけど、必要な数学ができて、計算できるようになったという側面もあったのだろう



自由流線理論の歴史については、

Ideal free streamline flow over a curved obstacle
https://www.sciencedirect.com/science/article/pii/0377042795002723

のIntroductionに、少し書かれているけど、適当に補うと、以下のような感じ。

流れの中に物体があるばあい、物体が流体から受ける力を計算するには、物体の表面に働く圧力を求め積分すればいい(今井功流体力学(前編)』§30参照。翼に働く表面摩擦の場合は、せん断応力を積分しないといけないけど)。このことは、18世紀には既に分かっていて、この方針で計算した結果、ダランベールパラドックスが導かれた(今井功流体力学(前編)』§28参照)。

実験的には、(水中でも空気中でも)抗力が流速の二乗に比例するという結果は何度も確認されたし、ダランベール自身も、水中で平板に働く抗力を測定している。ダランベールが理論的に考えたような流れが生じるのは、レイノルズ数が低い時のみ(今井功流体力学(前編)』§2、1-1図)で、一方、実測したような条件下では、後流ができ、抗力の主要な部分は圧力抵抗に由来する(低レイノルズ数で正しいStokes抵抗則は1851年に発表されている)

自由流線理論の開拓者は、Helmholtz(1868)らしい。今井功流体力学(前編)』§52に、噴流への適用が書かれていて、このへんの仕事っぽい(?)。"free streamline"という名前は、Kirchhoffによるとされる。


1876年に、Lord Rayleighは、平板の後流部を死水域と仮定して、一様流に対して(大きな)迎角を持つ平板に働く抗力を計算した。但し、剥離は、平板の前端と後端で起こると仮定している(この流れは、Helmholtz流と呼ばれることがある)。この結果は、今井功流体力学(前編)』§53に説明されている。論文は、以下にある

On the resistance of fluids
https://www.tandfonline.com/doi/abs/10.1080/14786447608639132?journalCode=tphm16

抗力の計算値は、測定値の半分弱でしかなかったけど、(1+4/π)倍すると、1798年のVinceの測定結果(水中の測定らしい)とよく一致すると書いてある。この謎の定数倍因子について、論文では"The fourth column represents the law of resistance according to the formula now proposed, a factor being introduced so as to make the maximum value unity"と書かれている。補正がなければ、流れに対して垂直に立てた平板の抗力係数は、C_D = 2\pi/(4 + \pi) \approx 0.88で、補正因子を掛けると、C_D=2.0を得る。こうして、(ニュートン流近似は別として)非粘性流体力学の範囲で、初めて、抗力が流速の二乗に比例するという結果が得られた

補正因子が必要になるのは、死水域の圧力を、静圧に等しいとしたためで、実際は負圧になっているとされる。適当な背圧を外部パラメータとして与えてやれば、正しい結果になるというのが、論文で与えられた謎の補正因子の現代的解釈だと思う。今井功流体力学(前編)』$55では、背圧を無次元の形で表した空洞係数(cavitation number)を導入して(Q=0は、背圧が静圧と等しい場合で、上記の論文の補正は、Q=4/πに相当する)、"C_D(Q)=C_D(0)(1+Q)によって、物体の受ける抗力を近似できると考えられるが、一般的な証明は、いまのところ与えられていない"と書かれている。

こうして、外部パラメータが必要ではあるけど、一様流中の平板に働く(迎角が大きい時の)圧力抵抗が計算できるようになった。同じ方法で、揚力も計算できるけど、論文には何も書かれていない。前縁で剥離が起こるようになって以後、大迎角に対する揚力係数の実験値を、見かけることが、あまりない気がするのだけど、

Aerodynamic characteristics of seven symmetrical airfoil sections through 180-degree angle of attack for use in aerodynamic analysis of vertical axis wind turbines
https://doi.org/10.2172/6548367

に、いくつかの対称翼に対して、迎角、レイノルズ数ごとのデータが記載されている。迎角と揚力係数の関係を、この論文の実験データの一つと、Helmholtz流(空洞係数の補正あり)に対するC_L=(1+4/\pi)\dfrac{2\pi\sin(\alpha)\cos(\alpha)}{4+\pi\sin(\alpha)}、Kutta-Jouwkoski流に対するC_L = 2\pi \sin(\alpha)(今井功流体力学(前編)』式(33.15))を、プロットしてみる。対称翼なので、負の迎角は、見る必要がない。

f:id:m-a-o:20190824140049p:plain
迎角と揚力係数

元々のダランベールパラドックスは、円柱に対する計算だったけど、Rayleighは、円柱に自由流線理論を適用した計算もしてない。代わりに(?)、1879年の論文で、Magnus効果(1852年、ドイツのHeinrich Magnusが報告していた)について考察し、揚力の計算を試みている。論文では、円柱に対するKutta-Joukowskiの定理が書かれているけど、循環の具体的な値は決めてない(ので定量的な予測はされてない)。

On the Irregular Flight of a Tennis-Ball
https://doi.org/10.1017/CBO9780511703966.054

この説明は、今でもよく見かけるけど、飛翔するテニスボールを、このようにモデル化するのは(多分)正しくない。飛翔するテニスボールの場合、レイノルズ数は、10^5オーダー(空気の動粘性係数の逆数くらい)で、ボール自身が回転するので、厄介であるけど、無回転の時と同様、背後に、後流が発生しているらしい。この後流領域を、死水領域としてモデル化することが妥当かは知らないけど、どっちにしろ、何も断らずにKutta-Joukowskiの定理を適用していい問題ではないと思う

cf)マグヌス効果の物理的メカニズムについて
https://www.jstage.jst.go.jp/article/jjsass/57/667/57_667_309/_article/-char/ja

【Magnus効果の応用】1920年代に、ドイツのAnton Flettnerという人が、風を受ける回転円筒を備え、Magnus効果で推進するローター船を作ったとされている。回転する円筒で飛行する航空機(Rotor airplane/Flettner airplane)も、試作されたことはあるようだけど、あまり記録がない。Flettnerは、その後、初期のヘリコプター開発にも関わったようだけど、Flettnerが、Flettner airplaneを作ったということは、ないっぽい。ローター船の場合、風速は、数m/sだろうけど、直径は約4メートルだったらしいので、レイノルズ数は、テニスボールの場合と同程度か、それより大きいと考えられるので、単純にKutta-Joukowskiの定理で、推力を説明するのは、正しくない結果を与えるだろうと思われる



特に重要な進展があったわけではないようだけど、1890年に、Joukowskiも自由流線理論に関する論文を書いている。

1907年には、イタリアのLevi-Civitaが、平板以外のcurved obstalceに適用可能な一般論を展開したとされている(多分、今井功流体力学(前編)』§54にあるような内容だと思うけど嘘かもしれない)。Levi-Civitaの論文は、以下で読めるっぽいけど、(イタリア語が)読めない。Levi-Civitaといえば、「接続」の枕詞くらいの認識だけど、こういうこともやっていたらしい

Scie e leggi di resistenza [英訳:Trails and laws of resistance]
http://doi.org/10.1007/BF03013504



1923年になって、イギリスのBrodetskyが、円柱に対して、自由流線理論を適用した計算を発表した。

Discontinuous fluid motion past circular and elliptic cylinders
https://royalsocietypublishing.org/doi/10.1098/rspa.1923.0014

また、ドイツのCurt Schmiedenは、1929年、1932年、1934年に、同じく円柱に自由流線理論を適用して、論文を書いている。例によって、ドイツ語なので読めない
Die unstetige Strömung um einen Kreiszylinder
https://link.springer.com/article/10.1007/BF02079711

Über die Eindeutigkeit der Lösungen in der Theorie der unstetigen Strömungen
https://link.springer.com/article/10.1007%2FBF02080809

円柱の場合は、一般に解析的に解くことはできないけど、複素対数速度が、単位円板上で正則なので、そこでTaylor展開することができ、係数を数値的に決定できる。Schmiedenが得た面白い結論の一つとして、物理的に許される剥離点の剥離角(円柱前部のよどみ点から剥離点までの角度を円柱中心周りに測った大きさ)は、概ね55度〜125度の範囲にあるというのがある(らしい)。これらの角度の範囲外では、2つの自由流線が交わるということが起きる。実際の流れにおいて、どこで剥離するかは、Reynolds数に依存して異なるが、可能な流れは、ポテンシャル流の計算で、ある程度制限される。

これは、定常性を仮定した議論なので、剥離点の位置が時間的に変動する場合は、これらの限界を超える可能性はある。実験的に剥離位置を調べた結果は、
円柱を過ぎる流れの剥離
https://www.jstage.jst.go.jp/article/jjsass1969/20/226/20_226_610/_article/-char/ja/
の第10図に、レイノルズ数の範囲6.0e4〜5.0e6での剥離角の測定結果が載っている。レイノルズ数が3.5e5〜1.5e6では、剥離角が140度強に達しているけど、それより低い領域では80度前後、高い領域では120度程度と、Schmiedenの計算値内に留まっている。レイノルズ数が3.5e5〜1.5e6では、円柱に剥離泡が付着する(第11図の超臨界領域)らしく、そもそも、流れが、自由流線理論でモデル化されてるものとは異なる。


1933年に、ChaplyginとLaurentievは、平板上面真ん中で剥離する流れに働く揚力と抗力を計算したらしいが、論文が、どれか分からない。冒頭の論文の計算に近いものか、実質的に同じものかもしれない。冒頭の論文と同じ計算は、Schmiedenが1941年に行っている。

[NACAーTM-961] Flow around wings accompanied by separatioai of vortices
https://ntrs.nasa.gov/search.jsp?R=19930094455

Schmiedenの論文と冒頭の論文は、複素ポテンシャルの取り方に定数差があるけど、それ以外は、同じ。今井功流体力学(前編)』式(54.1)に合わせて、定数項を0とするなら、
w = \dfrac{C}{2}(\zeta^{2} + \zeta^{-2}) - 2C \cos(\theta_S) (\zeta+\zeta^{-1})
となる(\zetaは、Schmiedenの論文では\tauになっている。それ以外の記号は、論文通り。冒頭の論文の記号では\theta_S=\sigma_0で、C=\dfrac{a^2}{2})。

そもそも(私の興味)は、Helmholtz流とKutta-Joukowski流を繋ぐ解があるなら面白いと思ったのだけど、残念ながら、そうはなってない。Schmiedenは、2つの極限状況を扱っていて、一つは、Helmholtz流になり、Helmholtz流から、前縁の剥離点を後端へ向かってう動かしていく途中で、2つの自由流線の交差が起きるようになるっぽく、それを、もうひとつの極限的ケースと位置づけている。つまり、迎角の大きさに応じて、上面途中の剥離点が動ける範囲が制限を受け、(迎角0度の場合を除いて)後端まで動かすことはできない。

数学的には、Schmiedenの論文に於けるk^2は、0を超えて、-1まで動くことができ、冒頭の論文では、「ただあまり剥離点を後緑近くにもつていくと,上下の剥離流線が交るようになると思われるが,そのときにもこの理論はなお近似的な意味は持ち得るであろう」として、さらっと、限界を超えているけど、Schmiedenは、自由流線が交差するべきでないという考えから、k^2は正としている。名前がないと不便なので、k^2=0の時の流れを、Schmieden流と呼んでおく

Schmieden流の場合は、(冒頭の論文では、より簡潔に、Poisson-Schwartzの積分公式で計算されている通り)
\dfrac{dz}{dw} = \dfrac{ 1+ \zeta^2 e^{-\alpha}}{1+\zeta^2 e^{i\alpha}}
なので、dz/d\zetaを決定できる。積分すると、z=z(\zeta)が求まる。Schmieden流に限らない一般の場合も、解析的に計算することができる。ついでに、Schmiedenの論文の式(3)の分母内に符号の誤植がある

Schmieden流では、迎角α>0の時には、上面途中で、剥離が生じている。迎角が10度の時、剥離位置は、後端から、L/5ほど進んだ所(Lは平板の長さ)にある。Helmholtz流でも揚力はあったけど、Schmieden流の場合は、標準的な揚力計算の結果と類似した力が出る。平板翼のSchmieden流に対する揚力係数は、迎角をαとした時
C_L = 2 \pi \sin(\alpha) \cdot f(\alpha)
の形で書ける。αが小さい時、f(α)は1より少しだけ小さい。αが10度の時、f(α)は、0.908で、通常の揚力計算によるより10%ほど小さい。

圧力抵抗は、複素対数速度をTaylor展開した一次の項の係数に比例する(今井功流体力学(前編)』式(54.13)にある通り)ので、Schmieden流の圧力抗力係数は、剥離があるにも関わらず、0となる。


Schmiedenの論文冒頭には、以下のように書かれている。

The flow around wings computed by the usual method leads in the case of a finite trailing edge to a stagnation point in the trailing edge due to the Kutta-Joukowsky condition of flow governing this region. As a result, the theoretical pressure distribution differs substantially from the experimental values in the vicinity of the trailing edge.

Kutta条件を課して計算すると、後縁角が正の場合、後端の流速は0なので、後端近傍の圧力分布は、計算値と実験値が大きくずれると言っている。上の方であげた守屋富次郎の報告中の第10図に、NACA-4412翼型(迎角6.4度)に対する実験値と計算値が比較してあって、後端近くまで実験値がプロットしてあるけど、確かに合ってない

自由流線理論ならいいかというと、自由流線理論では、デフォルトの(Schmiedenの計算でも使われている)Kirchhoffの死水モデル(今井功流体力学(前編)』§54)に従う限り、(平板翼に限らず)後端の流速は一様流速に等しくなるはずだけど、NACA-4412翼型の実験値を見ると、これも成立していなそう。これは、背圧を静圧と等しいとしているせいで、cavitation numberを決めれば解決するかもしれない

レイノルズ数に於ける、平板に対する揚力係数の実験値が、意外と見つからなくて困るのだけど、最初の方に書いたように、かなり小さい迎え角でも剥離泡が形成されて、理論値とずれ始めるらしい(要出典)


とりあえず、Schmiedenの計算と標準的な揚力計算で、平板翼に対して、揚力係数が近いだけでなく、翼表面の流速分布が近いことも確認しておく。自由流線理論でも、死水域の外では、流速分布からベルヌーイの定理で、圧力分布を決めていい

標準的な揚力計算の結果は、今井功流体力学(前編)』式(33.1)と(33.9)を使い(a=1/2とした)、自由流線理論の結果は、Schmiedenの論文の式(10)及び(12)を使う(zの積分定数は、zの範囲が、-1から+1になるように決定した)

横軸は、平板の位置座標(-1から+1の範囲。普通は、0〜1になるように取る)で、縦軸は、平板上の流速を一様流速で割った無次元量(守屋富次郎の報告内の式(14)で、\dfrac{v}{V}と書いてある量。裏と表の速度2つがある)。自由流線理論の計算は、途中で剥離が生じ、一部は、死水域に入っているので、そこの流速分布は決まらない

以下は迎角10度の時のグラフで、青が標準的な揚力計算から出る流速分布で、赤は自由流線理論から出る流速分布。異なる考えから出発して、式が全然違うにも関わらず、よく一致している。迎角を小さくすると、両者の区別が殆どつかなくなるので、多少は、差が見える値として、この迎角を選んだ。

f:id:m-a-o:20190824140147p:plain
Kutta-Joukowski流とSchmieden流の流速分布

自由流線理論の論文は、その後も出ているけど、実際の翼型に適用したという類の話は少ないように見える。

任意物体のまわりの自由流線理論
https://repository.exst.jaxa.jp/dspace/handle/a-is/24280

では、自由流線理論の新しいモデル(自由流線上で、流速は一定でなく、無限遠で一様流速に等しくなる。もはや定義上の"自由流線"じゃない)を提案し、数値計算法を与えている。シリンダーと(比較的迎角の大きめな)NACA4412翼型で、圧力分布を再現できたと書いているけど、剥離位置と背圧を、圧力分布から推定したと書いてあって、ただのパラメータフィッティングなので微妙。


自由流線理論を、現在の標準的な流体力学のカリキュラムで習うことは、まずないだろうけど、翼の揚力計算について無力で誤った雑魚理論というわけではない