JP4695376B2 - 加熱又は冷却特性評価方法及び装置、反応容器の操業管理方法及び装置、コンピュータプログラム、並びにコンピュータ読み取り可能な記録媒体 - Google Patents

Description

本発明は、例えば鋼材などの各種固体材料を加熱又は冷却する場合（被評価材）に、その被評価材を加熱又は冷却処理するプロセスの加熱又は冷却分布特性を評価するのに用いて好適な加熱又は冷却特性評価方法及び装置、さらには、例えば高炉、燃焼による鋼材加熱炉、石炭ガス化反応炉などの高温のガス反応又は液体反応を伴う反応容器の操業を管理するのに用いて好適な反応容器の操業管理方法及び装置、コンピュータプログラム、並びにコンピュータ読み取り可能な記録媒体に関する。

例えば、鉄鋼製造プロセスにおいては、鋼材を加熱又は冷却するプロセスが多数存在する。熱延プロセスにおけるスラブ加熱炉、焼鈍プロセスにおけるストリップ連続焼鈍炉、過時効炉などは、加熱プロセスの代表プロセスであり、熱延圧延機上でのスプレー水冷却、熱延ホットランテーブルでのラミナー水冷却などは、冷却プロセスの代表プロセスである。これらの各種加熱・冷却プロセスにおいては、加熱バーナーの種類や容量、又は冷却ノズル種類やその冷却特性を好適に選定・配置して、鋼材（被評価材）の加熱又は冷却分布特性を制御する必要がある。

従来から、鋼材（被評価材）の内部に埋め込まれた熱電対により測定された温度情報から、そのプロセスの加熱又は冷却特性を推定することがなされている。通常の場合、鋼材厚み方向に配置された２つの熱電対の温度変化から、その２点の温度差の経時変化を捉え、その値に熱伝導度を乗じて、熱電対間距離で割り戻した値を、平均的な熱流束として評価する方法がある。

しかしながら、前記のような推定では、熱電対周辺の平均的な熱流束を捉えているに過ぎず、鋼材（被評価材）の表面での、空間的な加熱又は冷却分布特性を評価していることにはならない。鋼材（被評価材）内部のある点での平均的な熱流束が鋼材表面での熱流束と等しくならない理由は、実際の現象では、加熱又は冷却の効果は、鋼材（被評価材）表面から、鋼材内部（熱電対埋め込み位置）へ熱伝導によって伝わるが、有限の熱容量を有する熱伝導現象は、表面と内部の間で若干の時間遅れを有する（非定常性）からである。

これに対して、鋼材（被評価材）内の熱伝導現象を非定常１次元の熱伝導逆問題と考えて、１つの熱電対温度変化、又は、１次元方向に並んだ複数の熱電対温度変化から、鋼材表面の熱流束変化を推定する手法が提案されている。図３は、複数の熱電対「×」が埋め込まれた鋼材（被評価材）の２次元断面を示している。鋼材内部に破線で境界を示しているが、１次元とはこの破線に沿った方向の熱流れのみを考慮したことを意味している。すなわち、例えば、１ａ→１ｂ→１ｃや１ｄ→１ｅ→１ｆ方向の熱伝導を想定した場合に、鋼材上面の熱流束を推定する。このとき、鋼材下面の冷却条件を既知と仮定して、未知とした鋼材上面熱流束を求めることが一般的である。もちろん、既知と未知の境界条件を反対にすることも可能である。

また、高炉、燃焼による鋼材加熱炉、石炭ガス化反応炉などの高温のガス反応又は液体反応を伴う反応容器の操業を管理する場合、反応容器内の状況（例えば、燃焼挙動）を観測し、その状況を管理する必要がある。

この場合にも、反応容器の壁に埋め込まれた熱電対により測定された温度情報から反応容器内の状況を推定することがなされている。例えば、急激な温度上昇があれば、その熱電対周辺の反応容器内において異常な発熱が生じていると推定し、逆に極端な温度下降があれば、その熱電対周辺の反応容器内において発熱反応域の縮小などの発熱量低下が生じていると推定するなどの経験的な手法である。

しかしながら、前記のような推定では、実際の反応容器内での温度異常の発生タイミングと温度測定したタイミングとの間でタイムラグが発生することは避けられない。これは、反応容器内の温度異常が熱流束変化として反応容器の表面に伝わり、その後、反応容器の壁材料内部での熱伝導によって熱電対に温度変化をもたらすためであり、原理的に熱伝導現象は若干の時間遅れを有する（非定常性）。

これに対して、反応容器壁内の熱伝導現象を非定常１次元の熱伝導逆問題と考えて、１つの熱電対温度変化、又は、１次元方向に並んだ複数の熱電対温度変化から、反応容器の内表面における熱流束変化を推定する手法が提案されている。例えば図１４は、複数の熱電対「×」が埋め込まれた反応容器（加熱炉）の炉壁近くの２次元断面を示している。炉壁内に破線で境界を示しているが、１次元とはこの破線に沿った方向の熱流れのみを考慮したことを意味している。すなわち、例えば、１ａ→１ｂ→１ｃや１ｄ→１ｅ方向の熱伝導を想定した場合に、炉内表面における熱流束を推定する。このとき、炉外表面の冷却条件を既知と仮定して、未知とした炉内表面における熱流束を求めることが一般的である。もちろん、既知と未知の境界条件を反対にすることも可能である。

ここで、前述したような熱伝導逆問題による推定手法として、例えば、特許文献１では、高炉炉床に埋め込まれた熱電対から、非定常１次元熱伝導方程式の逆問題解析することにより、端点の熱流束を推定する手法について述べられている。この手法の一つは、１点の熱電対温度変化と、端点の冷却条件（既知と仮定）から、その反対側の端点の熱流束を推定する手法である。このような冷却条件は、熱伝達係数と冷却水温度で与えることになるが、特に、熱伝達係数は、冷却水の平均流速から経験相関式により推定することになるので、不確実な推定値になる場合があり、その値を使って逆問題推定した反対側端点の熱流束推定値の精度に、悪影響を及ぼす可能性がある。また、もう一つの手法として、２点の熱電対温度変化を用いた推定手法についても述べているが、２点の内、１点を固定温度境界条件として与えて解く手法であるので、２点の相対的な温度変化を捉えて推定することは難しい上に、固定温度境界条件上での熱流束の推定は可能であるが、固定温度境界条件に選んだ側の、その外側延長線上の端点熱流束は推定できないことになる。更に、前記いずれの方法においても、解析長さを固定して両端の熱流束を求める手法ではなく、耐火物表面に付着する炉内溶融物による厚みの変化と、熱流束変化を同時に推定する手法である。凝固・溶解現象によって付着量を増減するロジックを逆問題解析に導入すると、第一に、計算手続きが複雑になって計算が不安定化しやすくなるという問題がある。第二に、各時間ステップで解析長さを変化させる計算手続きが入ると、長さを変化させた前後の温度分布の推定方法に不確定な要素が混入する可能性があるので、熱流束の推定精度が悪くなる可能性も否定できない。

このように、従来の逆問題解析手法では、不十分な点が多く、複数の熱電対情報から、解析長さを固定して、その両端の熱流束を同時に推定する手法を新たに確立して、非定常な熱流束の変化を精度良く、安定的に推定する技術が重要となる。

これに対して、複数の熱電対の計測温度から、その温度変化を十分に表現できるように、試行錯誤的に温度分布を推定し、両端の温度分布を同時に推定する手法も考えられる。しかし、このような手法では、熱電対の数が増えると計算が複雑化して、全ての熱電対の計測温度変化を満たす温度分布解を得ることは、極めて難しくなる。また、それぞれの熱電対において、計測温度と計算温度の差の絶対値を何処まで小さくすべきかの基準を決めることが困難なので、計算手続きを一般化することが難しい。

この一つの例として、２つの熱電対温度から、特許文献１の手法を応用して、未知の熱流束を、２つの端点で交互に変えて計算し、見かけ上、同時に端点の熱流束を推定する方法が考えられる。即ち、固定温度境界条件とする計測温度を交互に変えて繰り返して計算し、両方の熱電対における計測温度と計算温度が、ある程度一致した時点で、その時間ステップでの両端の熱流束解とするものである。しかし、この手法においては、それぞれの熱電対において、計測温度と計算温度の差の絶対値が、どの程度まで小さくなった時点で解とすべきかを決めることが困難で、場合によっては、片方の熱電対温度を極めてよく表現するが、もう一方の熱電対温度はあまり表現できないような場合でも、解として認識してしまう危険性をはらんでいる。つまり、２つの熱電対位置において、計測温度と計算温度の差の絶対値を最小化するに際して、独立した２つの熱電対位置での最小化のバランスをどの程度にするべきかの基準を、適切に設定することが難しい。さらに、複数の熱電対の場合まで、この方法を拡張すると、解の判定が極めて難しくなることは言うまでもない。

次に、２次元逆問題解析の手法の例としては、例えば、本出願人が特許文献２に開示したものがあり、この手法はそのまま１次元逆問題解析へも適用できる。また、１次元逆問題解析の例として、Beckらにより提案された解析手法が知られている（非特許文献１参照）。

また、逆問題解析の最近の手法として、カルマンフィルター理論や、射影フィルタ理論などの確率的推定法を適用することも考えられる。この手法は、現状では、後述する式（１）の左辺をゼロと置いた、定常熱伝導方程式（観測方程式）への適用が検討されているが、非定常項を含めて適切に観測行列を構成できれば、同様の逆問題解析ができる可能性がある。この定常微分方程式への、確率推定法の適用例としては、非特許文献２に詳しい。

更に、特許文献３には、非定常熱伝導方程式の内挿関数マトリックスについての記述がある。ところが、ここで述べられている内挿関数は、有限要素法などに一般的に用いられている内挿関数の表式を示したものであり、本発明のような非定常熱伝導方程式を満たす内外挿関数とは、全く異なるものである。

特開２００１−２３４２１７号公報 特開２００２−２０６９５８号公報 特開平１０−１００６４号公報 J.V.Beck他著、「Inverse Heat Conduction」, 1985, Wiley, New York, P1−P280 登坂 他著、「逆問題の数理と解法・偏微分方程式の逆解析」、東京大学出版会、１９９９年、Ｐ１９１−Ｐ２８９

ところが、本来の非定常１次元熱伝導逆問題は、鋼材上面及び鋼材下面や炉内表面及び炉外表面での境界条件を同時推定することであり、片側の境界条件を既知と仮定した逆問題解法では、未知とした境界条件の近似的な答えしか得ることができない。例えば鋼材の場合であれば、ある熱電対の温度変動が、上述のような鋼材上面の熱流束変化によるものなのか、鋼材下面の熱流束変化によるものかを区別することはできないことになる。同様に、例えば反応容器の場合であれば、ある熱電対の温度変動が、上述のような反応容器内の熱流束変化によるものなのか、反応容器外に設置された冷却装置の接触不良などによって引き起こされるような反応容器外の熱流束変化によるものかを区別することはできないことになる。

また、より厳密に評価するには、熱伝導現象は、図３においては破線を跨いで左右方向、図１４においては破線を跨いで上下方向にも起こるはずであり、２次元での逆問題を解くことが必要となる。この場合には、図３の左右境界や図１４の上下境界が断熱と仮定した場合においても、上下境界や左右境界の細かな熱流束分布を推定する２次元逆問題を構成する必要があることになる。

本発明は前記のような点に鑑みてなされたものであり、例えば鋼材の上面及び下面や反応容器の内表面及び外表面における熱流束分布や温度分布を同時推定可能とすることを目的とする。更には、１次元だけでなく、２次元、３次元といった空間次元数にも容易に適用可能とすることを目的とする。

本発明は、方法の発明について述べると、被評価材内部の複数の温度情報測定点において測定された温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、前記被評価材の内部又は表面における温度分布又は熱流束分布を演算する加熱又は冷却特性評価方法、及び、反応容器の操業管理方法であって、
前記温度情報測定点は前記被評価材内部に２次元又は３次元的に複数存在し、
前記非定常熱伝導方程式は、密度ρ、比熱Ｃ_p、ｘ方向の熱伝導度ｋ_x、ｙ方向の熱伝導度ｋ_y、ｚ方向の熱伝導度ｋ_zとして、後述する式（１）で表され、
位置ベクトル（ｘ，ｙ，ｚ）、時間ｔとし、上記非定常熱伝導方程式を解析的に満足し厳密解を与えるｘ，ｙ，ｚ，ｔの関数Ｆを用いて、パラメータα_j,i、基準位置ベクトル（ｘ^j，ｙ^j，ｚ^j）、基準時間ｔⁱ、基準位置ベクトルの数Ｎ_j、基準時間の数Ｎ_iとして、Ｆをｊ，ｉに対応してＦ_j,iとし、前記非定常熱伝導方程式の解を、該非定常熱伝導方程式を解析的に満足し厳密解を与える関数として、後述する式（５）により表現するものであって、
ｋを温度情報測定位置、ｌを温度サンプリング時間とし、温度情報測定点において測定された温度情報ａ_k,lとして、前記パラメータα_j,iを、後述する式（６）の連立方程式を用いて決めることを特徴とする。

本発明によれば、鋼材（被評価材）内部の温度情報測定点において測定された温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、前記鋼材（被評価材）の両端面（上面と下面、又は、左面と右面など）における温度分布や熱流束分布を同時推定することができる。従って、例えば、ある温度情報測定点における温度変動が、どちらの端面の熱流束変化によるものなのかを区別するようなことも可能となる。

また、本発明によれば、反応容器の壁内の温度情報測定点において測定された温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、前記反応容器の内表面及び外表面における温度分布や熱流束分布を同時推定することができる。従って、例えば、ある温度情報測定点における温度変動が、反応容器内の熱流束変化によるものなのか、反応容器外の熱流束変化によるものかを区別するようなことも可能となる。

さらに、内外挿関数を用いた簡便な手法により逆問題解析を行うことができ、１次元だけでなく、２次元、３次元といった空間次元数にも容易に適用することができる。

以下、図面を参照して、本発明の実施形態を説明する。
（第１の実施形態）
図１には、本実施形態の加熱又は冷却特性評価装置の概略構成を示す。同図に示すように、加熱又は冷却特性評価装置は、鋼材（被評価材）内部に埋め込まれた熱電対（図３を参照）により測定された温度情報が入力される入力部１０１と、入力部１０１に入力される温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、鋼材（被評価材）表面における温度分布又は熱流束分布を演算する演算部１０２と、演算部１０２により演算された鋼材（被評価材）表面における温度分布又は熱流束分布を、例えば図示しないディスプレイに表示などするための出力部１０３とを備えている。

ここで、非定常熱伝導方程式は、下記の式（１）により表される。

本実施形態では、前記非定常熱伝導方程式を満たす適当な内外挿関数を用いることにより、より簡便な手法としての鋼材（被評価材）表面における境界条件の変化を推定するようにしている。内外挿関数とは、測定点での温度を結んで、その点以外の領域、例えば、解析領域全体又は一部を表現する関数である。外挿のできない内挿関数としては、１次関数近似やスプライン補間などが知られているが、非定常熱伝導方程式を満たしながら、外挿も可能な関数は知られていない。内挿とは既知点に囲まれた内部の未知点を推定することをいい、外挿とは既知点の外側や周囲を含めて推定することをいう。

以下、図２のフローチャートを参照して、前記演算部１０２において行われる演算処理について説明する。演算部１０２では、まず、所定の内外挿関数及びパラメータを用いて非定常熱伝導方程式の解を表現する（ステップＳ２０１）。

本発明者らは、鋭意研究を重ねた結果、下記の式（２）で表現される非定常熱伝導方程式を満たす内外挿関数形を用いることで、より物理的に意味のある内外挿が施せることを見出した。

上式（２）のｔは時間を表し、また、ｘ、ｙ、ｚは位置ベクトル要素を表し、一般の３次元座標系にも適用可能である。τ_x、τ_y、τ_z、Ａ_x、Ａ_y、Ａ_z、Ｘ、Ｙ、Ｚは、適当な任意定数を表し、対象とする系によって、最適な値は変化する。これらの任意定数の値の選択には、注意する必要がある。

別の関数形として、例えば、下式（３）のように表現することも可能である。

同様に、τ_xy、τ_z、Ａ_xy、Ａ_z、Ｘ、Ｙ、Ｚは、適当な任意定数を表し、対象とする系によって、最適な値は変化する。前記式の座標系ｘ，ｙ、ｚの関係は、相互に交換可能であることは言うまでもない。

また、別の関数形として、例えば、下式（４）のように表現することも可能である。

これも同様に、τ_xyz、Ｘ、Ｙ、Ｚは、適当な任意定数を表し、対象とする系によって、最適な値は変化する。

これらの関数Ｆ(ｘ，ｙ，ｚ，ｔ)は、自動的に非定常熱伝導方程式（１）を満たす。この関数Ｆ(ｘ，ｙ，ｚ，ｔ)を用いて非定常熱伝導方程式の解を一般的に表現すると、下記の式（５）として表現される。

上式（５）のｘ^j、ｙ^j、ｚ^jは、任意の基準位置ベクトルの各要素、ｔⁱは任意の基準時間を表し、ｘ、ｙ、ｚ及びｔは、温度を推定しようとしている点の位置ベクトルの要素及び時間である。また、Ｎ_j、Ｎ_iは、それぞれ基準位置ベクトルの数、及び、時間方向の基準時間の数である。これらの数は、それぞれ、温度情報測定点の数、即ち、熱電対による温度測定点の数、及び、測定温度の時間方向のサンプリング数と一致させることが多いが、必ずしも一致させる必要はない。そして、α_j,iはパラメータであるが、この値が決まれば、任意の位置ベクトル（ｘ，ｙ，ｚ）、時間ｔでの温度分布Ｔ(ｘ，ｙ，ｚ，ｔ)を決めることができるのである。

次に、上式（５）により表現される非定常熱伝導方程式の解中のパラメータα_j,iの値を、熱電対により測定された温度情報を用いて決める（ステップＳ２０２）。このパラメータα_j,iの値は、下記の連立方程式（６）を解くことで決めることができる。

上式（６）のａ_k、_lは熱電対により測定された温度Ｔ(ｘ^k，ｙ^k，ｚ^k，ｔ^l)を示しており、上付き文字のｋは測定位置（ｘ^k，ｙ^k，ｚ^k）、上付き文字のｌはサンプリング時間ｔ^lを表す。

以上述べた手法を用いることで、空間及び時間方向に離散的な温度測定データがあれば、非定常熱伝導方程式に支配される鋼材（被評価材）内部領域全体（任意の時空間位置）での温度推定値が得られることになる。

ここで、熱伝導逆問題というのは、計算領域を支配する非定常熱伝導方程式を基にして、領域内部の温度情報を既知として領域境界での温度や熱流束などの境界条件又は初期条件を推定する問題を指す。これに対して、熱伝導順問題というのは、既知である境界条件を基にして、領域内部の温度情報を推定する問題を指す。

前記手法においては、鋼材（被評価材）境界面（表面）での温度分布も同時に推定していることとなり、間接的ではあるが、熱電対により測定された温度情報から鋼材（被評価材）表面の境界条件を決める逆問題となっている。

また、鋼材（被評価材）境界の温度分布だけではなく、その近傍の温度分布から境界での温度勾配が推定できるので、結果的には鋼材（被評価材）境界位置での熱流束変化も推定できることになる。本発明の大きな特徴の一つは、解析領域全体の温度分布の経時変化を簡便に推定できることである。通常の熱伝導逆問題では、熱電対の温度（離散測定点の温度）だけではなく、その他の解析領域における、ある時間断面での温度分布（一般には初期温度分布）が既知であることを前提として、定式化していることが普通である。ところが、実際問題として、鋼材内部の温度分布は、どの時間軸を取っても、不明であることが一般的である。従って、通常の逆問題解法を採用した場合、いろいろな工夫を施して、実際の温度分布を探索・推定しながら、安定的に解を探索する手法を、適宜付加していくことが求められる。ところが、本発明の手法は、原理的には、離散測定点での温度変化さえあれば、解析領域全体の温度分布の経時変化を簡便に推定できるのである。

また、この手法では空間次元数の制約はないので、空間１次元、２次元、３次元の逆問題解析手法としてそのまま適用することができる。

（実施例）
本発明の実施例について説明する。２次元非定常熱伝導方程式系について検討した。２次元の場合、支配方程式は、下記の式（７）のように簡略化される。また、内外挿関数も、同様にして下記の式（８）のように簡略化される。

図４には、ある鋼材内部に埋め込まれた熱電対（４点）を用いて、２次元の非定常熱伝導を仮定して、本発明により、冷却プロセスの逆問題解析を試みたモデルを模式的に示す。鋼材は幅0.1ｍ、高さ0.1ｍの正方形（太線）の形鋼であり、紙面奥行き方向に長く伸びている。４つの熱電対ＴＣ１〜ＴＣ４は、図４に示すように、幅位置（ｘ方向）と高さ位置（ｙ方向）が、それぞれ重ならないように設置している。即ち、正方形左面側の左隅位置を原点（０で図示）として、それぞれの熱電対位置をｘ，ｙ成分で表示すると、ＴＣ１＝（0.075,0.025）、ＴＣ２＝（0.015，0.045）、ＴＣ３＝（0.085，0.055）、ＴＣ４＝（0.025，0.075）である。仮定した鋼材の熱物性値は、比熱Ｃ_p＝434Ｊ／（ｋｇ／Ｋ）、密度ρ＝7500ｋｇ／ｍ³、熱伝導度（等方性）ｋ_x＝30.0Ｗ／（ｍ・Ｋ）、ｋ_y＝30.0Ｗ／（ｍ・Ｋ）である。

図５（ａ）〜５（ｄ）に示す各実線は、図４におけるｑ１、ｑ２、ｑ３、ｑ４の４面に人工的に付与した熱流束境界条件を示している。図５では、熱流束が正（図５（ａ））と負（図５（ｂ）、５（ｃ）、５（ｄ））で示しているが、正は鋼材に熱が流入していることを示し、負は熱が流出していることを示している。それぞれの熱流束は、４面の各全域に均一に付与した。この順問題解析においては、通常の陰解法（時間ステップ6.0秒）を用いた。分割数は、１０×１０分割である。

図６は、初期温度を１７℃均一状態から計算を開始した場合の各熱電対位置の温度変化を示している。本発明による逆問題解析方法によって、これらの熱電対温度を用いて、熱流束ｑ１、ｑ２、ｑ３、ｑ４を推定した。２次元の逆問題の場合も、本発明では正方形領域全体の温度分布を推定しているので（式（５）参照）、境界面の熱流束は、この温度分布を使って推定しなければならない。従って、ｑ２面、ｑ３面に対しては、式（９）を用い、ｑ１面、ｑ４面に対しては、式（１０）を用いて、熱流束を推定している。即ち、推定点の位置ベクトル（ｘ_p又はｙ_p）における温度（Ｔ_p）と、推定点極近傍の位置ベクトル（ｘ_pより3.0ｍｍ内側：ｘ_p,c又は、ｙ_pより3.0ｍｍ内側：ｙ_p,c）における温度（Ｔ_p,c）を式（５）により推定し、式（９）又は式（１０）に代入して、熱流束を計算している。

その結果が、図５のプロットである。順問題で与えた境界条件である実線と一致し、良好に元の境界条件を推定していることが分かる。逆問題解析の時間ステップは12.0秒であり、基準位置ベクトルは、熱電対位置に一致させ（４点）、時間方向の基準時間とその数（３点）は、熱電対の既知とした測定時間とその数に一致させた。即ち、時間方向の既知温度の数は、一つの熱電対に対して３点である（合計１２点）。基準時間を少しずつ前に進めながら経時変化を解析していったが、その最後の点（現在時間）での各端面での熱流束ｑ１、ｑ２、ｑ３、ｑ４を推定している。この時のパラメータの値は、それぞれＸ＝0ｍ、Ｙ＝0ｍ、Ａ_x＝1.0、Ａ_y＝1.0、τ_x＝48000秒、τ_y＝48000秒である。

本実施例のフィッティングでは、４つの各辺の熱流束を平均的に再現できる値を選定しているが、例えば、ｑ１面の精度を優先して、ｑ４の精度の悪化を無視する場合などでは、Ｘ、Ｙ、Ａ_x、Ａ_y、τ_x、τ_y、の設定をｑ１面の精度が向上するように適正化することがある。このように、これらのパラメータの選択は、そのニーズに合わせて変化する。

このようにして求めたパラメータＸ＝0ｍ、Ｙ＝0ｍ、Ａ_x＝1.0、Ａ_y＝1.0、τ_x＝48000秒、τ_y＝48000秒を使って、鋼材冷却時の実温度を用いた逆問題解析を実行した。即ち、約１０００℃に均熱化した鋼材を４つの面より（図４参照）ガスジェット冷却した。図７、図８に、各熱電対の温度経時変化を示している。図７の温度変化を用いて解析したが、４つの熱電対温度がほとんど重なってしまい、その差のイメージが掴めないので、図８に時間を区切った拡大図を示している。鋼材温度が低下するに従って、それぞれの熱電対温度差が縮小していくが、ＴＣ１が一番高く、ＴＣ４が一番低く、ＴＣ２とＴＣ３の温度が中間で近接している傾向はほとんど変化しない。

これらの熱電対温度変化を用いて逆問題解析をした結果を図９に示した。熱流束の推定位置は、４つの各辺の中央位置（図４に○で図示）としている。図９から判断すると、上面ｑ１の冷却能力が最も高く、下面ｑ４の冷却能力が一番低いことが分かる。側面に関しては、ほぼ同等の冷却能力であるが、左面ｑ２の方が、若干冷却能力が高いことが推察できる。このような解析結果を参考にして、冷却ガス量や、ノズル分布を修正して、冷却分布の均一化を図ることが可能となる。２次元逆問題解析では、４面の熱流れを同時に考慮できるので、より精度のよい各面表面での非定常熱流束変化を推定することが可能となる。

ここで、熱電対の設置位置として、例えば、ＴＣ１＝（0.075，0.025）、ＴＣ２＝（0.025，0.025）、ＴＣ３＝（0.075，0.075）、ＴＣ４＝（0.025，0.075）のような人工的な境界条件に対して対称の位置に設定することは望ましくない。この場合、ＴＣ１とＴＣ２、及び、ＴＣ３とＴＣ４の温度変化は、ほぼ等しくなり、ｘ方向の熱流れを推定することが困難となるからである。４つの熱電対ＴＣ１〜ＴＣ４は、図４に示したように、幅位置（ｘ方向）と高さ位置（ｙ方向）が、それぞれ重ならないように設置することが最も望ましいのである。

さらに、我々は、図４と全く同じ問題を、目的関数Ω_K,Lを用いた方法により、最適なパラメータを決めるべく検討した。この場合、２次元問題であるので、目的関数は、パラメータτ_x、τ_y、Ｘ、Ｙの関数として下式（１１）のように表現できる。

式（１１）中の下付文字「given」により、順問題で与えた熱流束を表し、「predict」により、逆問題で推定した熱流束を表している。Ｋは、比較する熱流束の座標位置の数を表し、２次元では、少なくとも４点あることが好ましい。本実施例の場合は、図１０の「○」で示した位置での熱流束を比較したので、Ｋ＝１６点とした。図のように、「○」は、熱流束を付与した各四辺の４等分割中心位置に、それぞれ配置している。また、Ｌは、時間方向の総サンプリング数を表し、少なくとも現在の１点および過去の１点は必要である。本実施例では、１２０秒を起点として９点（６２０秒間隔でサンプリング）としたので、Ｌ＝９点としている。

２次元の場合は、式（１１）に示すように、目的関数Ω_K,Lは、τ_x、τ_y、Ｘ、及びＹの４つのパラメータの関数となる。厳密には、Ａ_x、Ａ_yの関数でもあるが、Ａ_x＝1.0、Ａ_y＝1.0とした。これらの４つのパラメータを全て振って、式（１１）の最小値を求めることとした。変化させた範囲は、τ_x、τ_y、については、12000秒〜2400000秒の間を等間隔で２０分割し、Ｘ、Ｙについては、-3ｍ〜3ｍの間を等間隔で６分割している。即ち、４つのパラメータについて、（２０×２０×６×６）通りの組み合わせについて、式（１１）の値を全て比較して、式（１１）を最小とするパラメータ値を決めた。その結果、τ_x＝228600秒、τ_y＝23400秒、Ｘ＝-3ｍ、Ｙ＝0ｍの時に、目的関数の最小値が得られることが分かった。

全く同様の条件の検討を、熱流束ではなく、温度を用いた下式（１２）で示す目的関数を用いて、比較した。

その結果、式（１２）の最小値を与えるパラメータは、τ_x＝171600秒、τ_y＝12000秒、Ｘ＝1ｍ、Ｙ＝1ｍであり、この実施例の場合には、式（１１）を使った場合とは異なることが分かった。

式（１１）で決定したパラメータを用いて、図５と同様のプロットをした結果が、図１１である。同様に、式（１２）で決定したパラメータを用いて、図５と同様のプロットをした結果が、図１２である。両者ともに、図５の場合よりも、更に良好に順問題で与えた境界条件（実線）を推定していることが分かる。このように、目的関数（１１）、又は、（１２）を用いて、パラメータ決定を行うと、より正確に、順問題で与えた熱流束分布、ひいては解析領域の温度場を推定することができる。目的関数（１１）、又は、（１２）を用いて決定したパラメータにより、図７、図８の熱電対温度変化を用いて逆問題解析を行ったところ、図９とほぼ等しい解析結果が得られた。

以上の考え方は、このような冷却特性の評価のみならず、加熱特性の評価へも容易に適用可能であることは言うまでもない。また、３次元問題に対しても、容易に拡張可能であることも言うまでもない。

（第２の実施形態）
図１３には、第２の実施形態の反応容器の操業管理装置の概略構成を示す。同図に示すように、反応容器の操業管理装置は、反応容器の壁内に埋め込まれた熱電対（図１４を参照）により測定された温度情報が入力される入力部１０１と、入力部１０１に入力される温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、反応容器の内表面及び外表面における温度分布又は熱流束分布を演算する演算部１０２と、演算部１０２により演算された反応容器の内表面及び外表面における温度分布又は熱流束分布を、例えば図示しないディスプレイに表示などするための出力部１０３とを備えている。

ここで、非定常熱伝導方程式は、上式（１）により表される。本実施形態では、前記非定常熱伝導方程式を満たす適当な内外挿関数を用いることにより、より簡便な手法としての反応容器の内表面及び外表面における境界条件の変化を推定するようにしている。内外挿関数とは、測定点での温度を結んで、その点以外の領域、例えば、解析領域全体又は一部を表現する関数である。外挿のできない内挿関数としては、１次関数近似やスプライン補間などが知られているが、非定常熱伝導方程式を満たしながら、外挿も可能な関数は知られていない。内挿とは既知点に囲まれた内部の未知点を推定することをいい、外挿とは既知点の外側や周囲を含めて推定することをいう。

なお、前記演算部１０２において行われる演算処理については、上記第１の実施形態で説明したものと同様であり（図２及び式（２）〜式（６））、ここではその詳細な説明は省略する。

以上述べた手法を用いることで、空間及び時間方向に離散的な温度測定データがあれば、非定常熱伝導方程式に支配される反応容器壁領域全体（任意の時空間位置）での温度推定値が得られることになる。

ここで、熱伝導逆問題というのは、計算領域を支配する非定常熱伝導方程式を基にして、領域内部の温度情報を既知として領域境界での温度や熱流束などの境界条件又は初期条件を推定する問題を指す。これに対して、熱伝導順問題というのは、既知である境界条件を基にして、領域内部の温度情報を推定する問題を指す。

前記手法においては、反応容器の壁境界での温度分布も同時に推定していることとなり、間接的ではあるが、熱電対により測定された温度情報から反応容器の内表面及び外表面の境界条件を決める逆問題となっている。

また、反応容器の壁境界の温度分布だけではなく、その近傍の温度分布から境界での温度勾配が推定できるので、結果的には反応容器の壁境界位置での熱流束変化も推定できることになる。

本発明の大きな特徴の一つは、反応容器の内表面及び外表面だけでなく解析領域全体の温度分布の経時変化を簡便に推定できることである。通常の熱伝導逆問題では、熱電対の温度（離散測定点の温度）だけではなく、その他の解析領域における、ある時間断面での温度分布（一般には初期温度分布）が既知であることを前提として、定式化していることが普通である。ところが、実際問題として、反応容器壁内の温度分布は、どの時間軸を取っても、不明であることが一般的である。従って、通常の逆問題解法を採用した場合、いろいろな工夫を施して、実際の温度分布を探索・推定しながら、安定的に解を探索する手法を、適宜付加していくことが求められる。ところが、本発明の手法は、原理的には、離散測定点での温度変化さえあれば、解析領域全体の温度分布の経時変化を簡便に推定できるのである。

また、この手法では空間次元数の制約はないので、空間１次元、２次元、３次元の逆問題解析手法としてそのまま適用することができる。

（実施例１）
本発明の実施例について説明する。１次元非定常熱伝導方程式系について検討した。１次元の場合、支配方程式は式（１３）のように簡略化される。また、内外挿関数も、同様にして式（１４）のように簡略化される（Ａ_x＝1.0）。

図１５、図１６には、高炉炉底の底盤（図１５）、側壁（図１６）に設置されているカーボン煉瓦内に埋め込まれた熱電対を用いて、１次元の非定常熱伝導を仮定して、本発明により逆問題解析を試みたモデルを模式的に示す。いずれもの場合においても、熱電対は、冷却面側に偏って２本埋め込まれており、ＴＣ１は高温側熱電対、ＴＣ２は低温側熱電対を示す。これら熱電対ＴＣ１、ＴＣ２で測定される温度変化から、本発明に示した逆問題解析により、高温熱流束面での非定常熱流束ｑ１と、冷却面での非定常熱流束ｑ２とを同時に推定する。図１５では、高温熱流束面位置は冷却面から4.0ｍ奥の定点（ｘ＝4.0ｍ）、図１６では、高温熱流束面位置は冷却面から2.0ｍ奥の定点（ｘ＝2.0ｍ）である。仮定したカーボン煉瓦の熱物性値は、比熱Ｃ_p＝712Ｊ／（ｋｇ／Ｋ）、密度ρ＝2300ｋｇ／ｍ³、熱伝導度ｋ_x＝21.2Ｗ／（ｍ・Ｋ）である。

本発明は、定数τ_x、Ｘを適切に設定することが重要である。そこで、まずは、図１７Ａ又は図１８Ａのような人工的な境界条件ｑ１、ｑ２を付与した順問題解析を実行し、ＴＣ１、及び、ＴＣ２温度の経時変化を求めた（図１７Ｂ又は図１８Ｂ）。ここで、順問題解析では、通常の差分近似の陰解法（時間ステップ28800秒）を用いた。そして、図１７Ｂ又は図１８Ｂの各熱電対温度を既知条件として、本発明による逆問題解析を実行し、最適なτ_x、Ｘを探索する。図１５のモデルに対しては、図１７の条件を付与し、図１６のモデルに対しては、図１８の条件を付与した。本実施例では、熱電対間隔は不変であると考えることができるので、Ｘ＝0ｍ一定が妥当と判断して、好適なτ_xを選定することにした。τ_xを選定するに際しては、基準位置ベクトルや、基準時間をどこに設定するかも重要な要素である。図１９、図２０に、基準点（基準位置ベクトルと基準時間）、温度既知点（温度測定している熱電対位置ベクトルと温度既知時間）、推定点（温度推定する位置ベクトルと推定時間）の関係を模式的に示している。図１９、図２０の例では、基準点と温度既知点を一致させて、図１９では、時間方向５点、図２０では、時間方向３点の温度既知点を用いていることを示している。図１９、図２０のいずれの場合も、推定点は、時間方向には、最も現在に近い点（「現在」と図示）を選び、位置ベクトルとしては、ｑ１、及び、ｑ２位置に設定した。図１９及び図２０中では、推定点を「×」で示している。図１５のモデルに対しては、図１９の設定で推定点での熱流束ｑ１、ｑ２を求め、図１６のモデルに対しては、図２０の設定で推定点での熱流束ｑ１、ｑ２を求め、それぞれ、実際に与えた境界条件である図１７、図１８のｑ１、ｑ２と比較して、比較的良好に境界条件を再現できるτ_xの値を決定することにした。パラメータα_j,iを決めるに際しては、各時間ステップにおいて、式（６）の連立方程式を解いて更新・決定している。本発明での逆問題解析では、式（５）で示したように、領域全体の温度分布を推定することになるので、ｑ１、ｑ２の値は、直接には求めることができない。従って、推定点の位置ベクトル（ｘ_p）における温度（Ｔ_p）と、推定点極近傍の位置ベクトル（ｘ_pより3.0ｍｍ内側：ｘ_p,c）における温度（Ｔ_p,c）を推定し、上式（９）で計算して、推定点での熱流束qであると仮定した。

まずは、図１５のモデルに対して、順問題解析で与えたｑ１及びｑ２と、逆問題解析で求めたｑ１及びｑ２を比較した結果を図２１に示す。実線が順問題で与えたｑ１（太線）、ｑ２（細線）であり、プロットが逆解析で求めたｑ１（◆）、ｑ２（■）の推定値である。逆解析の時間ステップは28800秒であり、このときのτ_xの値は1800000秒であった。このτ_xの値は、解析の単位系によっても変わるが、今回の解析ではMKS単位系を用いている。ｑ１、ｑ２、共に、比較的良好に順問題の境界条件を再現していることが分かる。τ_xの選定値には比較的余裕幅があり、τ_x＝1800000秒±200000秒程度では、熱流束推定値に大きな精度悪化は見られなかった。ただし、τ_x＝1000000秒、τ_x＝2400000秒などに設定すると、明らかに推定精度が低下することが確認されている。

同様にして、図１６のモデルに対して、順問題解析で与えたｑ１及びｑ２と、逆問題解析で求めたｑ１及びｑ２を比較した結果を図２２に示す。実線が順問題で与えたｑ１（太線）、ｑ２（細線）であり、プロットが逆解析で求めたｑ１（◆）、ｑ２（■）の推定値である。逆解析の時間ステップは28800秒であり、このときのτ_xの値は1000000秒であった。ｑ１、ｑ２、共に、ほぼ一致し、順問題の境界条件を再現していることが分かる。τ_x＝1800000秒±200000秒にすると、推定精度が極端に悪化する。即ち、図１６のモデルに対しては、図１５の場合とは異なったτ_xの値で、推定精度が良好となることが分かった。このように、好適なτ_xは、基準点、温度既知点、推定点の相対的な位置関係や、それらの数によって、異なった値になる傾向にあることが分かる。

以上のようにして求めたτ_xの値を使って、高炉の実温度を用いて逆問題解析した結果（４月１日から約４ヶ月間の熱流束ｑ１及びｑ２経時変化）が、図２３及び図２４である。図２３は、高炉炉底底盤の図１５のモデルの場合（τ_x＝1800000秒）であり、図２４は、高炉炉底側壁の図１６のモデルの場合（τ_x＝1000000秒）である。これらの図は、もっと長期に渡った解析結果から抜粋したものである。

いずれの図においても、図２３Ｃ、図２４Ｃのグラフの熱電対ＴＣ１温度、図２３Ｄ、図２４Ｄのグラフの熱電対ＴＣ２温度データを用いて、図２３Ａ、図２４Ａのグラフの高温側熱流束ｑ１と、図２３Ｂ、図２４Ｂグラフの低温側熱流束ｑ２を同時推定している。図２３Ａ、図２４Ａのグラフ中の「定常法」は、従来より用いられている貫流熱流束推定方法であり、ＴＣ１とＴＣ２の温度差の絶対値に、熱伝導度ｋ_xを乗じて、ＴＣ１位置とＴＣ２位置の距離の絶対値で割り戻した値である。特に図２４に顕著に観察されるが、本発明のｑ１推定より数日遅れて、定常法のｑ１が同じ様な形状で推移をしていることが分かる。これは、高炉カーボン煉瓦の熱容量（ρＣ_p）の効果による遅れ時間の影響である。定常法では、この遅れ時間が考慮できないので、平均的な熱流束しか評価することができないが、本発明の方法では、この効果も考慮した推定ができるので、炉内の熱流束変化を早く検知することが可能である。

次に、我々は、好適なτ_xを決めるべく検討を重ね、解析領域の両端点（２点）における熱流束ｑ１及びｑ２について、下式（１５）で表される目的関数Ω_K,Lの極小点、又は、最小点になるように、τ_x、Ｘを設定することが望ましいことを見出した。

式（１５）中の下付文字「given」により、順問題で与えた熱流束を表し、「predict」により、逆問題で推定した熱流束を表している。Ｋは、比較する熱流束の座標位置の数を表し、１次元では最低２点は必要である。本実施例の場合は、両端点（２点）の熱流束を比較している。比較する座標位置の数Ｋは、本実施例のように、２点である必要はなく、解析領域中央の適当な位置での熱流束も加えて、３点以上にしても問題ないことは言うまでもない。また、Ｌは、時間方向の総サンプリング数を表し、少なくとも現在の１点および過去の１点は必要である。１次元の場合は、目的関数Ω_K,Lは、τ_x及び、Ｘの関数となる。本実施例の場合をより具体的に書き下すと、下式（１６）のようになる。

ここで、これまでの表記ｑ１は、ｑ_given,1,l、又は、ｑ_predict,1,lに対応し、ｑ２は、ｑ_given,2,l、又は、ｑ_predict,2,lに対応している。

図１５のモデルを用い、Ｘ＝0ｍ一定として、前記目的関数の特性をτ_xの関数として調査した結果を図２５に示した。ここで、目的関数Ω_K,Lを求めるに際して、時間方向の総サンプリング数Ｌは、80hrを起点として２８点（８８時間間隔でサンプリング）で、図１８Ａのような１８０時間周期の台形状の熱流束ｑ１と、周期状のｑ２について、式（１６）を適用した。図２５の縦軸は、目的関数Ω_K,Lであるが、値の変化が急激であるので、対数目盛で示している。複数のラインが存在するのは、順問題で与えた熱流束として、図１８Ａの他に、図１８Ａの台形状の熱流束の最大値や最小値、及び、その周期を、適当に変化させて式（１６）を適用してτ_xの変化を調査した結果もプロットしているためである。

図２１を用いて視覚的な判断で決めたτ_xの値（Ｘ＝0ｍ一定）は、τ_x＝1800000秒であったが、図２５では、τ_x＝1250000秒付近にも極小値が存在することが分かる。その他にも極小値を示す値が存在し、極小値のτ_xを使用して実温度データを用いて解析することにより、極小値以外でのτ_xにより解析した場合よりも好適となるのであるが、例えば、最初の極小値であるτ_x＝700000秒を用いて、実温度データを用いて解析してみると、感度が高くなりすぎて、推定値の振動が激しくなることが判明した。また、２番目の極小値（最小値）であるτ_x＝940000秒は、極小値の幅が極めて狭いので、条件の変化によって、精度が悪化することが予想される。従って、本実施例では、τ_x＝1250000秒が、好適τ_xであると判断した。このように、いくつかの極小値が存在している場合に、どのτ_xの値を採用するかについては、別途、実温度データを用いた解析評価や、その極小形状による判断を必要とする場合もあるが、目的関数Ω_K,Lを使って、どのτ_xの値が採用候補になり得るかの当たりを付けることにより、より好適なτ_xの値を簡単に探索することが可能となるのである。また、この例では、Ｘ＝0ｍで固定して好適値を探索したが、τ_x及びＸの値を同時に変化させて極小値を探索してもよいことは言うまでもない。

図２６Ａには、τ_xの好適値を決める際に目的関数Ω_K,Lを用いた場合（τ_x＝1250000秒）と用いない場合（τ_x＝1800000秒）の高温側熱流束ｑ１の推定結果の違いを示している。この時解析に用いた２つの熱電対温度変化は、図２６Ｂ、図２６Ｃに示している。目的関数Ω_K,Lを用いない場合の結果は、図２３の中から、４月２９日より約１０日間を抜き出してプロットしている。両者に大きな違いは観察されないが、目的関数Ω_K,Lを用いた場合（τ_x＝1250000秒）の方が、若干早めに熱流束の動きを推定していることが分かる。

（実施例２）
図２５のモデルを用いて、「未来の温度変化」を推定する例についても説明する。図２７に、基準点、温度既知点、推定点の関係を模式的に示した。本実施例の場合では、推定点の位置をＴＣ１の位置、及び、ＴＣ２の位置として、時間方向には、未来の温度変化を推定することを意味する。温度既知点（過去）、及び、推定点（未来）の単位時間ステップは、28800秒として、図２７のように、異なった位置ではあるが、等間隔に並んでいる。また、基準点は、半時間ステップ14400秒ずらして設定し、基準点の時間ステップは、２倍の57600秒としている。また、本実施例では、上式（６）の内外挿関数を用いて、τ_x＝6000000秒、Ｘ＝0ｍとして、解析を実施している。その他の設定は、実施例１と同様である。

図２８には、解析結果の一例を示す。本発明方法の特徴を示す為に、人工的な境界条件であるｑ１及びｑ２（太実線右軸）を与えた時の、ＴＣ１位置及びＴＣ２位置での温度変化（細実線左軸）を用いて、その推定精度を検証した。この順問題解析は、通常の差分近似の陰解法（時間ステップ28800秒）を用いた。図２８のプロット（◆、■）は、本発明の方法による推定結果を示すが、そのプロットの直前にある６つの既知温度点（ＴＣ１が３つ、ＴＣ２が３つ）を用いて、６つの推定点（未来点）を推定した結果である（図２７参照）。全体的には、ほぼ一致し、非常によい推定性能を有していることが分かる。ただし、一部の点において、全く違った方向の推定値を示している場合が散見できる。その代表例を楕円で囲って示しているが、本発明の推定値が「低下」しているのに対して、実際には、温度上昇となっている。これは、既知温度点での温度が、ステップ状の大きな境界条件の変化を跨いだ点に位置しており、本発明で用いた既知温度点が、熱流束が低い時（1000Ｗ/ｍ²）の影響を多大に受けているためであると考えられる。このように、極端に大きく境界条件が変化する場合は、本発明方法の推定精度は低下するが、境界条件に大きな変化がない場合には、極めて良好な推定精度を有することが分かる。

全く同じような検討結果として、図２９に、台形状の熱流束によって、ＴＣ１、ＴＣ２を生成させて、本発明により未来点を推定した場合の結果を示す。この場合には、図２８に比べると、滑らかな境界条件の変化であるため、本発明による未来温度の推定精度は、極めて良好となることが分かる。

（実施例３）
別の本発明の実施例について説明する。２次元非定常熱伝導方程式系について検討した。２次元の場合、支配方程式は上式（７）のように簡略化される。また、内外挿関数も、同様にして上式（８）のように簡略化される。

図３０には、ある鋼鉄製反応容器壁面内部に埋め込まれた熱電対（４点）を用いて、２次元の非定常熱伝導を仮定して、本発明により逆問題解析を試みたモデルを模式的に示す。厚み0.1ｍの鋼鉄壁は、図の上下方向に続いているが、本モデルでは、高さ0.2ｍ、厚み0.1ｍの長方形領域（太線破線）を解析領域として抜き出して、その４辺を貫流する熱流束を推定するモデルを構成した。４つの熱電対ＴＣ１〜ＴＣ４は、図３０に示したように、厚み位置（ｘ方向）と高さ位置（ｙ方向）が、それぞれ重ならないように設置している。即ち、長方形領域の左隅位置を原点（０で図示）として、それぞれの熱電対位置をｘ，ｙ成分で表示すると、ＴＣ１＝（0.075，0.05）、ＴＣ２＝（0.015，0.09）、ＴＣ３＝（0.085，0.11）、ＴＣ４＝（0.025，0.15）である。仮定した鋼鉄壁の熱物性値は、比熱Ｃ_p＝434Ｊ／（ｋｇ／Ｋ）、密度ρ＝7500ｋｇ／ｍ³、熱伝導度（等方性）ｋ_x＝30.0Ｗ／（ｍ・Ｋ）、ｋ_y＝30.0Ｗ／（ｍ・Ｋ）である。

図３１（ａ）〜（ｄ）に示す実線は、図３０におけるｑ１、ｑ２、ｑ３、ｑ４の４面に人工的に付与した熱流束境界条件を示している。図３１では、熱流束が正（図３１（ａ））と負（図３１（ｂ）〜（ｄ））で示しているが、正は破線長方形領域に熱が流入していることを示し、負は熱が流出していることを示している。それぞれの熱流束は、４面の各全域に均一に付与した。この順問題解析においては、通常の陰解法（時間ステップ6.0秒）を用いた。分割数は、１０×１０分割である。

図３２は、初期温度を１７℃均一状態から計算を開始した場合の各熱電対位置の温度変化を示している。本発明による逆問題解析方法によって、これらの熱電対温度を用いて、熱流束ｑ１、ｑ２、ｑ３、ｑ４を推定した。２次元の逆問題の場合も、本発明では長方形領域全体の温度分布を推定しているので（上式（５）参照）、境界面の熱流束は、この温度分布を使って推定しなければならない。従って、ｑ２面、ｑ３面に対しては、上式（９）を用い、ｑ１面、ｑ４面に対しては、上式（１０）を用いて、熱流束を推定している。即ち、推定点の位置ベクトル（ｘ_p又はｙ_p）における温度（Ｔ_p）と、推定点極近傍の位置ベクトル（ｘ_pより3.0ｍｍ内側：ｘ_p,c又は、ｙ_pより3.0ｍｍ内側：ｙ_p,c）における温度（Ｔ_p,c）を上式（５）により推定し、上式（９）又は上式（１０）に代入して、熱流束を計算している。

その結果が、図３１のプロットである。順問題で与えた境界条件である実線と一致し、良好に元の境界条件を推定していることが分かる。逆問題解析の時間ステップは12.0秒であり、基準位置ベクトルは、熱電対位置に一致させた。また、時間方向の既知温度の数は、一つの熱電対に対して３点であり（合計１２点）、基準時間を少しずつ前に進めながら経時変化を解析していったが、その最後の点（現在時間）での各端面での熱流束ｑ１、ｑ２、ｑ３、ｑ４を推定した。時間方向の基準点と温度既知点は一致させている。この時のパラメータの値は、それぞれＸ＝0ｍ、Ｙ＝0ｍ、Ａ_x＝1.0、Ａ_y＝1.0、τ_x＝55000秒、τ_y＝55000秒である。本実施例のフィッティングでは、４つの各辺の熱流束を平均的に再現できる値を選定しているが、例えば、ｑ１面の精度を優先して、ｑ４の精度の悪化を無視する場合には、Ｘ＝0ｍ、Ｙ＝0ｍ、Ａ_x＝1.0、Ａ_y＝1.0、τ_x＝72000秒、τ_y＝72000秒などに設定する方法も考えられる。このように、これらのパラメータの選択は、そのニーズに合わせて変化する。

このようにして求めたパラメータＸ＝0ｍ、Ｙ＝0ｍ、Ａ_x＝1.0、Ａ_y＝1.0、τ_x＝55000秒、τ_y＝55000秒を使って、反応炉鋼鉄壁の実温度を用いた逆問題解析を実行した。図３３に、各熱電対の温度経時変化を示している。これらの熱電対温度変化を用いて逆問題解析をした結果を図３４に示した。熱流束の推定位置は、４つの各辺の中央位置（図３０に○で図示）としている。２次元逆問題解析では、反応容器面から冷却面への熱流れだけでなく、上下方向の熱流れも考慮できるので、反応炉内部のより精度のよい非定常熱流束変化を推定することが可能である。図３４の例では、３時間３０分頃より、反応容器内部の熱流束ｑ２が低下しているが、それに合わせて冷却面側の冷却能力も徐々に緩めたこと（ｑ３の絶対値低下）が、はっきりと観察できる。

ここで、熱電対の設置位置として、例えば、ＴＣ１＝（0.080，0.05）、ＴＣ２＝（0.020，0.05）、ＴＣ３＝（0.080，0.15）、ＴＣ４＝（0.020，0.15）のような人工的な境界条件（図３０参照）に対して対称の位置に設定することは望ましくない。この場合、ＴＣ１とＴＣ２、及び、ＴＣ３とＴＣ４の温度変化は、ほぼ等しくなり、ｘ方向の熱流れを推定することが困難となるからである。４つの熱電対ＴＣ１〜ＴＣ４は、図３０に示したように、厚み位置（ｘ方向）と高さ位置（ｙ方向）が、それぞれ重ならないように設置することが最も望ましいのである。

さらに、我々は、図３０と全く同じ問題を、目的関数Ω_K,Lを用いた方法により、最適なパラメータを決めるべく検討した。この場合、２次元問題であるので、目的関数は、パラメータτ_x、τ_y、Ｘ、Ｙの関数として上式（１１）のように表現できる。

上式（１１）中の下付文字「given」により、順問題で与えた熱流束を表し、「predict」により、逆問題で推定した熱流束を表している。Ｋは、比較する熱流束の座標位置の数を表し、本実施例の場合は、図３５の「○」で示した位置での熱流束を比較したので、Ｋ＝１６点とした。図のように、「○」は、熱流束を付与した各四辺の４等分割中心位置に、それぞれ配置している。また、Ｌは、時間方向の総サンプリング数を表し、本実施例では、１２０秒を起点として９点（６２０秒間隔でサンプリング）としたので、Ｌ＝９点としている。

２次元の場合は、上式（１１）に示すように、目的関数Ω_K,Lは、τ_x、τ_y、Ｘ、及びＹの４つのパラメータの関数となる。厳密には、Ａ_x、Ａ_yの関数でもあるが、Ａ_x＝1.0、Ａ_y＝1.0とした。これらの４つのパラメータを全て振って、上式（１１）の最小値を求めることとした。変化させた範囲は、τ_x、τ_y、については、12000秒〜2400000秒の間を等間隔で２０分割し、Ｘ、Ｙについては、-3ｍ〜3ｍの間を等間隔で６分割している。即ち、４つのパラメータについて、（２０×２０×６×６）通りの組み合わせについて、上式（１１）の値を全て比較して、上式（１１）を最小とするパラメータ値を決めた。その結果、τ_x＝171600秒、τ_y＝46200秒、Ｘ＝-2ｍ、Ｙ＝-2ｍの時に、目的関数の最小値が得られることが分かった。

全く同様の条件の検討を、熱流束ではなく、温度を用いた上式（１２）で示す目的関数を用いて、比較した。

その結果、上式（１２）の最小値を与えるパラメータは、τ_x＝171600秒、τ_y＝46200秒、Ｘ＝-2ｍ、Ｙ＝-2ｍであり、この実施例の場合には、上式（１１）を使った場合と一致することが分かった。

これらのパラメータを用いて、図３１と同様のプロットをした結果が、図３６である。図３１の場合よりも、更に良好に順問題で与えた境界条件（実線）を推定していることが分かる。このように、目的関数（１１）、又は、（１２）を用いて、パラメータ決定を行うと、より正確に、順問題で与えた熱流束分布、ひいては解析領域の温度場を推定することができる。目的関数（１１）、又は、（１２）を用いて決定したパラメータにより、図３３の熱電対温度変化を用いて逆問題解析を行ったところ、図３４とほぼ等しい解析結果が得られた。以上の考え方は、３次元問題に対しても、容易に拡張可能であることは言うまでもない。

上述した実施形態の加熱又は冷却特性評価装置や反応容器の操業管理装置は、コンピュータのＣＰＵ或いはＭＰＵ、ＲＡＭ、ＲＯＭなどにより構成され、ＲＡＭやＲＯＭに記憶されたプログラムが動作することによって実現される。従って、コンピュータに対し、前記実施形態の機能を実現するためのプログラム自体が上述した実施形態の機能を実現することになり、そのプログラム自体は本発明を構成する。

また、前記プログラムをコンピュータに供給するための手段、例えばかかるプログラムを格納した記録媒体は本発明を構成する。かかるプログラムコードを記録する記録媒体としては、例えばフレキシブルディスク、ハードディスク、光ディスク、光磁気ディスク、ＣＤ−ＲＯＭ、磁気テープ、不揮発性のメモリカード、ＲＯＭなどを用いることができる。

また、コンピュータが供給されたプログラムを実行することにより、上述の実施形態の機能が実現されるだけでなく、そのプログラムがコンピュータにおいて稼働しているＯＳ（オペレーティングシステム）或いは他のアプリケーションソフトなどと共同して上述の実施形態の機能が実現される場合にもかかるプログラムは本発明の実施形態に含まれることはいうまでもない。

更に、供給されたプログラムがコンピュータの機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに格納された後、そのプログラムの指示に基づいてその機能拡張ボードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部又は全部を行い、その処理によって上述した実施形態の機能が実現される場合にも本発明に含まれることはいうまでもない。

なお、前記実施形態において示した各部の形状及び構造は、何れも本発明を実施するにあたっての具体化のほんの一例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその精神、又はその主要な特徴から逸脱することなく、様々な形で実施することができる。例えば、本発明をネットワーク環境で利用すべく、全部或いは一部のプログラムが他のコンピュータで実行されるようになっていてもかまわない。

第１の実施の形態の加熱又は冷却特性評価装置の概略構成を示すブロック図である。 評価装置における演算処理を説明するためのフローチャートである。 複数の熱電対が埋め込まれた被評価材の２次元断面を示す図である。 実施例における熱電対ＴＣ１、ＴＣ２、ＴＣ３、ＴＣ４の配置関係を示す図である。 実施例における解析結果を説明するための図である。 実施例における解析に用いたデータを説明するための図である。 実施例における解析に用いた実データを説明するための図である。 実施例における解析に用いた実データ詳細を説明するための図である。 実施例における実データによる解析結果を説明するための図である。 実施例における熱流束比較位置を説明するための図である。 実施例における別の方法１による解析結果を説明するための図である。 実施例における別の方法２による解析結果を説明するための図である。 第２の実施形態の反応容器の操業管理装置の概略構成を示すブロック図である。 複数の熱電対が埋め込まれた反応容器（加熱炉）の炉壁近くの２次元断面を示す図である。 実施例１における熱電対ＴＣ１、ＴＣ２の配置関係を示す図である。 実施例１における他の条件の熱電対ＴＣ１、ＴＣ２の配置関係を示す図である。 実施例１における解析に用いたデータを説明するための図である。 実施例１における解析に用いたデータを説明するための図である。 実施例１における他の条件の解析に用いたデータを説明するための図である。 実施例１における他の条件の解析に用いたデータを説明するための図である。 実施例１における解析方法を説明するための図である。 実施例１における他の条件の解析方法を説明するための図である。 実施例１における解析結果を説明するための図である。 実施例１における他の条件の解析結果を説明するための図である。 実施例１における実データによる解析結果を説明するための図である。 実施例１における実データによる解析結果を説明するための図である。 実施例１における実データによる解析結果を説明するための図である。 実施例１における実データによる解析結果を説明するための図である。 実施例１における他の条件の実データによる解析結果を説明するための図である。 実施例１における他の条件の実データによる解析結果を説明するための図である。 実施例１における他の条件の実データによる解析結果を説明するための図である。 実施例１における他の条件の実データによる解析結果を説明するための図である。 実施例１における他の解析方法の解析結果を説明するための図である。 実施例１における他の解析方法の実データによる解析結果を説明するための図である。 実施例１における他の解析方法の実データによる解析結果を説明するための図である。 実施例１における他の解析方法の実データによる解析結果を説明するための図である。 実施例２における解析方法を説明するための図である。 実施例２における解析結果を説明するための図である。 実施例２における他の条件の解析結果を説明するための図である。 実施例３における熱電対ＴＣ１、ＴＣ２、ＴＣ３、ＴＣ４の配置関係を示す図である。 実施例３における解析結果を説明するための図である。 実施例３における解析に用いたデータを説明するための図である。 実施例３における解析に用いた別のデータを説明するための図である。 実施例３における別の解析結果を説明するための図である。 実施例３における熱流束比較位置を説明するための図である。 実施例３における別の解析結果を説明するための図である。

符号の説明

１０１ 入力部
１０２ 演算部
１０３ 出力部

Claims (16)

1. 被評価材内部の複数の温度情報測定点において測定された温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、前記被評価材の内部又は表面における温度分布又は熱流束分布を演算する加熱又は冷却特性評価方法であって、
   前記温度情報測定点は前記被評価材内部に２次元又は３次元的に複数存在し、
   前記非定常熱伝導方程式は、密度ρ、比熱Ｃ_p、ｘ方向の熱伝導度ｋ_x、ｙ方向の熱伝導度ｋ_y、ｚ方向の熱伝導度ｋ_zとして、下式
   

   

   で表され、
   位置ベクトル（ｘ，ｙ，ｚ）、時間ｔとし、上記非定常熱伝導方程式を解析的に満足し厳密解を与えるｘ，ｙ，ｚ，ｔの関数Ｆを用いて、パラメータα_j,i、基準位置ベクトル（ｘ^j，ｙ^j，ｚ^j）、基準時間ｔⁱ、基準位置ベクトルの数Ｎ_j、基準時間の数Ｎ_iとして、Ｆをｊ，ｉに対応してＦ_j,iとし、前記非定常熱伝導方程式の解を、該非定常熱伝導方程式を解析的に満足し厳密解を与える関数として、下式
   

   

   により表現するものであって、
   ｋを温度情報測定位置、ｌを温度サンプリング時間とし、温度情報測定点において測定された温度情報ａ_k,lとして、前記パラメータα_j,iを、下式の連立方程式
   

   

   を用いて決めることを特徴とする加熱又は冷却特性評価方法。
2. 前記内外挿関数は、Ｘ、Ｙ、Ｚ、τ_x、τ_y、τ_z、Ａ_x、Ａ_y、Ａ_zを任意の定数として、下式
   

   

   の関係を有することを特徴とする請求項１に記載の加熱又は冷却特性評価方法。
3. 前記内外挿関数は、Ｘ、Ｙ、Ｚ、τ_xy、τ_z、Ａ_xy、Ａ_zを任意の定数として、下式
   

   

   の関係を有することを特徴とする請求項１に記載の加熱又は冷却特性評価方法。
4. 前記内外挿関数は、Ｘ、Ｙ、Ｚ、τ_xyzを任意の定数として、下式
   

   

   の関係を有することを特徴とする請求項１に記載の加熱又は冷却特性評価方法。
5. 熱流束又は温度の境界条件を与えた順問題解析を実施し、その解析で得られた数点の温度の時間変化情報に基づいて、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行い、順問題解析で付与した境界条件又は温度分布と一致するように、Ｘ、Ｙ、Ｚ、τ_x、τ_y、τ_z、Ａ_x、Ａ_y、Ａ_zの値を決定することを特徴とする請求項２に記載の加熱又は冷却特性評価方法。
6. 熱流束又は温度の境界条件を与えた順問題解析を実施し、その解析で得られた数点の温度の時間変化情報に基づいて、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行い、順問題解析で付与した点数Ｋの位置及び総サンプリング数Ｌの時間での熱流束ｑ_given,K,L又は温度Ｔ_given,K,Lと、逆問題解析で求めた対応する位置及び時間での熱流束ｑ_predict,K,L又は温度Ｔ_predict,K,Lが、下式
   

   

   又は、下式
   

   

   を極小化、又は、最小化するように、Ｘ、Ｙ、Ｚ、τ_x、τ_y、τ_z、Ａ_x、Ａ_y、Ａ_zの値を決定することを特徴とする請求項２に記載の加熱又は冷却特性評価方法。
7. 熱流束又は温度の境界条件を与えた順問題解析を実施し、その解析で得られた数点の温度の時間変化情報に基づいて、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行い、順問題解析で付与した境界条件又は温度分布と一致するように、Ｘ、Ｙ、Ｚ、τ_xy、τ_z、Ａ_xy、Ａ_zの値を決定することを特徴とする請求項３に記載の加熱又は冷却特性評価方法。
8. 熱流束又は温度の境界条件を与えた順問題解析を実施し、その解析で得られた数点の温度の時間変化情報に基づいて、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行い、順問題解析で付与した点数Ｋの位置及び総サンプリング数Ｌの時間での熱流束ｑ_given,K,L又は温度Ｔ_given,K,Lと、逆問題解析で求めた対応する位置及び時間での熱流束ｑ_predict,K,L又は温度Ｔ_predict,K,Lが、下式
   

   

   又は、下式
   

   

   を極小化、又は、最小化するように、Ｘ、Ｙ、Ｚ、τ_xy、τ_z、Ａ_xy、Ａ_zの値を決定することを特徴とする請求項３に記載の加熱又は冷却特性評価方法。
9. 熱流束又は温度の境界条件を与えた順問題解析を実施し、その解析で得られた数点の温度の時間変化情報に基づいて、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行い、順問題解析で付与した境界条件又は温度分布と一致するように、Ｘ、Ｙ、Ｚ、τ_xyz、の値を決定することを特徴とする請求項４に記載の加熱又は冷却特性評価方法。
10. 熱流束又は温度の境界条件を与えた順問題解析を実施し、その解析で得られた数点の温度の時間変化情報に基づいて、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行い、順問題解析で付与した点数Ｋの位置及び総サンプリング数Ｌの時間での熱流束ｑ_given,K,L又は温度Ｔ_given,K,Lと、逆問題解析で求めた対応する位置及び時間での熱流束ｑ_predict,K,L又は温度Ｔ_predict,K,Lが、下式
    

    

    又は、下式
    

    

    を極小化、又は、最小化するように、Ｘ、Ｙ、Ｚ、τ_xyz、の値を決定することを特徴とする請求項４に記載の加熱又は冷却特性評価方法。
11. 反応容器の壁内の複数の温度情報測定点において測定された温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、前記反応容器の内表面及び外表面における温度分布又は熱流束分布を演算する反応容器の操業管理方法であって、
    前記温度情報測定点は前記反応容器の壁内に２次元又は３次元的に複数存在し、
    前記非定常熱伝導方程式は、密度ρ、比熱Ｃ_p、ｘ方向の熱伝導度ｋ_x、ｙ方向の熱伝導度ｋ_y、ｚ方向の熱伝導度ｋ_zとして、下式
    

    

    で表され、
    位置ベクトル（ｘ，ｙ，ｚ）、時間ｔとし、上記非定常熱伝導方程式を解析的に満足し厳密解を与えるｘ，ｙ，ｚ，ｔの関数Ｆを用いて、パラメータα_j,i、基準位置ベクトル（ｘ^j，ｙ^j，ｚ^j）、基準時間ｔⁱ、基準位置ベクトルの数Ｎ_j、基準時間の数Ｎ_iとして、Ｆをｊ，ｉに対応してＦ_j,iとし、前記非定常熱伝導方程式の解を、該非定常熱伝導方程式を解析的に満足し厳密解を与える関数として、下式
    

    

    により表現するものであって、
    ｋを温度情報測定位置、ｌを温度サンプリング時間とし、温度情報測定点において測定された温度情報ａ_k,lとして、前記パラメータα_j,iを、下式の連立方程式
    

    

    を用いて決めることを特徴とする反応容器の操業管理方法。
12. 加熱又は冷却特性を評価するための装置であって、
    被評価材内部の複数の温度情報測定点において測定された温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、前記被評価材の内部又は表面における温度分布又は熱流束分布を演算する演算部を備えた加熱又は冷却特性評価装置であり、
    前記温度情報測定点は前記被評価材内部に２次元又は３次元的に複数存在し、
    前記非定常熱伝導方程式は、密度ρ、比熱Ｃ_p、ｘ方向の熱伝導度ｋ_x、ｙ方向の熱伝導度ｋ_y、ｚ方向の熱伝導度ｋ_zとして、下式
    

    

    で表され、
    位置ベクトル（ｘ，ｙ，ｚ）、時間ｔとし、上記非定常熱伝導方程式を解析的に満足し厳密解を与えるｘ，ｙ，ｚ，ｔの関数Ｆを用いて、パラメータα_j,i、基準位置ベクトル（ｘ^j，ｙ^j，ｚ^j）、基準時間ｔⁱ、基準位置ベクトルの数Ｎ_j、基準時間の数Ｎ_iとして、Ｆをｊ，ｉに対応してＦ_j,iとし、前記非定常熱伝導方程式の解を、該非定常熱伝導方程式を解析的に満足し厳密解を与える関数として、下式
    

    

    により表現するものであって、
    ｋを温度情報測定位置、ｌを温度サンプリング時間とし、温度情報測定点において測定された温度情報ａ_k,lとして、前記パラメータα_j,iを、下式の連立方程式
    

    

    を用いて決めることを特徴とする加熱又は冷却特性評価装置。
13. 反応容器の操業を管理するための操業管理装置であって、
    反応容器の壁内の複数の温度情報測定点において測定された温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、前記反応容器の内表面及び外表面における温度分布又は熱流束分布を演算する演算部を備えた反応容器の操業管理装置であり、
    前記温度情報測定点は前記反応容器の壁内に２次元又は３次元的に複数存在し、
    前記非定常熱伝導方程式は、密度ρ、比熱Ｃ_p、ｘ方向の熱伝導度ｋ_x、ｙ方向の熱伝導度ｋ_y、ｚ方向の熱伝導度ｋ_zとして、下式
    

    

    で表され、
    位置ベクトル（ｘ，ｙ，ｚ）、時間ｔとし、上記非定常熱伝導方程式を解析的に満足し厳密解を与えるｘ，ｙ，ｚ，ｔの関数Ｆを用いて、パラメータα_j,i、基準位置ベクトル（ｘ^j，ｙ^j，ｚ^j）、基準時間ｔⁱ、基準位置ベクトルの数Ｎ_j、基準時間の数Ｎ_iとして、Ｆをｊ，ｉに対応してＦ_j,iとし、前記非定常熱伝導方程式の解を、該非定常熱伝導方程式を解析的に満足し厳密解を与える関数として、下式
    

    

    により表現するものであって、
    ｋを温度情報測定位置、ｌを温度サンプリング時間とし、温度情報測定点において測定された温度情報ａ_k,lとして、前記パラメータα_j,iを、下式の連立方程式
    

    

    を用いて決めることを特徴とする反応容器の操業管理装置。
14. 加熱又は冷却特性を評価するための処理をコンピュータに実行させるコンピュータプログラムであって、
    被評価材内部の複数の温度情報測定点において測定された温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、前記被評価材の内部又は表面における温度分布又は熱流束分布を演算する処理をコンピュータに実行させるコンピュータプログラムであり、
    前記温度情報測定点は前記被評価材内部に２次元又は３次元的に複数存在し、
    前記非定常熱伝導方程式は、密度ρ、比熱Ｃ_p、ｘ方向の熱伝導度ｋ_x、ｙ方向の熱伝導度ｋ_y、ｚ方向の熱伝導度ｋ_zとして、下式
    

    

    で表され、
    位置ベクトル（ｘ，ｙ，ｚ）、時間ｔとし、上記非定常熱伝導方程式を解析的に満足し厳密解を与えるｘ，ｙ，ｚ，ｔの関数Ｆを用いて、パラメータα_j,i、基準位置ベクトル（ｘ^j，ｙ^j，ｚ^j）、基準時間ｔⁱ、基準位置ベクトルの数Ｎ_j、基準時間の数Ｎ_iとして、Ｆをｊ，ｉに対応してＦ_j,iとし、前記非定常熱伝導方程式の解を、該非定常熱伝導方程式を解析的に満足し厳密解を与える関数として、下式
    

    

    により表現するものであって、
    ｋを温度情報測定位置、ｌを温度サンプリング時間とし、温度情報測定点において測定された温度情報ａ_k,lとして、前記パラメータα_j,iを、下式の連立方程式
    

    

    を用いて決める演算処理をコンピュータに実行させることを特徴とするコンピュータプログラム。
15. 反応容器の操業を管理するための処理をコンピュータに実行させるコンピュータプログラムであって、
    反応容器の壁内の複数の温度情報測定点において測定された温度情報から、非定常熱伝導方程式を満たす内外挿関数を用いた逆問題解析を行うことにより、前記反応容器の内表面及び外表面における温度分布又は熱流束分布を演算する処理をコンピュータに実行させるコンピュータプログラムであり、
    前記温度情報測定点は前記反応容器の壁内に２次元又は３次元的に複数存在し、
    前記非定常熱伝導方程式は、密度ρ、比熱Ｃ_p、ｘ方向の熱伝導度ｋ_x、ｙ方向の熱伝導度ｋ_y、ｚ方向の熱伝導度ｋ_zとして、下式
    

    

    で表され、
    位置ベクトル（ｘ，ｙ，ｚ）、時間ｔとし、上記非定常熱伝導方程式を解析的に満足し厳密解を与えるｘ，ｙ，ｚ，ｔの関数Ｆを用いて、パラメータα_j,i、基準位置ベクトル（ｘ^j，ｙ^j，ｚ^j）、基準時間ｔⁱ、基準位置ベクトルの数Ｎ_j、基準時間の数Ｎ_iとして、Ｆをｊ，ｉに対応してＦ_j,iとし、前記非定常熱伝導方程式の解を、該非定常熱伝導方程式を解析的に満足し厳密解を与える関数として、下式
    

    

    により表現するものであって、
    ｋを温度情報測定位置、ｌを温度サンプリング時間とし、温度情報測定点において測定された温度情報ａ_k,lとして、前記パラメータα_j,iを、下式の連立方程式
    

    

    を用いて決める演算処理をコンピュータに実行させることを特徴とするコンピュータプログラム。
16. 請求項１４又は１５に記載のコンピュータプログラムを記録したことを特徴とするコンピュータ読み取り可能な記録媒体。

Images