線形安定性理論と固有値問題 - アーチェリー矢の側面境界層流れの不安定性と乱流遷移機構の解明

を得る．境界条件は，撹乱は表面および無限遠で消失するとして

u=0 (9.7)

とする．

さらに，撹乱uの振幅が微小である場合は式（9.5）の2次の微小非線形項(u· ∇)uを無視する．この場合、式（9.5）は線形撹乱方程式，

∂u

∂t +(U¯ · ∇)u+(u· ∇)U¯ = −∇p+ 1

Re∇²u (9.8)

となる．線形安定性理論において撹乱の初期振幅は微小であり、時間発展は式（9.6）-（9.8）で決定する．

9.1.2

線形撹乱のノーマルモード

撹乱は、次に示すnormalモードの組み合わせとして書ける．

u=v(x)exp(−iωt), p=q(x)exp(−iωt) (9.9) ここでxは一般座標を示し，ω^{は複素数である．式（}9.9）を式（9.8）、（9.6）に代入し，

iωv+(U¯ · ∇)v+(v· ∇)U¯ = −∇q+ 1

Re∇²v (9.10)

∇ ·v=0 (9.11)

である．物体表面と無限遠での境界条件は

v=0 (9.12)

である．したがって，固有関数の非自明解は式（9.10^）でω,v,q の固有値問題を解いて得られる．

複素数ω =ωr +iωi に対し，撹乱のノーマルモード（9.9^{）の時間因子は} exp(−iωt)= exp(ωit)

| {z }

gr owth/decay

(cosωrt−i sinωrt)

| {z }

oscillations

(9.13)

と書ける．ωr は撹乱の角周波数をωi は振幅の成長率を表す．ωi > 0^{の時は撹乱の振幅は} 時間とともに増幅し撹乱が成長することを意味する．一方、ωi < 0の時、撹乱は時間減衰するため基本流は線形安定である。一般にωi は撹乱の成長を表すため撹乱の成長率(growthrate) と呼ばれる．

る。基本流が並進対称性を持つ場合異なるFourier成分を独立に取り扱うことが可能である．

本研究において，線形安定性解析は円筒座標系を導入する．基本流U¯ ^はθ^{方向に対して不} 変であるため、撹乱は次の式のようにθ方向にスペクトル分解される。

« u_z u_r u_θ p

ª®®®

= ©

v_z(r,z) v_r(r,z) v_θ(r,z) q(r,z)

ª®®®

exp[i(mθ−ωt)] (9.15)

整数mは周回方向波数である．

平行流近似に基づく場合，基本流はz方向に対しても不変である．したがって，z方向に対してもスペクトル分解可能であり、

« u_z u_r u_θ p

ª®®®

= ©

« v_z(r) v_r(r) v_θ(r) q(r)

ª®®®

exp[i(αz+mθ−ωt)] (9.16)

のように記述できる。固有関数はr方向のみの関数となり，α,mは任意のパラメータである，

平行流近似に基づくためU_r = 0とする．式（9.16）を円筒座標系のNavier-Stokes方程式に加えて線形化すれば，円筒座標系の線形撹乱方程式を得る．

i(αU_z −ω)v_z+U_z^′v_r +iαq = 2 Re

[

v_z^′′+ 1 rv_z^′ −

(

α²+ m² r²

) v_z

]

(9.17) i(αU_z −ω)v_r +q^′= 2

Re [

v_r^′′+ 1 rv_r^′−

(

α²+ m²+1 r²

)

v_r − 2im r² v_θ

]

(9.18) i(αU_z −ω)v_θ + im

r q= 2 Re

[ v_θ^′′+ 1

rv_θ^′ − (

α²+ m²+1 r²

)

v_θ+ 2im r² v_r

]

(9.19) iαv_z +v_r^′ + 1

rv_r + im

r v_θ =0 (9.20)

ここで、^′は半径方向の微分を意味する．境界条件は，

v_z = v_r =v_θ =0, (r =0, r → ∞) (9.21)

である．

9.1.4

相似変数の導入と

Chebyshev

多項式による離散化

円筒座標系(r, θ,z)に対して以下の相似変数を導入する．

ζ =2 ( z

Re )¹₂

, σ = 2r−1

ζ r ∈ [0.5,∞) (9.22)

したがって，z,r 方向の偏微分は次のように変換される．

∂

∂z = ∂σ

∂z

∂

∂σ + ∂ζ

∂z

∂

∂ζ = 2 ζRe

( ∂

∂ζ − σ ζ

∂

∂σ )

(9.23)

∂

∂r = ∂σ

∂r

∂

∂σ + ∂ζ

∂r

∂

∂ζ = 2 ζ

∂

∂σ (9.24)

ここで，N−1次のChebyshev多項式に対して，N 個の極値をとる点（以下Chebyshev点

と書く）x_n は次の様に定義される．

x_n = cos ( nπ

N−1 )

0≤ n ≤ N −1 x_n ∈ [−1,1] (9.25)

σ ∈ [0,∞]^{であるため式（}9.26）により、Chebyshev点の区間[-1,1]と半無限区間[0,∞)^を対応させる．

σn = σ(ˆ 1+ x_n)

1−x_n σn ∈ [0,∞] (9.26)

x_n = −1→ σn =0，すなわちシャフト表面に対応し，無限遠では x_n = 1→ σn = ∞^に対応する．ここでσˆ は点の分配の疎密を決めるパラメータである．

任意の点σn における値 f(σn)^は，σ ^および f(σ)^{の値を使用して}spline^{補間により決定} する．

h_n(x)=

k=0,k,n x_n −x_k = δnk 0 ≤ n≤ N−1 (9.28) と書ける．式（9.27）をxで微分し，

dV(x) dx =

N−1∑

k=0

D_nkV^k, D_nk = dh_k(x)

dx (9.29)

とする．

2^階微分D⁽²⁾^，4^階微分D⁽⁴⁾^はD^{をそれぞれ}2^乗，4乗することで導出される．

D⁽_nk²⁾ = D_niD_ik (9.30)

D⁽⁴⁾_nk = D⁽²⁾_niD⁽²⁾_ik (9.31)

物体表面の圧力は定義されないため圧力項のChebyshev補間多項式Pの次数はN−1次となる．

P = ^N−2∑

k=0

P^khˆ_k(x) (9.32)

hˆ_n(x)= x_n+1

x+1 h_n(x) 0≤ n ≤ N −2 (9.33)

微分行列をDˆ として，

Dˆ_nk = x_k +1

x_n +1D_nk − 1

x_n+1δnk 0 ≤ n,k ≤ N−2 (9.34)

と定義される，

9.1.6 Chebyshev

選点法による離散化

速度場の撹乱vは無限遠と物体表面で0になるためChebyshev点の両端を省略しN−1点となる．したがって、微分はN−2次の多項式展開で表される．

dv dx

_x=x

N∑−2 k=1

D_nkv^k (9.35)

d²v dx²

x=xn

N∑−2 k=1

D²_nkv^k. (9.36)

圧力場の撹乱qに関しても無限遠に対応する点を省略するためN −3次の多項式展開，

N∑−2 k=1

q^kh˜_k(x) (9.37)

h˜_k(x)= x_n +1

x+1 · x_n −1

x−1 ·h_n(x) 1 ≤ n≤ N−2 (9.38)

で表される．よってその微分は

dq dx

x=xn

N−2

∑

k=1

D˜_nkq^k (9.39)

D˜_nk = 1

x_n²−1{(x²_k−1)D_nk −2x_nδnk} 1≤ n,k ≤ N−2 1≤ n,k ≤ N−2 (9.40) で表される．

したがって軸対称問題における線形撹乱方程式（9.17^）-^（9.20^{）は次式で書ける．}

i(α¯U_z −ω)¯ v_z +U_z^′v_r +i ¯αq = 1

Re_δ{v_z^′′+ξv_z^′− (α¯²+m²ξ²)v_z} (9.41) i(α¯U_z −ω)¯ v_r +q^′= 1

Re_δ[v_r^′′+ξv_r^′ − {α¯²+(m²+1)ξ²}v_r −2imξ²v_θ] (9.42) i(α¯U_z −ω)¯ v_θ +imξq= 1

Re_δ{v_θ^′′+ξv_θ^′ − (α¯²+(m²+1)ξ²)v_θ+2imξ²v_r} (9.43)

i ¯αv_z +v_r^′ +ξv_r +imξv_θ =0 (9.44)

ここで

Re_δ = ζRe/2, ξ= ζ

1+ζσ, (9.45)

α¯ = ζα, ω¯ = ζω. (9.46)

（9.49）で書ける．

ω¯AX¯ = BX¯ (9.49)

単位行列をI として、

A=©

« I

I I

0 ª®®®

, B= ©

b₀₀ b₀₁ 0 b₀₃ 0 b₁₁ b₁₂ b₁₃ 0 b₂₁ b₂₂ b₂₃ b₃₀ b₃₁ b₃₂ 0

ª®®®

. (9.50)

各小行列のn行k 列成分およびσ^{に関する微分行列を，}

(b₀₀)nk = α¯U_zⁿδnk + i Re_δ

(

D⁽²⁾_nk +ξnD⁽¹⁾_nk − (α¯²+m²ξn²)δnk

(b₀₁)nk = −iU_z^′ⁿδnk, (b₀₃)nk =αδ¯ nk, (b₁₁)nk = (b₂₂)nk = α¯U_zⁿδnk + i

Re_δ [

D⁽²⁾_nk +ξnD_nk⁽¹⁾ − (α¯²+(m²+1)ξn²)δnk

], (b₁₂)nk = 2ξn

Re_δξnδnk, (b₁₃)nk =−i ˜D_nk⁽¹⁾, (b₂₁)nk = − 2m

Re_δξn²δnk, (b₂₃)nk =mξnδnk,

(b₃₀)nk = αδ¯ nk, (b₃₁)nk =−i(D⁽¹⁾_nk +ξnδnk), (b₃₂)nk =mξδnk

D⁽_nk¹⁾ = 1

2 ˆσ(1−x_n)²D_nk, D⁽_nk²⁾ = 1

4 ˆσ²(1−x_n)³((1− x_n)D²_nk −2D_nk), D˜⁽¹⁾_nk = 1

2 ˆσ(1−x_n)²D˜_nk, D˜_nk = 1

x_n²−1((x_k²−1)D_nk −2x_nδnk),

と定義する．境界条件より1 ≤ n,k ≤ N −2の範囲で定義される．各小行列は N−2^次の正方行列であるから A,Bは4N−8次の正方行列となる．

9.1.7

圧力項の消去

圧力項を消去する．これにより，圧力に対する境界条件を削除でき，また行列のサイズを縮小し，固有値問題の計算時間を節約できる．式（9.50^）を，

ω¯v_z = b₀₀v_z+b₀₁v_r +b₀₃q (9.51)

ω¯v_r = b₁₁v_r +b₁₂v_θ+b₁₃q (9.52)

ω¯v_θ = b₂₁v_r+b₂₂v_θ +b₂₃q (9.53)

0= b₃₀v_z+b₃₁v_r+b₃₂v_θ (9.54)

と書き直す．v_z = (v_z¹, . . .,v_z^N⁻²)であり，他の成分についても同様である．式（9.54^）より，

係数行列を左からかけて，

ω¯b₃₀v_z = b₃₀b₀₀v_z+b₃₀b₀₁v_r +b₃₀b₀₃q, (9.55) ω¯b₃₁v_r = b₃₁b₁₁v_r +b₃₁b₁₂v_θ+b₃₁b₁₃q, (9.56) ω¯b₃₂v_θ = b₃₂b₂₁v_r +b₃₂b₂₂v_θ +b₃₂b₂₃q. (9.57) したがって，

q= P_zv_z +P_rv_r +P_θv_θ (9.58)

と書ける．ただし．

P_z = −(b₃₀b₀₃+b₃₁b₁₃+b₃₂b₂₃)⁻¹b₃₀b₀₀, (9.59) P_r =−(b₃₀b₀₃+b₃₁b₁₃+b₃₂b₂₃)⁻¹(b₃₀b₀₁+b₃₁b₁₁+b₃₂b₂₁), (9.60) P_θ = −(b₃₀b₀₃+b₃₁b₁₃+b₃₂b₂₃)⁻¹(b₃₁b₁₂+b₃₂b₂₂) (9.61) である．したがって，固有値問題は式（9.62^），^（9.63^{）のように}3N−6^{次となる．}

ω¯Y¯ = B¯Y¯, Y¯ = (v_z,v_r,v_θ) (9.62)

B¯ =©

b₀₀ b₀₁ 0 0 b₁₁ b₁₂ 0 b₂₁ b₂₂

ª®®

¬ +©

b₀₃P_z b₀₃P_r b₀₃P_θ b₁₃P_z b₁₃P_r b₁₃P_θ b₂₃P_z b₂₃P_r b₂₃P_θ

ª®®

(9.63)

9.1.8

線形安定性解析結果

線形安定性解析結果を図9.1^{に示す．〇は長谷川ら}^[17] ^の結果でRe =10000^，m= 1,2^，* はRe = 10000の非軸対称流のm= 0モードの流れにm= 1,2の撹乱を挿入して線形安定性

図9.1: 線形安定性解析結果，ωi −z依存性，Re= 10000，流線形鏃．

ドキュメント内アーチェリー矢の側面境界層流れの不安定性と乱流遷移機構の解明 (ページ 57-65)