JP5034194B2 - シミュレーション方法 - Google Patents

Description

本発明は、ガス、その他雰囲気中での加熱による処理において、材料内の粒子の挙動をシミュレーションする方法に関する。

近年、浸炭などの表面処理では、外部雰囲気からの拡散により母相を硬化するのみならず、拡散層内において炭化物などの硬質粒子を形成、成長させて、はるかに高い硬さを達成している。拡散層の硬さ制御のためには、侵入元素の濃度プロファイルとともに、拡散粒子の粒径や体積率に関する情報がきわめて重要である。

実用鋼など多成分系では、各元素の相互拡散の扱いが複雑であり、熱力学平衡と収支バランスとを同時に満たす収束計算が困難である。現状では市販ソフトＤＩＣＴＲＡが唯一、上記問題を１次元でのみ数値解析することができる。

ＤＩＣＴＲＡは多成分系の１次元拡散を数値解析するソフトであり、例えば浸炭雰囲気を境界条件として、単相の合金中へ一方向に拡散するＣ（炭素）の挙動をシミュレーションできる。また、２相以上が隣接するモデルにおいて、異相界面の熱力学平衡および収支を計算することができ、温度変化に伴う固−液変態、固−固変態など拡散律速現象のシミュレーションをその対象としている。

ＤＩＣＴＲＡの単相の一方向シミュレーションでは、拡散層内の各深さで全元素濃度が刻々計算されるため、それに対応して熱力学平衡計算を行えば、粒子の平衡モル分率のみ求めることができる。しかし、実際にはその瞬間に平衡モル分率に達する訳ではなく、また、あくまで母相のみのモデルなので、粒子の粒径等に関しては計算できない。

一方、２相以上を扱うモデルでは、基本的に温度以外の変化がない閉鎖熱力学系が、その時々の温度で熱力学的に最も安定な相比率と濃度分配に向かって、拡散による界面移動を起こす様をシミュレーションできる。これは、例えば、室温で安定であった粒子が高温で量的に不安定となり、安定な相比率となるべく母相中へ固溶していく過程など、材料内の自発的な拡散現象を扱うものである（例えば、非特許文献１参照）。

すなわち、ＤＩＣＴＲＡでは、１次元でモデル化できるシンプルな現象に限って、それぞれ別個のモデルが存在しているが、各モデルを連携させた解析は想定されていない。このため母相と粒子の２相があり、さらに外部からの拡散流速が作用するような場合には、該粒子について妥当な境界条件を知る手段がなく解析できなかった。

さらに高度な解析としては、表面から深さ方向へ複数の粒子を配置し、それぞれの地点での流速変化をバランスさせるシミュレーション方法も考えられる。しかし、実用金属材料の多成分系においては計算収束の困難さから現実的でなく、ＤＩＣＴＲＡに限らずその他のソフトウェアを用いても実施できなかった。
Ａ．Ｂｏｒｇｅｎｓｔａｍ，Ａ．Ｅｎｇｓｔoｍ，Ｌ．Ｈoｇｌｕｎｄ ａｎｄ Ｊ．Aｇｒｅｎ"Ｄｉｃｔｒａ，ａ ｔｏｏｌ ｆｏｒ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｄｉｆｆｕｓｉｏｎａｌ ｔｒａｎｓｆｏｒｍａｉｔｉｏｎ ｉｎ ａｌｌｏｙｓ" Ｊ．ｐｈａｓｅ Ｅｇｕｉｌｉｂｒｉａ，ｖｏｌ．２１，Ｎｏ．３，２０００，ｐ２６９

本発明は、上記従来の課題を解決することを目的とする。すなわち、本発明は、材料を加熱処理する際に、拡散層内の所定深さにある粒子の粒径の経時変化を精度よく予測できる粒子のシミュレーション方法を提供することを目的とする。

上記課題は、下記本発明のシミュレーション方法により解決することができる。

すなわち、本発明のシミュレーション方法は、材料を前加熱を行った後に行う表面処理で、粒子の挙動として、前記材料の拡散層内の所定の深さにおける、当該粒子の粒径の経時変化を予測するシミュレーション方法であって、前記表面処理前の状態で、予め存在する粒子とその粒子の周辺にある母相とからなる球状モデルを作成する第１の過程と、前記材料表面から侵入する侵入物質の活量と時間との活量経時関係を求める第２の過程と、前記第２の過程で得られた活量経時関係を境界条件として、前記第１の過程により作成された前記球状モデルから、粒子の挙動を予測することを特徴とする。

拡散層内のある深さに位置する粒子の粒径等を計算予測するには、（１）マクロな流速と粒子まわりの局所流速との関係、および（２）母相との界面での熱力学平衡と収支バランスとを考慮する必要がある。そこで、本発明の粒子のシミュレーション方法では、まず、第１の過程として、所望の深さに存在する粒子を考え、単一の粒子とその周辺の母相からなる球状モデルを作成する。球状モデル内の濃度の初期条件は、例えば、粒子内および母相内で濃度均一として行うのが良い。

また、第２の過程として、まず外部雰囲気を境界条件に、表面から侵入する元素の全体的な拡散流速を、平面的な１次元モデルで計算する。その結果をもとに、所望の深さを通過する流速の時間変化をとらえ、侵入元素の活量の時間変化として出力する。結果は、数値データとして保存するか、適当な近似関数（活量経時関係）を与えておく。

そして、第２の過程における平面１次元モデルによる拡散解析から抽出した活量経時関係を境界条件として、第１の過程における球状モデルに当てはめることで、拡散層内の所望の深さにある粒子と周辺母相が受ける拡散流速とを間接的ながら再現することが可能となり、拡散層内の所定深さにある粒子の粒径等の経時変化を精度よく予測することができる。

本発明のシミュレーション方法においては、前記第１の過程における前記球状モデルが、前記前加熱を行った直後の前記粒子の濃度分布をシミュレーションして作成されたものであることが好ましい。すなわち、加熱前の粒径や均一な濃度をそのまま用いるよりも、前加熱が完了して外部から拡散流速が侵入する直前の状態をシミュレーションした結果を用いることで、より現実的に、粒子のシミュレーションを行うことができる。

また、本発明のシミュレーション方法においては、前記粒子の挙動が、前記拡散層内の所定の深さにおける、当該粒子の粒径の経時変化である。

さらに、本発明のシミュレーション方法においては、前記表面処理が浸炭処理であり、前記粒子がセメンタイト粒子であることが好ましい。

本発明のシミュレーション方法によれば、材料を加熱処理する際に、拡散層内の所定深さにある粒子の粒径の経時変化を精度よく予測できる。

本発明のシミュレーション方法は、加熱雰囲気中での処理によって材料表面に形成される拡散層内の粒子の挙動（粒径の経時変化等）を予測するものである。そして、種々の材料および種々の加熱処理の場合に適用することができる。例えば、各種炭化物粒子の挙動や、窒化処理における窒化物の変化、あるいは、浸漬、ＰＶＤなどの特定の元素をコーティングした後、加熱により拡散浸透させる場合の材料内部粒子の変化を精度よく予測することができる。

以下、高濃度の浸炭（ＣＤＣ）処理におけるセメンタイト（以下、「θ」と表すことがある）粒子の挙動に、本発明の手法を適用した例を説明する。

本発明のシミュレーション方法によるＣＤＣ処理の解析にあたり考慮した材料の相変態およびθの変化を図１に示す。図１に示すように、まず、所定温度への前加熱によって材料の母相がα／γ変態を起こすとともに、成分によってはθが固溶し始める。次に、浸炭性雰囲気が整うと浸炭が始まり、拡散していきたＣが粒子に到達するとθは逆に成長をし始める。

このような現象に対し、図２及び図３に示す手順（第１の過程および第２の過程）で、浸炭層内のθ粒径の変化をシミュレーションする。個別のモデル作成、拡散の計算等にはＤＩＣＴＲＡを用いることが好ましい。以下、上記手順について説明する。

＜第１の過程＞
第１の過程では、前加熱を行う前の状態で、侵入物質（炭素）を含む粒子とその粒子の周辺にある母相とからなる球状モデルを作成する。すなわち、球状化処理材などの球状炭化物を想定して、図２に示すような、中心にθ、外周に母相であるフェライト（以下、「α」と表す）を、その間に高温で相変態により安定となる母相であるオーステナイト（以下、「γ」と表す）を潜在的に配置した球状モデルを用いる。γは、α／θ界面近傍の計算点が、１×１０^-5Ｊ／ｍｏｌの熱力学的駆動力を持ったときに生成すると仮定する。加熱速度・保持時間など前加熱時の熱履歴を設定して、母相のαからγへの変態、炭化物の固溶をシミュレーションし、最終状態を保存する。

具体的な計算の流れは図１１に示すとおりで、具体的には下記（１）〜（６）のようになる。

（１）対応する母相（α）の半径計算：
まず、球状モデルは、材料中に予め存在する分散粒子の半径ｒθと体積率ｖθとから、１個の粒子θと周囲の母相αからなるモデルとして作成する。ここで、モデル全体の半径をｒとすると、ｒ＝ｒθ＋ｒαで、ｖθ＝（ｒθ／ｒ）³の関係が成立する。従って、対応する母相の半径は、「ｒα＝ｒθ／（ｖθ）^1/3−ｒθ」の式で計算し、シミュレーション領域のｒαおよびｒθを求める。

（２）θとαとの平衡計算：
θとαとの初期成分濃度を求める。実測しても良いが、θとαとに熱力学的平衡が成り立つと仮定し、各成分の化学ポテンシャルμが両相で等しくなる（すなわち、θ＋αの２相混合系の化学エネルギーＧが最小となる）濃度を数値的に解析する。成分ｉのμ_iは各相のＧのｉ濃度に関する偏分であり、各相のＧは世界標準のデータベースから供給される。ここでは、Ｎｅｗｔｏｎ−Ｒａｐｓｏｎ法によりＧの最小化を行う市販ソフトＴｈｅｒｍｏ−Ｃａｌｃを用いて、５００℃での平均濃度を計算した。

（３）温度および時間ステップの設定：
具体的には、前加熱過程の温度変化を時間の関数Ｔ（ｔ）で設定する。本実施例の被処理材は、温度が上昇すると母相がαからγへ相変態するため、球状モデルのθとαとの界面の濃度から上記の熱力学平衡計算を逐次行い、界面でγが出現するための熱力学的駆動力を持つかどうかを判定する。

（４）θとαとの界面、およびαとγとの界面の成分濃度探索：
γが出現するとθ／γ／αの３相モデルとなり、θ／γ、γ／αの２つの界面が存在する。界面において隣りあう２相にそれぞれ熱力学的平衡が成り立つとし、θ側／γ側およびγ側／α側について成分ｉの界面濃度を（２）の方法で計算する。このとき（２）と違うのは、界面だけの局所的な平衡なので材料全体の平均濃度の制約がなく、平衡する界面濃度は何通りもある。いずれにしても、界面に位置する計算点のＧおよび濃度から、界面のμ_iが求まり、３相について各相内で各成分のμ_iの勾配を計算できる。

（５）θ、α、γ内の成分流量の計算：
各相内での成分ｉの流量Ｊ_iは、μ_iの勾配に比例する。比較係数である拡散係数は、全成分の相互作用を含むデータベースから供給される。

（６）θ、γおよびαのそれぞれの間の物質保存：
各相内を流れ、単位時間あたりに界面を通過した成分の流量Ｊ_iが、θ／γ、γ／αの２つの界面の移動ｖを起すはずなので、Ｊ_iとｖ×（成分ｉの界面濃度の差）とのバランスを計算する。差分が大きければ（４）にもどって平衡する界面濃度を再探索する。この流量バランスが全成分について成立するように、Ｎｅｗｔｏｎ−Ｒａｐｓｏｎ法で収束させる。差分が収束すれば界面の移動量が決定され、新しい分散粒子の半径ｒθが求まる。次の時間ステップに進む。このようにして、加熱によって変化する母相およびθの状態がシミュレーションされるので、加熱が完了した最終状態を浸炭直前の状態として保存する。

＜第２の過程＞
第２の過程では、材料表面から侵入する侵入物質の活量と時間との活量経時関係を求める。すなわち、まず、平面モデルを用いて等温浸炭による一方向の拡散をシミュレーションし、拡散層全体にわたる最終的な濃度プロファイルを得る。これを参考に粒径を予測したい深さ位置を決め、その位置でのＣ活量の時間変化を出力する。出力したデータを数値解析し、指数関数もしくは多項式でフィッティングして、活量経時関係としての式を導き出す。

具体的な計算の流れは図１１に示すとおりで、具体的には下記（７）〜（１０）のようになる。

（７）境界条件の計算：
実際の操業で外部のガス雰囲気は、例えば浸炭処理ではカーボンポテンシャルＣＰで制御される。Ｆｅ−Ｃ系のＣ濃度、温度と活量との関係図から、ＣＰを熱力学量であるａ_cに変換する。

（８）温度、時間ステップ設定：
ａ_cを境界条件として被処理材の化学成分を均一に有する平面モデルを作成する。母相γと粒子θがあることを想定しているが、θ内の拡散は無視する。処理（浸炭）温度を時間の関数Ｔ（ｔ）で設定する。被処理剤の表面と内部との化学ポテンシャルμ_iが、左記と同様にＧのデータベースから計算されるので、その勾配を求める。

（９）成分流量の計算：
μ_iの勾配が、ガス成分ｉの流量Ｊ_iとなって、内部へ拡散していく。（５）と同様に流量を求め、表面からの深さＺに対するｉの濃度プロファイル、ａ_cなどを計算する。

（１０）各深さでのγとθとの平衡計算：
内部で成分ｉが増加すると、熱力学的にはｉを主成分とするθが増加する。各計算点の成分濃度に対し、（２）の熱力学平衡計算を行い、θに対するγの平衡濃度を算出する。これが、新しい母相の濃度で、次の時間ステップに進む。

次に、図３に示すように、先に保存した前加熱終了時の球状モデルを再び読み込み、これに境界条件として上記のＣ活量変化を時間関数で与え、所定の深さでの炭化物成長をシミュレーションする。

具体的な計算の流れは図１１に示すとおりで、具体的には下記（１１）〜（１５）のようになる。

（１１）所望深さのＣ活量の時間変化：
上記の結果から、拡散層内でシュミレーションしたい粒子の位置Ｚにおける、ａ_cの時間変化を出力し、簡便に入力できるように関数近似を行う。

（１２）温度、時間ステップ設定：
当該設定は、上記（３）と同様にして行う。

（１３）θ／γ界面の成分濃度探索：
当該探索は、上記（４）と同様にして行う。但し、上記（４）と違うのは、前加熱過程で処理温度に達し、母相がすべてγへ相変態した状態を引き継いでいるため、球状モデルはθとγの２相のみで、界面もθ／γの１つだけであることである。

（１４）θ、γ内の成分流量の計算：
θ、γ内の成分流量の計算は上記（５）と同様にして行う。

（１５）θとγとの間の物質保存：
θとγとの間の物質保存は上記（６）と同様にして行う。

（シミュレーション例）
一般軸受鋼であるＳＵＪ２（Ｆｅ−１．５Ｃｒ−１．０Ｃ）球状化処理材についてシミュレーションした。なお、球状化処理材でのθ半径は０．２２μｍで既知とし、体積率から計算した球状モデルの半径は０．４０２μｍであった。

図４に、上記モデルに加熱前処理過程として５００℃より１℃／ｓで８５０℃に昇温、０．５ｈ保持する温度条件を与えたときのθ半径の変化を、また、図５に、θ粒子を中心とした最終的なＣおよびＣｒの濃度プロファイルを示す。

これから、前加熱を終えたときθはわずかずつ固溶しつつある段階で、θ粒子／母相の境界付近でＣｒが大きな濃度勾配を持っていると予想される。浸炭の直前、θ粒子まわりはこのような遷移的状態であり、これをそのまま浸炭過程の初期状態とすることで、より現実に即したシミュレーションが可能となる。

図６に、平面モデルによるＳＵＪ２の５ｈ浸炭後のＣ濃度プロファイルのシミュレーション結果を、実際の浸炭材での分析結果と比較して示す。モデル長さは２ｍｍで、左端を材料表面として実際の浸炭条件にならい、境界条件は８５０℃、ＣＰ（カーボンポテンシャル）＝１．２９である。なお、ＣＰ＝１．２９は、グラファイトを基準としたＣ活量が１．０になる条件でＣＤＣ処理の典型的なＣＰ値である。

計算したプロファイルはほぼ分析波形の平均的な値になっているので、拡散層内の全体にわたって実際に近いＣ流速が再現できていると思われる。ここで、Ｃ濃度＝１．３ｗｔ％の位置を深さ約２００μｍの位置と見なし、浸炭中における２００μｍ位置のＣ活量の時間変化を１００μｍ、３００μｍ位置とともに出力したのが、図７である。２００μｍ深さの結果を解析し、活量と時間の関係として、下記式（多項式）を得た。なお、下記式中の「ｔ」は、時間（秒）を表す。

（式）：「Ｃの活量」＝０．７２−２．０２×１０^-14×ｔ³−４．５９×１０^-10×ｔ²−３．３７×１０^-6×ｔ

図８および図９は、先に保存した前加熱完了後の球状モデルに、上記のＣ活量に関する境界条件を与えたシミュレーションによるθ粒子径、および体積率の変化を示す図である。実際の浸炭材で、分析にもとづいてＣ濃度＝１．３ｗｔ％の位置を切断し、測定した値を合わせて示す。測定値は２００μｍ深さでの計算結果に最も近い。切断位置や分析誤差を考慮すれば、このように実用的な精度で、所定の位置の粒子径および体積率が予測可能である。

一方、比較として、図１０は、図６で計算された各深さでのＣ濃度に対して、時間ごとに熱力学計算を行い、得られたθの平衡モル分率の変化を体積率に変換したものである。実際の浸炭材で見られたθの量は、この平衡体積率よりもはるかに小さい。すなわち、濃度が刻々変化する浸炭過程において、θが熱力学的に最も安定な量に達することは到底ないのであって、平面モデルと熱力学計算だけの従来法では、浸炭層内のθ粒子径や体積率を予測することはできない。

高濃度の浸炭処理の解析で考慮した材料中の変化を示す図である。 本発明のシミュレーションの流れを説明する図である。 本発明のシミュレーションの流れを説明する図である。 前加熱処理におけるセメンタイトの粒径変化を示す図である。 前加熱終了直後（浸炭直前）における炭素とクロムの濃度プロファイルを示す図である。 浸炭層内の各深さにおける炭素活量の時間変化を示す図である。 浸炭層内の各深さにおける炭素活量の時間変化を示す図である。 浸炭層内の各深さにおけるθの粒子径変化を示す図である。 浸炭層内の各深さにおけるθ粒子体積率変化を示す図である。 炭素濃度変化に対応したθ粒子の平均体積率を示す図である。 本発明のシミュレーションの流れを説明する図である。

Claims (3)

1. 材料を前加熱を行った後に行う表面処理で、粒子の挙動として、前記材料の拡散層内の所定の深さにおける、当該粒子の粒径の経時変化を予測するシミュレーション方法であって、
   前記表面処理前の状態で、予め存在する粒子とその粒子の周辺にある母相とからなる球状モデルを作成する第１の過程と、
   前記材料表面から侵入する侵入物質の活量と時間との活量経時関係を求める第２の過程と、
   前記第２の過程で得られた活量経時関係を境界条件として、前記第１の過程により作成された前記球状モデルから、前記粒子の挙動を予測することを特徴とするシミュレーション方法。
2. 前記第１の過程における前記球状モデルが、前記前加熱を行った直後の前記粒子の濃度分布をシミュレーションして作成されたものであることを特徴とする請求項１に記載のシミュレーション方法。
3. 前記表面処理が浸炭処理であり、前記粒子がセメンタイト粒子であることを特徴とする請求項１又は２に記載のシミュレーション方法。

Images