博士論文

粉体流動層の粒子ダイナミックス

東北大学大学院理学研究科物理学専攻市來健吾平成８年

概要

本研究では粉体流動層の数値解析を行う。本研究の特徴は、これまでに開発されて来た手法のほとんどが現象論的に流体と粒子の相互作用を導入したのに対し、この流体力学的相互作用を正確に扱いモデルの構成を行った点である。このようなモデルを構成する為に、粒子集団を連続体近似で扱うのではなく、個々の粒子の運動を直接扱う。また現実の流動層から最小限の基本的なメカニズムのみを取り出して「理想流動層」を定義し、この系をモデル化する。これは、現実の粉体流動層は一般に粒子の大きさや系に流し込む流体の流速の範囲が広く、これら全てを包括的に記述する第一原理的なモデルの構成は不可能である為である。

粒子と流体の相互作用は本来流体の粘性による摩擦の効果である為、理想流動層では粘性が支配的な低レイノルズ数の極限で流体を記述する。粒子間の直接相互作用は剛体的で、衝突は弾性とした。また粒子は全て同じ大きさの球とし、粒径や質量、形状の分布は無視した。低レイノルズ数を仮定している為、流入流速が大きな場合の挙動、例えば希薄流動状態や輸送状態などを記述するのは不可能である。従って現実の流動層の示す現象の中でも静止状態から流動状態への流動化現象と、気泡流動状態を含む比較的穏やかな流動状態が直接の考察対象となる。

モデル化で最も重要な低レイノルズ数での粒子間相互作用は、支配方程式の線形成の為に、粒子速度と流体に及ぼす力とを関係づける行列によって表される。この行列は粒子濃度が希薄な場合や２粒子のみの場合は求められている。今考察しようとする粒子濃度が高く多くの粒子が存在する場合は、希薄な場合と２粒子の場合の解をセルフコンシステントに結合する方法が良い結果を示す事が、コロイド粒子系の数値シミュレーションに於いて示された。本研究ではこの近似の正当性を確認する意味から、沈澱現象での平均沈降速度の解析にこの近似を用い、既存の解析が改善される事を示す。

粉体はコロイド粒子と異なり粒子の慣性が重要となる。しかし低レイノルズ数での粒子間相互作用を用いると、この相互作用の近接効果の持つ特異性 (抵抗の発散) の為に、モデルに有効に慣性を導入できないことが分かる。本研究ではこの特異性を繰り込んだ慣性を用いる事でモデルを構成した。この繰り込みが現実の系の振舞を良く再現する事は、現実の系では流体の連続性の限界や粒子表面での凹凸などにより流体力学的相互作用の特異性が隠されることを意味する。

このモデルを特徴付けるパラメータは、流入速度と無次元化された粒子質量であるストークス数である。このパラメータ空間で系統的に数値シミュレーションを行う。この結果、流入速度に流動化相転移が生じる臨界流速が存在する事が分かる。流動層の膨張率と流入速度の関係も実験結果を定性的に再現する。また流動状態はチャンネル状態と気泡状態の２種類があり、ストークス数がこれらの状態を特徴付けることが分かる。

粒子の運動エネルギーの振舞から、流入速度が「温度」の役割を持つ事が示唆される。この「温度」を用いて揺動散逸定理 (アインシュタインの関係) を理想流動層に適用し、自己拡散係数から有効粘性率を定義する。この有効粘性の流入速度依存性が、現実の実験結果と無矛盾である事が示される。また有効粘性率はストークス数依存性を含めて速度分布の非ガウス性の強さと極めて良く一致している事が分かる。

これらの結果は液体論で用いられる空孔模型で説明する事が出来る。つまり理想流動層の定常状態と液体の間に強い類似性が示される。

目次

1 はじめに 1 1.1 何を研究するのか . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 研究の背景 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 現実の流動層と理想流動層 4

2.1 現実の流動層 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.1 流動層の構成要素 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.2 流動層の示す現象 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.3 流動状態の分類 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2 流動層の研究 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 理想流動層 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.3.1 流体力学的相互作用 . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3.2 粒子間衝突 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.3 理想流動層の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3 シミュレーション方法 11 3.1 粒子の運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 流体力学的相互作用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3.2.1 ストークス近似 (低レイノルズ数近似) . . . . . . . . . . . . . . . . 12 3.2.2 希薄な場合の相互作用行列 . . . . . . . . . . . . . . . . . . . . . . . 16 3.2.3 ２体問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2.4 高濃度での相互作用行列 . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2.5 周期境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.3 粒子慣性の導入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3.1 粒子慣性の無い場合 . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3.2 粒子慣性の導入の困難 . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3.3 lubrication の繰り込み . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3.4 無次元化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3.5 繰り込みの物理的意味 . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.3.6 固定粒子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4 計算結果と解析 39

4.1 シミュレーションの設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1.1 パラメータ空間 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1.2 シミュレーションのサイズについて . . . . . . . . . . . . . . . . . . 41

4.2 動的挙動 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.2.1 振動的挙動と気泡 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.2.2 気泡状態とチャンネル状態 . . . . . . . . . . . . . . . . . . . . . . . 44 4.2.3 気泡生成周期 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2.4 気泡の上昇速度の流入速度依存性 . . . . . . . . . . . . . . . . . . . 46

4.3 定常状態 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.3.1 平均運動エネルギー . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.3.2 流動層の平均重心 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.3.3 自己拡散係数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.3.4 有効粘性率 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.3.5 速度分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

5 議論 56

5.1 シミュレーション結果の解釈 . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.1.1 ストークス数依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.1.2 流入速度依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.1.3 液体とのアナロジー . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.2 一様状態の欠如 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

6 結論と課題 63 6.1 結論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 6.2 課題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

Chapter 1

はじめに

1.1 何を研究するのか

花瓶か水槽のような容器の中に、粒の揃った砂か小さなガラスビーズを入れてみる。これだけでは砂は容器の底に堆積するだけである。この砂は静止状態にあると言える。容器の底が、砂は通さないが水や空気が通る位の小さな穴が沢山開いている多孔質性の物質で出来ているとし、この容器の底から上向きに流体を流し込む。流れの速さがある値を越えると、重力で底に堆積し互いに動く事の出来なかった砂は吹き上げられ、ある程度自由に動き回れるようになる。この砂は流動状態にあると言える。ここで扱った砂やガラスビーズは「粉体」と呼ばれ、粉体が静止状態から流動状態に変化する現象は「流動化」と呼ばれる。この「粉体」の「流動化」が本研究の研究対象である。

近年、粉体は統計物理学の研究対象として盛んに研究されるようになってきた。粉体を構成する粒子はマクロな大きさを持つ。従って変形や振動などの散逸効果が重要となって来る。その為例に示したように、粉体は外部からエネルギーを供給しないと流動化しない。エネルギーの供給の仕方には様々な方法がある。典型的なものとしては、本研究で考察する流体の流れによって駆動される「粉体流動層」の他に、容器の振動によって駆動される「粉体振動層」と呼ばれる系もある。また地震の時、振動と地下水によって生じる地盤の液状化現象として広く知られている現象も、粉体の流動化現象の１つと考える事が出来る。

本研究では「粉体流動層」を現象論を極力排し粒子スケールでのモデルを構成する。またこのモデルを用いて「流動化現象」を数値的に解析する。この研究により、これまで看過されて来た流動層の粒子スケールでの振舞と、粉体現象に於ける流体の寄与の一般的な考察を目指す。

以下本章では本研究の背景として、粉体全般に対する最近の物理的関心を中心に紹介する。

1.2 研究の背景

粉体とは、例に示した砂のような粒 (つぶ) 状の物質の集団を指すが、「粒径が何センチまでで、密度が...」と言うような定義はない。おそらく最も正しい粉体の定義は、「散逸を持つ粒子の集団」である。その意味で、例に挙げた「地盤」はもとより、台所の砂糖や机上の鉛筆、また街を流れる人間や車など、我々の身近にある物の多くは「粉体」と見倣す事ができる [1]。

そもそも「粉体」が物理学者の関心を集めたのは、砂山崩しの数値模型が自己組織化臨界現象を示した事が発端であろう [2]。それを契機に、改めて粉体の示す現象を物理的に見て行くと、他の物質形態とは異なる特異な振舞を示す事が分かって来た [3, 4]。

粉体の最も特徴的な点は、通常の意味での「温度」が効かない事と、「散逸」が最も重要なメカニズムである事である。これは考えている粒子がマクロな大きさをもつ事による。なぜならば、この場合は熱的なエネルギーよりも重力などの外力によるエネルギーが十分大きく、また粒子の変形や内部自由度が無視できないからである。粉体には温度が効かない為、通常の気体や液体が従う熱力学とは全く異なった振舞を示す。例えば、２種類の砂を混ぜた系を外から振動を加えて流動化させると、大きな粒子と小さな粒子が分離する現象を示す。これは通常の熱力学に反する振舞である。

このように粉体を入れた容器を振動させる事により系を流動化させた系は「粉体振動層」と呼ばれ、古くはファラデーによる研究に遡る事が出来る [5]。振動層の示す特徴的な現象には、振動加速度を上げる事により生じる流動化とそれに伴う粒子対流がある [6, 7]。また薄い振動層では、液体の表面波に非常に良く似たパターンが観測される [8, 9]。

最近１０年の間に様々な「特異な粉体の振舞い」に対する各論的研究が成され、それぞれの現象に対する理解は深まって来たと言える。最近、これらの研究成果から「粉体の物理学」として何か普遍的な物理を引き出そうとする機運が高まって来ている。しかし粉体の示す現象は幅広く、理解が進んだのはほんの一部に過ぎない。実際、一般に粉体現象は大気中で行われているにも拘らず、これまでのほとんどこの解析では大気の効果が無視されて来た、或は大気の効果が無視できる状況が主に研究されて来た、と言える。粉体流動層はまさにこの粒子を取り囲む流体の効果が主役になる系である。流動層に対する物理的理解が振動層のそれに比べると進んではいない事がこの状況を反映している。

流体が絡んで来る粉体現象が避けられて来た理由は、流体の効果が複雑である点にある。粉体現象の時間、空間スケールは広く、そのすべてを包括的に扱おうとすると、流体力学の非線形性に正面から取り組まなくてはならない。これは簡単に粉体に応用出来る問題ではない。

1.3 構成

本論文の構成は以下の通りである。

２章では実際の流動層に見られる現象を概観し、既存の「流動層」の解析手法を紹介する。既存の手法と本研究の根本的な違いは、流動層で最も重要なメカニズムである流体と粒子の相互作用の取り扱いにある。本研究ではこれまでのようにこの効果を現象論的に導入するのではなく、流体力学の境界値問題として評価する事でモデルの構成を行う点にある。この目的の為現実の流動層から本質的であると思われる最小限のメカニズムを取り出し「理想流動層」を構成する。

３章では、「理想流動層」の数値モデルの具体的な定式化を行う。本研究の成果の１つとして、粉体のような粒子慣性が重要となる系に流体力学的相互作用を第一原理的に導入したモデルを構成した事が挙げられる。また高濃度での流体力学的相互作用の計算方法を応用して平均沈降速度の理論解析を行い、実験結果と良く一致する結果を得る事を報告する。

４章では構成したモデルを用いて、「流動化」現象の系統的なシミュレーションを行う。その結果から「静止状態」から「流動状態」へ移る流動化転移が存在する事を示す。さらに流動状態を定量的に解析する事で、「気泡状態」と「チャンネル状態」に対応する２つの状態が存在する事が分かる。またこれらの転移現象と粒子の運動エネルギーや拡散係数、有効粘性のパラメータ依存性について結果を示す。

５章は議論として、シミュレーションで得られた定常状態での系の振舞いに対する解釈を行い、液体との類似性を議論する。また、現実の系では観測される一様状態が、理想流動層のシミュレーションに見られなかった点について、その原因を議論する。

最後に６章で結論と今後の課題をまとめる。

現実の流動層と理想流動層

本章では流動層が示す実際の現象を紹介し、これまで研究された解析手法をまとめる。後半では、現実の流動層から本質的なメカニズムのみを取り出して「理想流動層」を構成する。

2.1 現実の流動層

2.1.1 流動層の構成要素

そもそも流動層は化学工業の反応炉として工業的に極めて重要な装置である [10]。そのため２０世紀前半に実用化されて以来、常に経験が (工学的) 理論に先行して来た歴史を持つ。

実際の流動層はそれぞれの目的に応じて、流体の流入機構や粒子の再循環装置など特別な仕組みが付加される。しかし流動層の基本構成は、粒子は通さないが流体は流し入れる事のできる分散器が底に取り付けられた縦長の容器と、その中に入れられた粉体からなる。この分散器としては、例えば小さな穴やノズルを持った平板や多孔質性のフィルターなどが用いられる。

この装置の大きさ (容器の断面の大きさ) は、実験的な装置では数センチ程度のものから、実際に工業的に使われるプラントでは数メートルのものまで様々である。

一方用いられる粒子は、粒径が数十ミクロンから数ミリ程度、粒子と流体の密度差が 0.1 ∼ 10[g/cm^{3}] 程度のものである。実際に幾つかの典型的な物質の密度を示しておく。

\begin{matrix}物質 & 密度 [g/cm^{3}] \\ 砂 & 1.4 ∼ 1.7 \\ 石英 & 2.22 \\ ガラス & 2.4 ∼ 2.6 \\ アルミ & 2.69 \\ 銅 & 8.93 \\ 鉛 & 11.34 \\ 水 & 1.0 \\ 空気 & 1.2 × 10^{−3} \end{matrix}

流体として用いられる典型的な物質は、気体としては空気、液体としては水である。

2.1.2 流動層の示す現象

流動層は流体の流入流速を変化させる事により、様々な状態を示す事が実験的に知られている。

流速が十分に小さい時は、粒子は重力により容器の底に堆積した「静止状態」にある。流速がある値 (流動化開始速度) を越えた所で粒子は流動化を起こす。流動層の流動状態

には様々なものがある (図 2.1)。

Figure 2.1: 様々な流動状態。

流体として液体を用いた流動層では粒子が一様に膨張した「一様状態」が観測される。一方空気による流動層の流動状態には「気泡状態」と呼ばれる気泡 (空気の塊) が層の中を上昇して行く状態も観測される。この気泡は流速を上げて行くと次第に大きくなり、容器のサイズにまで大きくなった「スラグ状態」と呼ばれる状態に移行する。スラグとは横方向に広がった気泡の事である。流速の上昇に従い一般に粒子濃度は低くなって行き、十分に大きい流速では粒子が流体によって輸送される「輸送状態」に至る。極めて粒径の小さい粒子 (数ミクロン程度) を用いた場合は粒子間の凝集力が重要になる為、流動化しにくくなる。このような系では「チャネル状態」と呼ばれる流体の通る流路が形成される状態が観測される。

2.1.3 流動状態の分類

流動層の示すいろいろな流動状態間の遷移は流速の変化によって生じるが、どの状態を変化して行くかは流体や粒子の物性に依存する。流動状態の分類の１つの試みとして、物質パラメータの違いをアルキメデス数 Arで代表させ、流速をレイノルズ数 Re で表し、 (Re, Ar) 空間での相図を書くことができる (図 2.2)。ここでアルキメデス数 Arとレイノルズ数 Re は以下のように定義する。 (本章での見積りでは、粉体工学の分野での慣習に従って、長さのスケールに粒子半径 a ではなく直径 d を用いる。 )

Ar = \frac{d^{3}ρ_{f}(ρ_{p} − ρ_{f})g}{μ^{2}}  (2.1)
Re = \frac{dρ_{f}u^{∞}}{μ}  (2.2)

ここで g, μ, ρ_{p}, ρ_{f} はそれぞれ重力加速度、流体の粘性率、そして粒子と流体の密度である。アルキメデス数は、以下で出て来る本研究で用いるモデルを特徴付ける無次元パラメータ

Figure 2.2: 流動化状態図。 (文献 [11] より転載。 )

であるストークス数 St と以下の関係がある。

St ≡ \frac{mU_{0}}{6πμa^{2}} = \frac{1}{162} \frac{ρ_{p}}{ρ_{f}} Ar  (2.3)

ここで m は粒子の質量、 U_{0}は１粒子の沈降速度である。

U_{0} = \frac{mg}{6πμa} \frac{ρ_{p} − ρ_{f}}{ρ_{p}}  (2.4)

つまりストークス数は粒子の慣性エネルギーと流体の粘性エネルギーの比である。この相図で最も注目したいのは、流動化現象がレイノルズ数の小さい領域でも実験的に観測されている事である。この事実は本研究での流体力学的相互作用の定式化の大きな拠り所である。

但し現実の流動層の振舞いはこの相図に示されるように物質パラメータが１つだけで特徴付けられる様な簡単なものではない。端的に言えば粒子の大きさも重要なパラメータであり、実際粒径と、粒子と流体の密度差によって粒子を分類する試みも成されている [12]。

2.2 流動層の研究

「流動層」を研究する場合、粉体と流体の２種類の構成要素を扱う必要がある。この系を記述するスケールの違いにより、幾つかのアプローチが考えられる。

スケールモデル特徴
巨視的２流体モデル粒子を連続体とする
運動論モデル粒子の挙動を確率的に記述
微視的離散要素法粒子と流体の相互作用は粗視化
本研究理想流動層の第一原理モデル

最も巨視的な立場では、粉体も流体と同様に連続体と見倣して系を記述する。この場合流動層を２成分流体系として、幾つかの基本的な仮定のもとに支配方程式を構成していく [13, 14, 15, 16, 17]。この立場の長所は、モデルが解析的に扱いやすい事や、微細な構造を直接扱わない為に流動層に限らない普遍的な性質を議論できる点などがある。しかし半面、そもそも離散系である粉体を連続体として記述する事の妥当性に対する疑問は残る。仮に連続体として記述するとしても、粉体は通常の流体とは著しく異なった流体として振舞う (例えば粘弾性を示す) 為に、モデルの構成時に現象論を導入せざるを得ない。この現象論、特に粒子相の応力に対する構成方程式が自明ではなく、多くのモデルが存在している点が問題である。

一方より微視的な立場でのモデルも幾つかある。この立場では粉体を離散的な粒子として直接記述する。

運動論モデルは、気体の分子運動論の手法を粉体系に応用した手法で、ボルツマン方程式のように粒子の挙動を確率的に表現する [18]。このアプローチはもともと希薄気体に対して開発された為、流動層のような極めて高い粒子濃度での適用するという方法論自身に困難がある。また流体と粒子の相互作用の導入は、運動論の枠内では成功しておらず、２流体モデルと同様現象論を別に用意しなければならない。

離散要素法 (DEM) [19, 20, 21] は粒子間の直接相互作用を、弾性衝突を表すバネ、散逸を表すダッシュポットやスライダーなどの力学的なモデルで記述する。一方流体との相互作用は粒子相の粗視化を行い、２流体モデルと同様の現象論的な流体相との摩擦を導入し、また流体相は粘性を無視した方程式を用いている。この手法は現象を非常に良く再現しており、工学の分野での要求に最も叶った手法である。しかし導入されるパラメータの数は多く、用いる値も実験結果を使うなど、実験と独立な解析とは言えない。

このように、これまでの解析は流動層で最も本質的で重要な流体と粒子の相互作用に曖昧な現象論を用いている。本研究ではこの流体と粒子の相互作用を流体力学の境界値問題として正確に取り扱う事でモデルの構成を試みる。この目的の為、現実の流動層から本質的ではないと思われる複雑性を極力排除して「理想流動層」を構成し、この理想流動層を直接の研究対象とする。

2.3 理想流動層

2.1節で見たように流動層の現象は多岐に渡り、全体を包括する普遍的なモデルは一般には構成不可能である。本研究では、現実の流動層から本質的であると思われる基本的なメカニズムのみを取り出して「理想流動層」を構成し、その系を数値的に解析する。その際、流体の流入速度が比較的小さい場合の「流動化」現象を念頭に置きモデルを構成する。

流動層を構成する必要不可欠で考慮すべき効果には、流体を介しての粒子間相互作用、粒子の衝突、粒子の慣性の寄与、重力が挙げられる。

2.3.1 流体力学的相互作用

はじめに流動層に於いて最も重要である流体力学的相互作用について考える。着目する粒子の挙動の時間スケールが流体の音速に比べて遅い場合、流体は非圧縮と考えられる。この時流体はナビエ・ストークス方程式に従い、相互作用は粒子表面での境界値問題を解いて得られる。流動層に於ける流体力学的相互作用は、流体の流れで粒子が引きずられる効果である。これは流体の流れから粒子が受ける抵抗であり、流体の粘性が示す効果である。従って２流体モデルや有限要素法で仮定されるように、流体を非粘性とするのは正しくない。むしろ粘性効果こそ本質的と言える。本研究では、流体力学的相互作用を粘性が支配的である低レイノルズ数の極限で扱う。この理想化された条件は現実の流動層で実現するのは難しいと思われるかも知れない。しかし 2.1.3節の図 2.2からも分かる通り、レイノルズ数が比較的小さい流入速度でも系は流動化し、気泡の発生も実験的に観測されている [10, 11]。

流体力学的相互作用の実際の定式化に関して考える。粉体流動層は容器の中に粒子を入れ、流体を容器の底から流し込むことで構成される。従って系は束縛された流体中で相互作用する多粒子系である。ここで本質的なのは、流そうとする流体が粒子を避ける事が出来ない事である。この状況は空間的に束縛された流体中の粒子の挙動を考えるか、流体、粒子とも空間に無限に分布している系を仮想的に考えなければならない。本研究では後者の方法で系を構成する。これは、実際に数値シミュレーション可能な定式化が現在のところ後者の状況に限られてしまう為である。この選択の持つ意味は、現実の系に存在する容器 –特に側壁–の効果が効かない現象を見ようとしていると言える。低レイノルズ数の仮定から本研究の対象は流速の比較的遅い時の流動層の示す挙動に限られる。この場合、観測される気泡は容器の大きさに比べて十分小さい。この為側壁の効果を無視する事は正当化される。流体と粒子がともに空間に無限に分布する系は周期境界条件を用いて記述する。一方粉体を支える容器の底は、空間的に固定された粒子を用いてモデルに導入する。

2.3.2 粒子間衝突

次に粒子の衝突に関して考える。弾性衝突する剛体粒子系は古くから気体の分子模型として用いられて来た。数値シミュレーション法である分子動力学 [22] の結果から、高い粒子濃度での固相-液相相転移である Alder 転移を除き、系は特異な振舞は示さない。一方粉体は巨視的な大きさをもつ為、衝突により粒子のエネルギーは散逸する。つまり粉体の衝突は非弾性である。衝突が非弾性である系は、弾性衝突する系と異なり興味深い振舞を示す。実際、粉体振動層は非弾性粒子系と振動する容器によってモデル化する事ができるが、この系の数値シミュレーションから、ある振動加速度以上で流動化相転移が観測される事が示された [6, 7]。また最も単純な例である真空中での非弾性粒子系、つまり空気など粒子を取り巻く流体の効果を一切無視した系の研究から、非弾性を特徴付ける反発係数の大きさにより一様状態が不安定になる事が分かった [23]。

流動層は散逸系である。つまり流体の流れによりエネルギーが系に供給され、系自身の持つ散逸のメカニズムによって与えられたエネルギーが消費される。両者がバランスする事で系は定常状態に至る。流動層の場合、散逸のメカニズムとして第一に粒子表面での粘性境界条件を通しての流体との相互作用が考えられる。現実の流動層では当然それに加え非弾性衝突も考えられる。実際これら２つの効果のカップリングは興味ある問題であるが、本研究の趣旨から非弾性効果は二次的な役割であると仮定し、理想流動層での粒子間の衝突は弾性的であるとする。

2.3.3 理想流動層の定義

以上の議論から、本研究では理想流動層を、低レイノルズ近似で周期境界条件を用いて定式化した流体力学的相互作用、粒子間の弾性衝突、粒子の慣性の寄与、重力で構成する。

ここで本質的でないとの仮定により無視された効果としては、次のようなものが考えられる。先の議論によって排除された有限レイノルズ数効果、流体を束縛する境界としての容器の効果、衝突に於ける非弾性効果。加えて単純化の為本研究では、粒子の形状と多分散性、粒子の回転、分子間力など粒子間の (衝突以外の) 直接相互作用を無視した。

これらの効果は、より現実的な状況を目指す場合重要となる効果である。先に指摘した通り、高レイノルズ数効果は流入速度の比較的速い領域では無視する事は出来ない。現実の系では粒径が小さくなると粒子間の静電気力や分子間力による凝集力が無視できなくなる。

しかしながら、これらの副次的な効果を排除した理想流動層を研究する事は、流動層に対するより本質的な理解を進めるものと期待できる。

シミュレーション方法

本章では前節で定義した理想流動層の具体的な数値モデルを構成する。

3.1 粒子の運動方程式

本研究では粒子スケールでの数値モデルを構成する。従って粒子の運動方程式が支配方程式である。理想流動層の粒子の従う運動方程式は以下のように書ける。

m\frac{dU}{dt} = F_{f} + F_{g} + F_{c}  (3.1)

ここで m, U, F_{f}, F_{g}, F_{c} はそれぞれ粒子の質量、粒子の速度、粒子が流体から受ける力、重力、衝突力である。ここでボールド体はベクトルを表し、またこのベクトルは粒子の次元も含むとする。つまり３次元空間で N粒子の速度ベクトル U は以下のような 3N次元ベクトルである。

U = [ U^{(1)} ; U^{(2)} ; ⋮ ; U^{(N)} ]  (3.2)

上付添字は粒子を表す。運動方程式 (3.1) 式の各項の内、重力 F_{g}は以下のように書ける。

F_{g} = −E_{z}mg \frac{(ρ_{p} − ρ_{f})}{ρ_{p}}  (3.3)

ここで g, ρ_{p}, ρ_{f}は重力加速度と、粒子と流体の密度である。また E_{z}は (3.2) 式の形に拡張した z方向の単位ベクトルで、重力は−z方向に働くとした。

衝突力 F_{c}はハードコアとし、また完全弾性衝突を接触時の運動量交換の関係式で表現する事で、実際の数値シミュレーションに於いては力としては導入しない。

運動方程式 (3.1) 式の残りの２つの項に関しては以下の節で詳細に議論する。

3.2 流体力学的相互作用

本節では運動方程式 (3.1) の粒子が流体から受ける力 F_{f} を議論する。以下簡単に内容をまとめる。

まず最初に理想流動層の構成の際に仮定した低レイノルズ数近似を具体的に検証し、この近似のもとで流体から受ける力 (粒子間相互作用) がレジスタンス行列と呼ばれる量で記述できる事を示す [24]。次に粒子濃度が希薄な場合 (粒子間隔が大きい場合) に対してこの行列を具体的に書き下す。２粒子間の相互作用は厳密に計算する事ができる [25]。この結果を次に示す。その後、粒子濃度が高い場合の行列の近似的な計算手法である Stokesian dynamics の方法を紹介する [26, 27, 28]。シミュレーション法として開発されたこの方法を平均沈降速度の問題に応用し、既存の理論の高濃度での振舞を改善する事を示す [29]。理想流動層では容器を導入する代わりに周期境界条件を課すが、このときの行列の計算を次に紹介する [30, 31, 32]。

3.2.1 ストークス近似 (低レイノルズ数近似)

非圧縮性

粒子の特徴的な速度が流体の音速に比べて十分小さい時、流体を非圧縮として扱う事が許される。本節では流動層で用いられる２つの典型的な流体である水と空気について、この非圧縮性の仮定を検証する。

水と空気の密度、粘性率、音速はそれぞれ以下のような値である。

\begin{matrix}流体 & 密度 [g/cm^{3}] & 粘性率 [g cm^{−1}s^{−1}] & 音速 [cm/s] \\ 空気 & 1.2 × 10^{−3} & 1.82 × 10^{−4} & 3.31 × 10^{4} \\ 水 & 1.0 & 1.0 × 10^{−2} & 1.5 × 10^{5} \end{matrix}

流動層での特徴的な速度は、流体の流入速度である。この流入速度は流体と粒子の相対速度の平均値と見倣す事ができる。１粒子が単独で流体中を落下 (沈降) する速度以上の流速を流動層に与えると、粒子は全て吹き飛ばされてしまう^{1}。このため、流入速度にはこの１粒子の沈降速度 U_{0} と言う上限が存在する。レイノルズ数が小さい場合 U_{0}は次のように書ける。

U_{0} = \frac{mg}{6πμa} \frac{(ρ_{p} − ρ_{f})}{ρ_{p}}  (3.4)

ここで m = (4/3)πa^{3}ρ_{p}は粒子の質量、 ρ_{p}, ρ_{f} はそれぞれ粒子と流体の密度、 μは流体の粘性率、 a は粒子の半径である。レイノルズ数が小さいという条件は本研究の解析では常に仮定される条件である。粒子のスケールでのレイノルズ数は一般に次のようになる。

Re = \frac{ρ_{f}aU}{μ}  (3.5)

ここで Uは特徴的速度である。低レイノルズ条件 Re < 1 から、粒子半径 a、特徴的速度 Uに上限が定められる。特徴的速度 Uの上限として先の１粒子の沈降速度 U_{0}を用いると、この条件下での粒子半径の上限が決まる。従ってこの最大粒子半径 a_{max}での１粒子の沈降速度 U_{0}(a_{max}) は、低レイノルズ数の限界での最大の特徴的速度である。この速度が流体の音速に比べて十分小さいことは、実際に計算すると分かる。この見積もりで、粒子の密度を典型的なガラスの値 ρ_{p} = 2.5[g/cm^{3}] を用いた。

流体 a_{max}[cm] U_{0}(a_{max})[cm/s]  音速 [cm/s]
空気 3.70 × 10^{−3}  4.10 × 10^{1}  ≪ 3.31 × 10^{4}
水 6.74 × 10^{−3}  1.49 × 10^{0}  ≪ 1.5 × 10^{5}

従って低レイノルズ数近似の範囲では流体の圧縮性は無視できる。

線形化

非圧縮粘性流体は一般にナヴィエ・ストークス方程式で記述できる。

ρ_{f} {∂_{t}u + (u · ∇)u} = μ∇^{2}u − ∇p  (3.6)
∇ · u = 0  (3.7)

(3.6) 式を、特徴的速度 U、特徴的長さ L、特徴的時間スケール T で無次元化する。

\frac{ρ_{f}LU}{μ} \{\frac{L}{UT} ˆ∂_{t}û + (û · ˆ∇)û\} = ˆ∇^{2}û − ˆ∇p̂  (3.8)

1流速が大きい場合、粒子は吹き飛ばされ輸送される。この輸送状態の粉体の挙動も (特に工業上) 興味深い研究対象であるが、本研究では触れない。

ここで圧力は p = (L/μU)p̂と無次元化した。 (3.8) 式に現れる２つの無次元量はそれぞれ、レイノルズ数 Re、ストローハル数 Sh と呼ばれる量である。

Re = \frac{ρ_{f}LU}{μ}  (3.9)
Sh = \frac{L}{UT}  (3.10)

レイノルズ数 Re は、流体の慣性と粘性の比を表している。粒子間相互作用を考える場合、特徴的長さスケールは L ∼ a なので前節で考察したレイノルズ数 (3.5) に一致する。一方ストローハル数 Sh は、外部から与えられる時間スケールと、特徴的速度と特徴的長さから決まる時間スケールの比を表す。一定流速下での流動層の定常状態を考える場合、この外部から与える時間スケールは特に無い。従って低レイノルズ数 Re ≪ 1 の場合、流体は次の方程式に従うと近似できる。

0 = μ∇^{2}u − ∇p  (3.11)
∇ · u = 0  (3.12)

この近似をストークス近似と呼ぶ。この時流体の方程式は定常となるが、今の問題では粒子表面での境界条件を通して時間発展が入って来る。流体の方程式が定常である事は、流体の慣性が小さい為、粒子の移動に伴う境界の変化に比べて流体が追随する時間スケールが速い事を反映している。

行列表現^{2}

ストークス近似のもとでの粒子の流体力学的相互作用を考える。その為に無限に広がった流体中に N 個の粒子が分散している系を考える。境界条件は、粒子表面 S_{α}と無限遠方 S_{∞}で以下のように与える。

\begin{aligned}u & = & U^{α} & on & S_{α} \\ u & = & 0 & on & S_{∞} \end{aligned}  (3.13)

支配方程式の線形性から、速度場 u と圧力場 p を次のように分解する。

u = ∑_{α=1}^{N} u^{α}  (3.14)
p = ∑_{α=1}^{N} p^{α}  (3.15)

2本節の内容は以下の文献に依る [24]。

(u^{α}, p^{α}) は以下の境界条件を満たす解とする。

u^{α} = \{ \begin{matrix}U^{α} & on S_{α} \\ 0 & on S_{∞}, S_{β≠α} \end{matrix}  (3.16)

つまりこれらの解は、 α 粒子のみ運動し、残りの粒子は全て静止している時の解を表す。当然これらの解の重ね合わせは元の問題の境界条件 (3.13) を満たす。この解の分解は形式的なものであり、 N 体問題を１体問題に帰着させた訳ではない点は注意するべきである。

更に線形性から、 (u^{α}, p^{α}) は U^{α}に比例するはずである。

u^{α} = V^{α} · U^{α}  (3.17)
p^{α} = μP^{α} · U^{α}  (3.18)

ここで比例係数として導入した V^{α}は２階テンソル、 P^{α}は１階テンソルである。これらは、もとの方程式に代入する事で U^{α}, μに依存しない事が分かる。

一般に速度場 u と圧力場 p から流体のストレステンソルΠが構成できる。

Π = −pI + μ \{\frac{∂u}{∂x} + (\frac{∂u}{∂x})^{†}\}  (3.19)

ここで I は単位テンソル、^{†}は転置を表す。このストレステンソルをβ 粒子表面で積分する事で粒子が流体に及ぼす力 F^{β} が求められる。

F^{β} = − ∫_{S_{β}} dS Π · n  (3.20)

ここで n は表面 Sの流体を向いた法線ベクトルである。解 (u^{α}, p^{α}) に同様の操作を行って力を計算してその結果をまとめると、粒子が及ぼす力と粒子の並進速度の間に以下のような行列関係式が書き下せる。

[ F^{1} F^{2} ⋮ ] = μ [ \begin{matrix}A^{11} & A^{12} & ⋯ \\ A^{21} & A^{22} & ⋯ \\ ⋮ & ⋮ & ... \end{matrix} ] · [ U^{1} U^{2} ⋮ ]  (3.21)
A^{αβ} = − ∫_{S_{α}} Σ^{β} · dS  (3.22)

ここで、 Σ^{β}は３階テンソルで、粘性率μ、粒子の並進速度 U に依存しない量である。

Σ^{β} = −P^{β}I + \{\frac{∂V^{β}}{∂x} + (\frac{∂V^{β}}{∂x} )^{†}\}  (3.23)

粒子が流体に及ぼす力 F と粒子の並進速度 U を結び付ける行列をレジスタンス (抵抗) 行列 R と呼ぶ。

F = R · U  (3.24)

この行列は粘性流体中での粒子間の流体力学的相互作用を表し、粒子の位置にのみ依存するような量である事が分かる。この線形方程式 (3.24) は逆に解く事も可能である。つまり速度を力に関係づける行列も作る事が出来る。

U = M · F  (3.25)

この行列はモビリティ(易動度) 行列 M と呼ばれ、レジスタンス行列の逆行列である。

M = R^{−1}  (3.26)

ストークス方程式のガリレイ不変性から、 (3.24) 式は一様流 u^{∞}の中の粒子系に対しても成立する。但し粒子速度は一様流 u^{∞}との相対速度に代わる。

F = R · (U − u^{∞})  (3.27)

u^{∞}は粒子が存在しない時の流体の流速であり、流動層で一般に使われる「空塔速度」と同一である。

本節での議論は形式的なものであり、流体力学的相互作用の計算、つまり粒子表面での境界値問題が、ストークス近似の下ではレジスタンス行列を計算する事と等価である事が示されたに過ぎない。以下ではこのレジスタンス行列を具体的に計算する方法を議論する。

3.2.2 希薄な場合の相互作用行列

本節では粒子濃度が希薄な場合のモビリティ行列を導出する。その為に最初ポイントフォースに対するストークス方程式の解としてオセーンテンソルを示し、粒子の作る速度場をポイントフォースの多重極展開を用いて表す。粒子が流体に及ぼす力を速度場から求める為にファクセンの法則を示して、最後にモビリティ行列の形にまとめる。

オセーンテンソル

ストークス方程式の基本解であるポイントフォースに対する解 (グリーン関数) を求める。ストークス方程式は次のように書ける。

−μ∇^{2}u(r) + ∇p(r) = δ(r)f  (3.28)

∇ · u(r) = 0  (3.29)

速度場 u(r) と圧力場 p(r) をフーリエ変換する。

u(r) = \frac{1}{(2π)^{3}} ∫ dk ũ_{k}e^{ik·r}  (3.30)
p(r) = \frac{1}{(2π)^{3}} ∫ dk p̃_{k}e^{ik·r}  (3.31)

(3.32)

従って波数空間でのストークス方程式は次のように書ける。

μk^{2}ũ_{k} + ikp̃_{k} = f  (3.33)
ik · ũ_{k} = 0  (3.34)

これは ũ_{k}, p̃_{k} について以下のように解ける。

ũ_{k} = \frac{1}{μk^{2}} (I − \frac{kk}{k^{2}}) · f  (3.35)
p̃_{k} = −i\frac{k}{k} · f  (3.36)

ũ_{k}, p̃_{k} を実空間に戻すと原点に存在する強さ fのポイントフォースに対するストークス方程式の解が以下のように書ける。

u(r) = \frac{1}{8πμr} (I + \frac{rr}{r^{2}}) · f  (3.37)
p(r) = \frac{1}{4π} \frac{r}{r^{3}} · f  (3.38)

この解はストークス方程式の最も基本的な解であり、特に以下のテンソルをオセーンテンソルと呼ぶ。

J(r) = \frac{1}{r} (I + \frac{rr}{r^{2}})  (3.39)

多重極展開

粒子の作り出す速度場は粒子表面でのストレスΠの分布で表す事ができる。

u(r) = −\frac{1}{8πμ} ∫_{S} dSJ(r − ξ) · Π · n  (3.40)

ここで Sは粒子表面を、 ξは表面の位置を表す。 n は流体を向いた法線ベクトルである。粒子から遠く離れた場所での速度場は、粒子の中心ξ_{0}での多重極展開を用いて表す事が出来る。

u_{i}(r) = −\frac{1}{8πμ} ∑_{n=0} \{\frac{1}{n!} [ Π_{n≠0,m=1}^{n} \frac{∂}{∂ξ_{l_{m}}} ] J_{ij}(r − ξ_{0}) ∫ dS [ Π_{n≠0,m=1}^{n}(ξ_{l_{m}} − ξ_{0}_{l_{m}})] Π_{jk}n_{k}\}
= −\frac{1}{8πμ} (1 + \frac{1}{6}a^{2}∇^{2}) J_{ij}(r − ξ_{0}) · F_{j}  (3.41)

インデックス i, j, k, l_{m}はベクトル (テンソル) の成分を表し、同じものは縮役をとる。ここでは簡単の為展開は２次までをとり、またストレスΠの球面でのモーメントは０次の (粒子が流体から受ける) 力 F = ∫ dSΠ · n のみを考慮した。

相反定理^{3}

ストークス方程式に従う２つの解 (u_{1}, p_{1}, μ_{1}), (u_{2}, p_{2}, μ_{2}) を仮定する。この２つの解は流体の占める領域 V が同じ (境界 Sが同じ) だが、境界条件は違うとする。変形テンソル e とストレステンソルΠは以下のように定義される。

e = \frac{1}{2} (∇u + (∇u)^{†})  (3.42)
Π = −pI + 2μe  (3.43)

このとき対称性と非圧縮性から以下の関係が成り立つ。

μ_{2}Π_{1} : e_{2} = μ_{1}Π_{2} : e_{1}  (3.44)

ここで

Π_{1} : e_{2} = ∇ · Π_{1} · u_{2} − u_{2}∇ : Π_{1}  (3.45)

を用いて、次の関係が成り立つ。

μ_{2} (∇ · Π_{1} · u_{2} − u_{2}∇ : Π_{1}) = μ_{1} (∇ · Π_{2} · u_{1} − u_{1}∇ : Π_{2})  (3.46)

この式を流体の占める体積 V で積分すると、次の関係が示される。

μ_{2} [∫_{S} dS n · (Π_{1} · u_{2}) + ∫_{V} dr u_{2}∇ : Π_{1}] = μ_{1} [∫_{S} dS n · (Π_{2} · u_{1}) + ∫_{V} dr u_{1}∇ : Π_{2}]

(3.47)

この関係式を (ストークス方程式の) 相反定理と呼ぶ。これはこの関係を使う事でレジスタンス行列やモビリティ行列が対称行列である事が証明される事による。ファクセンの法則^{4}

3本節の内容は以下の文献に依る [24, 33]。

相反定理を用いてファクセンの法則を示す。ファクセンの法則とは、「ある速度場の中に粒子を置く場合、その粒子がどのような力を受けるか (流体に及ぼすか)」を表すものである。相反定理で考える２つの状況は、１つは粒子が並進速度 U で動いている解 (u_{1}, p_{1})、もう１つは静止している粒子と場所 yに強さが F のポイントフォースがある時の解 (u_{2}, p_{2}) とする。考える瞬間での粒子の位置は同じとする。２つの解の条件を式で書くと次のようになる。

u_{1} = U on S  (3.48)
∇ · Π_{1} = 0 in V  (3.49)
u_{2} = 0 on S  (3.50)
∇ · Π_{2} = −Fδ(x − y) in V  (3.51)

相反定理 (3.47) を用いると、

0 = U · ∫_{S} dS n · Π_{2} − F · u_{1}(y)  (3.52)

となる。ここで u_{1}(y) は１粒子が作り出す速度場で、並進速度 U を用いて次のように書ける。

u_{1}(y) = 6πμa (1 + \frac{1}{6}a^{2}∇^{2}) \frac{1}{8πμ}J(y − ξ_{0}) · U  (3.53)

ここでξ_{0}は粒子の中心の位置である。従って、状況２で粒子が流体に及ぼす力 F_{2} = − ∫_{S} dS n· Π_{2}は次のように書ける。

F_{2} = −6πμa (1 + \frac{1}{6}a^{2}∇^{2}) \frac{1}{8πμ}J(y − ξ_{0}) · F  (3.54)
= −6πμa (1 + \frac{1}{6}a^{2}∇^{2}) u′_{2}(ξ_{0})  (3.55)

ここでオセーンテンソルの性質 J(y − ξ_{0}) = J(ξ_{0} − y) を用いて、「ξ_{0}にあるポイントフォースが yに作る速度場」を「yにあるポイントフォースがξ_{0}に作る速度場」と読みかえた。従って u′_{2}(ξ_{0}) は「状況２で、粒子が無い時の粒子の中心での速度場」である。

任意の速度場はポイントフォースの重ね合わせで記述できるので、任意の速度場中に一般化出来る。また、ストークス方程式のガリレイ不変性から並進運動する粒子への一般化もできる。従って、ある速度場 u′中に並進速度 U を持つ粒子を場所ξ_{0}に置いた時、粒子が流体に及ぼす力 F は次に様に書ける。

4本節の内容は以下の文献に依る [33]。

F = 6πμa (1 + \frac{1}{6}a^{2}∇^{2}) (U − u′(ξ_{0}))  (3.56)

これをファクセンの法則と呼ぶ。証明には近似などが無いので、この関係は厳密である。ファクセンの法則を使うと、粒子の受ける力をその粒子が無い時の速度場から直接計算できる。つまり、 N粒子問題を N − 1 粒子問題に帰着できるように思われる。しかし、それは間違いである。実際、粒子が無い時の速度場 u′はある境界条件を満たしている。仮想粒子 (速度場 u′に置こうとしている粒子) の受ける力を計算する際は粒子表面での境界条件はその粒子の速度で決められている。しかし、その時生じる速度の乱れは元の速度場 u′を決める境界条件を乱しており、その補正はこの考察では採り入れられていない。 N粒子問題を考える時、ここに誤差が入る。この誤差は粒子間隔が大きい時無視できるので、希薄極限では良い近似である。

モビリティ行列

(3.41) 式の粒子の作る速度場の多重極展開と (3.56) 式のファクセンの法則を用いて、希薄極限での流体力学的相互作用を表す行列を構成する。 N個の粒子が流体中に分散している系で、１つの粒子αに着目する。粒子αが流体に及ぼす力 F^{α}はファクセンの法則を用いて次のように書ける。

F^{α} = 6πμa (1 + \frac{1}{6}a^{2}∇^{2}) (U^{α} − u′^{α}(x^{α}))  (3.57)

U^{α}は粒子αの並進速度、 u′^{α}はα以外の粒子の作り出す速度場、 x^{α}は粒子αの中心である。 u′^{α}を N − 1 粒子の多重極展開の重ね合わせで表すと次のように書ける。

u′^{α}(r) = ∑_{β≠α} \frac{1}{8πμ} (1 + \frac{1}{6}a^{2}∇^{2}) J(r − x^{β}) · F^{β}  (3.58)

(3.57) 式に代入して、粒子速度 U^{α}に対してまとめると、以下の関係を得る。

U^{α} = \frac{F^{α}}{6πμa} + ∑_{β≠α} \frac{1}{8πμ} (1 + \frac{1}{6}a^{2}∇^{2})^{2} J(x^{α} − x^{β}) · F^{β}  (3.59)

N粒子についてまとめると、モビリティ行列 M^{∞}を使って以下のように書ける。

U = M^{∞} · F  (3.60)

ここでモビリティ行列は次のような小行列から構成される。

(M^{∞})^{αα} = \frac{1}{6πμa}I  (3.61)
(M^{∞})^{αβ} = \frac{1}{8πμ} (1 + \frac{1}{6}a^{2}∇^{2})^{2} J(x^{α} − x^{β}) (α ≠ β)  (3.62)
M^{∞} = [ \begin{matrix}(M^{∞})^{11} & (M^{∞})^{12} & ⋯ & (M^{∞})^{1N} \\ (M^{∞})^{21} & (M^{∞})^{22} & ⋯ & (M^{∞})^{2N} \\ ⋮ & ⋮ & ... & ⋮ \\ (M^{∞})^{N1} & (M^{∞})^{N2} & ⋯ & (M^{∞})^{NN} \end{matrix} ]  (3.63)

3.2.3 ２体問題

前節では速度場の多重極展開とファクセンの法則を用いて希薄極限でのモビリティ行列 M^{∞}を計算した。

しかし、厳密に相互作用行列を計算する事は一般には不可能である。最も単純な場合である流体の中に２粒子のみ存在する場合の行列も、厳密に解かれたのは 80 年代の事である [25]。

ここでは簡単にその内容に触れ、シミュレーションで用いる結果をまとめておく。 2 体問題の対称性から、簡単のため 2 次元空間で議論を行う。このことにより一般性を

失うことはない。また実際にシミュレーションで使う並進速度と力の部分だけを考え、トルクや回転は無視する。粒子が無いときの流体の流れは一様流 u^{∞}とする。このとき、以下のようなレジスタンス問題が定義できる。

[ \begin{matrix}F^{1} \\ F^{2} \end{matrix} ] = μ [ \begin{matrix}A^{11} & A^{12} \\ A^{21} & A^{22} \end{matrix} ] · [ \begin{matrix}U^{1} − u^{∞} \\ U^{2} − u^{∞} \end{matrix} ]  (3.64)

ここで A^{αβ}は 2 階テンソルである。系に存在するベクトルは２つの粒子の作る境界条件を特定する中心間ベクトル r = x^{2} − x^{1}のみであることから、系に存在する２階テンソルは以下のように構成されることが分かる。

A_{ij}^{αβ} = X_{αβ}e_{i}e_{j} + Y_{αβ}(δ_{ij} − e_{i}e_{j})  (3.65)

ここで e = r/rとした。 X, Y は rに依存するスカラー関数である。従って、このスカラー関数 X(r), Y (r) を決めれば、一般のレジスタンス行列を構成する事が出来る。つまり、特別な場合の問題を解く事でこれら２つの関数を定義すればよい。実際 Xは軸対称な運動に関する部分である。中心間ベクトルが以下の場合、

e = [ \begin{matrix}0 \\ 1 \end{matrix} ]  (3.66)

レジスタンス行列は以下のようになる。

[ \begin{matrix}A_{xx}^{11} & A_{xy}^{11} & A_{xx}^{12} & A_{xy}^{12} \\ A_{yx}^{11} & A_{yy}^{11} & A_{yx}^{12} & A_{yy}^{12} \\ A_{xx}^{21} & A_{xy}^{21} & A_{xx}^{22} & A_{xy}^{22} \\ A_{yx}^{21} & A_{yy}^{21} & A_{yx}^{22} & A_{yy}^{22} \end{matrix} ] = [ \begin{matrix}Y_{11} & 0 & Y_{12} & 0 \\ 0 & X_{11} & 0 & X_{12} \\ Y_{21} & 0 & Y_{22} & 0 \\ 0 & X_{21} & 0 & X_{22} \end{matrix} ]  (3.67)

従って速度が軸対称な場合、

U^{α} = [ \begin{matrix}0 \\ U^{α} \end{matrix} ]  (3.68)

X_{αβ}のみ現れ、一方速度が非軸対称な場合、

U^{α} = [ \begin{matrix}U^{α} \\ 0 \end{matrix} ]  (3.69)

Y_{αβ}のみが現れる。

解くべき問題が限定されたとは言え、これらの関数を定義するには具体的に境界値問題を解かなければならない。 Jeffrey & Onishi (1984) は極座標系で書かれたストースク方程式の基本解を各粒子中心に対して用い、２つの極座標の変換を用いて境界条件を書き下し、スカラー関数 X(r), Y (r) を中心間距離の逆数 1/rの展開の形で与えた。

本論文で実際に用いる結果を以下にまとめておく。

X_{11} = \frac{1}{4} \frac{1}{1 − 4/s^{2}} − \frac{9}{40} ln(1 − 4/s^{2}) − \frac{3}{112}(1 − 4/s^{2}) ln(1 − 4/s^{2}) + f_{0}^{XA} − \frac{1}{4}
+ ∑_{m=2,4···}^{∞} \{(\frac{1}{4})^{m} f_{m}^{XA} − \frac{1}{4} − \frac{9}{20} \frac{1}{m} + \frac{12}{112} \frac{1}{mm_{1}} \} (\frac{2}{s})^{m}  (3.70)
X_{12} = −\frac{1}{2} \frac{1}{s(1 − 4/s^{2})} − \frac{9}{40} ln \frac{s + 2}{s − 2} − \frac{3}{112}(1 − 4/s^{2}) ln \frac{s + 2}{s − 2} − \frac{12}{112} \frac{1}{s}
− ∑_{m=1,3···}^{∞} \{(\frac{1}{4})^{m} f_{m}^{XA} − \frac{1}{4} − \frac{9}{20} \frac{1}{m} + \frac{12}{112} \frac{1}{mm_{1}} \} (\frac{2}{s})^{m}  (3.71)
Y_{11} = −\frac{1}{6} ln(1 − 4/s^{2}) + f_{0}^{Y A} + ∑_{m=2,4···}^{∞} \{(\frac{1}{4})^{m} f_{m}^{Y A} − \frac{1}{3} \frac{1}{m}\} (\frac{2}{s})^{m}  (3.72)
Y_{12} = −\frac{1}{6} ln \frac{s + 2}{s − 2} − + ∑_{m=1,3···}^{∞} \{(\frac{1}{4})^{m} f_{m}^{Y A} − \frac{1}{3} \frac{1}{m}\} (\frac{2}{s})^{m}  (3.73)

\begin{matrix}n & f_{n}^{X} & f_{n}^{Y} \\ 0 & 1 & 1 \\ 1 & 3 & 3/2 \\ 2 & 9 & 9/4 \\ 3 & 19 & 59/8 \\ 4 & 93 & 465/16 \\ 5 & 387 & 2259/32 \\ 6 & 1197 & 14745/64 \\ 7 & 5331 & 89643/128 \\ 8 & 19821 & 570017/256 \\ 9 & 76115 & 4451395/512 \\ 10 & 320173 & 33678825/1024 \\ 11 & 1178451 & 266862875/2048 \end{matrix}

Table 3.1: 二粒子問題の厳密解の展開係数。 Jeffrey & Onishi (1984)[25] に依る。

ここで、 m_{1}は次のように定義される。

m_{1} = −2δ_{m2} + (m − 2)(1 − δ_{m2})  (3.74)

また、 s = r/a は粒子半径 a でスケールされた粒子間距離 rである。

3.2.4 高濃度での相互作用行列

前の２つの節では限られた状況下での相互作用行列の具体的な形を示した。それらは希薄極限でのモビリティ行列 M^{∞}と２粒子の厳密な解に対応するレジスタンス行列である。今求めたい行列は多粒子が比較的高い粒子濃度で分布している場合に対してである。 Brady と共同研究者達は Stokesian dynamics 法と呼ばれるシミュレーション法の中で、この行列を精度良く計算する方法を開発した [28, 32, 31]。本節ではこの方法を説明し、その物理的な意味を平均沈降速度の計算への応用の過程で示す [29]。

Stokesian dynamics の方法

希薄極限でのモビリティ行列 M^{∞} の展開次数を上げる事で高濃度粒子系に適用しようとすると、行列の正定値性が破れるという問題に直面する。この問題は高濃度での抵抗係数が負の値となる事と等価である。これは自分の運動と反対の方向に摩擦力を受ける状況を意味し、全く非物理的である。

Brady らは当初この行列の正定値性が破れの困難が、２粒子問題の厳密なレジスタンス行列の重ね合わせを用いる事で回避される事を示した。しかし得られた行列の精度は満足いくものではなかった [26, 27]。

彼らは更に改良を試み、希薄極限でのモビリティ行列 M^{∞}と２粒子問題の厳密なレジスタンス行列をセルフコンシステントに結合した行列を計算した [28]。具体的には、まず２粒子問題の厳密なレジスタンス行列 R_{2B}から希薄極限でのモビリティ行列に既に含まれている寄与をあらかじめ引いた行列 R_{2B}^{lub}をつくり、

R_{2B}^{lub} = R_{2B} − (M_{2B}^{∞})^{−1}  (3.75)

全てのペアに対して重ね合わせた N粒子の lubrication^{5} 行列 R^{lub}を計算する。この行列と希薄極限でのモビリティ行列 M^{∞}を結合した行列で近似的なレジスタンス行列 R を構成する。

R = (M^{∞})^{−1} + R^{lub}  (3.76)

この行列を用いた計算はコロイド粒子系に対して数多く用いられ、その正当性は実際に示されている。

平均沈降速度

Stokesian dynamics 法は近似であり、その理論的な証明はなされていない。ここでは早川尚男氏と筆者による Stokesian dynamics 法の平均沈降速度の計算への応用を紹介し、この方法の有効性と物理的な意味を示す [29]。

平均沈降速度は、流体中に一様に分布する粒子が重力を受けて沈降する時の平均速度の事である。低レイノルズ数のもとでは、流体中にある粒子の速度 U と流体に及ぼす力 F は、流体力学的相互作用を表すモビリティ行列 M を用いて次のように書ける。

U = M · F  (3.77)

5lubrication とは粘性流体の近接物体間の相互作用の事で、「滑べり」や「潤滑」の意味を持つ。紙片を落した時床の上を滑べって行く現象がこの効果の典型的な例である。この現象では「僅かな隙間があって滑べる」部分に注意が行ってしまうが、「接触を阻む効果」 (今の場合床になかなか接触しない効果) という側面が lubrication 効果にはある事を注意しておく。

粒子の慣性が無視できる場合、 F は粒子に流体以外から働く力 (今の場合重力) に等しい。従って平均沈降速度< U >は次のように書ける。

< U >= \frac{1}{N} ∑_{α}^{N} ∑_{β}^{N} M^{αβ} · F_{g}  (3.78)

ここで F_{g}は１粒子に働く重力で３次元の定数ベクトルである。平均沈降速度はこのようにモビリティ行列の要素の平均であり、流体力学的相互作用を直接反映した物理量である。

粒子濃度が希薄な場合、２粒子相互作用の重ね合わせで平均沈降速度を近似できる。 g(r) を２体分布関数として、平均沈降速度は次のように書ける [34]。

\frac{< U >}{U_{0}} = 1 − 5ϕ + 3ϕ ∫_{2}^{∞} dr r \{g(r) − 1\}
+ϕ ∫_{2}^{∞} dr r^{2}g(r) \{x_{11} + 2y_{11} + x_{12} + 2y_{12} − 3 (1 + \frac{1}{r})\} + O(ϕ^{2}) (3.79)

ここで U_{0} = F_{g}/6πμa は１粒子の沈降速度、 ϕは粒子の体積分率で、 rは粒子半径 a で規格化してある。 x_{11}, y_{11}, x_{12}, y_{12}は２体相互作用を表すモビリティ行列 M_{2B}を構成するスカラー関数で、以下のように定義される。

M_{2B} = [ \begin{matrix}M^{11} & M^{12} \\ M^{21} & M^{22} \end{matrix} ]  (3.80)
M^{αβ}(r) = x_{αβ}(r)\frac{rr}{r^{2}} + y_{αβ}(r) (I − \frac{rr}{r^{2}})  (3.81)

ここで r は粒子半径 a で規格化された粒子間ベクトル、 x(r), y(r) は２体のモビリティ行列に対する厳密解である [25]。希薄極限での分布関数として一様分布を仮定すると、

g(r) → \{ \begin{matrix}0 & (r < 2) \\ 1 & (r ≥ 2) \end{matrix}  (ϕ ≪ 1)  (3.82)

平均沈降速度は以下のようになる。

\frac{< U >}{U_{0}} = 1 − 6.55ϕ  (3.83)

この導出は希薄極限で行われており、沈降速度の濃度展開の１次までとったものである。従って高濃度では使えない。実際ϕ ∼ 0.15 で負の値をとってしまう (図 3.1の Batchelor を参照)。

(3.79) 式の２体分布関数 g(r) に剛体粒子系のパーカス・イェビックの動径分布関数を用いることで、高濃度での沈降速度を見積もろうとする試みがある [35]。しかし結果はほとんど改善されない (図 3.1の P-Y 参照)。しかし高次の寄与である (3.79) 式の右辺第４項を無視すると高濃度でも正の沈降速度を与え、結果ももっともらしい事が分かった (図 3.1の Brady & Durlofsky 参照)。パーカス・イェビックの動径分布関数のフーリエ変換を用いてこの場合の沈降速度は解析的に求まる。

\frac{< U >}{U_{0}} = \frac{(1 − ϕ)^{3}}{(1 + 2ϕ)}  (3.84)

しかし、この結果は２つの問題点を含んでいる。１つは高次の寄与が無視できる明確な物理的理由が不明である点、もう１つは希薄極限で Batchelor の結果 (3.79) を再現していない点である。また定量的に見ると高濃度での値が実験の２倍程度大きくなっていることも問題である。

ここで Stokesian dynamics の方法を応用してみる [29]。実際にこの概念を沈降の問題に適用してみると、

F_{g} = {< (M^{∞})^{−1} > + < R^{lub} >} < U >  (3.85)

と書ける。 < M^{∞} >は多重極展開の低次の寄与のモビリティ行列の平均である。従って (3.84) 式の結果と同じものである。一方、 < R^{lub} >は (3.75) 式から次のように書ける。

< R^{lub} >= n ∫_{2}^{∞} dr g(r) [R_{2B}^{11} + R_{2B}^{12} − {(M_{2B}^{∞})^{−1}_{11} + (M_{2B}^{∞})^{−1}_{12}}]  (3.86)

ここで n = N/V は粒子数密度、 R_{2B}は二粒子問題の厳密解で、 (3.70) から (3.73) 式のレジスタンス関数 X(r), Y (r) を用いて次のように書ける。

R_{2B}^{αβ} = X_{αβ}(r) + 2Y_{αβ}(r)  (3.87)

一方希薄極限でのモビリティ行列 M_{2B}^{∞}は (3.61) 式と (3.62) 式で定義される小マトリックスから構成できる。ここではストークス因子 6πμa で無次元化して書く。

M_{2B}^{∞} = [ \begin{matrix}I & a^{∞} \\ a^{∞} & I \end{matrix} ]  (3.88)
a^{∞} = [\frac{3}{2} \frac{a}{r} − (\frac{a}{r})^{3}] \frac{rr}{r^{2}} + [\frac{3}{4} \frac{a}{r} + \frac{1}{2} (\frac{a}{r})^{3}] (I − \frac{rr}{r^{2}})  (3.89)
x^{∞}(r) = \frac{3}{2} \frac{a}{r} − (\frac{a}{r})^{3} , y^{∞}(r) = \frac{3}{4} \frac{a}{r} + \frac{1}{2} (\frac{a}{r})^{3}  (3.90)

希薄極限でのモビリティ行列の逆行列 (M_{2B}^{∞})^{−1}は、次のように求められる。ここで議論する沈降問題の場合、粒子に作用する力は全て同じである。 F_{1} = F_{2} = F またモビリティ行列 M2B^{∞}の対称性から U_{1} = U_{2} = U となることが分かる。従ってこの逆行列を求める際には、次の形の 1 粒子問題で考えて良い。

U = M_{1B}^{∞} · F  (3.91)

ここで M_{1B}^{∞} = (I + a^{∞}) で、 2 × 2 行列である。（２粒子問題は粒子の配置の対称性から、一般性を失うことなく２次元で考えられる。） M_{1B}^{∞}の逆行列は、明らかに M_{2B}^{∞}と次の関係がある。

(M_{1B}^{∞})^{−1} = (M_{2B}^{∞})^{−1}_{11} + (M_{2B}^{∞})^{−1}_{12}  (3.92)

この行列も一般に rの関数で書ける。

(M_{1B}^{∞})^{−1} = X^{∞}(r)e_{i}e_{j} + Y^{∞}(r) (δ_{ij} − e_{i}e_{j})  (3.93)

と書ける。更に r = (0, 0, r) の場合を考えると、各行列の要素が次のようになる。

M_{1B}^{∞} = [ \begin{matrix}1 + y^{∞}(r) & 0 \\ 0 & 1 + x^{∞}(r) \end{matrix} ]  (3.94)
(M_{1B}^{∞})^{−1} = [ \begin{matrix}Y^{∞}(r) & 0 \\ 0 & X^{∞}(r) \end{matrix} ]  (3.95)

従って、 (M_{1B}^{∞})^{−1}は求まった。

X(r) = \frac{1}{1 + x^{∞}(r)}, Y (r) = \frac{1}{1 + y^{∞}(r)}  (3.96)

従って、 (3.86) 式は次のように書ける。

< R^{lub} > = n ∫ dr g(r) (R_{2B}^{11} + R_{2B}^{12} − (M_{1B}^{∞})^{−1})  (3.97)
= ϕ ∫ dr r^{2}g(r) \{X_{11} + X_{12} + 2Y_{11} + 2Y_{12} − \frac{6r^{3}(−2 + 5r^{2} + 4r^{3})}{(−2 + 3r^{2} + 2r^{3})(2 + 3r^{2} + 4r^{3})}\}

g(r) にパーカス・イェビックの分布関数を用いて数値積分を行うと、積分はほとんどϕに依存しない事が分かる (図 3.1参照)。希薄な場合の一様分布 (3.82) 式を使うと以下のような結果を得る。

< R^{lub} >≃ 1.49ϕ  (3.98)

以上から、求めようとする平均沈降速度< U > /U_{0}は次のようになる。

\frac{< U >}{U_{0}} = < M >
= \frac{1}{< (M^{∞})^{−1} > + < R^{lub} >}  (3.99)
= \frac{< M^{∞} >}{1+ < R^{lub} >< M^{∞} >}  (3.100)

< M^{∞} >に (3.84) 式の結果を用い、 < R^{lub} >に (3.98) 式を用いると、以下の結果を得る (図 3.1参照)。

\frac{< U >}{U_{0}} = < M >
∼ \frac{(1 − ϕ)^{3}}{1 + 2ϕ + 1.49ϕ(1 − ϕ)^{3}}  (3.101)

この結果は、希薄極限で

\frac{< U >}{U_{0}} = 1 − 6.49ϕ  (3.102)

となる。この結果 (3.101) は、希薄極限での (3.79) 式への一致は完全ではないが、 (3.84) 式

Figure 3.1: 平均沈降速度の理論と実験結果。 Ours は (3.101) 式、 Ours with P-Y はパーカス・イェビックの分布関数を用いた結果、 Batchelor は (3.83) 式、 P-Y は (3.79) 式にパーカス・イェビックの分布関数を用いたもの、 Brady & Durlofsky は (3.84) 式、その他のマークは各々以下の実験結果である。 exp 1 は Bacri et al.(1986)[36]、 exp 2 は Buscall et al.(1982)[37]、 exp 3 は de Kruif et al.(1987)[38]、 exp 4 は Paulin et al.(1990)[39]、 exp 5 は Xue et al.(1990)[40]。

よりも改善されている。また、実験値との一致も (3.84) 式よりも改善されている。このことは、 Stokesian dynamics の方法の正当性の間接的な証明でもある。

以上平均沈降速度の理論を紹介して来たが、図 3.1に実験結果と合わせてプロットする。この図からも分かる通り、 Batchelor の希薄極限での理論 (3.79) とその理論にパーカス・イェビックの分布関数を用いた結果は薄い濃度で負の値をとってしまう。 Brady らの近似 (3.84) は高濃度でも定性的には良い振舞いを示すが、中間の濃度での実験からのずれは明らかである。我々の結果 (3.101) は Brady らの結果の不一致を明らかに改善しており、実験との一致も良いことが分かる。物理的意味

前節で示したように Stokesian dynamics の方法は、高濃度コロイド系の挙動を正しく再現し、平均沈降速度の高濃度での振舞いも実験と良く合う近似である。

Stokesian dynamics 法は「本来多体問題である流体力学的相互作用を表す行列を２体効果の重ね合わせで近似できる方法はないか」という問への１つの回答である。

今の場合重ね合わせの方法には２種類ある。１つは速度場の重ね合わせ、つまりモビリティ形式での重ね合わせで、もう１つは力の重ね合わせ、つまりレジスタンス形式での重ね合わせである。モビリティ形式での重ね合わせは 3.2.2節で示した通り希薄極限では正しいが、粒子間隔が小さくなると速度場の境界条件を合わせる為に高次の展開が必要となり、２体効果の重ね合わせで近似できなくなる。

２体効果の重ね合わせは、前節の沈降速度の理論で見た通り主要な寄与は体積分率ϕに比例する補正を与える。従ってモビリティ形式での重ね合わせを行えば、高濃度では必然的に負の沈降速度を与えてしまう。逆の問題であるレジスタンス形式で重ね合わせは、粒子濃度の増加に伴って抵抗が増加して行くという自然な傾向を表し、非物理的な負の抵抗 (負の沈降速度) を与えるような事は起こらない。

このモビリティ形式での重ね合わせとレジスタンス形式で重ね合わせは、電流回路の抵抗の並列結合と直列結合の関係に似ている。

U^{α} = ∑_{β} M^{αβ} · F^{β} モビリティ形式 (並列結合) I = ∑_{i} \frac{V}{R_{i}}  (3.103)
F^{α} = ∑_{β} R^{αβ} · U^{β} レジスタンス形式 (直列結合) V = ∑_{i} R_{i}I  (3.104)

このように Stokesian dynamics 法は、希薄極限で正当化される寄与 (遠方場の寄与) は速度場の重ね合わせで表現し、この近似が正当性を失う近接場の寄与 (lubrication 効果) は、物理的に妥当な結果を導く力の重ね合わせで表現する事で全体として破綻の無い近似法であると言える。

3.2.5 周期境界条件

Stokesian dynamics 法では M^{∞}が長距離的であるので、周期境界条件を用いる時は取扱に注意が必要である。実際のシミュレーションで周期境界条件を持つ系の流体力学的な相互作用を計算する時、 Ewald の方法を用いると効率的に行う事が出来る [30]。 3.2.2節の結果から、 M^{∞}によって以下のモビリティ問題が定義できる。

U^{α} = \frac{1}{6πμa} ∑_{γ} ∑_{β=1}^{N} M̂^{∞}(x^{α} − x^{β} + r_{γ}) · F^{β}  (3.105)

ここで、 M̂^{∞} = 6πμaM^{∞} はストークス因子で無次元化した希薄極限でのモビリティ行列で、 r_{γ}は周期セル γを表す格子ベクトルである。 M̂^{∞}を具体的に書き下すと、以下のようになる。

M̂^{∞}(r) = \frac{6πμa}{8πμ} (1 + \frac{a^{2}}{6} ∇^{2})^{2} J(r)
= \frac{3}{4} \frac{a}{r} (I + \frac{rr}{r^{2}}) + \frac{1}{2} \frac{a^{3}}{r^{3}} (I − 3\frac{rr}{r^{2}})  (3.106)

これは微分を使って表すと、以下のように書ける。

M̂^{∞} = (\frac{3}{4}a + \frac{1}{4}a^{3}∇^{2}) (∇^{2}I − ∇∇) r  (3.107)

格子和の収束を早める方法として、 Ewald の方法を使う。 M̂_{ij}^{∞}(r) を、次のように誤差関数を用いて２つの部分に分ける。

M̂^{∞}(r) = M^{(1)}(r) + M^{(2)}(r)  (3.108)
M^{(1)}(r) = (\frac{3}{4}a + \frac{1}{4}a^{3}∇^{2}) (∇^{2}I − ∇∇) \{r erfc(ξ r)\}  (3.109)
M^{(2)}(r) = (\frac{3}{4}a + \frac{1}{4}a^{3}∇^{2}) (∇^{2}I − ∇∇) \{r erf(ξ r)\}  (3.110)

ここで、誤差関数は次のように定義する。

erfc(x) = \frac{2}{\sqrt{π}} ∫_{x}^{∞} exp(−z^{2})dz  (3.111)
erf(x) = \frac{2}{\sqrt{π}} ∫_{0}^{x} exp(−z^{2})dz  (3.112)
erf(x) + erfc(x) = 1  (3.113)

従って、次の関係が得られる。

6πμaU^{α} = F^{α} + ∑_{γ} \begin{matrix}N \\ ∑ & ′ \\ β=1 \end{matrix} M^{(1)}(x^{α} − x^{β} + r_{γ}) · F^{β}
+ ∑_{γ} ∑_{β=1}^{N} M^{(2)}(x^{α} − x^{β} + r_{γ}) · F^{β} − M^{(2)}(r = 0) · F^{α}  (3.114)

この式の始めの格子和∑′は、セルフパート (同一セル内の同一粒子の寄与) を除くことを表し、二番目の格子和は、全てとるものとする。右辺第一項が、セルフパートの寄与を表している。最後の項は、二番目の格子和に含まれるセルフパートの寄与を取り除いている。始めの格子和は、誤差関数 erfc(ξr) の為、実空間で早く収束する。一方、二番目の格子和は、逆空間で早く収束する形を持っている。分離パラメータξは、この二つの格子和の両方が早く収束するように選ぶ。

M^{(1)}は直接計算する事が出来て、以下のように書ける。

M^{(1)}(r) = erfc(ξr) \{I (\frac{3}{4} \frac{a}{r} + \frac{1}{2} \frac{a^{3}}{r^{3}} ) + \frac{rr}{r^{2}} (\frac{3}{4} \frac{a}{r} − \frac{3}{2} \frac{a^{3}}{r^{3}} )\}

+\frac{e^{−ξ^{2}r^{2}}}{\sqrt{π}} [I (4a^{3}ξ^{7}r^{4} + 3aξ^{3}r^{2} − 20a^{3}ξ^{5}r^{2} − \frac{9}{2}aξ + 14a^{3}ξ^{3} + a^{3}ξ \frac{1}{r^{2}} )

+\frac{rr}{r^{2}} (−4a^{3}ξ^{7}r^{4} − 3aξ^{3}r^{2} + 16a^{3}ξ^{5}r^{2} \frac{3}{2}aξ − 2a^{3}ξ^{3} − 3a^{3}ξ \frac{1}{r^{2}} )] (3.115)

二番目の格子和についてみる。これは、波数空間で早く収束する。格子ベクトル r_{γ}は、整数 n_{1}, n_{2}, n_{3}とセルサイズ L_{1}, L_{2}, L_{3}を用いて以下のように書ける。

r_{γ} = [ \begin{matrix}n_{1}L_{1} \\ n_{2}L_{2} \\ n_{3}L_{3} \end{matrix} ]  (3.116)

ここで波数ベクトル k とセルの体積 V を以下のように定義する。

k_{λ} = ( \frac{2πm_{1}}{L_{1}} , \frac{2πm_{2}}{L_{2}} , \frac{2πm_{3}}{L_{3}} )

V = L_{1}L_{2}L_{3}

波数空間での格子和は、ポワソンの公式を用いて以下のように変形できる。

∑_{γ} ∑_{β=1}^{N} M^{(2)}(x^{α} − x^{β} + r_{γ}) · F^{β} = \frac{1}{V} ∑_{λ} ∑_{β=1}^{N} e^{−ik_{λ}·(x^{α}−x^{β})}M^{(2)}(k_{λ}) · F^{β}
= \frac{1}{V} ∑_{λ} ∑_{β=1}^{N} cos(k_{λ} · (x^{α} − x^{β}))M^{(2)}(k_{λ}) · F^{β}(3.117)

ここで、 M^{(2)}(k) は M^{(2)}(r) のフーリエ変換である。

M^{(2)}(k) = ∫ dr M^{(2)}(r) exp(ik · r)
= 6πa (I − \frac{kk}{k^{2}}) \frac{1}{k^{2}} (1 − \frac{1}{3}a^{2}k^{2}) (1 + \frac{k^{2}}{4ξ^{2}} + \frac{k^{4}}{8ξ^{4}} ) e^{− \frac{k^{2}}{4ξ^{2}}}  (3.118)

最後の項 M^{(2)}(r = 0) は M^{(2)}(k) を用いて計算できる。

M^{(2)}(r = 0) = \frac{1}{(2π)^{3}} ∫ dk M^{(2)}(k)
= \frac{I}{\sqrt{π}} (6aξ − \frac{40}{3}a^{3}ξ^{3})  (3.119)

最終的に周期境界条件のもとでの希薄極限のモビリティ行列は以下のように書ける。

6πμaU^{α} = F^{α} + ∑_{γ} \begin{matrix}N \\ ∑ & ′ \\ β=1 \end{matrix} M^{(1)}(x^{α} − x^{β} + r_{γ}) · F^{β}
+\frac{1}{V} ∑_{λ} ∑_{β=1}^{N} cos(k_{λ} · (x^{α} − x^{β}))M^{(2)}(k_{λ}) · F^{β} − M^{(2)}(r = 0) · F^{α}(3.120)

ここで M^{(1)}(r), M^{(2)}(k), M^{(2)}(r) はそれぞれ式 (3.115),(3.118),(3.119) により与えられる。

(3.120) 式の形から、周期境界条件の下でもセル内の N粒子に対してモビリティ行列が構成できることが分かる。従って 3.2.4節で示した Stokesian dynamics 法は M^{∞}を (3.120) 式に置き換える事で周期境界条件でも使う事が出来る [32, 31]。

3.3 粒子慣性の導入

本節では、前の 3.2節で示した流動層で重要となる粒子と流体の相互作用を用いて粒子の運動方程式を構成する際の、粒子慣性の取扱を議論する。

前節の議論から相互作用が既に分かった訳であり、運動方程式を書く事は一見自明の事と思われる。しかし粒子の慣性を有効に導入するように運動方程式を書くことは、想像される程単純ではない。本節では、まず粒子慣性が無い時の運動方程式を考察し、この式に単純に粒子慣性を導入する事が出来ない事を示す。次に、繰り込まれた粒子慣性の導入によりこの困難が取り除かれることを示し、最後にこの「繰り込まれた粒子慣性」の物理的意味を議論する。

3.3.1 粒子慣性の無い場合

3.2節の結果から粒子が流体から受ける力 F_{f} はストークス近似のもとで以下のように書ける。

F_{f} = −R · (U − u^{∞})  (3.121)

ここで R はレジスタンス行列で、流体による相互作用を表し、粒子の配置に依存する量である。 U,u^{∞}は粒子速度、流体の流入速度を表す。

最初に粒子の慣性が無視できる場合を考える。コロイド粒子系などのような系がこの場合に対応する。この時粒子に働く力は流体からの力 F_{f}と重力 F_{g}のみを考えると、運動方程式は次のように書ける。

0 = −R · (U − u^{∞}) + F_{g}  (3.122)

この方程式は未知数である粒子の速度 U について解く事が出来る。

U − u^{∞} = R^{−1} · F_{g}  (3.123)

従って初期条件として粒子の配置が与えられれば、その瞬間の粒子の速度が (3.123) 式を用いて解ける。この速度を積分する事で微小時間後の粒子の位置が計算できる。この操作を逐次的に行うことで粒子の軌道が得られる。見掛け上 (3.123) 式は定常のようであるが、得られる粒子速度によって粒子の配置 (境界条件) が時間と共に変化する。

粒子間衝突は今の場合そもそも導入する必要は無い。実際 (3.123) 式で決まる速度を粒子が持つ場合は、レジスタンス行列の持つ (或は低レイノルズ数流体力学が持つ) lubrication 効果の為に、粒子間隔が小さくなるとその粒子間の相対速度も小さくなり、原理的に粒子は接触する事は無い。これは粒子の中心間距離 rが小さくなると、抵抗係数が 1/(r − 2a) の形で発散する事に対応する。

このことは必ずしも粒子配置が一様になる事を意味しない。それは抵抗が常に反発力ではないからである。実際に粒子に働く力は抵抗係数と粒子速度の積であるため、 (他の粒子の効果など) 何かのはずみである粒子ペアの粒子間隔が一旦小さくなると、そのペアは近付きにくくなると同時に離れにくくなる。この結果、粒子慣性が無視できる系では粒子はクラスターを組み、一様状態は得られない。このことは、流体力学的相互作用が粒子間の相対速度を減少させる摩擦として作用している事、粒子濃度が高い程抵抗が大きい事を考えると自然である。

粒子慣性が無視できるような系は通常シアー流れの中での挙動など、流体力学的相互作用が作り出すクラスターを壊すメカニズムが存在する状況での振舞が主な研究対象となっている。

3.3.2 粒子慣性の導入の困難

粒子慣性を (3.122) 式に単純に導入すると、運動方程式は次のようになる。

m\frac{d}{dt}U = −R · (U − u^{∞}) + F_{g}  (3.124)

しかしこの粒子慣性の効果 mdU/dt は、 R に含まれる lubrication 効果の為にほとんど効かない。つまり、慣性の無視できる時の状況と同じように粒子は衝突する事は無い。実際に、ある接近した粒子の挙動を取り出してみる。簡単の為粒子の動径方向のみ考えると、粒子は次の方程式に従う。

m\frac{dU}{dt} = −\frac{1}{r − 2a}U  (3.125)

ここで粒子速度を U = dr/dt とした。このとき、粒子の負の方向への初速 (接近する方向) がどんなに大きくても粒子は接触点 r = 2a に至ることはできない。つまり粒子の初期の運動エネルギーがいかに大きくても、抵抗力によって散逸するエネルギーが接触点 r = 2a に至る過程で無限になり得るからである。

3.3.3 lubrication の繰り込み

現実の流動層では粒子間衝突は重要な寄与をもつ。しかし前節で見たように、慣性の無視できる運動方程式 (3.122) に普通に慣性を導入しても、 lubrication 効果によって粒子間衝突を引き起こすような粒子の相対速度が全て緩和してしまう。つまり、粒子の運動方程式 (3.1) に lubrication をもつ流体力を用いると慣性の効果が隠されてしまう。これは lubrication のもつ特異性による。

この矛盾を解決する為に、 lubrication の特異性を打ち消すように繰り込まれた粒子慣性を (3.122) 式に導入する。つまり、運動方程式を次のように書く。

m′R̂ · \frac{d}{dt}U = −R · (U − u^{∞}) + F_{g}  (3.126)

ここでR̂ = R/6πμa は無次元化されたレジスタンス行列、 m′は繰り込まれた質量である。この方程式に左からR̂^{−1}を掛け、次式で与えられる速度 Vを用いると、

V − u^{∞} = R^{−1} · F_{g}  (3.127)

次のように変形できる。

m′ \frac{d}{dt}U = −6πμa(U − V)  (3.128)

この式は粒子速度 U が時間 m′/6πμa で Vに向かって緩和する事を意味する。このことから Vを終端速度と呼ぶ。 Vの定義式 (3.127) は (3.123) 式と同じである。つまり粒子慣性が無い時の粒子速度に一致する。よって粒子速度 U の Vからのずれは時間 m′/6πμa の間残る。このずれが粒子衝突を引き起こす。この結果、 (3.126) 式は現実的な振舞を示すモデルとなっている事が分かった [41]。

3.3.4 無次元化

後の議論との関係を明確にし、実際のシミュレーションを行う時の便宜の為に、ここで運動方程式を無次元化しておく。１粒子の沈降速度 U_{0}と粒子半径 a で無次元化すると、 (3.128) 式は次のように書ける。

St′ \frac{d}{dt̂}Û = −(Û − V̂)  (3.129)

ここで無次元化された量はˆを用いて表す。 St′は (繰り込まれた) ストークス数で以下のように定義した。

St′ = \frac{m′U_{0}}{6πμa^{2}}  (3.130)

ストークス数は粒子の慣性と流体の粘性の比を表す無次元量である。無次元化された終端速度V̂は次のように決まる。

V̂ − û^{∞} = −R̂^{−1} · E_{z}  (3.131)

ここで E_{z}は (3.2) 式の形に拡張した z方向の単位ベクトルである。後の議論の為に普通に慣性を導入した運動方程式 (3.124) も同様に無次元化しておく。

StR̂^{−1} · \frac{d}{dt̂}Û = −(Û − V̂)  (3.132)

ここで St はストークス数である。

St = \frac{mU_{0}}{6πμa^{2}}  (3.133)

ここでは (3.129) 式との関係を見やすくする為に、両辺に R^{−1}を左から掛けて整理した形で書いた。

本論文では以後無次元化された方程式のみ扱うので、簡単の為、無次元化を表すˆは省略する。

3.3.5 繰り込みの物理的意味

3.3.2節の議論のように、粒子の慣性があっても流体力学的相互作用を考慮すると粒子間衝突が生じ得ないのではないか、との疑問はこの分野の大きな問題である。

最近この疑問に直接取り組んだ研究が Koch らによって行われた [42]。彼らの主張は「粒子が接近した時、間隔が周りの流体分子の平均自由行程程度になったとき流体力学的な記述が破綻し、従って lubrication は効かなくなる」というものだ。彼らはこの考えを具体的に用いて、有限な粒子慣性を持つ粒子系のシアー流れの中でのシミュレーションを行った [43]。彼らのシミュレーションは普通に粒子慣性を導入した運動方程式 (3.132) を直接計算したものである。しかし上述の通り、この式の直接計算では粒子の衝突が lubrication の効果によって隠されてしまう。ここで彼らは lubrication の効く粒子間の隙間にカットオフ (その値は粒子半径で規格化して 0.01) を導入し、隙間がそれ以下の場合はカットオフでの値を用いて流体力学的相互作用を計算した^{6}。興味深い事に、彼らのシミュレーション結果から St_{s}/ < R_{c} >というパラメータが系の振舞いを特徴付ける事が示されている。ここで St_{s} = mγ/6πμa はシアー流れ中でのストークス数で、 γはシアー率である。 St_{s}の物理的意味は流動層中でのストークス数 (3.133) と同じで、粒子の慣性と流体の粘性の大きさの比である。また< R_{c} >は固定粒子中でトレーサー粒子が受ける抵抗率で、シミュレーションで導入されたカットオフに依存する量である。この量はレジスタンス行列の平均と解釈できる。

カットオフを導入するモデルと、繰り込まれた慣性を用いたモデルの比較を行う。重複するが２つの方程式をここに書く。

StR_{c}^{−1} · \frac{d}{dt}U = −(U − V)  (3.134)
St′ \frac{d}{dt}U = −(U − V)  (3.135)

ここで R_{c}はカットオフを導入したレジスタンス行列を意味する。この２つの式の違いは左辺の慣性項にある。最も大きな違いはスカラーとテンソルの違いである。具体的には、 U の Vへの緩和がカットオフを導入するモデルの場合、 R_{c}^{−1}を通して多体効果 (多分散性) が入って来るが、繰り込まれた慣性を用いたモデルでは粒子依存性は無い。これは繰り込まれたストークス数 St′が R_{c}^{−1}のある種の平均値を用いて以下のように近似されたものと解釈できる。

6(3.132) 式は直接差分化した場合数値的に不安定である。彼らは少し変わった差分化を行ってシミュレーションしている。つまり速度の発展は中央差分を用い、その為に粒子の位置とは半ステップずらして計算している。

St′ ∼ St < R_{c}^{−1} >∼ StR_{c}^{−1}  (3.136)

従って、カットオフを導入するモデルのシミュレーションを特徴付ける St/ < R_{c} >が繰り込まれた慣性を用いたモデルの St′に対応する事を意味することが分かる。

以上の考察から、本研究での繰り込まれた慣性を用いたモデル化は、 lubrication にカットオフを導入するモデルの近似に対応する事が分かる。しかし繰り込まれた慣性を用いたモデルには長所がある。それは、カットオフを導入するモデルが当然カットオフの値に依存するのに対し、粒子慣性を繰り込んだモデルではカットオフというパラメータが無い点にある。 lubrication の発散を抑制するメカニズムは Koch らが指摘した流体の連続的な記述の破綻のみではない。例えばより現実的に考えられるものとして、粒子表面の乱れ (細かな凹凸) がある。この乱れを考えた場合、粒子の衝突は平均的な粒子表面間の隙間がゼロになる前に生じてしまう為、 lubrication の発散は抑制される。つまり lubrication の発散を隠すメカニズムはいろいろなものが考えられる。繰り込まれた慣性を用いたモデル化は lubrication の発散を隠すメカニズムの詳細に依らない定式化になっている。

しかしカットオフを導入するモデルに存在する慣性項の多体効果 (多分散性) が、繰り込まれた慣性を用いたモデルには無いことは留意するべきであろう^{7}。

3.3.6 固定粒子

これまでの議論では空間に固定された粒子が無い場合を議論して来たが、モデルの定式化の最後に系の中に固定粒子を導入した場合の運動方程式を議論する。

レジスタンス行列は粒子の配置にのみ依存するので、自由粒子と固定粒子の区別はない。まず、自由粒子と固定粒子を含めた系のレジスタンス行列が計算できているとする。

[ \begin{matrix}F_{m} \\ F_{f} \end{matrix} ] = [ \begin{matrix}R_{mm} & R_{mf} \\ R_{fm} & R_{ff} \end{matrix} ] · [ \begin{matrix}U_{m} − u^{∞} \\ −u^{∞} \end{matrix} ]  (3.137)

ここで、自由粒子 (mobile particle) の部分は添字 m を、固定粒子 (fixed particle) の部分は添字 fを用いて表す。また固定粒子の速度 U_{f}はゼロとした。自由粒子の運動方程式で必要となるのは、 F_{m}である。

7この多分散性の欠如は、本研究で示す繰り込まれた慣性を用いたモデルの抱える問題の１つである「一様流動層が見られていない」原因ではない事は 5.2節で議論する。

F_{m} = R_{mm} · (U_{m} − u^{∞}) − R_{mf} · u^{∞}  (3.138)

この場合修正を受けるのは終端速度 Vだけである。実際、固定粒子がある場合の終端速度は次のように書ける。

V − u^{∞} = R_{mm}^{−1} · (−E_{z} + R_{mf} · u^{∞})  (3.139)

この終端速度を使うと運動方程式 (3.129) は全く修正の必要は無い。

計算結果と解析

4.1 シミュレーションの設定

4.1.1 パラメータ空間

実際のシミュレーションを行う粒子の支配方程式 (3.129) に含まれるパラメータにはストークス数 St ^{1}と流入速度 u^{∞}がある。この他に、粒子数 Nや周期境界のセルの大きさ L_{x}, L_{y}, L_{z}がある。

本研究では粒子の運動を (x, z) の鉛直面内に限定した２次元的な状況でのシミュレーションを行った。ここで L_{y} = 2 とし、 y方向にも周期境界条件は適用する。従って流体力学的相互作用は３次元的なものである。 (これは主に計算機の能力の限界による。 )

粒子数 N = 266 で周期セルの大きさは (L_{x}, L_{y}, L_{z}) = (34, 2, 100) の場合を中心にして、上記のコントロールパラメータ u^{∞}, St を系統的に変化させてシミュレーションを行い、系の振舞いを定量的に解析する。その他のセルも用いたが、本研究では上の状況を補う為のものである。実際に行ったシミュレーションの粒子数と周期セルについては以下の表にまとめておく。

1実際には「繰り込まれたストークス数 St′」であるが、本章では単に「ストークス数 St」と呼ぶ事にする。

粒子数周期セルサイズ
N  (L_{x}, L_{y}, L_{z})
260  (14, 2, 250)
◎ 266  (34, 2, 100)
266  (34, 2, 75)
266  (34, 2, 50)
276  (70, 2, 70)

Figure 4.1: パラメータ空間 (u^{∞}, St)。粒子数 N = 266、周期セル (L_{x}, L_{y}, L_{z}) = (34, 2, 100) の結果に対する。図中の印は dead が「静止状態」、 bubbling が「気泡状態」、 channelling が「チャンネル状態」、また線は 4.3.1節で定義する流動化開始速度である。

粒子数 N = 266 で周期セルの大きさが (L_{x}, L_{y}, L_{z}) = (34, 2, 100) の場合のシミュレーションを、図 4.1にまとめておく。縦軸のストークス数 St に対しては対数目盛をとっている。ここでは、実際にシミュレーションで得られた状態によって３種類に分類し印を付した。また同時に 4.3.1節での議論で見積もる流動化開始速度 u_{c}も合わせてプロットした。

ここで示した状態の分類は、見ための分類である。定量的な分類は以下の節で行うが、ここではその状態の典型的な粒子の配置を紹介する。流速が小さい場合は粒子は底に堆積したまま完全に運動を停止してしまう (図 4.2)。この状態を「静止状態」と呼ぶ。

Figure 4.2: 静止状態の例。 St = 10.0, u^{∞} = 0.05 での計算結果。時間は左から右に 50 step 毎に示した。 1 step は無次元化された時間で 1 である。

流速が大きくなると流動化開始速度 u_{c} (定義は 4.3.1節で与える。 ) で流動化する。シミュレーション結果からストークス数の大きさによって流動形態が２種類観測された。ある程度ストークス数が大きい時は、底の部分で形成された気泡が粒子層の中を上昇して吹き抜ける挙動が観測された (図 4.3)。この状態を「気泡状態」と呼ぶ。

一方ストークス数が小さい場合、流体が流れる流路が形成される「チャンネル状態」が得られた (図 4.4)。

しかし二つの図 4.3と 4.4を見ても分かる通り、この違いはかなり微妙なものである。以下の節ではこの点も含めてシミュレーションで得られた結果の定量的な解析を試みる。 Figure 4.3: 気泡状態の例。 St = 10.0, u^{∞} = 0.15 での計算結果。時間は左から右に 50 step 毎に示した。 1 step は無次元化された時間で 1 である。

Figure 4.4: チャンネル状態の例。 St = 0.5, u^{∞} = 0.15 での計算結果。時間は左から右に 20 step 毎に示した。 1 step は無次元化された時間で 1 である。

4.1.2 シミュレーションのサイズについて

今回主に計算を行った設定は、 N = 266, (L_{x}, L_{y}, L_{z}) = (34, 2, 100) である。この設定を選んだ理由は大きく２つある。１つはセルの横幅の問題、もう１つは層高の問題である。

横幅に関しては、狭い場合系の挙動が１次元的 (鉛直成分のみ) になってしまう。２次元的振舞と１次元的振舞ではその挙動が大きく変わってしまう事が分かる。この違いを流動状態で言うと、２次元的挙動で特徴的な状態は気泡状態であり、１次元的挙動はスラグ状態である。現実の系では流速を上げる事で気泡状態からスラグ状態に転移する事が知られている。本研究が流動化現象の解析を目指している事からも、２次元的挙動が観測できる程の横幅は必要である。

層高に関しては、もし気泡の大きさと同程度の層高であれば、気泡状態とチャンネル状態を識別する事が困難になる。多様な流動状態を考察する為にもある程度の高さは必要である。

もちろん粒子数には実際の計算のコストから決まる上限がある為、セルの横幅や層高を無闇に大きく取る事は不可能である。現在のところ粒子数は 300 を少し下回るくらいが限界である。従って、本研究での解析の中心に上述の設定を据えた。

本研究でのモデル化では流体力学的に容器の効果を導入する代わりに周期境界条件を用いている。現実の系との対応を考えた場合、特に鉛直方向の周期性が奇異に感じられる^{2}。ここで鉛直方向のセルサイズ L_{z}のみを変化させた場合のシミュレーションを２つ行った。ここでは基本設定での値 L_{z} = 100 に対して、 L_{z} = 75, 50 について計算を行った。図 4.5 に３つの場合の粒子の挙動を示す。いずれも１つの気泡が底の辺りから吹き抜けるまでの

2最近、１方向は周期境界条件を取らないで２次元的に周期境界条件を課した場合の流体力学的相互作用の計算が行われている事を、御本人より教えて頂きました [44]。従って将来的には鉛直方向の周期境界を取らない状況での解析ができるだろう。

Figure 4.5: 系の振舞の鉛直方向のセルサイズ L_{z}依存性。パラメータは上から L_{z} = 100, 75, 50、その他は N = 266, St = 10, u^{∞} = 0.2 での計算結果である。時間は左から右に、シミュレーション開始後 300step から 20 step 毎に表示した。 1 step は無次元化された時間で 1 である。この結果から、気泡の大きさ、形、上昇速度もほとんど同じ事が分かる。

挙動である。この図から、セルの高さを除いて、気泡の大きさ、形、上昇速度などはほとんど変化が無い様に見える。定量的な物理量に関する比較として、図 4.6に粒子の運動エネ

Figure 4.6: 系の振舞の鉛直方向のセルサイズ L_{z}依存性。ここでは運動エネルギーの時間変化を示した。パラメータは L_{z} = 100, 75, 50、その他は N = 266, St = 10, u^{∞} = 0.2 での計算結果である。この結果から、変動幅、変動周期もほとんど同じ事が分かる。

ルギーの時間変化を示した。この図から分かるように、エネルギーの変動幅、変動周期もほとんど同じ事が分かる。このように想像されるよりも系の挙動の L_{z}依存性は無い事が分かった。今回の L_{z} = 100 という設定は大きな流速 u^{∞} = 0.8 程度でも、流体に吹き上げられた粒子が天井に届くか届かないか、と言う事を目安に決めた。

4.2 動的挙動

4.2.1 振動的挙動と気泡

粒子は底に配置された固定粒子により支えられているため、粒子の重心は平均的には運動していない。しかし流動化した系は粒子間の相対速度は持っている。１粒子当たりの粒子の運動エネルギー Eを考える。

E(t) = \frac{1}{N} ∑_{α=1}^{N} U^{α} · U^{α}  (4.1)

この物理量 Eは流動状態と静止状態を区別する為の良い観測量である。実際にこの量の時間変化を計算して見ると、気泡の生成に伴ってきれいなピークが観測される事が分かった。 Figure 4.7: 運動エネルギーの時間変化。 N = 266, St = 20, u^{∞} = 0.3 での結果に対する。気泡の存在に対応して振動的な挙動が見られる。

典型的な気泡流動状態での E(t) の振舞を図 4.7に示す。ここで通常用いられる速度の分散 < |U− < U > |^{2} >を使わない理由を述べる。本研究の場合固定粒子を導入している為この粒子との衝突では運動量が保存しない。従って分散の時間変動には衝突に伴う不連続点が数多く現れる。一方 Eは図 4.7に見られるように連続的に変動する為、後の解析に便利である。この連続性は、結果を解釈する時も良い見通しを与える。

流動層の状態を特徴付ける物理量として、粒子の重心 H も大切な量である。

H(t) = \frac{1}{N} ∑_{α=1}^{N} z^{α}  (4.2)

ここで z^{α} は粒子 α の位置の z成分である。固定粒子の位置を z = 0 とする。図 4.7に示した同じシミュレーションに対して重心の時間変化を計算した結果を図 4.8に示す。当然、

Figure 4.8: 重心の時間変化。 N = 266, St = 20, u^{∞} = 0.3 での結果に対する。気泡の存在に対応する振動的挙動が見られる。

系の中に気泡が出来ると Hは大きくなり、気泡が抜けると小さくなる。

粒子の運動エネルギーと気泡の存在の関係を直接的に見るには、同じシミュレーションに対してエネルギー E(t) と重心 H(t) を同時にプロットすれば良い。この結果を図 4.9に示す。このシミュレーションでは初期条件を底に堆積し、粒子がほとんど静止した状態か

Figure 4.9: 運動エネルギーと重心の関係。 N = 266, St = 20, u^{∞} = 0.3 での計算結果。このシミュレーションでは気泡の生成が観測されている。この気泡の吹抜けに伴って、この図で反時計周りに状態変化して行く。

らはじめている。この図から、シミュレーションの始めの部分ではエネルギーはそれほど大きくならず、重心だけが大きくなって行く様子が見られる。これは粒子は最初の堆積した状態のまま流れによって持ち上げられる過程に対応する。持ち上げられた粒子は下の部分から崩れて行き、底との間に出来た空隙が最初の気泡となって粒子の中を上昇して行く。気泡は粒子を巻き上げながら上昇して行き、抜け切ったあと吹き上げられた粒子が再び下降して層の中に戻って行く。重心の最大値は気泡が粒子層から抜け切る時に生じる。エネルギーの最大値は気泡が粒子層表面の粒子を大きく吹き上げて、その粒子が層に再び戻る過程の中で生じる。実際にエネルギーの最大値は重心の振動の頂点ではなく、ほとんど振動の中央の辺りで得られている。吹き上げられた粒子が再び堆積すると、粒子濃度が高くなって流体の摩擦が大きくなり、再び吹き上げられる過程がこの図から読みとれる。

4.2.2 気泡状態とチャンネル状態

粒子の重心 Hの変動は、気泡状態とチャンネル状態を区別する為の比較的良い指標である。つまり、流路が形成されれば、仮にその流路が不安定で断続的に崩壊と吹上を繰り返したとしても、重心の時間変動は比較的小さいはずである。 (この場合は純粋なチャンネル状態と言うよりも、むしろ気泡列と呼ぶべき状態かも知れない。 ) 実際に N = 266, u^{∞} = 0.3 の場合の重心の分散< (H− < H >)^{2} >のストークス数依存性を図 4.10に示す。ここで横軸のストークス数 St に対しては対数目盛をとっている。このように明らかにあるストークス

Figure 4.10: 重心の分散のストークス数依存性。 N = 266, u^{∞} = 0.3 での計算結果。 St_{c} = 5 ∼ 10 辺りで振舞いが大きく変わっている。 St < St_{c}では分散は小さく、従ってほとんど重心の変動が無く、チャンネル状態の性質を反映している。一方 St > St_{c}では分散が大きく、重心の振動が大きく、気泡状態の特徴を表している。このことから、 St_{c}がチャンネル状態と気泡状態の転移点に対応している事が分かる。

数 St_{c} = 5 ∼ 10 を境にし、ストークス数が小さい領域 St < St_{c}では重心の分散が小さい、つまり重心の変動が小さい事が分かる。このことは、系がチャンネル状態である事を反映している。一方ストークス数が大きい領域 St > St_{c}では分散が大きく、重心の変動が大きい事が分かる。実際系は気泡状態にあり、気泡の生成に伴って重心は大きく振動している。

気泡状態に対して行ったように St < St_{c}で系がチャンネル状態にある時の重心とエネルギーの関係を図 4.11に示す。この図と気泡状態の図 4.9を比較すると、振動を表す円の形が明らかに縦長になっており、重心の変動が小さい事を表している。 Figure 4.11: チャンネル状態での運動エネルギーと重心の関係。 N = 266, St = 2, u^{∞} = 0.3 での計算結果。図 4.9と比較すると、時間変化のループが縦長になっている事が分かる。

本節での考察から St_{c} 付近がチャンネル状態から気泡状態への流動状態の転移に対応する事が分かった。これから示す解析に於いていろいろな物理量のパラメータ依存性を見て行くが、殆どの結果のストークス数依存性に St_{c} = 5 ∼ 10 辺りに特異な変化が見られる (図 4.19,4.22, 4.24,4.27,4.29を参照)。これらの特異性は流動状態の転移との関係は 5.1.1節で議論する。

4.2.3 気泡生成周期

前節の議論から、運動エネルギーや重心の振動が気泡の存在に対応する事が分かった。このことから、気泡の生成周期 T_{bub}を運動エネルギー E(t) の時間的な振動周期から計算することが出来る。

具体的な気泡の生成周期 T_{bub}の計算は次のように行った。まず初期のトランジェントな振舞いを除く為、エネルギーの時間変化から最初のピークを越えた後の極小値をとる時間から以後を考察の対象とする。この領域でエネルギーの平均値を求め、エネルギーが平均値を (それよりも小さい値から) 越えた時間を調べ、その時間間隔を気泡周期 T_{bub}とした。この解析を流動状態 (気泡状態とチャンネル状態を含む) が得られたシミュレーションに対して行った。

流入速度依存性

図 4.12に気泡生成周期 T_{bub} の流入速度 u^{∞}依存性を示す。この結果から、一旦流動化す

Figure 4.12: 気泡生成周期の流入速度依存性 N = 266, St = 10 の計算結果。気泡の生成周期は流入速度に敏感ではない事が分かる。

ると気泡の生成周期は流入速度に敏感ではないことが分かる。図からは流動化開始速度 (臨界流入速度) 付近で幾分周期の増加傾向が読みとれる。しかし、誤差が大きくなっていることが示すように、この領域では１つのシミュレーションの内に気泡が小数しか観測されなかった。これは気泡の周期が長くなると言うよりは、むしろ初期の静止状態から流動化するまでに時間が掛かる事に依っている。

ストークス数依存性

図 4.13に気泡生成周期 T_{bub} のストークス数 St 依存性を示す。ここで横軸のストークス

Figure 4.13: 気泡生成周期のストークス数依存性。シミュレーション結果 (N = 266, u^{∞} = 0.3) とフィッティング結果 (A_{T} = 5.33 ± 0.15, B_{T} = 30.9 ± 4.4)。ストークス数に線形の関係が示唆される。

数 St に対しては対数目盛をとっている。この振舞は単純な線形依存性を示している。実際以下の式でうまくフィッティング出来る。

T_{bub} = A_{T} St + B_{T}  (4.3)

このフィッティングは、支配方程式 (3.129) に見られるようにストークス数 St が系の (速度変化の) 時間スケールを特徴付けている事によっている。但しこのフィッティングでのパラメータ B_{T} の値と物理的意味が不自然である。シミュレーションからは St = 0 の極限では、系はチャンネルを作り粒子は完全に静止してしまう。つまり T_{bub} = 0 に対応するからである。ここで仮定したフィッティング (4.3) 式は 4.2.2節で議論した St_{c}での相転移の特異性は T_{bub}には無いとする立場である。図 4.13から St_{c}辺りで特異な振舞いが見られるか見られないかは、判断の分かれるところである。

4.2.4 気泡の上昇速度の流入速度依存性

粒子濃度のプロファイルの時間発展から、気泡の上昇速度を計算する事が出来る。その結果を図 4.14に示す。流動化開始速度で気泡の速度が立ち上がり、流入速度の上昇にともなって気泡速度も大きくなることが分かる。巾的振舞を仮定し、以下の式でフィッティングを行う。

U_{bub}(u^{∞}) = A_{U}(u^{∞} − u_{c})^{α}  (4.4)

Figure 4.14: 気泡の上昇速度の流入速度依存性。シミュレーション結果とフィッティング結果 (α = 0.441 ± 0.040)。

このとき巾の指数はα = 0.441 ± 0.040 であり、 u_{c}は 4.3.1節で定義される流動化開始速度である。

4.3 定常状態

前節では系の示す動的な挙動について解析して来た。本節では動的な側面ではなく定常状態の系の平均的な振舞に着目する。ここで定常状態とは、先の気泡周期の計算の際に定義した様に、エネルギーの最初のピークを越えた後の最初の極小値をとる時間以後とする。これは一回気泡が系の中を吹き抜ける事で、ほとんどトランジェントな振舞が消えて系が定常状態になることが、実際のシミュレーションから見られる事による。得られた結果に対する解釈は 5.1節でまとめて行う。

4.3.1 平均運動エネルギー

運動エネルギーの定常状態での時間平均をĒとする。

Ē = \frac{1}{T} ∫_{T} dt E(t)  (4.5)

ここで T を上述のように定義した定常状態の時間とした。平均エネルギーは、コントロールパラメータである流入速度 u^{∞}とストークス数 St に依存する。

流入速度依存性

平均エネルギーĒの流入速度 u^{∞}依存性を図 4.15に示す。ある流入速度まではエネルギー

Figure 4.15: 平均運動エネルギーの流入速度依存性。 N = 266, St = 10 の結果に対する。臨界流速 u_{c}で平均エネルギーが流速に線形に立ち上がる振舞いが見られる。

はゼロである。つまり粒子は静止状態にある事が分かった。流速がそれ以上になると運動エネルギーは流速に比例する振舞を示す事が分かる。この流速を臨界流速 u_{c}とする。臨界流速付近でのエネルギーの振舞を拡大したものを図 4.16に示す。この図から臨界流速 u_{c}で

Figure 4.16: 臨界流速 u_{c}付近での平均運動エネルギーの流入速度依存性。 N = 266, St = 10 の結果に対する。臨界流速 u_{c}で平均エネルギーが流速に線形に立ち上がる振舞いが見られる。

の平均エネルギーの立上りは線形である事が分かる。これらの図は粒子数とストークス数が N = 266, St = 10 の結果である。その他のストークス数でのシミュレーション結果に対しても同様の振舞を示す。

平均運動エネルギー Ē の流入速度 u^{∞} 依存性を、次のような関数でフィッティングする。

Ē(u^{∞}) = A_{E}(u^{∞} − u_{c})  (4.6)

ここで、フィッティングパラメータとして A_{E}, u_{c} を導入した。これらはストークス数 St に依存する量である。 A_{E}(St) は運動エネルギーの流速に対する増加率に対応し、一方 u_{c} は先に導入した臨界速度、あるいは流動化開始流速に対応する。それぞれのストークス数 St 依存性を図 4.17、図 4.18 に示す。

Figure 4.17: 平均運動エネルギーĒ のフィッティングパラメータ A_{E}(St) のストークス依存性。 St_{c} ∼ 5 の辺りでピークを示しているようである。

Figure 4.18: 臨界流速 u_{c}のストークス数依存性。大きなストークス依存性は見られない。

A_{E}(St) はあるストークス数 St_{c} = 1 ∼ 5 でピークを持つ様に見える。つまり St ∼ St_{c} で系の平均運動エネルギーは u^{∞} に敏感になっていることが分かる。一方臨界流速 u_{c} はストークス数にほとんど依存せず、ほぼ u_{c} ∼ 0.1 である。この解析で得られた臨界流速が、本章の最初で示した図 4.1のパラメータ空間 (St, u^{∞}) での分類図に示した値である。ストークス数依存性

平均エネルギーĒのストークス数 St 依存性を図 4.19に示す。この図から、明らかにあ

Figure 4.19: 平均運動エネルギーのストークス数依存性。 N = 266, u^{∞} = 0.3 の結果に対する。明らかに St_{c} = 3 ∼ 4 でピークを示している。

るストークス数 St_{c} = 3 ∼ 4 でĒがピークを持つ事が分かる。 Ēの流入速度依存性に対して式 (4.6) が成り立つとすると、ある流速 u^{∞} での運動エネルギーのストークス数依存性 Ē(St) は、本質的に AĒ(St) と同じ振舞をしめす。実際、この図に見られるように、 Ē(St) と AĒ(St) は類似の振舞いをしている。 (AĒ(St) はフィッティングによる誤差を含んでいる為多少ばらついている。 )

4.3.2 流動層の平均重心

運動エネルギーに対して行った解析を粒子の重心 H(t) に対し行う。定常状態での平均重心をH̄とする。

H̄ = \frac{1}{T} ∫_{T} dt H(t)  (4.7)

流入速度依存性

平均重心H̄の流入速度 u^{∞}依存性を図 4.20に示す。ここでは St = 10 での結果を示す。

Figure 4.20: 流動層の平均重心の流入速度依存性 N = 266, St = 10.0 の計算結果。流速に線形に振舞う領域が見られる。これは実験で観測された結果に類似している。

平均エネルギーĒの解析で得られた臨界流速 u_{c} ∼ 0.1 で転移を示すことが分かる。転移点近傍での立ち上がりはエネルギーの時と異なり u^{∞} の線形ではない。しかし流入速度が u_{c} より少し大きくなると線形領域 0.15 ≤ u^{∞} ≤ 0.7 が存在するように見える。このような流入速度に比例して層が膨張することは実験的にも観測されている [45]。彼らの実験結果を図 4.21に示す。縦軸は層の膨張率で、横軸は流体の流入速度である。流体は空気 (密度は Figure 4.21: Furukawa & Ohmae (1958) の流動層の層膨張の流入速度依存性の実験結果。ここで流体は空気 (密度は 0.0012[g/cc])、粒子は Poly(vinyl acetate) bead (sphere) 粒子密度は 1.3[g/cc] である。図中のラベルは粒子の直径をミクロンで示してある。

0.0012[g/cc])、粒子は Poly(vinyl acetate) bead (sphere) 粒子密度は 1.3[g/cc] である。この実験結果の中でも粒径の小さい場合が本シミュレーションの結果に良く一致している事は、モデルの仮定 (低レイノルズ数) を考えると自然である。

ストークス数依存性

平均重心H̄のストークス数 St 依存性を図 4.22に示す。この図から、あるストークス数

Figure 4.22: 流動層の平均重心のストークス数依存性。シミュレーション結果 (N = 266, u^{∞} = 0.3) とフィッティング結果 (C_{H} = 3.087 ± 0.055, D_{H} = 16.51 ± 0.14)

St_{c} = 2 ∼ 3 以下では重心はストークス数に依らないことが分かる。一方 St_{c} を越えると重心がストークス数の増加に伴って高くなる。増加の程度は対数的である。実際以下の式でフィッティングを行うと良い一致を示す。

H̄(St) = C_{H}(u^{∞}) log(St) + D_{H}(u^{∞})  (4.8)

フィッティングに用いたシミュレーションは N = 266, u^{∞} = 0.3 で、その時のフィッティングパラメータの値は C_{H} = 3.087 ± 0.055, D_{H} = 16.51 ± 0.14 である。

4.3.3 自己拡散係数

実際のシミュレーションでは粒子の位置と速度の時系列が得られる。そのため粒子速度の時間揺らぎに対するパワースペクトルは容易に計算する事が出来る。 D_{ij}(ω) を粒子速度のパワースペクトルとし、以下のように定義する。

D_{ij}(ω) = \frac{1}{N} ∑_{α=1}^{N} ∫_{−∞}^{∞} dt e^{−iωt}U_{i}^{α}(t)U_{j}^{α}(0) =< Ũ_{i}^{α}(ω)Ũ_{j}^{α}(−ω) >_{α}  (4.9)

ここでŨ(ω) は粒子速度自身のフーリエ変換である。また粒子に対する平均を<>_{α}の形で略記した。

Ũ_{i}^{α}(ω) = ∫_{−∞}^{∞} dt e^{−iωt}U_{i}^{α}(t)  (4.10)

このパワースペクトル D(ω) を用いて自己拡散係数 D_{p}を求める事が出来る。

D_{p} = ∫_{0}^{∞} dt < U^{α}(0) · U^{α}(t) >_{α}= \frac{1}{2} (D_{xx}(0) + D_{zz}(0))  (4.11)

ここで今考察しているシミュレーションでは粒子の運動が x, z成分に限られている事を用いた。実際に行ったシミュレーションで得られた時系列がそれほど長くなく、またあるシミュレーションではトランジェントな部分があまりに長い為、データの収束性がそれほど良くないが、系の振舞を特徴付ける重要な統計量である。

流入速度依存性

自己拡散係数 D_{p}の流入速度 u^{∞}依存性を図 4.23に示す。この図から平均エネルギーや重

Figure 4.23: 粒子速度のパワースペクトルから求めた自己拡散係数の流入速度依存性。臨界流速 u_{c}で転移が見られる。

心と同様に臨界流速 u_{c} で転移を示していることが分かる。転移の立ち上がりは線形ではない。 u^{∞}の比較的大きな領域では線形的な依存性を示すようにも見える。しかし 4.3.4節での議論で考察するように、拡散係数の線形性のずれが有効粘性に重要な意味を持つ事が分かる。

ストークス数依存性

自己拡散係数 D_{p}のストークス数 St 依存性を図 4.24に示す。この振舞もエネルギーや重

Figure 4.24: 粒子速度のパワースペクトルから求めた自己拡散係数のストークス数依存性。 St_{c} = 5 ∼ 10 でピークを持つ事が分かる。心と同様にあるストークス数 St_{c} = 5 ∼ 10 で特異な振舞を示している。

この拡散係数に関しては、次の 4.3.4節で再び議論する。

4.3.4 有効粘性率

シミュレーションで観測された物理量に、速度のパワースペクトルから求めた拡散係数 D_{p} がある。拡散係数の意味は、熱平衡系での揺らぎの緩和を特徴付けるパラメータである。「揺動-散逸定理」により、この熱平衡での揺らぎの緩和が平衡からのずれの小さい線形応答の範囲での揺らぎの振舞と同じであるとの仮定から、非平衡での系の応答である抵抗係数が平衡系を特徴付ける拡散係数を用いて表現できる事が示される。最も簡単な「揺動-散逸定理」の表式は、アインシュタイン関係式として知られる以下の式である。

\frac{1}{ζ} = \frac{D}{kT}  (4.12)

ここで D, k, T, ζ はそれぞれ拡散係数、ボルツマン定数、温度、抵抗係数である。一般に抵抗係数は粘性率に比例する。ここではアインシュタインの関係式の抵抗係数を有効粘性率の定義として、粉体の有効粘性率を以下のように定義する。

μ_{e} = \frac{u^{∞}}{D_{p}}  (4.13)

ここで通常の温度の代わりに流入速度 u^{∞}を用いた。このことは平均エネルギーが流入速度に対して図 4.15に見られるような転移を示す事による。図 4.25に有効粘性率μ_{e}の流入速度 u^{∞}依存性を示す。この図から臨界流速 u_{c}で拡散係数 D_{p}がゼロから立ち上がる事に対応

Figure 4.25: 有効粘性率の流入速度依存性。シミュレーション結果 (N = 266, St = 10.0) とフィッティング結果 F = 13.7 ± 1.9, E_{f} = 0.113 ± 0.017。

し、有効粘性が無限大 (静止状態=固体状態) から減少して行く事が分かる。

Fukukawa&Ohmae(1958)[45] は実験的に流動層の有効粘性率を測定している。実験は膨張率の測定と同じ流体と粒子を用いている。実験方法は流動層の中に円柱状のドラムを入れ、一定速度で回転させた時ドラムにかかる負荷を測定することで粘性率を求めている。彼らの測定結果を図 4.26に示す。シミュレーション結果は彼らの実験結果と定量的には一 Figure 4.26: Furukawa & Ohmae (1958) の有効粘性率の流入速度依存性の実験結果。ここで流体は空気 (密度は 0.0012[g/cc])、粒子は Poly(vinyl acetate) bead (sphere) 粒子密度は 1.3[g/cc] である。図中のラベルは粒子の直径をミクロンで示してある。

致していないが、定性的な傾向は一致していることが分かる。彼らは実験結果に対して以下のアレニウス型の関数形でフィッティングしている。

μ_{e} = F e^{E_{f} /u^{∞}}  (4.14)

図 4.25の結果をフィッティングした結果は F = 13.7 ± 1.9, E_{f} = 0.113 ± 0.017 である。この結果は 5.1.3節の議論で示すように、液体の粘性率の温度依存性と類似した振舞いである。

一方、有効粘性率μ_{e}のストークス数 St 依存性を図 4.27に示す。この図から、あるストー

Figure 4.27: 有効粘性率のストークス数依存性。 N = 266, u^{∞} = 0.3 の計算結果。 St_{c} = 5 ∼ 10 で極小値を取っている事が分かる。

クス数 St_{c}付近で粘性率が減少しており、ストークス数が小さくなると粘性が増加する傾向が見られる。

4.3.5 速度分布

粉体を含め、散逸を持つ多粒子系の示す速度分布は、散逸を持たない系でのガウス分布と事なり、強い非ガウス性を示す事が知られている。

シミュレーション結果の粒子速度の分布関数 P (U_{x}) を考察する。水平成分を選んだのは対称性が良いからである。実際鉛直成分の速度分布は吹き上げる過程と沈降する過程の非対称性が大きく現れており、ガウス分布からのずれを議論するのに適当ではない。

図 4.28に分布関数の流入速度 u^{∞}依存性を、図 4.29にストークス数 St 依存性を示す。ここで、横軸の粒子速度は各分布の標準偏差で規格化してある。図 4.28からは、分布の流入速度 u^{∞}に対する大きな依存性は見られないようである。一方図 4.29からはストークス数 St が小さい程、非ガウス性が強くなる傾向がみられる。 Figure 4.28: 粒子速度分布関数の流入速度依存性。 N = 266, St = 10 の計算結果。横軸の速度は標準偏差で規格化している。

Figure 4.29: 粒子速度分布関数のストークス数依存性。 N = 266, u^{∞} = 0.3 での計算結果。横軸の速度は標準偏差で規格化している。

分布の非ガウス性を定量的にみる為に、４次のキュミュラント C_{4}(U_{x}) を計算した。

C_{4}(U_{x}) =< U_{x}^{4} > −3 < U_{x}^{2} >^{2} −4 < U_{x} >< U_{x}^{3} > +12 < U_{x} >^{2}< U_{x}^{2} > −6 < U_{x} >^{4}

(4.15)

この４次のキュミュラント C_{4}(U_{x}) を分散 (２次のキュミュラント)C_{2}で規格化した量 C_{4}/(C_{2})^{2} を、流入速度 u^{∞}とストークス数に対してプロットした (図 4.30と図 4.31)。この量はガウ

Figure 4.30: 速度分布の非ガウス性の流入速度依存性。 N = 266, St = 10 の計算結果。臨界流速 u_{c}の付近で非ガウス性が強くなっている事が分かる。

ス分布の時は当然ゼロとなり、指数分布の時は C_{4}/(C_{2})^{2} = 6 となる。この結果からは分布の形からははっきりしなかった流入速度依存性に対して、臨界流速 u_{c}に近い程非ガウス性が強くなっている傾向が分かる。一方ストークス数依存性に関しては、これまでの多くの結果と同様に St_{c} = 5 ∼ 10 付近で特異な振舞いを示している。今の場合、 St_{c}で非ガウス性が小さくなり、特にストークス数の小さい方に対して非ガウス性の明らかな増加傾向が見られる。

この２つの図を 4.3.4節の有効粘性のパラメータ依存性を示した図 4.25と図 4.27と比較すると、興味深い事に、非ガウス性の強さと粉体の有効粘性のパラメータ依存性が極めて良く一致している事が分かる。つまり、臨界流速 u_{c}から流入速度の増加に伴い有効粘性は急速に減少するが、速度分布の非ガウス性も同様に小さくなり、ストークス数 St_{c}で有効粘性は小さくなるが、その点で非ガウス性も小さくなり、ストークス数が減少すると有効粘性と非ガウス性が同様に大きくなっている。このような散逸の強さと非ガウス性の関連は非 Figure 4.31: 速度分布の非ガウス性のストークス数依存性。 N = 266, u^{∞} = 0.3 の計算結 ; 果。ある St_{c} = 5 ∼ 10 の付近で非ガウス性が弱くなっており、それより小さいストークス ; 数では非ガウス性強くなっている事が分かる。弾性多粒子系に於いても見られる [46, 47]。

55

議論

5.1 シミュレーション結果の解釈

5.1.1 ストークス数依存性

状態

ストークス数を変化させるとある値 St_{c} = 1 ∼ 10 を境に２つの流動状態–チャンネル状態と気泡状態–が得られる。これらの状態の区別は重心の分散の大きさで決める事が出来ることを 4.2.2節で示した。これは、チャンネル状態では系は流路を形成し、その付近で流動化している為重心の時間変動はほとんど見られないが、一方気泡状態では系の中を気泡が上昇し、それに伴って重心も変動するという事実による。実際図 4.10に示したように、 St < St_{c}では重心の分散は小さくチャンネル状態にあり、 St > St_{c}では重心の分散は大きく気泡状態にあることが分かる。このように St = St_{c}はチャンネル状態と気泡状態の境界に対応するが、定常状態での系の統計量–平均エネルギーや平均重心など–も St = St_{c}で振舞いが変化する事が分かる。

平均運動エネルギー

平均運動エネルギーのストークス数依存性は図 4.19に見られるように、 St < St_{c}ではストークス数の増加に伴ってエネルギーも増加し、 St > St_{c}ではストークス数の増加に伴ってエネルギーは減少する。

最初にチャンネル状態 (St < St_{c}) でのエネルギーの振舞いを議論する。 St = 0 では粒子の運動が無くなり、１つの安定な流路が形成される (図 5.1)。 St = 0 の時粒子速度は流体 Figure 5.1: ストークス数 St = 0 の時の安定なチャンネル状態。 N = 266, u^{∞} = 0.3, St = 0 の計算結果。時間は左から右に進み、 20 step 毎に示した。 1 step は無次元化された時間で 1 である。

力学的相互作用で決まる終端速度となり、この場合ほとんど流動化しない。ストークス数は粒子速度が終端速度に緩和する無次元化された時間であり、終端速度からのずれの大きさを表す。このずれが流路の崩壊、つまり流動化を引き起こす。従って、ストークス数が増加すると流動化が大きくなり、エネルギーも増加する事が理解できる。

Figure 5.2: 気泡に伴う粒子の対流。 N = 266, u^{∞} = 0.3, St = 10 の計算結果。周期セルを横に２つ並べたスナップショットである。

Figure 5.3: 気泡に伴う粒子の対流。 N = 266, u^{∞} = 0.3, St = 100 の計算結果。図 5.2と比較すると周期セルを横に２つ並べたスナップショットである。気泡の輪郭が曖昧になっている事が分かる。

気泡状態 (St > St_{c}) では気泡に伴って生じる粒子の二重対流がエネルギーの大半を担っている。この二重対流は、粒子が気泡の中では吹き上げられ、気泡の上部から気泡の周囲を沿って下降することにより形成されている (図 5.2)。この粒子の対流は、粒子の濃度差 (=気泡) により生じると理解できる。つまり流体は障害物である粒子の濃度が小さい領域を流れようとし、そのため気泡の中での流体の流れは相対的に強く、そこに飛び込んだ粒子は吹き上げられる。一方粒子濃度の高い領域 (気泡の外縁の両側) では相対的に流体の流れが弱く、粒子は下降する。従って、粒子の濃度差が小さい場合はエネルギー (対流の強さ) もまた小さくなる事が期待される。ストークス数 St は終端速度への緩和時間である。終端速度は流体力学的相互作用から決まり、粒子の配置に依存する。従ってストークス数は粒子の濃度分布の変化に対する粒子速度の応答時間である。ストークス数が増加すると濃度変化に対する応答が遅くなり、気泡と粒子相の境界が曖昧になって来る (図 5.3)。このため気泡状態ではストークス数の増加に伴ってエネルギーが減少すると考えられる。

以上の議論から、チャンネル状態から気泡状態に転移する点 (St = St_{c}) で平均運動エネルギーがピークを持つ事が分かる。平均運動エネルギーの大きさは流動化の大きさを表すから、 St = St_{c}で系は最も流動化していると言える。このことは自己拡散係数のストークス数依存性に見られるピークや有効粘性率の極少に反映している。

平均重心

平均重心のストークス数依存性は図 4.22に示した。

チャンネル状態 (St < St_{c}) で平均重心がほとんど変化していない。この事は、チャンネル状態が流体の流れる流路のまわりが流動化している状態であり、気泡状態とは異なり平均重心は流動化の大小にあまり依存しない。

一方気泡状態 (St < St_{c}) での平均重心の振舞いは、式 (4.8) のようにストークス数に対して対数的に増加する。これは次のように書き換える事も可能である。

St = exp (\frac{H̄ − D_{H}}{C_{H}} )  (5.1)

D_{H}は St = 1 での平均重心H̄である。したがって式 (5.1) はストークス数 St が平均重心の

Figure 5.4: 平均重心のストークス数依存性のフィッティングパラメータ C_{H}の流入速度依存性。流入速度に線形の依存性が見られる。

St = 1 からの増加分H̄ − D_{H} に対して指数的に増加する事を意味する。系の横幅 L_{x}, L_{y} は固定されているので、平均重心の増加分は層膨張率に対応する。 C_{H}の流入速度依存性を図 5.4に示す。この図から C_{H}は流入速度 u^{∞}に線形的に依存する事が分かる。今の場合ストークス数は系の振舞いを特徴付ける時間スケールである。実際に運動方程式 (3.129) に見られるようにストークス数は粒子速度の緩和時間を表し、また図 4.13に示した様に気泡の生成周期もストークス数に線形依存する。つまり系の特徴的時間 T が次のように書ける事を意味する。

T ∼ exp (C \frac{ΔV}{u^{∞}} )  (5.2)

ここで Cはある定数である。この式の意味は 5.1.3節で議論する。

5.1.2 流入速度依存性

流動化転移

流入速度 u^{∞}を増加させた時、ある流速 u_{c}を越えると系は静止状態から流動状態に転移する。図 4.18に示したように流動化開始速度 u_{c}はストークス数 St には依存していない。これは u_{c}が静止状態の粒子配置での平均沈降速度に対応し、平均沈降速度は粒子の配置にのみ依存する為である。このことは以下のように説明する事が出来る。

平均運動エネルギーの大きさは系の流動化の大きさを表す。粒子は大きさを持ち互いに重なり合う事は出来ない為に、流動化する為にはそれだけの隙間が必要となる。流動層は流体の一様な流れにより駆動されている。理想流動層はガリレイ不変であり、系に一様な流れを与える事 (系は静止した状態で) と、系を一定速度で下降させる事 (流体の流れは無い状態で) とは等価である。今底の役割をする固定粒子を導入しているが、後者の立場に立つと固定粒子の速度は−u^{∞}である。もし u^{∞}が静止状態での平均沈降速度よりも小さい場合は底の移動に対して粒子は追随する事が出来る。しかしもし u^{∞}が静止状態での平均沈降速度を越えると底と粒子層の間に隙間が生じる。この隙間が種となって系の流動化が生じると考えられる。つまり流動化開始速度 u_{c}は静止状態での平均沈降速度に対応する事が分かる。

底と粒子層の間に出来た隙間は一般に気泡となって層の中を上昇する。これは重い物質 (粒子) の中に軽い物体 (気泡) が置かれた時生じる自然な運動である。この意味で沸騰する水の中の気泡の振舞いと同じであり、気泡流動状態の挙動が沸騰を連想させる事は当然である。

平均運動エネルギー

平均運動エネルギーの流入速度依存性は図 4.15に示した。この結果から流動化開始速度 u_{c}よりも小さい流速 u^{∞}では平均運動エネルギーはほとんどゼロで、流動化開始速度を越えると流入速度に線形に立ち上がる事が分かる。流動化開始速度以下では系は流動化できない事は既に議論した。ここでは平均運動エネルギーの転移点での流入速度に対する線形性を議論する。

今議論しているのは個々の気泡の運動ではなく、この運動を平均した時の運動エネルギーである。流動化開始速度の議論と同様にガリレイ不変性を考慮し後者の立場に立つと、系に供給されるエネルギーの平均量は簡単に見積もる事が出来る。この場合流体は静止しており、固定粒子および自由に動ける粒子とも平均として−u^{∞}の速度で沈降している。従って１粒子当たり単位時間に mgu^{∞}のエネルギーを受け取っている事が分かる。つまりエネルギー供給は流入速度 u^{∞}に比例することが分かる。

定常状態でのエネルギーを決める為には、このエネルギー供給とエネルギー散逸をバランスさせなければならない。エネルギー散逸を見積もるには個々の気泡の挙動を正確に取り扱わなければならず、現時点では完全な議論は出来ていない。

以上の議論から、 u^{∞} < u_{c}で平均運動エネルギーĒ = 0 となり、 u^{∞} > u_{c}では (少なくとも転移点近傍では) 平均運動エネルギーĒが流入速度に対して線形の依存性を示す事は理解できたと言える。

5.1.3 液体とのアナロジー

本節では散逸系であり非平衡状態にある理想流動層の定常状態と熱平衡系、特に液体との間に見られるアナロジーを議論する。

空孔模型

理想流動層の定常状態で平均された系の振舞いは、流入速度 u^{∞}を有効温度と見倣す事で温度依存性が液体と極めて良い一致を示す事が分かった。例えば、平均運動エネルギー (内部エネルギー) が流入速度 (温度) に線形的に依存する事、流入速度 (温度) が大きくなると層膨張率 (体積) が大きくなる事、粘性率が流入速度 (温度) に対してアレニウス型の依存性を示す事など。

粘性率の温度依存性に対する次の式は Andrade により経験的に示された [48]。

μ ∼ exp (\frac{E_{f}}{kT} )  (5.3)

この式は次のように解釈する事が出来る。流動性は分子間に隙間 (自由体積) が出来るにつれて大きくなる。この隙間を作るのに必要なエネルギーが E_{f}である時、そのような隙間が出来る確率は exp(−E_{f}/kT) に比例する。粘性率μは流動性 (移動度) の逆数であるから (5.3) 式が成り立つ。この考え方は液体論の空孔模型で使われたものである [49]。

4.3.4節で示した理想流動層の有効粘性率 (4.14) 式は、流入速度 u^{∞}を有効温度としてこのメカニズムが成り立っている事を示している。 5.1.1節で議論した平均重心のストークス数依存性も同じ解釈で理解できる。つまり 5.1.1節の議論から系の特徴的時間スケールが (5.2) 式のようにΔV/u^{∞}に指数的に依存する事が分かったが、これは流入速度 u^{∞}を温度と見倣した時の、膨張率ΔV が生じるのに必要な時間と解釈する事が出来る。

このように定常状態での理想流動層の性質は液体論の空孔模型でコンシステントな説明が出来る。

揺動散逸定理

本研究で見積もった有効粘性率μ_{e}は流入速度 u^{∞}が温度と見倣すことが出来る事を用い、熱平衡状態と熱平衡近傍の非平衡状態に対して成り立つ揺動散逸定理の１つであるアインシュタインの関係を理想流動層に適用し、自己拡散係数から求めた量である。

この量が、実験で実際に系にシアー流れ (速度勾配) を作り出し応力を測定し求めた粘性率と定性的に一致する事は、理想流動層に於いてアインシュタインの関係が成立する事を意味する。また、粒子の速度分布に見られる非ガウス性の大きさがこの有効粘性率の大きさに対応している事も、アインシュタインの関係から求めた有効粘性率が妥当である事を示している。

この事は驚くべき結果である。なぜならば、理想流動層は本来散逸的な系であり、ここで揺動散逸定理が成立するべき理由は無いからである。この結果はまた (拡張された) 揺動散逸定理が、より普遍的に (ある種の非平衡状態でも) 成立している可能性を示唆している。非平衡統計力学の構築という文脈に於いても、理想流動層で揺動散逸定理が何処まで普遍的に成立するかと言う事を明らかにする事は重要である。その為には理想流動層に於いて直接シアー流れに対する応答を測定し粘性率を計算し、アインシュタインの関係が成立しているかどうかを調べる事が必要かつ重要である。

5.2 一様状態の欠如

今のところ理想流動層のシミュレーションからは「一様状態」が観測されていない。この事は、現実の流動層と比較した場合の理想流動層の最も大きな疑問であり問題点である。本節ではこの「一様状態の欠如」の理由を議論する。

実験結果から予想される「一様状態」が見られるパラメータ領域は、流速が流動化開始度を僅かに上回る程度で、アルキメデス数 (従ってストークス数) が小さい領域である。これは「固液流動層では気泡状態は見られず、一様状態のみ観測される」という経験的な傾向と矛盾しない。

最初に 3.3.5節で議論した、シミュレーションに用いたモデルには「粒子速度の緩和に多体効果が入っていない」ことが、一様状態の欠如に関係するかどうかを考える。その為に、本研究で示したモデルに「粒子速度の緩和の多体効果」が導入された時の振舞の変化をみる。慣性項が St′dU/dt → StR_{c}^{−1}dU/dt と変化する時、粒子間隔の接近した部分ではレジスタンスは大きくなるので、慣性項は小さくなり、一方粒子間隔の広い部分では慣性項は大きくなる。慣性項が小さくなると、粒子速度は流体力学的相互作用から決まる終端速度に近い値をとるようになる。終端速度に従う系はコロイドのような系であり、クラスタリングを起こす傾向がある。従って、慣性の多体効果を無視した近似が「一様状態」を消すものではないことが分かる。

現在考えられる最も有力な可能性は、考察したシミュレーションが鉛直面内に限られた運動しかできない事によるものである。粒子の運動を鉛直面内に限り、この面に垂直な方向の周期境界の間隔を粒子の直径としている為、静止状態 (最密充填) での粒子の周りの隙間が鉛直方向に揃った傾向になる。もし運動の制限が無くなると静止状態 (最密充填) での隙間は等方的になる。この違いは、正方格子と六方最密格子の関係に似ている。流動層では流体を鉛直方向に流すので、隙間が鉛直方向に揃っている場合は流体は流れやすくなる。このことは、シミュレーションで得られた流動化開始速度 u_{c} ∼ 0.1 が現象論の与える値よりも大きくなっている事からも分かる。実際同じ体積分率ϕ ∼ 0.60 での平均沈降速度の値を (3.101) 式から見積もると、 < U > /U_{0} ∼ 0.028 となりシミュレーションの結果が３倍も大きい事が分かる。理想流動層での流動化開始速度が、気泡が発生する程大きなものになっているとすると、一様状態が得られない事は説明できる。この描像の検証は、３次元的に粒子が運動できる状況でのシミュレーションを行う事で出来る。これは今後の課題である。

この他にも、そもそも理想流動層を構成する理想化の過程で無視した効果、例えばある程度大きい流速で粒子の周りに生じるカルマン渦による粒子の振動など、が「一様状態」を作る為に重要である可能性も考え得る。

結論と課題

6.1 結論

本研究では、これまで行われて来なかった流体を介して行われる粒子の多体相互作用を正確に扱った流動層のモデル化を行った。この目的の為に、現実の流動層から必要最小限の基本的なメカニズムを取り出した「理想流動層」を構成し、この系をモデル化した。

実際の計算で必要となる高濃度での粒子間の流体力学的な相互作用の計算の為に、粒子の慣性が無視できる場合に開発された Stokesian dynamics 法を拡張し、有限の粒子慣性をもつ粒子系に応用する為の数値手法を開発した。粒子慣性は流体力学的相互作用がもつ lubrication の特異性を繰り込んだ形で導入した。最近の研究との比較からこの繰り込みが、現実の系では様々な理由–例えば流体の連続体記述の破綻や粒子表面の微小な凹凸など–によってこの特異性が隠されてしまう事を意味する事が明らかになった。

またシミュレーションに用いた高濃度での流体力学的相互作用の計算方法を応用して平均沈降速度の理論解析を行い、実験結果と良く一致する結果を得た。

本研究ではこの計算手法を用いた「理想流動層」の「流動化」に関する数値シミュレーションを系統的に行い、系の示す相転移に伴う物性の変化を初めて系統的に調べた。

このシミュレーションの結果から、流動化転移が臨界流速 u_{c} で起こる事を数値的に示した。また流動層の膨張率の流入速度依存性が実験結果と定性的に良い一致を見た。これらの結果から流入速度 u^{∞}が「温度」の役割を持つ事が示唆された。「温度」として流入速度 u^{∞}を用いて揺動散逸定理の最も簡潔な表式であるアインシュタインの関係を理想流動層に適用し、自己拡散係数から有効粘性率を定義した。その結果この有効粘性率の流入速度依存性は実験で得られた結果と定性的に同じ振舞を示す事が分かった。

粒子の速度分布の解析からは、非ガウス性の強さの流入速度依存性とストークス数依存性を調べた。この結果速度分布の非ガウス性の強さは流入速度依存性だけでなくストークス数依存性も含め、有効粘性の強さと同一の傾向を示す事が分かった。

一方粒子の慣性の大きさを特徴付けるストークス数依存性については、ある特徴的な点 St_{c} が存在し、その点を境に St > St_{c}では気泡状態、 St < St_{c}ではチャンネル状態 (気泡列状態) というような系の流動形態の転移がある事を発見した。平均エネルギーのストークス数依存性は流動状態の違いによって説明され、 St_{c}でピークを持つ事が分かった。このことが St = St_{c}での自己拡散係数のピークや有効粘性率の極少に反映している。重心のストークス数依存性は、系の時間スケールが層膨張率と流入速度 (有効温度) の逆数の積に指数的に依存する事が分かった。

これらの結果は、流動層の流入速度依存性と液体の温度依存性との間の類似性を強く示唆する結果である。実際系の定常状態の振舞いは液体論の空孔模型のように系の流動性が粒子のまわりに出来る隙間 (自由体積) により生じると言う解釈で理解できる。流動層は散逸が支配的な非平衡系であるのに対し液体は通常熱平衡系として扱われるが、この両者の振舞いが似ている事は興味深い。

6.2 課題

本研究では熱平衡状態から遠く離れた散逸系である「理想流動層」の定常状態に関して、流入速度 u^{∞}を温度と見倣して揺動散逸定理 (アインシュタインの関係) を仮定し有効粘性率を見積もり、実験結果と定性的に一致する事を示した。この仮定の検証の為には、「理想流動層」で実際にシアー流れに対する系の応答を計算し粘性率を直接求め、揺動散逸定理により得られた結果と比較する事が必要である。

本研究では初期条件として主に静止状態を用いたシミュレーションを行った。これは固体の融解のような低温から高温に変化した際の系の変化を見たと言える。初期状態を静止状態ではなくランダムに分布する一様状態を用いたシミュレーションも少数 (転移点近傍は行っていない) ではあるが行った。これまでの所特に初期状態に依存する結果は得られていない。しかし特に転移点近傍での系の振舞いにおいてヒステリシスが存在するのか、或は転移が連続的であるのかを特定する事は、流動化転移を理解する上で重要であり、初期状態依存性を精密に考察する必要がある。

また本研究では計算時間の問題から、シミュレーションのシステムサイズ依存性に関しては限られた状況についてしか行っていない。層高、周期セルの横幅に関する系統的な解析は今後行う予定である。また一様状態の欠如との関連からも３次元的な状況での系の振舞いも考察するべき課題である。

謝辞

早川尚男先生 (現京都大学) には修士の時以来指導頂き、また 1996 年４月よりは京都大学の研究室での研究環境を提供して頂きました。佐々真一先生 (東京大学)、田口善弘先生 (東京工業大学)、西森拓先生 (茨城大学) には、粉体の研究会で多くの助言を頂きました。土井正男先生 (名古屋大学) には 1996 年の９月に１週間程研究室に滞在させて頂き、また同じ期間滞在された John F. Brady 教授との議論の機会を与えてくださり、感謝致します。両教授のアドバイスは、この研究をまとめる上で大きな助けになりました。早川先生、佐々先生と京都大学の水口毅先生、武末真二先生、佐野光貞先生、黒田明義君、大信田丈志君には本論文に関して貴重な御意見を数多く頂きました。都築俊夫先生 (東北大学) には、修士以来東北大学の研究室でのセミナーを通じ多く御指導を頂き、また京都大学で研究活動を行うに当たり様々な便宜を計って頂き、感謝致します。最後に東北大学理学部物性理論研究室と京都大学総合人間学部物性理論研究室の先生方、院生の方々、秘書の方々には、特に D3 の１年間私が京都で研究活動を行った事もあり、大変お世話になりました。

本研究は日本学術振興会の援助を受けて行われました。本研究で示したシミュレーションは東京大学物性研究所のスーパーコンピュータを用いて行いました。

PUBLICATIONS

[1] Kengo Ichiki and Hisao Hayakawa, “Simulation of granular particles in flow by the Stokesian-Dynamics method”, Int. J. Mod. Phys. B,7, 1899–1911 (1993).

[2] 市來健吾, “粘性流体中における分散多体系の動的挙動”, 東北大学大学院理学研究科修士論文, (1994).

[3] 市來健吾, “粘性流体中における分散多体系の動的挙動”, 物性研究, 62, 732–792 (1994).

[4] Hisao Hayakawa and Kengo Ichiki, “Statistical Theory of Sedimentation of Disordered Suspensions”, Phys. Rev. E, 51, R3815–R3818 (1995).

[5] Kengo Ichiki and Hisao Hayakawa, “Dynamical simulation of fluidized beds: Hydrodynamically interacting granular particles”, Phys. Rev. E, 52, 658–670 (1995).

[6] Kengo Ichiki, “Hydrodynamic effects in density waves of granular flows”, (submitted to Phys. Rev. E.)

Bibliography

[1] 田口善弘. 砂時計の七不思議. 中公新書, 1995.

[2] P. Bak, C. Tang, and K. Wiesenfeld. Phys. Rev. A, 38:364, 1988.

[3] H. M. Jaeger, S. Nagel, and R. P. Behringer. Phys. Today, April:32, 1996.

[4] H. M. Jaeger, S. Nagel, and R. P. Behringer. Rev. Mod. Phys., 68:1259, 1996.

[5] M. Faraday. Philos. Trans. R. Soc. London, 52:299, 1831.

[6] Y-h. Taguchi. Phys. Rev. Lett., 69:1367, 1992.

[7] J.A. C. Gallas, H.J.Herrmann, and S. Sokołowski. Phys. Rev. Lett., 69:1371, 1992.

[8] F. Melo, P.B. Umbanhower, and H.L. Swinney. Phys. Rev. Lett., 72:172, 1994.

[9] F. Melo, P.B. Umbanhower, and H.L. Swinney. Phys. Rev. Lett., 75:3838, 1995.

[10] 鞭巌, 森滋勝, 堀尾正靭. 流動層の反応工学. 培風館, 1984.

[11] 千葉忠俊, 吉田邦夫編著. 流動層概論. 朝倉書店, 1996.

[12] D. Geldart. Powder Technol., 7:285, 1973.

[13] D. A. Drew. Annu. Rev. Fluid Mech., 15:261–291, 1983.

[14] G. K. Batchelor. J. Fluid Mech., 193:75–110, 1988.

[15] S. Sasa and H.Hayakawa. Europhys. Lett., 17:685–689, 1992.

[16] M. F. Goz. J. Fluid Mech., 240:379–404, 1992.

[17] T. S. Komatsu and H. Hayakawa. Phys. Lett. A, 183:56–62, 1993.

[18] D. Gidaspow. Multiphase Flow and Fluidization. Academic Press, London, 1994.

[19] Y.Tsuji, T.Tanaka, and T.Ishida. Powder Technol., 71:239, 1992.

[20] Y. Tsuji, T. Kawaguchi, and T. Tanaka. Powder Technol., 77:79–87, 1993.

[21] T. Tanaka, T. Kawaguchi, and Y. Tsuji. Int. J. Mod. Phys. B, 7:1889–1898, 1993.

[22] B. J. Alder and T. E. Wainwright. J. Chem. Phys., 31:459–466, 1959.

[23] I. Goldhirch and G. Zanetti. Phys. Rev. Lett., 70:1619, 1993.

[24] J. Happel and H. Brenner. Low Reynolds number hydrodynamics. Martunus Nihhoff, Dordrecht, 1973.

[25] D. J. Jeffrey and Y. Onishi. J. Fluid Mech., 139:261–290, 1984.

[26] G. Bossis and J. F. Brady. J. Chem. Phys., 80:5141, 1984.

[27] J. F. Brady and G. Bossis. J. Fluid Mech., 155:105, 1985.

[28] L. J. Durlofsky, J. F. Brady, and G. Bossis. J. Fluid Mech., 180:21, 1987.

[29] H. Hayakawa and K. Ichiki. Phys. Rev. E, 51:R3815, 1995.

[30] C. W. J. Beenakker. J. Chem. Phys., 85:1581–1582, 1986.

[31] J. F. Brady and G. Bossis. Annu. Rev. Fluid Mech., 20:111–157, 1988.

[32] J. F. Brady, R. J. Phillips, J. C. Lester, and G. Bossis. J. Fluid Mech., 195:257–280, 1988.

[33] S. Kim and S. J. Karrila. Microhydrodynamics. Butterworth-Heinemann, Boston, 1991.

[34] G. K. Batchelor. J. Fluid Mech., 52:245–268, 1972.

[35] J. F. Brady and L. J. Durlofsky. Phys. Fluids, 31:717–727, 1988.

[36] J. C. Bacri, C. Frenois, M. Hoyos, R. Perzynski, N. Rakotomalala, and D. Salin. Europhys. Lett., 2:123, 1986.

[37] R. Buscall, J. W. Goodwin, R. H. Ottewill, and T. F. Trados. J. Colloid Interface Sci., 85:78, 1982.

[38] C. G. de Kruif, J. W. Jansen, and A. Vrij. A sterically stabilized silica colloid as a model supramolecular fluid. pages 315–346. Wiley-Interscience, New York, 1987.

[39] S. E. Paulin and B. J. Ackerson. Phys. Rev. Lett., 64:2663, 1990.

[40] J. Z. Xue, E. Herbolzheimer, M. A. Rutgers, W. B. Russel, and P. M. Chaikin. Phys. Rev. Lett., 69:1715, 1992.

[41] K. Ichiki and H. Hayakawa. Phys. Rev. E, 52:658, 1995.

[42] R. R. Sundararajakumar and D. L. Koch. J. Fluid Mech., 313:283–308, 1996.

[43] A. S. Sangani, G. Mo, H-K. Tsao, and D. L. Koch. J. Fluid Mech., 313:309–341, 1996.

[44] K. Ishii. J. Phys. Soc. Jpn., 46:675–680, 1979.

[45] J. Fukukawa and T. Ohmae. Ind. Eng. Chem., 50:821–828, 1958.

[46] Y-h. Taguchi and H. Takayasu. Europhys. Lett., 30:499–504, 1995.

[47] J. J. Brey, F. Moreno, and J. W. Dufty. Phys. Rev. E, 54:445–456, 1996.

[48] E. N. da Andrade. Phil. Mag., 17:497,698, 1934.

[49] 戸田盛和, 松田博嗣, 樋渡保秋, 和達三樹. 液体の構造と性質. 岩波書店, 1976.