微小液滴の乱流衝突成長に対する Large-Eddy Simulation(乱流と輸送現象:コーヒーカップから宇宙まで)

微小液滴の乱流衝突成長に対する Large-Eddy

Simulation

地球シミュレータセンター 大西 領 (RyoOnishi) ,

高橋桂子 (Keiko Takahashi)\star 1

京都大学大学院工学研究科機械理工学専攻 松田景吾 (KeigoMatsuda),

小森 悟 (Satoru Komori) \star 2

\star 1

_{Earth Simulator}

Center, Japan Agencyfor

Marine-Earth

Science

andTechnology

\star 2Department

ofMechanical Engineering andScience, Kyoto University

1

緒言

気相乱流中で微小液滴粒子が互いに衝突するという現象は,

蒸気発生器, 噴霧器, サイクロン分離 器などの工業装置内[山1の流れにおいてだけでなく, 鉛直方向に発達する対流雲中においても見られる. 従って, その衝突挙動, 衝突成長過程を正確に予測することは, 工業装置の最適設計・操作, および 対流雲とそれに伴う豪雨の高精度予測に関連して, 機械工学的, 環境工学的に極めて重要である. 本研究では特に, 環境工学分野での重要性に注目した. 日本では1990年代以降, 1 時間あたりの雨 量が 100mmを超える局地的豪雨による被害が増加傾向にある. 特に都市部での局地的集中豪雨は徒大 な被害をもたらしている. 例えば,

2000

年

9

月に発生した東海豪雨では愛知県内だけでも約

8650

億 円もの被害が出たと言われている [$8$ 】. このため, 局地的豪雨の発生機構の解明と発生予測が急務の課題 とされている. _{この課題を解決するためには, 豪雨をもたらす対流雲の発生・発達をできるだけ正確} に予測できる数値シミュレーション手法を開発する必要がある

.

対流雲中では, 強い上昇流に伴って 発生する乱流のために雲粒の衝突成長が促進されることが指摘されている$M\mathfrak{l}\epsilon 116$ ]. そのため, 冒頭で述 べた乱流中での微小液滴粒子の衝突現象を正確に考慮できる数値シミュレーション手法の開発が, 対 流雲の発生発達予測. つまりは豪雨予測の高精度化に欠かせな$Aa$

.

既に,

液滴径までを考慮した数値シミュレーション手法は開発されており

,

液滴の大きさまでを考

慮して衝突成長を計算することが可能である.

しかし, その数値シミュレーションでは液滴粒子の衝

突に及ぼす乱流の影響が十分に考慮されていない

.

例えば,

Geresdi

and

$Ra\epsilon mu\epsilon 8en$ (2004)[剛による

山岳降雨の計算では, 乱流の影響が全く考慮されていない衝突頻度モデル

\iota 8]

が使用されている

.

また,

$Iffnn$

et al.

$(2006)^{[f]}$_{によるスコールラインの計算では, 乱流効果係数 (定数) を導入することによっ}

て大雑把に乱流の影響を考慮しただけである.

以上の状況を考慮し, 本研究では, 微小液滴 (雲粒) _{の乱流衝突成長を正確に考慮できる数値シミ} $=$レーション手法を確立することを目的とした. 具体的には, _{$Large\cdot Eddy$}

Simulation

(LES)_による数 値シミ$=$ レーション手法を開発することを目的とした.

第 2 節では, 粒子の乱流衝突頻度を予測するために必要な乱流衝突頻度因子モデル開発について記

述し, 第3節では, _{LES による乱流衝突成長計算手法の開発について記述する}. 第 4 節では, 気相乱

流中での液滴の乱流衝突成長に対して_{, 液滴運動をラグランジュ法によって追跡する直接数値計算}

(Direct

_Numerical

Simulation,DNS)と本LESの両方を適用し, その結果を比較することによって本

LES

手法の信頼性を確認した.

2.

乱流衝突頻度因子モデル

2

$\cdot$

1

粒子の衝突成長式

粒予の運動および成長をオイラー法によって計算する際には, 衝突による粒子数密度関数$n_{p}(r,x_{*},t)$ の時間変化率を次式により計算する. $\frac{\partial n(r,l)}{\partial\prime}=\int\kappa_{\iota}\langle r’’,r’$)

$n_{p}(r’’,t)n_{p}\langle r’,t$)$dr’-f^{K\langle r,r’)n_{p}(r,t)n\{r’,t)\theta’}$ (1)

ここで, r”=(P-rs)l\hslash である. $K_{\iota}\langle r_{l}$,

r\acute )

は衝突頻度因子と呼ばれ

,

半径

$r_{1}$ の粒子と半径りの粒子が衝突す

$N.(r_{1},r_{2})=K.(r_{1}.r_{2})N_{p}(r_{1})N_{p}(r_{2})$ (2) この衝突頻度因子モデルの代表的なものとして, 次式で表される

Hydrodynamic

Kernelモデルがある 18]. $\langle\kappa_{c.h\rangle^{*}}.\rangle=R_{12}^{2}|V_{p.\infty}(r_{1})-l_{p.\circ}’\langle r_{2}f$ (3) ここで, R,,は衝突半径と呼ばれ2粒子の半径の和(rl+r2)を表し, $V_{p}$.g\langle )は半径$r$の粒子の終末速度を表 す. 式(3) 中には, 気相の流れ場の状態を表す量 (例えば, レイノルズ数) が含まれていないことから も明らかなように, この Hydrodynamic

Kernel

では衝突に及ぼす気相乱流の影響は全く考慮されてい ない. 気相乱流の影響を考慮した衝突頻度因子モデルとして, 我々は既に同一径粒子間の乱流衝突頻度因 子モデルを開発した [9]. しかし, 異径粒子間の乱流衝突頻度因子に関しては, 信頼性の高いモデルは提 案されていなかった. そこで, 本研究では, まず異径粒子間の衝突頻度因子モデルの開発を行った. 本節では, そのモデルに関して記述する.

2

$\cdot$

2

重力による沈降を考慮しない揚合の異径粒子間乱流衝突頻度因子モデル

重力による粒子の沈降を考慮しない場合には, _{乱流の効果に起因する衝突頻度因子 (Kc.\omega lb)} は次のよ うに表される$\mathfrak{l}10$]. $\langle K_{c.1w4}\rangle=2rR_{\mathfrak{l}2}^{2}(|w_{r}|\rangle g_{11}(R_{12})$ (4) ここで, 添字turbは乱流の効果により生じた項であることを表しており, $<|\iota\nu_{r.1_{-\sqrt{}>}}$ およびg’2(R) はそれ

ぞれ, 半径$r_{1}$ の粒子と半径$r_{2}$ の粒子との接触時相対接近速度 (radial

relative

vclocity atcontact) および近 接粒子分布関数(radial

distribution function

$al\bm{m}n$_論ct) である. なお, 以降では$w_{r,12}$ を$w_{r}$ と略記する.

\leq wA>は 2 粒子が接触

(衝突) する際の相対速度の絶対値の時間平均値である

.

gl2(R) は粒子が不均一

に分布することによる衝突促進効果を表す.

Zhou

et

al.

$(2001)^{[11]}$_は,

Williams

&Crane(1983)2] および

Kruis&Kusters(1997)1\iota 3]

の解析モデルを改良して

<|wA>

に関する次のモデルを提案した

.

$\langle|w_{r}|\rangle_{\backslash ub}=(\frac{2}{\pi}\langle w_{r}^{2}$

)

$)^{1\prime 2}$

(5)

$=[ \frac{2}{\pi}((|V_{r.\infty}||\rangle+(|1l_{r.. n}’|\rangle)]^{1\prime l}$

ここで,$w_{r..ocd}$ は接触時相対接近速度のうち2粒子間相対速度の変動強度に起因する成分を表し,$u_{r..h-}$

は流体の局所速度勾配に起因する成分を表す. これら 2 成分は, それぞれ次のようにモデル化された.

$(|w_{\infty\epsilon d}^{2}|)=C_{*}.\langle\ell$)$u^{\prime 2} \frac{r}{\gamma-1}\{(\theta_{1}+\theta_{2})-\frac{\theta_{1}\theta_{2}}{\theta_{1}+\theta_{l}}[\frac{1+\theta_{1}+\theta}{(1+\theta_{l}X]+\theta_{2})}]^{\iota/1}\}[\frac{11}{([+\theta_{I}K1+\theta_{2})(1+w_{1}X^{[+\gamma\theta_{2})}}]$ (6)

$\Downarrow w\grave{i}m|\rangle=\frac{1}{15}(\frac{R}{l_{n}})^{2}\lambda^{2}$ (7)

式 (6) において, $\theta_{i}$は積分スケール$l,$, 流体の速度変動のRMS値u

および粒子の緩和時間\mbox{\boldmath $\tau$}p.’ を用いて,

$\theta$

,

$=2.5\tau_{p}.’\cdot u’/l_{l}$ のように定義され, $\phi=\max(\theta_{2}/\theta_{1}, \theta_{1}/\theta_{-}),$ $\gamma=\phi\cross 0.183u^{2}/(ev)^{1\prime 2}$ である. また, 右辺第一項

$C,.(\phi)$ _{は経験パラメータであり},

Zhou

at

al.

(2001)[11] は次のように定義した.

$C_{l’}(\phi)=1.0+0.6\exp[-(\phi-1)^{1.5}]$ (8)

中の異径粒子間の近接粒子分布関数

(g12(R)turb)

は以下のようにモデル化された[1}].

$g_{12}(R_{12})=1+\rho_{I2}[g_{I1}(R_{II})-1]^{I\prime 2}k(R,,)-1]^{1/}-$ (9)

ここで, $\rho_{12}\text{は_{}\psi}\tau_{p.1}/\tau_{p,2},$$\tau_{p},/\tau_{p.1}$) を用いて次のようにモデル化された.

$\rho_{12}=2.7\exp(-\psi\}+0.205\exp\{\triangleleft.0206\psi)^{x}0.5[1+\tanh(\psi-3)]$ (10)

なお,

Zhou

$\epsilon t$

al. (2001)111]

の論文中においては

,

右辺第一項の係数は 26 である.

しかし, $\psi\approx 1$ のとき に

\mbox{\boldmath$\rho$}12

$=1$ _{であるべきことを考慮すると, 27の方がふさわしい. 同一径粒子間の近接粒子分布関数}

{

$gdR\rangle_{urb}$) に関しては, 既に我々の過去の研究 [9] において,_{信頼性の高いモデルが開発されているので}, 以上により$g_{12}(R)_{turb}$

.

さらには衝突頻度因子を求めることができる

.

2

$\cdot$

3

重力による沈降を考慮する揚合の異径粒子間乱流衝突頻度因子モデル

前節では, 式 (4) が議論された. しかし, 式(4)では, 式(3)で表される, 2粒子間の終末速度差に起因 する衝突は考慮されていない. 現実の流れ場では重力により, 粒子は沈降するので, 沈降速度が速い 粒予が遅い粒子を捕集するという現象が起こる

.

この衝突を表現したのが式(3)に示された Hydrodynamic Kemel (Kchvdr) である. _{乱流だけでなくこの重力の効果までを考慮した全衝突頻度因子}

($K_{c}$.bUr)は,

重力の効果までを考慮した接触時相対接近速度

(<|w4>t

。

d

および近接粒子分布関数

\omega I2(R\mbox{\boldmath $\lambda$}r\infty 1)

を用いて式 (4)と同様に次のように表される.

$K_{c,\iota\alpha 1}|(r_{1}, r_{2})=2\pi R^{2}<|wt>|\alpha d$ $g_{1_{-}}\cdot(R\lambda_{ad}$ (11)

式中の$<|_{1}\nu\sqrt{}>_{t\propto\lrcorner}$ に関しては, w, の確率分布がガウス分布に従うと仮定すると

$<|\#\cdot\sqrt{}>_{1\alpha 1}l=(<|w\lambda>_{\iota urb^{2}}+<|w\sqrt{}>bydr^{\sim})^{l\prime 2}$ (12)

という関係が成り立つので, 式(5)および$<|\nu 4>hydr=|t_{p.\infty}’\{r_{\iota}$) $-V_{p.\infty}(r_{2})|/2$ の関係を用いて算出することが できる. しかし, g’2(R\a’に関しては全く明らかになっていない. そのため, $\mathcal{K}_{1A}$を厳密に予測する

ことは現時点ではできない. このような状況の中, Falkovichetal.

{

$2002)^{[l]}$_は$k_{t\circ I-}’$を_{$K_{1}$}と$K_{hyi1}$の単純

な和で表現した.

$K_{c.1\propto l1}\approx t_{turb}+k_{hy4r}’$ (13)

しかし, 上式を用いると, $g_{1_{-}},(R_{12})\approx 1$ の場合, $<|w’\sqrt{}> ad\approx<|w\lambda>_{t\iota rb}+\triangleleft w\sqrt{}>hydr$となる. これは式 (12)と整

合しない. この問題を解決するために, 本研究では次式を提案する [14].

$K_{c.\mathfrak{i}oIat}\approx(K_{c.t\ovalbox{\tt\small REJECT}}+K_{bdr})^{1/2}$ (14)

本式を用いると, $g_{1_{-}}(R_{1_{-}})\approx 1$ の場合に式(12)を満足する.

2

$\cdot$

4

直接数値計算

(D-irect

_Numerical

Simulation,

DNS)

_{による衝突頻度因子の算出}

モデルの予測結果を評価するために用いる衝突頻度因予に関するデータを

DNS によって求めた. 気

相乱流場の運動をフーリエ_{. ガラーキン法を用いた擬スペクトル法により計算した. 支配方程式は次} 式で示される渦度方程式である.

$\frac{\partial\omega,}{\partial t}+U_{j}\frac{\partial\omega_{J}}{a_{j}}=\omega_{\text{ノ}}\frac{\partial U_{l}}{a_{\text{ノ}}}+V\frac{\partial_{\mathcal{O}1_{j}}^{2}}{a_{j}a_{J}}+F_{;}$ (15)

渦度の時間発展には

4

次精度のルンゲクッタ法を用

$Aa$

) また, 3分の2法を用いてエイリアジング誤差

て, 計算領域中に $u\cross 64\cross 64$ 個の計算格子数を等間隔に配置した. 側面境界条件として周期境界条件 を用い, また, 波数が1である渦の強さが一定になるように強制外力$F_{j}$を与えることによって, 定常等 方性乱流場を形成させた. なお, 本数値計算法を用いることにより, 理想的な定常等方性乱流場が形 成されることは, 我々の過去の研究で既に確認されている[16].

DNS

においては粒子運動をラグランジュ法を用いて計算した. このとき, 粒子の運動方程式は次式 で表される.

$\frac{\partial a)}{\partial t}+U_{j}\frac{\partial a\gamma_{J}}{b_{j}}=\omega_{j}\frac{\partial U,}{a_{J}}+\nu\frac{\partial^{2}\omega_{/}}{a_{j}a_{j}}+F_{t}$ (16)

重力による粒予の沈降を考慮する場合には右辺中の$g$を98$[m/s^{-}]$に設定して計算を行った. 雲粒を想定

して, 微小水滴を混入粒子とした. 上述の定常等方性乱流場の中に, 半径$r_{1}$および r\tilde ’ を持った粒子をそ

れぞれ2048個ずつ, 合計4096個混入させた. 2 粒子間の距離が衝突半径以下になれば, 衝突したとみ

なした. 衝突による運動量交換によって引き起こされる余分な衝突を避けるために, 衝突後は2粒子 のうち一方を消滅させた. 雲粒の代表的な数密度$n$は$n\approx 10^{1}$

[l/m3]

であり

,

これと代表的な雲粒半径$r\approx$

$10[\mu m]$より, 粒子体積率\phi、は$\phi_{v}\approx 10^{\prec}\sim 10^{-5}$, _{粒了質量率 \phi ,は}$\phi_{n}\approx 10^{-2}\sim 10^{-1}$程度であると見積もら

れる. これは粒子濃度が十分低いことを意味しているので, 本研究では粒子による乱流変調を無視し

た. 衝突頻度因子$K_{c}$の算出法$n$時間ステップ目において, 半径$r_{1}$の粒子と半径$r_{2}$の粒子が衝突した回数

を$N_{\infty 1.12}^{\hslash}$ とする. そのときの衝突頻度(Nc’$\prime 01\sim 2$) は次式のように求まる.

$N_{\epsilon.12}^{n}= \frac{N_{Y}^{n_{V.t2}}}{t_{d}’\Delta t}$ (17)

ここで, $V_{d}$は計算領域の体積を表し, 捜は計算時間間隔を表す. このとき, $n$時間ステップ目の衝突頻

度因子K訟は$K_{c.1_{\vee}}^{n}$

.

$=N_{\epsilon}’’./n.n\iota z_{P^{1}’\cdot-}$から求まる. 時間ステップ毎に求めた K 翫を平均することによって, 時

間平均した衝突頻度因子K,2を求めた. $\frac{Table1T.u.ltu1entstatistics.ofairi1owinDNS}{\frac{u’ l,I10^{3}m1l,[10^{4}m1l_{\lambda}f10^{5}m1}{0.1409043644.64}}$

2

$\cdot$

6 衝突頻度因子最に関する結果と考察

表1に, DNS から得られた気相乱流場に関する乱流特性値を示す. 表1に示した乱流特性値を持っ た気相乱流場中での, 重力による粒子の沈降を考慮しなかった場合のK\epsilon に関する本モデルの予測結果 および DNSの結果を図1に示す. なお, 一方の粒子半径を$r_{1}=30[\mu m]$として, もう一方の粒子半径$r_{2}$を グラフの横軸にとった. また,

縦軸は認で割ることによって無次元化された邸を表している

.

図 1 よ り, $r_{2}$が $25[\mu m]$よりも小さい領域では本研究で提案したモデルの結果とDNSの結果の間に差具が見ら れるものの, その差異は最大でも 3 割程度であり, 本モデルの結果はDNSの結果と良好に一致すると 言える. 続いて, 重力による粒子の沈降を考慮した場合の瓦の結果を図 2 に示す. 比較のために, 乱流の効

果を考慮しない衝突頻度因子モデルである

Kt.hyd’

の結果も同じグラフに示されている

.

図2より,

DNS

の結果は$K_{c,hydr}$の結果よりも大きい. これは, 乱流によって衝突が促進されることを意味する. また, 式 (13) を用いた予測縞果は, r-,が大きくなり, 濃度ムラの効果が小さくなる領域, つまり$g_{12}(R\succ 1$ とな る領域でDNSの結果よりも大きな予測となっている. これに対して, 式 (14) を用いたモデルの予測結果 は全領域でDNSの結果とよく一致する. つまり, 重力の効果を考慮する場合には, 本研究で提案した 式(14)を用いることが妥当であることがわかる.

Fig.

1 Comparison ofDNS

and

model

predictions$ofK_{c}$

in no-gravity

case.

$(r_{1}\approx 30\mu m)$

Fig.

2

Comparison

ofDNS

andmodel

Predictions

$ofK_{c}$ingravity

case.

$(r_{l}=30\mu m)$

3.

粒子の衝突成長に対する

_Large-Eddy

Simulation(LES)

3

$\cdot$

1

気相乱流場の計算方法

気相乱流場を計算する際の支配方程式は連続の式およびN-S 方程式である. これらにフィルタ操作 を施すことにより

LES

での支配方程式を得ることができる. このようにして得られた支配方程式を以 下に示す. $\frac{\partial\overline{U_{i}}}{a,}=0$ (18) $\frac{\partial\overline{U_{i}}}{\partial\prime}+\overline{U_{\dot{\text{ノ}}}}\frac{\partial\overline{U_{i}}}{a_{J}}=-\frac{1}{\rho}$$\frac{a\overline{P}}{e_{l}}+(\nu+\nu_{SCS})\frac{\partial^{2}\overline{U_{j}}}{a_{J}a_{j}}$ (19)

419)中のvsos はサブグリッドスケール (sub grid scale, SGS) 渦粘性係数であり, 静的スマゴリンスキー

$\nu_{scs}=(C_{\backslash }.\Delta)^{2}\ulcorner D|$ (20) ここで, \Delta は計算格子幅であり, $C_{s}$はスマンゴリンスキー定数である. 本研究では, 等方性乱流場にお ける理論値である$C_{s}=0.173$ を用いた. また, $|-D|$は計算格子スケール (grid scale, GS)のひずみ速度テン ソルの大きさであり,

ひずみ速度テンソルー

D,j

$= \frac{1}{2}(\frac{\alpha\prime_{/}-}{\alpha_{j}}+\frac{\delta\overline{l\prime_{j}}}{a_{l}})$ を用いて次のように表される. $\pi=\sqrt{2\overline{D_{j}D_{j}}}$ (21)

8

$\cdot$

2

粒子成長の計算方法

気相乱流中での粒子運動の支配方程式は粒子の質量密度関数$\sqrt{}$r) の輸送方程式である. これにフィ ルタ操作を施すと次式が得られる.

$\frac{\partial\overline{q_{P}(r)}}{\ }+\overline{U_{J}}\frac{\partial\overline{q_{p}(r)}}{a_{j}}=(\frac{\partial\overline{q_{p}(r)}}{\partial\iota})_{ca\prime}+(\frac{\partial\overline{q_{p}(r)}}{\partial t})_{\epsilon a\prime d}+V_{p,\infty}\frac{\partial\overline{q_{p}(r)}}{\partial x_{2}}+K_{m}\frac{\partial^{2}\overline{q_{p}(r)}}{a_{\text{ノ}}a_{J}}$ (22)

式(22)の右辺第1項および第2項は粒子質量密度の衝突による増減および凝縮・蒸発による増減をそれ

ぞれ表している. 本研究では衝突成長を主な対象としたので, 凝縮蒸発の項を考慮しなかった. こ

れは凝縮成長の速さに比べて衝突成長の速さが非常に大きい場合を対象としていることを意味する. 式 (22) 中の衝突成長項は, 粒子数密度関数の時間変化率と次の関係がある.

$( \frac{\partial\overline{q_{\rho}(r)}}{\partial t})_{\infty l}=m(r)(\frac{\partial\overline{n_{p}(r)}}{\partial\prime})_{co1}$ (23)

この式と式 (2)を用$Aa$ さらに衝突頻度因子$K_{c}$に対してモデルを与えることによって, 粒子の衝突成長

を計算することができる. また, 式(22)中の右辺最終項に現れる陥は$q_{p}(r)$のSGS渦拡散係数であり, 本 研究では

\sim

$=v_{scs}$と仮定した$[17][18]$

.

2 節において開発した衝突頻度因子モデルは, 入力として気流の速

度変動のRMS値$u’$, 積分スケールみ,

_{テイラーマイクロスケール}

1\mbox{\boldmath $\lambda$}

_{およびコルモゴロフスケール}

l\eta

_を必要

とする. これらのうち, みおよびu’は流体の大スケール運動に支配されるので, $l_{l}\approx\overline{l,}$および$u’\approx\overline{u’}$ と考

え, 流体のGS成分から算出した. —方,

_{l\mbox{\boldmath $\lambda$}}

_および

_T,,

_{は流体の小スケール運動に支配される}

_.

まず,

_らを

求めるために, 粘性消散率\epsilon を次式のようにGS成分(句s)とSGS成分(\epsilon s$\propto$)で近似した.

$\epsilon\approx\epsilon_{CS}+\epsilon_{SGS}$

$= \nu\frac{\partial\overline{U,}}{\partial x_{j}}\frac{\partial\overline{U_{(}}}{\ }+(C_{s} \Delta)^{2}ff^{3}$

(24)

なお, 実際の気象計算では格子解像度が低いためv $\epsilon_{GS}<<\epsilon_{S0S}$と予想されるが, 本研究では\epsilon os も無視す

ることなく考慮した. 式 (24) から得られた\epsilon を$l_{\eta}=(v^{3}/e)^{1\prime 4}$に代入することにより,

-を算出した

.

さらに,

\epsilon および u’ から$l_{\lambda}=u’\sqrt{15\nu/\epsilon}$の関係を用いて,

l\mbox{\boldmath $\lambda$}を算出した.

_っまり, _{局所等方性を仮定することによ}

って, $l_{\eta}$と

T\mbox{\boldmath $\lambda$}

を算出した

.

オイラー法を用いたLES によって粒子成長を計算するためには, 粒子半径を離散化しなければならな

$\{\begin{array}{l}r_{1}=0[\mu \mathfrak{n}]r_{k}=r_{2}\cdot C^{k-2}\end{array}$

(25)

ここで, $r_{2}=5[\mu m],$ $C=1.03$ として, クラス 1からクラス

100

($r_{1\infty}=90$6[\mbox{\boldmath $\mu$}m]) までを計算対象とした.

4.

乱流中での粒子の衝突成長に対する数値計算

本節では, 4節で開発した

LES

手法を用いて乱流中での粒子の衝突成長を計算した. 同様の条件で

DNS も行い,

LES

の結果と DNSの結果を比較することにより, 本LES手法の信頼性を検証した.

$4\cdot 1$

LES

の計算条件

代表長さ (k) を$0.005[m]$として, 計算領域の大きさを$2_{y}\cdot cA$ x2d0x2y-c恥に設定した. そして, 計算領域

中に $16xl6xl6$_{個の計算格子を等間隔に配置した}

.

壁面境界条件として周期境界条件を用$A\backslash$ また, 波 数が

1

である渦の強さが一定になるように強制外カ$F_{l}$を与えることによって, 定常等方性乱流場を形成 させた. フィルタ関数としてはスペクトルカットオフフィルタを用$Aa$ _{カットオフ波数を}8_とした. 各計算格子点における粒予径の初期分布として

,

図 3 中に太実線で示される分布を用いた. この初期 分布は, 次のようにして得られた. まず,

SeeBelberg

et

al. (19%)[19]

が用いた凝縮核の数密度分布を持っ た粒子群を, その平均半径が$30[\mu m]$になるまで凝縮成長させた. _次に, _{計算領域内の総粒了数が荀%} 個になるように, 得られた粒了径分布に係数を掛け合わせた. なお, 半径 $30[\mu m]$_{という大きさは対流} 雲中で見られる代表的な雲粒の大きさである. 以上の計算条件の下で, 4節で説明したLESを実行した. なお, 衝突頻度因子モデルとして,

2 節で開発したモデルを用いた場合と Hydrodynamic

Kcmel モデルを 用いた場合の両方で計算を行った. $4\cdot 2$

DNS

の計算条件

LESと同じ大きさの計算領域に対して, 24 節で示した方法を用いて, 気相乱流場および粒子運動を 計算した. ただし, 定常な流れ場を形成させるために強制外力を加える際には, 波数が1である渦度 の大きさが前節のLESにおける大きさと同じになるように設定した.衝突成長を計算する前準備として,

異なった径を持った 4096 個の粒子を計算領域内に等間隔に配置した.

この4096個の粒子の粒子径分 布は,

LES

で用いた初期粒予径分布と同一になるように設定した

.

計算開始後しばらくの間, 弾性衝突 を仮定して粒予位置の時間発展を計算した

.

十分に時間が経過し, 粒子分布に及ぼす初期粒子位置の 影響が無視できると考えられる時間 (本研究では, 積分時間スケール$T_{e}(=l7u)$の$3$倍を目安にした) 経 過した時点を$t-T$ _{とした. そして}, $t>0$ _{では衝突による粒予成長を計算した}. その際, 衝突した 2 粒子 は直ちに合体して1つの粒子になると仮定して, 質量保存と運動煮保存の関係を満たすように衝突後 の粒子の半径速度を決めた. $4\cdot 3$

LES

と

DNS

の比較結果および考察

表 2 に, LESおよび_{DNS から得られた気相乱流場に関する乱流特性量を示す}. 両者の結果の差異は 最大でも1割程度であることから, 本LES手法を用いることにより乱流特性量を精度良く評価するこ とができることがわかる. 図 3 に, $t/T_{e}=40$_{における質量密度関数の結果を示す}. _{実線で示されたDNSの結果に注目すると,}

明瞭な粒子分布の極大がいくつか見られる

.

例えば,「$-31$[\mbox{\boldmath $\mu$}m]付近に見られる極大は\sim 3O[\mbox{\boldmath $\mu$}m]の粒予

同士が衝突してできた粒子を表しており,「$-43$ [\mbox{\boldmath $\mu$}m]付近に見られる極大は「-30 [\mbox{\boldmath $\mu$}m]の粒子と「$-3l[\mu m]$

の粒子が衝突してできた粒子を表している

.

図3 より, 衝突頻度因子として,

_{Hydrodynamic Kemel}

(Kc.hydr) を用いた LESの結果では,

DNS

の結果に比べて,「$-38$[pm]の粒子が極端に少ない. これは, $K_{hyb}$

を用いると, _{粒子の成長速度をかなり過小評価することを意味する}

.

つまり, 乱流中での粒子の衝突

成長を予測する場合には, $A_{by}$,を用いることは不適切である. 一方, 本研究で開発した衝突頻度因子

モデルを用いた LES の予測結果はDNSの結果と非常に良く一致しており, 本LES手法が乱流中での粒子

$V^{-}\Xi$

嘗

$\wedge*_{\theta}^{\llcorner}$

$r[\mu m]$

Fig.

3

Predictions

ofthe

drop

mass

densityfunction at$t\Pi e=40$

.

5.

結

$z$

本研究では, 乱流中での粒子の衝突成長を正確に予測可能な

Large-Eddy

Simulation(LES)手法を開発

することを目的とした. そのために, まず 2 粒子間の衝突頻度因予モデルを開発した. その際, 乱流

によって引き起こされる衝突と, 粒子の終末速度差によって引き起こされる衝突の両方を考慮する方

法として式 (14) を提案した. そして, 開発した衝突頻度因子モデルを LES に組み込んで乱流中での粒

子の衝突成長を実際に計算した. その結果を直接数値計算(Direct

Numerical

Simulation, DNS)による予測

結果と比較することにより, 本LES 手法の儒頼性を確認した. 一方, 従来のメソスケールの気象計算

で用いられてきたHydrodynamic Kernel を用いた

LES

は, 対流雲中での粒子の衝突成長を正しく予測で

きないことを明らかにした. 本研究で開発した, 衝突頻度因子モデルおよびLES手法は, 機械工学の 分野だけでなく, 環境工学の分野にも応用されることが期待される.

謝辞

本研究は京都大学大学院高木一氏 (現, (株) 神戸製鋼所) の協力を得た. 記して深く謝意を表する. また, 本数値計算の一部は独立行政法人国立環境研究所地球環境研究センターのスーパコンピュータ $(NEC:SX\cdot 6)$を用いて行われた. センターの関係者各位に深く感謝する. 本研究は科学技術振興機構

(JST)

戦略的創造研究推進機徹CREST)マルチスケール マルチフィジックス現象の統合シミュレーション領 域「災害予測シミュレーションの高度化プロジェクト」 の援助を受けて行われた.

参考文献

[1] Williams,J.J.E.

and

Cmne, R.l., Dmp Coagulation

in Cross-Over

Pipe Flows

of Wet

Steam,

Journal

_of

[2] Abrahamson, J., Collision Rates of Small

Particles

in

a

Vigorously

Turbulent

Fluid,

Chemical

EngineeringScience,Vol.30(1975),

pp. 1371-1379.

[3]

_{Foundation ofRiver&Basin}

Integrated

Communications

ed., PORTAL,Vol.15(2002),Pp.2-8.

[4] Falkovich, $G$ et al.,

Acceleration of Rain Initiation

by Cloud Turbulence, Nature, Vol.419 _(2002), pp.151-154.

$\int 5]$ Lynn, B.H., Khain, A.P,, Dudhia, J., Rosenfeld, D., Pokrovsky A.

and

Seifert, A.. “Spectral (Bin)

Microphysics Coupled with

a

Mesoscale Model(MM5).Pan II:

SimuIation

ofa CaPE Rain Event with

a

SquallLine, $\Lambda tonthly$WeatherReview,Vol. 133(2005),Pp.59-71.

[6] Houze, R.A.Jr. $a\mathfrak{n}d$Medina, S.,

Turbulence

as a

_{$M\epsilon chanism$} for

$ot\mathfrak{v}graphiC$

Precipitation

Enhancement,

Journal

_of

the

_AtmosPheric

Sciences,Vol.

62

(2005),

pp. 3599-3623.

[7] Geresdi,

1. and

Rasmussen, R., Freezing Drizle

Formation in

Stably

Sffalifid

$Lay\epsilon r$

CloMs.

Pan II:

Role

of

Giant

Nuclei and Aerosol Panicle Size

$Disn\cdot ibution$

and Solubility,

Journal

of

$tk$

Amospkric

Sciences,

Vol.

62

(2004),

pp.2037-2057.

[81 Hall, W.D., $A\propto tailedM\ddagger c\iota ophysical$

Model Within

a

$\mathfrak{m}0- Dim\epsilon nsio\mathfrak{n}d$ Dynamic

Framework:

Model

$\propto scription$

and

$Pr\epsilon li\min_{\Psi}$Results,

Journal

ofthe

AlmosphericSciences,

Vo1.37

$(19lO)$

.pp.2486-2507.

[91 Onishi,R. andKomori, S.,

Modeling of the

collision

kemel

of$monodisp\epsilon rsepar\dot{\mathfrak{n}}cl\epsilon s$

in

turbulentflows,

Transactions

_ofthe

$J\varphi an$Sociely

ofMechanical

$Engineers_{1}$Series$B$.(in print)

[10] Wang, L.-P. et al.,

Statistical Mechanical

Description and Modelling of Turbulent

ColIision

of Inenial Particles,

Journa’

_of

Fluid

Afechanics,$Vot.415$\langle 2000),

pp. 117-153.

[11] Zhou,

Y.

et al., Modelling$Turbul\epsilon nt$

Collision of

$Bidisp\epsilon oee$lnertialParucles,Journal

ofFluu

Mecknics,

Vo1.433

$(2W1),$Pp.77-104.

[12] Williams, J.J.E. and Crane, R.I., _{Pmicle Collision Rate}

_in

_Turbulenl Flow, Imemational Joumal

of

MultiPhase

$Fl\alpha\nu$

.

Vo1.9,No.4(1983),Pp.421-435.

[13] Kruis,F.E. andKusters,K.A.,$lhe$_Collision$Rat\epsilon$

of

$Puticl\epsilon s$

in

$Turbul\mathfrak{n}\epsilon t$Flow,

Ckmical

Enginceling

Science,$VoI.158(1997),$$pp.201\cdot 230$

.

[141 $0\mathfrak{n}iShi$, R.,

Takahashi.

K. and

Komori, S., _{$Laoe\triangleright Bddy$}

Simulation

for Particle

Collision

_{$c_{IOW}m$}

in

TurbulentFlows, Transactions

_of

theノa\mbox{\boldmath$\mu$}n Society

_of

MechanicalEngineers. Series$B$, Vo1.72,No.722

(2006).

_pp.244

$1\cdot 2448$

.

[15] lshioka,K.,

2005:

isPack-0.71,$ht\phi;//www.gfd- dennou.oryarch/ispacW$,

_GFD

$\propto nnou$

Club.

[16] Onishi,

R. and

Komori, S.,

The

gravity $eff\infty ts$

on

the col[ision _{$\Re qucncy$}

of

$monodisp\epsilon\iota sePaHicleS$

in

turbulent

flows

(in Japanese), Transactions

_of

the $Ja\mu n$ Socie’y

of

Mechanical

$Engi’\kappa ers$,

&ries

$B$,

$Vol.71$,No.705\langle 2005),

pp. 1225-1230.

[17] Kogan, $Y.L.$,The $Simulat\ddagger on$of

aConvective

Cloud

in

a3-D Model with Explicit$Mic\iota ophysics$

.

Pan1;

Model Description and Sensitivity $Exp\epsilon riments$, Journal

of

the Atmo\varphi &ric

Sciences.

Vo1.48

(1990).

$pp.1160- 1189$

.

[18] $Koga\mathfrak{n},$ $Y.L$

.

$\epsilon t$ al., Modeling of Sratocumulus Cloud Layers in aLaIge Eddy Simulalion Model with

Explicit Microphysics, Joumal of the

_ArnosPheric

$sCie\mathfrak{n}CCs$

.Vo1.52

(1995).pp.2923-2940.

[191 $Se\epsilon O\epsilon 1berg_{*}$M. et

al..

$St\propto hastic$

Simulations

as

aBenchmark

for Muhcmatical

$M\epsilon thods$

Solv

$ing$ the

Coalescence Equation,$\Lambda$tmospheric_{$Resea\kappa h$},Vo1.40(1996),Pp.33-48.

[20] Pinsky, M. et al.,

Scochastic

Effects of Cloud Droplet Hydrodynamic Interaction in

a

$Turbule\mathfrak{n}t$ Flow,