アーチェリー矢の側面境界層流れの不安定性と乱流遷移機構の解明

  電気通信大学大学院  

情報理工学研究科 機械知能システム学専攻

1932022

小野口 貴仁

指導教官：宮嵜 武 教授

副指導教官： 守

裕也准教授

令和

3

年

1

月

25

日

第 1 章 緒言 1 第 2 章 無次元数の定義 5 2.1 使用する無次元量・用語の定義 . . . 5 第 3 章 基本流の数値計算 6 3.1 軸対称基本流の数値計算 . . . 7 3.2 非軸対称基本流の数値計算 . . . 13 第 4 章 ２次不安定性解析 18 4.1 ２次不安定性 . . . 18 4.2 基礎方程式. . . 18 4.3 固有関数の周回方向の離散化 . . . 20 4.4 Chebyshev選点法による半径方向の離散化 . . . 21 4.5 固有値問題への帰着 . . . 22 第 5 章 基本流の設定 24 5.1 有限振幅撹乱の挿入 . . . 24 5.2 非軸対称流のフーリエ変換 . . . 25 5.3 流れの重ね合わせ . . . 28 5.4 主流の補間. . . 28 5.5 基本流の比較 . . . 29 第 6 章 解析結果 31

6.1 収束の確認. . . 31 6.2 非軸対称流のm = 0モードの流れに同m = 1モードを挿入した流れに対す る2次不安定性 . . . 34 第 7 章 実験結果との比較と議論 44 第 8 章 結言 47 第 9 章 付録 49 9.1 線形安定性理論と固有値問題 . . . 49 9.2 Bi-global安定性解析 . . . 57 付録 A Bi-global 安定性解析 62 A.1 基礎方程式. . . 62 A.2 主流の離散化 . . . 64 A.3 固有値問題. . . 64 A.4 妥当性の検証方法 . . . 65 参考文献 67  

3.1 鏃形状. . . 6 3.2 鏃（実寸mm）のCAD図，上: 流線形鏃，下: 椎型鏃．長谷川[30]_{より再生． 7} 3.3 流線形鏃を装着したベアーシャフトの計算格子，長谷川[30]より再生． . . . 9 3.4 主流の速度分布. Red=2.0 × 104,SP= 0.015, Uz. . . 10 3.5 鏃ごとの主流速度分布Uz. Red = 2.0 × 104, 1.0 × 104, SP = 0.025.左から z = 0, 50, 100, 150, 200． . . . 11 3.6 周回方向速度分布 U_θ. Red = 2.0 × 104, 1.0 × 104, SP = 0.025. 左から z = 0, z = 50, z = 100, z = 150, z = 200. . . 12 3.7 周回方向速度分布U_θ. Red = 2.0 × 104. SP= 0.015, 0.025, 0.05,左からz=0, z=100, z=200. . . . 13 3.8 迎角0.5 度がついた矢周りの直接数値計算結果のz 方向速度分布の全体図， (a)椎型鏃，(b)流線形鏃． . . . 15 3.9 迎角 0.5 度がついた矢周りの直接数値計算結果の z 方向速度分布，Re = 15000，流線形鏃． . . . 17 3.10 迎角 0.5 度がついた矢周りの直接数値計算結果の _θ 方向速度分布，Re = 15000，流線形鏃． . . . 17 5.1 upz の空間構造，Re= 20000, z = 5 . . . 24 5.2 主流のフーリエ変換，z方向速度，Re=15000，N=64，流線形鏃． . . . 26 5.3 主流のフーリエ変換，θ方向速度，Re=15000，N=64，流線形鏃． . . . 27 5.4 軸対称流れの補間および微分結果，Re=20000，N=64，流線形鏃． . . . 28 5.5 m= 0モードの流れの補間および微分結果，Re=15000，N=64，流線形鏃． . 29

5.6 軸対称流の比較と非軸対称流のm = 0モード流れ． . . . 30 6.1 固有関数，Re = 15000, z = 50，α = 1，N = 64，M = 9，流線形鏃，迎角0.5 度． . . . 33 6.2 固有値分布，Re = 15000，z= 20，α = 1.7，N = 64，迎角0.5度． . . . 34 6.3 α − ωi 依存性，Re= 15000，N = 64，流線形鏃，迎角0.5度．. . . 34 6.4 z− ωi 依存性，N = 64，迎角0.5度 . . . 35 6.5 固有関数，Re = 15000, z = 2，α = 3.4，N = 64，ω_i = −0.039765453，流線 形鏃，迎角0.5度． . . . 36 6.6 z− ωi 依存性．青，赤，黄色の実線は，m= 0モードの流れにm = 1モードの 流れをそれぞれ1倍，1.2倍，1.6倍して挿入した結果を表す(Re= 15000，流 線形鏃)．赤線に黒丸，緑線に黒丸は伊藤ら[22]_{らの結果を表す(Re = 20000，} 椎型鏃)． . . . 38 6.7 ωi − α依存性，Re= 15000，流線形鏃，迎角0.5度 . . . 39 6.8 固有関数，Re = 15000, z = 2，α = 3.8，N = 64，ωi = −0.035940，A= 1.6， 流線形鏃，迎角0.5度． . . . 40 6.9 固有関数，Re= 15000, z = 50，α = 1.7，N = 64，ω_i = 0.029636，A= 1.6， 流線形鏃，迎角0.5度． . . . 41 6.10 固有関数，Re= 15000, z = 68，α = 2.4，N = 64，ω_i = 0.048560，A= 1.6， 流線形鏃，迎角0.5度． . . . 42 6.11 ωi − α依存性，流線形鏃 Re = 20000の軸対称流に流線形鏃 Re = 15000の 非軸対称流のm= 1モードの流れを挿入した． . . . 43 7.1 CD 値の Re数依存性, A/C/E矢，椎型鏃．[30]より引用． . . . 44 7.2 CD のγ 依存性，Re = 10000.[30] より引用． . . . 45 7.3 矢にかかる加速度の時間変化，初速 56.4[m/s]，A/C/E 矢，椎型鏃．Ortiz et al.[29] より引用． . . . 46 7.4 矢にかかる加速度の時間変化，初速60.0[m/s]，A/C/E矢，流線形鏃．大場 [31]_{より引用． . . . . 46} 9.1 線形安定性解析結果，ωi − z依存性，Re= 10000，流線形鏃． . . . 57

表目次

矢の歴史は長く，古くは狩猟や戦争の道具として用いられ，現代ではアーチェリーや弓道 などの競技として幅広い年代に親しまれている．アーチェリー競技では新素材や新機構の導 入により弓具の性能が飛躍的に向上している．それに伴い，勝利者の得点は年々向上してお り，競技で求められる矢の命中精度も年々高まっている．オリンピックのアーチェリー競技 では 70 m先の直径12.2 cmの円を常に射抜くことが要求される。競技者の放つ矢の速度は およそ60m/sで、矢の直径を代表長としたRe数は約2.0 × 104に達する．飛翔中の矢に働く 空気抵抗は重力の0.5倍を超えるが，その半分以上はシャフトに働く摩擦抵抗からの寄与で あり、残りは矢羽に働く空気抵抗である．矢側面の境界層流れが層流である場合と乱流であ る場合では空気抵抗が有意に異なり，鉛直方向に20 cm程度のずれが生じるため，乱流遷移 の詳細なメカニズムを解明し，それを制御することはアーチェリー競技のさらなる発展のた めに重要である． Park et al.[1]は水槽実験を行い，Re = 2.74 × 104（代表長さは矢の直径）において先端から 乱流境界層になることを示した。 Denny[2]は矢の形状から推定される抗力係数や揚力係数といったを求め，飛翔中の矢の軌 道を数値的に算出した．Zanevskyy[3]_{は、矢の姿勢と迎角の関係を取り入れ，矢の空力特性} を算出する手法を提案した．Liston[4]，Park[5]は矢を構成している鏃・シャフト・矢羽・ノッ クにおけるそれぞれの抵抗の総和が矢の抵抗であるとし，矢の側面境界層流れが層流または 乱流であるかにより抗力係数が有意に異なることを指摘した． 鈴木ら[6] _{は日本航空宇宙研究機構（JAXA）の磁力支持天秤装置付き風洞（MSBS風洞）} を使用し，クロスボウ（細長比57）の矢の抗力係数CD、揚力係数CL、ピッチングモーメン

ト係数CM を測定した．測定が行われた104 ≤ Re ≤ 5 × 104の範囲では，境界層の状態は鏃

形状に依存することが示された．流線形鏃の場合では層流境界層となり凹凸型鏃を装着した 場合，乱流境界層となった。椎型鏃の場合，Re = 9.8 × 103付近で境界層は乱流に遷移した．

アーチェリー競技者が使用するEaston社製のA/C/E矢の空力特性はMiyazaki et al.[7, 8]の

MSBS風洞実験と飛翔実験といった，支持干渉のない 2種類の実験により詳細に調べられ た．MSBS風洞実験では鏃形状によらず矢を気流に対して平行に置く場合に，Re ≤ 1.3 × 104 で境界層は層流となる．一方，椎型鏃を装着した矢を気流に対して，0.75◦傾けた場合には Re= 1.1 × 104付近で乱流遷移することが示された．また，圧縮空気によって矢を発射する飛 翔実験結果は，椎型鏃を装着した場合は1.3 × 104 < Re < 1.8 × 104の領域でCD が層流境界 層に相当する値（CD ∼ 1.5）と乱流境界層に相当する値（CD ∼ 2.7）に二値化し，流線形鏃 の場合ではこの乱流遷移領域がより高Re 数領域（1.3 × 104 ≤ Re ≤ 2.4 × 104）に広がった． Miyazaki et al. [8]はCD の二値化の原因が矢の迎角のばらつきにあることを指摘し，矢の迎 角が気流に対して_±0.6度以上つくと境界層が乱流遷移すると指摘した．実際，迎角がつかな いように放たれた椎型鏃の矢のCD は，遷移領域のRe = 1.7 × 104でも層流値（CD = 1.34） となることが確認されている．矢周りの流れが乱流遷移する機構解明のためには，これらの 結果を理論的に説明することが不可欠である． 乱流遷移の機構の1つには，層流の時間的に定常な流れが微小撹乱に対し不安定になり， その微小撹乱が時間的に発展すると遷移が始まるものがある．撹乱が発達する場合は，その 流れは不安定性であり，そうでない場合は安定な流れと呼ぶ． アーチェリー矢のような円筒側面に発達する境界層流れの安定性の研究はいくつか存在す る．Seban&Bond[9]やKelly[10]，Glauert&Lighthill[11]は細長い円柱を過ぎる軸対称境界層流

れの境界層方程式を導出し、べき級数展開により近似解を求めた． Tutty et al.[12] は非回転円柱の側面境界層流れについての数値計算を行い，その線形安定性 解析を平行流近似に基づいて行った．最初に線形成長する撹乱は周回方向波数m= 1のモー ドであり，臨界Re 数は2120であることが示された．次いで，m = 2, 3, 0, 4の順に不安定と なり臨界 Re 数はそれぞれ12140, 20204, 24878であることが示された．Herrada et al.[13] は 非回転および回転円柱の側面境界層流れの線形安定性を調べた．その結果，非平行流効果に より臨界Re数はわずかに低下し，回転効果により大きく低下することを示した．Muralidhar et al.[14, 15] は軸周りに高速回転（無次元回転数Sp = 0.1）する円柱の側面境界層流れの線形 不安定性解析を平行流近似に基づいて行なった．線形化 Navier-Stokes方程式の離散化には

アーチェリー矢は細長比はO(102) 程度であり先端には鏃が装着されている．そのため， アーチェリー矢の側面境界層流れの安定性解析は，鏃の影響も加味した基本流に対して行う 必要がある． 長谷川ら[17] _{はアーチェリー矢の側面境界層流れの線形安定性を行った．その結果，矢周} りの境界層流れはm = 1, 2のモードによる線形不安定な流れであることが示された。特に， m = 1のモードの撹乱は1.0 × 104 ≤ Re ≤ 2.0 × 104の範囲で，成長することが確認された． また，流線形鏃の場合の線形成長率が椎型鏃の場合に比べてわずかに大きくなった． 不安定性の乱流遷移への寄与を調べるため，eN 法[18] _{により解析結果からN 値を評価し} た．風洞実験の場合，外乱強さT u評価し乱流遷移するN 値が推定される．平板境界層流れ において，Mack[19] の式 N = −2.4 ln(Tu) − 8.43が知られている．しかし，円筒面に対する 同様の評価式が見出せなかったため，平板境界層流れの場合の式を使用し N 値を計算した． その結果，N 値の最大値は2.5 程度であり，MSBS風洞実験において迎角0 度の場合に層 流を観測した事実と整合している．さらに，Miyazaki et al.[8]の圧縮空気による飛翔実験で， 初期迎角を制御すれば Re ≤ 1.8 × 104でCd は層流値となっていることから，迎角が0◦で Re ≤ 1.8 × 104の場合，矢周りの境界層流れが層流を維持したことを理論的に裏付けた． 小野口ら[20]_は田中[21]_{の軸周りの回転（Sp= 0.015, 0.025, 0.05）効果を含む，軸対称境界} 層流れの線形安定性解析を平行流近似に基づいて行なっている．その結果，非回転の場合に 主要な役割を果たしていたm= 1, 2 の撹乱の成長率は回転効果により有意に増加した．一方 で周回方向逆向きの撹乱を加えた場合（m = −1, −2），成長率は非回転の場合よりも減少，回 転効果によって撹乱の成長が抑制されることも示している．線形成長するm > 0の撹乱のN 値を計算したところ，SP = 0.025の場合N 値でN = 9.44となった．実際の矢も矢羽の効果 により矢軸周りにSp = 0.015程度で回転しているがMSBS風洞実験結果[8]_{は層流を観測し} ており．線形安定性解析結果とも整合する．したがって，実際の矢の回転数よりもはるかに

大きい回転数に到達しなければ乱流遷移しないと推測される． 一方で，Miyazaki et al.[8]の実験から，気流に対する迎角が乱流遷移の重要なファクターで あることが指摘されている．気流に対する迎角がつくと、矢の先端からm= 1のモードの撹 乱が誘起される．伊藤ら[22] _{はアーチェリー矢の境界層の線形安定性解析におけるm = 1の} モードの固有関数を有限振幅撹乱として田中[21] _{の基本流に重ね合わせ，迎角がついた場合} の流れを再現し，その新たな基本流に対する2次不安定性解析を行った．その結果，挿入し た有限振幅撹乱の影響は先端付近で強く出現し，m = −1のモードが最も不安定となった．ま た，eN 法により解析結果からN 値を評価したところ，有限振幅撹乱の振幅強度が一様流速 の5%に達すると境界層の乱流遷移の可能性が高まることが示された． 椎型鏃および流線形鏃を装着したアーチェリー矢に迎角がついた場合の主流に対し2次不 安定性解析を行う．得られた解析結果を伊藤ら[22]_{の結果や宮嵜ら}[7, 8]_{の結果と比較する．} 本論文は次のように構成される．2章では使用する無次元数の定義を行う．3章では2次 不安定性解析に使用した非軸対称流および軸対称流の数値計算方法および数値計算結果を述 べる．4章では 2次不安定性解析の理論および解析手法を述べる．5章では2次不安定性解 析結果を述べる．6章では得られた2次不安定性解析結果を我々のグループが行っている実 験結果と比較・議論する．7章では本論文の結果を要約し，今後の展望を述べる．

本研究では以下の無次元量を使用して解析,考察を行う.

2.1

使用する無次元量・用語の定義

レイノルズ数: Re 物体まわりの流速をレイノルズ数を用いて,無次元化する. レイノルズ数を以下の式で定義 する. Red = V d ν (2.1) V :矢の速度または主流速度[m/s], d:矢の直径[m],ν:動粘性係数[m2/s] また,半径を代表長さとしたレイノルズ数Rer とRe数は次式の関係がある. Red = 2Rer (2.2) スピンパラメーター: SP 矢の回転数を以下の式で無次元化する. SP = πdf V (2.3) f :矢の毎秒回転数[rps], V :矢の速度[m/s]

第

₃

章

基本流の数値計算

本章では，非軸対称基本流の数値計算方法および軸対称基本流の数値計算方法を述べる． 装着する鏃は競技で使用されている椎型鏃と当研究室で作成した流線形鏃の2種類を用いる． 図3.1に鏃の概形を示す． (a)椎型鏃. (b)流線形鏃. 図3.1: 鏃形状. 鏃形状の数値計算への反映はCAD データを用いる．使用したCAD データを図3.2に示 す．鏃のCADデータは，適当な数の点(100点程度)の座標を取り,スプライン補間を施して 作成された．

図3.2: 鏃（実寸mm）のCAD図，上: 流線形鏃，下: 椎型鏃．長谷川[30]_{より再生．} 数値計算には円筒座標系を用いる．シャフト方向をz，シャフトに垂直かつ離れる方向に r，周回方向をθ とする．z，rは半径で規格化する．原点は鏃とシャフトの境目におけるシャ フトの中心とした．また，矢の進行方向に対して逆向きにz軸の正をとる．

3.1

軸対称基本流の数値計算

流線形鏃と椎型鏃を装着した矢羽なしのシャフトをすぎる，軸対称基本流の数値計算方法 [17]_{について述べる．}

3.1.1

数値計算方法

気流に対して平行な矢の周りの流れは軸対称流とみなせるため，主流はzとr の関数にな る．Navier-Stokes方程式をデカルト座標系から円筒座標系(z, r, θ)に変換し,それぞれの速度 成分を_(Us z, Urs, U_θs)とすると, Navier-Stokes方程式と連続の式は次式のようになる. z軸方向

の原点はシャフトと鏃の継ぎ目とし，r 軸方向の原点はシャフト中心とした． ∂Us z ∂t + Uzs ∂Us z ∂z + Urs ∂Us z ∂r = − ∂p ∂z + 1 Red ( ∂2_Us z ∂z2 + 1 r ∂Us z ∂r + ∂2_Us z ∂r2 ) (3.1) ∂Us r ∂t + Uzs ∂Us r ∂z + Urs ∂Us r ∂r − U_r2s r = −∂p_∂r + 1 Red (_∂2 U_rs ∂z2 + 1 r ∂Us r ∂r − U_rs r2 + ∂2_Us r ∂r2 ) (3.2) ∂Us θ ∂t + Uzs ∂Us θ ∂z + Urs ∂Us θ ∂r + U_rsU_θs r = 1 Red ( ∂2_Us θ ∂z2 + 1 r ∂Us θ ∂r − U_θs r2 + ∂2_Us θ ∂z2 ) (3.3) ∂Us z ∂z + ∂Us r ∂r + U_rs r = 0 (3.4) 流れの基礎式中の未知数はUz, Ur, pの3つであるが,流れ関数・渦度法は流れ関数ψ と方 位角方向の渦度_φ を導入することでNavier-Stokes方程式の未知数を2つに減らすことがで きる. φ = ∂Urs ∂z − ∂Us z ∂r (3.5) U_zs = 1 r ∂ψ ∂z, Urs = − 1 r ∂ψ ∂r (3.6) 流れ関数・渦度法は2次元流れの計算のみに限定されるが,計算処理が短く, 連続の式を厳 密に満たせることが利点である. 式(3.1)と式 (3.2)から圧力項を消去し, 式(3.5) の関係式を使って置き換えを行うと,渦度 輸送方程式(3.7)を得る. ∂φ ∂t + Uzs ∂φ ∂z + Urs ∂φ ∂r = 1 Red ( ∂2_φ ∂z2 + 1 r ∂φ ∂r + ∂2_φ ∂r2 ) (3.7) 式(3.6)を式(3.4)，(3.7)に代入すると次式を得られる． ∂2_ψ ∂r2 − 1 r ∂ψ ∂r + ∂2_ψ ∂z2 = −rφ (3.8) ∂φ ∂t + 1 r ∂ψ ∂r ∂φ ∂z − 1 r ∂ψ ∂z ∂φ ∂r = 1 Red ( ∂2_φ ∂z2 + 1 r ∂φ ∂r + ∂2_φ ∂r2 ) (3.9)

1201× 301とした．シャフト部分の長さは半径を基準とした無次元長さで240になるよう設 定した．図3.3に作成された計算格子を示す．

3.1.3

数値計算結果

数値計算された主流を示す． (a)椎型鏃 (b)流線形鏃 図3.4: 主流の速度分布. Red=2.0 × 104,SP= 0.015, Uz. 図3.4aがRed = 2.0 × 104における椎型鏃を装着した場合，図3.4bが流線形鏃を装着した 時の主流である．椎型鏃では, z=0において,逆流が現れている. これは流れが剥離しているた めである. z > 100鏃の種類による境界層厚さの違いはみられない．

(a)椎型鏃. (b)流線形鏃. 図 3.5: 鏃ごとの主流速度分布 Uz. Red = 2.0 × 104, 1.0 × 104, SP = 0.025. 左から z = 0, 50, 100, 150, 200． 図3.5は主流の速度プロファイルの Red 数毎の比較である. Red 数が大きくなると境界層 が厚くなる. また, Red=2.0×104 における椎型鏃の先端付近の逆流速度は Red=1.0×104の3 倍程度である．

(a)椎型鏃. (b)流線形鏃. 図3.6: 周回方向速度分布U_θ. Red = 2.0 × 104, 1.0 × 104, SP = 0.025. 左からz= 0, z = 50, z = 100, z = 150, z = 200. 図3.6はθ 方向速度プロファイルのRed数毎の比較である. 椎型鏃は鏃の影響が先端で現 れており，U_θ のr 方向速度勾配が流線形鏃と比較して小さい．

(a)椎型鏃. (b)流線形鏃 図3.7: 周回方向速度分布U_θ. Red = 2.0 × 104. SP = 0.015, 0.025, 0.05,左から z=0, z=100, z=200. 図3.7はRed=2.0 × 104におけるθ方向速度プロファイルのSP毎の比較である. 粘性の影 響により，下流になるにつれてr方向速度勾配が小さくなっている.SPの上昇に伴い，r 方向 速度勾配は小さくなる．

3.2

非軸対称基本流の数値計算

非軸対称流の数値計算は岩津ら[25] _{によって行われた．}

3.2.1

数値計算法

非軸対称流の数値計算を行う．非圧縮のNavier-Stokes方程式と連続の式を基礎方程式とす る．計算は一般座標を用いた空間および時間はともに2次精度差分を使用した．時間積分に はJameson[26]の3段低用量ルンゲクッタ法を使用し，Chorinタイプの射影法を適用した．

3.2.2

計算格子の生成

計算格子は代数型格子生成法で格子を生成した. 最小格子幅はz方向は1.0 × 102，r方向は 1.0 × 103とした．格子数は z, r, θ 方向にそれぞれ2073× 171 × 83とした．シャフト部分は 円柱とみなし，長さは半径を基準とした無次元長さで258になるよう設定した．

3.2.3

数値計算結果

図3.8にRed = 15000で，それぞれ椎型鏃および流線形鏃を装着したケースの z方向速度 分布Uz を各 z において2次元カラープロットしたものの全体図を示す．流れは紙面裏側か ら表側にかけて迎角0.5 度がついている．流線形鏃においては定常となったが，椎型鏃では 非定常となった．椎型鏃を装着した場合，z ≥ 10から流れが乱れ，境界層が乱流になった可 能性が考えられる．流線形鏃を装着した場合，z ∼ 50より後端では π₂ < θ < 3₂π にかけて速 度が遅くなる領域があり，迎角の影響が顕著に表れている．

(a)

(b)

図3.8: 迎角0.5 度がついた矢周りの直接数値計算結果のz 方向速度分布の全体図，(a)椎型 鏃，(b)流線形鏃．

図3.9に数値計算結果の z方向速度Uz を示す．流れは紙面の左側かつ紙面の裏面から表面 の方向に迎角0.5度でシャフトに当たっている．カラーマップのスケールが−1 ≤ U_z ≤ 1で あることに注意されたい．境界層厚さは99%境界層厚さ(δ99)で，δ99|z=2 = 1.24，δ99|z=50 = 1.84である．Uz はシャフトを回り込むことによって減速してしまうため，境界層は図の左側 よりも右側の方が厚くなる． 図3.10に数値計算結果の_θ 方向速度分布を示す．Re = 15000で流線形鏃を装着した場合 である．流れは紙面の左側かつ紙面の裏面から表面の方向に迎角0.5 度でシャフトに当たっ ている．カラーマップのスケールが_{−5 × 10}−3 _{≤ Uθ ≤ 5 × 10}−3であることに注意されたい． また，迎角がついているために境界層の外側でもU_θ , 0であることにも注意されたい．迎角 を_γ とすると無限遠でのU_θ は次式の関係がある． U_θ|r=∞ = sin γ sin θ (3.10) 反時計回りの流れを負，時計回りの流れを正としてプロットしている．流れは線対称で，粘 性の影響はシャフトのごく表面でしか表れない．z = 50においては，周り込んだ流れが剥離 し，渦を生成している．非回転の矢であるため，U_θ の大きさは小さい．z = 2の場合，境界 層内においてU_θ ∼ O(10−5)，z = 50の場合，境界層内においてU_θ ∼ O(10−5)である．しか しながら，将来的には回転する矢についても同様の解析を行うことができるよう，U_θ , 0と する．

(a) z= 2. (b) z= 50. 図3.9: 迎角0.5度がついた矢周りの直接数値計算結果のz方向速度分布，Re = 15000，流線 形鏃． (a) z= 2. (b) z= 50. 図3.10: 迎角0.5度がついた矢周りの直接数値計算結果のθ 方向速度分布，Re = 15000，流 線形鏃．

第

₄

章

２次不安定性解析

本章では２次不安定性の基礎理論および解析手法を述べる．平行流近似を適用し，主流をz 方向に独立して解析を行う．また，固有関数をθ方向にFourier級数展開することにより,異 なるモードの撹乱の相互作用を考慮することが可能になる．

4.1

２次不安定性

飛翔中の矢に迎角がついたとする．この時矢側面の流れは，迎角がついていない状態の流 れにさらに撹乱が加わったものとみなすことができる．これを新たな基本流とみなすと，この 基本流に対する不安定性を考えることができる．これを2次不安定性(Secondary Instability) と呼ぶ． ¯ U = U₀+ Au (4.1) ここでU₀は基本流，uは加わった撹乱，Aは加わった撹乱の振幅，U¯ は新たな基本流である．

4.2

基礎方程式

図3.9∼図3.10より，境界層内においてUr のオーダーはUz のオーダーよりも2桁程度小 さく，U_|θのオーダーと比較しても2桁程度小さいため，Ur = 0とみなす．よって新しい基 本流U0は次式のようにかける． ¯ Uz(r, θ, t) = Uz(r) + uz(r) cos θ (4.2)

­­ ­ « u_θ p ®® ® ¬ ­­ ­ « v_θ(r, θ) q(r, θ) ®® ® ¬ ここで軸方向波数_αは実数パラメータで，_ω は複素未知数である．平行流近似を適用したこ とで撹乱は軸方向の波数ごとに分解できる． 式(4.5)をNavier-Stokes方程式に代入し線形化すると線形撹乱方程式(4.6)-(4.9)を得る. i(α ¯Uz− ω)vz + ¯ U_θ r ∂vz ∂θ + v_θ r ∂ ¯Uz ∂θ + d ¯Uz dr vr + iαq = 2 Red[ ∂2_v z ∂r2 + 1 r ∂vz ∂r + 1 r2 ∂2_v z ∂θ2 − α 2 vz] (4.6) i(α ¯Uz− ω)vr + ¯ U_θ r ∂vr ∂θ − 2 ¯U_θ r vθ + ∂q ∂r = 2 Red[ ∂2_v r ∂r2 + 1 r ∂vr ∂r + 1 r2 ∂2_v r ∂θ2 − (α 2₊ 1 r2)vr − 2 r2 ∂vθ ∂θ ] (4.7) i(α ¯Uz− ω)vθ + ¯ U_θ r ∂vθ ∂θ + v_θ r ∂ ¯Uθ ∂θ + ( d ¯U_θ dr + ¯ U_θ r )vr + 1 r ∂q ∂θ = 2 Red[ ∂2_v θ ∂r2 + 1 r ∂vθ ∂r + 1 r2 ∂2_v θ ∂θ2 − (α 2₊ 1 r2)vθ+ 2 r2 ∂vr ∂θ ] (4.8) iαvz+ ∂v_∂rr + vr r + 1 r ∂vθ ∂θ = 0 (4.9) 境界条件は, vz = vr = vθ = 0 (r = 0, r → ∞) (4.10) である. 式(4.6)-(4.9)は偏微分方程式であるため容易には解けない. そのため，式(4.5) の固 有関数をθ方向にFourier級数展開しすると常微分方程式の標準固有値問題 ωAX = B (4.11) に帰着する．Aは単位行列，X は固有関数となる．

4.3

固有関数の周回方向の離散化

線形撹乱(4.5)をθ方向に−M ∼ M まで複素フーリエ級数展開する. ©­ ­­ ­ « vz(r, θ) vr(r, θ) v_θ(r, θ) q(r, θ) ª® ®® ® ¬ = M ∑ m=−M ©­ ­­ ­ « Am(r) Bm(r) Cm(r) Dm(r) ª® ®® ® ¬ eimθ (4.12) 複素数で展開することにより,回転のパラメータを１つの変数のみで扱えるようになる. m番目のFourier成分を式(4.6)-(4.9)に代入すると常微分方程式

i(αUz− ω) Am+ iαuzAm−1+ iαu∗zAm+1

+ ( U_θ r − Ω ) imAm+ u_θ r i(m − 1)Am−1+ u∗_θ r i(m + 1)Am+1 + i ruzCm−1− i ru ∗ zCm+1 + dUz dr Bm+ duz dr Bm−1+ du∗_z dr Bm+1+ iαDm = 2 Red [ d2Am dr2 + 1 r dAm dr − ( α2₊ m2 r2 ) Am ] (4.13)

i(αUz− ω) Bm+ iαuzBm−1+ iαu∗zBm+1

+ ( U_θ r − Ω ) imBm+ u_θ r i(m − 1)Bm−1+ u∗_θ r i(m + 1)Bm+1 − 2Uθ r Cm− 2u_θ r Cm−1− 2u∗_θ r Cm+1+ dDm dr 2 Red [ d2Bm dr2 + 1 r dBm dr − ( α2₊ m2+ 1 r2 ) Bm− 2im r2 Cm ] (4.14) i(αUz− ω) Cm+ u_θ r i(m − 1)Cm−1+ u∗_θ r i(m + 1)Cm+1 + ( U_θ r − Ω ) imCm+ u_θ r i(m − 1)Cm−1+ u∗_θ r i(m + 1)Cm+1 + i ruθCm−1− i ru ∗ θCm+1 + ( dU_θ dr + U_θ r ) Bm+ ( du_θ dr + u_θ r ) Bm−1+ ( du∗_θ dr + u∗_θ r ) Bm+1+ im r Dm = 2 Red [ d2Cm dr2 + 1 r dCm dr − ( α2₊ m2+ 1 r2 ) Cm+ 2im r2 Bm ] (4.15)

円筒座標系(z, r)を座標系(ζ, σ)で表す. 田中[21]_{の数値計算は直径を代表長さにしている} ため, Srikanth et al.[14, 15] の導入した相似変数に修正を加え,次のように定義する. ζ = 2( z Re )1 2 , σ = 2r− 1 ζ r ∈ [0.5, ∞] (4.17) σ の導入でr の定義区間が(1, ∞)から(0, ∞)となる. 微分も以下のように変換される． ∂ ∂z = ∂σ ∂z ∂ ∂σ + ∂ζ ∂z ∂ ∂ζ (4.18) ∂ ∂r = ∂σ ∂r ∂ ∂σ + ∂ζ ∂r ∂ ∂ζ (4.19) したがって,式(4.17)をz, r で微分することで ∂ ∂z = 2 ζ Re ( ∂ ∂ζ − σ ζ ∂ ∂σ ) (4.20) ∂ ∂r = 2 ζ ∂ ∂σ (4.21) を得る.

4.4.2

半径方向の離散化と補間

σ方向の離散化にはChebyshev選点法を採用する．選点数N = 64とし,点の分配を制御す るパラメータσ = 5ˆ とする. ˆσは小さいほどシャフト表面で点が密に分布し,大きいほど表面 付近で疎になる. x = σ − ˆσ σ + ˆσ σ ∈ [0, ∞] → x ∈ [−1, 1] (4.22) xn = cos ( _nπ N− 1 ) 0 ≤ n ≤ N − 1 (4.23)

4.4.3

微分行列の作成

Chebyshev選点法における微分の演算はChebyshevの微分行列で行うことができる. 本研 究ではLloyd[27]の方法で微分行列を作成した．N+ 1点のChebyshev点に対し微分行列Dは N+ 1次の正方行列である．N 次の微分行列 DN のi行j列成分(DN)i j について (DN)00 = 2N2+ 1 6 , (DN)N N = − 2N2+ 1 6 (DN)j j = −xj 2(1 − x2_j), j = 1, . . ., N − 1 (DN)i j = ci cj (−1)i+j (xi− xj), i, j, i, j = 0, . . ., N wher e, ci = { 2, i = 0 or N 1, otherwise (4.24) である.

4.5

固有値問題への帰着

式(4.17)の変換を導入することで, 周回方向に離散化された線形撹乱方程式(4.13)-(4.16) は以下のように変形される.

i( ¯αUz− ¯ω) Am+ i ¯αuzAm−1+ i ¯αu∗zAm+1

+ i(ξUθ − Ω)mAm+ iξuθ(m − 1)Am−1+ iξu∗_θ(m + 1)Am+1

+ iξuzCm−1− iξu∗zCm+1 + U′ zBm+ u′zBm−1+ u∗′z Bm+1+ i ¯αDm = 1 Re_δ [ A′′_m+ ξ A_m′ − ( ¯α2+ ξ2m2)Am ] (4.25)

i( ¯αUz− ¯ω) Bm+ i ¯αuzBm−1+ i ¯αu∗zBm+1

+ i(ξUθ − Ω)mBm+ iξuθ(m − 1)Bm−1+ iξu∗_θ(m + 1)Bm+1

− 2ξUθ− 2ξuθC_m−1− 2ξu∗_θC_m+1+ D_m′

= 2 Re_δ [ B_m′′ + ξB_m′ −{α¯2+ ξ2(m2+ 1)}Bm− 2iξ2mCm ] (4.26)

m ただし, Re_δ = ζ Red, ξ = ζ 1+ ζσ, ¯ α = ζα, ω = ζω¯ (4.29) とした. プライムはσ に関する微分を表す. 境界条件は, Am = Bm= Cm = 0 (σ = 0, σ → ∞) (4.30) である. 以上より, ¯ω について4(N − 2) × 2(M + 1)次の固有値問題 ¯ ω ¯A ¯X = B ¯X (4.31) に帰着する. さらに,圧力項を消去することで3(N − 2) × 2(M + 1)次の固有値問題となり,計 算時間を短縮することができる. 式(4.28)より Am = Bm′ + ξBm+ iξmCm (4.32) となるため, Amを消去することも可能である. 本研究では圧力項の消去のみ行った. 固有値計算はMATLABのeig関数によるQR法を用いた．

第

₅

章

基本流の設定

本章では2次不安定性解析に用いる基本流の作成方法について述べる．

5.1

有限振幅撹乱の挿入

伊藤ら[22] _{が2次不安定性解析に使用した流れの設定手法について述べる．(4.1)式におい} て長谷川ら[17] _{が使用した軸対称流U}s 0に，同じく長谷川ら [17]_{の結果から得たα ∼ 0の固有} 関数(m = 1)を u として用いる．_{α × δ99 = const} となるように_αを設定する．本研究では α × δ99 = 0.01565とする． 図5.1 にz 方向の有限振幅撹乱 |upz| の空間構造を示す．θ 方向の有限振幅撹乱upθ は， upz より 3 桁程度小さい．(図略) なお，up の振幅強度は |upz|max と定義し，目安として 2.5%, 5.0%の場合について解析を行う． 図5.1: upz の空間構造，Re= 20000, z = 5

uc_z = 1 π 2π 0 Uz(r, θ) cos mθdθ (5.2) us_z = 1 π ∫ 2π 0 Uz(r, θ) sin mθdθ (5.3) uc_θ = 1 π ∫ 2π 0 U_θ(r, θ) cos mθdθ (5.4) us_θ = 1 π ∫ _2π 0 U_θ(r, θ) sin mθdθ (5.5) (5.6) フーリエ変換した結果のz 方向速度を図5.2 に，θ 方向速度を図5.3に示す．なお，us_z と u_θc はO(10−5)であるため図表記を省略した．Uz，すなわち基本流の m = 0成分の流れは境 界層となり，r − 1 ∼ 0.3で一様流となる．また，m = 1 の成分の流れは z = 2 において， r− 1 ∼ 0.12で一様流と比較して5%の強度であるのに対し，z = 30においてはr − 1 ∼ 0.2 で15%程度の強度である．また，m = 2, 3の流れも後端になるにつれ発達するが，3%程度 の強度であり，m= 1モードの流れと比較すると無視できるほどに小さい． us_θ は迎角がついているため，m = 1モードの流れは無限遠で定常にはなるが 0にならな いことに注意されたい．z = 2において，m = 1モードの流れはr ∼ 0.6で流れは定常とな る．m= 2モードの流れもみられるがその大きさはm= 1モードと比べると1桁程度小さい． z= 30において，m= 1モードの流れはr ∼ 0.6で流れは定常となる．m= 2モードの流れも みられるがその大きさはm= 1モードと比べると1桁程度小さい．

(a) z=2．

(b) z=30．

(c) z=60．

(a) z=2．

(b) z=30．

(c) z=60．

5.3

流れの重ね合わせ

Re = 20000において非軸対称流の結果が存在しないため，長谷川ら[17] が使用した軸対称 流Us₀に，非軸対称流をフーリエ変換したm = 1モードの流れを uとして挿入する．

5.4

主流の補間

長谷川ら[17] _{が使用した流れの補間結果および微分結果を図5.4に示す．シャフト表面付} 近で微分の振動がみられる． (a) z=2． (b) z=50． 図5.4: 軸対称流れの補間および微分結果，Re=20000，N=64，流線形鏃． 非軸対称流のm= 0モードの流れの補間結果および微分結果を図5.5に示す.

(a) z=2． (b) z=50． 図5.5: m = 0モードの流れの補間および微分結果，Re=15000，N=64，流線形鏃． 青線は補間された主流速度，赤点は補間された点を，紫の点は補間点における微分値を表 す．端点においても微分値は振動しない．

5.5

基本流の比較

長谷川ら[17] _{が使用した主流と，非軸対称流のm= 0モードの流れを図5.6に示す．m = 0} モードの流れは，99%境界層厚さで比較するとz = 2では非軸対称流の方が軸対称流よりも 2.7%厚くなり，z = 200では47%厚くなる．

本章では2次不安定性解析結果について述べる．半径方向補間点数N と周回方向展開項数 M を決定し，2次不安定性解析を行う． 2次不安定性解析は次の3パターンの流れについて行った；迎角のついた基本流をフーリエ 変換し，m= 0モードの流れにm= 1モードの流れを加えた流れ，同様のm= 0モードの流 れにm = 1モードの流れの振幅を変えた流れ，長谷川ら[17] が使用した軸対称流れにm = 1 モードの流れを加えた流れである．

6.1

収束の確認

半径方向補間点数N と周回方向展開項数 M を決定するため，固有値の収束確認を行う．

6.1.1

N

の決定

表6.1 に， N と成長率 ωi の関係を示す．N が2倍になると計算時間は10倍になる．ま た，N ≥ 64で固有値は4桁一致した．計算時間も考慮にいれ，解析ではN = 64とした．

表6.1: Nと成長率ωi の関係. N ωi t[s] 32 -0.22225 2.2 64 -0.22238 23 128 -0.22239 203 256 -0.22239 1748

6.1.2

M

の決定

図6.1に，固有関数を示す．M ≥ 7で固有関数の振幅はO(10−3)以下であるため，解析は M = 7で行う．

(e) m= −5 (f) m= −4 (g) m= −3 (h) m= −2

(i) m= −1 (j) m= 0 (k) m= 1 (l) m= 2

(m) m= 3 (n) m= 4 (o) m= 5 (p) m= 6

(q) m= 7 (r) m= 8 (s) m= 9

6.2

非軸対称流の

m

= 0

モードの流れに同

m

= 1

モードを挿入

した流れに対する 2 次不安定性

6.2.1

固有値解析結果と

α

依存性

2次不安定性解析の結果を図6.2に示す．固有値は連続モードと離散モードが存在し，連 続モードは流れ場が無限領域にわたっているときに現れ，安定である．連続モードはσˆ の値 に大きく依存するが，離散固有値は依存しない．上部の緑線で囲われた箇所の固有値が離散 モードの第一固有値である．また，固有値をσ = 5, 10ˆ で比較したところ5桁の一致を確認 した． 図 6.2: 固有値分布，Re = 15000，z = 20， α = 1.7，N = 64，迎角0.5度． 図6.3: α − ωi 依存性，Re = 15000，N = 64， 流線形鏃，迎角0.5度． αのω_i 依存性を図6.3に示す．z< 100においてはα − ω_i のグラフは上に凸となり極値を とる．先端であるz = 2においてはα = 3.4でωi が最大となる．下流にいくほどωi を最大 にするαの値は小さくなり，z = 90ではα = 0.5である．z ≥ 98では離散固有値は見つから なかった．撹乱の波長 _λは_{λ = 2π/α}で表される．_{α = 0.5}では _{λ ∼ 12}となる．そのため， 不安定になる軸方向波長が長くなる後端においては平行流近似を用いた不安定性解析は適当 ではない可能性がある．

図6.4: z− ωi 依存性，N = 64，迎角0.5度

6.2.3

固有関数

図6.5に z = 2，_{α = 3.4} の場合における第1固有値に対応する固有関数を示す．横軸は 99%境界層厚さで規格化した半径方向距離，縦軸は固有関数の振幅を表す．m = ±3の固有 関数が最も発達している．m= ±4, ±2は振幅が境界層の半分のところにおいて最も大きくな るが，境界層内で0に収束する．

(a) m= −7 (b) m= −6 (c) m= −5 (d) m= −4 (e) m= −3 (f) m= −2 (g) m= −1 (h) m= 0 (i) m= 1 (j) m= 2 (k) m= 3 (l) m= 4 (m) m= 5 (n) m= 6 (o) m= 7 図6.5: 固有関数，Re = 15000, z = 2，_{α = 3.4}，N = 64，_{ωi = −0.039765453}，流線形鏃，迎 角0.5度．

の場合におけるα << 1の撹乱を挿入した． 伊藤らの結果において，5%の場合は先端で撹乱の成長率はO(10−1)となり，不安定性が強 まっている．2.5%の場合は z= 10で撹乱の成長率は最大となり，以降は単調減少する． 本研究の結果は，振幅が1倍の場合，z < 100でωi < 0 となり，安定である．振幅が1.2 倍の場合，z > 36でω_i > 0となり，不安定性が表れた．z= 64でω_i は最大となる．振幅が 1.6倍の場合，不安定性はz > 20で表れる．以降はz = 50までωi は大きくなり，z = 52で さらに増大し，伊藤ら[22]_{の結果よりも大きな値となる．振幅を増大させていくとω} i の値は 大きくなり，流れの不安定性が強まる． 本研究と伊藤ら[22]_{の結果はそのふるまいが大きく異なる．図5.6が示すように，軸対称流} と非軸対称流のm= 0モードの流れは大きくことなっている．この違いが2次不安定性解析 結果に反映されたものと考えられる．

図6.6: z− ωi 依存性．青，赤，黄色の実線は，m= 0モードの流れにm= 1 モードの流れを それぞれ1倍，1.2倍，1.6倍して挿入した結果を表す(Re= 15000，流線形鏃)．赤線に黒丸， 緑線に黒丸は伊藤ら[22]_{らの結果を表す(Re= 20000，椎型鏃)．} 図6.7にm= 1モードの流れの振幅を1.2倍，1.6倍した時のω_i−α依存性を示す．A= 1.2 の場合，_ωi は最大の最大値はz= 2においてα = 3.4，z= 20においてはα = 2.3，z = 50に おいてα = 1.5，z = 70においてα = 1.3である．A = 1.2, 1.6において，ω_i − αグラフは z= 50, 70においてピークが2つ存在する．また，下流になるにつれωi を最大にするαの値 は小さくなっている． 図6.8 にA = 1.6，z = 2におけるωi の最大値に対する固有関数を示す．縦軸は固有関数 の振幅，横軸は_{(r − 1)}を99%境界層厚さで規格化した半径方向距離である．m= ±1が最も 振幅が大きく，続いてm = ±2が大きい．不安定性もm = 1のモードから生じたものと考え られる． 図6.9に A= 1.6，z= 50におけるω_i の最大値に対する固有関数を示す．m= 0，m= ±1， m = ±2，m = ±3，m= ±4の順に固有関数が発達している．そのため，固有値もこれらのい ずれかのモードから生じていると思われるが，そのモードまでは特定できなかった．

(a) A= 1.2 (b) A= 1.6 図6.7: ωi − α依存性，Re= 15000，流線形鏃，迎角0.5度

(a) m= −7 (b) m= −6 (c) m= −5 (d) m= −4 (e) m= −3 (f) m= −2 (g) m= −1 (h) m= 0 (i) m= 1 (j) m= 2 (k) m= 3 (l) m= 4 (m) m= 5 (n) m= 6 (o) m= 7 図6.8: 固有関数，Re= 15000, z = 2，_{α = 3.8}，N = 64，_{ωi = −0.035940}，A= 1.6，流線形 鏃，迎角0.5度．

(e) m= −3 (f) m= −2 (g) m= −1 (h) m= 0

(i) m= 1 (j) m= 2 (k) m= 3 (l) m= 4

(m) m= 5 (n) m= 6 (o) m= 7

図6.9: 固有関数，Re = 15000, z = 50，_{α = 1.7}，N = 64，_{ωi = 0.029636}，A= 1.6，流線形 鏃，迎角0.5度．

(a) m= −7 (b) m= −6 (c) m= −5 (d) m= −4 (e) m= −3 (f) m= −2 (g) m= −1 (h) m= 0 (i) m= 1 (j) m= 2 (k) m= 3 (l) m= 4 (m) m= 5 (n) m= 6 (o) m= 7 図6.10: 固有関数，Re = 15000, z = 68，_{α = 2.4}，N = 64，_{ωi = 0.048560}，A= 1.6，流線形 鏃，迎角0.5度．

挿入した流れの _{ωi − z}依存性を示す．z = 2ですでにωi > 0 であり，不安定性が表れる． z= 20で急激に撹乱の成長率は増大し，z = 78で最大(ωi = 0.258)となる． 撹乱の成長率は伊藤ら[22]_{らの結果と比較すると最大値で 3倍程度高くなり，不安定性が} 非常に強まっている． 図6.11: ωi − α依存性，流線形鏃Re = 20000の軸対称流に流線形鏃Re= 15000の非軸対称 流のm= 1モードの流れを挿入した．

第

₇

章

実験結果との比較と議論

2次不安定性解析結果を実験結果と比較する．図7.1はMSBS風洞実験結果である．矢羽 は非装着のA/C/E矢に椎型鏃を装着した場合のCD のRe依存性を示す．

0

0.5

1

1.5

2

2.5

Re

₁₀

4

0

1

2

3

4

5

C

D

Seban-Bond-Kelly

= 0

= 0.5

= 1.0

= 2.0

図7.1: CD 値のRe数依存性, A/C/E矢，椎型鏃．[30]より引用． 図中の₊は迎角を0度に設定した場合である．黒線で示す理論値 [9, 10] _{に沿うことから，} 層流境界層のCD値に相当する．▲は迎角を0.5度 つけた場合である．12000≤ Re ≤ 14000

よってCD 値は異なるが，定性的には同じ傾向がみられる．|γ| > 1◦ でCD 値は急激に上昇 し，境界層が乱流遷移したものと考えられる．|γ| < 0.6◦では層流のCD 値である．これは2 次不安定性解析で Re= 10000，迎角0.5度の椎型鏃を装着した場合の流れが安定であったこ とと一致する． 図7.2: CD のγ 依存性，Re = 10000.[30] より引用． Ortiz et al.[29] は A/C/E 矢に 3 軸加速度センサを挿入し，矢の回転数を低くするため， Straight-vaneを装着し圧縮空気による飛翔実験を行った．図7.3 に飛翔中の矢にかかる加速 度の変化を示す．t < 0.4sで矢にかかる加速度はおよそ0.65Gで，t > 0.45sではおよそ0.3G である．また，我々のグループが2020年に行った同様の実験結果を図7.4に示す．流線形鏃 を装着したA/C/E矢を使用した．初速は60.0m/sであり，Re = 20000に相当する．流線形鏃 の場合，矢は発射後から層流の境界層で飛翔した．このことは Re = 15000でも 2次不安定 性解析で流れが安定だったことと一致する． 飛翔中の矢の迎角と境界層の状態を時間的に追跡できれば乱流遷移メカニズムの解明の糸

口となる．

図7.3: 矢にかかる加速度の時間変化，初速56.4[m/s]，A/C/E矢，椎型鏃．Ortiz et al.[29]よ

り引用．

図7.4: 矢にかかる加速度の時間変化，初速60.0[m/s]，A/C/E矢，流線形鏃．大場[31]_より引 用．

迎角がついたアーチェリー矢の側面境界層流れに対する2次不安定性を調べた． 迎角0.5度の場合において，Re = 5000, 10000は椎型鏃および流線形鏃ともに定常流であっ た．Re = 15000では流線形鏃は定常流となったが，椎型鏃は非定常流となった．Re = 20000 では椎型鏃および流線形鏃ともに非定常流となった． 非軸対称流をフーリエ変換して得られたm = 0モードの流れにm = 1モードの流れを挿 入した場合，Re = 15000ではシャフト全域において流れは安定であった．安定化の要因は m = 0モードの流れの境界層厚さが軸対称の数値計算結果と比較して厚くなったことだと思 われる．同Re数でm= 1モードの振幅を1.2倍にした場合ではz > 38で不安定性が表れた． 同様に振幅を1.6倍にした場合ではz > 18で流れは不安定となり，撹乱の成長率はO(10−2) であった． Re = 20000の軸対称流れに Re = 15000のm = 1モードの流れを挿入した流れの場合， z > 2で流れは不安定になり，z = 20で撹乱の成長率は急激に上昇する．z = 70で撹乱の成 長率は最大となり，_{ωi ≈ 0.25} であり，撹乱の成長率はO(10−2)となり，強い不安定性が表 れた． 2次不安定性解析結果をMiyazakiet al.[7, 8] _{の実験結果と比較した．MSBS風洞実験結果で} は迎角が0.5度程度では層流境界層のCD値であり，飛翔実験では流線形鏃を装着した場合， Re > 15000でも層流境界層のCD値を観測した．したがって，迎角が0.5度のような微小な 迎角の場合においてはアーチェリー矢周りの流れは層流であることが考えられる．この結果 は2次不安定性解析結果で流れが安定であったことと一致する． したがって，アーチェリー矢の乱流遷移の原因は矢に迎角がつくことだと考えられる．ま

た，導入したm = 1モードの振幅を大きくしたとき，撹乱の成長率が急激に上昇した．この

ことから，遷移のメカニズムは矢に迎角がつくことで境界層に大きな外乱が導入されること により急激に境界層が乱流遷移するバイパス遷移である可能性を検討することが必要である．

9.1

線形安定性理論と固有値問題

本章では線形不安定性解析と解析結果について述べる．

9.1.1

線形撹乱方程式

無次元形のNavier-Stokes方程式と連続の式 ∂U ∂t + (U · ∇)U = −∇P + 1 Re∇ 2_{U (9.1)} ∇ · U = 0 (9.2) から得た定常流( ¯U, ¯P)に微小な撹乱(u, p)が加わったとする． U= ¯U + u, P = ¯P + p (9.3) 式（9.1）から、定常流が満たす方程式、 ( ¯U · ∇) ¯U = −∇ ¯P + 1 Re∇ 2_{U¯ (9.4)} を引くと撹乱方程式 ∂u

∂t + ( ¯U · ∇)u + (u · ∇) ¯U + (u · ∇)u = −∇p + 1 Re∇ 2 u (9.5) を得る．連続の式（9.2）についても同様に定常流の場合の連続の式を引くと， ∇ · u = 0 (9.6)

を得る．境界条件は，撹乱は表面および無限遠で消失するとして u= 0 (9.7) とする． さらに，撹乱uの振幅が微小である場合は式（9.5）の2次の微小非線形項_{(u · ∇)u}を無視 する．この場合、式（9.5）は線形撹乱方程式， ∂u

∂t + ( ¯U · ∇)u + (u · ∇) ¯U = −∇p + 1 Re∇ 2_{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∇ 2_{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¯ は_θ方向に対して不 変であるため、撹乱は次の式のようにθ方向にスペクトル分解される。 ©­ ­­ ­ « uz ur u_θ p ª® ®® ® ¬ = ©­­­ ­ « vz(r, z) vr(r, z) v_θ(r, z) q(r, z) ª® ®® ® ¬ exp[i(mθ − ωt)] (9.15) 整数mは周回方向波数である． 平行流近似に基づく場合，基本流はz方向に対しても不変である．したがって，z方向に対 してもスペクトル分解可能であり、 ©­ ­­ ­ « uz ur u_θ p ª® ®® ® ¬ = ©­­­ ­ « vz(r) vr(r) v_θ(r) q(r) ª® ®® ® ¬ exp[i(αz + mθ − ωt)] (9.16) のように記述できる。固有関数はr方向のみの関数となり，α, mは任意のパラメータである， 平行流近似に基づくためUr = 0とする．式（9.16）を円筒座標系のNavier-Stokes方程式に 加えて線形化すれば，円筒座標系の線形撹乱方程式を得る． i(αUz − ω) vz+ Uz′vr + iαq = 2 Re [ v_z′′+ 1 rv ′ z − ( α2₊ m2 r2 ) vz ] (9.17) i(αUz − ω) vr + q′= 2 Re [ v_r′′+ 1 rv ′ r − ( α2₊ m2+ 1 r2 ) vr − 2im r2 vθ ] (9.18) i(αUz − ω) vθ + im r q= 2 Re [ v_θ′′+ 1 rv ′ θ− ( α2₊ m2+ 1 r2 ) v_θ+ 2im r2 vr ] (9.19) iαvz + vr′ + 1 rvr + im r vθ = 0 (9.20)

ここで、′は半径方向の微分を意味する．境界条件は， vz = vr = vθ = 0, (r = 0, r → ∞) (9.21) である．

9.1.4

相似変数の導入と Chebyshev 多項式による離散化

円筒座標系_{(r, θ, z)}に対して以下の相似変数を導入する． ζ = 2( z Re )1 2 , σ = 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点 と書く）xn は次の様に定義される． xn = cos ( _nπ N− 1 ) 0≤ n ≤ N − 1 xn ∈ [−1, 1] (9.25) σ ∈ [0, ∞]であるため式（9.26）により、Chebyshev点の区間[-1,1]と半無限区間_{[0, ∞)}を対 応させる． σn = ˆ σ(1 + xn) 1− xn σn ∈ [0, ∞] (9.26) xn = −1 → σn = 0，すなわちシャフト表面に対応し，無限遠では xn = 1 → σn = ∞に対応 する．ここでσˆ は点の分配の疎密を決めるパラメータである． 任意の点_σn における値 f(σn)は，σ および f(σ)の値を使用してspline補間により決定 する．

hn(x) = k=0,k,n xn − xk = δnk 0 ≤ n ≤ N − 1 (9.28) と書ける．式（9.27）をxで微分し， dV(x) dx = N−1 ∑ k=0 DnkVk, Dnk = dhk(x) dx (9.29) とする． 2階微分D(2)，4階微分D(4)はDをそれぞれ2乗，4乗することで導出される． D(2)_nk = DniDik (9.30) D(4)_nk = D(2)_niD(2)_ik (9.31) 物体表面の圧力は定義されないため圧力項のChebyshev補間多項式Pの次数はN− 1次と なる． P = N−2_∑ k=0 Pkˆhk(x) (9.32) ˆhn(x) = xn+ 1 x+ 1 hn(x) 0 ≤ n ≤ N − 2 (9.33) 微分行列を_Dˆ として， ˆ Dnk = xk + 1 xn + 1 Dnk − 1 xn+ 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 = N∑−2 k=1 Dnkvk (9.35) d2v dx2   x=xn = N∑−2 k=1 D2_nkvk. (9.36) 圧力場の撹乱qに関しても無限遠に対応する点を省略するためN − 3次の多項式展開， q= N∑−2 k=1 qk˜hk(x) (9.37) ˜h_{k(x) =} xn + 1 x+ 1 · xn − 1 x− 1 · hn(x) 1 ≤ n ≤ N − 2 (9.38) で表される．よってその微分は dq dx   x=xn = N−2 ∑ k=1 ˜ Dnkqk (9.39) ˜ Dnk = 1 x_n2− 1{(x 2 k− 1)Dnk − 2xnδnk} 1 ≤ n, k ≤ N − 2 1 ≤ n, k ≤ N − 2 (9.40) で表される． したがって軸対称問題における線形撹乱方程式（9.17）-（9.20）は次式で書ける． i( ¯αUz − ¯ω)vz + Uz′vr + i ¯αq = 1 Re_δ{v ′′ z + ξvz′− ( ¯α2+ m2ξ2)vz} (9.41) i( ¯αUz − ¯ω)vr + q′= 1 Re_δ[v ′′ r + ξvr′ − { ¯α2+ (m2+ 1)ξ2}vr − 2imξ2vθ] (9.42) i( ¯αUz − ¯ω)vθ + imξq = 1 Re_δ{v ′′ θ + ξvθ′ − ( ¯α2+ (m2+ 1)ξ2)vθ+ 2imξ2vr} (9.43) i ¯αvz + vr′ + ξvr + imξvθ = 0 (9.44) ここで Re_δ = ζ Re/2, ξ = ζ 1+ ζσ, (9.45) ¯ α = ζα, ω = ζω.¯ (9.46)

（9.49）で書ける． ¯ ωA ¯X = B ¯X (9.49) 単位行列をI として、 A=©­­­ ­ « I I I 0 ª® ®® ® ¬ , B =©­­­ ­ « b00 b01 0 b03 0 b11 b12 b13 0 b21 b22 b23 b30 b31 b32 0 ª® ®® ® ¬ . (9.50) 各小行列のn行k 列成分および_σに関する微分行列を， (b00)nk = ¯αUznδnk + i Re_δ ( D(2)_nk + ξnD(1)_nk − ( ¯α2+ m2ξn2)δnk ) , (b01)nk = −iU ′_n z δnk, (b03)nk = ¯αδnk, (b11)nk = (b22)nk = ¯αUznδnk + i Re_δ [ D(2)_nk + ξnD_nk(1)− ( ¯α2+ (m2+ 1)ξn2)δnk ] , (b12)nk = 2ξn m Re_δξnδnk, (b13)nk = −i ˜D (1) nk, (b21)nk = − 2m Re_δξ 2 nδnk, (b23)nk = mξnδnk, (b30)nk = ¯αδnk, (b31)nk = −i(D(1)_nk + ξnδnk), (b32)nk = mξδnk D(1)_nk = 1 2 ˆσ(1 − xn) 2_D nk, D(2)_nk = 1 4 ˆσ2(1 − xn) 3_{((1 − x} n)D2_nk − 2Dnk), ˜ D(1)_nk = 1 2 ˆσ(1 − xn) 2_D˜ nk, ˜ Dnk = 1 x_n2− 1((x 2 k − 1)Dnk − 2xnδnk), と定義する．境界条件より1 ≤ n, k ≤ N − 2の範囲で定義される．各小行列は N− 2次の正 方行列であるから A, Bは4N− 8次の正方行列となる．

9.1.7

圧力項の消去

圧力項を消去する．これにより，圧力に対する境界条件を削除でき，また行列のサイズを 縮小し，固有値問題の計算時間を節約できる．式（9.50）を， ¯ ωv_z = b00v_z+ b01v_r + b03q (9.51) ¯ ωv_r = b11v_r + b12v_θ+ b13q (9.52) ¯ ωv_θ = b21v_r+ b22v_θ + b23q (9.53) 0= b30v_z+ b31v_r+ b32v_θ (9.54) と書き直す．v_z = (v_z1, . . ., v_zN−2)であり，他の成分についても同様である．式（9.54）より， 係数行列を左からかけて， ¯ ωb30v_z = b30b00v_z+ b30b01v_r + b30b03q, (9.55) ¯ ωb31v_r = b31b11v_r + b31b12v_θ+ b31b13q, (9.56) ¯ ωb32v_θ = b32b21v_r + b32b22v_θ + b32b23q. (9.57) したがって， q= Pzvz + Prvr + Pθvθ (9.58) と書ける．ただし． Pz = −(b30b03+ b31b13+ b32b23)−1b30b00, (9.59) Pr = −(b30b03+ b31b13+ b32b23)−1(b30b01+ b31b11+ b32b21), (9.60) P_θ = −(b30b03+ b31b13+ b32b23)−1(b31b12+ b32b22) (9.61) である．したがって，固有値問題は式（9.62），（9.63）のように3N− 6次となる． ¯ ω ¯Y = ¯B ¯Y, _Y¯ _{= (v} z, vr, vθ) (9.62) ¯ B=©­_­ « b00 b01 0 0 b11 b12 0 b21 b22 ª® ® ¬ +©­_­ « b03Pz b03Pr b03Pθ b13Pz b13Pr b13Pθ b23Pz b23Pr b23Pθ ª® ® ¬ (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，流線形鏃．

9.2

Bi-global 安定性解析

本章ではBhoraniya and Vinod[32] _{の方法を基にbi-global線形不安定性解析について述べ} る．軸対称基本流をz方向とr 方向に展開し，流れの’全体的な’安定性を調べる．なお，簡 単のために軸対称基本流は非回転とする．

9.2.1

基礎方程式

bi-global線形不安定性解析には非圧縮のNavier-Stokes 方程式を用いる．流れに微小な撹 乱が導入されるとすれば，流れの速度成分および圧力成分は ¯ U = U + u, ¯V = V + v, ¯W = 0 + w, ¯P = P + p,  (9.64) と書ける．さらに微小な撹乱は次の式で表される．