― 22 ―SENAC Vol. 54, No. （4 2021.10）［共同研究成果］

乱流パフの時空間欠性と DP 普遍性に関する一考察：
環状クエット・ポアズイユ流の直接数値シミュレーション
竹田 一貴†・森松 浩隆†・塚原 隆裕‡†

東京理科大学大学院‡
東京理科大学

理工学研究科 機械工学専攻
理工学部 機械工学科

乱流の発生・維持に関する根本的な問いの一つとして，壁乱流の亜臨界乱流遷移問題がある．
古くは 19 世紀より，パイプ流の乱流遷移が実験的に研究されてきたが，そのメカニズムは難解で，流体運動の非線形性ゆえに遷移の普遍則における境界条件依存性については理論的アプローチによる解釈も難しい．近年，大規模なベクトル型計算機の並列計算の高効率化が実現し，Navier-Stoke
式を直接解く直接数値シミュレーション（DNS）で，巨大流路実験をもってしても観測困難なほ
どの長時間にわたる乱流を追跡した観測が可能となった．
本研究では，パイプ流を模擬するような同軸二重円筒間で駆動方法が複合的なせん断流である
環状クエット・ポアズイユ流を対象として，乱流－層流共存の間欠パターン形成に関する内外円筒の半径比依存性やパターンの時空間欠性に関する普遍性を調査した．本解析では大規模領域を対象としており，局在乱流（乱流パフ）の確率的挙動を捉え，パイプ流における典型的な乱流パ
フとの比較，有向浸透現象（directed percolation, DP）の普遍性との特徴量を評価した．

1.

序論

流れには異なる性質を持つ「層流」と「乱流」があり，これらは流動の複雑具合によってエネ
ルギ散逸率が大きく異なる．流動していても渦の生成に乏しく定常的な流れとなる「層流」は流れの自由度が低く，制御が比較的容易であるが，エネルギ散逸率が低い．一方で，大小様々な渦が混在し絶え間なく新たな渦を生成する「乱流」状態はエネルギ散逸率が高い．さらに，物質の高い拡散性や活発な運動量・熱輸送が特徴である．この乱流状態は，工学的側面として原料を十分均一に混ぜ込む必要のある化学薬品生成，薬剤をはじめ化粧品や食品の攪拌，また高い伝熱特性を要求する熱交換器で利用されている．こういった利点の反面，乱流は強非線形・非定常で流動を予測しにくい．また，航空機の翼面で発生する摩擦抵抗を低減し省エネルギ化を目指すうえで，翼面上の流れを層流に保ち乱流化を防ぐ技術開発などは工学的に重要である．他にも，自然，
社会，産業においても流体の流動現象は切っても切れない関係にあり，気象予報，自然災害軽減，感染症拡大，製品性能向上において，流れの再現と予測が鍵となっている．近年，それらの解析
では数値流体力学（CFD）シミュレーションが大きな役割を果たしており，身の回りに溢れる空
気や水など様々な流体の流れについて我々は須らく CFD の恩恵を受けて生活している．この層流と乱流の切替わり，いわゆる乱流遷移に関する研究の歴史は古く，15 世紀の Leonardo da Vinci による渦のスケッチ（流路から流れ出た層流が池に流入し無数の渦を生成する様子），19 世紀のOsborne Reynolds によるパイプ流乱流遷移の実験的観測と定量評価を本格的な研究の始まりとしている．本研究は，乱流遷移の一つである亜臨界乱流遷移に着目したものである．
壁面せん断流の亜臨界乱流遷移では，レイノルズ数 Re（e.g., Rebulk = 管内平均流速 × 管直径 ÷
作動流体動粘度）が，自発的な層流の乱流化を考える線形安定性理論から求まる臨界値 ReL よりもはるかに低いにもかかわらず，乱流維持することが知られている．さらに，この亜臨界遷移域
において，特に完全な層流化を引き起こす直前のレイノルズ数で，特徴的な層流と乱流の共存構
造（いわゆる大規模な乱流間欠構造）を形成することが分かっている．この共存構造の形成は低乱流パフの時空間欠性と DP 普遍性に関する一考察：
環状クエット・ポアズイユ流の直接数値シミュレーション― 23 ―いレイノルズ数において起こるにも関わらず，層流状態に比べ高い伝熱特性・物質拡散性を有しており，亜臨界遷移状態の積極的な利用が期待されている．例えば，前述の馴染みあるパイプ流
では，流れの駆動方向に層流と乱流が間欠的に分布する 1 次元間欠性として，「乱流パフ」が形成される[1]．しかし，これら局在乱流は時空間的に大規模なスケールで現れるため，その形成メカ
ニズム，時空間欠性や遷移の普遍性については，理論的なアプローチが難しく未解明問題の一つであった．近年では，並列計算技術や大型ベクトル並列計算機の利用によって，乱流モデルを一切用いない直接数値シミュレーション（direct numerical simulation, DNS）を超大規模スケールで行うことが可能になり，Hof グループを中心に大規模な実験と DNS が行われ，乱流パフの分裂と減衰に要する時間のレイノルズ数依存性を発見し，その大域安定性から永久的に乱流維持可能な下限値の下臨界値 Reg（= 2040）を見出した[2,3]．また，チャネル流（平行平板間流）でも Reg 近傍
では，チャネル幅の数百倍のスケールを有した帯状の局在乱流，それらが複数整列したストライ
プ構造が形成される[4,5]．他のせん断流でも，それぞれ特徴的な局在乱流構造が確認され始めた[6-8]．さらに，これらの局在乱流はいずれも，決定論的な支配方程式に従いながらも確率的な分裂減衰挙動を呈し，一度完全に層流化すると自発的に乱流化することはないという吸収状態転移の観点から，有効浸透現象（directed percolation，DP）との関連性が活発に議論されている．いく
つかのせん断流では，それぞれ特徴的な局在乱流構造が確認されると共に時空間欠性の統計量評
価がなされており，Pomeau [9]が予言した DP 普遍クラスとの関連性が実証されつつある[10–16]．
具体的には，乱流パフの様に間欠構造が空間 1 次元的に広がりを持つ(1+1)-D DP，または空間 2 次
元的な(2+1)-D DP の特徴が報告されているが，2 つの DP 普遍クラスを決定する境界条件依存性などはいまだ不明である．
同軸二重円筒間の環状流路において，亜臨界遷移で形成される局在乱流構造の流路形状依存性，流れの駆動方法（クエット／ポアズイユ）依存性は興味深く，流路形状を決定する幾何パラメー
タの円筒比 η = rin / rout [0:1]によって間欠構造は次のように整理される．低円筒比（η ≤ 0.3）の流路では乱流パフに準ずる 1 次元間欠構造が，高円筒比（η ≥ 0.5）ではストライプ構造に準ずる空間 2 方向に間欠性を持つ螺旋構造が形成される [14,16–20]．つまり，円筒比に応じてパイプ流もしくはチャネル流における亜臨界遷移の特徴が，この環状流路において表出する．この系を研究対象とすることで，チャネル流～パイプ流の繋がり（もしくは不連続性）が調査できる．しかし，
環状クエット流と環状ポアズイユ流のいずれも，円筒比 η → 0 で厳密にパイプ流とは一致せず，
チャネル流～パイプ流を滑らかに接続して考えることは不可能であった．そこで，本研究対象として，パイプ流の基本流を模擬するべく，環状クエット・ポアズイユ流に設定した．
本研究では，環状クエット・ポアズイユ流の亜臨界遷移域の下臨界値（大局的臨界値 Reg）近傍において，クエット成分・ポアズイユ成分の強さを調整することでパイプ流の流速分布を模擬したうえで，局在乱流（乱流パフなど）の形状，挙動に関する流路円筒比依存性を調査した．また，大規模計算機を活用することで，数値計算特有の周期境界条件の影響を極力排除すべく，可能な限り巨大な主流方向領域を設定し，乱流パフの時空間欠性について，より大域的に起こる局在乱流の広がりと維持に関する本質と DP 普遍性の関連性に迫った．

図１ 環状クエット・ポアズイユ流の模式図と座標系．本研究では 2 つの駆動方向は一致する．

― 24 ―

2.

SENAC Vol. 54, No. （4 2021.10）

計算対象と方法

非圧縮性ニュートン流体の環状クエット・ポアズイユ流（図 1）が解析対象であり，軸方向（x）には周期境界を課して無限の管路を表現している．ただし，計算領域は有限長であり Lx = 409.6h
または 1500h とした．本解析系は同軸二重円筒系で，軸方向への内円筒のスライドと圧力勾配によって駆動する複合的な壁面せん断流である．また，流路形状を決定する幾何パラメータである
円筒比 η（= rin / rout）について，1 の極限をとると無限に広いチャネル流，0 の極限をとると太さ
無限小の内円筒を外円筒の中心に持つパイプ形状流路となる．前報[21]では内円筒のスライド駆動のみであったため，η → 0 でも流速分布はパイプ流のそれには漸近しない．しかし，環状クエット・ポアズイユ流では，流れの駆動方法のクエット成分とポアズイユ成分（圧力勾配）の強さを調整することでパイプ流と同等の流れを解析することができる．本報告では，パイプ形状に近い
低円筒比 η = 0.1，0.3 の流路における結果を紹介する．
流体の速度 u = (ux, ur, uθ) と圧力 p に関する支配方程式は，円筒座標系の連続の式及び Navier–Stokes 方程式である．式 (2) の最終項はポアズイユ成分，つまり圧力勾配駆動を表す外力項であり，Fp に関する圧力 P は場に依存しない定数である．また，式 (1)–(3) の上添え字 * は代表速度の内円筒スライド速度 Uw，内外円筒の半径差（ギャップ幅）h = rout − rin による無次元化を表す．
∇∗ � �∗ = 0,∂� ∗ �∗ + ��∗ � � ∗ ��∗ = −∇∗ �∗ +

(1)
��
1 ∗ ∗
∆� +
� ,
��� �
���

(2)

d�∗
ℎ� d�
(3)
�
��
=
−
.
�
d� ∗
��� d�
ここで，t は時間，μ は粘度である．レイノルズ数は，動粘度 ν を用いて Rew = Uwh / ν と定義した．
これは，平面クエット流での一般的な定義の 4 倍の値となる[22]．
�� ≡ −

基本流について，クエット成分とポアズイユ成分の調整は Fp で制御する．例えば，基本流である層流状態の流速分布に限って図 2 に紹介する．このとき，層流流速分布が Hagen-Poiseuille flow（HPf）と一致する（内円筒壁面で速度勾配がゼロとなる）には，Fp が次式を満たす必要がある．
�1 − ���
(4)�� ��� ��� = −4 �� �1 − 2 ln �� − 1
円筒比を問わず，その層流解が図 2 (a,b) の黒破線のような HPf の流速分布として得られている．
図 2 の環状クエット・ポアズイユ流（aCPf）の層流状態における主流方向流速分布は，クエット
成分（紫色）とポアズイユ成分（緑色）の線形和で表現される．しかし，本研究では乱流を対象としているため，その非線形性から亜臨界遷移で見られる乱流－層流共存構造についても，単純に環状クエット流（aCf）と環状ポアズイユ流（aPf）の要素の線形和として理解することはできない．
また，亜臨界遷移下での流れは層流と乱流が共存しているため，Fp の値は図 2（d 実線）よりも高く設定してある．例えば，Rew = 1600 の場合，Fp = 6.5 を採用しており，内円筒壁面で摩擦抗力係
数は Cf in = 0 の要件を厳密には満たしていない．ただし，内円筒近傍での乱流生成が支配的でない観点からすると，ナイーブな問題にはならない．
本解析は一切の乱流モデルを用いない直接数値解析（direct numerical simulation, DNS）を採用した．境界条件は円筒壁面で滑り無し条件，軸・周方向には周期境界条件を課した．本 DNS では，スタガード格子を採用し，周方向に 256 分割の格子数，主流方向の格子数は 4096 もしくは 16384，
r 方向には不等分割格子にて 64 分割とした．時間刻み Δt は CFL 条件を満たす可能な限り大きな値を用いている．式 (1)，(2) のカップリングには Fractional-step 法を用いた．時間進行は粘性項
半径方向成分について 2 次精度 Crank–Nikolson 法を，その他の項に 2 次精度 Adams–Bashforth法を用いた．空間的離散化には有限差分法を用い，主流及び周方向に 4 次精度中心差分を，半径方向に不等分割格子の 2 次精度中心差分をそれぞれ適用した．非圧縮性流体における圧力ポアソ

乱流パフの時空間欠性と DP 普遍性に関する一考察：
環状クエット・ポアズイユ流の直接数値シミュレーション― 25 ―

図 2 環状クエット流（aCf）と環状ポアズイユ流（aPf）の層流流速分布の線形和で表される環状クエット・ポアズイユ流（aCPf）の流速分布．(a) 低円筒比流路 η = 0.1，(b) 高円筒比流路 η= 0.9 では圧力勾配関数 Fp の適切な設定により，円管内流れの層流流速分布（Hagen-Poiseuilleflow, HPf）を忠実に再現している．HPf の模擬だけでなく，例えば，低円筒比流路 η = 0.1 で Fp
= −2 とすると，(c) の様に外円筒側（y* > 0.5）で逆流するような基本流を得ることもできる．
(d) は環状流路において HPf を模擬する Fp HPf の円筒比依存性（式 (4)）と純粋なクエット流（Fp= 0）を表している．
ン方程式には，高速フーリエ変換で波数空間にて解く解法を採用しているが，MPI 並列計算における計算速度のボトルネックとなりやすい箇所である．MPI，OpenMP ハイブリッド並列計算の並列数増加に伴う損失については第 4 章で報告する．
本解析対象における線形安定の臨界値 ReL は，低円筒比（η ≤ 0.1）においても流れ場が HPf と一致する場合は任意のレイノルズ数に対し線形安定（ReL → ∞）である．さらに，流れ場が HPf と一致しない場合でも，本解析で扱うレイノルズ数 Rew = O(103) と比較して ReL は極めて高い[23]．よって，いずれの流路系に対しても Rew ≪ ReL の，まさに亜臨界レイノルズ数領域の調査であるため，自発的に層流から乱流へ遷移することはない吸収状態転移を考えることになる．よって，一連の DNS 解析にて，最初に比較的高い Rew および初期擾乱を与えた場のシミュレーションを実施
し，十分発達した乱流場を初期条件として準静的に Rew および Fp を下げていき，大域的層流化に至るまでを調査した．また，亜臨界特有の乱流－層流共存場から成る時空間欠性の統計的な特徴を得るためにレイノルズ数の減少後には十分時間をかけ，主に統計的定常状態を議論する．ただし，臨界値近傍での局在乱流の動的挙動に関する議論についてはこの限りではなく，レイノルズ数を急減させた際の乱流の減衰，いわゆる臨界クエンチングを行った結果について議論する．

― 26 ―

SENAC Vol. 54, No. （4 2021.10）

図 3 環状クエット・ポアズイユ流における乱流間欠構造．同軸二重円筒の，ある(x-r)断面上半分を半径方向速度変動で可視化してあり，黒が乱流，白が層流を示す．Rew = 1600, Fp = 6.5 でパ
イプ流の流速分布を模擬した低円筒比（η = 0.1）において，典型的な軸方向に延びた赤破線で示す釣り鐘型の乱流パフを形成する．ただし，見やすさのため半径，主流方向長さは縮尺が異なる． 時間経過で乱流パフは移流すると共に，減衰（青矢印）と維持（赤矢印），分裂を繰返す．

3.

結果と考察

本章では，大規模 DNS で得られた結果として，環状クエット・ポアズイユ流で形成される大規模間欠構造を紹介する．特に，低円筒比流路で形成される 1 次元間欠性に対して，過渡的・統計的な時空間欠性に関するいくつかの最も基本的な冪的スケーリングの観点から DP 普遍性との関連を議論する．ただし，以降に示す本結果から，すべての壁乱流の遷移現象について DP 普遍性の有無を議論するものではなく，固有の特徴に関する境界条件依存性を抽出しているものといえる．

3 . 1 環 状 クエ ッ ト ・ ポ ア ズ イ ユ 流 で 形 成 さ れ る 大 規 模 間 欠 構 造
過去の筆者らの研究で，環状クエット・ポアズイユ流においても，環状クエット流と環状ポア
ズイユ流で形成される大規模な層流－乱流共存構造の円筒比依存性と同等の流路形状依存性を持つことが確認された[14,16–20]．例えば，円筒比 η = 0.8 の高円筒比流路では，内円筒に巻き付く螺旋状の局在乱流，いわば，チャネル流におけるストライプ構造に準ずる大規模間欠構造を形成した．螺旋やストライプ構造は，主流方向と周（スパン）方向に対し間欠的で，2 次元パターンを呈する[5]．一方，低円筒比流路（η ≤ 0.3）では，図 3 に示すような複数の局在乱流（乱流パフ）を形成する．この局在乱流は周方向に間欠性がコヒーレントであり，層流・乱流共存の繰り返しは
主流方向のみに限られた 1 次元間欠性をもつ．したがって，この局在乱流構造はまさにパイプ流における 1 次元間欠構造を構成する乱流パフの定性的特徴[2,3]と一致している．
環状クエット・ポアズイユ流における乱流パフの主流方向の長さは，η ≤ 0.3 では円筒比に依らず，いずれの流路系においても大きさは Lx = 30h で概ね同じある．これは，同様な流路系である
環状クエット流では円筒比が小さくなるにつれて，乱流パフの長さが短くなる結果と一致しない[14,16–20]．さらに，乱流パフの分裂・消滅の挙動についても流れの駆動方法に大きく依存する．
図 3 は，η = 0.1 の環状クエット・ポアズイユ流内の乱流パフの時間発展を示す．アンサンブル平
均した内円筒壁面での摩擦抗力係数は<Cf in> ≈ 0 を満たすように，コントロールパラメータを Rew= 1600, Fp = 6.5 と設定しパイプ流の流速分布を模している． 乱流パフは主流方向に概ねバルク速度と一致した移流速度 Uf で移動しながら，分裂・減衰を大域的に繰り返す．この際，臨界値 Reg以上であり半永久的に間欠性のある乱流を維持する．
本研究では，環状流路で形成される乱流パフの特徴を定量的に調査する趣旨のもと，以降 η =0.1，0.3 の結果に着目して，その時空間欠性に関する議論を進める．

乱流パフの時空間欠性と DP 普遍性に関する一考察：
環状クエット・ポアズイユ流の直接数値シミュレーション― 27 ―

図 4 臨界値近傍（Rew ≈ Reg）の環状クエット流（aCf）と環状クエット・ポアズイユ流（aCPf）
で形成される局在乱流（乱流パフ：黒色）の時空間ダイアグラム．円筒比は，(a,b) η = 0.1，(c,d)η = 0.3．いずれも横軸は主流方向の空間を表しており，乱流パフと一緒に移動するような座標系である．乱流パフの移動速度 Uf は概ねバルク速度と一致している．また，縦軸は時間経過を示しており，乱流パフの挙動を追跡している．(a) の赤矢印は，時空間欠性の指標となる層流間隔（乱流パフ間距離）ΔT，ΔLx を示す．

3 . 2 乱 流 パフ の 時 空 間 欠 性 の 統 計 量 解 析
本節では，低円筒比の環状クエット流と環状クエット・ポアズイユ流で形成される乱流パフの
時空間欠性について定性的，定量的に議論を行う．時空間欠性については，図 4 に示す乱流パフの時空間ダイアグラムを用いて説明する．時空間ダイアグラムは，主流方向空間と時間の 2 軸を持つ．ただし，乱流パフの移流を考え，空間軸は乱流パフと共に動く移動座標である．つまり，主流方向に間欠的に配列される乱流パフの分裂・合体・減衰の時間進展を観測することになる．
このとき，乱流パフの主流方向への移流速度 Uf は，概ね式（5）に示すバルク速度 ubulk と一致している．図 4 の描写に関する詳細な定義は既報[13]を参照されたい．
����� =

���� �� ���� ��1�� �� �� ��, �, �, ��d�d�d�d����� ℎ2� �
�
�
���

(5)

まず，図 4 を用いて，それぞれの流れで起こる分裂・合体・消滅の特徴を定性的に比較する．
環状クエット流では円筒比の低下により乱流パフの空間長さが短く，より明確に分裂の頻度が増えており，さながら(1+1)-D DP の時空間欠性に酷似してくる[16]．また，図には掲載していないが，円筒比 η = 0.1 の環状クエット流では，その確率的挙動のレイノルズ数依存性が Lemoult etal.[10]のクエット流における 1 次元間欠性のものと酷似する．一方，環状クエット・ポアズイユ流
では，乱流パフの空間長さは円筒比依存性をあまり持たないものの，分裂挙動については，η = 0.1に比べ η = 0.3 の流路におけるパフの方が活発な分裂がより確率性を呈していると言える．パイプ流では分裂と減衰の頻度は，臨界点から遠ざかるにつれて二重指数関数的に増加すると言われている[2]．図 4 に示す各流れの臨界値はそれぞれ異なり，環状クエット流では円筒比が小さい順に
Reg = 1628，1555 [16]，環状クエット・ポアズイユ流では Reg = 1540，1465 付近である[24]．この点で，図 4 の各間欠性を形成したレイノルズ数と臨界値との距離（|Rew − Reg| / Reg）がそれぞれ異なるため分裂減衰頻度が異なることは自明であるが，特に η = 0.1 の環状クエット・ポアズイユ
流における臨界値近傍での乱流パフは長時間分裂することがなく，非(1+1)-D DP 的な現象といえる．定量的な比較評価は環状クエット・ポアズイユ流のもの（図 4（b,d））に限り図 5 で議論する．

― 28 ―

SENAC Vol. 54, No. （4 2021.10）

0(b) 10-1(a) 10

Re w = 1600
1575
1550
1540
(1+1)-D DP**

PDF ( L x )-3

PDF ( T )

-2

10
10

-4

10

0

10

1

10

L x

*

1

2

10

-4

4

10

10

Re w = 1500
1475
1470
1465
(1+1)-D DP

-2

10

3

10
T *

-1

*

-3

10

0

10

10
PDF ( T )

-2

-6

0
(d) 10

Re w = 1500
1475
1470
1465
(1+1)-D DP

10
*

-4

10

10

2

10

-1

PDF ( L x )

-3

10

-5

(c) 10

10

-3

10

-4

10

-5

10

-5

10

-6

10

-6

10

-2

10

10

-5

10

Re w = 1600
1575
1550
1540
(1+1)-D DP

-1

10

0

10

1

10

2

L x

*

10

-7

10

0

10

1

10

2

10

3

T

*

10

4

10

5

10

図 5 環状クエット・ポアズイユ流における乱流パフの時空間欠性に関する有効浸透現象（DP）との関連性の可否：(a, b) 円筒比 η = 0.1，(c, d) η = 0.3．統計的定常状態において，図 4 (a)の赤矢印で定義される層流間距離 ΔT，ΔLx の確率密度について時空間平均をとった分布である．レイノルズ数は低くなるにつれて Reg に近づく．黒破線はまさに臨界点における(1+1)-D DP モデ
ルによる理論的な空間・時間間欠性の確率密度分布で，いずれも冪則に従う．だだし，空間・時間のスケーリング指数はそれぞれ，μ┴1D = 1.748，μ║1D = 1.84 である．
図 4 で時空間欠性が観測されたケースについて，代表速度をバルク流速 ubulk，代表長さを水力
等価直径（= 2h）とするバルクレイノルズ数 Rebulk に換算すると，パイプ流でのレイノルズ数とよく比較できる．環状クエット・ポアズイユ流において η = 0.3 では，臨界値近傍（Rew ≈ Reg）の値
が Rebulk = 2145–2300 程度となり，η = 0.1 では Rebulk = 2045–2190 程度となる．円筒比を低下させると，パイプ流で乱流パフが形成されるレイノルズ数，臨界レイノルズ数（Reg,bulk = 2040）と概ね一致するようになる．また，円筒比 η = 0.1 に対し環状クエット流で Reg,bulk ≈ 1000，環状ポアズイ
ユ流で Reg,bulk ≈ 3000 となることが知られており[16,24]，この観点からも低円筒比の環状クエット・
ポアズイユ流によるパイプ流模擬の妥当性が示された．
次に，時空間欠性については乱流パフ同士の距離もしくは時間（層流間隔）について，それぞれ長い／短い層流間隔の含まれる存在割合，いわば確率密度分布を議論する．層流間隔は臨界値以上（Rew > Reg）のケースに限り，図 4(a)の間欠構造内にある赤矢印 ΔT，ΔLx で定義されており，確率密度について統計的定常状態において約 T* = O(104)の長時間平均を行った．層流間隔は臨界値以上で ΔL ~ exp(−ΔL/ξ) に従うことが古くから知られている．さらに，DP 普遍性に従う時空間
欠性の場合，臨界値に近づくにつれ ξ → ∞となるように，図 5 の黒破線で示すような冪則でスケーリングされるようになることが理論的解析で分かっている．また，時定数 ξ は相関長と呼ばれ，
ξ ~ |Rew − Reg|/Reg でスケーリングされることも知られているが，今回は取り扱わない．層流距離の確率密度分布について，η = 0.1 の空間間欠性は図 5 (a)のように Reg 近傍でも DP に従わず，黒乱流パフの時空間欠性と DP 普遍性に関する一考察：
環状クエット・ポアズイユ流の直接数値シミュレーション(a) Lx = 409.6h

0― 29 ―(b) Lx = 1500h

0

10

10

Sustain

-1

10

1

10

Re w = 1575
1550
1540
1530
1525
(1+1)-D DP

Ft

Ft

Re w = 1575
1550
1540
1530
1525
(1+1)-D DP

-1

10

2

10

3

t

*

10
= t Uw / h

4

10

1

10

2

10

3

t

*

10
= t Uw / h

410

図 6 円筒比 η = 0.1 の環状クエット・ポアズイユ流で，レイノルズ数を Reg 近傍に急減させた
際の局在乱流の減衰過程，いわゆる臨界クエンチング．初期状態はいずれも同じ間欠的な乱流状態である．臨界値以上では乱流維持し，乱流間欠率 Ft は統計的定常値を示す．一方，Reg 以
下では，指数関数的に減衰し有限時間内に層流化する．黒破線は(1+1)-D DP モデルによる理論
的な Ft の減衰線であり，まさに Rew = Reg では時間の冪乗則（Ft ~ t−α） に従う．指数 α の値は(1+1)-D DP 固有の値 α1D = 0.159 を用いており，乱流パフの挙動が(1+1)-D DP 的であることを示唆している．また，領域サイズが大きい(b)では Ft のゆらぎが小さくなる．
矢印で示すような特定の長距離間隔（ΔLx* ≈ 100，160）を顕著に多く含む．乱流パフが特定の間隔で配列する原因として主流方向領域の不足の可能性も排除できないが，非 DP 的な時空間欠性は特出すべき点である，ただし，無限長さを持つパイプ流の実験においては，非 DP 的でありながらも特定の長距離間隔を持つことはなく，指数則に従うことが分かっている[3]．一方，η = 0.3 では
Reg に近づくにつれ，長距離の層流間隔についても DP に従うようになる典型的な空間間欠性である（図 5 (c)）
．時間間欠性は，どちらの円筒比，いずれのレイノルズ数においても，短–長時間の層流間隔に対して DP に従う傾向が一致している（図 5 (b,d)）．しかし，明確なレイノルズ数依存性が見えていないことから，領域を拡大した系を対象に，さらに長時間の観測が必要であることも示唆されたところである．

3 . 3 乱 流 パフ の 動 的 臨 界 特 徴 に つ い て － 臨界 ク エ ン チ ン グ の 観 点 か ら
本節では，局在乱流の減衰過程について調査することで，臨界値近傍での動的特徴を議論する．
臨界現象を測る基本的な手法である臨界クエンチングを行った．臨界クエンチングはある乱流場から臨界値近傍にめがけてレイノルズ数を急減させた際の過渡的な乱流減衰を議論するものである．今回は領域内に占める乱流の割合，いわゆる乱流間欠率の時間変化を測定した．
図 6 は円筒比 η = 0.1 の環状クエット・ポアズイユ流での臨界クエンチングによる過渡減衰のレ
イノルズ数依存性を乱流間欠率 Ft の時間変化として示しており，(a,b) は主流方向の計算領域が約 4 倍異なる．また，両対数表示しているので，横軸は右に進むにつれて指数的に長時間化しており，時間が経つにつれ，Ft の変化に要する時間は指数的に増加する．初期状態はいずれも流路一面に乱流パフが広がった同じ状態であり，そこから減衰し臨界値以上では Ft が統計的定常状態に収束し，臨界値以下では指数的に減衰し層流に至る．また，まさに臨界値では，Ft の減衰が時間の冪乗（Ft ~ t−α）に従っており，指数 α の値は(1+1)-D DP 固有の値 α1D = 0.159 であった．したがって，非 DP 的な空間間欠性を持ちながらも DP 的な減衰を起こすことが示唆された．また，図6 (a) では Reg 近傍の Rew = 1525–1550 でゆらぎが顕著に大きいが， 図 6 (b) ではより多くの乱流パフを捉えることができるため，統計的により良い Ft の減衰を捉えている．ただし，計算領域の増大は計算負荷の増大に直結するため長時間の解析が困難になる．本解析では SX-Aurora 8VE を用いた MPI，OpenMP ハイブリッド並列計算により実現できた．

― 30 ―

SENAC Vol. 54, No. （4 2021.10）

以上の観点から，低円筒比の環状クエット・ポアズイユ流における乱流パフを伴った亜臨界遷移について，その時空間欠性から部分的に DP 普遍性に従う特徴を抽出することができた．しかし，DP モデルの間欠性が持つ確率的な分裂・減衰挙動から乖離した時空間欠性からも DP 的な冪則が現れており，さながら擬 DP 的な要素も抽出されたといえる．この点で，DP 普遍クラスとの部分的な一致については，環状流，ひいては壁面せん断流固有の特徴が抽出されたともいえる．
時空間欠構造の普遍性に関する特徴のより包括的な理解においては，流れの境界条件を変えた解析を実施する必要がある．また，厳密な DP 普遍クラスとの評価には，さらにいくつかの時空間ス
ケーリングが必須になるため，本解析に比べ大領域・長時間の DNS が必要になってくる．

4.

計算性能

本研究では， MPI，OpenMP による並列化プログラミングを施した自作の Fortran コードを用いて高負荷の DNS を実施している．また，チューニングによってベクトル化効率を向上させてある（表 1）．よって，SX シリーズなどのベクトル型計算機でその実行効率は発揮されるが，以下に新たに導入された SX-Aurora TSUBASA を用いた計算性能を記す．また，旧機 SX-ACE で実施した同等の負荷となる数値解析との比較も行う．
本研究の解析条件の一つとして格子数 Nx × Nr × Nθ = 16,384 × 64 × 128 で 100,000 ステップ分のプログラム性能を抜粋して表 1 に示す．このとき，東北大学サイバーサイエンスセンターの SXAurora にて実行し，ノード内 8 コア OpenMP 並列化および 8 ベクトルエンジン（VE，いわゆるノード）間 MPI 並列化を施したコードでの検証結果である．また，前述の計算条件とは異なるが，同一コードによる DNS を東北大学サイバーサイエンスセンターSX-ACE，SX-Aurora を用いた同計算で比較検証した．ただし，本研究で使用した MPI 化 DNS コードでは，x-z 面において FFT を行う都合上，y 方向に要素を各ノードに対して割り振っている．そのため，最大並列ノード数は y方向格子数に依存しており，本解析では Ny = 64 であるため，最大並列は 64VE である．まず，並列運用時の性能として特筆すべき点は表 2，図 7 に示すように，本解析で最大規模の解析（8,192
× 64 × 4,096）に対して，MPI による使用メモリの増大は顕著であるものの，計算性能については
MPI による顕著な時間ロスは発生していない点である．さらに，SX-Aurora と SX-ACE を比較して，3 倍以上の高速化が成し遂げられた．
SX-Aurora TSUBASA 1VE 内の OpenMP 性能については既報[21]を参照されたい．

26

1024 × 64 × 1204 (= 2 )272048 × 64 × 1204 (= 2 )282048 × 64 × 2048 (= 2 )294096 × 64 × 2048 (= 2 )304096 × 64 × 4096 (= 2 )318192 × 64 × 4096 (= 2 )

350
300
250

2

200

×

150

Real Time26

 Ntot. grids [sec.]

400

100
50
0

1

10

100

Number of VE

図 7 MPI，OpenMP ハイブリッド並列計算性能．並列 VE 数を 1–64 ノードにした際の，1,000
ステップにかかる時間．縦軸の実計算時間は，解析対象の総格子数を計算規模として規格化してある．4VE 以上ではいずれの計算ケースでも，概ね同じ並列計算性能が発揮されている．

乱流パフの時空間欠性と DP 普遍性に関する一考察：
環状クエット・ポアズイユ流の直接数値シミュレーション

表1

プログラム性能（SX-Aurora TSUBASA）．100,000 ステップ分の計算結果．
実行時間 [hh:mm:ss]
11:14:01
平均ベクトル長
CPU 時間 [hh:mm:ss]
89:37:09
ベクトル演算率9
総 Flops 値（real）
36.2 × 10
MOPS (real)― 31 ―

212.2
94.7%
60,173

表2

主要な MPI, OpenMP ハイブリッド計算のプログラム性能．1,000 ステップ分の計算結果．
使用機
格子数
VE 数
コア数／VE 数 総メモリ
実行時間
1
16GB
240 sec2
19GB
210 sec
4
24GB
215 sec
1024 × 64 × 10248834GB
150 sec
16
54GB
120 sec32
98GB
100 sec2
70GB
1,360 sec
4
80GB
1,280 sec
SX-Aurora2048 × 64 × 2048
8
8
100GB
570 sec
16
145GB
415 sec
32
236GB
290 sec
8369GB2,870 sec
4096 × 64 × 4096
16
8
511GB
2,110 sec
32
790GB
960 sec
32
1232GB
2,500 sec
8192 × 64 × 4096
8
642630GB
1,333 sec
4096× 64 × 4096
32
4
667GB
2,075 sec
SX-ACE
8192 × 64 × 4096
32
41311GB
3,900 sec

5.

結論

パイプ流を模擬した低円筒比の環状クエット・ポアズイユ流の亜臨界遷移において，円筒比の異なる流路を対象に DNS を実施した．複数の乱流パフからなる大規模間欠構造のレイノルズ数，
圧力勾配関数依存性を調査し，パイプ流で形成される典型的な乱流パフの性質，有向浸透現象（DP）の特徴と評価することで，模擬パイプ流の乱流遷移に潜む普遍性の有無を探った．
本解析では低円筒比流路（η = 0.1, 0.3）において，パイプ流や環状流で観測されている乱流パフと概ねサイズや形状の一致した局在乱流構造が形成されることを発見した．各流路で形成される
乱流パフはそれぞれ異なる分裂・衝突・減衰挙動を示し，時空間欠性に関する統計量からも差異が抽出された．ただし，いずれも臨界点である層流化の直前では，(1+1)-D DP 普遍性に従うような特徴の一部が得られた．しかし，η = 0.1 の流路における乱流パフは DP のような分裂は起きておらず，得られた DP 普遍性の一部は擬 DP 普遍性の特徴ともいえる．また，下臨界値は両流路でいずれも Rebulk = 2000–2300 前後であり，パイプ流の下臨界値 Reg = 2040 と非常に近い値となった．
本研究結果を含め，壁面せん断流では亜臨界遷移において乱流パフに限らず多様な間欠パターンを形成し，それぞれが特有の分裂挙動を呈すことが分かってきている．現在，主にパイプ流や環状流における局在乱流の挙動について，低次元モデルが提案され[25]，筆者らも DP モデルを基にした確率モデルを検討している．確率モデルの設計には膨大な現象観察が重要で，引続き大規
模高負荷な 3 次元 DNS を長時間実施する必要があるため， SX-Aurora TSUBASA のように高性能
スーパーコンピュータ利用が不可欠である．将来的に，開発中の確率モデルを CFD に組み込むことで，工学的応用に向けた亜臨界遷移下での乱流の活用を模索できると期待する．
本研究結果の詳細については既報[26,27]を参照されたい．

― 32 ―

SENAC Vol. 54, No. （4 2021.10）

謝辞
本研究は，東北大学サイバーサイエンスセンターのスーパーコンピュータ SX-ACE，SX-Auroraを利用することで実現することができた．また，MPI 並列計算にあたっては同センター関係各位に有益なご指導とご協力をいただいた．本研究は，JSPS 科研費（基盤研究 (B) JP19H02071）の助成を受けて行われたものである．また，筆頭著者は，東京理科大学マテリアル人材フェローシッ
プ制度（MEXT 事業）の支援を受けた．ここに記し，感謝の意を表する．
参考文献[1] Wygnanski, J. and Champagne, F. H., “On transition in a pipe. Part 1. The origin of puffs and slugs andthe flow in a turbulent slug,” J. Fluid Mech., 59 (1973), pp. 281–335.
[2] Avila, K., Moxey, D., Lozar, A. de, Avila, M., Barkley, D. and Hof, B., “The onset of turbulence in pipe
flow,” Science, 333 (2011), pp. 192–196.
[3] Mukund, V. and Hof, B., “The critical point of the transition to turbulence in pipe flow,” J. Fluid Mech.,
839 (2018), pp. 76–94, 2018.
[4] 塚原隆裕, 石田貴大, “平面ポアズイユ流の亜臨界遷移における下臨界レイノルズ数”,ながれ:日本流体力学会誌, 34.6 (2015), pp. 383–386.
[5] Tsukahara, T., Seki, Y., Kawamura, H. and Tochio D., “DNS of turbulent channel flow at very low
Reynolds numbers,” In: Proc. 4th Int. Symp. on Turbulence and Shear Flows Phenomena (2005), pp.
935–940.
[6] Andereck, C. D., Liu, S. S. and Swinney, H. L., “Flow regimes in a circular Couette system withindependently rotating cylinders”, J. Fluid Mech., 164, pp. 155–183, 1986.
[7] Cros, A. and Le Gal, P., “Spatiotemporal intermittency in the torsional Couette flow between a rotating
and a stationary disk,” Phys. Fluids, 14, 3755, 2002.
[8] Prigent, A., Grégoire, G., Chaté, H., Dauchot, O. and van Saarloos, W., Large-scale finite-wavelength
modulation within turbulent shear flows, Phy. Rev. Lett., 89 (2002), 014501.
[9] Pomeau, Y., “Front motion, metastability and subcritical bifurcations in hydrodynamics,” Physica D:
Nonlinear Phenomena, 23 (1986), pp. 3–11.
[10] Lemoult, G., Shi, L., Avila, K., Jalikop, S. V., Avila, M. and Hof, B., “Directed percolation phasetransition to sustained turbulence in Couette flow,” Nature Phys., 12 (2016), pp. 254–258.
[11] Sano, M. and Tamai, K., “A universal transition to turbulence in channel flow,” Nature Phys., 12 (2016),
pp. 249–253.
[12] Shimizu, M. and Manneville, P., “Bifurcations to turbulence in transitional channel flow,” Phys. Rev.
Fluids, 4 (2019), 113903.
[13] Chantry, M., Tuckerman, L. S. and Barkley, D., “Universal continuous transition to turbulence in aplanar shear flow,” J. Fluid Mech., 824 (2017), R1.
[14] Kunii, K., Duguet, Y. and Tsukahara T., “Laminar-turbulent coexistence in annular Couette flow,” J.
Fluid Mech., 879 (2019), pp. 579–603.
[15] Hiruta, Y. and Toh, S., “Subcritical laminar-turbulent transition as nonequilibrium phase transition in
two-dimensional Kolmogorov flow,” J. Phys. Soc. Japan, 89 (2020), 044402.
[16] Takeda, K., Duguet, Y. and Tsukahara T., “Intermittency and critical scaling in annular Couette flow,”
Entropy, 22 (2020), 988.
[17] T. Ishida, Y. Duguet, and T. Tsukahara, Transitional structures in annular Poiseuille flow depending on
radius ratio, J. Fluid Mech., 794 (2016), R2.
[18] Ishida, T., Duguet, Y. and Tsukahara T., “Turbulent bifurcations in intermittent shear flows: from puffs
to oblique stripes,” Phys. Rev. Fluids, 2 (2017), 073902.
[19] 塚原隆裕，石田貴大，“直接数値解析による環状流路ポアズイユ流の亜臨界遷移の研究”，東北大学サイバーサイエンスセンター大規模科学計算システム広報 SENAC, 52(2) (2019), pp. 9–17.

乱流パフの時空間欠性と DP 普遍性に関する一考察：
環状クエット・ポアズイユ流の直接数値シミュレーション― 33 ―[20] 竹田一貴，塚原隆裕，“環状クエット流の亜臨界遷移における局在乱流構造と DP 普遍性の円筒比依存性”，日本流体力学会年会 2020 講演論文集，(2020).
[21] 塚原隆裕，竹田一貴，“直接数値解析を用いた仮想拡張した環状流路における亜臨界遷移の研究”，東北大学サイバーサイエンスセンター大規模科学計算システム広報 SENAC，53(4) (2020),
pp. 90–98.
[22] Duguet, Y., Schlatter, P. and Henningson, D. S., “Formation of turbulent patterns near the onset oftransition in plane Couette flow,” J. Fluid Mech., 650 (2010), pp. 119–129.
[23] Gittler, P. “Stability of axial Poiseuille-Couette flow between concentric cylinders,” Acta Mech., 101(1993), pp. 1–13.
[24] 森松浩隆，“低円筒比環状流路の亜臨界遷移域において生じる一次元間欠性に関する研究”，
東京理科大学大学院修士論文，(2021).
[25] Barkley, D., “Theoretical perspective on the route to turbulence in a pipe,” J. Fluid Mech., 803 (2016), P1.
[26] 森松 浩隆，塚原隆裕，“低円筒比環状流路の亜臨界遷移に生じる間欠構造に関する調査”，第34 回数値流体力学シンポジウム 講演論文集，A05-1 (2020).
[27] Morimatsu, H. and Tsukahara T., “Laminar-turbulent intermittency in annular Couette-Poiseuille flow:
Whether a puff splits or not,” Entropy, 22 (2020), 1353.