論

J. Soc. Powder Technol., Japan, 61, 258–267 (2024)
doi: 10.4164/sptj.61.258

文

Original Paper

湿式ボールミル中の媒体ボール衝突速度および角度が砕料粒子粉砕挙動におよぼす影響の解析
久志本 築*，加納 純也

Analysis for Effects of Collision Velocity and Angle of Grinding Balls
on Particle Grinding Behavior during Wet Ball Milling
Kizuku Kushimoto* and Junya Kano
Received 22 December 2023; Accepted 16 February 2024

The simulation using the Advanced Discrete Element Method-Computational Fluid Dynamics (ADEM-CFD)model analyzed the motion and breakage behavior of particles during wet ball milling in order to investigate the effectsof the collision velocity and angle between two grinding balls on the particle grinding behavior. The analyses showedthat the particles were more damaged when the grinding balls collided with a faster collision velocity in closer to anormal collision direction. The degree of the particle damage was quantified by the release energy, which was definedwith the potential energy released by breaking the internal bonds of the particles. The release energy was strongly
related to the normal components of the impact energy calculated from the collision velocity between two grinding
balls.
Keywords: Wet ball milling, ADEM-CFD model, Simulation, Collision angle, Collision velocity.

1. 緒

言

らは，湿式遊星ボールミルにおける媒体ボールの材質および直径が珪砂の粉砕におよぼす影響について実験的に

湿式ボールミルは，砕料粉体が懸濁するスラリー中で解析し，ボール径とボール密度が小さくなると粉砕平衡

媒体ボールを衝突させることで砕料粒子を微細化する粉径も小さくなることを報告した[4]。また，横田らは，湿砕方法の一種であり，粉砕速度が高速かつ大規模化もで式転動ボールミルによるアルミナ粉体の粉砕速度におよき，砕料粒子を数 μm 程度の粒子径まで微細化可能であぼす媒体ボール充填量とポット回転数の影響を実験的にることが特徴である。そのため，セラミックス[1]，鉱業

解析し，粉砕速度が最大となる回転数は媒体ボール充填[2]，医薬品[3]など幅広い産業製品の原料，中間製品，最量により異なることを報告した[5]。一方で，実際の湿式

終製品の製造に用いられている。

ボールミル中の粉砕を考えると，媒体ボールの衝突速度

しかしながら，湿式ボールミルは，媒体ボール，砕料

や衝突角度といった衝突状況も粉砕に影響すると考えら

粒子，分散媒や装置条件など粉砕に影響をおよぼす因子れる。しかしながら，媒体ボールの衝突状況が粉砕におが多く，粉砕の制御が難しいことが課題である。とりわ

よぼす影響に関する報告は見当たらず，これは，特定のけ，湿式ボールミルにおける媒体ボールは，粉砕を直接速度，特定の角度で媒体ボールを衝突させたときの砕料担う物体であるため，その物性や特性あるいは衝突状況

粒子の粉砕挙動を実験から解析することが難しいためとが粉砕におよぼす影響を把握することは，湿式ボールミ考えられる。このような理由から，媒体ボールの衝突状

ルによる粉砕を制御する上で必要である。

況が粉砕におよぼす影響については十分に理解できてい媒体ボールの材質や直径，充填量などの条件が，湿式ないのが現状である。

ボールミルによる粉砕におよぼす影響については，これこうした実験では解析困難な現象に対し有効な手法とまでにも数多くの実験的な報告がある。たとえば，横山

してシミュレーションによる解析がある。シミュレーショ
ンはコンピュータ上で実際の現象を模擬する手法であり，東北大学 多元物質科学研究所（〒 980-8577 宮城県仙台市青葉区片平 2-1-1）
Institute of Multidisciplinary Research for Advanced Materials,
Tohoku University（2-1-1 Katahira, Aoba-ku, Sendai, Miyagi 980-8577, Japan）* Corresponding Author kizuku.kushimoto.d2@tohoku.ac.jp

258

ボールミル内の媒体ボール，分散媒，砕料粒子の挙動を可視化できるだけでなく，媒体ボール間や砕料粒子に働く瞬間的な力や積算したエネルギーを抽出することもできる。特に，各構成要素に作用する合力から並進と回転の運動方程式をその要素ごとに立て逐次的に解くことで粉体工学会誌構成要素群の運動を記述する Discrete Element Method

Method（FDM）を用いる。B) 媒体ボールの運動は Discrete（DEM）と，流体の運動を記述する連続の式と運動方程Element Method（DEM）[9]により追跡し，C) 砕料粒子の式を近似的に解き解析する Computational Fluid Dynamics

運動および破壊挙動は Advanced DEM（ADEM）[11]を適（CFD）を連成させた手法は湿式ボールミル内の媒体ボー用する。なお，D) 砕料粒子と媒体ボールの相互作用はル挙動の解析に頻繁に用いられている。たとえば，Beinert

DEM により計算する。E) 分散媒と媒体ボールの相互作

らは DEM と CFD の連成シミュレーションにより，湿式用は体積力型埋め込み境界法[14,15]により求め，F) 分散

ボールミル中の媒体ボール挙動を表現し，媒体ボール同媒と砕料粒子の相互作用は DEM-CFD モデル[8]をベース

士の衝突点で生じる応力エネルギーを Normal impact，に推算する。以下では，媒体ボールの運動，砕料粒子のTorsion，Shearing と Rolling の 4 つに分類する手法を開発

運動と破壊，分散媒の運動の順に，シミュレーションモ

し，それら指標により粉砕機の種類，運転状態，媒体ボーデルの詳細を説明する。

ルの衝突状況が特徴づけられることを報告した[6]。ま2.1 媒体ボールの運動

た，Carvalho らは，DEM と CFD の連成シミュレーショ

媒体ボールの運動は DEM[9]により追跡する。DEM でンにより，連続式の縦型湿式ボールミル中の媒体ボール

は，媒体ボール 1 個や砕料粒子 1 粒といった各構成要素の運動を解析し，ミルの上部と下部で媒体ボールの衝突に作用する力を求め，並進と回転の運動方程式を要素ごの傾向が異なることを見出した[7]。このようにシミュとに立て逐次的に解くことで構成要素群の運動を記述す

レーションを用いることで実験では取得困難な情報を得る。今回の場合，媒体ボールに作用する力は，媒体ボーることができる。

ル同士の接触力，媒体ボールと砕料粒子の接触力と，媒一方，既存のシミュレーションの多くは媒体ボールの体ボールが分散媒から受ける流体力を考える。媒体ボー

運動のみを表現しているため，今回対象としている砕料

ルに作用する流体力の推算に体積力型埋め込み境界法

粒子の粉砕挙動については，直接的な解析は難しい。こ[14,15]を適用する。体積力型埋め込み境界法では，空間

うした中，著者らは，スラリー中で媒体ボールに捕獲さを区切る計算格子（流体セル）における，固体状物質のれる砕料粒子の運動および破壊挙動を表現するシミュ

体積割合を求め，その体積割合に応じて，流体に強制力

レ ー シ ョ ン モ デ ル と し て ， Advanced Discrete Elementを作用させ，その反作用力を固体状物質に与えることで，

Method-Computational Fluid Dynamics（ADEM-CFD）モデ

流体と固体状物質の相互作用を計算する手法である。固ルを開発し，さらに，液中で接近する媒体ボール周りの体状物質を媒体ボールへ，流体を分散媒へ置き換え，媒

砕料粒子の運動と破壊挙動について実験と比較すること

体ボールが分散媒から受ける流体力の計算に体積力型埋で，ADEM-CFD モデルの妥当性も確認した[8]。したがっめ込み境界法を適用すると，媒体ボールの並進および回て，ADEM-CFD モデルを用いることで，これまでほとん

転の運動方程式はそれぞれ次式で与えられる。

ど解明されてこなかった媒体ボールの衝突状況が砕料粒
子の粉砕挙動におよぼす影響を解析できると考えた。
そこで本研究では，ADEM-CFD モデルを適用したシ
ミュレーションにより，スラリー中で衝突する 2 つの媒
体ボール間の砕料粒子挙動を表現し，このときの媒体ボールの衝突速度や衝突角度といった衝突状況をパラメータにすることで，媒体ボールの衝突状況が砕料粒子の粉砕
挙動におよぼす影響を解析する。
2. 媒体ボール間の砕料粒子挙動を表現するシミュ
レーション方法
湿式ボールミルでは分散媒，媒体ボール，砕料粒子が相互作用しながら共存している。ADEM-CFD モデルでは，湿式ボールミル中の媒体ボール間で粉砕される砕料
粒子挙動をシミュレーションにより表現するために，次の 6 つの運動と相互作用をモデル化する。
A) 分散媒の運動
B) 媒体ボールの運動
C) 砕料粒子の運動および破壊挙動
D) 砕料粒子と媒体ボールの相互作用
E) 分散媒と媒体ボールの相互作用
F) 分散媒と砕料粒子の相互作用
A) 分 散 媒 の 運 動 は Computational Fluid Dynamics（CFD）により表現し，差分近似には Finite Difference

Vol. 61 No. 5（2024）ki
d mball, i uball, i= ∑ Finj + Fit j − ρf
dt
V
j=1

f ball dV

(1)

ball

ki
d(I ball, i · ωball, i)= ∑ rball, i × Fit j + Rr − ρf
dt
V
j=1

rf × f ball dV
ball

(2)
ここで，添え字i, jは着目する要素の番号とその近傍に存在する要素の番号をそれぞれ表し，k は接触している要素数を表す。なお，ここでの要素は媒体ボールと後述する
砕料粒子を構成する構成粒子のことを意味する。mballは

媒体ボールの質量， I ballは媒体ボールの慣性モーメントテンソル，uballは媒体ボールの並進速度ベクトル，ωballは

媒体ボールの角速度ベクトル，F nとF tはそれぞれ法線方
向と接線方向の接触力，Rrは転がり抵抗，rball, iは媒体ボール中心から接触点に向かう相対位置ベクトル，rf は媒体

ボール中心から分散媒の計算格子（流体セル）中心に向かう相対位置ベクトル，V ballは媒体ボールの体積， ρf は分散媒の密度， f ballは媒体ボールが分散媒から受ける流

体力をそれぞれ表す。なお，媒体ボールの衝突状況が砕料粒子の粉砕挙動におよぼす影響を抽出するために，重力の影響は無視している。

法線方向の接触力F nはバネとダッシュポットを並列接続した Voigt モデルを適用し，接線方向F tにはバネとダッ

259

シュポットに加えスライダーを追加したモデルを適用す解くことで，その集合体である砕料粒子の運動を記述することで摩擦の影響も考慮している。このようにモデルる。このように砕料粒子を扱うと，構成粒子間の連結バ化すると，F nとF tはそれぞれ次式で与えられる。
F n = K n Ln + ηn dL
dt

ネの変形により砕料粒子の変形が表現できるとともに，

n

(3)

t
dφ
F t = min μ F n , K t Lt + rball φ + ηt dL + rball
dt
dt

Lt
Lt

作用力とそのほか外力から運動方程式を時々刻々解くこ

ここで，添え字n，tはそれぞれ法線と接線を表す。K はバ
ネ定数，ηはダッシュポットの粘性係数，Lは並進変位ベ
クトル，φは角変位ベクトル，μは摩擦係数をそれぞれ表
す。なお，法線方向のバネ定数 K n は媒体ボールのヤン
グ率Eとポアソン比νから Hertz の弾性接触理論[10]に基づき与え，接線方向のバネ定数K tはラメ定数の定義式とK nから与えると，K nとK tはそれぞれ次式で計算できる。
−1

n
K = K
2 1 + νit

うに設定することで砕料粒子の破壊も表現できるようになる。したがって，構成粒子個々に働く連結バネからの(4)

rir j 1 − νi2 1 − ν2j+
Kn = 4E jπ
3π ri + r j Eiπ

この連結バネがある閾値以上に変形した時に破断するよ(5)とで砕料粒子の運動，変形，破壊が表現される。今回の場合，外力には分散媒からの流体力と接触力が働くことを考慮すると，構成粒子 1 個の運動方程式は次式で与えられる。ここで，構成粒子は砕料粒子を構成する一部であることから構成粒子自身の回転は無視できるものとし
ている。
ki
si
d mp, i up, i
S, t
S, n
+ Fim
= ∑ Finj + Fit j − ∑ Fim+ Ff p, i (11)
dt
j=1
m=1

ここで，siは連結バネで接続された構成粒子数，mは着目する構成粒子と連結バネで接続された構成粒子の番号を意味する。F S, nとF S, tはそれぞれ連結バネから受ける力の(6)

さらに，粘性係数ηはバネ定数を元に臨界減衰となるよう

法線成分と接線成分を表し，それぞれ次式で与えられる。
F S, n = KS,n lS,n +ηS,n ΔlS,n
Δt

(12)

F S, t = KS, t lS, t +ηS, t, ΔlS, t
Δt

(13)

に次式で与える。
η = 2 mK(7)

また，転がり抵抗 Rrは次式で与えられる[11,12]。

ここで，KSは連結バネ定数，ηSはダッシュポットの粘性(8)係数，ΔlSは連結バネの変位ベクトルを表す。連結バネ定

ここで，μrは転がり抵抗係数，bは接触面の半径をそれぞ縮試験から得られる応力ひずみ線図の勾配が合うようにり，統一されたモデルは確立されていないため，本論文

る値を持つはずであるが，それらバネ定数の比率についφRr = −μrb F
φnれ表す。転がり抵抗を表現するモデルは複数存在しておでは湿式ボールミル中の媒体ボール挙動を表現した実績[12]のあるモデルを採用した。
媒体ボールと分散媒の相互作用 f ballは，体積力型埋め込み境界法[14,15]に基づき次式から計算する。
f ball = −α vball − u fΔtvball = uball + ωball × rball

数KSは砕料粒子の弾性に相当する係数であり，単粒子圧決定される。連結バネ定数の法線成分と接線成分は異なては決定方法が確立されていない。一方で，接線成分の連結バネ定数を法線成分の連結バネ定数の 1/1000 程度まで小さくすると安定した計算が実施できることが経験的にわかっているため，今回は接線方向の連結バネ KS, tの法線方向の連結バネ定数 KS, n に対する比（KS, t /KS, n）を(9)

1/1000 に固定することにした。なお，ダッシュポットは(10)

合，ダッシュポットを考慮しなくても十分に安定した計

ここで，αは着目する流体セル中で媒体ボールが占める
体積割合，vballは分散媒の計算格子（流体セル）の中心における媒体ボールの速度ベクトル，u f は分散媒の速度ベ

クトル，Δtはタイムステップをそれぞれ表す。

媒体ボールの運動は，Eqs. (3)～(10) から求められる媒体ボールへの作用力とモーメントをそれぞれ Eq. (1) とEq. (2) へ代入し解き，媒体ボールの位置と速度を逐次更新することで追跡する。
2.2 砕料粒子の運動および破壊挙動
砕料粒子の運動および破壊挙動は ADEM[11]により表計算の安定化のために導入されたものであり，今回の場算が実施できたため無視している。
連結バネの破断は，連結バネの自然長 l0に対する法線方向の伸びが構成粒子径に対しある一定のひずみεmaxに到達するときに起こるものとして，次式で定義する。
lS, n − l0 /dp > εmax(14)

ここで，l0は自然長であり初期配置における 2 つの構成

粒子中心間距離から求め，dpは構成粒子径を表す。εmaxは砕料粒子の脆性に相当し，単粒子圧縮試験から得られる
応力ひずみ曲線の破断位置が合うように決定される。

現する。ADEM では，1 個の砕料粒子を複数の微細な構媒体ボールに作用する流体力については体積力型埋め成粒子が連結バネと呼ばれる仮想的なバネで結合した集込み境界法を用いて推算した。しかしながら，この方法合体であるとみなし，構成粒子ごとの運動方程式を逐次を ADEM で表現される砕料粒子にそのまま適用しようと

260

粉体工学会誌すると，構成粒子よりも十分に小さい流体セルが必要とめり込み量δにある任意の閾値δmaxを与え，この閾値を超なるため，計算負荷が高く現実的な時間で計算を完了す

えた接触が生じた DEM 粒子を ADEM 粒子に変換し，一ることが難しい。一方，ADEM-CFD モデルでは構成粒子度変換された砕料粒子は常に ADEM 粒子として扱うこと

群を充填層とみなすことで，DEM-CFD カップリングモ

で，変形と破壊を考慮する必要のある砕料粒子について

デル[16]の枠組みを応用できるようにしているため，構のみ ADEM 粒子として扱うようにする。なお，DEM に成粒子よりも大きな流体セルサイズを用いても構成粒子

おいてめり込み量δの許容される最大値は粒子径の 0.1～に作用する流体力を計算することができ，大幅に計算時0.5%と経験的に知られていることから[19]，δmaxは砕料粒間を削減できる。そこで，構成粒子に作用する流体力Ff p

は DEM-CFD カップリングモデルに従い次式から推算

2.3 分散媒の運動
分散媒の運動は，Anderson と Jackson が提案した局所する。
Ff p = −εf βV pu − up − V p∇P
1 − εf f(15)

ここで，V pは構成粒子の体積，uf は分散媒の速度ベクト

ル，εf は分散媒の体積割合，Pは分散媒の圧力である。また ， β は 相 間 交 換 係 数 で あ り ， 川 口 ら [16] と 同 様 に ，Wen–Yu[17]と Ergun[18]の式をεf に応じて使いわけ，次のように計算する。
1 − εf2rpεf2150N
1 − εf μf+ 1.75ρf εf N1p ∑i =p 1 uf − up
2rpβ=when εf ≥ 0.8(16)

ここで，rpは構成粒子半径，μf は分散媒の粘度，ρf は分散媒の密度，N pは同一流体セル中に存在する構成粒子の数である。CDは流体抵抗係数であり次式より求める。
24 1 + 0.15Rep0.687
Rep
CD =
0.43

when Rep ≤ 1000

(17)

ここで，Repは粒子レイノルズ数であり，次式で計算する。
2rp ρf εf•μfNp

連続の式と運動方程式に対し近似解を求めることで表現する。以下の 2 式はそれぞれ，局所体積平均を施した連続の式と運動方程式である。
∂εf+ ∇ ⋅ εf uf = 0
∂t

(19)

μ
∂εf uf+ ∇ ⋅ εf uf ⋅ uf = εf − 1 ∇P + f ∇2uf + f f − p + f ball
ρf
ρf
∂t

ここで，tは時間，εf は分散媒の体積割合，uf は分散媒の速度ベクトル，Pは分散媒の圧力，μf は分散媒の粘度，ρfは分散媒の密度である。f f pは砕料粒子から分散媒への反作用力であり，同一の流体の計算格子（流体セル）中の構成粒子が分散媒から受ける作用力Ff pの合計値から計算される。 f ballは Eq. (9) で計算される流体力の反作用力であり，ここでは，分散媒の速度を媒体ボールの速度に修正する強制力と捉えられる。また，連続の式（Eq. (19)）

when Rep > 1000

N
∑i =p 1 uf − up

体積平均化処理された二流体モデル[20]の基礎式の内，(20)

when εf < 0.8
Npρf 1 − εf −2.7 1
3N p ∑i = 1 uf − up
4 CD 2rp εf

Rep =子径 0.44%（構成粒子径の 3%）に設定した。

と運動方程式（Eq. (20)）は Hirt らの Highly Simplified
Marker and Cell（HSMAC）法[21]により連成し，これら微分方程式の近似解は計算空間を有限の大きさを持つ直行格子（流体セル）で区切り差分近似することで求める。
分散媒の計算領域は Fig. 1 に示した領域と同じとして，(18)

前後，左右，上下面はすべてスリップ壁境界とした。流

砕料粒子の運動および破壊挙動は，Eqs. (12)～(18) から求められる構成粒子への作用力を Eq. (11) へ代入し解き，構成粒子の位置，速度，連結バネの結合状態を逐次更新することで表現する。
Filled up with dispersant (water)

一方で，砕料粒子の変形や破壊を考慮する必要があるの

Grinding ball
(Zirconia)
Φ 250 μm

は，砕料粒子がほかの物体と強く接触した場合だけで十分である。そこで今回は，強い接触が検出されるまで砕710 μm

ここで，ADEM は砕料粒子を複数の構成粒子の集合体（ADEM 粒子）として扱うため計算負荷が大きくなる。

料 粒 子 を DEM で 表 現 さ れ る 1 個 の 真 球 形 状 の 粒 子（DEM 粒子）として扱うモデルを導入する。初期状態にGrinding ball
(Zirconia)
Φ 250 μm

おいて，すべての砕料粒子は DEM 粒子として扱われ，このときの運動方程式や，その方程式を解くのに必要な接触力や分散媒から受ける流体力は 2.1 節で述べた媒体

Particles placeat random in 3-D space.

ボールと同様にして計算する。ここで，DEM 粒子では砕710 μm料粒子の変形や破壊を直接的に表現することは難しいものの，微小な変形であれば，DEM 粒子同士のめり込みにより表現される。そこで本モデルでは，DEM 粒子同士のVol. 61 No. 5（2024）

Particles(Alumina)
Φ 17.0 μm

Fig. 1

Calculation system for the analysis of particle grinding
behavior during wet ball milling261

体セルサイズについては，砕料粒子の粉砕挙動に影響をRelative ball velocity
vr,ball

およぼす可能性があるため，その粉砕挙動を解析できる
程度に十分小さくする必要がある。ここで，媒体ボール

Grinding ball

間の砕料粒子挙動を考えると，媒体ボール間に捕獲されLine throughtwo ball centers

Collision angle θc

るまでの砕料粒子は，媒体ボール間の流体流れに沿って移動し，その後，砕料粒子が媒体ボール間に捕獲され破壊されると考えられる。すなわち，捕獲後の砕料粒子はCollisionpoint

媒体ボール間によってその移動が拘束されるといえる。
したがって，接近する 2 つの媒体ボール間の砕料粒子挙Collision plane

動は，砕料粒子径（17.0 μm）よりも小さい流体セルサイズを用いていれば十分と考えられる。そこで本計算では，
流体セルサイズを砕料粒子径の半分程度の 7.5 μm とGrinding ballした。
3. 媒体ボール間の砕料粒子挙動解析におけるシミュ
レーション条件

Fig. 2

Definition of collision angle

突する 2 つの媒体ボールの中心間を結ぶ直線に対し垂直本シミュレーションでは，スラリー中で異なる衝突角な面を衝突面，2 つの媒体ボールの速度ベクトルの差を度で衝突する 2 つの媒体ボール間の砕料粒子挙動を解析

相対速度ベクトルとし，媒体ボールの衝突角度は，相対する。Fig. 1 に解析に使用する計算系を示す。計算系は速度ベクトルが衝突面に対し成す角度と定義する。すな

分散媒で満たされており，媒体ボールの初期表面間距離

わち，衝突角度が 90°のとき法線方向の衝突を，0°のとは 100 μm で一定なるように衝突角度に応じてずらしてき接線方向の衝突をそれぞれ意味している。

配置する。また，砕料粒子は計算空間中にランダムに配置する。

Table 1 にシミュレーションで用いた各種物性・特性値
および実験との比較から決定されるパラメータの内，主Fig. 2 に媒体ボールの衝突角度の定義を図示する。衝要なものをまとめる。媒体ボールはジルコニアボール，

Table 1 Simulation conditions

Dispersant (Water)
Density
Viscosity
Fluid cell size

kg/m3
Pa·s
μm

1.0×103
1.0×10−3
7.5

kg/m3
μm
%
-

3.95×103
17.0
25.0
0.5
1.0

N/m
μm
-

1.2×104
2.5
177
982
0.007

Particles (Alumina)
Density
Diameter Dp
Volume fraction
Friction coefficient
Rolling friction coefficient
ADEM parameters
Joint spring coefficient KS,nDiameter of constituent particle dp
The number of constituent particles in one particle
The number of joint springs in one particle
Maximum strain εmax-

Grinding balls (Zirconia)
Density
Diameter DB
Friction coefficient
Rolling friction coefficient
Young’s modulus
Poisson’s ratio
Normal spring coefficient between two balls
Initial velocity
Collision angle θc

262

kg/m3
μm
MPa
N/m
m/s
deg

6.0×103250
0.5
1.0
1.0×1020.30
5.8×104
0.50, 0.75, 1.00
0–90
in 15-degree increments

粉体工学会誌

Exp. (N=10) 24)

Sim.(εmax=0.007, KS,n=1.2×104 N/m)

The number of broken joint springs [-]90

Grinding ball

Grinding ball80

Grinding ballStress [MPa]

70
60
40

20 µm

(a) 0 deg

30

Grinding ball
20 µm(b) 45 deg(c) 90 deg

Tangential

20

NormalFig. 4 Broken behaviors of particles at different collision angles;

1000.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16
Strain [-]Fig. 3

Grinding ball

Grinding ball
20 µm

50

Comparison of the simulated stress-strain curve with theexperimental curves: the red curve is the stress-strain curvefrom the simulation using the best-fit condition, the green
and blue curves only change the joint spring constantand the maximum strain from the best-fit condition,
respectively(a), (b), and (c) are the collision angles 0 deg, 45 deg, and90 deg, respectively (Initial velocity of grinding balls is
100 cm/s.)

4. 媒体ボール衝突角度が砕料粒子挙動におよぼす影
響の解析
媒体ボールの衝突角度が砕料粒子の粉砕挙動におよぼす影響を観察するために，異なる媒体ボール衝突角度における砕料粒子の粉砕状態を可視化し観察する。Fig. 4

砕料粒子はアルミナ粒子，分散媒は純水をそれぞれ想定に初速 100 cm/s，衝突角度 0°，45°，90°で衝突した瞬間している。媒体ボールと砕料粒子に関する摩擦係数を実の媒体ボール間の砕料粒子の粉砕状態を衝突角度ごとに測することは容易ではないため，摩擦係数の中間的な値示す。図中の色は，構成粒子ごとに数えた破断した連結として 0.5 を仮定し与えている。粉砕機は湿式媒体撹拌

バネの本数を意味し，0 本のときは白，1 本以上のときは型ミルを想定し，媒体ボールの初速は湿式媒体撹拌型ミその本数が増えるごとに青から赤色に変化させている。

ル中の媒体ボール挙動を解析した報告[22,23]を参考に与衝突角度が法線方向に近づくと，より多くの連結バネがえており，2 つの媒体ボールそれぞれに同じ速さで逆向破断していることから，砕料粒子が大きく損傷している

きになるように与えた。ADEM のパラメータである連結

ことがわかる。Fig. 4 の砕料粒子の粉砕状態の観察から，

バネ定数および最大ひずみは，砕料粒子であるアルミナ

衝突角度と砕料粒子が受ける損傷の大きさには何らかの粒子の応力ひずみ線図が合うように決定される。今回は，関係があることが示唆される。そこで，この損傷の程度

微小圧縮試験機（島津製作所，MCT-510）によりアルミを次のように定量化し衝突角度との関係を調査する。まナ粒子（住友化学，A-26（50%粒子径 50 μm））一粒の加ず，実際に砕料粒子が破壊されるときを考えると，その重－変位線図を文献[24]より引用し，試験ごとの粒子径

破断面は結合により拘束されていた状態から解放されるから真球を仮定することで応力ひずみ線図を取得した。

と解釈することができ，ここで解放されたエネルギー分，一方，シミュレーションについては，平行平板間に砕料

砕料粒子は損傷を受けたと考えることができる。このこ

粒子を一つ静置し，片方の平板を 1 cm/s で移動させ圧縮とを ADEM の連結バネに置き換えると，連結バネの破断することで応力ひずみ線図を求めた。Fig. 3 は，実験かにより，連結バネのポテンシャルエネルギーが解放され，ら得られた応力ひずみ線図と，連結バネ定数と最大ひず

その分のエネルギーが砕料粒子の損傷の程度を表すと考みを調節したシミュレーションから得られた応力ひずみ

えられ，この解放されたエネルギーを解放エネルギーE R

線図を実験と比較したものである。図中の赤色の曲線がとすると，次のように定式化できる。

今回採用した条件（εmax = 0.007, KS, n = 1.2 × 10 N/m）を3用いたシミュレーションから得られた応力ひずみ線図であり，そのほかの黒色で示した曲線は砕料粒子 10 個をそ
れぞれ圧縮試験して得られた応力ひずみ線図を文献[24]より引用したものである。シミュレーションから得られ

N

b
2
2
E R = ∑ 1 KSn,i dp, iεmax, i = 1 N bKSn dpεmax
2
2
i=1

(21)

ここで，N bは破断した連結バネの本数を表す。Fig. 5 に，

この Eq. (21) から求めた解放エネルギーと衝突角度の関

た応力ひずみ線図は実験の応力ひずみ線図の群の中央付係を媒体ボールの初速ごとに示す。なお，図中のプロッ近に位置していることが見てとれる。ここで，応力ひずトは，砕料粒子の初期配置をランダムに変更し 6 回行っ

み線図の傾きが砕料粒子の剛性を，破断ひずみがその靭た平均値であり，エラーバーはその標準誤差を表す。媒性をそれぞれ表していることを加味すると，実際の砕料

体ボールの衝突角度が大きくなり法線方向に近づくほど粒子の平均的な剛性や靭性をシミュレーション上の砕料

解放エネルギーは増加しており，Fig. 4 の観察結果と一粒子に反映できていることがわかる。

致している。また，媒体ボールの初速の増加に対しても解放エネルギーは増加しているものの，衝突角度の増加に対し解放エネルギーも増加する傾向は一貫している。

Vol. 61 No. 5（2024）

263

Release energy ER [J]

4×10-93×10-92×10-9Initial velocity of grinding balls
1.00 m/s
0.75 m/s
0.50 m/s

Grinding ball

N=6

Shear

Particles

1×10-90Compression
030
60
Collision angle [deg]

90

Shear

Grinding ballCompression

Fig. 5 Relationship between collision angles and release energy:

Shearwhere the plots are average values, and the error bars are
standard errors

Fig. 7

Schematic diagram of the particle grinding behavior
between two grinding balls

6a に示した媒体ボールが法線方向に衝突する場合の砕料
粒子挙動を見てみると，媒体ボール間で複数の砕料粒子により形成される構造ごと圧縮され，砕料粒子が粉砕されていることがわかる。一方，Fig. 6b のように媒体ボールが接線方向に衝突する場合，砕料粒子が作る構造は破壊されているものの，砕料粒子自体の粉砕は法線方向のときほど起きていないことがわかる。Fig. 7 に媒体ボール間に捕獲される砕料粒子の構造と働く応力の概念図を示す。媒体ボールが砕料粒子群を捕獲するとき，砕料粒
子群が作る構造が，巨大な媒体ボールの表面から応力を受け変形させられている状態にあると考える。このように考えると，砕料粒子群が作る構造の変形は，媒体ボールが法線方向に衝突すると圧縮が支配的となり，接線方(a) 90 deg(b) 15 deg

The number of broken joint springs [-]Fig. 6

Snapshots of particle motion and breakage behavior when
the collision angles are (a) 90 and (b) 15 degrees,
respectively (Initial velocity of grinding balls is 100 cm/s.)

向に衝突するとせん断が支配的になる。ここで本計算では，スラリー中の砕料粒子の体積濃度は 25.0 v.%である
ため，砕料粒子が自由に移動できる空間が十分にあるといえる。こうした粒子濃度が疎な条件でせん断方向に応力が加えられたとき，砕料粒子群が作る構造は破壊できても，個々の砕料粒子は応力を逃れるように移動できるため，砕料粒子を破壊するだけの応力を砕料粒子自体にこれらのことから，媒体ボールの衝突角度が法線方向に加えることは難しい。一方で，圧縮方向への変形は，砕近づくほど，媒体ボールの衝突速度によらず砕料粒子が料粒子が存在する空間自体が押し縮められるため，媒体受ける損傷の程度は大きくなることがわかる。

ボールが圧縮方向に対して持つエネルギーが，捕獲され媒体ボールの衝突角度が砕料粒子の損傷の程度に影響る砕料粒子群の構造的な強度と砕料粒子自体の強度を超をおよぼした原因を調査するために法線，接線方向の衝えた分だけ破壊することができるものと考えられる。ま突それぞれの場合について，砕料粒子の粉砕過程を可視とめると，媒体ボールの衝突が法線方向に近くなるほど化する。Fig. 6 にシミュレーションから得られた媒体ボー

砕料粒子の損傷が大きくなった要因は，圧縮方向への変ル接近時の砕料粒子挙動のスナップショットを示す。Fig.

形の方が砕料粒子を破壊する観点では応力を砕料粒子自264

粉体工学会誌体に伝えやすいためといえる。なお，今回想定した湿式

度ベクトルであり，今回は初期の媒体ボールの速度ベク

媒体撹拌型ミルでは，質量濃度で 5 から 75 w.%程度で運トルから求める。Fig. 8 にEn，Et，En + Etと Fig. 6 で示し用されることが多く[25]，今回の砕料粒子の割合を質量た解放エネルギーERの関係をそれぞれ示す。法線方向成濃度に直すと 56.8 w.%であることから，実際の運用を考分の衝突エネルギーと解放エネルギーERの関係を表した

えても濃度は比較的高い部類に位置しており，今回のシ

Fig. 8a においてすべてのプロットが一本の曲線上に集中

ミュレーション条件における砕料粒子の空間的な自由度していることから，砕料粒子の解放エネルギーは，法線が特別高いわけではない。

方向の衝突エネルギーと強い相関があることがわかる。

以上の考察が正しいとすれば，砕料粒子の損傷の程度，

また，Fig. 8b から接線方向の衝突エネルギーが増加する

すなわち解放エネルギーは，砕料粒子群に強い法線方向

ほど解放エネルギーは減少していることから，接線方向の応力を作用させるほど大きくなると考えられる。したの衝突は砕料粒子に損傷を与えることは難しいことを表がって，解放エネルギーは，衝突する 2 つの媒体ボールしている。さらに Fig. 8c より，法線と接線を合計した衝の圧縮方向のエネルギーで整理できる可能性がある。こ突エネルギーは衝突角度が砕料粒子の損傷に与える影響

うした媒体ボールが衝突するときのエネルギーには，衝を考慮することは難しいこともわかる。したがって，媒突エネルギーが考案されており[26]，これを媒体ボール

体ボール衝突時に砕料粒子が受ける損傷は，衝突エネル

が一回衝突するときの法線方向成分 En，接線方向成分 Et

ギーの法線方向成分に支配される可能性が示唆される。

について衝突角度 θcを用いて記述すると，それぞれ次式で与えられる。

En = 1 mball vr,ball sin θc
2

2

(22)

2
E = 1 mball vr,ball cos θc
2
t

(23)

ここで，vr,ballは 2 つの衝突する媒体ボール表面の相対速

今回のシミュレーションで示されたような，法線方向の衝突エネルギーが砕料粒子の変形や状態の変化に影響をおよぼすことは実験的にも報告されている。たとえば，

4×10-91.00 m/s
0.75 m/s
0.50 m/s
N=61×10-90に損傷を与えることができる可能性があるといえる。

Initial velocity of grinding balls

3×10-92×10-9に衝突するように設計できれば，より効率的に砕料粒子

Release energy E R [J]Release energy E R [J]

4×10-9すなわち，同じ衝突エネルギーが生じる条件で湿式ボールミルを動かしたとしても，媒体ボール同士が法線方向

0.0

4.0×10-8
8.0×10-8
Normal impact energy E n [J]Initial velocity of grinding balls
1.00 m/s
0.75 m/s
0.50 m/s

3×10-9N=62×10-9
1×10-90

1.2×10-7

0.0(a) Normal component

Release energy E R [J]

4×10-93×10-92×10-94.0×10-8
8.0×10-8
Tangential impact energy E t [J]

1.2×10-7(b) Tangential component

Initial velocity of grinding balls
1.00 m/s
0.75 m/s
0.50 m/s
N=61×10-90

0.0

4.0×10-8
8.0×10-8
Impact energy E n+ E t [J]

1.2×10-7(c) Total
Fig. 8 Relationship between the impact energy and the release energy: where (a), (b), and (c) are normal component,

tangential component, and summation of normal and tangential components of the impact energy, respectively. The
plots are average values, and the error bars are standard errors

Vol. 61 No. 5（2024）

265

Fujita らは，乾式ボールミルにより鉄のアトマイズ粉をた。その結果，媒体ボールが高速かつ法線方向に近い衝処理し，鉄粉の塑性変形後の組織の特徴から，ボールミ

突角度で衝突すると砕料粒子がより大きな損傷を受ける

ルの変形は圧延に近いことを示した[27]。また，われわ

ことがわかった。また，砕料粒子の損傷の程度を ADEM

れの既往の報告では，湿式ボールミル中のメカノケミカ

で用いられる連結バネが破断するときに解放されるポテル反応で合成されるチタン酸リチウム水和物（LHTO）

ンシャルエネルギー（解放エネルギー）で定量化した。

の生成速度と，DEM シミュレーションを用いて計算され求めた砕料粒子の解放エネルギーは，媒体ボールの衝突

るミル中の媒体ボール挙動から求めた衝突エネルギーをエネルギーの内，法線方向成分に対し強い相関関係があ

比較し，衝突エネルギーの法線方向成分が LHTO の生成ることも明らかにした。この相関関係から，湿式ボール

速度に強く寄与していることを示した[28]。したがって，

ミル中の砕料粒子の破壊に対しては，衝突エネルギーの今回のシミュレーション結果は既往の実験事実とも符合

法線成分が強く寄与していることが示唆された。

する。
5. 結

［謝辞］本研究の一部は，ホソカワ粉体工学振興財団

言の HPTF21110 の 助 成 ， JSPS 科 研 費 JP20K22457 と湿式ボールミル中の砕料粒子の粉砕挙動におよぼす媒
体ボール衝突角度と衝突速度の影響について ADEM-

JP22K14525 の助成を受けたものです。ここに記して謝意を表明する。

CFD モデルを適用したシミュレーションにより解析しNomenclatureb: Radius of contact area

[m]

t

CD

: Drag coefficient[-]uball : Velocity of grinding balldp

: Particle diameter

[m]

uf

: Velocity of dispersant

[m/s]

E

: Impact energy

[J]

up

: Velocity of particle[m/s]: Release energy[J]vball : Grinding ball velocity at the center of fluid cell [m/s]ERf ball : Interaction between grinding ball and dispersant[m/s2]

f fp: Time[s][m/s]vr,ball : Relative velocity between two grinding balls[m/s]

V ball : Volume of grinding ball[m3]

: Interaction between constituent particle and dispersant

Vp: Volume of constituent particle[m3][m/s2]α: Volume ratio of grinding ball in the fluid cell

[-]
3

F

: Contact force[N]β: Interchange coefficient

Ff p: Fluid force acted on constituent particle

[N]

εf

: Volume ration of dispersant in the fluid cell[N]εmax : Maximum strain

FS: Joint spring force

I ball : Inertia moment of grinding ball[kg·m2]k: The number of contact elements

K: Elastic coefficient

[N/m]
[N/m]

[kg/( s·m ) ]
[-]
[-]

η

: Damping coefficient between two contact elements

ηS: Damping coefficient between two jointing particles[-][N·s/m]KS: Joint spring coefficientl0: Natural length of joint spring

[m]

θc

: Collision anglelS: Relative displacement vector of joint spring[m]μ: Coefficient of friction[-]

L: Relative displacement vector between contact objects

μf: Viscosity of dispersant[Pa·s][m]μr: Coefficient of rolling resistance

[kg]

ρf

: Density of dispersant: Rotational displacement vector of grinding ball [rad]mball : Mass of grinding ball[N·s/m]mp: Mass of constituent particle

[kg]

φ

Nb

: The number of broken joint springs[-]

ωball : Angular velocity of grinding ball

Np

: The number of constituent particles in fluid cell

[-]
[Pa]

P

: Pressure of dispersant

rball: Relative position vector from grinding ball center tocontact pointrp: Particle radiusrf: Relative position vector from grinding ball center to

Rep: Particle Reynolds number

Rr: Rolling resistance

fluid cell center.

266

[deg]

[-]
[kg/m3][rad/s]Subscript

i: Interest particle[m]j: Particle contacting the interest particle[m]n: Normal directiont: Tangential directionm: Particle connecting the interest particle with joint

[m]
[-]

sprint[N·m]

粉体工学会誌References[1] Dokuritsu gyoseihojin koonseramikkuzairyo dai124iinkai,Senshin seramikkusu no tsukurikata to tsukaikata, The
Nikkan Kogyo Shimbun, Ltd. (2010) p. 28.
[2] M. He, Y. Wang, E. Forssberg, Slurry rheology in wet ultra‐
fine grinding of industrial minerals: A review, Powder Tech‐nol. 147 (2004) 94–112.
[3] M. Malamatari, K. M. G. Taylor, S. Malamataris, D.
Douroumis, K. Kachrimanis, Pharmaceutical nanocrystals:
Production by wet ball milling and applications, Drug Dis‐
cov. Today 23 (2018) 534–547.
[4] T. Yokoyama, Y. Taniyama, G. Jimbo, Q. Zhao, The grinding
equilibrium size of in-water grinding of silica sand by a plan‐etary ball mill, J. Soc. Powder Technol., Japan 28 (1991)
751–758.
[5] K. Yokota, H. Watanabe, H. Tomino, Y. Kondo, Influence ofthe charged amount of ball media on grinding speed in wet
rotational ball milling, J. Jpn. Soc. Powder Powder Metal‐
lurgy 45 (1998) 680–684.
[6] S. Beinert, G. Fragniere, C. Schilde, A. Kwade, Analysis andmodelling of bead contacts in wet-operating stirred media
and planetary ball mills with CFD–DEM simulations, Chem.
Eng. Sci. 134 (2015) 648–662.
[7] R. M. de Carvalho, A. L. R. Oliveira, H. A. Petit, L. M.
Tavares, Comparing modeling approaches in simulating acontinuous pilot-scale wet vertical stirred mill using PBMDEM-CFD, Adv. Powder Technol. 34 (2023) 104135.
[8] K. Kushimoto, S. Ishihara, J. Kano, Development of
ADEM–CFD model for analyzing dynamic and breakage
behavior of aggregates in wet ball milling, Adv. Powder
Technol. 30 (2019) 1131–1140.
[9] P. A. Cundall, O. D. L. Strack, A discrete numerical method
for granular assemblies, Geotechnique 29 (1979) 47–65.
[10] S. P. Timoshenko, J. N. Goodier, Theory of Elasticity,
McGraw-Hill Book Company (1970).
[11] S. Natsui, S. Ueda, M. Oikawa, Z. Fan, J. Kano, R. Inoue, T.
Ariyama, Optimization of physical parameters of discreteelement method for blast furnace and its application to theanalysis on solid motion around raceway, ISIJ international
49 (2009) 1308–1315.
[12] A. Kondo, S. Ishihara, K. Kushimoto, T. Kozawa, J. Kano,
M. Naito, Correlation between grinding results in a tumblingball mill with liquid media and the analysis of ball motions
using DEM simulation, J. Soc. Powder Technol., Japan 56(2019) 148–155.
[13] S. Ishihara, Q. Zhang, J. Kano, ADEM simulation of particle
breakage behavior, J. Soc. Powder Technol., Japan 51 (2014)
407–414.

Vol. 61 No. 5（2024）[14] T. Kajishima, S. Takiguchi, H. Hamasaki, Y. Miyake, Turbu‐
lence structure of particle-laden flow in a vertical planechannel due to vortex shedding, JSME Int. J. Ser. B Fluids
Therm. Eng. 44 (2001) 526–535.
[15] Y. Yuki, S. Takeuchi, T. Kajishima, Efficient immersedboundary method for strong interaction problem of arbitrary
shape object with the self-induced flow, J. Fluid Sci. Tech‐nol. 2 (2007) 1–11.
[16] T. Kawaguchi, T. Tanaka, Y. Tsuji, Trans. JSME Ser. B 58
(1992) 2119–2125.
[17] C. Y. Wen, Y. H. Yu, Mechanics of fluidization, Chem. Eng.
Prog. Symposium Series 162 (1966) 100–111.
[18] S. Ergun, Fluid flow through packed columns, Ind. Eng.
Chem. 41 (1949) 1179–1184.
[19] P. W. Cleary, DEM prediction of industrial and geophysicalparticle flows, Particuology 8 (2010) 106–118.
[20] T. B. Anderson, R. Jackson, Fluid mechanical description of
fluidized beds. Equations of motion, Ind. Eng. Chem. Fun‐damen. 6 (1967) 527–539.
[21] C. W. Hirt, J. L. Cook, Calculating three-dimensional flows
around structure and over rough terrain, J. Comput. Phys. 10(1972) 324–340.
[22] Y. Yamamoto, R. Soda, J. Kano, F. Saito, DEM simulation ofbead motion during wet bead milling using an enlarged parti‐
cle model, Int. J. Miner. Process. 114–117 (2012) 93–99.
[23] R. Soda, J. Kano, F. Saito, Analysis of effect of pin configu‐
ration on beads motion in a stirred mill by DEM, J. Soc.
Powder Technol., Japan 46 (2009) 180–186.
[24] K. Kushimoto, Study for dynamic and breakage behavior ofparticles in wet ball milling., Ph. D. Thesis, Graduate Schoolof Environmental Studies, Tohoku University (2019) p. 162.
[25] T. Nakayamam, Cho biryushi·nanoryushi wo tsukuru bizu‐
miru, Kogyo Chosakai Publishing Co., Ltd. (2010) 106–112.
[26] H. Mori, H. Mio, J. Kano, F. Saito, Ball mill simulation in
wet grinding using a tumbling mill and its correlation togrinding rate, Powder Technol. 143–144 (2004) 230–239.
[27] S. Fujita, H. Sato, S. Motozuka, Deformation type during theball milling process: A comparative study of the microstruc‐tures formed by ball milling, uniaxial compression, androlling, Powder Technol. 426 (2023) 118598.
[28] T. Kozawa, K. Fukuyama, K. Kushimoto, S. Ishihara, J.
Kano, A. Kondo, M. Naito, Effect of ball collision direction
on a wet mechanochemical reaction, Scientific reports 11(2021) Article number 210.

267