管内サスペンション流れにおける浮遊粒子の集中現 象

象

著者 山下 博士

発行年 2020‑03‑31

学位授与機関 関西大学

学位授与番号 34416甲第768号

URL http://doi.org/10.32286/00021315

学位授与 2020 年 3 月  関西大学審査学位論文 

管内サスペンション流れにおける浮遊粒子の集中現象 

理工学研究科  総合理工学専攻  基礎・計算物理 

17D6001  山下博士 

理工学研究科  総合理工学専攻  基礎・計算物理  17D6001  山下博士   

管内サスペンション流れにおける浮遊粒子の集中現象   

  粒子サスペンションを一様な断面をもつ直管に流すと，浮遊粒子が慣性の影響により周 囲の流体から揚力（流れに対し垂直方向の力）を受けて，下流の管断面で特定の位置に集中 する現象が起こる．この粒子集中現象は医療や生体流体の分野で，細胞をその媒体から外力 なしに分離するマイクロ流体デバイスを実現させると期待されており，近年活発に研究さ れている．分離デバイス開発に向けてこれまでに多くの研究が行われているが，その設計に 必須である粒子の集中位置予測の理論は未だ不十分であり，現状の流路設計は試行錯誤に よるものが多い．これは粒子集中位置が管断面形状，管断面に対する粒子サイズ（サイズ 比），流量（レイノルズ数）といったさまざまなパラメータに依存して大きく変化すること が原因である．これまでの研究により，円管内流れでは粒子は環状に集まり，矩形管内流れ では管断面内の数点に集まることがわかっている．特に正方形管の場合，流量の増加につれ て粒子集中位置およびその数が変化することが示されている． 

  本研究では，粒子集中位置の複雑な変化を予測するために，パラメータ毎に粒子集中現象 が如何な変化をするのかを系統的に明らかにすることを目的とする．基礎的な管断面形状 である円管と正方形管の場合について，管内サスペンション流れにおける浮遊粒子の集中 現象を数値解析と実験の両面から調べた． 

  揚力が零となる粒子平衡位置を数値的に探索するとともにその安定性を評価することで，

いずれの管断面形状の場合においても，粒子集中位置およびその数の変化が平衡位置の分 岐現象として説明できることが示された．この見通しを基に，サイズ比による粒子集中現象 の変化は統一的に理解されることがわかった．数値解析による予測は，本研究の観察実験に よって確認された． 

  第 1 章「緒言」では，管内サスペンション流れで見られる粒子集中現象についての先行 研究を引用してその複雑さを説明する．円管内流れでは低レイノルズ数の場合に一つ，中レ イノルズ数で二つ，高レイノルズ数で一つの粒子集中位置が管断面内に出現することが報 告されている．正方形管の場合では，いずれのレイノルズ数の場合にも断面軸上に現れる面 心集中点，高レイノルズ数において対角線上の四隅付近に出現する対角集中点，および中レ イノルズ数領域においてそれらの間に現れる中間集中点の存在が示されている．また各々

集中現象における複雑さを明らかにする動機（本研究の目的）を示す． 

  第 2 章「方法」では，数値解析手法と実験方法の詳細を述べる．球に働く揚力の管断面内 分布を求める数値解析，およびさまざまなパラメーターの場合に出現する揚力が零となる 粒子平衡位置を探索してその安定性を判断する手法を説明する．マイクロメーターサイズ の管断面をもつ流路内を浮遊する球形粒子が，管下流断面で通過する位置を観察する方法 を示す．また，用意する作動流体の物性値の制御に関する説明を行う．

  第 3 章「結果と考察」では正方形および円管の場合で得られた結果を説明する．管の断 面形状が正方形の場合に実験的に示された，レイノルズ数による粒子集中パターンの移り 変わりは，断面内に出現する粒子平衡位置に関する分岐現象として捉えられることがわか った．観察される粒子集中点は安定な粒子平衡位置である．特定のレイノルズ数より大きい 場合に中間集中点や対角集中点が観察される現象は，サドル・ノード分岐や亜臨界ピッチフ ォーク分岐としてそれぞれ説明されることがわかった．この見通しを基にサイズ比依存性 を調べると，それぞれの臨界レイノルズ数の変化によって分岐構造自体が変化することが 捉えられ，これまでに報告されていない粒子集中パターンの遷移をとるサイズ比が示唆さ れた．また，先行研究で数値的に予測された比較的サイズ比が大きい場合にのみ低レイノル ズ数領域で出現する対角集中点の存在を，本数値解析でも確認した．さらに，レイノルズ数 の増加に伴ってこの粒子集中パターンから標準的なパターン（面心集中点のみ）へと移り変 わるまでの変化が明らかにされた．本研究で実行された観察実験はそれらの変化を初めて 捉えた． 

  管断面形状が円管の場合で見られた粒子集中位置の数の変化に対しても，粒子平衡位置 の分岐現象として見通せることがわかった．レイノルズ数の増加によって粒子集中位置が 一つから二つとなる変化，二つから一つとなる変化ともにサドル・ノード分岐として捉えら れることが示された．正方形管の場合と同様で，サイズ比による分岐構造の変化が起きるこ とがわかった．粒子集中位置が二つとなるレイノルズ数領域は比較的小さいサイズ比の場 合に現れ，比較的大きいサイズ比となると平衡位置の分岐現象がなくなることが数値的に 示された． 

  第 4 章は「結言」であり，本研究を総括した． 

以上 

1 緒言 1

2 ^方法 5

2.1 数値解析による粒子平衡位置の探索 . . . . 5

2.1.1 球に働く揚力の数値解析 . . . . 5

2.1.2 揚力分布と粒子平衡位置 . . . . 8

2.2 浮遊粒子の観察実験による粒子集中位置測定 . . . . 9

3 結果と考察 12 3.1 正方形管内サスペンション流れで見られる粒子集中現象 . . . . 12

3.1.1 数値解析手法の検証 . . . . 12

3.1.2 100^{オーダーの}Re数で見られる粒子平衡位置の分岐現象 . . . . 15

3.1.3 10^{オーダーの}Re数で見られる粒子平衡位置の分岐現象. . . . 30

3.2 円管内サスペンション流れで見られる粒子集中現象 . . . . 38

3.2.1 数値解析手法の検証 . . . . 38

3.2.2 共存する二つの粒子集中位置 . . . . 40

4 ^結言 44

1 緒言

一般に，流体に密度差のない粒子を混ぜて流すと，粒子は流れに追従して運動する．流れの 可視化に用いられるトレーサー粒子はその一例である．ただし特定の条件下では，流体との相 互作用によって粒子が流れに追従せずに浮遊する場合がある．例えば，粒子の回転に起因する

Magnus効果が挙げられるだろう．一様流れ中に回転する球がある場合，球は主流と球の角速

度に直交する方向の力を受けて流れを横切ろうとする．投球におけるバックスピンやトップス ピンによる球軌道の変化はMagnus効果として知られる．剪断流れ中の粒子に関するSaffman 効果という現象も知られている．剪断流れの方向に粒子が先行または遅れる際，粒子は主流と 垂直の方向に移動しようとする．これらの効果は，粒子の運動が流れ場を大きく変えることに よって現れる．本研究では，管内サスペンション流れにおける浮遊粒子の集中現象として知ら れる，Segr´e–Silberberg効果を考える．中立浮遊粒子が流れに身を任せているにもかかわら ず，流れを横切る方向の力が粒子に作用する．

一様な断面をもつ直管内の流れに流体と同じ密度の球形粒子を混ぜると，粒子の大きさが管 断面幅の1/10程度以上の場合，十分下流の管断面内において粒子が特定の位置だけを通過する．

この粒子配置の集中現象はReynolds^数(Re =U D/ν^{，平均流速}U^，管幅D^{，流体の動粘度}ν) が数百程度までの範囲の円管内流れを用いた実験でSegr´e & Silberberg[1]^{により初めて観察さ} れた．浮遊粒子の通過位置は管断面内の半径の0.6倍程度の動径位置に集中する．以後，発見 者に因みSegr´e–Silberberg効果と呼ばれるようになった．Asmolov[2]は平板間流れに中立浮遊 する球形粒子の運動に対して，粒子サイズが平板間距離に対して十分小さいとした接合漸近展 開法による理論解析を行い，Segr´e & Silberbergが観察した粒子集中現象のシナリオを示した．

遅い粘性流れ，つまりNavier–Stokes方程式にある慣性(非線形)項を無視した流れの場合，管 内流れ中の球形粒子に作用する揚力(主流に対して垂直方向の力)は生じない．一方，慣性(非 線形性)が効く場合，管内流れ中の粒子には流体から揚力が作用し，粒子は主流を横切る運動を 行う．管壁付近を移動する粒子には管中央に向かう揚力が作用し，粒子は壁から離れようとす る．これはwall effectと呼ばれる．一方管軸付近に位置する場合には，粒子は管壁に向かう (速度勾配の絶対値が大きい方向の)揚力を受ける．この揚力はshear-gradient induced lift と呼ばれる．管断面内には，これら二つの効果が均衡する粒子運動の落着位置が現れる．この 位置は粒子に作用する揚力が零になる粒子平衡位置であり，Segr´e & Silberberg^{が観察した粒} 子集中位置に対応する．

円管内サスペンション流れ中の浮遊粒子の集中位置はRe数による複雑な変化を呈する．Matas et al.[3]^は広範なRe^数領域(Re⩽1700)における集中位置を実験によって調べ，Re^数が600^程 度以上ではSegr´e & Silberbergが観察した粒子集中位置(^外側環)に加えて管中央側にもう一つ

別の集中位置(^内側環)が現れることを見出した．彼らが同時に行った理論解析では，外側環が 粒子平衡位置である一方で，内側環の位置で揚力が零にならない．そのため，内側環の観察結 果が粒子の管内運動の過渡的な位置を示したものであるのか，あるいは粒子サイズを無視した 理論解析では捉えられないのかが不明であった．内側環の真仮という問題に対して，Morita et al.[4]とNakayama et al.[5]はさまざまな管長の場合での浮遊粒子の下流断面内通過位置を実験 的に調べ，内側環に対応する位置を通過する粒子が管長が長くなるにつれ減少傾向にあるケー ス(Re^数領域)^{を示した．また}Nakayama et al.[5]は管長を管幅の千倍とした場合でも消えずに 残る二つの粒子集中位置が現れるRe数領域の存在を明らかにし，観察された二つの粒子集中位 置が共に粒子平衡位置であることを粒子サイズを考慮した数値解析[5]^{によって示した，}

Segr´e–Silberberg効果によって見られる粒子集中現象は近年，粒子や細胞を流体から選抜す

るマイクロ流体デバイス開発を可能にする技術として，生体流体の分野で注目されており医療 等への応用が期待されている[6, 7, 8, 9]．マイクロメーターサイズの矩形断面をもつ流路を用

いたDi Carlo et al.[10]の顕微計測実験によって，矩形管内流れに浮遊する粒子の集中位置は，

円管の場合が特定の動径位置であったのに対して，下流断面内の特定の点に集中することが示 された．矩形管の場合に見られる粒子集中点を利用すれば，強い外力による分離方法(^例えば遠 心分離)では損傷してしまうような細胞を，その懸濁液を流路に流すだけで流体から分離できる 可能性がある．

矩形管内流れに浮遊する粒子の慣性集中現象に対して多くの研究がなされてきた．Chun &

Ladd[11]は，正方形管内流れに中立浮遊する複数の球形粒子の運動を数値的に解くことで，粒

子集中点には断面の各辺中央付近と対角線上の四隅付近の二種類が存在することを示した．本研 究では前者を面心集中点，後者を対角集中点と呼ぶ．Re⩽300では浮遊粒子の多くが面心集中 点に落着し，Re⩾500では対角集中点に到達することが数値解析により予測された．Di Carlo

et al.[10]は一辺50µmの正方形断面をもつ流路に粒子または赤血球のサスペンションを流す実

験を行い，低Re数領域(Re⪅50)において両者が面心集中点に到る観察結果を報告した．Choi

et al.[12]はホログラフィーを利用した顕微計測によってガラス管内を浮遊する粒子の断面内通

過位置を捉え，Re数が大きくなる程，より多くの粒子が面心集中点を通過することを示した．

またDi Carlo et al.[13]は正方形管内を浮遊する球形粒子に作用する揚力を求める数値解析を

実行し，実験で観察される粒子集中点はその揚力が零となる粒子平衡位置であることを示した．

Amini et al.[14]の数値解析結果は浮遊粒子が粒子集中点に対して管壁側か中央側に位置するか

によって，管内に生じる小さな二次流れの特徴が変化することを示した．Hood et al.[15]は摂 動展開した上で正方形管の場合の揚力を数値的に求め，粒子直径と管幅のサイズ比(β =d/D， 粒子径d)によって粒子平衡位置の動径距離が変化することを説明した．

Miura et al.[16]は粒子間相互作用を無視するために希薄な粒子サスペンション(^{粒子の体積} 分率は10⁻²%のオーダー)を用い，ミリメーターサイズの正方形断面をもつ管内流れに浮遊す る粒子を管出口正面から観察し，下流断面における粒子の通過位置を解析した．幅広いRe数領 域(Re⩽1200)を調べ，高Re数領域(Re⪆500)において面心集中点と対角集中点が共存する ことを実験的に示した．Nakagawa et al.[17]はMiura et al.[16]が行った希薄サスペンション流 れを模擬するような，正方形管内流れに中立浮遊する単一の球形粒子の運動を数値的に解いた．

高Re数領域での粒子軌道の落着位置は面心集中点と対角集中点であり，それらが共存すること を示した．これらの観察結果および数値解析結果はChun & Ladd[11]の数値解析による予測と 異なる結果であり，粒子間相互作用が粒子集中位置に影響を及ぼすことを示唆している．Shichi

et al.[18]はマイクロメータースケールの管断面をもつガラス管を用いて無次元管長を大きくと

り，管内に流入した全ての浮遊粒子がいずれかの粒子平衡位置に到達した状態の観察を試みた．

サイズ比β≈1/10において，正方形管内流れに浮遊する球形粒子の集中位置には，これまでに 確認されていた面心集中点と対角集中点の他に，それらの間に位置する粒子集中点が存在する ことが報告された．これを本研究では中間集中点と呼ぶ．この新たな集中点はRe数が大きい 場合になるにつれて対角集中点に近付く特徴をもつことが実験的に示され，対応する数値解析 [19]によっても確認された．さらにそれぞれの粒子集中点が現れるRe数領域が異なるために，

低Re数領域では面心集中点のみ，中Re数領域では面心集中点と中間集中点，高Re^数領域で は面心集中点と対角集中点が現れるといった粒子集中パターンのRe数による移り変わりが起こ ることが示された．

矩形管の場合で見られる粒子集中現象はRe数だけでなく，サイズ比β^{や管断面のアスペク} ト比によっても大きく変化することが報告されている．サイズ比β依存性に関して，正方形管 内流れを浮遊するサイズ比β ≈1/3以上の比較的大きな球形粒子の集中位置が，低Re^数領域 (⪅40)において対角集中点だけとなることを予測した数値解析結果が報告されている[20]^．こ の粒子集中パターンを実験で確認した研究はない．管断面のアスペクト比に関して，つまり長 方形断面をもつ流路の場合の粒子集中位置は，低Re数領域で断面の各長辺の中央付近[21, 22]， Re数がより大きい場合となると短辺中央付近が加わること[23, 24]が示されている．

このようにSegr´e–Silberberg効果が起こす管内サスペンション流れにおける浮遊粒子の集中 現象には，管断面形状，粒子径と管断面のサイズ比β^，Re数といった多くのパラメーターに依 存する多様な粒子集中位置や集中パターンが存在することがわかってきている．しかし，それら のパラメーターによる粒子集中位置の変化を予測する理論は未だ不十分であり，そのため，懸濁 粒子を分離するマイクロ流体デバイスの最適な設計は現状困難となっている．本研究では，粒 子集中位置の予測のために，パラメーター毎に粒子集中現象が如何な変化をするのかを系統的 に明らかにすることを目的とする．管内流れに中立浮遊する球形粒子の運動を数値解析と実験

の両面から調べる．粒子集中現象を大きく変えるパラメーターの一つである管断面形状に関し て，基礎である円形断面と正方形断面の場合に注目する．Segr´e–Silberberg効果によって現れ る粒子集中パターンの変化を粒子平衡位置の分岐現象として説明し，サイズ比βやRe数によ る変化の見通しを提示する．

2 方法

2.1 数値解析による粒子平衡位置の探索

管内流れに浮遊する球に作用する慣性に起因する揚力の管断面内分布を数値解析により求め，

揚力が零となる断面内位置(粒子平衡位置)やそのパラメーター変化を詳細に調べる．本研究で は管の断面形状として，正方形と円形の二通りを選んだ．ここでは正方形管の場合を主軸とし て，定式化やデータ解析手法を説明する．

2.1.1 球に働く揚力の数値解析

正方形断面をもつ直管内流れに中立浮遊する単一の球の運動を考える．管内を満たす流体は

非圧縮Newton流体とし，連続の式とNavier–Stokes方程式に支配される．流れに浮遊する球

はNewtonの第二法則に従う．ただし球に作用する揚力を求めるために，球の並進運動につい

ては管軸方向のみとして，管断面方向の運動は制限する．このとき球はどの方向にも自由に回 転運動を行う．Fig. 1に描いたように，座標系は球の運動とともに管軸(x^軸)^{方向に並進移動} する．管断面中央を原点とし，断面各辺と平行にy^およびz^{軸をとる．}

O x

y z

ω_p wall velocityvpx

D

periodic repeat length32d

Figure 1. 正方形管内流れに浮遊する球．座標系は球の並進運動とともに移動する．

正方形管流れに中立浮遊する球に働く揚力を数値的に解く手法として，Kajishima et al.[25]

の埋め込み境界法を選んだ．流体と球の運動をつないだ支配方程式は

∇ ·u= 0， (1)

ρ [∂u

∂t + (u· ∇)u ]

=−∇p+µ∇²u+f_p−ρd

dt(vpxex)^， (2) m_pdv_px

dt =−

∫

V

f_pxdV， (3)

Ip

dω_p dt =−

∫

V

rp×f_pdV (4)

となる．V は計算領域である．ここでのuは液相と固相を統一的に表現した速度で

u= (1−α_p)u_f+α_pu_p (5)

とする．u_f^{は流体の速度である．}upは球表面および内部の速度であり

u_p=ω_p×r_p (6)

として，角速度ωpと球中心からの位置ベクトルrpで記述される．vpxexは球の並進速度であ り，そのために管壁の速度は−vpxexとなる．式(5)^に現れるαpは，計算セル内の液相と固相 を判断するためのスカラー関数であり，液相ではα_p = 0，固相ではα_p= 1となるように

αp = 1 2

(

1−tanh|r_p| −0.5d δα

)

(7)

と定義する．δ_αは液相と固相の界面に関するパラメーターであり，0< α_p<1の幅を決定する．

本研究ではNakagawa et al.[17]を引用しδ_α= 0.6δ_x = 0.03dとした．δ_xは計算領域の差分格子 幅である．式(2–4)でのf_pは，埋め込み境界法における液相と固相の相互作用であり，

f_p= ρ δt

αp(up−u)ˆ (8)

と表される．δ_tは式(2)左辺第一項を前進差分で離散化したときの時間刻み幅である．uˆ は

uˆ =uⁿ+δt

[

−1

ρ∇p−(uⁿ· ∇)uⁿ+µ

ρ∇²uⁿ− d

dt(vpxex) ]

(9)

であり，計算領域すべてが液相として得られる仮速度である．上付き添字nは離散的な時刻に 対応する．式(8)は式(2)の離散式

uⁿ⁺¹ = ˆu+δt

f_p

ρ (10)

が液相(α_p = 0)でuⁿ⁺¹ = ˆu，固相(α_p= 1)でuⁿ⁺¹ =u_pとなるように与えられている．

境界条件は管壁での滑りなし条件と管上流断面と下流断面に対する周期境界条件である．そ れらは

u|_y=±^D

2 =−v_pxe_x，u|_z=±^D

2 =−v_pxe_x， (11)

u|_x=−16d= u|_x=16d^， ∂u

∂x

x=−16d

= ∂u

∂x

x=16d

(12)

とそれぞれ表される．またこの定式化では液相と固相の界面(球表面)における滑りなし条件が

満たされる．流れは管の上流断面と下流断面の間の圧力差δpによって駆動し，

p|_x=−16d= p|_x=16d+δ_p (13)

と表される．

空間に関して二次中心差分法で離散化し，定常状態，つまり球に作用する抗力とトルクが零 であるときの揚力

F =−

∫

V

f_pdV =Fyey+Fzez (14)

を求める．

正方形管の場合にはFig. 2aのように球のみを埋め込み境界法で表現したが，円管の場合では

Fig. 2bのように球だけでなく円管も埋め込み境界法で表現する．体積重み付き平均速度(5)と

Navier–Stokes方程式(2)をそれぞれ

u= (1−αp−αw)uf+αpup+αw(−vpxex)^， (15) ρ

[∂u

∂t + (u· ∇)u ]

=−∇p+µ∇²u+f_p+f_w−ρd

dt(v_pxe_x) (16) と変更する．αwは円管を表現するためのスカラー関数であり，f_wはそれに対応する液相と固 相間の相互作用を表す体積力である．その他は前述の正方形管の場合の同様である．

(a) (b)

x y z

x y z

Figure 2. 埋め込み境界法．黒色の枠が計算領域である．(a)正方形管内流れに浮遊する球．等

値面はα_p = 0.5(球表面)を表す．(b)円管内流れに浮遊する球．等値面はα_p = 0.5(球表面) お よびα_w = 0.5(円管壁)を表す．青線は静止座標系での流線を表す．(a)β= 1/5，Re = 100．(b) β= 1/10，Re = 100．

2.1.2 揚力分布と粒子平衡位置

2.1.1項の方法で球が管断面内のさまざまな位置にある場合に受ける揚力を計算し，揚力の断

面内分布を求める．本研究では粒子平衡位置の見積りの精度とそのための計算コストのバラン スと取るために，計算点YijとしてGauss–Lobatto^点ξnを用いた

Y_ij =Y ξ_ie_y+Y ξ_je_z，ξ_n= cosπn

N ^，n= 0,1, . . . , N (17) をとる．Y ^は(D−d)/2より小さい定数とする．また対称性から0⩽i⩽j⩽N/2^{の場合のみ} を解く．具体的には一つのパラメーターセット(β,Re)^{に対して，}Fig. 3a^のようにN = 14^とし て計36ヶ所での揚力を計算する．選点数はYamashita et al.[26]を引用して決定された．この 選点上で計算された揚力F_ij =F(Y_ij)を

F˜(y, z) =

∑N i=0

∑N j=0

F_ijh_i(y/Y)h_j(z/Y) (18)

とLagrange補間する．補間のために用いる多項式h_k(ξ)は

h_k(ξ) = 2 N

∑N

l=0

1 ckcl

T_l(ξ_k)T_l(ξ)，k, l= 0,1, . . . , N (19)

としてGauss–Lobatto^点ξnとChebyshev^多項式Tnを用いて定義される．多項式(19)^中のcn

は，0< n < N^のときcn= 1^，n= 0, N ^のときcn= 2^{とする．以降，式}(18)^のF˜ ^{を簡単のた} めにF ^{と記述する．}

管断面形状が矩形である場合，その断面に関する周方向対称性の欠如のために，揚力に周方 向成分が生じる得る．この考えから本研究では，内挿によって得られた揚力分布F ^{の動径方向} および周方向成分が零である等値線(Fig. 3bにおけるそれぞれ青線と緑線)^{に注目する．}Fr = 0 の等値線の内側(^管中央側)^ではFr >0^{であり，外側}(^管壁側)^ではFr <0^{である．そのため，}

この曲線(^青線)は円管内流れの場合での環状の粒子平衡位置に相当すると考えられる．Fr = 0 の等値線上でF_θ= 0の位置，つまりF_r = 0の等値線(青線)とF_θ = 0の等値線(緑線)の交点 が粒子平衡位置r_eである．r_eの詳細な座標はF(r_e) =0をNewton–Raphson法で解くことで 求める．また管内流れに浮遊する球は計算された揚力に沿って管断面内運動を行うと仮定して，

Fig. 3bの赤の矢印付き曲線のような球中心の断面内軌道を計算することで，それぞれの粒子平

衡位置の安定性を見積る．球の軌道の落着点となる粒子平衡位置は安定である(Fig. 3bの○印)． 平衡位置近傍の球軌道がそれに近づくものと離れるものがある場合，その粒子平衡位置はサド ル点である(Fig. 3b^{の断面対角線上の×印})．粒子平衡位置周りのどの軌道もそれから離れる場 合，それは不安定な平衡位置である(Fig. 3b^{の断面中央の□印})^．

(a)

y z

0 D/4 D/2

D/2

D/4

0

0.04

|F|/(ρU2d2)

0

Y

(b)

y z

0 D/4 D/2

Figure 3. (a) 球に作用する揚力の管断面内分布Fij(^矢印)．灰色の領域は球中心が存在できる 断面内領域を表す．(b)揚力分布を内挿することで見積もられた正方形断面内の粒子平衡位置．

○印は安定な平衡位置，×印はサドル点，□印は不安定な平衡位置である．青の曲線は揚力分 布の動径方向成分F_rが零となる等値線を表す．軸上，対角線上，管壁の沿うような曲線の緑線 は揚力分布の周方向成分F_θが零である等値線を表す．赤の矢印付き曲線は，計算された揚力に 沿って球が側方移動するとして見積もられた，球中心の断面内軌道を表す．球は矢印の方向に 運動する．

円管の場合では式(17, 18)を

Ri=Y ξi，n= 0,1, . . . , N/2 (20)

F˜_r(r^′) =

∑N/2

i=0

F_r(R_i)h_i(r^′/Y) +

∑N

i=N/2+1

F_r(R_N−i)h_i(r^′/Y)，−Y ⩽r^′ ⩽Y (21) として離散的な計算結果から連続的な揚力の動径方向分布Fr(r) = ˜Fr(r)^{を用意し，粒子平衡} 位置や揚力分布のパラメーターによる変化を評価する．

2.2 浮遊粒子の観察実験による粒子集中位置測定

Shichi et al.[18]が開発した正方形管内流れに中立浮遊する球形粒子の正面観察実験を行い，

管下流において粒子が通過することができる断面内位置を解析することで，数値解析によって予 測された粒子集中現象との比較を行う．Fig. 4aは実験装置の概略図である．粒子サスペンショ ンは，密度1.05×10³kg/m³のポリスチレン製球形粒子(Thermo Science)^{とそれと同密度のグ} リセリン水溶液を混合することで用意される．このとき粒子間相互作用を無視するために，サ スペンション中の粒子の体積分率を0.02%程度と希薄にした．用意された粒子サスペンション

をシリンジを介して，一辺が数百µmの正方形断面をもつガラス管(VitroCom)^{内にシリンジポ} ンプによる陽圧で流す．Fig. 4aに描いたように，長作動距離対物レンズを取り付けた高速度カ

メラ(Photron FASTCAM mini AX100)をガラス管出口の断面に対して正面に設置する．浮遊

粒子が管出口を通過する際に撮影されたFig. 4bのような画像からImageJ (National Institutes

of Health)を用いて，粒子中心の管断面内位置の解析を行う．

(a)

(b)

6–600mm

Flow direction

Square glass duct Tub

Adapter

Backlight Sylinge pump

Dilute suspension

(∼0.02% vol. frac.) Objective lens

High-speed camera Focal plane

300µm

100µm

Figure 4. (a) 正方形管内流れに中立浮遊する球形粒子の正面観察実験．(b) 管出口を通り過ぎ

る粒子の撮影画像．

各々の実験のReynolds数を精度良く決定するために，分散媒であるグリセリン水溶液の密度 と粘度はその溶液温度を考慮する経験式[27, 28]によって見積られる．粒子の中立浮遊条件を 満たす密度のグリセリン水溶液を用意するために，グリセリンの重量比W =Mg/(Mg+Mw) を求める．MgとMwはそれぞれグリセリンと純水の重量である．グリセリン水溶液の密度ρgw

をグリセリンの密度ρ_gと純水の密度ρ_wで表すと

ρ_gw= κρ_gρ_w (1−W)ρg+W ρw

(22)

となる．グリセリンと純水を混合することによる体積収縮率κ^は

κ= V_g+V_w Vgw

= 1 +c[

sin(W^1.31π)]0.81

，c= 1.78×10⁻⁶T²−1.82×10⁻⁴T+ 1.41×10⁻² (23) というVolk and K¨ahler[27]^{の経験式で溶液温度}T ^と重量比W ^{から得られる．}Vgはグリセリ ンの体積，V_wは純水の体積，V_gwはグリセリン水溶液の体積である．Volk and K¨ahler[27]は グリセリンと純水の密度の温度に関する経験式も提案しており，本研究ではそれらを採用した．

ρ_gw= 1.05×10³kg/m³ となる重量比W は，サーミスタ温度計(AS ONE WT-100)により測定 された溶液温度の実測値Tを用いて式(22, 23)を反復計算することで得られる．経験式から得ら れたW ^{を用いて作成された，}19.1–23.2^◦Cでのグリセリン水溶液の密度測定(Kyoto Electronics DA-650)^結果は1.0497×10³±4.3955×10⁻¹kg/m^{3 であった．}

グリセリン水溶液の粘度µgwはCheng[28]^の経験式

µ_gw=µ¹_g^−Aµ^A_w，A= 1−W+ abW(1−W)

aW +b(1−W)^，a= 0.705−0.0017T，b= (4.9 + 0.036T)a^2.5 (24) によって，実測した溶液温度T ^とVolk and K¨ahler[27]の経験式を用いて得られた重量比W から見積もられる．グリセリンと純水の粘度，µgとµwについても，本研究ではCheng[28]^に よって提案された温度に関する経験式を採用した．この経験式の検証としてUbbelohde^粘度

計(Sibata)による測定結果との比較を行い，差が0.5%より小さいことを確認した．式(24)^か

ら得られた19.1–23.2^◦Cでのグリセリン水溶液の粘度は，主にその温度に依存して，µ_gw = 1.67–1.87×10⁻³Pa·s となった．

3 結果と考察

3.1 正方形管内サスペンション流れで見られる粒子集中現象

3.1.1 数値解析手法の検証

2.1節で説明した本研究の数値解析手法およびデータ解析の検証のために，正方形管内流れに浮 遊する球に作用する揚力の分布について先行研究との比較を行った．Fig. 5はサイズ比β = 1/4.5

の場合でRe = 40のときの，それぞれの研究で得られた管断面y軸上の揚力分布である．本研

究で計算された揚力分布(Fig. 5の◆印)はDi Carlo et al.[13]の有限要素法による結果(■印) やNakagawa et al.[17]の差分法での埋め込み境界法で得られた結果(○印)と一致することが確 認された．また内挿によって連続関数として得られた揚力分布(赤の曲線)は先行研究より少な い計算点であるにもかかわらず，その分布の特徴をよく捉えている．またFig. 5^{から，本研究で} 得られた粒子平衡位置がNakagawa et al.[17]の結果と非常によく一致していることが見てとれ る．本研究で用いた数値解析手法(差分法での埋め込み境界法)は彼らと同じものであるが，採 用した座標系が異なる．本来，結果は座標系に依存しないが，数値解析であるために差が生じ る．Nakagawa et al.[17]やYamashita et al.[26]の研究でとられた静止座標系では，埋め込み境 界法で表現された球が差分格子を横切る際に計算される力が数値振動する．そのため，定常に なったと考えられる数百から数千ステップの時間データを平均することで球が受ける揚力は評 価された．一方，本研究ではガリレイ変換を施すことで球に固定した座標系を採用したために，

球が格子を横切ることがなく，揚力の時間発展における数値振動はほとんど起こらない．した がって揚力の評価に関して，最終状態の時刻の計算値をそのまま採用すれば十分である．これ ら二種類の座標系で得られる数値結果がよく一致することは，本研究の数値解析が妥当である ことを示唆する．

本数値解析の妥当性を確認するために，内挿で得た粒子集中点，および初期位置からそこに到 るまでの側方軌道をShichi et al.[18]の実験結果と重ね合わせた．Fig. 6^{の黒点で描かれた分布} 図は，サイズ比β = 1/8^でありRe = 100の場合に，それぞれ異なる管長L^{で観察された浮遊粒} 子の下流断面内における通過位置の分布である．Fig. 6aはこの中で最も管入口近く(^{無次元管長}

L/D= 50)で観察された分布図であり，浮遊粒子が少し側方運動を行った場合である，粒子集

中現象の過渡を表す分布図である．Fig. 6b,cから，さらに管長が長い場合には，浮遊粒子が通 過した分布域が狭まっていくことが見てとれる．Fig. 6dは粒子が管内を長距離(L/D= 750)浮 遊したの場合の観察結果である．Shichi et al.[18]は，L/D= 1500とした場合でもL/D= 750 の場合と同様の粒子分布図を得た．Fig. 6dの粒子分布は粒子の側方運動が完了した状態の観察 結果と判断できる．Fig. 6a–dの粒子分布図の推移は，Fig. 7に描いた本数値解析結果を用いて，

次のように説明される．正方形管内流れに浮遊する球形粒子の側方運動には二つの段階が存在す

-1.2e-01 -1.0e-01 -8.0e-02 -6.0e-02 -4.0e-02 -2.0e-02 0.0e+00 2.0e-02 4.0e-02

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

y/(D/2) Fy/(ρU2d2)

Di Carlo et al. 2009 Nakagawa et al. 2015 Present study

Figure 5. 管断面の軸上に球が位置する場合に受ける揚力の分布．サイズ比β = 1/4.5，Re = 40．

■印: Di Carlo et al.[13]^，○印: Nakagawa et al.[17]^，◆印: 本研究で得られた数値結果．◆印 をつなぐ曲線は2.1.2項で示した内挿方法により見積もられた揚力分布である．

る．第一段階(Fig. 7の緑の軌道)において，管に流入した浮遊粒子は，その初期位置が管壁付近 であればそこから離れ，管中央付近であれば管壁側に向かうといった側方移動を行い，pseudo Segr´e–Silberberg (pSS) ring[18]と名付けられた環状の位置(赤の破線)上に到る．この特徴は，

管入口での一様な粒子分布からFig. 6a^を経てFig. 6bとなる変化と一致する．第二段階(Fig. 7 の青の軌道)^{では，浮遊粒子は}pSS ringに沿った断面内運動を行い，最終的に粒子集中点に到達 する．これはFig. 6c^を経てFig. 6dの粒子分布となることに対応する特徴である．本数値解析 により予測された粒子の側方軌道(Fig. 7^またはFig. 6^{の赤の矢印付き曲線})^はShichi et al.[18]

の示した粒子分布の変化を定性的に説明するだけでなく，それらの位置に関しても定量的によ く一致する．Fig. 6a–cにおいて，粒子分布の輪郭と数値解析で得られたpSS ring (破線)の位置 が重なり，またFig. 6dから両者で得られる粒子集中位置が一致している．

(a) L/D= 50 (b) L/D= 125 (c) L/D= 250 (d) L/D= 750

D/2

0

−D/2

Figure 6. Shichi et al.[18]の実験により観察された粒子通過位置の管長に伴う推移と本研究で 得られた球の側方軌道との比較．サイズ比β = 1/8^，Re = 100．それぞれの黒点は通過した粒 子の中心位置を表す．いずれの粒子分布も250点以上の粒子数を含んでいる．L/D^{は無次元の} 管長である．赤の矢印付き曲線は本研究の数値解析によって予測された浮遊粒子の側方軌道を 表し，その中で破線で描いた粒子軌道はpSS ring[18]に対応する．○印は数値解析から予測さ れた粒子集中点を表す．

pS Srin

g 1st stage

1ststage 2nd stage

2nd stage

Figure 7. 正方形管内流れに浮遊する粒子の運動における二段階の側方移動．第一段階(緑の軌

道)^{において，浮遊粒子は}pSS ring (^赤の破線)^{に向かう．第二段階}(^青の軌道)^{では，粒子は}pSS ring^{に沿って粒子集中点}(^○印)^{に到達する．}

3.1.2 100^{オーダーの}Re数で見られる粒子平衡位置の分岐現象

正方形管で見られるSegr´e–Silberberg効果による浮遊粒子の下流断面内集中パターンについ て，Shichi et al.[18]^{はサイズ比}β ≈ 1/10の場合を実験的に調べ，Fig. 8a–d^{に描いたような，}

Re数により変化する四種類の粒子集中配置を発見した．Re⪅200で観察される粒子集中位置 はFig. 8aのように面心集中点のみであり，200⪅Re⪅300では面心集中点と中間集中点(Fig.

8b)が観察される．より高いRe数の場合には，限られたRe数領域で面心集中点，中間集中点，

および対角集中点の三種類が共存する粒子集中現象(Fig. 8c)が観察される．さらに高いRe数 領域では中間集中点は消滅し，面心集中点と対角集中点が管断面内に現れる(Fig. 8d)．実験的 に示されたこのような粒子集中現象のRe数変化の背後にひそむメカニズムを明らかにするため に，粒子と管断面のサイズ比1/10⩽β ⩽1/3^{の範囲の場合に，}100⩽Re ⩽400^{で現れる粒子} 平衡位置を数値解析によって調べた．

面心集中点 中間集中点 対角集中点

(a) (b) (c) (d)

Figure 8. Shichi et al.[18]が報告した，正方形管内流れに中立浮遊する球形粒子の下流断面に おける粒子集中現象のRe数変化．(a–d)は観察される粒子集中パターンの模式図であり，右ほ ど，Re数が大きい場合である．(a)面心集中点が現れる粒子集中パターン，(b) 面心集中点と 中間集中点のパターン，(c)面心集中点，中間集中点，および対角集中点が共存するパターン，

(d)面心集中点と対角集中点のパターン．

2.1.2項で示した粒子平衡位置の探索とその安定性の評価をパラメーター毎に行うことで，管

断面内には複数種類の粒子平衡位置が現れること，それらは不安定，サドル，安定な平衡位置 と区別できること，粒子平衡位置の安定性はRe数によって変化する場合があることを捉えた．

100⩽Re⩽400のRe数領域は，平衡位置の変化に関する複数の臨界Re数の存在のために，安 定な平衡位置の組が異なるRe数領域毎に分類できることがわかった．数値解析から予測され た安定な粒子平衡位置は実験で観察される粒子集中点に，臨界Re数は粒子集中点のRe数によ る出現や消滅に関連する．つまり実験で見られる粒子集中パターンの変化(Fig. 8a–d)^は粒子平 衡位置の分岐として説明できることがわかった．また粒子直径と管断面のサイズ比β^{によって，}

粒子平衡位置毎の臨界Re数はそれぞれ複雑な変化をし，そのために異なる粒子集中パターン遷

移(^例えばFig. 8a,d,c,d^{という遷移})が存在し得ることが明らかとなった．以降，その詳細を数

値解析結果の代表例を用いて説明する．

i.^サイズ比β = 1/8^の場合

粒子直径と管幅の比β = 1/8の場合に100⩽Re⩽400の範囲で管断面内に出現する粒子平 衡位置とその安定性をFig. 9^{に描いた．横軸は}Re数，縦軸は粒子平衡位置の方位角θeである．

管断面中央にはどのRe数の場合でも不安定である粒子平衡位置が存在するが，以降の議論に関 与しない．そこで，本研究では縦軸に方位角をとったり，正方形断面の対称性から管断面の第 一象限(0⩽θe⩽π/2)に現れる粒子平衡位置を描いている．赤矢印はpSS ring^{に対応する粒子}

軌道(Fig. 7赤破線を参照)の進行方向を表す．粒子の側方運動の特徴において二段階目に相当

するこの軌道は粒子平衡位置の安定性の判断に直結する．Fig. 9にある×印を結ぶ破線はサドル 点である粒子平衡位置を表す．一方，○印をつなぐ実線は安定な平衡位置を表す．Fig. 9から，

調べられたどのRe数においても管断面軸上(θ_e= 0, π/2)および対角線上(θ_e=π/4)に粒子平 衡位置が存在することが見てとれる．また特定のRe数を超えると(Re>Re_I1)，軸と対角線の 間(0< θe< π/4^またはπ/4< θe< π/2)^{に複数の粒子平衡位置}(Fig. 9^のA^，B^，およびC)^が 出現する．軸上平衡位置は安定，平衡位置A^{はサドル点，平衡位置}B^{は安定，平衡位置}C^は サドル点である粒子平衡位置であることが見てとれる．この内，平衡位置B^とC^はRe>ReI2

となると消滅する．一方，対角線上平衡位置は特定のRe^数(Re_C1)を境にサドル点から安定な 平衡位置に変化する．ここでRe_I1,Re_I2,Re_C1は粒子平衡位置の分岐が起きる臨界Re数である．

Re<Re_I1，Re_I1<Re<Re_C1，Re_C1<Re<Re_I2，およびRe_I2<ReのRe数領域ではそれ ぞれ異なる安定な粒子平衡位置の組(粒子集中パターン)となる．

Re<ReI1において管断面内に出現する粒子平衡位置は，断面中央を除くと，軸上(θe= 0, π/2) と対角線上(θe =π/4)^{の二種類である．この}Re^{数領域の代表例として}Re = 200^{の場合の数} 値解析結果をFig. 10aに描いた．断面対角線上の四隅付近に粒子平衡位置(Fig. 10aの×印)が 存在することがわかる．赤破線で描いたpSS ringに対応する粒子軌道がその位置から離れる軌 道であるため，この平衡位置はサドル点である．一方，軸上の粒子平衡位置(Fig. 10aの○印) は側方軌道の落着点となっていることが見てとれ，この位置が安定な平衡位置であることがわ かる．Re<ReI1では軸上の平衡位置のみが安定であるため，実験ではFig. 8a^{のような粒子集} 中パターンが観察されることになる．赤の一点鎖線で描いた粒子軌道はセパラトリクスであり，

この軌道に挟まれる領域に初期断面内位置をもつ浮遊粒子は対応する安定な粒子平衡位置に落 着する．例えばFig. 10a^では，θ =−π/4, π/4上にセパラトリクスが存在することが見てとれ る．浮遊粒子の初期位置がこれらに挟まれた領域内(Fig. 10a^{の橙色領域})^{であれば，}θ = 0^上 の平衡位置に集中する．Fig. 11は数値解析から得られた安定な粒子平衡位置とShichi et al.[18]

および改良された観察実験で得られた粒子集中点を比較した図である．Fig. 11aのRe = 100の

0 0.1 0.2 0.3 0.4 0.5

100 150 200 250 300 350 400

Re θe/π

ReI1 ReC1 ReI2

A B

C

Figure 9. サイズ比β= 1/8の場合の粒子平衡位置の分岐図．横軸にRe数，縦軸に平衡位置の

方位角θ_eがとられる．○印をつなぐ実線は安定な粒子平衡位置を表す．×印をつなぐ破線はサ ドル点である．赤の矢印はpSS ring上での側方軌道の方向を模式的に表す．

場合の比較から，数値解析の予測と観察実験結果がよく一致することがわかる．

ReI1 <Re <ReC1では，前述したRe <ReI1で見られた軸上と対角線上の粒子平衡位置に 加えて，それらに挟まれた領域(0< θ < π/4^等)^{に二種類の平衡位置}(Fig. 9^のA^とB)^が現れ ることがわかった．このRe^{数領域の代表例として}Re = 270の場合に数値解析から得られた浮 遊粒子の側方軌道と粒子平衡位置をFig. 10bに描いた．それぞれの平衡位置の周りの粒子軌道 から，軸上の平衡位置とBは安定な平衡位置であり，平衡位置Aと対角線上の平衡位置にサド ル点が捉えられる．セパラトリクスで挟まれた橙色領域に初期断面内位置をもつ浮遊粒子は軸 上平衡位置に落着し，緑色領域内であれば平衡位置Bに落着する．安定である軸上平衡位置は 実験で観察される面心集中点に，Bは中間集中点に対応し，Re_I1<Re<Re_C1ではFig. 8bの 粒子集中パターンとなることを本数値解析は示唆する．Fig. 11bにRe = 280の場合の数値解析 結果と実験結果の比較を描いた．軸上平衡位置およびBの位置は，浮遊粒子が通過する観察結 果と一致する．Re = 200 (Fig. 10a)の場合と比較すると，軸上および対角線上の平衡位置の安 定性は変化していないことが見てとれる．Re = 270^の場合のpSS ring^{を成す粒子軌道}(A^か ら軸上平衡位置への軌道，A^からBへの軌道，および対角線上平衡位置からB^への軌道)^を，

Re = 200の場合の対角線上平衡位置から軸上平衡位置への軌道と比較すると，AおよびBの

粒子平衡位置がpSS ring上に現れることがわかる．

第三のRe^{数領域である}ReC1 <Re<ReI2では，第二のRe^数領域(ReI1 <Re <ReC1)^で 存在した軸上平衡位置，対角線上平衡位置，平衡位置A^とBに加えて，さらにもう一つ別の平

(a) ^{Re = 200} (b) ^{Re = 270}

(c) ^{Re = 300} (d) ^{Re = 350}

y D/2 z

D/4

0

0 D/4 D/2

A B

A B C

A

Figure 10. 管断面内に現れる粒子平衡位置．β = 1/8^．(a) Re = 200^，(b) Re = 270^，(c) Re = 300^，(d) Re = 350^．□印: 不安定な平衡位置，×印: ^{サドル点，○印}: ^{安定な平衡位置} (粒子集中点)．赤破線はpSS ringに対応する粒子の側方軌道である．赤の一点鎖線はサドル点 である粒子平衡位置に到達する粒子軌跡である．赤矢印はその位置での揚力の方向を示す．青 線はF_r = 0の等値線であり，緑線はF_θ = 0の等値線である．

衡位置Cが対角線上平衡位置とBの間に出現する．Fig. 10cに描いたRe = 300の場合の結果 を見ると，軸上平衡位置，平衡位置A，およびBの安定性は第二のRe数領域の場合と同じで あることがわかる．一方，対角線上平衡位置は安定な平衡位置に変化していることが見てとれ る．新たに出現した平衡位置Cはサドル点である．この第三のRe数領域において安定な粒子 平衡位置は，軸上平衡位置，平衡位置B，および対角線上平衡位置の三種類である．それぞれ が実験での面心集中点，中間集中点，および対角集中点に対応する．三種類の集中点が共存す る粒子集中パターン(Fig. 8c)が起こり得ることを数値解析によって確認した．Re = 325^の場合 に数値解析で得られた安定な粒子平衡位置と実験結果の比較をFig. 11c^{に描いた．実験で観察} された面心集中点と中間集中点の位置は，それぞれ軸上平衡位置および平衡位置B^{とよく一致}

(a) ^{Re = 100} (b) ^{Re = 280} (c) ^{Re = 325}

(d) ^{Re = 350} (e) ^{Re = 400}

Figure 11. 数値解析により予測された安定な粒子平衡位置(^○印)^とShichi et al.[18]^および本 研究の観察実験で得られた粒子集中点(粒子分布)の比較．それぞれの赤点は通過した浮遊粒子 の中心位置である．β= 1/8．(a) Re = 100，(b) Re = 280，(c) Re = 325，(d) Re = 350，(e)

Re = 400．(c–e)に描いた赤破線は観察実験で用いた管幅D= 400µmの正方形管の断面形状を

表す．断面四隅の丸みの半径はおよそ0.1D^である．

している．一方，対角線上平衡位置を通過する浮遊粒子は観察されず，数値結果と実験結果に 差異があることがわかった．この差に関しては後述する．

第四のRe^{数領域である}ReI2<Reにおいて得られた数値解析結果の代表例として，Re = 350

の場合をFig. 10dに描いた．安定である軸上平衡位置，サドル点である平衡位置A^{，および安}

定である対角線上平衡位置の三種類が，第三Re数領域の場合と同様に，管断面内に存在する ことが見てとれる．一方，第三のRe数領域では存在した安定な平衡位置Bとサドル点である Cは消滅していることがわかる，安定な粒子平衡位置は軸上平衡位置と対角線上平衡位置であ り，そのためこのRe数領域ではFig. 8dのような粒子集中パターンが実験で観察されることに

なる．Fig. 11dおよびeに，第四Re数領域で観察された実験結果と，対応する数値解析結果を

重ねて描いた．Re = 400^の場合(Fig. 11e)は，観察された面心集中点および対角集中点と軸上 平衡位置および対角線上平衡位置がそれぞれよく一致する．Re = 350^の場合(Fig. 11d)^は，面 心集中点と軸上平衡位置，対角集中点と対角線上平衡位置はそれぞれ一致していることがわか

るが，実験ではそれらに加えて中間集中点が観察されていることが捉えられる．この観察結果 は数値解析で予測された第三Re数領域の結果に近いものである．これについても後述する．

これまではそれぞれのRe数領域特有の粒子集中現象を個別に説明してきた．つぎにそれらの Re数領域間の変化，つまりどんな粒子平衡位置の分岐現象であるのかに焦点をあてた説明を行

う．2.1.2項で説明したように，粒子平衡位置は揚力が零である位置であり，揚力の動径方向成

分Fr = 0^の等値線(粒子の動径運動に関するヌルクライン)^{と周方向成分}F_θ = 0^の等値線(^周 方向運動に関するヌルクライン)の交点である．これらの等値線のRe数変化を追うことで，安 定な粒子平衡位置の組がそれぞれの臨界Re数を境に移り変わることは粒子平衡位置の分岐現象 として説明できることがわかった．Fig. 12にそれぞれのRe数領域でのF_r = 0およびF_θ = 0 の等値線を模式的に描いた．ここではF_r = 0の青の曲線とそれに沿うようなF_θ = 0の緑の曲 線の二つに注目する．そのため，以降それぞれをr-nullcline^およびθ-nullcline^{として表現する．}

Fr = 0^およびFθ = 0の等値線の数値解析結果はFig. 10^やFig. 11においてそれぞれ青線と緑 線で描いている．

(a)100≤Re<Re_I1 (b) Re_I1<Re<Re_C1 (c)Re_C1<Re<Re_I2 (d)Re_I2<Re≤400

B

A

C B

A

A

θ-nullcline r-nullcline

Figure 12. F_r= 0およびF_θ = 0の等値線のRe数変化の模式図．β = 1/8，100⩽Re⩽400．青 の曲線(r-nullcline)はF_r= 0の等値線である．緑で描いた管壁に沿うような曲線(θ-nullcline)， 各軸と対角線はFθ = 0の等値線に相当する．緑矢印はその付近のFθの符号を表しており，例え ば0< θ < π/4の範囲について，緑の曲線を境に管中央側ではFθ<0^{であり管壁側では}Fθ >0 である．

Fig. 12aに描いたように，第一Re数領域(100⩽Re<Re_I1)ではr-nullclineはθ-nullclineの 内側に位置している．Re数が大きい場合となるとr-nullclineとθ-nullclineはそれぞれ変形し，

Fig. 12bのように断面の軸と対角線の間で二度交差する．これらが平衡位置AおよびBである．

この変化の過渡にはr-nullclineとθ-nullclineが接する場合が存在し，そのRe数がRe_I1である．

第二Re数領域で見られる平衡位置A^とB^はReI1<Reにおいて同時に出現する平衡位置であ ることがわかり，第一と第二Re数領域の間の変化がサドル・ノード分岐として捉えられること がわかった．

Re数がさらに増加するにつれてr-nullclineとθ-nullclineはより複雑な変化を示す．Fig. 12c

に描いたように第三Re^{数領域では，第二}Re数領域での特徴を残したまま，断面四隅付近の部 分でr-nullclineがθ-nullclineより管壁側に位置するようになる．Fig. 12bとFig. 12cの比較か ら，臨界Re数Re_C1でのr-nullclineとθ-nullclineの様子は，対角線上でそれらが接した状態で あることがわかる．Re>Re_C1となり，わずかでもr-nullclineがθ-nullclineより壁側に位置す るとヌルクラインの交差が増加して平衡位置Cが現れる．また，この変化は対角線上平衡位置 の安定性を変える．対角線上平衡位置がθ-nullclineより管中央側に位置するとき(Fig. 12b)，そ の平衡位置周りのFθはそれから離れる向き(^符号)となっているため，この平衡位置がサドル 点であることが捉えられる．Re^数がReC1を超えて対角線上平衡位置がθ-nullcline^{より壁側に} 位置するとき(Fig. 12c)^{，その周囲の}F_θの符号からこの位置が安定であることがわかる．この ように第二から第三Re数領域への集中パターンの遷移は，粒子平衡位置の亜臨界ピッチフォー ク分岐として説明できることがわかった．

(a) (b) (c)

D/2

0

−D/2

Figure 13. (a)^{正方形管での}Poiseuille^{流れ．背景の等値線}: ^{速度の大きさ}upoi，矢印: −∇upoi． いずれの場合においても，暖色ほど値が大きく寒色ほど値が小さい．(b) ^{管内流れに浮遊する} 球に作用する揚力の特徴を考慮して用意した，模擬的な揚力場(r²−cr)∇upoi．r =c^において 揚力が零となる．ここではc= 0.4Dとした．rは管断面中央からの動径距離とする．矢印の色 は，(a)の場合と同様で，大きさを表す．(c) 模擬的な揚力場の周方向成分の等値線．色は符号 を表しており，赤は正(反時計回りの方向)，青は負(時計回りの方向)である．

θ-nullclineを境に管中央側か壁側かどうかで揚力の周方向成分F_θの符号が変化するという特

徴(例えば0 < θ < π/4領域においてθ-nullclineの内側ではF_θ < 0，外側ではF_θ > 0)に関 して補足する．Fig. 12a–dのどの場合でも見られるこの特徴は，shear-gradient induced liftと

wall effectによって説明される．球に作用する揚力はその断面内位置でのPoiseuille流れの速度

勾配の方向に対して概ね平行であり反対の方向(管壁方向，Fig. 13a矢印を参照)となることが shear-gradient induced liftから考えられる．そのため，管断面0< θ < π/4^{領域の管中央付近} ではFθ <0となる．揚力の方向は管壁に付近になるにつれて管壁から離れる方向となる(wall effect)．その方向の周方向成分に注目すると，管壁付近でのF_θ^の符号は0< θ < π/4^領域では 正となることがわかる．視覚的な理解のために，Poiseuille流れの速度勾配を管中央付近でその