離散要素法によるフェロコークス製造用竪型乾留炉の炉体大型化に伴う粒子の降下挙動変化の解析

Synopsis:

With the aim to decrease the reducing agent rate in blast furnace, the use of ferro coke containing metallic iron as a catalyst was proposed. This new burden material is expected to lower thermal reserve zone temperature. For the production of ferro coke, the use of a vertical coke oven is examined to realize the optimum carbonizing process. In this type of the oven, the burden materials have a distribution in residence time due to the friction effect on the wall and the interaction between the burdens. It has a great influence on the properties and quality of the ferro coke. Increasing the size of the vertical coke oven is essential for producing a large amount of ferro coke efficiency with stable quality.

In this study, the burden descending behavior in an enlarged vertical coke oven was calculated by Discrete Element Method (DEM). The effects of the shape of the discharging part and the discharging method on the residence time of the burden were analyzed based on the calculated results. In addition, considering the reinforcement of the enlarged coke oven, the burden descending behavior in the vertical coke oven with reinforcing horizontal beam inside was also examined.

The calculation results showed that distribution variation in the residence time could be decreased by increase in the size of the vertical coke oven. Furthermore, the results showed the possibility of controlling the residence time by controlling the discharge method.

1. 緒言

高炉で用いる炭材の反応性の向上は，熱保存帯の温度を低下させ高炉での炭材の利用効率を向上させると考えられている。このため高反応性コークスの製造および利用技術に関する研究開発が行われている^1,2,3,4)。この中で，微粉の石炭とペレットフィードから構成され，触媒作用のある鉄分を多量に含むフェロコークス³⁾は高反応性炭材として，基礎検討結果並びに工業プロセス化に向けた研究が報告されている⁴⁾。フェロコークスの製造工程では，従来の室炉式コークスと異なり，ブリケット等に塊成化した原料を乾留するため，竪型乾留炉の利用により炉頂より装入した塊を炉内で加熱し，炉下部で回収する連続式プロセスが従来プロセスに比べ効率的であると考えられる。

具体的な乾留方法を検討するために，前報⁵⁾ではフェロコークス製造用の竪型乾留炉における粒子降下特性を，冷間実験による観察，また個別粒子の運動解析が可能な離散要素法(以下，DEM:Discrete Element Method)⁶⁾を用いて定量的に検討した。前報では，フェロコークス製造に用いる竪型乾留炉内の固体流動，滞留時間分布などを解析し，適切な炉体設計，操業指針を得ることを目的とし，降下に伴う粒子拡散のモデル化を行った。さらに，そのモデルを用い，実機プロセスでの乾留時間の分布について検討した。この結果，粒子と壁の摩擦により壁に接した粒子の降下速度が遅く，粒子の滞留時間分布が拡大する原因になることが明確に示された。高炉内の装入物降下挙動については，壁からの一定距離の領域で壁との摩擦によって粒子混合が生じることが報告されている^7,8)。粒子の滞留時間を一定に制御するためには壁面の影響を小さくする必要がある。

ここで，移動する粒子量に対して壁面の効果を相対的に低下させるためには，炉体の大型化が有効な手段である。また大型化はフェロコークス実用化のための生産量の確保に必要な技術である。石炭の乾留には熱の供給が必要であり，壁面からの外熱や内部に高温ガスを流通することが考えられる。炉径の拡大は水平方向の温度分布を増大し，乾留の制御を困難にする可能性がある。そこで外熱による乾留を効果的に行うために矩形の炉体を用いることが提案できる。これを用いた場合，矩形の長辺より熱あるいは高温ガスを導入することにより温度分布を短辺の範囲に維持することが出来る。また，長辺を延長することにより乾留条件に影響を与えずにスケールアップを行うことが出来る点に長所がある。

本研究では，スケールアップを想定した矩形竪型乾留炉の形状での粒子の炉内降下挙動をDEMを用いて解析を行った。炉内の粒子は重力により自然に降下するものとして，炉底の形状および一度の排出量を変化させ，炉内粒子の縦方向の分散の変化を算出した。また，実炉体のスケールアップ時の強度保持を考慮し，炉内に円柱状の補強材を水平に設置した場合の，支柱の粒子運動への影響も考察した。

2. DEMによる炉内固体流れのシミュレーション

2・1 DEMの基礎方程式

DEMは対象とする全ての粒子の動きを接触力により決定する数値計算である。粒子に働く外力として，重力と，その粒子に直接触れる他の粒子，および壁の影響を考慮した。ある粒子に働く外力がわかれば，ニュートンの運動の第2法則に基づく運動方程式を解くことで，加速度が求められる。これを時間積分することにより粒子の変位が求まり，その運動軌跡を表現することが可能となる。

DEMについては既報⁵⁾で述べたためここでは概略の説明にとどめる。粒子間の応力にはVoigtモデルを用いた。このモデルでは接触点の法線，せん断方向にはそれぞれ2粒子の相対変位によって弾性変形し反力を発生するバネと相対速度に比例した抗力を発生する粘性ダッシュポットが並列に配置される。また，摩擦力による相互作用を考慮し，せん断方向には摩擦スライダーを導入した。このとき，2粒子間に働く力による接触粒子の並進の変位u_iは次式で表される。   

$$m_i\frac{d^2{u}_i}{d{t}^2}+\eta \frac{d{u}_i}{dt}+K{u}_i=0$$(1)

また，回転変位ω_iは次式で表される。   

$$I_i\frac{d^2{\omega }_i}{d{t}^2}+\eta \frac{d{\omega }_i}{dt}+K{R}_i{}^2{\omega }_i=0$$(2)

m，R，Iおよびωはそれぞれ粒子質量，半径，慣性モーメント，回転角速度を表す。また，K，ηはばね定数とダッシュポットの粘性係数である。一般の粉体では，1個の粒子の周りに多数の粒子が接触しているので，それら全ての接触点に対して(1)式，(2)式が成立する。

ある時刻Tにおける法線方向nの接触力F_nおよびせん断方向sの接触力F_sは，ばねによる弾性力eとダッシュポットによる粘性抵抗力dの和として，それぞれ(3)式，(4)式で表される。   

$$F_{n,T}=\sum {}_j\left(e_{n,T}+d_{n,T}\right)$$(3)

  

$$F_{s,T}=\sum {}_j\left(e_{s,T}+d_{s,T}\right)$$(4)

(3)式中のe_n,T，d_n,Tは(5)式で表される。   

$$e_{n,T}=e_{n,T+\Delta T}+K_n\Delta u_n , d_{n,T}={\eta }_n\frac{\Delta u_n}{\Delta t}$$(5)

ばね定数と粘性抵抗力は(6)式である。   

$$K_n=2D_b{\left(\frac{1-{\nu }_i^2}{E_i}+\frac{1-{\nu }_j^2}{E_j}\right)}^{-1},{\eta }_n=2\sqrt{m{K}_n}$$(6)

ここでEとνは粒子のヤング率とポアソン比，D_bは接触面半径を示す。また，(4)式中のe_s,T，d_s,Tは(7)(8)式で表される。   

$$e_{s,T}=e_{s,T+\Delta T}+K_s\Delta u_s,d_{s,T}={\eta }_s\frac{\Delta u_s}{\Delta t}$$(7)

  

$$K_s=8D_b{\left\{\frac{2\left(2-{\nu }_i\right)\left(1+{\nu }_i\right)}{E_i}+\frac{2\left(2-{\nu }_j\right)\left(1+{\nu }_j\right)}{E_j}\right\}}^{-1},{\eta }_s=2\sqrt{m{K}_s}$$(8)

ΔTはタイムステップである。粒子間せん断方向について，クーロンの摩擦法則から得られるせん断力より大きな力は作用しない。せん断応力が粒子間摩擦力より大きい場合はすべりが生じるため，摩擦係数をμとして(9)式が成り立つ。   

$$F_s=\mu \left|F_{n,ij}\right|$$(9)

DEMにおいて個々の粒子は完全な真球として取り扱われる。しかし，実際の粉体は不定形状であり，真球に比べて回転しにくい。このような粒子形状を表現する為に，回転方向のトルクTと逆向きの回転抵抗Mを導入した。   

$$I_i\frac{d{\omega }_i}{dt}=\sum {}_j\left(T_{ij}+M_r\right)$$(10)

  

$$M_r=\frac{3}{8}\alpha \cdot D_b\left|F_n\right|$$(11)

T_i,jは粒子iが受けるトルク，M_rは接触面でのせん断応力を示す。αは回転抵抗係数であり，粒子表面物性に基づいて決定した⁹⁾。

2・2 DEMの計算条件

パイロットプラントスケールの竪型乾留炉をスケールアップすることを想定してDEMにより竪型炉の解析モデルを作成した。Fig.1に示す(a)のように上部の幅1 m奥行1 mの正四角柱と下部四面が斜面して狭くなるテーパー部により構成される高さ16 mを一つのユニットとし，これを並列に並べ大型の竪型乾留炉を構成した。各ユニット間には壁面は存在しない。(b)は本計算の対象とした五連ユニットの形状を示す。下部は個々のユニットが独立したテーパー部を有する形状(c)と，五連のユニット全体が一つのテーパー部を有する形状(d)を想定し，炉下部形状が与える影響について解析を行った。また，炉体補強用の支柱を用いることを想定し，二番目と四番目のユニットの中央，下から8.0 mの位置にユニットの列と直交するように支柱を設置した計算も行った。支柱は円柱形で任意に径を設定した。

Fig. 1.

 Schematic diagram of furnace calculated in present study.

粒子の移動による降下の表現のために，炉底部に粒子が消失する領域を設定し，その領域に存在する粒子を一定時間間隔で消去することで粒子の排出を表現した。炉底から上方に15 mの位置にストックラインを設け，充填層の高さがこのストックラインを超えるように逐次粒子を上部より落下し，供給を行った。移動層の降下は，炉底部での消失領域の体積と粒子消去間隔を変化させることで調整した。粒子径は，フェロコークスの実粒子径を想定して6.0 cmとした。粒子の密度は，フェロコークスの値⁵⁾を用いた。フェロコークスは乾留前後によってその物性が異なるが，乾留前後の測定値の平均値を用いた。粒子間および粒子−壁間の接触摩擦係数はフェロコークス粒子，また耐火物壁を用いたせん断試験により実測して求めた値⁵⁾を用いた。計算に用いた計算条件を，Table 1に示す。

Table 1.  Calculation conditions.

Particle radius [cm]	3.00
Particle Young modulus [GPa]	0.54
Particle Poisson's ratio [-]	0.22
Particle density [g/cm3]	1.05
Friction coefficient (particle-particle) [-]	0.25
Particle rolling friction coefficient [-]	0.30
Friction coefficient (particle-wall) [-]	0.45
Time step [s]	1.00 × 10–4
Particle number [-]	80000 (1unit) 400000 (5units)

フェロコークス製造プラントでの原料の滞留時間は数時間に及ぶため，実際の時間スケールで計算を行うことは極めて困難である。本研究では，シミュレーションでの排出時間間隔の短縮により計算時間の短縮を行った。計算結果を用いて，滞留時間の換算を行い実炉での滞留時間分布を求めた。シミュレーションでの降下挙動と実炉での降下挙動との関係をFig.2に示す。一定の時間間隔で粒子の排出を行うため，降下挙動は図中step1~3の様に排出，降下，停止に分けられる。シミュレーション(a)では計算時間短縮のため，排出間隔を小さくしている。実際の降下挙動(b)ではstep1,2に対して，step3にかかる時間が極めて長いため，step3にかかる時間と排出間隔はほぼ等しくなる。このとき，一回の排出で原料が降下する距離は排出量によって決まる。このため，同一の排出量の場合，原料がある距離を降下するのに必要な排出回数は，シミュレーションと実炉で同じ回数となる。シミュレーションで求めた排出回数と任意に設定した排出時間間隔から，実炉での滞留時間分布を求めることができる。

Fig. 2.

 Comparison between residence times of particle discharging in simulation and actual behavior.

3. 炉内充填層降下挙動の解析結果

3・1 炉下部形状の影響

粒子の降下挙動が把握しやすいように，炉底から3 m，6 m，9 m，12 mの位置に他の粒子と物性が同一であるトレーサー粒子を設定した(Fig.3)。排出部の下部0.5 mに存在する粒子を1.0秒ごとに消去し，充填層の降下を制御すると17秒後にFig.3(b)のように壁面近傍の粒子の降下が遅れ，中心部の粒子が先行する。Fig.4に20秒後と27秒後のトレーサーの分布を示す。壁面近くの粒子の降下が遅いことが判る。図の左右の壁付近の粒子は図の手前および奥の壁の影響を受けるため，特に降下が遅れることが判る。この二面の壁付近の粒子の遅れは前報の結果とほぼ一致する。左右の壁は図の左右方向に1ユニット程度の範囲で粒子の降下に影響を与えるため，中央の3ユニット内の粒子は手前と奥の壁の影響は明確に受けているが，左右の壁の影響をほとんど受けない。この壁面影響の減衰はFig.5に示すように実験でも確認できる。この実験は実炉の1/20スケールの容器にφ2 mm×3 mmの円筒形プラスチックペレットを充填し，炉下部より往復プラウを用い排出した際の，降下挙動の観察を行ったものである。初期条件(a)では炉頂部に位置する白色のペレットが，乾留が終了する位置まで到達した際(b)には，左右の炉壁面でのみペレットの遅れが確認でき，炉中心部ではほぼ均一に降下している。このことから炉のサイズの拡大は降下挙動の均一化に有効であることが確認できる。

Fig. 3.

 Positions of tracer at 0 s and 17 s.

Fig. 4.

 Tracer particle distribution at 20 s and 27 s.

Fig. 5.

 Descent behavior observation by the cold model.

次に27秒後の下部のテーパー部を各ユニット毎に設置した場合(a)と，全体が一体化したテーパー部を設置した場合(b)での比較を行う。テーパー部壁からの粒子に加わる応力の影響で，(a)では中央のユニットの粒子の降下が遅く，(b)では中央が速い降下を示した。これは，下部のテーパーによる縮流と摩擦が(a)ではより大きく影響し，粒子の降下が妨げられ，粒子の降下速度に差が生じたためであると考えられる。

応力分布をFig.6に示す。(a1, b1)と(a2, b2)はそれぞれ炉内上部と下部(位置)において，粒子が降下し応力付加が起きた時の応力分布を示している。粒子には排出時の降下に合わせて上方に応力が伝搬するようにかかることが分かった。カラースケールはある1つの粒子に作用する粒子間応力σを以下の式から導出した値を示す。   

$${\sigma }_{n,T}=\frac{{\sum }_j\left(e_{n,T}+d_{n,T}\right)}{4\pi R^2}$$(12)

Fig. 6.

 Influence of discharge area shapes on stress distribution.

ここでRは対象粒子の半径を示す。下部に個別のテーパー部を設置した(a)では局所的にテーパー部から粒子に応力がかかるため，炉下部の粒子のほとんどは応力の低い値を示す。一方，(b)では手前と奥の傾斜した壁から粒子全体に応力がかかるが，全体の粒子に均一な力がかかり，高い応力値を示す粒子が存在しないことが判る。

炉内での粒子の圧潰を防止するためには，最大応力の小さい(b)の方が有利である。また，フェロコークスは乾留により強度が上がるため，炉下部の応力がかかる部位では，強度が十分に上昇した製品が通過するように調整を行うことが必要となる。

3・2 粒子の排出方法の滞留時間分布への影響

乾留工程において高温場での乾留物保持時間は製品の性状に大きく影響するため，竪型乾留炉の縦方向中央部分を乾留温度域として，DEMにより粒子の滞留時間を算出した。Fig.7に濃灰色で示すように，炉下部より9 mの位置に0.2 m厚のトレーサー層を配置し，炉下部より4 mの位置をトレーサー粒子が通過する時間あるいは割合を算出した。炉下部の排出は掻き出し装置による間欠排出を再現するために，炉下部の0.6 m(a)，0.3 m(b)厚の粒子層を同時に消去し，粒子層を降下した。このような間欠排出ではFig.2の模式図に示すように，掻き出しにより粒子が排出され，その直後に粒子が重力により自由落下し，全体の粒子が落下後に動きが止まり，充填層が安定する。このサイクルを繰り返すことで充填層全体が降下する。シミュレーション上では1.0秒ごとに排出を行った。実炉での装入物全体の平均降下速度が2.5[m/h]になる様に(a)，(b)の排出間隔1.0 sを，それぞれ実炉での6.2 min，2.85 minに換算して滞留時間分布を求めた。この時，最初に排出されるトレーサーの滞留時間は2時間となる。

Fig. 7.

 Initial position of tracer and region for residence time calculation.

上記の条件で炉中間部を縦方向に5 mトレーサーが移動した割合を，Fig.8に示す。縦軸は各排出の後，充填層が安定した時に高さ4 mの位置を通過したトレーサーの割合を示す。横軸は最初に通過したトレーサーを基準(0 min)とした経過時間を示す。プロットは各排出時での通過割合を示している。通過割合が増加するに従い，一回当たりの排出での通過割合の増加は低減する。トレーサーの99%が通過するには最初のトレーサーの通過から，40 minほどの時間を要する。一回の排出当たりの粒子の排出量が多く，排出回数の少ない(a)の条件では最初のトレーサー検出位置に達してから2回の排出で8割程度，3回の排出で9割以上のトレーサーが高温部を通過することが判る。(b)では通過割合が8割，9割に達するまでの排出回数はそれぞれ4回，7回である。一回当たりの排出間隔を考慮すると，トレーサーが8割通過するまでにかかる時間の差は約10分，9割通過するまでにかかる時間の差は約15分である。これらの時間差は基準となる2時間に対してそれぞれ約8.3，12.5%の差となる。

Fig. 8.

 Effect of discharging region and frequency on relationship between residence time and passed tracer ratio.

トレーサーの上下方向の分布は壁の摩擦によるため，Fig.4のように壁面近傍の粒子が遅れ上側に，中央付近の粒子は下側に分布する。特に両端の粒子の遅れが大きいため，両端のユニットからの排出速度を高め粒子の降下を調整することを試みた。

Fig.9に示すように，(a)均一に炉下部30 cm，(b)両端の炉下部40 cm，中央の炉下部30 cm，(c)両端の炉下部50 cm，中央の炉下部30 cmを1回あたりの粒子排出の領域とした。Fig.10は各排出条件での計算開始から実炉換算で100 minが経過した時点の粒子の分布を示す。両端部の排出量を増やした(b)，(c)では，テーパー部でのトレーサー(黄色)の分布に大きな差がある。これに対して炉体中部のトレーサー(赤色)の分布はテーパー部に比べ，変動が小さいことが確認できる。これらより，狙い通り排出方法により壁面部での遅れを緩和することが可能であることが分かった。本研究では，炉底より4 mから9 mまでの区間を乾留温度域として設定したため，炉下方テーパー部でのトレーサー分布の変動は乾留時間に影響を与えない。

Fig. 9.

 Variation of discharging region for controlling residence time of particle.

Fig. 10.

 Tracer particle distributions for discharging conditions (a), (b), and (c).

Fig.11に炉中間部を縦方向に5 mトレーサーが移動した割合と時間の関係を示す。Fig.9(a)の条件では，初期の粒子排出が均一であるが，通過割合が0.96以上で両端のユニットの角の粒子の停滞が顕著に表れる。(b)では0.96以上での粒子の通過割合の上昇が(a)より高く壁の分布拡大への影響を緩和していることが判る。(c)では中心部ユニットから端のユニット領域への粒子の流れ込みが顕著になり全体にトレーサーの分布が大きくなることが判る。しかし，最初のトレーサーが5 m降下してから9割9分が通過するまでの時間は条件(a)よりも短く，極端に長時間滞留する粒子割合が低下することが判る。

Fig. 11.

 Effect of discharging condition on relationship between residence time and passed tracer ratio.

粒子の滞留時間を均一に制御するためには(b)に示すように並列したユニットの両端の降下速度を中心部より若干高め，壁の摩擦による降下分布を減少することが有効であると考えられる。

3・3 炉内支柱の粒子降下挙動への影響

3・3・1 支柱径の粒子への影響

矩形に炉体を拡張する際に，短辺方向への炉体の変形を防ぐために，炉内に支柱を設置した際の粒子降下挙動への影響を計算した。支柱の太さは直径20 cmと40 cmの円柱とし，2・1項で述べたように，炉底から8 mの位置に2本設置した。ここではトレーサー粒子を炉下部から3 m，6 m，9 m，12 mの位置に設置し，粒子の降下による分布の変化を計算した。シミュレーション上では炉底から40 cmの位置を粒子の排出範囲として1.0 s毎に排出を行った。Fig.12に降下開始直後，4.0 s経過後，8.0 s経過後を示す。トレーサー粒子が支柱の周辺で動きが変化し，分布に乱れが生じることが判る。8.0 s経過後の分布から，20 cmの支柱では支柱の位置のトレーサー粒子が遅れ，40 cmの支柱では支柱位置のトレーサー粒子が先行することが確認できる。

Fig. 12.

 Tracer particle distribution at 0 s, 4 s and 8 s with beam.

Fig.13にトレーサー粒子が支柱の上部に達した時と，下部に達した時の速度ベクトルを示す。支柱上部において20 cmでは支柱直上の粒子速度が低下し，支柱の間の粒子の速度が上昇しているが，40 cmでは支柱上部の粒子が相対的に速い速度で支柱の間に流れ込むことが判る。また支柱下部においても，速い速度で支柱の下に粒子が流れ込むことが判る。この挙動から，40 cmではもともと支柱の間に位置したトレーサー粒子が，支柱の下側に回り込み降下が先行する。また，支柱の上部に位置した粒子は支柱の間を遅れて降下したことが判る。

Fig. 13.

 Influence of diameter of beam on particle velocity at vicinity of beam.

Fig.14に支柱周辺における粒子降下時の応力分布を示す。色は応力の強度を示す。両条件も支柱の直上の粒子に大きな応力がかかっていることが判る。40 cmでは支柱間への流れ込の影響で，支柱の間での応力が大きくなっていることが判る。また，支柱直下に応力が低下する領域が確認できる。粒子にかかる応力の観点から，径の小さい支柱が有利であることが明確に示された。

Fig. 14.

 Influence of beam diameter on stress distribution.

3・3・2 支柱径の調整による粒子降下の均一化

支柱の径により支柱の位置でのトレーサーの遅れや先行が確認された。ある径の支柱を用いれば，支柱により粒子の遅れと先行が相殺される可能性がある。そこで，支柱の直径を1 cm単位で変化させて粒子の降下挙動を算出した。Fig.15に支柱直径35 cm，36 cm，37 cmにおける，支柱位置を通過後約4 m降下したトレーサー粒子の分布変化を示す。通過後のトレーサー粒子の分布を35 cmでは支柱の影響で支柱位置の粒子の遅れが確認できる。また，37 cmでは40 cmでの条件に比べると分布が小さいが，支柱位置の粒子の先行が確認できる。36 cmでは粒子の遅れ，先行の両方が見られるが，他の条件に比べ上下方向の分布が小さくなることが判る。支柱周り粒子の遅れと先行の境界は36 cm付近にあると考えられる。また，このような支柱を選択することで，粒子の滞在時間の分布を縮小できる可能性がある。しかし，支柱の影響による粒子流れの乱れは上下方向の分布を拡大する。この乱れを無くすことは困難であるため，粒子の乾留時間の均一化には支柱を用いないことが望ましい。

Fig. 15.

 Effect of beam diameter on tracer particle distribution.

4. 結言

パイロットプラントスケールのフェロコークス竪型乾留炉を並列に配置し，矩形竪型乾留炉にスケールアップすることを想定し，炉内での粒子の挙動をDEMにより解析を行った。炉下部の形状および排出方法と粒子の滞留時間および粒子にかかる応力の関係を導出した。また，実炉体のスケールアップ時の強度保持を考慮し，炉内の補強材の粒子運動への影響も考察した。以上より以下の知見を得た。

1)炉体を並列化することで，2面の炉壁により降下が遅れる角の影響が小さくなり，全体の降下速度が均一化する。

2)ユニット毎に排出口を設置する場合，粒子の滞留時間への影響は小さいが，炉下部の縮流により粒子へ応力が集中する。

3)並列化した端のユニットの粒子降下が遅くなるため，端の排出速度を中央部に比べ上昇することで，粒子の滞留時間の均一化が図れる。

4)炉内に支柱を設置した場合，支柱の直径により支柱位置の粒子の降下が促進あるいは遅滞する。その境界の支柱径を用いることで上下方向の粒子の分布拡大を緩和できる可能性がある。

謝辞

本研究はNEDO「資源対応力強化のための革新的製銑プロセス要素技術開発事業」の一環として行われたものであることを記すとともに，実験結果の提供を頂きましたJFEスチール研究所の関係者の方々に謝意を表します。

文献

• 1)   S.  Nomura,  H.  Ayukawa,  H.  Kitaguchi,  T.  Tahara,  S.  Matsuzaki,  M.  Naito,  S.  Koizumi,  Y.  Ogata,  T.  Nakamura and  T.  Abe: ISIJ Int., 45(2005), 316.
• 2)   S.  Nomura,  H.  Kitaguchi,  K.  Yamaguchi and  M.  Naito: Tetsu-to-Hagané, 92(2006), 481.
• 3)   T.  Yamamoto,  T.  Sato,  H.  Fujimoto,  T.  Anyashiki,  M.  Sato and  K.  Takeda: Tetsu-to-Hagané, 96(2010), 683.
• 4)   T.  Yamamoto,  H.  Fujimoto,  T.  Anyashiki,  K.  Fukada,  M.  Sato,  K.  Takeda and  T.  Ariyama: Tetsu-to-Hagané, 97(2011), 501.
• 5)   T.  Kon,  K.  Kojima,  S.  Natsui,  S.  Ueda,  R.  Inoue and  T.  Ariyama: Tetsu-to-Hagané, 98(2012), 459.
• 6)   P.A.  Cundall and  O.D.L.  Strack: Geotechnique, 29(1979), 47.
• 7)   M.  Ichida and  H.  Biausser: Tetsu-to-Hagané, 69(1983), S61.
• 8)   M.  Ichida,  K.  Tamura,  Y.  Okuno,  K.  Yamaguchi,  M.  Nakamura and  M.  Nakayama: Tetsu-to-Hagané, 78(1992), 58.
• 9)   S.  Natsui,  S.  Ueda,  M.  Oikawa,  Z.  Fan,  J.  Kano,  R.  Inoue and  T.  Ariyama: ISIJ Int., 49(2009), 1308.