微小粒子の乱流中の衝突特性に対する 比重の影響
横嶋 哲
1・島田 佳昭
21正会員 静岡大学准教授 学術院工学領域 数理システム工学系列 （〒432-8561浜松市中区城北3-5-1）
E-mail: [email protected]
2非会員 静岡大学 大学院総合科学技術研究科（修士課程） 工学専攻 数理システム工学コース （同上）
流体運動に完全追従するラグランジュ粒子や十分に重い慣性粒子について，乱流中の粒子衝突特性を調べた 例は多く存在する反面，粒子比重の影響を調べた例は存在しない．本研究は，比重が一様等方乱流中の微小粒 子のクラスタ形成に及ぼす影響を検討した先行研究（土木学会論文集A2, 72(2), I 459, 2016）の発展とみなせ，
比重が1.005から1000の範囲内で粒子の衝突特性を調べた．流体中の粒子挙動は，Basset履歴項を無視した
Maxey-Riley方程式によって表現した．比重の増加とともに粒子衝突頻度は大幅に高まること，その主要因はク
ラスタ形成に伴う粒子数密度の局所・瞬間的な高まりにあることが示された．粒子衝突頻度を相対接近速度と 動径方向分布関数で推定するモデルの妥当性も評価し，このモデルは粒子衝突頻度を非常に良好に予測できる ことが示された．
Key Words: particle clustering, preferential concentration, particle collision, effects of mass density ratio, homoge- neous isotropic turbulence, direct numerical simulation
1. はじめに
流体よりも十分に重い微小粒子（慣性粒子）群が乱流 中で，たとえその乱れが統計的に一様であっても，偏分 布（クラスタリング）することはよく知られている1)． 他方でそのクラスタリング特性が，粒子比重γを系統 的に変化させた場合にどのような影響を受けるかにつ いては，横嶋ら2)によって初めて検討された．彼らはス トークス数Sηの値を，慣性粒子の場合に視覚的にも容 易にクラスタ形成が認められる約0.7に固定して，比重 γの値を1.005から100の範囲で変化させた．彼らは，
(i)粒子がわずかに重いγ= 1.005の場合には，粒子速 度の初期条件に起因する過渡的なクラスタ形成3)が認 められるものの，その程度は非常に弱く，視覚的な確 認は困難であること，(ii)γの増加とともにクラスタ形 成が促進され，γ= 2.0前後で視覚的にもクラスタの存 在が容易に認識できること，および(iii)γ= 100程度に 達するとクラスタ分布はもはやγの値に依存せず，慣 性粒子とみなせること，を指摘した．
クラスタ形成は粒子の数密度を局所・瞬間的に高める ため，粒子間の接近/遭遇/衝突頻度を増加させることは 容易に想像できる．これは，局地的な集中豪雨をもた らし得る大気中での水滴の急激な衝突・合体速度4)や，
微生物と餌の遭遇・捕食頻度5),6),7)，あるいは粒子の凝 集効率8)などの多種多様な問題に密接に関連するため，
粒子比重が衝突特性に及ぼす影響について理解を深め
ることは重要である．しかしながら粒子の衝突特性に 関する検討例はこれまでのところ，流体に完全追従す るラグランジュ粒子（γ= 1.0，例えば文献9),10)）や慣 性粒子（γ→ ∞，例えば文献11),12)）に限られている．
本報では粒子比重が微小粒子の衝突特性に及ぼす影 響について検討する．具体的には粒子衝突頻度Γと，そ れに密接に関連するパラメータである衝突時動径方向 相対接近速度⟨|wr|⟩および動径方向分布関数g(R)につ いて詳しく調べた．検討対象とするデータは前報2)の 数値実験に基づくものの，以下の2点について改善を 図った：(a)前報2)では初期に粒子を一様ランダムに流 れ場に投入し，比重γの異なるそれぞれのケースにお いてコルモゴロフ時間の約60倍に渡って粒子挙動を追 跡した．粒子場はこの時間内に統計的定常状態に概ね 達したように見えるものの（例えば，前報2)の図–4），
その確認が十分になされたとは言えない．本研究では 積分時間をコルモゴロフ時間の180倍にまで延ばし，
完全に統計的定常な条件下で粒子の衝突特性を議論す る．(b)前報2)では比重γを100まで増加させ，上述の
(iii)の結論を得た．他方で大気中での水滴の振る舞い
を考える場合，粒子比重は800に達する．本報では比 重γ= 200,500,1000についても数値実験を新たに実施 し，γの値が十分に大きい場合の慣性粒子への収束の程 度について詳細に検討する．
2. 支配方程式と数値計算法
(1) 粒子の運動方程式
Squires and Yamazaki13)の再現実験を行った前報2)と 同じく，本研究でも微小粒子の運動方程式として，い
わゆるMaxey-Riley方程式14)（以下，MR方程式）を
採用する．MR方程式は一様でない流速場に投入され た，単一の微小球形剛体粒子に対する運動方程式で，粒 子に作用する力にはストークス抗力，圧力勾配力，（浮 力も含めた）有効重力，付加質量力，およびBasset履 歴力が含まれる．質点近似に基づくMR方程式が成立 するためには，粒径が流れ場の最小空間スケールより も小さくなければならない．Faxen補正項を付加する ことでこの制約はいくらか緩和される．乱流中の微小 粒子のクラスタリングや衝突特性を扱った先行研究の
大半1),11),12),15),16),17)では，大気中の液滴のように粒子比
重γが十分に大きな場合を扱うことを主たる理由とし て，圧力勾配，付加質量，およびBasset履歴の影響が，
Faxen補正項と併せて無視されてきた．他方でSquires
and Yamazakiでは，莫大な計算負荷を伴うBasset履歴
項は無視されたものの，圧力勾配と付加質量の寄与を考 慮したMR方程式が採用されており，本研究でも前報2) に引き続き，そのアプローチを導入した．同様の試み としてMaxey18)および杉山ら19)も挙げられる．Squires
and Yamazakiの表記に倣えば，本研究で扱うMR方程
式は次式で表される：
dxp
dt =vp, (1)
dvp
dt =α(u(xp)−vp−we3) +3
2R∗du(xp) dt . (2) ここでvpとxpは個々の微小粒子の速度と位置を，uは 背景流速を表し，重力加速度gは−x3方向に働くもの とする．eiはxi方向の単位ベクトルを示す．α−1, R∗ はそれぞれ速度差の緩和時間および粒子と流体の質量 密度比に関するパラメータ，wは静止流体中における 粒子の沈降速度を表し，時間スケールτ∗≡(2/9)a2/νf
（aは微小粒子の半径，νfは流体の動粘性係数），g，比 重γ≡ρp/ρf（ρp, ρfはそれぞれ粒子および流体の質量 密度）を用いて
α−1= (
γ+1 2 )
τ∗, R∗= 1
γ+12, w=τ∗(γ−1)g (3) と表される．比重γが十分に大きい場合には，式(2)は 粒子の緩和時間τp≡γτ∗を用いて
dvp
dt =u(xp)−vp
τp −ge3. (4)
のように，重力環境下での微小慣性粒子の運動方程
式17),20)に帰着する．
(2) 粒子の衝突と接触
本研究では同一粒径の多数の微小粒子が乱流中で衝 突する現象を扱う．粒子の衝突半径をRとすれば，衝 突を以下のように定義できる10)：ある時刻tにおける粒 子iと粒子jの中心間距離をdij(t)とすれば，時刻t1に 粒子iと粒子jが衝突したとは，dij(t1−∆t)> Rかつ dij(t1)≤Rが成り立つことを意味する．また，0.99R≤ dij(t)≤1.01Rを満たすとき，粒子iと粒子jは時刻t に接触状態にある，とする．
(3) 数値実験条件
数値実験の条件は前報2)とほぼ一致する．すなわちナ ヴィエ・ストークス方程式の直接数値シミュレーション
（DNS）によって，テイラー・レイノルズ数が約4213)の 統計的に定常な一様等方乱流場を再現した．多数の微 小粒子を一様ランダムな初期配置となるように投入し，
その衝突特性を調べた．ラグランジュ粒子や慣性粒子 を対象とした先行研究9),11),12)では，DNSから得られる ある瞬時速度場（凍結乱流場）を用いて粒子挙動が評 価された．凍結乱流場では乱流渦がその強度を保って 存在し続ける（渦の寿命が無限大）ため，衝突頻度が 過大評価されることが指摘されている11),12)．本研究で は，ナヴィエ・ストークス方程式の時間発展によって乱 流渦が成長/減衰する，いわば生きた乱流場をDNSで 再現し，粒子の振る舞いを追跡した．1.節で述べたよ うに，本研究では数値実験のシミュレーション継続時 間を，前報2)で用いられたコルモゴロフ時間τηの約60 倍から，その3倍の180τηに拡大し，統計的に定常な 状態で十分な量のサンプルデータを確保した．
境界条件の影響を強く受ける大スケールでは非定常 性や異方性が認められる場合にも，カスケード過程を 繰り返して達する小スケールな乱流場に着目すれば，そ の乱れ特性は概ね統計的に定常で一様等方と近似でき る．本研究では衝突半径Rがコルモゴロフ長ηに等し い場合12)を検討対象とする．衝突に最も影響をもたら すのは衝突半径と同程度のスケールの乱流変動である ので，背景乱流場に対して統計的に定常な一様等方乱 流を仮定することは，以下の議論とそこから得られる 知見の一般性を損なわない．
本研究では粒子比重γを1.005から1000の範囲で 変化させた．この際，慣性粒子のクラスタ形成を支配 するパラメータとして知られている16)，流体の最小ス ケールであるコルモゴロフ時間τηに関するストークス 数Sη≡τp/τηの値は0.668で一定に保った．各ケース における粒子の物理パラメータを表–1にまとめた．重 力加速度の影響を無視した．これは，前報2)において，
ストークス数を固定して比重を増加させると粒子の終 端速度パラメータwが飛躍的に増加した結果，クラス
表–1 比重γを系統的に変化させた場合の粒子の物理パラメー タ．ストークス数Sηの値が一定となるようにγa2の値 を(9/2)Sηη2に固定した．
γ τ∗
τη
gτη uη
ατη w uη
R∗ Sη≡ τp τη
1.005 0.6645 0 1.00 0 0.664 0.668 1.1 0.6071 0 1.03 0 0.625 0.668 1.2 0.5565 0 1.06 0 0.588 0.668 1.3 0.5137 0 1.08 0 0.556 0.668 1.4 0.4770 0 1.10 0 0.526 0.668 1.5 0.4452 0 1.12 0 0.500 0.668 2.0 0.3339 0 1.20 0 0.400 0.668 5.0 0.1336 0 1.36 0 0.182 0.668 10 0.0668 0 1.43 0 0.095 0.668 20 0.0334 0 1.46 0 0.049 0.668 30 0.0223 0 1.47 0 0.033 0.668 50 0.0134 0 1.48 0 0.020 0.668 100 0.0067 0 1.49 0 0.010 0.668 200 0.0033 0 1.49 0 0.005 0.668 500 0.0013 0 1.50 0 0.002 0.668 1000 0.0007 0 1.50 0 0.001 0.668
タ形成が阻害されることを確認したことによる．
(4) 数値計算法
3方向に周期境界条件を課した立方体領域内における 背景乱流場を483のフーリエ・ノードで離散表現し，擬 スペクトル法でナヴィエ・ストークス方程式を解くこと でその時空間発展を数値的に再現した．計算条件と得 られた乱流場の主要な統計量を表–2に示す．時間刻み 幅∆tには，粒子衝突を見過ごすことの無いように，約 0.03τηと非常に小さな値を与えた．Np= 1.5×483個 の微小粒子を初期に一様ランダムに配置し，粒子の存 在が流体運動に及ぼす影響を無視する，いわゆる1-way カップリング・シミュレーションを行った．粒子の初期 速度vp,0には粒子位置での流体速度u(xp,0)を与えた．
粒子速度の初期条件が粒子挙動に及ぼす影響について は前報2)を参照されたい．また，流体場および粒子場の 数値シミュレーションの詳細は横嶋ら10)を参照された い．利用した数値実験プログラムの有効性は付録Iに て詳細に確認された．
3. 数値実験結果
(1) 粒子クラスタリング
粒子の空間分布の偏りの程度を定量的に表す指標の ひとつである空隙率PV16)の比重依存性を図–1に示す．
空隙率PVは粒子を全く含まない微小体積が計算領域
全体に占める割合を表す量で，本研究では，(i)計算領 域を微小体積要素（ここでは背景乱流場の数値計算に 用いたフーリエ・ノード数と整合させ，483の微小立方 体セル）に分割し，(ii)微小粒子を1つも含まないセル 体積を足し合わせ，(iii)全計算領域の体積で規格化，す ることで評価した．
図–1(a)には全16ケースから得られたPVの時系列 分布を示す．時刻t= 0に一様ランダムに計算領域内に 投入された粒子はエネルギー保有渦の代表時間スケー ルT′ ≈11τηの数倍の時間をかけてクラスタを形成し，
準定常状態に達することがわかる．背景乱流場の振る 舞いは全ケースに共通であるため，PVの増減の様子が 似通っている．比重が増すとPVは増加し，ある分布 に収束するような振る舞いが認められる．図–1(b)には γ ≥ 20の7ケースの時系列分布のみを拡大して示す．
γ ≥100以上の4ケース間に差異を見出すことはこの スケールでも困難であり，瞬間的には，例えばγ= 200 のケースが最大となる場合も見受けられた．いずれの ケースにおいても系が十分に統計的定常状態に達して いる，区間120 ≤t/τη ≤180での時間平均値P¯v のγ 依存性をP¯v,γ=1.005で規格化して，図–1(c)にプロット した．γが増すにつれて，粒子を含まない計算セルの 割合はγ = 1.005の場合と比較して最大で2.6倍程度 増加し，P¯vは0.6弱程度の値に漸近する様子が見られ る．図–1(c)中の差し込み図にはP¯v,γ=1000からのズレ をP¯v,γ=1000で規格化した結果を示した．ズレはγが増 すにつれて単調に減少し，γ≥30ではズレは1%未満 に収まることがわかる．
クラスタ形成の様子を視覚的に確認するため，ある 瞬間（ここでは全てのケースに対してt= 180τη）の粒 子数密度np（各計算セルに含まれる粒子数）の2次元 平面内分布を図–2(a)-(p)に示す．比重γ= 1000の場合 にnpが最大値66をとった点を通る2次元平面を選び，
他の15ケースにおいても同一平面上でnpの分布を抽 出した．総計算格子数が483個に対して，計算領域内 の総粒子数Npが1.5×483個であるから，粒子分布が 完全に均一であれば平均的にnp= 1.5となることに留 意されたい．図–2中のccは，それぞれのケースのnp
と比重γ= 1000の場合のnpの空間分布の相関係数の 値を表す．図–2(q)には(a)-(p)と同一平面上での渦度ベ クトルω ≡ ∇ ×uの絶対値のカラーコンター図を示 した．流体に対して十分に重い粒子は慣性の影響で渦 度の小さな領域に集積することが指摘されている1),15)． 図–2(q)と(p)を比較すると，確かに渦度|ω|の大きな 領域には粒子は存在せず，それら高渦度領域に挟まれ た低渦度領域に選択的に粒子が集積するように見える．
粒子比重γ= 1.005と1.1の結果を比較すると，図– 1に示されたPv の値では両者の間に差異が認められ
表–2 背景乱流場の代表的な統計量．Nは空間離散化に用いた総スペクトルモード数，L, Tは計算領域（立方体）の一辺の長 さとDNS継続時間，∆xと∆tはDNSの時空間分解能，η, τη, uηおよびL′, T′, u′はそれぞれコルモゴロフ・スケー ルおよびエネルギー保有渦の長さ,時間,速度スケールを表す．
Reλ N L/η T /τη ∆x/η ∆t/τη L′/η T′/τη u′/uη 42.3 483 108.8 180.0 2.27 0.030 36.1 10.93 3.31
γ
100 101 102 103
¯ P
v/ ¯ P
v,γ=1.0051 2
(c)
γ
100 101 102 103
(¯Pv,γ=1000−¯Pv)/¯Pv,γ=1000 10-3 10-2 10-1 100
図–1 (a)空隙率Pvの時系列分布．γ= 1.005の場合にPvの値が最も低く，γが増すにつれてPvも増加．(b) (a)の拡大図．(c) 空隙率Pvの区間120≤t/τη≤180における時間平均値P¯vの粒子比重γに対する依存性を，P¯v,γ=1.005(= 0.223)で規 格化して表示．差し込み図はP¯v,γ=1000(= 0.584)からのズレをP¯v,γ=1000で規格化した結果を示す．
る（P¯v,γ=1.005= 0.223,P¯v,γ=1.1= 0.241）のに対して，
図–2(a),(b)から有意な差を見出すことは難しく，視覚的 にもクラスタ形成が認められるのはγ= 1.4（P¯v,γ=1.4= 0.340）程度から，と言える．他方でγ= 10の場合の γ= 1000の結果との相関係数は0.76に達し，クラスタ の骨格はγ= 10でも十分に再現される．
(2) 粒子衝突頻度
粒子の衝突因子ΓDNSの時系列分布を図–3(a)に示す．
衝突因子ΓDNS は，数値実験中のある時間ステップ中
（時間∆t）に観測された粒子衝突数Nc，計算領域の体
積V ≡L3，粒子数密度ρp ≡Np/V（図–2で扱った，
計算セル単位で評価される局所的な数密度npとは異な り，計算領域全体を対象として巨視的に評価される，単
位体積あたりの粒子数）を用いて
ΓDNS≡Nc/(∆t·V ·ρ2p/2) (5) のように求まる．すなわち衝突因子とは，単位時間内 に，単位体積に含まれる粒子（粒子数ρp）からなるペ ア（ペア数ρpC2 ≈ ρ2p/2）の衝突に至る割合，を意味 する．図–3(a)に示されるΓDNS時系列分布においても，
空隙率Pvの場合と同様に，比重とともに衝突頻度が大 きく変化するのはγ = 10以前であり，γ= 20以上の ケース間に有意な差を見出すことは，たとえ図–1(b)の ような拡大図によっても，難しい．
図–3(b)には，図–1(c)と同様に，120 ≤ t/τη ≤180 での時間平均値Γ¯DNSのγ依存性をΓ¯DNS,γ=1.005で規 格化して表示した．γが増すにつれて衝突頻度は増加 し，γ = 1.005の場合と比較して最大で約6.2倍に達
(a)γ= 1.005,cc = 4×10−4. (b)γ= 1.1,cc = 0.08. (c)γ= 1.2,cc = 0.14.
(d)γ= 1.3,cc = 0.19. (e)γ= 1.4,cc = 0.23. (f)γ= 1.5,cc = 0.27.
(g)γ= 2,cc = 0.39. (h)γ= 5,cc = 0.65. (i)γ= 10,cc = 0.76.
図–2 (a)-(p)粒子数密度npの2次元平面内分布のスナップショット．γ= 1000の場合にnpが最大値をとる点を通る平面を選
び，他の15ケースにおいてもその面上でnpの分布を描画．図中のccは，それぞれのケースのnpと比重γ= 1000の 場合のnpの空間分布の相関係数の値を表す．(q)は(a)-(p)と同一平面上での渦度ベクトルω≡ ∇ ×uの絶対値のモノ クロ・コンター図（コルモゴロフ・スケールで規格化）を表す．
する．
(3) 衝突因子モデル
前節で扱った衝突因子ΓDNSの平均量は次式で推定 できる12)：
Γsph= 2πR2⟨|wr|⟩g(R). (6) wrは2粒子の衝突時（dij =R）の相対接近速度w ≡ v(j)p −v(i)p と粒子中心の相対位置ベクトルR≡x(i)p −x(j)p
（|R| = R）を用いて次式のように評価できる（図–4 参照）：
wr≡w·R/R. (7)
ある瞬間に接触状態（2.(2)項で定義）にある全ての粒 子ペアに対してwrを評価し，その絶対値をアンサンブ ル平均することで⟨|wr|⟩を求めた．また，動径方向分
布関数g(R)はWanget al.12)に倣って次式で算出した：
g(R) = Ncontact·V
NpC2·Vs
. (8)
ここでVs = (4/3)π(1.013−0.993)R3は2粒子の接触 検出用の検査領域の体積，Ncontactはある瞬間に検出 された接触数を意味する．すなわちg(R)は，（0.99R≤ dij ≤1.01Rを満たし）接触状態にある粒子ペア数の，
粒子の偏分布による，粒子が均一に計算領域内に分布 した場合に対する増減の程度を表す．式(6)は，ある粒 子を中心に有する半径Rの球表面（面積4πR2）を通っ
(j)γ= 20,cc = 0.80. (k)γ= 30,cc = 0.81. (l)γ= 50,cc = 0.82.
(m)γ= 100,cc = 0.82. (n)γ= 200,cc = 0.82. (o)γ= 500,cc = 0.82.
(p)γ= 1000. (q) vorticity intensity.
図–2 続き．
て球内部に向かう流束と外部への流束が釣り合う，と の仮定に基づく．Wanget al.12)は内部に向かう流束が わずかに大きいこと，したがって式(6)は衝突因子を最 大で数%程度過少評価することを指摘した．
図–5はDNSから算出した⟨|wr|⟩，g(R)，およびそれ らに基づく衝突因子の推定値Γsphの時系列分布を示す．
比重が大きな（例えばγ ≥2）ケースにおいて，粒子 の衝突時動径方向相対接近速度⟨|wr|⟩が初期に急増す ることは，慣性の大きな粒子が乱流渦の作用によって 特定の領域に集積する過程で生じたと理解できる．比
重γ = 1.0のラグランジュ粒子では動径方向相対接近
速度⟨|wr|⟩は0.2(R/η)uηで良好に近似できることが知 られており10)，本数値実験においてもγ ≈1.0な場合 には上記の知見がよく当てはまることがわかる．比重 γ= 1.005のときg(R)は1前後で推移しており，この
場合には粒子の集積がほとんど生じないため，粒子の接 触が促進されないことを意味する．式(6)から求まった 衝突因子Γsphの経時変化の様子（図–5(c)）は，図–3(a) のΓDNSの挙動をよく再現しており，これは衝突因子推 定モデル（式(6)）が統計的定常時だけでなく，粒子の 集積形成過程でもよく機能することを意味する．
次に，⟨|wr|⟩，g(R)，Γsphの時間平均値の比重依存性 を図–6に示す．図–6(c)より，ΓDNSとΓsphのズレは全 てのケースで±0.2%程度に収まっており，粒子衝突頻 度は式(6)によって極めて精度よく推定できる．図–3(b) では，ラグランジュ粒子（γ= 1.005）と比較して慣性 粒子（γ = 1000）の粒子衝突頻度が約6.2倍高まるこ とが示された．式(6)の枠組みではこの内訳として，粒 子の衝突時の動径方向相対接近速度⟨|wr|⟩の高まりの 寄与が約1.27倍，クラスタ形成による局所的な粒子数
γ
100 101 102 103
¯ Γ
DNS
/
¯ Γ
DNS,γ=1.005
1 3 5
(b)
図–3 (a)式(5)に基づいて評価された衝突因子ΓDNSの時系列 分布．γ= 1.005の場合にΓDNSの値が最も低く，γが 増すにつれてΓDNSも増加．(b)衝突因子の時間平均値 Γ¯DNSの粒子比重γに対する依存性を，Γ¯DNS,γ=1.005(=
1.29)で規格化して表示．
i
j vp(i)
vp(j)
w=vp(j)-vp(i)
wr=w R/R xp(i)
xp(j)
O
図–4 動径方向粒子相対接近速度wrの定義．ベクトルRは それぞれの粒子の中心の位置ベクトルを用いてR ≡ x(i)p −x(j)p で表され，Rはその絶対値である．
密度の高まりの寄与が約4.86倍，と解釈できる．
(4) 考察
いわゆる流砂水理学のように，水中での土砂粒子の 振る舞いを扱う場合，その粒子比重は2.6程度である．
テイラー・レイノルズ数が42，衝突半径Rがコルモゴ ロフ長ηに等しい条件化で行われた本数値実験結果か らは，衝突因子ΓDNSはラグランジュ粒子の場合と比べ て約3.1倍（図–3(b)）となり，その内訳は粒子の衝突時 の動径方向相対接近速度⟨|wr|⟩が約1.1倍（図–6(a)），
図–5 (a)粒子の衝突時の動径方向相対接近速度⟨|wr|⟩，(b)動 径方向分布関数g(R)，(c)式(6)で推定された衝突因子 Γsph，の時系列分布．
動径方向分布関数g(R)が約2.9倍（図–6(b)）である．
空隙率Pvの値も約2.1倍（図–1(c)）となり，図–2にお いて視覚的にも明確なクラスタ形成が認められるケー スである．Nabiet al.21)は移動床流れの巨視的数値予測 における浮遊砂の取り扱いとして，質点近似に加えて，
粒子間相互作用と粒子混入による乱流変調を無視し，式 (2)と同様に簡素化されたMR方程式14)によって質点 粒子の挙動を追跡するアプローチを提案している．本 研究で得られた知見は，こういったアプローチに対し
γ
100 101 102 103
h|wr|i/h|wr|iγ=1.005 1.0 1.1 1.2 1.3
(a)
γ
100 101 102 103
¯g ( R )
1 2 3 4 5
(b)
γ
100 101 102 103
3
¯ Γ · τ /R
η 2 4 6(c)
Γ¯DNS
Γ¯sph
γ
100 101 102 103
¯ Γsp
h/
¯ ΓD
NS
0.99 1.00 1.01
図–6 (a)粒子の衝突時の動径方向相対接近速度⟨|wr|⟩，(b)動 径方向分布関数g(R)，(c)式(6)で推定された衝突因子 Γsph，の時間平均値の粒子比重γに対する依存性．(c) ではΓDNSも併せて示し，また，差し込み図には両者の 比Γ¯sph/Γ¯DNSを表示した．
てさらなる改良の方向性を示すものと言える．
他方で比重が1.1程度の水生微生物の場合には，そ の衝突頻度はラグランジュ粒子と比べて高々1.07倍に 留まる．粒子のクラスタ形成の程度を表す指標である 空隙率Pv，動径方向分布関数g(R)も共に1.07程度で あり，図–2において視覚的にクラスタの存在は認めら
れない．この分野では微生物の運動能力を，流体によ る移流の寄与と微生物が有する遊泳能力に分け，前者 については流体運動に完全追従するラグランジュ粒子 として取り扱われてきた5),6),7)．本研究結果はこのアプ ローチの有効性を支持するものである．
ただし，本研究では重力加速度の影響が無視されて いることに注意が必要である．2.(3)項に既に記された ように，前報では比重が2.0の場合に，Squires and Ya-
mazaki13)に倣って重力加速度の影響を加えた数通りの
条件下で検討したところ，いずれの場合にも終端速度 パラメータwが乱流速度変動強度を大きく上回った結 果，クラスタ形成が抑制されることが示された．よって 本報ではあくまでも乱流変動に起因する粒子の相対運 動にのみ着目して，粒子比重がクラスタ形成や衝突頻 度に及ぼす影響を議論した．他方で重力の存在が常に 粒子の偏分布や衝突頻度を抑制するわけではなく，そ の効果は乱流変動レベルや粒子パラメータに依って大 きく変化することが慣性粒子を対象とした研究で指摘 されている17)．乱流中では静水中よりも沈降速度が高 まること22)や，高粒子数密度条件下で乱流変調が生じ る場合には沈降がさらに促進されること23)も慣性粒子 に対しては報告されており，加えて混合粒径の場合に は，沈降速度の違いによってサイズの異なる粒子間の 衝突が促進されるなど，話はより複雑になる．これら の現象に対する粒子比重の影響について理解を深める には，任意の条件設定が可能となる数値実験の特長を 活かして，複雑化の要素をうまく切り離すアプローチ が有効と考えられ，今後の検討課題である．
最後に，前報2)では比重を100まで増加させたのに 対して，本報では比重が200, 500,1000のケースも検 討対象に加えた．その結果，粒子の偏分布の程度およ び衝突特性のいずれも，比重が100程度以上ではほぼ 不変であり，すなわち慣性粒子として式(4)によってそ の運動が記述できることが確かめられた．
4. おわりに
著者らは前報で2)，微小粒子の乱流中のクラスタ（偏 分布）形成に対する粒子比重γの影響を，直接数値シ ミュレーションによって検討した．本報ではその先行 研究を発展させ，クラスタ形成によって微小粒子の衝 突がどの程度促進されるかを検討した．
粒子が流体より高々数倍重い場合にもクラスタは明 瞭に形成され，100倍を超えると差異はほとんど認め られず，いわゆる慣性粒子（γ→ ∞）とみなせる．こ のとき粒子の衝突頻度ΓDNSは，流体運動に完全追従 するラグランジュ粒子とみなせるγ= 1.005の場合と 比べて，γ= 2,10,100,1000の場合にそれぞれ2.6倍，
t
0.0 0.5 1.0
hKi
0.8 1.6 1.2 1.4
(a)乱れエネルギー
t
0.0 0.5 1.0
hεi
0.4 0.6 0.8
(b)エネルギー散逸率
図–7 減衰一様等方乱流における乱れエネルギーとその散逸率の時系列分布．実線:本研究;プロット:参照DNS結果24)．
5.2倍，6.0倍，6.2倍に達する．式(6)の衝突因子モデ ルに基づけば，粒子の衝突促進は，衝突時の相対接近 速度の高まり（⟨|wr|⟩）の効果とクラスタ形成による局 所的な粒子数密度の増加（g(R)）の効果の2つに分け て考えることができる．数値実験結果から⟨|wr|⟩およ びg(R)を直接評価した結果，特に後者の寄与が前者を 数倍上回ることが確認された．また，式(6)は，ラグラ ンジュ粒子から慣性粒子まで広い条件下で，粒子の衝 突頻度を極めて良好に予測できることが示された．
水中の土砂粒子（γ≈2.6）の振る舞いを巨視的に扱 う場合，粒子のクラスタ形成やそれによる衝突促進の 影響が考慮されることはこれまでほどんど無かったも のの，本数値実験結果はその重要性を示唆するものと 言える．他方でγ ≈1.1の水生微生物に対しては，流 体運動に完全追従するラグランジュ粒子として近似的 に扱うことの有効性が支持された．
本数値実験は背景乱流場の乱流レベルや衝突半径，粒 子のストークス数をある値に固定して行われた．これ らのパラメータを系統的に変化させて，上述の知見の 普遍性を確かめることが望まれる．
付録 I 計算プログラムの精度検証
本数値実験で用いた計算プログラムの精度検証を目 的として，(i)一様等方乱流の減衰過程，および(ii)統 計的に定常な一様等方乱流中での微小慣性粒子の衝突 現象，に関する数値実験を行った．(i)は計算プログラ ム中の流体ソルバー部分の，(ii)は粒子挙動の追跡と衝 突判定に関する部分の，精度検証である．
(1) 一様等方乱流の減衰特性
本研究で用いた一様等方乱流ソルバーの精度評価と して，減衰性の一様等方乱流の数値シミュレーション 結果を，信頼性の高い参照DNSデータ24)と併せて図–7 と図–8に示す．643のフーリエ・モードを用いた擬ス ペクトル法でナヴィエ・ストークス方程式を解き，テ
イラー・レイノルズ数が約39から30まで減衰する過 程を再現した．乱れエネルギーおよびその散逸率の系 全体での空間平均値の時間変化率は参照DNSデータ24) と非常に良好に一致する．また，いくつかの時刻にお ける1次元パワースペクトルも参照データと良く一致 しており，特に高波数端領域において参照データが示 すデータの跳ね上がりは，本研究の計算結果では認め られない．これは本数値計算コードにおいてエイリア ス誤差の除去が良く機能しているためと理解できる．
(2) 統計的に定常な一様等方乱流中の微小慣性粒子の 衝突特性
流れ場はテイラー・レイノルズ数が約24で統計的に 定常な一様等方乱流場で，323のフーリエ・モードで速 度場を表現した．衝突半径Rが0.8∆x（ここで∆xは 計算格子サイズ）である微小慣性粒子1024個を初期に 流れ場にランダムに配置し，シミュレーション時間内 に生じた，スキームI11)に基づく粒子間の衝突回数を 記録した．Zhouet al.11)に倣って採用したこれらの条 件を本研究と，本研究とは完全に独立してプログラム 開発および計算実行がなされた渡邊25)で揃え，初期の 流速場も一致させた．ストークス数Sηを変化させなが ら数値実験を繰り返し，衝突頻度のストークス数依存 性を両者で比較した．計算プログラムに関する本研究 と渡邊の主要な相違点は，粒子位置で流体速度を求め る際に使用する補間法，および一様等方乱流を統計的 に定常に保つための低波数領域でのエネルギー注入法，
の2点のみである．渡邊の用いた数値計算法の詳細は Watanabe & Gotoh26)を参照されたい．
得られた衝突因子Γ≡Nc/(T V ρ2p/2)を図–9に示す．
ここでNc はシミュレーション中に観察された総粒子 衝突数，T はシミュレーション継続時間を示す．計算 プログラムに関する上述の相違点，およびシミュレー ション継続時間T や時間刻み幅∆tの相違を考慮すれ ば，両者から得られた衝突因子は極めて良好に一致し
k
100 101
1-D energy spectrum
10-4 10-3 10-2 10-1
t= 0.001
k
100 101
1-D energy spectrum
10-4 10-3 10-2 10-1
t= 0.500
k
100 101
1-D energy spectrum
10-4 10-3 10-2 10-1
t= 1.000
(a)乱れエネルギー
k
100 101
1-D dissipation spectrum
10-3 10-2 10-1
t= 0.001
k
100 101
1-D dissipation spectrum
10-3 10-2 10-1
t= 0.500
k
100 101
1-D dissipation spectrum
10-3 10-2 10-1
t= 1.000
(b)エネルギー散逸率
図–8 減衰一様等方乱流における乱れエネルギーとその散逸率の1次元パワースペクトル．実線:本研究;プロット:参照DNS 結果24)．
S
η≡ τ
p/τ
η0 2 4 6 8
Γ τ
η/R
30 2 4 6 8
present DNS Watanabe(2015)
図–9 一様等方乱流中の微小慣性粒子の衝突因子Γのストー クス数依存性．衝突因子Γはコルモゴロフ時間τηおよ び粒子衝突半径Rで規格化した．
ており，本研究で用いた計算プログラムの妥当性を支 持する．また，粒子緩和時間τpが増加してコルモゴロ フ時間スケールτηと同オーダーになるにつれて粒子衝 突が急激に促進される様子が認められる．
参考文献
1) Maxey, M. R.: The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields,J. Fluid Mech., Vol. 174, pp.441–465, 1987.
2) 横嶋 哲,藤井 秀太,宮原 高志:微小粒子の乱流中のクラ スタ特性に対する比重の影響,土木学会論文集A2, Vol.
72, No. 2, pp.I 459–I 465, 2016.
3) Jim´enez, J.: Oceanic turbulence at millimeter scales, Sci.
Mar., Vol. 61 (Supl. 1), pp.47–56, 1997.
4) Grabowski, W. W. and Wang, L.-P.: Diffusional and accre- tional growth of water drops in a rising adiabatic parcel: ef- fects of the turbulent collision kernel,Atmos. Chem. Phys., Vol. 9, pp.2335–2353, 2009.
5) Lewis, D. M. and Pedley, T. J.: The influence of turbulence on plankton predation strategies, J. Theor. Biol., Vol. 210, pp.347–365, 2001.
6) 馬場 謙二郎,横嶋 哲,益子 岳史,宮原 高志:乱流中のプ ランクトン間の遭遇率と摂食率,土木学会論文集B1, Vol.
69, No. 4, pp.I 817–I 822, 2013.
7) 横嶋 哲,安田 昌平,宮原 高志:乱流環境中を遊泳する水 生微生物の捕食頻度のモデリング,土木学会論文集A2, Vol. 71, No. 2, pp.I 713–I 718, 2015.
8) Cuthbertson, A. J. S., Dong, P. and Davies, P. A.: Non- equilibrium flocculation characteristics of fine-grained sed- iments in grid-generated turbulent flow,Coastal Engineer- ing, Vol. 57, pp.447–460, 2010.
9) Wang, L.-P., Wexler, A. S. and Zhou, Y.: On the collision rate of small particles in isotropic turbulence. I. Zero-inertia case,Phys. Fluids, Vol. 10, pp.266–276, 1998.
10) 横嶋 哲,益子 岳史,松坂 隆弘,宮原 高志:一様等方乱流中 のラグランジュ粒子の接触特性評価,京都大学数理解析 研究所講究録1822, pp.76–83, 2013.
11) Zhou, Y., Wexler, A. S. and Wang, L.-P.: On the collision rate of small particles in isotropic turbulence. II. Finite iner- tia case,Phys. Fluids, Vol. 10, pp.1206–1216, 1998.
12) Wang, L.-P., Wexler, A. S. and Zhou, Y.: Statistical mechan- ical description and modelling of turbulent collision of iner- tial particles,J. Fluid Mech., Vol. 415, pp.117–153, 2000.
13) Squires, K. D. and Yamazaki, H.: Preferential concentration of marine particles in isotropic turbulence,Deep Sea Res., Vol. 42, No. 11/12, pp.1989–2004, 1995.
14) Maxey, M. R. and Riley, J. J.: Equation of motion for a small rigid sphere in a nonuniform flow,Phys. Fluids, Vol.
26, pp.883–889, 1983.
15) Squires, K. D. and Eaton, K.: Preferential concentration of particles by turbulence, Phys. Fluids A, Vol. 3, pp.1169–
1178, 1991.
16) Yoshimoto, H. and Goto, S.: Self-similar clustering of in- ertial particles in homogeneous turbulence,J. Fluid Mech., Vol. 577, pp.275–286, 2007.
17) Onishi, R., Takahashi, K. and Komori, S.: Influence of grav- ity on collisions of monodispersed droplets in homogeneous isotropic turbulence,Phys. Fluids, Vol. 21, 125108, 2009.
18) Maxey, M. R.: The motion of small spherical particles in a cellular flow field,Phys. Fluids, Vol. 30, pp.1915–1928,
EFFECTS OF MASS DENSITY RATIO ON COLLISION STATISTICS OF SMALL PARTICLES IN A TURBULENT ENVIRONMENT
Satoshi YOKOJIMA and Yoshiaki SHIMADA
This is the first attempt to study the effects of mass density ratio on collision statistics in an isotropic particle-laden turbulent suspension. Previous studies have concentrated on the collisions of either the La- grangian particles following the fluid motions completely or the inertial particles which is much heavier than the surrounding fluid. The present study is an extension of Yokojima et al. (J. JSCE A2 72(2) I 459 2016), where the influence of mass density ratio on clustering/preferential concentration of monodisperse small particles in a turbulent environment was examined closely by direct simulations of turbulent flows.
It has been revealed from the present study that the collision frequency is increased with increase in the specific density, i.e., the ratio of the mass density of the particles on the fluid density, and that the primary factor is the local/instantaneous increase in particle number density owing to the particle clustering. The accuracy of the collision-frequency estimation model based on the radial relative velocity and the radial distribution function is also examined. The model prediction is found to be in excellent agreement with the simulation results.
1987.
19) 杉山 和靖,高木 周,松本 洋一郎:渦セル流れ中の気泡・
粒子の並進運動,日本機械学会論文集（B編）, 69巻680 号, pp.786–793, 2003.
20) Fallon, T. and Rogers, B.: Turbulence-induced preferential concentration of solid particle in microgravity conditions, Exp. Fluids, Vol. 33, pp.233–241, 2002.
21) Nabi, M., de Vriend, H. J., Mosselman, E., Sloff, C. J. and Shimizu, Y.: Detailed simulation of morphodynamics: 2.
Sediment pickup, transport, and deposition,Water Resour.
Res., Vol. 49, pp.4775–4791, 2013.
22) Wang, L.-P. and Maxey, M. R.: Settling velocity and con- centration distribution of heavy particles in homogeneous isotropic turbulence, J. Fluid Mech., Vol. 256, pp.27–68, 1993.
23) Bosse, T., Kleiser, L. and Meiburg, E.: Small particles in homogeneous turbulence: Settling velocity enhancement by two-way coupling,Phys. Fluids, Vol. 18, 027102, 2006.
24) de Bruyn Kops, S.:私信, 2004.
25) 渡邊 威:私信, 2015.
26) Watanabe, T. and Gotoh, T.: Hybrid Eulerian-Lagrangian simulations for polymer-turbulence interactions, J. Fluid Mech., Vol. 717, pp.535–575, 2013.
(2017. 5. 8受付)
