定常的な粒子輸送下での粒子-粒子層間の衝突過程
田邊章洋 1 , 新屋啓文 2 , 粟津暁紀 1 , 西森拓 1
1
広島大学大学院 理学研究科数理分子生命理学専攻
2
名古屋大学大学院 環境学研究科
概要
砂床や雪面において、風によって駆動される地表面付近の堆積粒子輸送は、粒子-流体相互作用、
粒子-粒子相互作用によって維持される。粒子輸送の発生初期では、風によって粒子が空気中に
取り込まれるが、十分に発達した粒子輸送においては、堆積粒子の取り込みは主に跳躍粒子と粒
子層の衝突によって行われることが知られている。そのため、跳躍粒子と粒子層の衝突過程を調
べることは、風に駆動される粒子輸送を理解する上で非常に重要である。本研究では、離散要素
法を用いて 1 粒子を粒子層への打ち込む計算機実験を行った。そして、得られた結果を既存の風
洞実験の結果で重み付けすることで、定常的な粒子輸送を定性的に再現した。
Grain-Bed Collision Process under Stationary Grain Transport
Takahiro Tanabe1 , Hirofumi Niiya 2 , Akinori Awazu1 , Hiraku Nishimori1
1
Department of Mathematical and Life Science, Graduate School of Science, Hiroshima University
2
Graduate School of Environmental Studies, Nagoya University
Abstract
On sand and snow fields, wind-induced stationary transport of grains near the bed surface
is sustained under complicated processes: grain-fluid interaction, grains midair collision, and
grain-bed collision. At first, grain entrainment into the air from the bed is caused by grainfluid interaction, but soon, a major part of grain entrainment shifts to be caused by grain-bed
collision processes. Therefore, grain-bed collision is considered as the most relevant issue in
studying the sand transport over the bed. In this paper, we simulate 1 grain collision with
bed using DEM (Discrete Element Method). Thereafter, we attempt to reproduce the grain
movements in the stationary grain transport by superimposing the numerical result at each
injection angle with the weight of the stationary distribution of injection angle obtained in
wind tunnel experiment.
1
はじめに
表面の凹凸などによって粒子が空中に飛び出すと、
運動モードは転がりから跳躍へ変化する。そして、
風による砂や雪粒子の運動形態は、粒子のサイズ
や風速に依存して転がり (Creep)、跳躍 (Saltation)、
浮遊 (Suspension) の3つに分類されることが知られ
跳躍粒子は、流体から運動量を受け取り、その後風
下の地表面に衝突する。その衝突によって、地表面
を構成する粒子は空気中に取り込まれる。この過程
ている [1]。堆積している粒子は、ある風速よりも強
はスプラッシュ過程と呼ばれる。
い風が吹く事で、運動を開始する。最初、運動を開
初期の粒子輸送は風による粒子の取り込みで発生
始した粒子は地表面に沿って転がるが、いったん地
するが、粒子輸送が十分に発達すると、風による粒
1
子の取り込みは極端に少なくなる。それに代わって、 で記述される。さらに、接線方向の粒子の衝突に関
スプラッシュ過程が粒子の取り込みの主要になり、 して、クーロンの摩擦則を実現するため、法線・接
定常的な粒子輸送は維持される [2]。そのため、スプ
線方向の力の比が摩擦係数 µ を超えた時に、滑りを
ラッシュ過程を調べることは、地表面付近 (特に跳
表現する摩擦スライダーを考慮する。具体的に、斥
躍が活発に起きている跳躍層) における粒子輸送メ
力は以下の形で計算される。
(
)(
)
cos θi,j − sin θi,j
Fri,j
i,j
n
Fr =
sin θi,j
cos θi,j
Fri,j
t
カニズムを理解する上で重要である。
これまでに、飛砂粒子の軌道や速度は高速度カメ
ラを用いた風洞実験などで観測されてきた [3] が、
{
地表面付近の粒子の運動は空間に占める粒子数が多
Fri,j
n
いため、測定が困難であった。ところが、近年の技
術の進歩により、レーザーを用いた位相差から地表
面付近の粒子の運動形態がとらえられつつある [4]。
また、粒子-粒子層間の衝突に対する実験的・理論的
Fri,j
t
=
(|r i − r j | >
−kn δni,j − ηn δ̇ni,j

 0
=
0
di +dj
2 )
(otherwise)
(|r i − r j | >
 min{−kt δti,j − ηt δ̇ti,j , µF
研究も行われており、放出粒子の入射角度・速さの
i,j
i,j δt
rn |δ i,j |
t
(2)
(3)
di +dj
2 )
} (otherwise)
ここで、n、t は粒子間の法線方向と接線方向、θi,j
統計的性質が議論されてきた [5–7]。
は粒子 i、j の中心を結ぶ直線と x 軸 がなす角を表
しかしながら、多くの先行研究で使用された粒子
し、sin θi,j = (zj − zi )/|r i − r j |、cos θi,j = (xj −
は、実際の砂粒子と比べ非常に大きい、または、単
xi )/|r i − r j | である。また、di は粒子 i の粒径、δni,j
一粒径である。そこで、本研究では、粒径に幅を持
つ粒子層を形成し、スプラッシュ過程について調べ
、δti,j は法線、接線方向に対する粒子 i、j の接触開
た。具体的には、粒子層に対して 1 粒子を打ち込むこ
始時からの変位の積算を表す。そして、k 、η は法線、
とで、粒子-粒子層間の衝突シミュレーションを行っ 接線方向の弾性、散逸に関するパラメータである。
4
2
2
2
た。そして、風洞実験で得られた粒子の特性 (入射 今回は kn = 10 [kg/s ]、kt = 3 × 10 [kg/s ]、ηn =
−3
−3
角度分布) を粒子の衝突過程に重み付けすることで、 4.1 × 10 [kg/s]、ηt = 2.0 × 10 [kg/s]、µ = 0.3 と
定常的な粒子輸送の様子を 1 粒子衝突の重ね合わせ した。
3
で表すことを試みた。
2
モデル
結果
上記モデルを用いて、我々は 2 次元鉛直断面で衝突
我々は 離 散 要 素 法 (DEM: Discrete Element
シミュレーションを行った。シミュレーション空間の
Method) を用いて、衝突過程の計算を行う。仮定と
底面 (粒子層の底面) に粒子を固定し、底面粒子との
して、粒子は球形であり、その形状は変化しないと
衝突は法線方向の斥力のみが考慮される。ここで、底
する。次に、粒子の回転運動は無視され、水平・鉛
面粒子は堆積粒子と比較して小さく (df = 100µm)、
直方向の 2 次元平面に沿った並進運動が考慮される。 衝突による衝撃が底面で反射して粒子層表面に影響
−2
最後に、全粒子の運動は Newton の運動方程式に従 を与えないよう、やわらかい (ηf = 1.0×10 [kg/s])
い、粒子に働く力として、重力と粒子間の接触時に
ものとする。また、水平方向の境界は周期境界とし
働く斥力を考える (式 (1))。
た。そして、粒子数は 2048 個で、各粒子の密度は
mi r̈ i = F i = mi gez +
∑
j
F i,j
r
ρp = 2500[kg/m3 ] で固定される。ただし、粒子の大
(1)
きさは単一ではなく、粒径 160µm、200µm、240µm、
280µm、320µm の粒子をそれぞれ 204、409、822、
409、204 個用いた。粒子のパラメータは砂を仮定し
量、位置、鉛直下向きの単位ベクトルを表す。また、 ており、これらの粒子のランダムな初期位置からの
g 、F i,j
r は重力加速度と粒子 i、j の接触による斥力 自由落下により、初期の粒子層は形成される。
を表す。
形成された粒子層に対して粒子を入射するが、入
ここで、i、j は粒子の番号、m、r 、ez は粒子の質
粒子間の斥力は 2 粒子の法線・接線方向にそれぞ
れ分けて考える。そして、これらの斥力は線形の粘
射粒子の大きさは粒子層の平均粒径 d¯ = 240µm を
¯ とする。こ
与え、初期位置は (x0 , z0 ) = (R, zb + d)
弾性で表現できると仮定し、Kelvin-Vogit モデルを
こで、R は [0, L] の一様乱数 (L=2.5cm:システムサ
用いる。つまり、衝突は反発と摩擦力の重ね合わせ
イズ)、zb は粒子層の最高点を表す。今回は、入射
2
(4)
図 1: 衝突過程の様子
入射粒子 : 黒、放出粒子 : 灰色 (入射粒子のリバ
ウンドも灰色)、破線は打ち込み粒子の初期高さ z0
速さ vi を 1.5、2、3、4m/s と固定し、入射角度を
5◦ から 90 °まで 5 °刻みで変化させ、放出粒子の角
度と速さ分布について調べた。ただし、入射速さ・
図 2: 放出角度分布 (vi = 2m/s): (青) 数値計算の結果
角度の各組み合わせ毎に粒子の初期位置 x0 を変え、
入射速さ
A
B
α
µ
σ
100 回のアンサンブルをとり、また、放出粒子は入
1.5m/s
2m/s
3m/s
0.55
0.83
2.2
14.6
17.3
25.8
0.096
0.089
0.087
ln 0.76
ln 0.64
ln 0.55
0.88
0.85
0.78
4m/s
2.6
28.6
0.087
ln 0.51
0.78
射粒子の初期の高さ z0 に到達した粒子と定義する
(図 1)。
数値計算で得られた放出角度と速さ分布を LiQiang
らの風洞実験の結果 [4] と比較した。風洞実験におい
表 1: 各入射速さに対する、放出角度 (式 (5))・速
さ (式 (6)) の Fitting パラメータ
て、粒子層は砂で構成されており、粒子層表面へ入射
する飛砂粒子および、衝突によって粒子層表面から
放出される粒子それぞれの角度と速さが測定された。
そして、入射・放出それぞれの角度と速さ分布に関す
る Fitting 関数が提案された。その結果として、入射・
放出に関わらず、角度分布は指数分布、速さ分布は対
する。入射の速さを変化させた場合も同様の結果が
得られた (表 1)。
次に、放出粒子の速さ分布について、風洞実験の
数正規分布を用いて近似された。しかしながら、[4]
結果は放出粒子の速さの平均値で無次元化されてい
では入射の角度と速さの相関については言及されて
る。従って、数値計算の結果もそれに合わせて、放
いない。従って、ここでは速さ一定の入射粒子が、[4]
( θi )
1
で得られた入射角度分布 f (θi ) = 7.4
に
exp − 36.3
出粒子の平均の速さで無次元化した。風洞実験の速
さ分布は対数正規分布
従って砂床に衝突してくると仮定し、定常状態にお
(
)
(ln ve∗ − µ)2
α
exp −
f (ve∗ ) = √
2σ 2
2πσve∗
ける粒子輸送と対応付けた。すなわち、数値計算の
結果をこの入射角度分布で重み付けすることで、放
出角度分布と放出速さ分布を見積もる。
(6)
でよく Fitting されている (α = 0.129、µ = ln 0.919
まず、放出角度分布について、風洞実験の放出角
、σ = 0.668)。入射の速さを 2m/s の場合、数値計算
度分布は指数分布
の結果は α = 0.089、µ = ln 0.64、σ = 0.85 でよく
(
)
θe
1
f (θe ) = exp −
A
B
Fitting された (図 3)。また、他の関数 (Maxwell 分
(5) 布、Gamma 分布、Levy 分布) を用いて Fitting を
試みたが、対数正規分布での Fitting が最も誤差を
でよく Fitting されている (A = 8.2、B = 40.6)。
小さく与えることがわかった。以上の結果から、速
入射の速さが 2m/s の場合の、数値計算の結果は一
さ分布に関しても、我々の計算結果は風洞実験と定
見すると指数分布に従っていない。ところが、放出
性的に一致する。入射の速さを変えた場合でも同様
角度 35 °以上の範囲に限定すると、指数的に減少
の結果が得られた (表 1)。
している (図 2)。35 °以下の浅い角度で放出される
最後に、放出角度と速さの相関について着目した
粒子を無視すると、Fitting パラメータは A = 0.83
(図 4)。図 4(a)、(b) 共に、低速の粒子は幅広い範囲
、B = 17.3 となり、風洞実験の結果と定性的に一致
に放出角度を持つのに対して、高速の粒子は狭い範
3
図 3: 放出速さ分布 (vi = 2m/s): ve∗ は平均放出
速さ (v̄e = 0.35m/s) で無次元化された速さ、(青)
数値計算の結果、(赤)Fitting 関数
囲に放出角度が集中している。しかしながら、(a)、
(b) を比較すると、高速粒子の放出角度に違いが生
じている。つまり、高速粒子の放出角度を考えると、
入射粒子のリバウンドは比較的浅い角度 (5 °∼30 °)
で多く発生する (図 4(a))。一方で、入射粒子自身で
なく、衝突により粒子層から飛び出す粒子は、比較
的深い角度 (30 °∼50 °) で飛び出すことが分かった
(図 4(b))。
4
図 4: 放出の角度と速さの関係 (vi = 2m/s): (a) 入
射粒子のリバウンドのみを考慮 (v̄e = 0.61m/s)、
(b) 粒子層から放出された粒子のみを考慮 (v̄e =
0.21m/s)
子の角度と速さは独立ではないと考えられる。従っ
まとめと今後の展望
て、その両者の相関を含んだ重みを考慮した上で、
今回、我々は堆積粒子の風による輸送過程のうち
放出角度と速さ分布は議論されるべきであり、これ
特に粒子衝突過程に注目し、離散要素法を用いて粒
は今後の課題である。
子-粒子層間の衝突に関する数値実験を行った。その
参考文献
結果、1 粒子の衝突過程の重ね合わせで定常的な粒
子輸送を定性的に再現した。具体的には、各入射の
[1] R. A. Bagnold, The Physics of Blown Sand
速さに対して風洞実験で得られた入射角度分布を重
and Desert Dunes. London: Methuen, 1941.
み付けし、放出角度と速さ分布を算出した。そして、
これらの分布は風洞実験で提案された Fitting 関数
[2] R. S. Anderson, and P. K. Haff, Science, 241,
820-823, 1988.
と定性的に一致した。
[3] Z. XueYong et al., Geomorphology, 36, 155165, 2001.
しかしながら、放出角度分布について、35 °以下
の浅い角度は風洞実験に基づく Fitting 関数から大
[4] K. LiQiang, G. LieJin, and L. DaYou, Sci
China Ser G-Phys Mech Astron, 51, 986-1000,
きく外れてしまう。その原因として、シミュレーショ
ンでは入射粒子の初期高さ z0 を超えた粒子のみを
2008.
放出粒子と定義しているためだと考えられる。つま
[5] B. T. Werner, and P. K. Haff, Sedimentology,
り、地表面を転がるような浅い角度で放出された粒
35, 189-196, 1988.
子は分布に含まれていない。また、表 1 に示すよう
[6] F. Rioual, A. Valance, and D. Bideau, Phys.
Rev. E, 62, 2450, 2000.
に、Fitting パラメータは風洞実験の結果と定量的に
一致しない。この理由として、我々の数値計算では
[7] M. Xing, and C. He, Geomorphology, 187, 94100, 2013.
風洞実験よりも速い粒子を入射していることが挙げ
られる。最後になるが、実際の系において、入射粒
4
ダウンロード

粒子層間の衝突過程 - Physics of Traffic