# Fileset

[669.pdf](https://mdr.nims.go.jp/filesets/77a3bd22-4132-4221-b67a-12546f192257/download)

## Creator

[石井 秋光](https://orcid.org/0000-0002-9261-4047)

## Rights

[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[High-fidelity Prediction of Solid-state Sintering Through the Integration of Experiments and Numerical Simulations](https://mdr.nims.go.jp/datasets/eb771987-cec1-4191-a172-f1b62e47fee8)

## Fulltext

物質・材料研究機構構造材料研究センター材料評価分野計算組織設計グループ研究員(〒3050047 つくば市千現 121)Highfidelity Prediction of Solidstate Sintering Through the Integration of Experiments and Numerical Simulations; Akimitsu Ishii(National Institute for Materials Science, Tsukuba)Keywords: sintering, simulation, phasefield method, data assimilation, in situ observation, nanoparticle2025年 6 月30日受理[doi:10.2320/materia.64.669]ま て り あMateria Japan第64巻 第10号(2025)実験と計算の融合による固相焼結挙動の高精度予測石 井 秋 光. は じ め に焼結は，粉末状の原料を加圧成型した後，融点以下の温度で熱処理することで多結晶材料を作製する技術である．焼結は比較的シンプルなプロセスであり，原料粉末さえ用意できれば様々な配合や形状で材料を作製できるため，粉末冶金やセラミックス分野のみならず，全固体電池，永久磁石，超伝導材料などの材料作製にも用いられている．近年では，金属3D プリンターの基盤技術としても注目を集めている．焼結によって作製された材料は多結晶であるため，材料の機械的，電気的，磁気的特性は内部のミクロ組織に強く依存する．焼結過程では，加圧成型体の相対密度が加熱時間の経過とともに上昇する「緻密化」と呼ばれる特有の現象が発生する．典型的には，加熱前の加圧成型体の相対密度は50～60程度であり，十分な時間の熱処理を行えば相対密度はほぼ100まで上昇する．しかしながら，材料の種類や組成，粒成長抑制とのトレードオフによる熱処理時間への制約などの要因により，相対密度が100まで上昇しない場合も多い．すなわち，材料内部に気孔が残留し，この気孔の有無や量は材料特性に強く影響を及ぼす．したがって，ミクロ組織の形成・変化および緻密化挙動を正確に予測・制御することは，焼結部材の高特性化を実現するために不可欠である．実験観察技術の進歩により，固相焼結挙動を 3 次元かつ非破壊でその場観察することが可能となっている(1)(2)．一方，実験のみでは明らかにすることが難しいメカニズムの解明や，任意の条件下での焼結挙動を予測するために，数値シミュレーションを用いた研究も多数行われている．特に，複雑な界面移動を比較的容易かつ精度良く予測可能なフェーズフィールド(Phasefield: PF)法(3)を用いた焼結シミュレーションは，この20年で活発に研究されてきた(4)(6)．しかしながら，PF 法は連続体モデルに基づく手法であるため，高精度な予測結果を得るためには正確な材料の物性値やパラメータ(原子拡散係数，界面エネルギー，界面易動度など，以下ではこれらを総称して材料パラメータと呼ぶ)が必要である．材料パラメータが不正確な場合，予測精度は低下し，そもそも計算に必要な材料パラメータがすべて揃っている材料は限られている．材料パラメータは，実験や第一原理計算，または分子動力学シミュレーションなどから同定することができるが，すべての未知の材料パラメータを実験や計算で取得するのは非常に時間のかかるタスクであり，現実的ではない．そこで筆者らの研究グループでは，固相焼結に「データ同化」(7)を適用することで，上記課題の克服を目指した．データ同化とは，実験データと数値シミュレーションをベイズ推論に基づいて融合することで，数値シミュレーションの予測精度の向上や未知パラメータの推定が可能なデータ科学技術である．筆者らは，Cu ナノ粒子を対象として固相焼結のその場観察および PF シミュレーションを実施し，独自に開発したデータ同化手法を用いて両者を融合した．その結果，複数の未知材料パラメータの同時推定を通じて，実験で観察された固相焼結挙動を PF シミュレーションで忠実に再現することに成功した．本稿では，この筆者の最新の研究成果(8)について，実験・シミュレーション・データ同化のそれぞれの手法と，得られた成果について紹介する．. Cu ナノ粒子固相焼結の実験観察データ同化を行うためには，時系列の実験観察データが必要となる．また，実現象を予測するためには，3 次元のデータを得ることが望ましい．そこで本研究では，走査透過電子顕微鏡(Scanning transmission electron microscopy: STEM)を用いて，約20個の Cu ナノ粒子からなるクラスターについての固相焼結中の 3 次元的な時系列変化をその場観察し図 1 (a) STEM を用いた観察により得られた Cu ナノ粒子焼結の時系列 3 次元データ．(b) データ同化用にノイズ除去を施した後の時系列データ． 　　　　　　新 進 気 鋭た．焼結が進行しない200°Cを待機温度とし，この温度下でTEM ホルダーを－62°から62°まで 2°ごと傾斜させながら観察することで，傾斜角ごとの高角度散乱暗視野(Highangleannular darkfield: HAADF)像を取得した．その後，傾斜シリーズデータを 3 次元再構成した．待機温度での観察と焼結の進行を交互に行い，焼結を進行させる際はターゲット温度まで急速加熱して 5 秒保持した後，待機温度まで急冷した．焼結は55秒行い，ターゲット温度は370°Cから420°Cまで段階的に上昇させた．さらに詳細な実験条件については文献(1)(8)を参考にされたい．図(a)には，3 次元再構成した時系列の実験観察データを示す．低温かつ短時間の加熱であるため，著しい緻密化や粒成長は見られないが，初期状態で存在した粒子間の空隙の縮小・消失が確認できた．観察された Cu ナノ粒子(コア粒子)の平均粒径は 150 nm であった．粒子の表面は，低温焼結を促進させる焼結助剤として機能する 10 nm 以下の微細粒子(シェル粒子)で覆われている(1)が，今回のその場観察では，コア粒子のみが観察された．データ同化には，画像処理技術を応用したノイズ除去を施したボクセル形式のデータ(図 1(b))を使用した．. 固相焼結フェーズフィールドモデル固相焼結の PF モデルはこれまでに数多く提案されているが，本研究では最も基礎的でシンプルながらも，焼結中に生じるミクロ組織変化と緻密化を同時に予測することが可能なWang(5)の固相焼結 PF モデルを使用した．使用した PF モデルについて以下に概説する．解析対象とする系内部の任意の座標の状態を表す変数として，焼結体の有無を表す変数 r と，多数の焼結粒子を区別するための変数 hi(i は粒子番号)を定義する．系内部では体積保存が成り立つことを仮定するため，r は積分値が時間変化しない保存量として扱う．一方で各粒子は焼結中に体積変化するため，hi は非保存量として扱う．粒子の内部では(r,hi)＝(1, 1)，焼結体の外部の気相では(r, hi)＝(0, 0)であり，焼結体の表面では r と hi がそれぞれ 0 から 1 に滑らかに変化する．また本モデルでは粒子は単結晶であることを仮定し，粒子同士の界面を粒界として扱う．粒界では r＝1 かつ複数の hi が 0 から 1 に滑らかに変化する．すなわち，rと hi の分布はそれぞれ焼結体形状と粒界分布を表現する．系の全自由エネルギー F を，r と hi を用いて次式で定義する．F＝fV{fbulk(r, hi)＋k2r2|:r|2＋∑ik2h2|:hi|2} dV ( 1 )ここで，fbulk はバルクの化学的自由エネルギー密度である．右辺被積分項の第 2, 3 項は，それぞれ r と hi についての勾配エネルギー密度であり，kr と kh は勾配エネルギー係数である．熱力学第二法則に基づき，全自由エネルギー F が時間と共に単調減少するように r と hi の時間発展方程式を導く．r の時間発展方程式は，CahnHilliard 型の発展方程式(9)に，固相焼結中に生じる粒子の剛体運動を考慮するための移流項を加えた次式で表す．&r&t＝:･(Mr:dFdr－r∑ini) ( 2 )ここで，ni は i 番目の粒子の剛体運動を表現するための移流速度である．この移流項が加わることで，緻密化が再現される．Mr は焼結体形状変化の速さを決定するモビリティであり，次式で決定される．Mr＝VmRT　Dvolh(r)＋Dvap{1－h(r)}＋Dsurfr2(1－r)2＋Dgb∑i∑j≠ihihj( 3 )ここで，Dvol, Dvap, Dsurf, Dgb がそれぞれ粒子内部・粒子外部(気相)・粒子表面・粒界での原子拡散係数である．Vm はモル体積，R は気体定数，T は絶対温度である．拡散係数はアレニウスの式で決定されることを仮定した．また，hi の時間発展方程式は，AllenCahn 型の発展方程式(10)に移流項を加えた次式で表される．&hi&t＝－MhdFdhi－:･(hini) ( 4 )ここで，Mh は粒界の易動度を決定するモビリティである．移流速度 ni は，以下の式で表されるように，並進速度と回転速度の和として定義する．ni＝mtrViFihi＋mroViTi×[r－rc, i]hi ( 5 )ここで，mtr と mro は粒子の並進運動と回転運動のモビリティである．Fi と Ti は i 番目の粒子に作用する力とトルクである．Vi と rc, i はそれぞれ i 番目の粒子の体積と重心位置を表す．Fi は次式で定義した．Fi＝fVkst∑j≠i(r－r0)〈hihj〉[:hi－:hj]d3r＋Fexi ( 6 )ここで，kst は粒子同士が接触した際に生じる力の大きさと，粒界での空孔過飽和を関連付ける剛性係数である．〈hihj〉は任意の座標が粒界かどうかを判定する演算子である．式( 6 )の右辺第一項は Wang によって提案されたものであり，粒子同士が接触している場合に粒子に働く力をモデル化したものである．一方で Fexi は本研究で独自に導入した図 2 非逐次データ同化手法の概念図．ま て り あMateria Japan第64巻 第10号(2025)力であり，粒子がわずかに離れている場合でも，互いに引き合うような微小な力を生じさせる．2 章で述べたように，実験で用いた粒子はコア粒子の表面を微細なシェル粒子が覆うように分布していたが，シェル粒子については実験データに含まれなかった．そのため，実験データではわずかに離れているように見える場所も，実際にはシェル粒子を介して粒子同士が接触している可能性があった．この可能性を考慮してFexi を導入した．本研究では，Fexi をファンデルワールス力の式の形で定義することで，数 nm 程度離れている場合にのみ微小な力が働く効果を再現した．式( 2 )と( 4 )は，有限差分法を用いて解いた．ここで，有限差分格子の位置を図 1(b)に示した実験データのボクセル座標と一致させることで，実験データと同座標のシミュレーション結果を取得した．比較的物性値が明らかな Cu であっても，400°C程度の低温での粒界拡散係数は未知であり，ni を計算するために必要なパラメータの値も不明である．そこで本研究では，特に焼結中の構造変化に強く寄与する，7 つの材料パラメータ(表面および粒界の拡散係数の前指数因子 D0surf, D0gb および活性化エネルギー Qsurf, Qgb と，Mh, mtr, kst)をデータ同化で推定した．. データ同化手法データ同化手法には，アンサンブルカルマンフィルタや粒子フィルタが代表的な逐次データ同化手法と，4 次元変分法やアジョイント法などの非逐次データ同化手法が存在する(7)．いずれの手法も長所と短所を持つが，焼結は熱処理時間の経過と共に段階が初期(ネック成長)・中期(緻密化)・後期(粒成長)と遷移するため，時系列データ全体の変化の過程を考慮する非逐次データ同化手法が適していると考えた．ベイズ推論に基づくデータ同化では，対象とする系の状態とその時間発展を表現するために，状態ベクトル xt と観測ベクトル yt を定義する．xt は時刻 t(0ttend)におけるシミュレーションのデータをまとめた状態ベクトルであり，本研究のように未知のパラメータを推定する場合には，推定対象のパラメータも xt に含む(このときの xt は拡大状態ベクトルと呼ばれる)．yt は実験データをまとめたベクトルである．本研究の場合は，xt は計算領域内すべての座標の r とhi，および推定対象の 7 種の材料パラメータ D0surf, D0gb, Qsurf,Qgb, Mh, mtr, kst で構成される．一方で yt は，図 1(b)に示したボクセルデータで構成される．xt には数値シミュレーションに由来する誤差が，yt には実験に伴う誤差が必ず含まれる．これらの誤差が独立でガウス分布に従う確率変数であると仮定すると，xt と yt をガウス分布型の確率密度関数(Probability density function: PDF)として扱うことができる．ここで図に，非逐次データ同化手法の概念図を示す．xtと yt が PDF で表される場合，0ttend における xt と yt の「データのずれ」は最尤推定法に基づいて次式で表される．J(x0)＝12(x0－xb0)TB－1(x0－xb0)＋tend∑t＝012(Ht(xt)－yt)TR－1t (Ht(xt)－yt) ( 7 )ここで xb0 は，データ同化を行う前に任意の方法で設定したシミュレーションの初期状態であり，第一推定値とも呼ばれる．B と R はそれぞれ，背景誤差共分散行列と観測誤差共分散行列であり，xb0 および yt に含まれる誤差を表す．Ht(･)は xt から yt と比較可能な量，つまりシミュレーション結果から実験データと比較可能なデータを取り出すための演算子である．本研究の場合は，実験で焼結体の形状データのみが得られているため，xt の要素の中から r のみを取り出す働きをする．式( 7 )は評価関数またはコスト関数と呼ばれ，J(x0)が小さいほど Ht(xt)と yt のずれが小さいことを意味する．したがって，この J(x0)が最小となる xa0 が最適な推定値である．言い換えれば，非逐次データ同化手法はJ(x0)の最小化問題を解くことで最適なシミュレーションの初期状態 xa0 を見つける手法である．一般的な非逐次データ同化手法(4 次元変分法・アジョイント法)は，J(x0)を最小化するために勾配 :J(x0)を導出し，勾配法などを用いて xa0 を探索する．しかしながら，:J(x0)の導出には非常に複雑な計算が必要であり，大きな計算コストと記憶容量を必要とする(11)．そこで筆者らは，勾配を使わず効率的に J(x0)を最小化するために，J(x0)の最小化計算にベイズ最適化(Bayesian optimization: BO)(12)を応用することを着想した．BO はブラックボックス最適化に使われる数理最適化手法であるため，目的関数の勾配を必要としない．また，比較的少ない試行回数で最適解が見つかり，局所解にも陥りにくい手法である．筆者らは BO を用いた非逐次データ同化手法 DMCBO (Data assimilationmethod Minimizing Cost function using BO)(13)を開発し，その有効性を実証した．本研究では，DMCBO のアルゴリズムに TPE(Treestructured Parzen Estimator)(14)を採用することで推定精度を向上させた発展版手法の DMCTPE(15)を使用した．図 3 評価関数の最小化反復過程における (a)(g) 材料パラメータの推定値および (h) 評価関数の変動．(a)(d) 表面拡散係数および粒界拡散係数の前指数因子および活性化エネルギー．(e) 粒界易動度を決定するモビリティ．(f) 並進運動モビリティ．(g) 剛性係数．(h) の実線は反復計算中に得られた J(x0)最小値の推移を示す．図 4 (a) データ同化で推定した材料パラメータを用いた PFシミュレーション結果．(b) 文献値(16)の表面拡散係数を用いた PF シミュレーション結果． 　　　　　　新 進 気 鋭. 固相焼結挙動の高精度予測時刻 t＝0 s から 55 s まで 5 s ごとに得られた実験データのうち，t＝0 のデータを固相焼結 PF シミュレーションの初期状態の作成に使用し，それ以外のデータをデータ同化に使用した．初期状態を実験データから作成したため，初期状態の粒子分布と粒界分布は既知であると仮定した．そのため，J(x0)の最小化のために最適化される初期状態の拡大状態ベクトル x0 の要素のなかで，実際に最適化された要素は推定した 7 種の材料パラメータのみである．初期状態 x0 を用いた PF シミュレーションの実行，J(x0)の評価，J(x0)の最小化計算を 1 サイクルとして，このサイクルを反復することで xa0 を得た．PF シミュレーションの詳細な設定は文献(8)を参照されたいが，空間的・時間的スケールは実験データと同じとした．図に，J(x0)を最小化する反復過程における，7 つの材料パラメータの推定値((a)～(g))，および J(x0)((h))の推移を示す．図 3(a)～(g)の縦軸の範囲が，本研究で設定したTPE による最適パラメータの探索範囲を表している．約150回の反復までの前半部では，比較的広い範囲で最適なパラメータの探索が行われているが，後半になると推定値のばらつきは小さくなり，推定値を連続的に変化させることで 7 つの材料パラメータの最適な組み合わせを探索している．このふるまいの変化は，TPE のアルゴリズムで自動的に行われる．反復回数が300回以上になると J(x0)の最小値の更新がほとんど見られなくなったため，反復は400回目で終了させた．J(x0)が最小となったのは308回目の反復であり，その時に推定された材料パラメータの値は，D0surf＝3.71×10－7 m2/s,Qsurf＝1.46×105 J/mol, D0gb＝2.26×10－6 m2/s, Qgb＝1.63×105 J/mol, mtr＝3.78×10－21 m5/(Js), Mh＝1.28×10－10 m3/(Js), kst＝1.74×106 N/m2 であった．図(a)に，これらの推定されたパラメータを使用した PF シミュレーションの結果を示す．この結果は，データ同化を行うことで PF シミュレーションにより実験データ(図 1)を忠実に再現できたことを示している．推定値がデータ同化対象の系に依存しない値であることを確かめるために，同じ実験条件で別の Cu ナノ粒子クラスターの実験データを取得し，PF シミュレーションと比較した(図)．この視野についてはデータ同化を行わず，シミュレ図 5 別視野での Cu ナノ粒子焼結の (a) 実験データと (b)PF シミュレーション結果．PF シミュレーションは図 4の系でのデータ同化から得られた材料パラメータを使用した．ま て り あMateria Japan第64巻 第10号(2025)ーションには前述の推定値を用いたが，焼結挙動は高精度に予測された．この結果から，本研究で得られた材料パラメータの推定値は，少なくとも同じ実験条件から得られる焼結挙動を Wang の PF モデルで忠実に再現することが可能な値であることが示された．推定した材料パラメータの中で，低温域での Cu の表面拡散係数については文献値が存在する(16)．そこで得られた推定値と文献値を比較すると，推定値は文献値よりも 3 桁程度小さかった．しかしながら，文献値を用いた PF シミュレーションの結果は，実験データと大きく乖離する(図 4(b))．このような差が表れた理由の一つは，実験データやPF モデルの不完全性にあると考えられる．本研究で使用した実験データはシェル粒子が含まれておらず，PF モデルも基礎的なモデルを採用している．今後さらに実験を改良すると同時に，PF モデルを発展させることで，より物理的妥当性が高い材料パラメータが得られるようになると期待できる．また一方で，推定値の妥当性よりも現象の再現や予測を重視するならば，不完全性を有する実験データ・PF モデルであっても，データ同化が有効に働くことも強調しておきたい．. お わ り に本稿では，データ同化を用いて固相焼結の時系列その場観察実験データと PF シミュレーションを融合することで，未知の材料パラメータを推定し，実験観察された焼結挙動を忠実に再現・予測可能なシミュレーションを実現した研究成果を紹介した．本研究で確立したデータ同化のアプローチを基盤とし，今後，実験と計算の両面でさらなる発展を続けることで，焼結体の構造やミクロ組織の高精度予測技術の発展に貢献できると期待される．本研究の成果は，東京農工大学 山中晃徳 教授，山本明保准教授，九州大学 波多聰 教授，斉藤光 准教授，ならびに関連研究室の学生の皆様との共働によって得られたものです．皆様の多大なるご協力とご尽力に，ここに記して謝意を表します．最後に，本研究は JSTCREST「超伝導インフォマティクスに基づく多結晶型超伝導材料・磁石の開発」(JPMJCR18J4)の支援受けて実施されました．関係各位に深く感謝いたします．文 献( 1 ) S. Ihara, M. Yoshinaga, H. Miyazaki, K. Wada, S. Hata, H.Saito and M. Murayama: Nanoscale, 15(2023), 1013310140.( 2 ) S. A. McDonald, C. Holzner, E. M. Lauridsen, P. Reischig, A.P. Merkle and P. J. Withers: Sci. Rep., 7(2017), 5251.( 3 ) 高木知弘，山中晃徳フェーズフィールド法数値シミュレーションによる材料組織設計，養賢堂，(2012), 113.( 4 ) M. Xue and M. Yi: Comput. Model. Eng. Sci., 140(2024),11651204.( 5 ) Y. U. Wang: Acta Mater., 54(2006), 953961.( 6 ) A. Ishii, K. Kondo, A. Yamamoto and A. Yamanaka: Mater.Today Commun., 35(2023), 106061.( 7 ) T. Tsuyuki and T. Miyoshi: J. Meteorol. Soc. Jpn., 85B(2007),331361.( 8 ) A. Ishii, A. Yamanaka, M. Yoshinaga, S. Sato, M. Ikeuchi, H.Saito, S. Hata and A. Yamamoto: Acta. Mater., 278(2024),120251.( 9 ) J. W. Cahn and J. E. Hilliard: J. Chem. Phys., 28(1958), 258267.(10) S. M. Allen and J. W. Cahn: Acta. Metall., 27(1979), 10851095.(11) A. Ishii, A. Yamanaka, E. Miyoshi, Y. Okada and A.Yamamoto: Model. Simul. Mater. Sci. Eng., 29(2021),065012.(12) B. Shahriari, K. Swersky, Z. Wang, R. P. Adams and N. deFreitas: Proc IEEE., 104(2016), 148175.(13) A. Ishii, A. Yamanaka, E. Miyoshi and A. Yamamoto: Mater.Today Commun., 30(2022), 103089.(14) J. Bergstra, R. Bardenet, Y. Bengio and B. K áegl: Adv. Neural.Inf. Process Syst., 24(2011), 25462554.(15) A. Ishii, A. Yamamoto and A. Yamanaka: Sci. Technol. Adv.Mater. Methods, 3(2023), 2239133.(16) H. G. Bowden and R. W. Balluffi: Philos. Mag., 19(1969),10011014.石井秋光★★★★★★★★★★★★★★★★★★★★★2022年 3 月 東京農工大学大学院工学府機械システム工学専攻博士後期課程修了2021年度 JSPS 特別研究員(DC2)2022年 4 月2025年 3 月 国立研究開発法人物質・材料研究機構 ICYS 研究員2025年 4 月現職専門分野計算材料科学，組織予測◎フェーズフィールド法を用いた数値シミュレーションによるミクロ組織予測技術の開発に従事．近年はデータ科学技術の援用も進めている．★★★★★★★★★★★★★★★★★★★★★