2017 Volume 103 Issue 12 Pages 755-762
Heat transfer analyses of casting processes offer an effective way of understanding solidification processes which is important for prevention of formation of casting defects and for control of segregations. The accuracy of the simulations depends on the accuracy of input parameters such as thermal conductivity and heat transfer coefficient. It is not straightforward to evaluate such parameters, especially thermal conductivities of bulks and heat transfer coefficient with high accuracy. It is very important to develop a reliable method for estimation of these parameters. In this study, particle filter, which is a method of data assimilation, is applied to estimation of thermal conductivity of melt and heat transfer coefficient during alloy solidification in mold casting and its applicability is systematically investigated on the basis of twin experiments. It is shown that a constant heat transfer coefficient can be accurately estimated in this method by utilizing a single cooling curve measured in the mold or in the melt near the mold, while the thermal conductivity can be accurately estimated from a single cooling curve measured in the melt. In addition, this method can be utilized for estimation of time-dependent heat transfer coefficient without any assumption and/or approximation on its time dependence. Importantly, the thermal conductivity and time-dependent heat transfer coefficient can be estimated simultaneously with high accuracy.
高品質の合金鋳塊の製造のためには，鋳造プロセスにおける溶質や不純物の偏析を制御し，鋳造欠陥の生成を防止することが必要である。偏析制御や鋳造欠陥の生成防止のためには，凝固過程を正しく理解することが重要であり，数値シミュレーションによる伝熱解析がその有効な手段となる1,2)。伝熱解析では，一般に以下の式を解く。
| (1) |
| (2) |
ここで，ρは密度，Cpは比熱，Tは温度，kは熱伝導率，ΔHは潜熱，fsは固相率，qmold/alloyは溶湯−鋳型間の熱流束，hは溶湯鋳型間の熱伝達係数，Talloyは鋳壁に接している溶湯の温度，Tmoldは溶湯に接している鋳型の温度である。式(1)は鋳型と溶湯における熱伝導方程式，式(2)は鋳型と溶湯の境界における熱伝達式を表す。適切な伝熱解析を行うためには，これらのパラメータに適切な値を入力する必要があるが，その入力値の精度が十分ではない場合があり，そのため伝熱解析が適切に行えないことがある。密度や比熱に比べれば，熱伝導率kと熱伝達係数hに関して精度の問題が生じることが多く，その詳細を以下に説明する。
合金の鋳造過程においては，固液界面で溶質分配が生じて濃度場が不均一になると共に液相中で流動が生じる。濃度場の不均一性と流動は凝固過程に大きな影響を及ぼすため，本来式(1)と式(2)に加えて，溶質拡散方程式と流体の運動方程式を連立して解く必要がある。しかしながら，計算コストの観点から濃度場と流体の解析を行わずに伝熱解析のみによって凝固過程を解析する場合が多々ある。そのようなシミュレーションを行う際には，濃度場の不均一性や流動の影響を含んだ見かけの熱伝導率kが必要となる。そのような見かけの値は，合金系はもちろん，鋳造条件の違いに応じても変化するため，kの見かけの値を高精度かつ簡便に求める手段が必要である。また，濃度場や流体の方程式を連立させたシミュレーションにおいては，kは通常の物性値に相当するが，その場合においても流体の物性であることから，合金によっては実測値が報告されていなかったり，実測値が報告されていても精度が低かったりする場合がある。そして，熱伝達係数hは合金系および鋳造条件に依存するとともに，一般には時間依存のパラメータである3,4,5)。hは温度の実測値と伝熱解析の結果が一致するように試行錯誤によって決定されることが多いが3)，試行錯誤の方法では常に十分な精度が得られるわけではない。さらには，hの時間変化を考慮する場合には，その時間変化を表す近似式を事前に仮定する必要があるが，その近似式が常に正しいとは限らないという懸念も残る。
したがって，kとhを同時かつ高精度に求める手法が必要である。その手法はhの時間変化に事前の仮定を必要とせず，試行錯誤に依らない方法であることが望まれる。また，特殊な鋳造条件のみに限定されずに，種々のプロセスに簡便に適用できること，さらには伝熱解析に濃度場や流体の方程式を連成させた場合にも適用可能な拡張性を有することが望まれる。伝熱問題におけるパラメータを求める試みは，いわゆる伝熱解析の逆問題に相当し，以前から様々な取り組みが行われてきた。例えば，感度解析と収束計算を併用することでhを推定する方法3)，ベイズ推定に基づいて定常状態のkとhを推定する方法6)，特殊な実験系で測定した熱物性から遺伝的アルゴリズムに基づいてkを推定する方法7)など様々な手法が提案されている。これらはそれぞれの適用範囲において優れた手法であると考えられるが，上記の拡張性と簡便さを兼備した手法の発展のためには，さらなる検討の余地がある。特に，kと時間依存のhを同時に推定することは一般に容易ではなく6)，新しい手法の開発が望まれる。本研究ではそのパラメータ推定の方法としてデータ同化8,9,10)に着目した。
データ同化は実測データをシミュレーションに取り入れるデータ科学の手法であり，シミュレーションの精度向上やパラメータの推定などに用いられる。この手法は気象学や海洋学の分野で応用が進んでいる手法である9)。本研究では，データ同化を鋳造プロセスの伝熱解析におけるパラメータ推定に適用し，その有効性を検証した。特に，データ同化の手法として，実装が容易な逐次型の手法である粒子フィルタ(もしくは逐次モンテカルロ法)を用いて，冷却曲線の情報のみから熱伝導率と熱伝達係数を推定することを試みた。伝熱問題はパラメータ推定や状態推定の好個の対象であることから，データ同化を伝熱解析に応用した例は既に報告されている11)。ただし，種々の鋳造プロセスに応用する際の実用的な観点から，(i)kとhの単独推定における温度測定箇所(熱電対の位置)と推定精度の関係，(ii)時間依存のhの推定，(iii)kと時間依存のhの同時推定に関する詳細を明らかにすることが必要であると考えられるが，これらに関して十分な調査は報告されていない。そこで本研究では，これら三項目の内容を詳細に調査し，粒子フィルタの有効性を調査した。
鋳造プロセスにおけるデータ同化の有効性を検証するためには，全てのパラメータが既知の実験系を用意して鋳造中の熱分析を行い，その冷却曲線からkとhの推定を行う必要がある。本研究では，パラメータが既知の実験系を用意する代わりに，双子実験と呼ばれる簡便な検証方法を用いた。この方法では，あらかじめ仮定したkとhの値を用いて伝熱解析を行い，その解析から算出された冷却曲線を実測データとみなしてkとhの推定を行う。そして，その推定結果とあらかじめ仮定したkとhの値を比較することで精度検証を行う。まずは伝熱解析の方法について以下に説明する。
本研究で対象とした鋳型と溶湯の模式図をFig.1に示す。Fig.1(a)は溶湯と鋳型の全体図を表し，Fig.1(b)は実際の計算で対象とした一次元システムを表す。y軸方向の伝熱は無視できるものと仮定し，x軸方向のみの一次元解析を行った。鋳型内と溶湯内の熱伝導は，式(1)を差分法により離散化し，エンタルピー法を用いることで解いた。溶湯の中心部，すなわちFig.1(b)の右端にはミラー条件を用い，鋳型と溶湯の境界温度は式(2)をもとに更新した。また，鋳型は一定温度の大気と接しているものとして，Fig.1(b)の左端の境界温度も熱伝達式により求めた。なお，鋳型は鋳鉄，溶湯はAl-4 mass%Cu合金を想定し，Table 1に示したパラメータを入力値とした。これらの中で，溶湯の熱伝導率kと鋳型−溶湯間の熱伝達係数hがパラメータ推定の対象である。続いて，双子実験における実測データの求め方について説明する。
Schematic illustrations of (a) whole area of mold and alloy and (b) one-dimensional calculation area.
| Parameter | Value | |
|---|---|---|
| Mold (cast iron) | Alloy (Al-4.0 mass%Cu) | |
| Initial temperature (T0/K) | 298.15 | 950 |
| Liquidus temperature (Tl/K) | – | 924.78 |
| Solidus temperature (Ts/K) | – | 905.91 |
| Density (ρ/kgm–3) | 7000 | 2800 |
| Thermal conductivity (k/Wm–1K–1) | 65.4 | Estimation |
| Heat transfer coefficient Mold/Alloy (h/Wm–2K–1) | Estimation | |
| Heat transfer coefficient Air/Mold (h/Wm–2K–1) | 120 | |
| Air temperature (T/K) | 298.15 | |
| Spatial grid spacing (Δx/mm) | 1 | |
| Time step (Δt/s) | 0.1 | |
例として，k=87 Wm−1K−1，h=600 Wm−2K−1としたときの伝熱解析の結果をFig.2に示す。図中の実線，点線，破線は溶湯における冷却曲線であり，それぞれ鋳型から34 mm，16 mm，1 mm離れた位置の温度変化である。本研究ではこれらの冷却曲線から1 s間隔で温度の値を抽出した。実際の鋳造プロセスにおける温度測定には誤差が生じる。その測定誤差を考慮するために，算出された冷却曲線に平均が0 K，標準偏差がσTの正規分布に従うノイズを与え，そのデータを実測値とした。σTの値はTable 2に示す三種類の値を検討した。σT=1 Kのときに得られた実測値の例を図中にプロットで示した。このような実測データを用いて，あらかじめ想定した真の値，つまりk=87 Wm−1K−1，h=600 Wm−2K−1が推定可能であるかを調査することが双子実験の内容である。3章では，hが時間変化する場合の推定について議論するが，その場合には時間依存のhを用いた伝熱解析をあらかじめ行い，Fig.2に示したような実測データを得てから推定を行ったことに注意されたい。
Calculated cooling curves in the alloy at 34 mm, 16 mm and 1 mm from the mold wall. The plots are the calculated values at 1 second interval with addition of noise. These plots correspond to measured values in this study.
| Parameter | Value |
|---|---|
| σT | 1 K, 3 K, 5 K |
| σh | 1% (Estimation of constant parameter) or 10% (Estimation of time-dependent parameter) |
本研究では粒子フィルタを用いてパラメータ推定を行った。粒子フィルタはベイズの定理6)に基づいた状態およびパラメータ推定をモンテカルロ法によって行う手法である。ここでは，kが既知であるとし，時間変化しないhを推定する場合を例に粒子フィルタの手順を示す。
まず，推定の開始において，hの真の値は分からないため，hの値をランダムに複数サンプリングする。Fig.3(a)にそのサンプリングの例を示す。縦軸はh，横軸が時間tの図であり，破線は真の値を表す。この例では，t=0において0~3000 Wm−2K−1の範囲に128通りの値が乱数によりサンプリングされている。そして，粒子フィルタにおいては，これらのhの値を用いた伝熱解析を実測データが得られる時刻まで行う。つまり，この場合は128通りの伝熱シミュレーションを1 sまで実施する。これらの伝熱解析のセットをアンサンブル，それぞれ解析を粒子と呼び，この例では128個の粒子を扱っている。その結果得られた冷却曲線をFig.3(b)に点線で示す。ここでは見やすさのため128本中20本の冷却曲線のみを示した。また，プロットは鋳型近傍の溶湯で測定した温度の実測データであり，計算結果も同じ場所の冷却曲線を表している。t=1 sにおいては，実測データが熱電対の数だけ存在し，そのそれぞれに対応する伝熱解析の結果が128通り存在する。このとき，粒子フィルタではフィルタリングという操作を行って，それぞれのシミュレーション結果の尤度を算出し，その尤度をもとにhの値を更新(リサンプリング)する。ここで，尤度とは実測データに対してシミュレーション結果がどの程度尤もらしいかを表す値である。尤度の計算方法を以下に説明する。
Schematic illustrations of (a) time change of particles (b) time change of cooling curves.
熱電対がm本の場合，t=1 sの時点で温度の実測データがm個ある。これをTobs=(Tobs,1，Tobs,2，…，Tobs,i，…，Tobs,m)Tとベクトル表示する。上付文字のTは転置の操作を意味し，Tobsは縦ベクトルである。同様に，粒子iの計算結果における熱電対の測定位置の温度をTicalc=(Ticalc,1，Ticalc,2，…，Ticalc,i，…，Ticalc,m)Tと表す。本研究で対象とするように，実測データの測定誤差が正規分布に従う場合，t=1 sにおける粒子iの尤度は下記の多変量正規分布に基づいて計算される。
| (3) |
ここで，Σは分散共分散行列(m×m行列)であり，各熱電対の測定誤差に相関がなく，全ての測定点の誤差が同じ分散σTに従う場合は，m×mの単位行列Eを用いて，
| (4) |
と与えられる。この尤度を全てのiに対して求め，下記の式で規格化する。
| (5) |
Nは粒子数であり，ここではN=128である。ここで，計算された粒子iの尤度はシミュレーションに用いたhの尤度に相当する。式(3)から明らかなように，実測データに近い計算結果の尤度は高い。
各粒子の尤度を計算した後は，その尤度に基づいてhのリサンプリングを行う。リサンプリングの目的は，尤度の低い粒子を破棄し，尤度の高い粒子を複製することにある。乱数を用いてリサンプリング後のhの頻度分布が式(5)のβtiに従うように求める。そのリサンプリング後のhの値をFig.3(a)のt=1 sにプロットした。ここで，例えば粒子i=55のh=2000の計算が破棄され，粒子i=123のh=700の計算が複製されるときには，粒子i=123の温度場の計算結果も含めて粒子i=55にコピーすることになる。ただし，粒子i=55と123が完全に同じ計算になることを防ぐため，粒子i=55のhにノイズを与えた。本研究では，そのノイズを平均0，標準偏差がσhの正規分布に従う乱数とした。Table 2に本研究で用いたσhの値を示す。なお，σhの値を変化させた場合に関しても調査を行ったが，本解析の範囲において推定精度はσhに大きく依存しなかったため，本論文ではTable 2に示した値の結果のみを示す。そして，リサンプリング後のhおよび温度場を用いて，伝熱解析を128通り行い，t=2 sまでの冷却曲線を得る(Fig.3(b)参照)。そこで，尤度を計算しリサンプリングを行う。この操作を実測データが存在する時刻で繰り返す。その結果，Fig.3(a)に示すように，hの値全体はこの操作に伴って真の値に近づいていく。
以上の議論では，hの単独推定を例にした。kとhを同時に推定する場合は，リサンプリングをhとkに対して行う点のみが異なり，他の操作は全て同じである。この手法では，伝熱解析のメインとなるアルゴリズムやコードを変更する必要はなく，式(3)−(5)による尤度の計算とリサンプリングの操作を追加するのみで実装可能である。また，伝熱解析に流動解析や濃度場の方程式を連成することも容易である。この点において，この方法は優れた拡張性を有しているといえる。また，実測データとして用いるのは，温度である必要はない。実測可能で計算から算出でき，対象とするパラメータの影響を受けるような量であれば，どのような量でも実測データとして使用できる。この点にもこの手法の魅力があると考えられる。ただし，粒子フィルタでは，粒子数が少ない場合，推定の途中からhの狭い範囲のみに粒子が集中して十分にサンプリングができないという問題(退化)が生じることがある。一般には粒子数の増加と共に推定精度は向上する。したがって，粒子数を十分に多くして推定を行う必要があるが，粒子数の増加は計算コストの増大を意味する。このため，問題に応じた適切な粒子数の設定が必要となる。後述するように，本問題においては，比較的少ない粒子数で高精度な推定が可能であることが示された。
まず，熱伝導率をk=87 Wm−1K−1，熱伝達係数hをh=600 Wm−2K−1の一定値であるとして，hのみの推定を行った。温度測定は一箇所のみ，粒子数は128とし，σT=1 Kの実測データから推定を行った。熱電対を溶湯の中心(鋳壁から34 mmの位置)に置いた場合の結果をFig.4(a)に示す。破線はhの真の値である。ここで，実線は最も尤度の高い粒子のhの値(最尤推定値)を表し，点線はhの最大値と最小値を表している。初期においては最尤推定値のhが大きく振動し，t=10 s付近で真の値に収束している。また，最大値と最小も初期においては幅が広いが，t=15 s付近で真の値の周辺に収束している。続いて，熱電対を鋳壁から1 mm離れた鋳型内に置いた場合の結果をFig.4(b)に示す。Fig.4(a)と比較して，非常に短時間で最尤推定値のh，最大値と最小値が真の値の付近に収束した。このように，熱電対を置く位置によって推定の際の粒子の挙動が著しく異なることが示された。
Time change of heat transfer coefficient estimated from cooling curve (a) at the center of melt and (b) at the mold near the mold wall.
そこで，鋳型や溶湯の様々な位置に熱電対を置いた場合について，それぞれ推定精度の検証を行った。その結果をFig.5に示す。ここでは，t=15 sからt=30 sにおいて1 sごとに得られた最尤推定値の時間平均を推定結果とした。さらに，同条件で5回推定を行い，それらの推定結果の平均値をプロットで表し，最尤推定値の最大値と最小値をエラーバーで示した。また，破線はhの真の値である。横軸は温度測定箇所を鋳壁からの距離で表した値であり，xが正のときは溶湯内部，負のときは鋳型内部で測定したことを意味する。この結果から，鋳型内に熱電対を置いたとき，エラーバーが真の値に対してほぼ10%以内の範囲に収まっており，高精度に推定を行えることが分かった。特に，鋳壁に近いほどより高精度な推定が可能であった。また，溶湯内に熱電対を置く場合であっても，鋳壁から5 mm程度の近い位置であれば比較的高精度に推定可能であることが示された。
Values of heat transfer coefficient estimated by using a cooling curve at different distance from the mold wall. The plot represents the average value of 5 times estimations and the error bars indicate the maximum and minimum values of the estimations.
同様の検証を溶湯の熱伝導率kの推定に対しても行った。真の値はk=87 Wm−1K−1とし，h=600 Wm−2K−1，粒子数を128，σT=1 Kとした。その結果をFig.6に示す。破線はkの真の値である。この図から，鋳型内に熱電対を置いた場合には高精度な推定が行えないこと，鋳壁近傍であっても溶湯内であれば推定が可能であること，および溶湯の中心部に熱電対を置いた場合に最も推定精度が高いことがわかった。ただし，鋳壁から3~5 mm離れた溶湯内の点では，その周囲と比較して推定精度が低下している。この特定の範囲で精度が低下する原因は今のところ不明であり，今後検討する必要がある。
Values of thermal conductivity estimated by using a cooling curve at different distance from the mold wall. The plot represents the average value of 5 times estimations and the error bars indicate the maximum and minimum values of the estimations.
上記の結果は，全てσT=1 Kの実測データを用いた結果である。σTの影響をFig.7と8に示す。Fig.7はh，Fig.8はkの推定結果であり，いずれもx=−3 mm(鋳型内)と3 mm(溶湯内)の結果を示している。また，横軸はσTである。いずれの場合も，σTが大きくなるにつれてエラーバーが長くなっている。ただし，hを推定する場合，熱電対が鋳型にあればσTが大きくても5回の平均値は真の値を良く推定できている。また，kの推定においては，熱電対を溶湯に置く限り，σTに大きく依存せずに高精度な推定が可能であることが示された。次から示す解析においても同様の傾向が示されたため，以降の議論ではσT=1 Kの結果のみを示す。
Dependence of heat transfer coefficient estimated from a cooling curve measured at x = –3 and 3 mm on the standard derivation of error in temperature measurement, σT.
Dependence of thermal conductivity estimated from a cooling curve measured at x = –3 and 3 mm on the standard derivation of error in temperature measurement, σT.
以上のように，データ同化によってhとkのそれぞれを推定可能であることが示された。特に，5回の推定のエラーバーが大きいときには，その平均値も真の値から逸脱するという一般的な傾向が見て取れた。つまり，真の値が不明であっても複数回のパラメータ推定からエラーバーの大きさを評価することで，温度測定箇所と推定精度の関係を予測することが可能である。したがって，実際の鋳造プロセスにおいて熱電対の設置箇所を決める前にこのような解析を実施することで，あらかじめ測定箇所と推定精度の関係を検証しておくことが可能であり，実用上有益であると考えられる。
3・2 時間依存の熱伝達係数の推定前節ではhが時間変化しない場合を想定したが，多くの場合hは時間依存の量である3,12)。通常，hの時間変化を求める際は，あらかじめその近似式を仮定し，その係数等を見積もる方法がとられる。このとき，その近似式がhの時間変化を適切に表現できているのか十分に注意を払う必要がある。ここで，Fig.4(b)の結果に再度注目すると，hを推定する際は鋳壁近傍の冷却曲線を用いることで，非常に短い時間で推定可能であることが分かる。この結果から，時間依存のhに対しても，最尤推定値のhは各時刻で速やかに真の値に収束することが期待される。つまり，データ同化に基づけば，hの時間変化に何ら仮定を導入せずとも推定が可能であると期待される。そこで，本研究では，最尤推定値のhの時間変化が真のhの時間変化を記述することが可能か否かを調査した。なお，各時刻の尤度からhの期待値の時間変化を求めることも可能であるが，本研究ではより簡便だと考えられる最尤推定値を用いた。
鋳造における熱伝達係数は，凝固収縮等の影響で時間とともに減少することが多々あり3,4)，例えば，aとmを定数としてh(t)=at−mの式が用いられている5)。そこで，本研究では，真の値をh(t)=3000 t−0.5 Wm−2K−1とした。このh(t)の式とk=87 Wm−1K−1を用いて再度伝熱解析を行い，そこから得られた冷却曲線を実測データとした。なお，ここではσT=1 Kの結果のみを示す。熱電対の測定箇所は一箇所，粒子数は128とし，σhはTable 2に示す通り元の値の10%とした。Fig.9に推定結果を示す。熱電対を鋳壁から1 mm離れた鋳型中に設置した場合(x=−1 mm)の結果である。破線がhの真の値の時間変化を表し，プロットが各時刻における推定結果である。推定結果は真の値と非常に良い一致を示したことが分かる。つまり，本手法では一つの冷却曲線のみからでも，時間変化するhを高精度に推定可能であることが示された。特に，hの時間変化に事前の仮定を必要としないことは強調すべき重要な事実と考えられる。
Estimated values (plots) of time-dependent heat transfer coefficient h(t) from a single cooling curve measured in the mold at 1 mm away from the mold wall. The dashed line represents the true value of h(t).
続いて，熱電対の位置を変えた場合の推定精度の変化をFig.10に示す。30秒間の冷却について同条件下で5回推定を行い，その推定誤差の平均値をプロットとして表した。エラーバーは推定誤差の最小値と最大値である。推定誤差ΔEは式(6)によって計算した。
| (6) |
Dependence of estimation error for time-dependent heat transfer coefficient on the position of a thermocouple.
ここで，htrue,jとhestimated,jは時刻jにおける真の値と最尤推定値，Mは1回の推定の中で粒子フィルタを適用した回数であり，ここではM=30である。鋳型内に熱電対を置いた場合，鋳型と溶湯の境界から3 mm程度の近い距離のとき推定誤差は10%程度であり，高精度な推定を行うことができることが分かった。また，溶湯内に熱電対を置く場合であっても，鋳壁から1 mm程度の非常に近い位置であれば同程度の精度で推定が可能であった。
3・3 熱伝導率と熱伝達係数の同時推定粒子フィルタによって，熱伝導率kと時間変化する熱伝達係数h(t)の同時推定が可能か否かを検証した。前節までの結果において，kおよびhを単独で推定する場合には，一本の熱電対で十分な精度の推定を行うことができた。ただし，熱電対の位置は，kの推定の際は溶湯，hの推定の際は鋳型もしくは鋳型近傍の溶湯内にすべきであることが示された。そこで，kとh(t)の同時推定においては，鋳型と溶湯の両方に熱電対を設置した。本研究では熱電対の設置を二通り試した。一つ目は，鋳型と溶湯のそれぞれに一本の熱電対を設置した計二本からの推定である。二つ目は，鋳型に二本と溶湯に一本を設置した計三本からの推定である。その熱電対位置の模式図をFig.11に示す。計二本からの推定においては，鋳型中の熱電対位置と溶湯中の熱電対位置は，鋳壁から同じ距離xとした。計三本からの推定においても同様に，鋳型と溶湯に鋳壁から同じ距離の位置に熱電対を設置した。ただし，三本目の鋳型内の熱電対は，もう一方の鋳型内の熱電対から10 mm離れた場所とした。なお，kとh(t)の真の値は，それぞれ，k=87 Wm−1K−1，h(t)=3000 t−0.5 Wm−2K−1とした。粒子フィルタにおいては，kとh(t)の両方をサンプリングする必要があるため，今までよりも多くの粒子を考慮する必要がある。そこで，粒子数は1282とした。
Schematic illustration of positions of thermocouples in simultaneous estimation of the thermal conductivity and time-dependent heat transfer coefficient. The circles (triangles) represent the positions of thermocouples in estimation using two (three) sets of thermocouples.
推定の結果，熱電対の位置を適切に調整することで，kとh(t)を同時かつ高精度に推定可能であることが示された。h(t)に関しては，熱電対が鋳壁近傍にあれば，Fig.10に示したように，各時刻でh(t)の真の値をほぼ予想できていた。その推定誤差をFig.12(a)に示す。推定誤差は式(6)によって計算した。推定は同条件下で5回行い，その推定結果の誤差の平均をプロットで表した。また，エラーバーは推定誤差の最大値と最小値である。熱電対と境界との距離が小さい場合(x≦4 mm)，計二本と計三本の推定ともに，推定誤差10%以内で推定が可能であることが示された。一方で距離が5 mmよりも大きいとき，熱電対が二本と三本の場合でやや差が現れた。xの値がいずれの場合も，熱電対を三本用いた推定の方が精度が高いことが分かる。
Results of simultaneous estimation of the time-dependent heat transfer coefficient and thermal conductivity. (a) Estimation error for the heat transfer coefficient with respect to the distance of thermocouples from the mold wall. (b) Estimated values of thermal conductivity with respect to the distance of the thermocouples from the mold wall.
また，Fig.12(b)にはkの推定結果を示す。これはh(t)と同時に推定した結果であることに注意されたい。破線はkの真の値である。単独で推定する場合と同様に，t=15 sからt=30 sの各時刻において最も尤度の高い粒子の熱伝導率の平均を推定値とした。そして，同条件下で行った5回の推定の推定結果の平均値をプロットで示した。エラーバーは最尤推定値の最大値と最小値である。h(t)の推定誤差の傾向と同様に，鋳壁からの距離が小さい場合は熱電対の本数に関わらず高精度な推定結果が得られた。しかし，熱電対が二本の場合，鋳壁からの距離が3 mmよりも大きいとき急に精度が低下するという結果が得られた。そして，kの推定結果においても，熱電対を三本用いることで推定精度の低下が抑制されることが示された。
本研究ではデータ同化の手法である粒子フィルタを鋳造プロセスにおける熱伝導率と熱伝達係数の推定に応用し，特に，(i)熱電対位置による推定精度の検証，(ii)時間変化する熱伝達係数の推定，そして(iii)熱伝導率と時間依存の熱伝達係数の同時推定に注目し，粒子フィルタの有効性の検証を行った。その結果，熱伝導率を単独で推定する場合は溶湯に，熱伝達係数を単独で推定する場合は鋳型に熱電対を置くことで高精度な推定が可能であることが示された。そして，熱伝達係数の尤度を利用することで，事前に関数形を仮定することなく時間変化する熱伝達係数を推定することができた。特に，鋳型に一本の熱電対を設置するだけで，高精度に推定可能であった。また，熱伝導率と時間依存の熱伝達係数を同時に推定することも可能であった。したがって，粒子フィルタは鋳造プロセスの伝熱解析におけるパラメータの推定に強力な手法として活用できることが示された。
既に述べたように，粒子フィルタの精度と計算コストは，粒子数と共に増加する。本研究では，上記(i)と(ii)の推定には128個という非常に少ない粒子を用いたが，いずれも十分高精度な推定が可能であった。上記(iii)の同時推定においては，1282の粒子数を用いることで高精度な推定が可能であった。いずれの場合においても，今回の伝熱シミュレーション自体が低計算コストであったため，ごく短時間で推定が可能であった。具体的には，Intel® CoreTM i7-4770K(3.50GHz)のCPUと32GBのRAMを搭載した汎用的な計算機でも，(i)と(ii)の推定は約1 s程度，(iii)の推定は約60 s程度であった。近年では，GPUをはじめとする並列計算技術が発展している。粒子フィルタは並列計算に向いた手法であるから，並列化によってより複雑な伝熱シミュレーションにおけるパラメータ推定も可能と考えられる。
