― 32 ―SENAC Vol. 55, No. （4 2022. 10）[ 共同研究成果 ]直接数値解析を用いた複合剪断流における亜臨界遷移現象の研究―直交した流れが局在乱流パターンに与える非線形相互作用―

松川 裕樹†・塚原 隆裕‡†東京理科大学大学院

理工学研究科 機械工学専攻
東京理科大学 理工学部 機械工学科‡

円管内や平行平板間などの壁乱流における亜臨界遷移現象のメカニズム・普遍則の解明に対しては理論的アプローチが困難で，従来は風洞や大規模な流路実験に頼ってきた．しかし近年の計算機技術の急速な進歩により，流れ場を数値的に取り扱う数値流体力学が発展するようになった．
本研究は大型ベクトル並列計算機（東北大学サイバーサイエンスセンター所有のスーパーコンピ
ュータ群 AOBA）を活用することで，同心二重円筒の内外円筒の周方向回転と軸方向圧力勾配駆動が合わさった複合剪断流（Taylor–Couette–Poiseuille 流）に対して大規模な直接数値解析を実施
し，亜臨界遷移現象の解明を目指した研究である．また，双方の駆動の大きさが拮抗した際に，従来観測されていた主流方向に対し傾斜した局在乱流構造とは異なる，輪状の局在乱流構造が現れることを発見した．

1.

序論

我々は日常生活を通してさまざまな形で流体の恩恵を受けている．例えば炊事・洗濯に用いる
生活用水に始まり，ヒートパイプや油圧機器といった工業製品として生活に利用しているほか，呼吸を通して空気を取り入れたり血液として酸素を体内に巡らせたりと，内的な利用も無意識のうちに行っている．そのため我々は古代から流体に親しみ生活に利用してきた．しかし，流体力学として確立されてきたのはここ数百年のことで，流体研究は従来，風洞実験や水路実験を用いて行われてきた．近年計算機技術の急激な進歩により流れ場を数値的に取り扱う数値流体力学（computational fluid dynamics, CFD）が発展し，流体現象の理解・予測・制御に大きく貢献しており，特に複雑な流体現象を再現するには膨大な格子数を要するためスーパーコンピュータ等の高性能な計算機が必要となる．流体流れは層流と乱流に大別することができ，層流では整った流れを呈し，その挙動は予測が比較的容易であるが，乱流では非線形性により流れが複雑となり，摩擦抵抗増大・伝熱特性向上といった特徴があるため工業的には乱流の促進・抑制が重要となる．
したがって乱流制御をする上で，層流・乱流の切り替わりである遷移域の研究は理学・工学の両側面から重要なものとなる．
壁面剪断流は線形安定性理論により求められる（これ以上では無限小撹乱に対して基本流が不安定となる）層流維持限界の線形安定臨界（局所安定臨界）レイノルズ数 ReL 以下の亜臨界遷移域でも乱流を維持することができ，時空間的な間欠性を持つ[1][2]．この間欠的な乱流構造は大域的安定臨界レイノルズ数 Reg まで維持することが可能であり，その形状から，乱流パフ（turbulent
puff），乱流縞（stripe），乱流帯（band）などと呼ばれる．例えばパイプ流は最も基本的な流路であ[3]
り，Reynolds の実験により層流・乱流・遷移状態のスケッチが行われたのは非常に有名である．
パイプ流は線形安定性解析によると ReL = ∞であり，非常に注意深く流体を流せば乱流へ遷移することはないはずである．しかし実際にはそのようなことは起きず，有限振幅撹乱に対して不安定となるため，線形安定性だけではなく非線形安定性も考慮する必要がある．パイプ流の亜臨界遷移域ではパフやスラグと呼ばれる局在的な乱流塊が現れることが知られており，パフはスラグに比べて発生するための擾乱が大きくなければならないが低 Re で発生し，スラグと異なり乱流塊が直接数値解析を用いた複合剪断流における亜臨界遷移現象の研究― 33 ―

拡大することはない．パフを十分長い流路内に発生させ，長時間観察すると分裂を繰り返す．しかしレイノルズ数が小さいとパフが消滅し完全な層流状態に至り，そこから自発的に乱流域が発生することはない．Avila et al.[4]はこのパフの分裂・消滅に着目し，生成されたパフが次の分裂までに要する時間とパフが消滅するまでの減衰時間を算出し，両者の観測結果を近似した二重指数関数が交わるレイノルズ数で統計的平衡状態が実現するとして Re = 2040 ± 10 を Reg と決定づけた．この発見は乱流遷移研究の大きなブレイクスルーとなり，Pomeau[5]が提唱した局在乱流構造の時空間的間欠性と非平衡統計物理における有向浸透現象（directed percolation, DP）との関連づけについて，改めて注目されるきっかけとなった．
このような亜臨界遷移の研究はパイプ流のような一次元系だけでなく二次元系でも報告されている．平行平板間内流れでは平面 Couette 流[6]や平面 Poiseuille 流[7]において主流に対して傾いたストライプ状の局在乱流構造（乱流縞）が観測される．それぞれの系の Reg に近づくにつれて乱流
縞はより断続的となり，局所的な乱流帯が不規則に分布した状態となる．この局在乱流の振る舞
いの臨界現象はパイプ流と同様 DP 普遍性に紐づけられる[8][9]．以上のことから，局在乱流の振る舞いを調べることは Reg 近傍の臨界現象の知見を深めるほか，Reg を決定するという点において重要であると言える．
亜臨界遷移の帯状の局在乱流は Taylor–Couette 流（Taylor–Couette flow, TCf）において Coles[10]が最初に発見した．TCf は同心二重円筒の内円筒と外円筒を周方向に回転させて流体を駆動する系
で，円筒の回転速度に応じてさまざまな構造が現れることが知られている[11]．TCf では，内円筒と外円筒の半径の比（以降は円筒比 η ≡ rin/rout と呼称）が大きい高円筒比において，内円筒と外円筒が対向回転すると亜臨界遷移となり，spiral turbulence（SPT）と呼ばれる乱流縞構造が現れる．
Andereck et al.[11]は η = 0.883 において実験により詳細な状態遷移図を作成し，以降の TCf 研究における重要な指標となっている．高円筒比 TCf は幾何的に 2 枚の平行平板をスライドさせて流体を駆動させる平面 Couette 流（plane Couette flow, pCf）に形状が近いため Prigent et al. [6]が実験によって pCf でも同様のストライプ構造が現れることを発見した．また，Tsukahara et al. [7]により平行
平板間に圧力勾配駆動を与える平面 Poiseuille 流（plane Poiseuille flow, pPf）においてもストライプ
構造が現れることが直接数値解析（direct numerical simulation, DNS）によって明らかになっている．
さらに大規模な DNS により pPf の乱流縞が Reτ ≤ 44 においても半永久的に維持することが可能であることを見出した[12][13]．さらに，pCf や pPf で現れる縞状／帯状構造の周囲に，大規模な二次流れが生じることがわかっている[12][14]．TCf のような環状流においても亜臨界遷移はよく研究されており，二重円筒の軸方向に圧力勾配を与えた環状 Poiseuille 流（annular Poiseuille flow, aPf）では η ≥ 0.5 において helical turbulence（HLT）と呼ばれる螺旋状の局在乱流構造が現れる[15]．SPT と
HLT はどちらも内円筒に巻き付いた局在乱流という点では共通しているが，TCf と aPf では主流
方向が周方向と軸方向で異なるため，乱流域を構成する縦渦が伸長する方向も異なっている．これら環状流は高円筒比（η → 1）で形状が pCf と pPf に近づくことから同様の縞構造が現れていると考えられる．しかしながら，これら二つの剪断流（Couette 流と Poiseuille 流）が直交する流れ場
において局在乱流パターンがどのように変化するかは不明である．
本研究対象の Taylor–Couette–Poiseuille 流（Taylor–Couette–Poiseuille flow, TCPf）は同心二重円筒間
内の内外円筒壁面周方向回転（TCf）と軸方向圧力勾配駆動（aPf）が組み合わさった複合剪断流である．このような複合剪断流では捩れた三次元境界層を形成している基本流となっており，壁面からの高さにより流れの方向が異なるため乱流構造が単一剪断流に比べて複雑になる．TCPf は実
験・数値計算を通して数多の研究がなされているがその多くは内円筒回転のみの条件を対象としており[16][17]，内円筒と外円筒が対向回転する亜臨界遷移を対象とした研究は数少ない[18]．本研究では軸方向に長い計算領域を持った TCPf に対して DNS を実施し，複合剪断流の亜臨界遷移現象を明らかにすることを目的としている．ちなみに DNS では，乱流に内在する主たる構造の全てを解像しつつ全容を十分に捉える必要があるため，微細な格子幅および広大な計算領域を要する．
特に，局在乱流が形成する空間パターンは大規模であるため，低レイノルズ数条件下とは言え，大負荷な計算実行を要するため，ベクトル型大型並列計算機の活用は不可欠である．

― 34 ―

SENAC Vol. 55, No. （4 2022. 10）

図 1 Taylor–Couette–Poiseuille 流の模式図．

2.

解析対象と支配方程式

本研究では作動流体を非圧縮性 Newton 流体とし，図 1 に示すような TCPf に対して DNS を実施した．乱流現象を忠実に再現できるよう，乱流モデルを一切用いていない．rin, rout, Ωin, Ωout はそれぞれ内円筒・外円筒半径，内円筒・外円筒回転角速度を表し，内円筒・外円筒回転速度はそれ
ぞれ uin = rinΩin, uin = routΩout となる．代表長さを両円筒間ギャップ幅 h = rout – rin，代表速度を両円筒の相対回転速度 u0 = uin – uout とすることにより，動粘性係数 ν を用いて
𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅in =𝑢𝑢𝑢𝑢in ℎ𝑢𝑢𝑢𝑢out ℎ
𝑢𝑢𝑢𝑢0 ℎ, 𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅out =, 𝑅𝑅𝑅𝑅𝑒𝑒𝑒𝑒 =
𝜈𝜈𝜈𝜈
𝜈𝜈𝜈𝜈
𝜈𝜈𝜈𝜈

(1)

と各レイノルズ数を定義する．支配方程式は，円筒座標系の連続の式と Navier–Stokes 方程式：
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑧𝑧𝑧𝑧∗ 1 𝜕𝜕𝜕𝜕(𝑟𝑟𝑟𝑟 ∗ 𝑢𝑢𝑢𝑢𝑟𝑟𝑟𝑟∗ ) 1 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜃𝜃𝜃𝜃∗
+
+ ∗
= 0,
𝜕𝜕𝜕𝜕𝑧𝑧𝑧𝑧 ∗ 𝑟𝑟𝑟𝑟 ∗ 𝜕𝜕𝜕𝜕𝑟𝑟𝑟𝑟 ∗𝑟𝑟𝑟𝑟 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑧𝑧𝑧𝑧∗𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕∗1 ∗2 ∗
1∗∗ )𝑢𝑢𝑢𝑢∗
(𝑢𝑢𝑢𝑢
+
⋅
∇
=
−
+
∇ 𝑢𝑢𝑢𝑢𝑧𝑧𝑧𝑧 +𝐹𝐹𝐹𝐹(𝑃𝑃𝑃𝑃),
𝑧𝑧𝑧𝑧∗∗
𝜕𝜕𝜕𝜕𝑡𝑡𝑡𝑡𝜕𝜕𝜕𝜕𝑧𝑧𝑧𝑧
𝑅𝑅𝑅𝑅𝑒𝑒𝑒𝑒
𝑅𝑅𝑅𝑅𝑒𝑒𝑒𝑒𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑟𝑟𝑟𝑟∗
𝑢𝑢𝑢𝑢𝜃𝜃𝜃𝜃∗2𝜕𝜕𝜕𝜕𝑝𝑝𝑝𝑝∗
1
𝑢𝑢𝑢𝑢𝑟𝑟𝑟𝑟∗
2 𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝜃𝜃𝜃𝜃∗∗∗ )𝑢𝑢𝑢𝑢∗
∗2 ∗
(𝑢𝑢𝑢𝑢
+
⋅
∇
−
=
−
+
�∇
𝑢𝑢𝑢𝑢
−
−
�,
𝑟𝑟𝑟𝑟
𝑟𝑟𝑟𝑟
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 ∗𝑟𝑟𝑟𝑟 ∗
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 ∗ 𝑅𝑅𝑅𝑅𝑒𝑒𝑒𝑒
𝑟𝑟𝑟𝑟 ∗2 𝑟𝑟𝑟𝑟 ∗2 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝜃𝜃𝜃𝜃∗𝑢𝑢𝑢𝑢𝑟𝑟𝑟𝑟∗ 𝑢𝑢𝑢𝑢𝜃𝜃𝜃𝜃∗
1 𝜕𝜕𝜕𝜕𝑝𝑝𝑝𝑝∗
1
𝑢𝑢𝑢𝑢𝜃𝜃𝜃𝜃∗
2 𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑟𝑟𝑟𝑟∗∗∗ )𝑢𝑢𝑢𝑢∗
∗2 ∗
(𝑢𝑢𝑢𝑢
+
⋅
∇
+
=
−
+
�∇
𝑢𝑢𝑢𝑢
−
+
�
𝜃𝜃𝜃𝜃
𝜃𝜃𝜃𝜃
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 ∗𝑟𝑟𝑟𝑟 ∗
𝑟𝑟𝑟𝑟 ∗ 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝑅𝑅𝑅𝑅𝑒𝑒𝑒𝑒
𝑟𝑟𝑟𝑟 ∗2 𝑟𝑟𝑟𝑟 ∗2 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕

(2)
(3)
(4)
(5)

である．式(3)右辺の外力項，軸方向一定圧力勾配駆動力の大きさを表す圧力勾配関数（定数）で，𝐹𝐹𝐹𝐹(𝑃𝑃𝑃𝑃) = −𝑅𝑅𝑅𝑅𝑒𝑒𝑒𝑒d𝑃𝑃𝑃𝑃∗
d𝑧𝑧𝑧𝑧 ∗

(6)と定義する．ここで，上付き添え字( )*は h または u0 による無次元化を表す．本研究では各レイノルズ数と圧力勾配関数 F(P)をコントロールパラメータとして解析を行った．また，本研究ではAndereck et al.[11]の研究に合わせ，円筒比を全て η = 0.883 で統一している．
式(2)–(5)のカップリングには Fractional-step 法を用いた．空間的離散化に有限差分法を用い，z方向と θ 方向には四次精度中心差分，r 方向は不等分割格子として二次精度中心差分を採用した．
時間的離散化については r 方向の粘性項に二次精度 Crank–Nicolson 法を，その他の項には二次精
度 Adams–Bashforth 法を用いた．計算領域は全ての条件で Lz × Lr × Lθ = 256h × h × 2πr であり，格子数は条件に応じて Nz × Nr × Nθ = 2048 × 64 × 256, 2048 × 64 × 512, 4096 × 64 × 256 としている．

直接数値解析を用いた複合剪断流における亜臨界遷移現象の研究― 35 ―

表 1 プログラム性能（SX-Aurora，1 ノード 8 コア並列時）．100,000 ステップ分の計算結果．
実行時間 [hh:mm:ss]
8:44:41
CPU 時間 [hh:mm:ss]
69:48:13
実行時間 1 秒あたりの MFLOPS 値
189,260
平均ベクトル長
255.4
ベクトル演算率
98.9 %MIPS [million instructions/sec]
4,127
CPU ポート競合時間 [sec]6487

図 2 純粋な aPf のギャップ中央 z–θ 面における半径方向速度 ur’/uτ,z の二次元可視化図：(a) Reτ,z =48, (b) Reτ,z = 52．[Matsukawa, Y. & Tsukahara, T., Phys. Fluids, 34, 074109 (2022); licensed under a CCBY license. Copyright © 2022 The Authors. Published by AIP Publishing.]
数値計算には東北大学サイバーサイエンスセンター所有のスーパーコンピュータ群 AOBA を利用した．使用した計算機は NEC 製のベクトル型スーパーコンピュータ SX-Aurora TSUBASA で，ノード内 OpenMP 並列化を施したコードでの検証である．表 1 に本研究で実行した計算における計算性能を示す．

3.

解析結果

3.1 Pure-aPf
TCPf の解析結果を示す前に，Couette 成分が存在しない純粋な aPf（pure-aPf）で実施した DNSの結果について議論する．DNS は軸方向摩擦レイノルズ数 Reτ,z ≡ uτ,zh/2ν = 48, 52 に対して，計算領域（Lz = 256h）を用いて亜臨界遷移現象を再現した．この計算領域は Ishida et al. [15]による aPfの DNS と比較しても十分大きく，大規模な間欠性を維持できる十分な格子解像度を使用している．また，初期条件を十分高いレイノルズ数の完全乱流場として，そこから目標のレイノルズ数に向けて長い計算時間をかけて準静的に変化させることにより層流・乱流共存場を再現している．
図 2 に pure-aPf の可視化結果を示す．z–θ 平面ギャップ中央で半径方向変動速度 ur’/uτ,z を可視化した二次元図で，緑色の箇所は速度変動が小さい層流域，赤と青の箇所は変動速度の大きいストリーク構造が入り乱れた乱流域を示している．図 2(a)の可視化図では Reτ,z = 48, (b)の可視化図で

― 36 ―

SENAC Vol. 55, No. （4 2022. 10）

図 3 Set A（Rein = 400, Reout = −1000）のギャップ中央 z–θ 面における半径方向速度 ur’*の二次元可視化図：F(P) = (a) 5, (b) 10, (c) 15, (d) 20．[Y. Matsukawa & T. Tsukahara, Phys. Fluids, 34, 074109(2022); licensed under a CC-BY license. Copyright © 2022 The Authors. Published by AIP Publishing.]

は Reτ,z = 52 を可視化しているがどちらにおいても局所的な乱流域（HLT）を確認できる．先行研
究[15]では一方向の乱流縞が確認されているが，これは本研究と比較して軸方向に短い領域（η ≥ 0.5
で Lz ≤ 74h）であったことが原因と考えられる．また，図 2(a)では乱流領域に端点を確認でき（z*≈ 100, 130）
，これは有限長の縞であることがわかる．

3.2 Set A
Set A（Rein = 400, Reout = −1000）は軸方向圧力勾配が存在しない状態（pure-TCf）の状態で層流状態を示し，特徴的な構造は現れないことが先行研究から明らかになっている[11]．図 3 は F(P)の増加に伴う流れ場の変化を示しており，ギャップ中央の z–θ 面内における半径方向変動速度 ur’* =ur’/u0 の二次元可視化である．図 3(a)では特徴的な乱流構造は現れず，完全な層流場となっている．
また，F(P) = 0 でも同様の可視化図となっており，先行研究[11]とも矛盾しない．(b)は F(P) = 5，(c)は F(P) = 10 の可視化図で，どちらも軸方向に長く伸びた局在乱流構造が現れていることが確認できる．(b)の状態における軸方向摩擦レイノルズ数は内円筒で Reτ,zin = 42.3，外円筒で Reτ,zout = 41.4
であり，これらの値は pPf の大域的安定臨界レイノルズ数 Reg[12][13]に近い値となっている．また，(b)では端点を持った有限長の縞であるのに対して，(c)では少なくとも当計算領域（Lz = 256h）では端点を確認できず無限長の縞となったことから，Reg 近傍においては縞が有限長となることが
pure-aPf の結果と合わせて確認できた．さらに，図 2(a), (b)と異なり，図 3(b), (c)では乱流縞の方向が一方向に揃っていることがわかる．Pure-aPf は軸方向の Poiseuille 成分のみを持ち，周方向のCouette 成分を持たなかった一方で，Set A の図 3(b), (c)は Couette 成分と Poiseuille 成分の両方を持つため，pure-aPf の持っていた対称性が Couette 成分によって破られたことで縞の向きが一方向に揃ったものと解釈できる．また，(b), (c)の状態では Couette 成分よりも Poiseuille 成分の方が支配的となっているため，これらの縞は HLT と見なせる．二点速度相関を基に算出した HLT の間隔
で，(c)の縞の z 軸正方向に対する傾き（α）を求めると α = 7°であった．(d)の F(P) = 20 では流れ場全体が完全乱流に近づいているものの，所々に間欠性が現れていることが確認できる．この構
造は環状 Couette 流（annular Couette flow, aCf）や偏心 aPf でも確認されており，laminar patch と呼ばれている[19][20]．
3.3 Set B
Set B（Rein = 800, Reout = −2000）は Andereck et al.[11]によると SPT が現れる領域に位置している．
図 4 は Set B のレイノルズ数に対して F(P)を大きくした際の流れ場の変化を示した可視化図で，
図 3 と同様，ギャップ中央 z–θ 面における半径方向変動速度を可視化している．図 4(a)は F(P) = 0の pure-TCf の可視化図である．
二方向を向いた乱流縞が流れ場全体に現れていることが確認でき，

直接数値解析を用いた複合剪断流における亜臨界遷移現象の研究― 37 ―

図 4 Set B（Rein = 800, Reout = −2000）のギャップ中央 z–θ 面における半径方向速度 ur’*の二次元可視化図：F(P) = (a) 0, (b) 2, (c) 4, (d) 6, (e) 8, (f) 10．[Y. Matsukawa & T. Tsukahara, Phys. Fluids, 34,
074109 (2022); licensed under a CC-BY license. Copyright © 2022 The Authors. Published by AIP
Publishing.]

先行研究[11]からこの縞は SPT であることがわかる．(b)は F(P) = 2 の可視化図で，(a)とは異なり，
SPT の方向が一方向に揃っていることがわかる．これは pure-aPf の対称性を Couette 成分が破ったのと同様，pure-TCf が持つ対称性を Poiseuille 成分が破ったことによるものであると解釈できる．
(c)は F(P) = 4 の可視化図で，(b)と同様一方向の SPT が現れたが縞の傾きが異なることがわかる．
Set A で現れた HLT は F(P) = 10, 15 の両方の条件で z 軸正方向に対して正の傾き（α > 0）を持っていたのに対して，図 4(b), (c)の SPT は負の方向（α < 0）に傾いていることがわかる．二点速度
相関を基に算出した SPT の間隔でこの角度を求めると，(b)では α = −52°であったのに対して(c)では α = −31°であった．(d)は F(P) = 6 の可視化図で，これまでの可視化図とは大きく異なり，乱流
縞が z 軸に対して垂直（α = 90°）に現れている．また，現れた乱流縞はそれぞれ独立しており，等間隔で軸方向に整列している．TCPf は環状流であるからこれらの垂直な乱流縞は内円筒に輪状に巻き付いているものであることがわかる．以降，この輪状の局在乱流を ring-shaped turbulence(RST)と呼ぶことにする．RST は軸方向の Poiseuille 成分が周方向の Couette 成分と拮抗する大きさとなったことにより現れたものと考えられる．(e)の F(P) = 8 ではまだ RST を確認できるものの，形状は(d)に比べて崩れている．しかし，この崩れた RST はこの形状を長時間保って下流に伝播し続ける．さらに(f)の F(P) = 10 では図 3(d)で見られたような laminar patch が現れた．

3.4 Set C
Set C（Rein = 1200, Reout = −3000）は Set B と同様，pure-TCf で SPT が現れるレイノルズ数域に位置している[11]．図 5 の可視化図は Set C のもので，見方は図 4 と同様である．(a)の F(P) = 0 は二方向の SPT が現れており，先行研究[11]の結果と矛盾しない．(b)の F(P) = 2 ではまだ SPT が二方向で現れており，(c)の F(P) = 4 でようやく一方向に揃っている．これは Set B の結果と異なっており，図 5(b)ではまだ Couette 成分の対称性を崩せるほど Poiseuille 成分が大きくないことが原因だと考えられる．さらに F(P)を大きくすると Set B のときと同様，(d)の F(P) = 6 で RST が現れた．
Set B で現れた RST と縞の間隔が異なるほか，図 5(a)–(d)の局在乱流が現れていない層流域ではノイズが入った流れ場となっていることも相違点として挙げられる．(e)の F(P) = 8 と(f)の F(P) = 10では軸方向に長く伸びた縞状構造が現れた．しかし，他の条件で観測された局所的な乱流構造とは異なり，層流領域と乱流領域の境界が不明瞭で，計算領域全体が乱れていることがわかる．

― 38 ―

SENAC Vol. 55, No. （4 2022. 10）

図 5 Set C（Rein = 1200, Reout = −3000）の z–θ 平面ギャップ中央における半径方向速度 ur’*の二次元可視化図：F(P) = (a) 0, (b) 2, (c) 4, (d) 6, (e) 8, (f) 10．[Y. Matsukawa & T. Tsukahara, Phys. Fluids, 34,
074109 (2022); licensed under a CC-BY license. Copyright © 2022 The Authors. Published by AIP
Publishing.]

3.5 Ring-shaped turbulence
Set B, C で見られた輪状の局在乱流構造，RST を二次元・三次元可視化したものを図 6 に示す．
図 6 で示したものは全て Rein = 800, Reout = −2000, F(P) = 6 で現れた RST である．(a)では RST を三次元可視化しており，赤が軸方向変動速度 uz’*の大きい領域，青が変動速度の小さいストリーク構造を示している．このような輪状の形状をした局在乱流パターンは他の系では確認されていない．
(b)–(d)は RST の z–θ 平面における軸方向変動速度を二次元可視化している（図 4(d)は半径方向変動速度の可視化）．(b)–(d)は可視化している位置が異なり，それぞれ(b)外円筒壁面付近（y*in ≡ r* −r*in = 0.9），(c)ギャップ中央（y*in = 0.5）, (d)内円筒壁面付近（y*in = 0.1）である．また，(b)外円筒
壁面付近と(d)内円筒壁面付近には壁面剪断応力の方向を矢印で併記している．TCPf は Couette 成
分と Poiseuille 成分が直交する複合剪断流であるため捩れた三次元的な速度分布を形成する．そのため壁面からの高さによって流れの方向が異なっている．また，図 6(b), (d)では軸方向変動速度を可視化しているため壁面付近で筋状に伸びたストリーク構造を確認することができる．これらのストリーク構造の向きは矢印で示した壁面剪断応力の方向と概ね一致している．したがって RSTを構成する縦渦構造（可視化しているのはストリーク構造だが低速・高速ストリークが並んでいることは縦渦構造がその場所に存在していることを示唆している）は高さによって向きが異なっているにもかかわらず，乱流縞としての大規模な構造自体は壁面からの高さによらず不変である
ことがわかる．これまで単一剪断流で確認されてきた乱流縞構造は主流に対して傾いて発現していたことを踏まえると，RST は円筒の軸に対して垂直ではあるものの，各高さにおける主流方向に対しては同様に傾いていると考えることができる．
3.6 流れ状態遷移図
前節までに，各レイノルズ数において F(P)に対する局在乱流パターンの依存性を調査し， 多様なパターンが存在することが判明した．改めて整理すると，Set A では有限長・無限長の一方向の乱流縞（HLT）が出現した．縞の長さは Poiseuille 成分の大きさに依存し，縞の方向選択性は Couette
成分が pure-aPf の対称性を破ったことにより変化したものである．Set B では F(P) = 0 の状態で Lz= 256h の軸方向に十分大きな計算領域を確保していれば二方向の SPT を捉えられ，軸方向圧力勾配の増加に伴い，一方向 SPT，RST を経由して laminar patch に至ることがわかった．二方向 SPTが一方向 SPT へと切り替わるのは Poiseuille 成分が pure-TCf の持つ対称性を破ることによるもの

直接数値解析を用いた複合剪断流における亜臨界遷移現象の研究― 39 ―

図 6 Rein = 800, Reout = −2000, F(P) = 6 における輪状局在乱流構造（ring-shaped turbulence, RST）の三次元可視化図と二次元可視化図．(a) 三次元可視化図（赤：uz’* ≥ 0.075，青：uz’* ≤ −0.075）．z–θ
面の軸方向変動速度の二次元可視化：y*in = (a) 0.9, (b) 0.5, (c) 0.1．[Y. Matsukawa & T. Tsukahara, Phys.
Fluids, 34, 074109 (2022); licensed under a CC-BY license. Copyright © 2022 The Authors. Published by
AIP Publishing.]
である．RST は内円筒に輪状に巻き付いた局在乱流構造で，それぞれの縞は独立して存在し，軸方向に等間隔に整列する．Set C は set B に似たパターン変化をするが，二方向 SPT から一方向に切り替わる F(P)の閾値が異なり，F(P) = 8, 10 では軸方向に長い縞状構造が現れるという特徴がある．これらの流れのパターンを，層流状態も含めて 8 種に分類し，図 7 の状態遷移図で整理した．

― 40 ―

SENAC Vol. 55, No. （4 2022. 10）

図 7 η = 0.883 の TCPf における流れの状態遷移図：(a) F(P) vs Re，(b) Reτ,z vs Reτ,θ．(a)と(b)はどち
らも Poiseuille 成分と Couette 成分の比を表している．下部に示したコンターは計算領域のギャッ
プ中央面全体における速度変動の瞬時場である．[Y. Matsukawa & T. Tsukahara, Phys. Fluids, 34,
074109 (2022); licensed under a CC-BY license. Copyright © 2022 The Authors. Published by AIP
Publishing.]

図 7 の(a), (b)はどちらも軸方向の Poiseuille 成分と周方向の Couette 成分の比を軸に取ったものだが，(a)では圧力勾配関数 F(P)と内外円筒の相対速度に基づくレイノルズ数 Re = Rein – Reout を比較
し，(b)では軸・周方向の摩擦レイノルズ数 Reτ,z, Reτ,θ を比較しているという点で異なる．状態遷移
図の全容と詳細で不明確な領域もあり，引き続き系統的調査が必要である．

4.

結言

本研究では η = 0.883 の Taylor–Couette–Poiseuille 流に対して DNS を実施し，亜臨界遷移現象を調査した．コントロールパラメータは軸方向圧力勾配関数 F(P)と，内円筒・外円筒の回転速度に基づくレイノルズ数 Rein, Reout である．TCf で層流となるレイノルズ数（set A: Rein = 400, Reout =−1000）と SPT が現れるレイノルズ数（set B: Rein = 800, Reout = −2000, C: Rein = 1200, Reout = −3000）に対して軸方向圧力勾配を付加することによって TCPf の局在乱流パターンを調査した．
TCPf において，RST 等を含む様々な流れパターンの存在が明らかになったが，当該系にはパラメータが多く存在しているため，本研究の結果はその一端を解明したものにすぎない．今後，パラメータをさらに変化させ，亜臨界遷移現象のメカニズム解明を目指すには高解像度の大規模計算を効率よく行う必要がある．引き続き東北大学サイバーサイエンスセンターとの共同研究を通して大規模科学計算システム AOBA の高性能スーパーコンピュータ群を積極的に活用することで本研究のさらなる発展が望まれる．
本研究結果の詳細については既報[21]を参照されたい．

謝辞
本研究は，東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができた．また，研究にあたっては同センター関係各位に有益なご指導とご協力をいただいた．本研究は JSPS 科研費（基盤研究(B) JP19H02071）の助成を受けて行われたものである．記して，感謝の意を表する．

直接数値解析を用いた複合剪断流における亜臨界遷移現象の研究― 41 ―

参考文献[1] P. Manneville, “Transition to turbulence in wall-bounded flows: Where do we stand?,” Mech. Eng. Rev.
3, 15-00684 (2016).
[2] L. S. Tuckerman, M. Chantry, and D. Barkley, “Patterns in wall-bounded shear flows,” Ann. Rev. Fluid
Mech. 52, 343–367 (2020).
[3] O. Reynolds, “An experimental investigation of the circumstances which determine whether the motionof water shall be direct or sinuous, and of the law of resistance in parallel channels,” Philos. Trans. R.
Soc. 174, 935–982 (1883).
[4] K. Avila, D. Moxey, A. De Lozar, M. Avila, D. Barkley, and B. Hof, “The onset of turbulence in pipe
flow,” Science 333, 192–196 (2011).
[5] Y. Pomeau, “Front motion, metastability and subcritical bifurcations in hydrodynamics,” Phys. D 23, 3–
11 (1986).
[6] A. Prigent, G. Grégoire, H. Chaté, O. Dauchot, and W. van Saarloos, “Large-scale finite-wavelength
modulation within turbulent shear flows,” Phys. Rev. Lett. 89, 014501 (2002).
[7] T. Tsukahara, Y. Seki, H. Kawamura, and D. Tochio, “DNS of turbulent channel flow at very low
Reynolds numbers,” in Proceedings of 4th International Symposium on Turbulence and Shear Flow
Phenomena, 935–940 (2005).
[8] M. Sano and K. Tamai, “A universal transition to turbulence in channel flow,” Nat. Phys. 12, 249–253(2016).
[9] G. Lemoult, L. Shi, K. Avila, S.V. Jalikop, M. Avila, & B. Hof, “Directed percolation phase transition
to sustained turbulence in Couette flow,” Nat. Phys., 12, 254–258 (2016).
[10] D. Coles, “Transition in circular Couette flow,” J. Fluid Mech. 21, 385–425 (1965).
[11] C. D. Andereck, S. S. Liu, and H. L. Swinney, “Flow regimes in a circular Couette system withindependently rotating cylinders,” J. Fluid Mech. 164, 155–183 (1986).
[12] 塚原，石田，“平面 Poiseuille 流の亜臨界遷移における下臨界レイノルズ数”，日本流体力学会
誌：ながれ 34, 383–386 (2015)．
[13] M. Shimizu and P. Manneville, “Bifurcations to turbulence in transitional channel flow,” Phys. Rev.
Fluids, 4, 113903 (2019).
[14] X. Xiao and B. Song, “Kinematics and dynamics of turbulent bands at low Reynolds numbers in channel
flow,” Entropy 22, 1167 (2020).
[15] T. Ishida, Y. Duguet, and T. Tsukahara, “Transitional structures in annular Poiseuille flow depending on
radius ratio,” J. Fluid Mech. 794, R2 (2016).
[16] J. M. Nouri and J. H. Whitelaw, “Flow of Newtonian and non-Newtonian fluids in a concentric annulus
with rotation of the inner cylinder,” J. Fluids Eng. 116, 821–827 (1994).
[17] M. Manna, A. Vacca, and R. Verzicco, “Reverse transition of a turbulent spiral Poiseuille flow at Ta =1500,” J. Fluid Mech. 941, A6 (2022).
[18] J. Langenberg, M. Heise, G. Pfister, and J. Abshagen, “Convective and absolute instabilities in counterrotating spiral Poiseuille flow,” Theoret. Comput. Fluid Dynamics, 18, 97–103 (2004).
[19] K. Kunii, T. Ishida, Y. Duguet, and T. Tsukahara, “Laminar–turbulent coexistence in annular Couette
flow,” J. Fluid Mech. 879, 579–603 (2019).
[20] K. Yatsushiro, Y. Horimoto, and T. Tsukahara, “Visualization of intermittent turbulent structure ineccentric annular pipe flow,” in Proc. the 13th Pacific Symp. on Flow Visualization and Image
Processing, 7 (2022).
[21] Y. Matsukawa and T. Tsukahara, “Subcritical transition of Taylor–Couette–Poiseuille flow at high radiusratio”, Phys. Fluids 34, 074109 (2022).