SENAC Vol. 52, No. （2 2019. 4）― 9 ―［共同研究成果］

直接数値解析による環状流路ポアズイユ流の亜臨界遷移の研究―

壁乱流遷移における局在乱流パターン形成の解明に向けて―塚原 隆裕†・ 石田 貴大‡†‡

東京理科大学 理工学部 機械工学科
宇宙航空研究開発機構 航空技術部門

円管内や平行平板間などの壁乱流に見る亜臨界遷移，そのメカニズムと普遍則の解明には理論的アプローチが困難で，従来は風洞や水路実験に頼ってきた．近年は，大型ベクトル並列計算機（東北大学サイバーサイエンスセンター保有の SX-ACE 等）により，直接数値シミュレーション（DNS）を用いたアプローチが可能となってきた．本研究では，圧力勾配駆動の二重円管内環状ポアズイユ流を対象にして，乱流維持の最小臨界レイノルズ数付近で観察される層流—乱流共存パターンに着目することで，パターン形成の有無と変化の連続性を明らかにすることを目指した．
特に，円筒比およびレイノルズ数依存性の調査により亜臨界遷移現象の普遍性を検討している．

1. 緒言現在の流体工学関連の“ものづくり”では，模型実験に先立って CFD（数値流体力学）シミュレーションにより設計指針を立て，機器形状や運転条件の最適化をある程度まで図ることが一般的な開発プロセスになっている．計算機環境と CG 技術の進歩によって益々，CFD 解析の担う役割・場面は増えてきたが，依然として予測困難の現象が層流—乱流遷移である．流れの状態は大別して，層流と乱流があり，前者では流れの時空間的変化が規則的または緩やかなため予測も比較的容易である．しかし，航空機の気流や発電プラントの冷却材循環系など，様々な流体流れの多くは乱流状態にある．微細な乱流渦が流体をかき乱すことで，運動量・熱・物質の混合が強められ，結果的に壁面摩擦・熱伝達・物質拡散が層流に比べて格段に促進される．これが機器性能に直結するため，問題となる流れが層流と乱流のいずれであるかを予測し，伝熱量などを見積もることは工学的に非常に重要である．よく教科書では，円管内流れの臨界レイノルズ数（一解釈として“層流⇔乱流が切り替わる無次元流速の閾値”）は 2000 や 2300 [1-4]と明記されているが，これは飽くまでも経験則であり，数学的または物理学的な裏付けに乏しい．実際，撹乱条件や入口形状の影響[5-7]を受ける上に，異なる流動形態（例えば矩形管やクエット流）になれば臨界レイノルズ数が大きく変わり未知であることも多い．最近では，後述する大規模な層流—乱流共存パターンを伴う場合には臨界値が下がることが，本研究を含め最新の実験や DNS によって明らかにされてきている[8-12]．しかし，そのパターン形成を十分に捉えるには巨大な水路や計算領域を要求されるため，大型並列計算機の利用は避けられない．
以上のように，乱流遷移の解明は未だに困難な課題であるが，より厳密な最小臨界レイノルズ数（あるいは乱流維持限界）の特定を目標に，直接数値シミュレーション（DNS）を用いた研究を行っている．本研究対象は二重円筒間を圧力勾配により駆動されて流れる，いわゆる環状ポアズイユ流とした．当該系は，規範的（カノニカル）な流れとして古くから乱流および遷移現象の研究対象として扱われている円管内流と平行平板間流において，その中間的な流路形態に相当するものである．内外円筒の半径比に応じて円管内流または平行平板間流に漸近することから，異なる流路形態に見るそれぞれの遷移過程を統一的に説明できる可能性を有した場である．以下に，円管内流と平行平板間流の先行研究と本研究目的を記す．

― 10 ―

SENAC Vol. 52, No. （2 2019. 4）

図 1 解析対象の環状ポアズイユ流．x 方向に一定の平均圧力勾配を課して，流れを駆動している．円筒座標系（x, r, θ）において，各方向速度成分を（ux, ur, uθ）と定義する．
主流方向（x）には周期境界条件を課している．

壁面せん断流れの乱流遷移は亜臨界過程（線形安定性理論に基づく臨界値よりも低いレイノルズ数で生じる，有限撹乱誘起の乱流遷移過程）でよく特徴づけられ，最小臨界レイノルズ数付近
では乱流と層流状態が空間的に同時に現れる層流—乱流共存パターンが発生し，最近では有向パーコレーションの普遍クラスとの関連も示唆されている[13-16]．パターンは流路形態に依存して，流れ方向に局在化した一次元的構造である乱流パフが円管内流[1, 7, 17, 18]で，層流—乱流界面が主流に対して傾きを有する縞状の二次元的パターンが平行平板間流[8-16]で生じることが知られている．いずれの層流—乱流共存パターンも，流路代表長さ（直径またはチャネル幅）に比べて数十倍～数百倍のスケールを有し，最小臨界レイノルズ数を下回るまで支配的な構造となっている．
各構造のメカニズムは乱流パフ[19]と縞状パターン[20]のそれぞれについて議論が進められている．円管と平行平板間の本質的差異はスパン（円周）方向自由度の有無と曲率の影響にあると考えられる．そして，スパン方向空間領域の制限度合によって乱流パフと縞状パターンの異なる形状間で変化するものと推測される．しかし，既存研究では個々の系に対する調査に限られており，パターンの維持機構の流動形態依存性及びパターン形状の変化過程が明らかではない．
本研究では，乱流パフと縞状パターンを繋げ得る同心二重円筒間（図 1 を参照）の環状ポアズイユ流を用いることで，円筒比依存のパターン形成（図 2）とパターン変化の連続性について調査した．層流―乱流共存パターンは大規模なスケールで生じるため，東北大学サイバーサイエンスセンター保有のベクトル計算機 SA-ACE を利用し，大規模計算領域による DNS を実施した．

2.

解析対象と計算手法

解析対象は，図 1 のような同心二重円筒間の環状ポアズイユ流れで，軸方向に一様な圧力勾配が付加されることで流れは駆動される．重要な幾何学パラメータである円筒比 η = ri/ro（ri：内円筒半径，ro:外円筒半径）を変えることで，周方向に閉じた系としてスパン方向長さが決定される．
特に η → 0 の場合は事実上，円管内流のものに似た流路となり，η → 1 では平行平板間流に漸近する．内外円筒半径差は d とする．
支配方程式は，円筒座標系での非圧縮性流体の連続の式とナビエ・ストークス方程式であり，
  u 0(1)u1
 (u   )u   p   2ut

(2)

ここで，u は速度ベクトル，ρ は密度，ν は動粘度，p は圧力である．空間的離散化に有限差分法を用いた．時間的離散化に二次精度 Crank-Nicolson 法と二次精度 Adams-Bashforth 法を用いた．

直接数値解析による環状流路ポアズイユ流の亜臨界遷移の研究― 11 ―(a) Helical turbulence at η = 0.5(b) Helical puff at η = 0.3(c) Straight puff at η = 0.1
図 2 環状ポアズイユ流における層流—乱流共存パターン．円筒比（η = 0.5，0.3，0.1）に応じて，螺旋乱流，螺旋パフ，（直線型）乱流パフを形成する．左側は三次元可視化図：(a, b) 周方向速度（赤 uθ > +1［x 下流方向に向かって時計回転］，青 uθ < −1［反時計回転］）；(c) 主流方向速度変動（赤 ux′ > +4，青 ux′ < −1）．右側は局在乱流領域の模式図．
円筒比に依存したパターンの変化を捉えるため，円筒比 η は 0.1 ⁓ 0.8 の範囲で，完全に発達した乱流場の解析からはじめ，摩擦レイノルズ数 Reτ = uτd/2ν（摩擦速度 uτ）を段階的に下げることで遷移レイノルズ数域の解析を行った（Reτ = 48 ⁓ 150）．計算領域はパターンを捉えるため，軸方向に十分長い領域（Lx/d = 51.2 ⁓ 180）を設け，周方向は全周（Lθ = 2π）としている．格子数は軸方向に最大 4096，周方向に最大 1024，壁垂直方向には不等間隔格子を用いて 128 点とした．

― 12 ―

SENAC Vol. 52, No. （2 2019. 4）

図 3 円筒比 η = 0.1 の環状ポアズイユ流における（直線型）乱流パフ，及び各円筒比の乱流間欠率のレイノルズ数依存性．三次元可視化図は瞬時の主流方向速度変動分布を示す．グラフの横軸 Rem は，バルク平均速度と d に基づくレイノルズ数．

3.

結果と考察

3.1

層流­乱流パターン変化

各極限と中間の円筒比におけるパターン発生の有無を調査し，図 2 のように円筒比に依存して(a) 螺旋乱流，(b) 螺旋パフ，(c) 乱流パフの異なった特徴をもつパターン形成が明らかになった．
さらに，乱流パフについてはレイノルズ数の減少に伴い，乱流パフの主流方向間隔が長くなることでより間欠的な乱れ場となる傾向が観測された．図 3 には円筒比 0.1 で見られた乱流パフのレイノルズ数依存性を示し，合わせて乱流間欠率とバルクレイノルズ数の関係を円筒比毎にプロットした．乱流間欠率は，計算領域に対する局在乱流領域の占有率を意味し，0 は完全に層流，1 に近いほど全域が乱流状態であることを示す．レイノルズ数が減少するにつれ，乱流間欠率は単調に減少するが，円筒比によって変化傾向に違いがある．特に，図に示したケースの内で唯一，乱流パフを形成する円筒比 0.1 では Rem = 1600 付近で急激な変化が見られた．これは，僅かなレイノルズ数の変化で乱流パフの個数が増減するためであり，有限な計算領域を用いている限りは連続的遷移と言うよりはむしろ階段状のグラフになってしまう．厳密な最小臨界レイノルズ数を特定するには主流方向に長い計算領域の設定が肝要である．今後，MPI 並列計算を用いた超長流路解析による連続的遷移および最小臨界レイノルズ数の解明が期待される．
本研究では，図 2 のようにパターン形成が円筒比に依存して変化する際の特徴と原因について，特に螺旋乱流で見られる大規模二次流れに注目して，さらなる調査を行った．

3.2

縞状パターンの傾斜性

螺旋乱流のように縞状パターンの維持機構は，空間的局在化した乱流領域の縁で生じる大規模二次流れによって説明され，小規模と大規模流れの明確なスケール分離が重要な仮定となる[20]．
小規模流れは，準秩序的な乱流運動であるストリークと縦渦構造に対応し，内層スケール（uτ/ν）

直接数値解析による環状流路ポアズイユ流の亜臨界遷移の研究― 13 ―

図 4 螺旋乱流における瞬時の x-θ 二次元速度分布：（上段）ギャップ中央の半径方向速（中）半径方向平度 ur 分布で赤・青は正・負の速度を示し，緑は層流領域に相当する；
均された x-θ 面内速度ベクトル（Ux, Uθ）；（下）ローパスフィルタで抽出された大規模二次流れを表す速度ベクトル（ℒUx, ℒUθ）．DNS の計算条件は η = 0.5，Reτ = 56．
によって特徴付けられ，パターンの乱流領域域内の微細構造に相当する．これに対して大規模流れは，外層スケール（d）で特徴付けられ，パターンの大域的な水平方向のサイズを表す．従って，両スケールの間には大きな隔たりが存在する．縞状パターンの（主流方向に対する）傾斜性の起源は大規模流れの質量平衡に由来するものと考えている．小規模流れと大規模流れ（Ui = ℒui，ℒ：ローパスフィルタ）が明確に分類でき，大規模流れのみを抽出できたと仮定すると，以下に示す大規模流れの連続の式が得られる．
U x 1 r U r 1  U0
r rr x(3)縞状パターンは半径方向に対してほぼ一様となり，主流 Ux とスパン方向 Uθ の二次元流れとみなせるため，半径方向平均を施すことで壁面すべり無し条件から以下の式のように変換できる．

Ro U

Rix

xdr  

Ro 1  U r  dr

Ri(4)パターン発生時における∫Ux dr は層流領域と乱流領域で値が異なり，その界面では局所的に，主流方向の勾配が生じる．

Ro U

Rix

xdr  0(5)これにより，式(4)でスパン方向速度成分と釣り合う必要が生じ，速度変化が顕著な層流−乱流界面においてスパン方向二次流れが誘起され縞状パターンは維持される．
図 4 は螺旋乱流の瞬時速度場に対してローパスフィルタを適用し，大規模二次流れを抽出した結果を示す．但し，速度ベクトルは式(4)で施しているように半径方向の空間平均された値を用いている．図 4 内の三図を比較すれば，局在乱流の傾斜帯に沿うような二次流れが上流側・下流側で対向して生じているのが分かる．

― 14 ―

SENAC Vol. 52, No. （2 2019. 4）

50
Lzo ()Lzi ()Lzi/d, Lzo/d40
3020
10000.1

0.2

0.3

0.40.5

0.60.70.80.91

図 5 円筒比に依存する周方向空間長さ（内円筒円周 Lzi と外円筒円周 Lzo）．

4.2

環状流におけるパターン形状の変化

環状流中のパターン形成変化の主要因は，円筒比 η に依存した周方向領域サイズ Lzi = 2πRi または Lzo = 2πRo と考えられる．図 5 に，円筒比に対する内円筒と外円筒壁面における周方向長さを示す．円筒比が大きい場合，周方向長さは無限大に漸近し，低円筒比では周方向長さは最大でも 6d程度となる．縞状パターンの維持機構は小規模と大規模流れの明確なスケール分類が重要となる．
しかし，円筒比が小さい場合，大規模流れ自体が維持可能な周方向領域を確保できないことが考えられる．図 2 からも推察して，円筒比 0.3 程度を下回ると螺旋形状が失われていくことから，おおよそ 10d が限界の外円筒円周の狭さと言える．
パターン形状の変化を定量的に捉えるため，大規模流れの局所速度の向き（主流方向に対する角度）α を抽出することで，その確率密度分布を得る．ここで大規模流れを得るために主流方向は，パターンのサイズを基にして 20d 以下の波長を除去した．また，周方向に関しては円筒比が小さい場合，大規模流れ相当の計算領域を保有できないため，ストリーク相当の小規模流れ（2d）のみを除去した．遷移レイノルズ数域における α の確率密度分布の結果を図 6 に示す．円筒比が大きい場合に α ≠ 0°でピークを有し，円筒比低下に従って α = 0°でシングルピークの凸分布となる．前者の傾向は大規模二次流れの傾斜性を，一方で後者は特徴的な二次流れの不在を示唆している．ここで，α = 0°付近の分布を偶関数，例えば二次曲線近似し，PDF( ) a 2  c

(6)で定義される係数 a(η, Reτ)に着目する．当係数が正ならば a ≠ 0 で，負ならば a = 0 でピークを有した確率密度分布であることに相当する．さらに a(η, Reτ) = 0 が得られれば，パターン変化（螺旋パフ⇔直線型乱流パフ）の臨界条件 ηcritical(Reτ)が特定できる．図 7 には，図 6 で示すような α の確率密度分布のピーク位置，いわゆる最頻角度 M(α)と，係数 a(η, Reτ)を示す．Reτ = 80 の比較的に高いレイノルズ数では発達した乱流状態（小規模流れ支配）であるため，円筒比に依存せずに M(α)= 0°となる．
これに対して，低いレイノルズ数では高円筒比で M(α) ≈ 10°の高い値を示すことから，周方向大規模流れが顕著に現れていることが分かる．一方，0.4 以下の低円筒比では M(α)が小さくなり，大規模流れの特徴が無くなっている．同様の結論が，係数 a(η, Reτ)の結果からも伺える．
特に注目すべきは，η = 0.2 ⁓ 0.4 の範囲で連続的な変化を示すことである．これは，円筒比に応じて螺旋乱流から乱流パフへの変化に連続性を有することを示唆している．
この結果から層流—乱流共存パターン形状の変化原因を考察する．円筒比が十分に大きい場合，平板間流れの理論を直接適用することが可能となり，主流方向大規模流れの変化が周方向大規模流れの変化と釣り合うことにより，環状流では螺旋乱流を形成する．しかし，円筒比が小さくな

直接数値解析による環状流路ポアズイユ流の亜臨界遷移の研究― 15 ―

PDF()

0.016
0.012

Re  = 80
Re  = 64
Re  = 56Re  = 520.0080246810

図 6 大規模流れ中の局所風向 α の確率密度分布．下図は η = 0.5 のものに限定．

図 7 大規模流れ中の局所風向 α の確率密度分布（図 6）に関する最頻角度 M(α)，及び式(6)で定義される係数 a の円筒比依存性．
ると，周方向空間が狭くて大規模流れが存在しなくなり，結果として周方向二次流れが消失した．
螺旋乱流から乱流パフへと，中間的な形態である螺旋パフを介して，層流—乱流共存パターンの形状変化が起きたと考えられる．

5.

結言

環状ポアズイユ流の亜臨界遷移過程について，特に層流—乱流共存パターンの円筒比依存性およびレイノルズ数依存性を調査した．大規模な計算領域による DNS 実施によって，螺旋乱流・螺旋パフ・乱流パフの発生を捉え，円筒比 0.3 ⁓ 0.4 で局在乱流領域の傾斜性を徐々に失い，低円筒比

― 16 ―

SENAC Vol. 52, No. （2 2019. 4）にて円管内流れのものと酷似した亜臨界遷移過程を経ることを見出した．これは異なる層流—乱流共存パターン形状の強い関連性を示唆するものである．
当該研究成果は文献[21, 22]にて詳細に報告している．さらに，摩擦係数や各種乱流統計量については Ishida & Tsukahara [23]で，駆動方法の異なる環状クエット流については Kunii et al. [24, 25]を参照されたい．
層流—乱流転移の連続性の議論には，まだ膨大な試行と巨大計算領域が必要とされる．本解析も含めて，そのような DNS の実現には高性能のスーパーコンピュータ利用が不可欠であり，MPI 並列プログラムの高効率化も必須である．引き続き，東北大学サイバーサイエンスセンターの共同研究などを通してプログラム高性能化および高度なコンピュータ環境利用が望まれる．

謝辞
本研究は，東北大学サイバーサイエンスセンターのスーパーコンピュータ SX-ACE を利用することで実現した．計算にあたっては同センター関係各位に有益なご指導とご協力をいただいた．
研究費の一部は科学研究費補助金として，新学術領域（研究領域提案型）「壁乱流亜臨界遷移の間欠乱流パターン形成の大規模 DNS 解析（16H00813）」および「亜臨界乱流遷移におけるグローバル安定性と大規模間欠構造の複雑流への研究展開（16H06066）」から支援を受けたものである．

参考文献[1] J. Wygnanski and F. H. Champagne, On transition in a pipe. Part 1. The origin of puffs and slugs andthe flow in a turbulent slug, J. Fluid Mech., 59, 281–335, 1973.
[2] F. M. White, Fluid Mechanics (4th edition), McGraw-Hill, New York, pp. 325–332, 2001.
[3] 日野幹雄，流体力学，朝倉書店，1992[4] JSME テキストシリーズ，流体力学，日本機械学会，2005．
[5] 小川益郎，河村洋，円管内ガス流の遷移域の圧力損失と熱伝達に及ぼす入口形状の影響，日本機械学会論文集 B 編，第 52 巻 第 477 号，2164–2169，1986．
[6] B. Hof, A. Juel, and T. Mullin, Scaling of the turbulence transition threshold in a pipe. Phys. Rev. Lett.,
91, 244502, 2003.
[7] J. Peixinho and T. Mullin, Finite-amplitude thresholds for transition in pipe flow, J. Fluid Mech., 582,169–178, 2007.
[8] A. Prigent et al., Large-scale finite wavelength modulation within turbulent shear flows. Phys. Rev. Lett.,
89, 014501, 2002.
[9] T. Tsukahara et al., DNS of turbulent channel flow at very low Reynolds numbers. In: Proc. 4th Int.
Symp. Turbulence and Shear Flow Phenomena, pp. 935–940, 2005; arXiv Preprint, 1406.0248.
[10] T. Tsukahara et al., Turbulence stripe in transitional channel flow with/without system rotation, In: Proc.
7th IUTAM Symp. Laminar-Turbulent Transition, Springer Netherlands, 421–426, 2010.
[11] 塚原隆裕，石田貴大，平面ポアズイユ流の亜臨界遷移における下臨界レイノルズ数，日本流体力学会誌「ながれ」，第 34 巻 第 6 号，383–386，2015.
[12] J. J. Tao, B. Eckhardt, and X. M. Xiong, Extended localized structures and the onset of turbulence inchannel flow, Phys. Rev. Fluids, 3, 011902, 2018.
[13] M. Sano and K. Tamai, A universal transition to turbulence in channel flow, Nature Phys., 12, 249–253, 2016.

直接数値解析による環状流路ポアズイユ流の亜臨界遷移の研究― 17 ―[14] G. Lemoult et al., Directed percolation phase transition to sustained turbulence in Couette flow. Nature Phys.
12, 254–258, 2016.
[15] M. Chantry, L. S. Tuckerman, and D. Barkley, Universal continuous transition to turbulence in a planar
shear flow, J. Fluid Mech., 824, R1, 2017.
[16] M. Shimizu and P. Manneville, Bifurcations to turbulence in transitional channel flow. arXiv
Preprint:1808.06479, 2018.
[17] K. Avila et al., The onset of turbulence in pipe flow, Science, 333, 192–196, 2011.
[18] D. Barkley et al., The rise of fully turbulent flow, Nature, 526, 550–553, 2015.
[19] M. Shimizu and S. Kida, A driving mechanism of a turbulent puff in pipe flow, Fluid Dyn. Res., 41,045501, 2009.
[20] Y. Duguet and P. Schlatter, Oblique laminar–turbulent interfaces in plane shear flows. Phsy. Rev.
Lett., 110, 034502, 2013.
[21] T. Ishida, Y. Duguet, and T. Tsukahara, Transitional structures in annular Poiseuille flow depending on
radius ratio, J. Fluid Mech., 794, R2, 2016.
[22] T. Ishida, Y. Duguet, and T. Tsukahara, Turbulent bifurcations in intermittent shear flows: from puffs tooblique stripes, Phys. Rev. Fluids, 2, 073901, 2017.
[23] T. Ishida and T. Tsukahara, Friction factor of annular Poiseuille flow in a transitional regime, Adv. Mech.
Eng., 9, 1, 2017. doi: 10.1177/1687814016683358[24] 國井康平，石田貴大，塚原隆裕，遷移域スライディング・クエット流の局在乱流構造と統計量に及ぼす壁面曲率の影響，日本流体力学会誌「ながれ」，第 35 巻 第 6 号，475–480，2016.
[25] K. Kunii et al., Laminar-turbulent coexistence in annular Couette flow, J. Fluid Mech., under review,
2019; arXiv Preprint, 1904.09160.