技術論文＞セルラーオートマトン法を用いたオーステナイト粒の異常粒成長シミュレーション

3

技術論文
Technical Paper

セルラーオートマトン法を用いたオーステナイト粒の異常粒成長シミュレーション
村田憲治＊1，2，福井ちひろ＊2，孫飛＊3，陳達徳＊3，足立吉隆＊3

Simulation of Abnormal Grain Growth Using the Cellular Automaton
MethodKenji MURATA, Chihiro FUKUI, Fei SUN, Ta-Te CHEN and Yoshitaka ADACHISynopsis
The abnormal grain growth of steel, which occurs during carburization, adversely affects properties such as heat treatment
deformation and fatigue strength. This study aimed to control abnormal grain growth by controlling the materials and processes.
Thus, it was necessary to investigate the effects of microstructure, precipitation, and heat treatment conditions on abnormal grain
growth. We simulated abnormal grain growth using the cellular automaton (CA) method. The simulations focused on the grainboundary anisotropy and dispersion of precipitates. We considered the effect of grain boundary misorientation on boundary
energy and mobility. The dispersion state of the precipitates and its pinning effect were considered, and grain growth simulations
were performed. The results showed that the CA simulation reproduced abnormal grain growth by emphasizing grain boundary
mobility and the influence of the dispersion state of the precipitate on the occurrence of abnormal grain growth. The study findings
show that the CA method is a potential technique for the prediction of abnormal grain growth.

1．緒言

大きなひずみの影響により，浸炭工程にてオーステナイト粒が異常粒成長を起こし，熱処理ひずみや，疲労強度自動車向けの動力伝達部品には歯車機構が用いられといった歯車に必要とされる特性に影響を与えてしまう

る．こうした歯車には高い表面硬さと高い靭性を両立さことが知られている 1）．従って，異常粒成長を予測するせることが求められるため，JIS-SCr420 などの表面処理ことは工業的に重要であるが，析出物の量，分散状態，用の機械構造用鋼に対して浸炭処理を行われることが一初期オーステナイト粒度といったいくつかの材料学的な般的である．近年，カーボンニュートラルの観点から，要因を実験的に定量化し，異常粒成長の過程を熱処理中自動車の製造時の CO2 排出量削減が求められている．に観察することは困難である．よって，異常粒成長を再CO2 を排出する自動車用歯車の製造工程の一つに，熱間現できるセルラーオートマトン（以下，CA という）法鍛造にて素形材を成型する工程がある．この工程を冷間を研究することとした．

鍛造とすることで，製造時の CO2 排出量の削減が可能となる．しかし，冷間鍛造とすることで，冷間鍛造時の粒成長には，例えば，モンテカルロ法 2）～ 5） や，CA法 6）～ 9），フェーズフィールド法 10），11） といったいくつ

2024年 3月 28日 受付＊ 1 大同特殊鋼㈱技術開発研究所（Corporate Research & Development Center, Daido Steel Co., Ltd.）＊ 2 名古屋大学 大学院工学研究科（Graduate School of Engineering, Nagoya University）＊ 3 名古屋大学 大学院工学研究科，工博（Dr. Eng., Graduate School of Engineering, Nagoya University）

4

電気製鋼 第 95 巻第 1 号 2024 年

かの方法が知られている．こうした方法は定常粒や異常の粒成長の駆動力は粒界の曲率半径を用いて式（2）のよ粒成長をシミュレーションするために用いられてきた．

うに示される．

例えば，早川らは，モンテカルロ法を用いて粒界方位差を導入した異常粒成長シミュレーションを行った 12），13）．

（2）

この中で異常粒成長は粒界の方位差が粒界エネルギーと易動度に影響を与えるときに起こると述べている．まここで，Rgbは曲率半径，σは粒界エネルギーである．また，Ye らは CA 法に粒界の異方性を導入することで，た，粒界の易動度は温度依存性があるとされ，式（3）の異常粒成長をシミュレーションした

ように示される．

14）．

一方で，析出物のピンニング効果は異常粒成長に重要な役割を果すことが知られている 15），16）．実用材におい（3）ては，こうした析出物は例えば，ミクロ偏析や，浸炭のここで，M0は易動度の係数，Q は活性化エネルギー，R影響でマトリクスの中に偏在する場合がある．浸炭によは気体定数，Tは温度である．

る析出物への影響を導入して粒成長のシミュレーション式（2）および式（3）に必要な定数についてはすで．彼らは NbCに多くの知見があるが，曲率を CA 法で扱う方法は一般フェーズフィールド法により異常粒成長をシミュレー扱うためである．本研究ではこの問題を解決するため，ションした．しかし，粒界の異方性と析出物の分散状態

Pierre らによる Box Count 法 19）を用いることで，曲率にを行った例として木下らの報告がある

17）のピンニング効果と浸炭における C 量の影響を導入し，をともに導入したシミュレーションの報告はない．
これまでの CA 法では，結晶粒の方位として，粒番号を割り当てる方法が用いられることが多かった 18）．こ的ではない．これは CA 法では空間を離散的に分割してよる粒成長の駆動力への影響を導入した．Box Count 法
では式（4）のように曲率 RL を表す．の方法はシンプルであるが，後述するように，実際の方（4a）位差分布と異なる方位差の分布となってしまう．こういった不確実性を排除すべきである．
本研究では，CA 法を用いて異常粒成長に対する粒界（4b）の易動度と析出物の分散状態の両者に着目した．その中で，粒番号よりも，結晶粒方位を割り当てる方法によりシミュレーションを行った．
フェーズフィールド法とは異なり，CA 法では離散的（4c）に界面を扱うため，粒界エネルギーの物理的な扱いが可能で，粒成長挙動を直感的に理解することが可能である．また，フェーズフィールド法よりも小さな計算コストで計算が可能である．こうした理由から本研究では

ここで，bは Fig. 1に示すような，セル x, yの周囲で考えるべきセル範囲を表す．例えば，b=5であれば第二近接のセルまでを対象として曲率を計算する．n othersは周囲[x, y][x+i, x+j]の異なる粒のセル数である．δ
はセル x, yとセルCA 法を用いた．

x+i, y+jとが異なる結晶粒であるとき 1，同一の結晶粒で2．計算方法

あれば 0である．grain[x, y]はセル x, yにおける粒の識別番号である．

オーステナイトの粒成長挙動を調査するため，以下に示すように，二次元の CA 法を用いた粒成長モデルを構

Range of cells to be considered

築した．
初めに粒成長の一般論を CA 法に落とし込むこととする．粒成長速度 v は単位体積当たりの粒成長の駆動力∆G と粒界の易動度 M の積として式（1）のように示される．

[ , ]

（1）
ここで x, yはセルの二次元座標である．単位体積当たりFig. 1. Schematic of the range of cells considered in thebox-counting method.

技術論文＞セルラーオートマトン法を用いたオーステナイト粒の異常粒成長シミュレーション

5

次に結晶粒の扱いについて述べる．CA 法やフェーズ
フィールド法ではそれぞれのセルには特定の粒番号を入間での粒界の方位差 θ x+i, y+jとして粒番号の差の絶対値をx, y最大の粒番号で割る式（5）に示す方法がある．

（5）Standardized frequency

れることが多い．この手法を用いた場合に二つのセルのここで，G maxは粒番号の最大値である．しかし，次の2点で実組織と異なってしまう．

Misorientation

1．最大の方位差が 90 °となってしまうこと．

(a)2．粒同士の方位差の分布をとると，高角側に向かって単調減少していく方位差分布となってしまうこと．
織と同様に，結晶方位を与えることとした．その際には，計算の都合上，方位を表す四元数を用いることとした．
四元数を用いて粒方位を表すことは，例えば高橋ら 20）によって提案されている．
立方晶の結晶 a,b の方位を表す四元数をそれぞれとq a, q b すると，この方位差 θa,b は式（6）のように，平易Standardized frequency

そこで，本研究では，各結晶粒に粒番号ではなく，実組に表される．
（6）
ここで，eiは立方晶の対称操作を表す 24種類の四元数である．また，（p, q）は四元数 p, qの内積を表す．オイ
ラー角や格子ベクトルを用いる場合には，四元数を用いる場合と比較してより複雑な操作が必要であり，四元数を用いることでより低い計算コストで方位差が求められる．また，四元数の 4つの実部は等価な扱いが可能なため，ランダム方位を与える際に，4つの実部に一様乱数を与え，規格化するのみで容易にランダム方位を与えらMisorientation(b)
Fig. 2. Distribution of grain boundary misorientation using(a) grain numbers and (b) random crystallographicorientations. Dashed line represents ideal random
misorientation distribution.21)
ギー，E0は高角粒界のエネルギー，θmはエネルギーが高角粒界と同じとみなせる閾値である．θmは 15 °を用い18）た．各セルの粒界エネルギー σ[x, y]は周囲のセルとの粒による方法で粒界方位差を定義した場合と，直接粒方位界エネルギーの和として，式（8）のように表し，周囲のをランダムに与える本手法で粒界方位差を定義した場合セルはムーア近傍のセルとした．

れる利点がある．結晶粒に粒番号を与え，Contieriらとで，方位差の分布をそれぞれ Fig. 2(a) および Fig. 2(b)に示す．直接粒方位を与える手法では，理論的な分布 21）に近い方位差の分布が得られることが分かる．

（8）

次に，粒界エネルギーの扱いについて述べる．粒界エ

粒界の易動度については，Humphry ら 24）による式（9）ネルギーには Read らによる式（7）を適用し，粒界エを適用することで，易動度の方位依存性を導入することネルギーの方位依存性を取り入れることとしたができる．さらに，易動度の不均一性を考え，特定の粒22），23）．

のみに高い易動度を割り当て，易動度を強調したモデル（7）
ここで，Gbx+i, y+jはセル x, yとセル x+i, y+jとの粒界エネル
x, yを考案しシミュレーションを行った．この易動度の強調は例えば，特定の粒のみが周囲の粒と異なる方位や対応粒界を持ち，易動度が高くなる場合などを想定している．

6

電気製鋼 第 95 巻第 1 号 2024 年

Start

（9）

= 0, cell[ , ; 0]

ここで，θmは易動度が大角粒界とほぼ同等になる角度を

LoopTime step< 1000示し，本研究では 15 °を用いた．

Loopcell [ , ]0 ≦ x ＜Xmax,0 ≦ y ＜Ymax（10）

Calculate[ , ]ここで，pは強調する場合の粒界易動度の係数である．
最後に析出物の影響について述べる．熱処理中に析出物は結晶粒界に対してピンニング効果を持ち，結晶粒粗大化を抑制するとされる．Zener によると 25），ピンニング効果を式（2）に対して導入すると，単位体積当たりの粒成長の駆動力は二次元の Zener の式として，式（11）のように示される．

min Δ,,,< ΔG x, y,>YesNo[ , ; + 1] =[ , ; ]

[ , ; + 1] =[ + , + ; ]Loop [ , ]
Loop

（11）
ここで，fとrはそれぞれ，析出物の体積分率，析出物の半径である．本研究ではこの Zenerの式を用いた．
さらに，組織の中で偏析などの影響により，ピンニング粒子に偏りがある場合を想定して，体積分率 f が勾配を持つよう，式（12）に示すように体積分率に分布を与End

Fig. 3. Overall flow of the CA simulation.
粒成長の CA 法による計算に与える初期組織は相変態による CA 法を基にしたもの 26） とし，平均粒径が異なるいくつかの初期状態からスタートして CA 法による計算を行った．今回は 1 ピクセルを 1 μm として計算し，初期の平均粒径は 3.19 mm である．境界条件は周期的えたときの粒成長挙動を確認した．

境界条件を用いた．上記シミュレーションは pythonコー（12）
ここで，f0は平均の体積分率，mはピンニング粒子の偏ドおよび C言語によるコードとして計算を行った．計算に用いたパラメータは Table 1のとおりである．

り度合いを表す整数，Ymaxは y方向のセルの数である．

3．結果および考察

本研究では式（13）に示すように，∆Gが低下するかつ，まず，易動度を強調せず，析出物がない場合のシミュ

粒成長速度 vが乱数で与えられるある z 値以上であれば，レーション結果を示す．初期組織，100，500，1000 ス

前述のように粒成長速度 vは式（1）のように表される．

∆Gが最も低下するように，隣接する粒と同じ相になるというルールで粒成長過程をシミュレーションする．

テップ後の組織をそれぞれ Fig. 4 に示す．結晶方位はIPF（Inverse Pole Figure）マップで示している．また，平均粒径の推移を Fig. 5 に示す．ここで，点線に示す放物則は式（14）のように表される．

（13）

（14）

ここで ∆G [x, y]-[x+i, x+j] ，はセル [x, y]のみがセル [x+i, x+j]

ここで tは時刻，t0は基準となる時刻，<Rg>は時刻 tでのafterの粒に置き換わったとしたときのセル x, yにおける ∆G

平均粒径，<R0>は初期の平均粒径，αは比例定数である．

ステップごとに更新される．計算全体の流れは Fig. 3にが開始する時刻を t0 として，最小二乗法によるフィッである．z は最大を z maxとする一様乱数であり，各セル，示すフローチャートのとおりである．

Fig. 5 では，初期粒径から計算が落ち着いて，粒成長ティングを行った直線を破線で示している．
平均粒径は Fig. 5 に示されるとおり，放物則に近い推7

技術論文＞セルラーオートマトン法を用いたオーステナイト粒の異常粒成長シミュレーション

Table 1. Parameters for cellular automaton.
SymbolValueQ

1.82 × 10R8.314

M04.67 × 10

E00.75TzmaxbXmax, Ymax

Initial state

6 27)−6 27)

28)

Unit

Description

J/mol

Activation energy

J/(K·mol)Gas constantm /(J·s)

Coefficient of grain boundary mobility

43J/m

Grain boundary energy

1173

KTemperature

5 × 10m/s

Maximum random number

7pixelRange of cells in the box count method

150, 150

pixel

Size of calculation cells

Step: 100−7Step: 500

Step: 1000111Fig. 4. Microstructures during specific steps.

100110移を示しており，定常粒成長であると考えられる．
次に，易動度を強調する粒の易動度を 10 倍，強調する粒の数を 10 個とした場合の粒成長挙動を Fig. 6 に示す．易動度を強調したことで，異常粒成長が認められる．易動度を強調したことは，対応粒界や偏析の影響を導入したことに相当すると考えることができる．例えば
Olmsted らは対応粒界の易動度が他の粒界よりも非常に高いことを報告している 29）．また，偏析により易動度が低下するとの報告もある 30）．
異常粒成長の発生の基準については，組織内での最大粒径 Rmax を平均粒径 <Rg> で割った Z g を基準とし，この値が 3 以上のとき異常粒成長とする．

この時の組織の一例を Fig. 7 に示す．Z g が 3 より大きい時に異常粒成長が発生していると考えてよい（式Fig. 5. Time evolution of average grain radius.
上で，異なる 10 の乱数シードを CA 法による計算の前に与えてそれぞれ CA 法による計算を行い，Z g をプロッ（15））．

トしてある．強調する粒の選択と，それぞれのステップ（15）ごとの変化は乱数に依存するため，結果が乱数シードに易動度を強調度合いを表す係数 p の影響を Fig. 8 に示

異常粒成長は易動度の係数が 5 以上の時に発生しはじより異なってくることとなる．

す．x 軸に強調の係数 p を，y 軸に 1000 ステップ後におめており，7 以上では乱数シードによらずすべての場合の数は 10 個である．次に，異常粒成長におよぼす強調

でも異常粒成長が発生するが，初期組織にある粒の数のける Z g を示している．このとき，易動度を強調する粒で発生している．また，易動度を強調する粒の数は 1 つする粒の数の影響を Fig. 9 に示す．x 軸に強調する粒の内，強調する粒の数が 20 % を超えると異常粒成長が発Fig. 8 および Fig. 9 では乱数シード以外の条件を揃えたと，強調された粒のみで定常粒成長となるためである．

数，y 軸に 1000 ステップ後における Z g を示している．

生しない．これは，易動度を強調する粒の数が多くなる

電気製鋼 第 95 巻第 1 号 2024 年

8Initial state

Step: 100

Step: 500

Step: 1000
111

110100Fig. 6. Microstructures during specific steps in the case of enhanced mobility.

p=3Zg = 2.3

p=4Zg = 2.8

p=5Zg = 3.6

p=6Zg = 3.8111

Normal grain growth

Abnormal grain growth100110Fig. 7. Effect of mobility coefficient p on the microstructure after 1000 steps.
Fig. 10 に 1000 ステップ後における計算組織に対するThe percentage of enhanced grains (%)

析出物の影響を示す．また，Z g に対する析出物の体積Random seed 0Random seed 1分率 f0 による影響を Fig. 11 に示す．これらの結果は析

Random seed 2Random seed 3

出物の体積分率が大きくなると，異常粒成長の発生を抑Random seed 4Random seed 5制することが可能であることを示している．ここでは析

Random seed 6Random seed 7

出物の半径 r は位置によらず，5 nm で一定とした．まRandom seed 8た，易動度を強調する粒の易動度を 10 倍，強調する粒Random seed 9の数を 10 個としている．析出物の体積分率に偏りがな

いとき，析出物の平均体積分率 f0 が約 6.0 × 10-5 以上であれば，異常粒成長を抑制することができる．
ここで，Gladman の式 31） との対応関係について議論する．Gladman によると異常粒成長発生の臨界析出物半The number of enhanced grains

The percentage of enhanced grains (%)

Random seed 0
Random seed 1
Random seed 2
Random seed 3
Random seed 4
Random seed 5
Random seed 6Zg

Random seed 7
Random seed 8
Random seed 9

Mobility coefficientFig. 8. Effect of mobility coefficient on Z g after 1000steps. (Dotted line is abnormal grain growthcriteria in Equation (15).)

The number of enhanced grains

Fig. 9. Effect of the number of enhanced grains on Zg
after 1000 steps. (Dotted line is abnormal grain
growth criteria in Equation (15).)

技術論文＞セルラーオートマトン法を用いたオーステナイト粒の異常粒成長シミュレーション

f0 = 2.4×10-4
Zg = 4.8

f0 = 4.8×10-4
Zg = 3.6

f0 = 7.2×10-4
Zg = 2.3

9f0 = 9.6×10-4
Zg = 1.9111

Abnormal grain growth100

Normal grain growth110

Fig. 10. Effect of the volume fraction of precipitates on the microstructures after 1000 steps.
径 r* は式（16）のとおり表される．

Random seed 0Random seed 1
Random seed 2（16）Random seed 3
Random seed 4
Random seed 5

ここで Z は初期の Zg，ピンニング粒子の体積分率を，Random seed 6Random seed 7

異常粒成長を抑制できる f0 ＝ 6.0 × 10-5 としたとき，r* ＝Random seed 80.99 nm と計算される．今回計算に用いた析出物粒径はRandom seed 95 nm であるから，Gladman の式より求められる値よりも約 5 倍大きくなったが，オーダーとしては近い値が得られており，実組織に近い挙動を再現できると考えられる．
析出物の偏りについて，析出物が Fig. 12 の上段のように，偏って存在しているとした場合の 1000 ステップ後の組織を Fig. 12 の下段に示す．m が 2 以下の時（偏りのない m=0 の時を除く）に，異常粒成長が発生した．

Volume fraction of presipitates f0Fig. 11. Effect of the precipitates volume fraction on Zg
after 1000 steps. (Dotted line is abnormal grain

加えて，析出物の体積分率が低い部分で異常粒成長が発

growth criteria in Equation (15).)生している．このように，析出物の体積密度に偏りがあれば，異常粒成長の発生が起こりやすくなる．特に，偏で考慮されてこなかった要因について導入した粒成長シりが大きいほど，異常粒成長が発生しやすくなる傾向にミュレーションモデルを提案した．こうした要因も導入ある．このことは，析出物分布の偏りが大きいほど異常したシミュレーションを行うことで，異常粒成長の予測粒成長の抑制には不利であることを示す．

が可能となる．組織制御が容易になることで材料開発へ

以上のように，易動度や析出物の偏りといったこれまm=0m=2m=4の知見が得られると考える．
m=6m=8f0×10-41.75
1.50
1.25
1.00
0.75
0.50
0.25m=0Zg = 2.4m=2Zg = 4.6m=4Zg= 4.0m=6Zg= 3.0m=8Zg= 2.8111

100Fig. 12. Effect of the precipitates volume fraction on the microstructures after 1000 steps.

11010電気製鋼 第 95 巻第 1 号 2024 年

4．結言

本研究では，結晶粒方位に実方位を用いた CA 法を用いて，異常粒成長の発生および，析出物による異常粒成長の抑制を再現することを試みた．その結果以下の結論が得られた．
1．CA法において，異常粒成長は粒界の易動度を強調することで発生し，強調の程度を 5倍以上とすると異常粒成長が発生する．このことから粒界の易動度の不均一性は異常粒成長発生の一要因であると示唆される．
2．粒番号ではなく，結晶方位を用いることで，粒界エネルギーと易動度の影響を正確に再現することが可能である．
3．異常粒成長は析出物のピンニング効果により抑制が可能である．しかし，組織の中で析出物の体積分率に偏りがあり，ピンニング効果が弱い場所がある場合には，異常粒成長の抑制が難しくなる．
4．以上の結果から，CA法は異常粒成長の再現に用いることが可能であり，浸炭処理などの熱処理中の異常粒成長挙動を再現するのに有効な手段となりうる．

16）紅林豊，中村貞行：電気製鋼，65（1994），67．

17）T. Kinoshita and M. Ohno：Comput. Mater. Sci.，177（2020），109558.

18）R. J. Contieri，M. Zanotello and R. Caram：Mater．Res.，
20（2017），688.

19）P. J. P. Pimienta，E. J. Garboczi and W. C. Carter：
Comput. Mater. Sci.，1（1992），63.

20）高橋裕，宮沢薫一，森実，石田洋一：日本金属学会
誌，50（1986），357.

21）D. C. Handscomb：Can. J. Math. ，10（1958），85.

22）W. T. Read and W. Shockley：Phys. Rev.，78（1950），275.

23）A. Mallick and S. Vedantam：Comput. Mater. Sci.，46（2009），21.
24）F. J. Humphreys：Acta Mater.，45（1997），4231.

25）C. S. Smith：Trans. Metall. Soc. AIME，175（1948），15.

26）大笹憲一，棗千修：鉄と鋼，102（2016），157.

27）S. Dépinoy，B. Marini，C. Toffolon-Masclet，F. Rochand A.-F. Gourgues-Lorenzon：Metall. Mater. Trans. A，48
（2017），2289．

28）H. S. Zurob，Y. Brechet and G. Purdy：Acta Mater.，49（2001），4183．
29）D. L. Olmsted，E. A. Holm and S. M. Foiles：Acta（文献）1）瓜田龍実，並木邦夫，飯久保知人：電気製鋼，59（1988），33．
2）O. M. Ivasishin，S. V. Shevchenko，N. L. Vasiliev andS. L. Semiatin：Acta Mater.，51（2003），1019.

3）Y．Saito and M．Enomoto：ISIJ Int．
，32
（1992），267.

Mater.，57（2009），3704.

30）J. W. Cahn，P. Fife and O. Penrose：Acta Mater.，45（1997），4397.
31）T. Gladman：Proc. R. Soc. Lond. A，294（1966），298.
本 論 文 は，「Matrials」17 巻， 記 事 番 号 138 に 掲 載4）H. V. Atkinson：Acta Metall.，36（1988），469.

さ れ た“Simulation of Abnormal Grain Growth Using the6）Y. J. Lan，D. Z. Li and Y. Y. Li：Metall. Mater. Trans.

Fei Sun, Ta-Te Chen and Yoshitaka Adachi（2024）（CC BY5）Q. Yu and S. K．Esche：Mater. Lett.，57（2003），4622.
B，37（2006），119．

7）F. Han，B. Tang，H. Kou，J. Li and Y. Feng：J. Mater.
Sci.，48（2013），7142.

Cellular Automaton Method”©Kenji Murata, Chihiro Fukui,
4.0, https://creativecommons.org/licenses/by/4.0/）[https://doi.org/10.3390/ma17010138] の内容を翻訳し，加筆修正したものである．

8）Y. He，H. Ding，L. Liu and K. Shin：Mater. Sci. Eng.
A，429
（2006），236.

9）J. Ogawa and Y. Natsume：Comput. Mater. Sci.，199（2021），110729.
10）L. -Q. Chen and W. Yang：Phys. Rev. B，50
（1994），15752.

11）C. E. Krill III and L. -Q. Chen：Acta Mater.，50（2002），3057.

村田憲治

福井ちひろ陳足立吉隆

12）Y. Hayakawa and J. A. Szpunar：Acta Mater.，45（1997），1285.

13）Y. Hayakawa and J. A. Szpunar：Acta Mater.，45（1997），4713.
14）L. Ye，B. Mei and L. Yu：Metals，12（2022），1717.

15）木下修司，上田武司，鈴木章：鉄と鋼，59（1973），446．

達徳

孫飛