「目的」
2次元翼理論にてよく出てくる複素速度ポテンシャル空間とジューコスフキー変換。概念が把握しずらい円中周りの循環Γの様子を詳しく観察する
「方針」
円柱周りの循環Γに関係する理論式を対象に、scilabにより複素計算を実行し、paraviewによりビジュアル化していきます。
１ 理論
１．１ 複素速度ポテンシャル空間(参考文献[1]4章)
複素速度ポテンシャルwは複素変数z(=x+iy)で微分すると実数部がvx、虚数部が-vyとなる性質があり、とても便利な関数となります。式で示しますと以下となります。
dw/dz=vx-ivy
そしてこの便利な複素速度ポテンシャルとして以下の3つが基本要素として教科書に出てきます。
A：一様流
w=Uz
U：一様流の速度
B：2重湧き出し
w=μ/z
μ：UR2
R：一様流との合体により速度「0」となる円半径
C：循環Γ(時計回り)
w=iΓ*logz/(2π)
Γ：4πUR
上記3つの関数の数値を計算し、ビジュアル化することでそのイメージを把握していきたいと考えます。
２．計算
２．１ 計算方法
これ以降検討する空間として以下としました(適当)。
・縦横-1～1[m]の正方形にて51×51の点群作成
この座標毎に複素速度ポテンシャルを微分することで速度vx、vyを求めます。
例えば一様流の場合、
dw/dz=d(Uz)/dz=U
これよりvx=U、vy＝0と簡単に計算出来ます。
しかしこれ以外は複素数計算となるため、これ以降はscialbにより計算しました(エクセルでも可)。プログラム骨子は以下です(2重湧き出しの場合を示します)。
z=x+%i*y;
VV＝-μ/z.^2;
Vx=real(VV);
Vy=-1*imag(VV);
更に下記paraview手順を実施することでビジュアル化が可能となります。
・paraviewにCSV(x,y,z,vx,vy,vz)として読み込む(z=vz=0)
・関数「table to point」でx、y、z座標を定義
・calculator「vx*iHat+vy*jHat+vz*kHat」で速度ベクトル生成
・Delaunay 2Dで各点を面で結ぶ
２．２ 複素速度ポテンシャルによる速度
（１）一様流
dw/dz=U
U：一様流の速度4[m/s]
予定通り一様流が表現出来ています。
（２）2重湧き出し
dw/dz=-μ/z2
μ：UR2＝4*0.52=1
R：0.5[m](適当)
教科書に出てくる蜘蛛のような形状が出てきました。
円頂上(R=0.5[m])での速度はvx=4[m/s]となっています(上部：時計回り、下部：半時計回り)。
（３）循環Γ
dw/dz=iΓ/(2πz)
Γ：4πUR＝8π
綺麗に循環Γにより各円周で一定速度の様子が観察できます。
円頂上(R=0.5[m])での速度はvx＝8[m/s]となっています(時計回り)。
（４）一様流+２重湧き出し
dw/dz=U-μ/z2
こちらも教科書に出てくる円を避ける一様流が観察されます。
円頂上(R=0.5[m])での速度はvx=8[m/s]で(1)と(2)の足し算になっています。
注目すべきは円の中にも流れが見られ、その流れが内部で完全に閉じている点です。
（本稿計算の動機はこの図でした。一様流に蜘蛛の形状をした2重湧き出しを足すことで円が形成されることがイメージ出来ないことが動機でした。）
（５）一様流+２重湧き出し+循環
dw/dz=U-μ/z2+iΓ/(2πz)
円頂上(R=0.5[m])での速度はvx=16[m/s]で(1)と(2)と(3)の足し算になっています。
逆に下方では速度が相殺され、そこが新しいよどみ点となっています。
尚円内は見たこともないような奇妙な図形となっています。
（６）一様流(迎角5度)+２重湧き出し+循環
前項までは一様流に迎角がついていませんでした。ここでは迎角を付けた場合の複素速度ポテンシャルによる速度の様子を観察してみます。
迎角αが点くことにより複素速度ポテンシャルは以下になるそうです。
w=U(e-iαz+R2eiα/z)+iΓ/(2π)*logz
但しΓ=4πURsinα
よって速度vxとvyはwを微分することで得られます。
dw/dz=U(e-iα-R2eiα/z2)+iΓ/(2πz)=vx-ivy
計算結果を以下に示します。
一様流はおよそ10度の傾き(5度ではありません)を持っており、これが循環Γによる吹きおろし(吹上げ？)を表現していると考えます(補足1)
念のため迎角60度の場合も計算してみました。 確かに循環Γの影響が出ています。よどみ点も90+60＝150度のところに発生しています。
「補足１」
一様流(迎角5度)+２重湧き出し+循環での複素ポテンシャルによる速度分布は以下(左図)でした。
そこで循環なしの場合を以下に示します(右図)。予想通り循環なしでは一様流が5度傾いた状態です。
ここから容易に推定されることとして迎角を90度に近づけていくと前後よどみ点が合流すると考えられます。実際の計算結果が以下で予想通り前後よどみ点が一か所に合流しました(全体が90度傾いています)。
ここで「一様流(迎角5度)+２重湧き出し+循環Γ」において発生する10度の偏向について観察してみます。
まずは前述「一様流(迎角5度)+２重湧き出し」(循環Γなし)での半径R＝0.5での周速分布を下記エクセルの図に示します。
上記流れの様子からは循環なしではθ=5度とθ=185によどみ点が出ていましたが、下記エクセルの図からも同じことが見て取れます。
ここに循環Γを足すわけですが、循環Γ=４πURsinαよりR＝0.5での周速vはv=Γ/2πR=2Usinαよりv=2*4*sin5=-0.697[m/s]となります。上図エクセルのA点(θ=0:最右端)ではvy=0.697であり、そこに周速vy=-0.697が加わるためvy=0となります。vx=0のためこの点が新しいよどみ点となるはずです。また上図エクセルB点(θ=190)ではv=-0.697となっており、そこに周速v=0.697が加わるためv＝0となり、この点も新しいよどみ点になると推察されます。
実際の計算結果「一様流(迎角5度)+２重湧き出し+循環Γ」を以下のエクセル図に示します。 予想通りθ=0とθ=190度のところによどみ点が移動していることを確認出来ました。
以上より循環Γを直接観察することは困難ですが、よどみ点の移動によりその存在を確認出来ました。
更にここで「一様流(迎角5度)+２重湧き出し+循環」の状況を揚力という観点で考察したいと考えます。
循環は揚力を発生する源として導入されています。その値は以下となります。
L=ρUΓ=4*2.19=8.76[N]
Γ=4πURsinα=8πsinα=2.19
一方円を囲む一様流は流管の方向が約10度(＝β)偏向しています。この偏向による作用力を計算してみます。流体の運動量保存則より
F=ρQUsinβ＝ρAU*vy
円形状前後のvyの様子は下図となります(横軸ｙ、縦軸vy)。
面積積分を行うことで上式におけるA*vyが計算されます。
面積計算よりA*vy=1.92(入口-出口)
∴F=UA*vy=4*1.92=7.67[N]
結果LとFはおおよそ一定しています。
この結果から循環Γは一様流の流れの向きを変えることで揚力を発生していることを確認出来ました(一般的な理解です)。
注目点として前述の複素速度ポテンシャル図から以下が推定されます。
・後縁では板に沿って流れが放出されるのに対し、前縁側で一様流よりも大きな迎角で流入してくるようです(別の言い方をすると変化は主に前縁で起きています)。
通常の感覚と違うので要注意です。但しこれはあくまで物理法則を考慮しない数学的解釈においてとなります。
３．ジューコスフキー変換
３．１ 理論
続いてジューコスフキー変換を観察してみます。
変換前座標z=x+iy、変換後座標Z=X+iYとしたとき以下の式がジューコスフキー変換となります。
Z=z+a2/z
この変換式でa=0.5のとき1[m]×1[m]、51×51点群(x,y))は以下の(X,Y)に変換されます。
色付き円[Ｒ=0.25(緑)、0.375(紫)、0.5(赤)、1(青)、1.25(黄)]の変換の様子を下図に示します。円が楕円形状に変換され、R=0.5(赤)では直線となっています。ジューコスフキー変換は複素ポテンシャル空間での循環理論を現実世界に転写する目的を持っており、ここでは円を平板(直線)に変換しています。
尚点が集中しているところが±1になります。
⇒
３．２ 複素速度ポテンシャルによる速度のジューコスフキー変換
ここからが本稿の本題で、平板板の周りに出来る速度分布を可視化してみます。
ジューコスフキー変換前の循環Γは綺麗な同心円流れ(同一半径上で一定値)を生成していました。この循環が変換後(実世界)でどのような挙動を示すのかを観察してみます。
・変換前複素ポテンシャルによる速度
w=U(e-iαz+R2eiα/z)+iΓ/(2π)*logz
dw/dz=U(e-iα-R2eiα/z2)+iΓ/(2πz)=vx-ivy
・変換後複素ポテンシャルによる速度
w=U(e-iαz+R2eiα/z)+iΓ/(2π)*logz
dW/dZ=dw/z*dz/dZ=[U(e-iα-R2eiα/z2)+iΓ/(2πz)]/[1-a2/z2]=Vx-iVy
(Vx、Vy：ジューコスフキー変換後空間での速度)
結果が以下です。循環Γがどこにいってしまったのかよくわかりません。
ここで前項で示したよどみB点190度がどこに変換されるかを確認します。
R＝0.5でのθ=190度の点はなんとX=-0.985に変換される結果となってしまいました。つまり複素速度ポテンシャル空間での10度は実空間では翼舷コード長の僅か0.75％(コード長では1/2倍になります)にしかならず、観測は事実上無理なレベルです(θ=185度はX＝-0.996に変換されました)。数学理論上出現する循環Γは実社会で計測することは非常に困難な結果となってしまいました。
2D実空間の計算結果を並べてみると複素ポテンシャル空間と実空間との類似点/相違点は以下となりました。
・前縁の速度増加領域(5.5[m/s]位のエリア)が対応(※1)。
・実空間でのよどみ点移動量はごくわずか(0.75%)のため、とても読み取れそうにありません。
・循環Γに対応する周運動はこの2D実空間の計算結果からは見て取れません。
（複素ポテンシャル空間ですら一様流と2重湧き出しを加えると打ち消しあってしまい、循環Γの姿を直接見て取れないため、実空間では見れないのは当然と言えます。）
※１：
ジューコスフキー変換後の薄板上部の増速部5.5[m/s](但しピーク値は18[m/s])と2D実空間での解析結果5.5m/sが一致する状況は類似する他形状でも起きており、偶然ではないと考えます。
一致する要因は以下と考えています。
・板形状では上面気流が早々に剥離し、円形状のような最大2倍の増速は期待出来ない。
・前述「一様流(迎角5度)+2重湧き出し+循環Γ」の様子で見た通り、気流変化は前縁が主で、かつ循環による誘導速度の2倍が作用する。
・ジューコスフキー変換前後で気流方向の空間は倍半分の変形を受けるが、前縁付近は実空間とちょうどバランスする領域が存在する。
以上から概算としては実空間での前縁速度増加分の半分が循環誘導速度と一致する見込みがあります。すると循環Γより以下の揚力計算が出来ます。
L=ρUΓ＝9.42[N](スパン方向1[m])
ρ：密度 今回は「1」とします。
U：一様流 4[m/s]
Γ：2πRｖ=2π*0.5[m]*(5.5-4)/2[m/s]=2.356
一方2D実空間計算結果CL=0.467(後述)より
L=1/2*ρU2*CL*ℓ=0.5*42*0.467*2[m]=7.47[N]
概算としては十分一致しています(CLが理論値に近づけば更に近い値になります)。
これは前縁流れの速度からおおよその揚力が簡単に計算出来ることを示しています。
[2D平板計算条件]
・解析空間 50×4×0.01[m]
・平面板形状 2×0.03[m]、端部R0.03[m]
(ジューコスフキー変換後の長さと一致させた)
・メッシュ
50×4×0.01[m]を100x8x1の合計800マスの2Dに分割。
平面板回りを2倍、4倍、8倍、16倍、32倍と段々に細分化。
・積分幅：adjustTimeStep=yesとしてクーラン数による自動設定。
・クーラン数 Co=0.5
・速度境界：入口4[m/s]、出口はinletOutlet。平面板表面はnoSlip、その他はslip。
・動粘度：空気 1e-5[m2/s]
・解析ソルバーはPimpleFoam(非圧縮・非定常乱流解析ソルバー)
・乱流モデルはRAS(kEpsilon)。
・k及びnutファイルの設定は適当(初期速度4[m/s]×1%で参考文献[2]4.5.5項に従う。)
・fvSchemes及びfvSolutionは適当。
・空力微係数計算結果：CL=0.467 / CD=0.0517 / L/D=9.0
４．まとめ
・複素ポテンシャル空間にて一様流と2重湧き出しにより生成される円の形状では外部と内部の流れが分離されることで円を表現していました。内部には幻想的な模様が生成されていました。
・揚力を生み出している循環Γは一定速度の円運動をしていますが、複素ポテンシャル空間ですら一様流と2重湧き出しを加えるとその姿を直接見て取れなくなっていました。当然実空間でも一定速度の周運動は観察できません。名残りとしてよどみ点が後方に移動するはずなのですが、あまりに僅かなためその存在を確認出来ませんでした。
参考文献
[1]BLUE BACKS 高校数学でわかる流体力学 竹内淳著 講談社
[2]OpenFoamによる熱移動と流れの数値解析 第2版 OpenFoamCAE学会編 森北出版
HA-technodefy 