JP4299554B2 - Inverse problem analysis method, apparatus, computer program, and computer-readable storage medium - Google Patents

Description

【００２２】
数１において、ρは反応容器の材料の密度、Ｃ_pは反応容器の材料の比熱、Ｔは反応容器内部の温度の計算値、ｔは時間、ｋは反応容器の材料の熱伝導度を表す。
【００２３】
熱伝導逆問題解析というのは、計算領域を支配する非定常熱伝導方程式を基にして、領域内部の温度を既知として、領域境界での温度や熱流束等の境界条件を推定することをいう。これに対して、熱伝導順問題解析というのは、既知である領域境界での温度や熱流束等の境界条件から領域内部の温度を推定することをいう。
【００２４】
２次元逆問題解析の手法の例としては、例えば、本願出願人が特開２００２−２０６９５８号公報に開示したものがあり、この手法はそのまま１次元逆問題解析へも適用できる。また、１次元逆問題解析の例として、Beckらにより提案された解析手法が知られている（Beckその他、Inverse Heat Conduction,1985,Wiley,New York）。
【００２５】
また、逆問題解析の最近の手法として、カルマンフィルター理論や、射影フィルタ理論等の確率的推定法を適用することも考えられる。この手法は、現状では、上記数１の左辺をゼロと置いた、定常熱伝導方程式（観測方程式）への適用が検討されているが、非定常項を含めて適切に観測行列を構成できれば、同様の逆問題解析ができる可能性がある。この定常微分方程式への、確率推定法の適用例としては、登坂その他、「逆問題の数理と解法・偏微分方程式の逆解析」（東京大学出版会（1999））に詳しい。
【００２６】
本実施の形態では、逆問題解析の手法として上記特開２００２−２０６９５８号公報に開示した考え方を用いる。すなわち、下記の数２に示すように、ある１次元方向（図２に示す１ａ→１ｂ→１ｃや１ｄ→１ｅ等）に配置された各熱電対で測定された温度Ｙと、反応容器の内表面及び外表面における熱流束の仮定値から非定常熱伝導方程式により算出された各熱電対位置での温度Ｔとの差の二乗の和が最小となる仮定値を反応容器の内表面及び外表面における熱流束として求める。なお、Ｊは熱電対の数を表す。
【００２７】
【数２】

【００２８】
このように複数の熱電対位置での温度Ｔ、Ｙを完全に一致させるような解（反応容器の内表面及び外表面における熱流束）を求めるのではなく、最小二乗的に満たすような解を求めることにより、現実的な熱流束変化の推定が可能となる。その理由は、測定温度には様々な測定誤差要因が含まれるため、完全に一致させることは実用的に意味がないといえるからである。
【００２９】
なお、計算を安定化させるために、正則化項を付加するようにしてもよい。下記の数３には、０次の正則化項の例を示す。ｐは推定熱流束の分割数の数であり、α₀は経験値から得られる正則化パラメータである。
【００３０】
【数３】

【００３１】
以下に、より具体的に、複数の熱電対位置での温度Ｙを既知として、反応容器の内表面及び外表面における熱流束を推定する定式化と、計算手続きの一例を示す。
【００３２】
下記の数４のＳ_mは全体の目的関数を表し、下記の数５は、実測温度Ｙと計算温度Ｔの偏差を表す目的関数を示す。下記の数６は、計算を安定化するために付加した目的関数であり、空間分割方向の値の急激な変化を抑える働きがある。数６中のα₀やα₁は、一定の経験値から得られる正則化パラメータである。
【００３３】
【数４】

【００３４】
【数５】

【００３５】
【数６】

【００３６】
上記数５では、ある熱電対で測定された温度Ｙと、熱流束の仮定値から熱伝導方程式モデルにより算出された温度Ｔの差の二乗が最小となるように目的関数を設定している。また、上記数６では、温度測定誤差があっても解が安定するように空間方向の正則化を施す目的関数を設定している。そして、数４を全体の目的関数として、下記の数７に示すように、未知である熱流束分割領域に対して極小点を探す。
【００３７】
【数７】

【００３８】
ここで、数８に示すように、解を安定させる目的で、各時間ステップの熱流束値が、一定の未来時間まで不変であると仮定する。時間ステップは、対象とする材料の熱物性・形状等によって変わる。数８のｑは熱流束を示し、ｍ時間ステップにおける熱流束ｑ_mから、将来時間ｍ＋ｒ−１時間ステップにおける熱流束ｑ_m+r-1が一定であると仮定している。
【００３９】
【数８】

【００４０】
そして、数７の極小化を、数８の仮定を用いて展開すると、数９に示すように、マトリクス形に展開することができる。
【００４１】
【数９】

【００４２】
数９のＸ^TＸは数４の右辺第１項から導かれ、Ｘ^TＸに続く２項（α₀Ｈ₀ ^TＨ₀＋α₁Ｈ₁ ^TＨ₁）は、数４の右辺第２項から導かれる（上付のＴは、転置行列を表す）。Ｘの構成は、補足式数１０として下部に、Ｘ_j,i,kとして示している。ここで、時間方向の分割数を示すｉは、最大Ｍ時間ステップまで変化し、熱電対の数を示すｊは、最大Ｊ個まで変化して、熱流束分布の分割数を示すｋは、最大ｐまで変化する。なお、数９の上付の＊は、繰り返し収束計算での参照値であることを示しており、Ｔ^*は温度参照値、ｑ^*は熱流束参照値である。１次元の場合は、両端の境界条件を推定するので、熱流束分布の分割数ｋは、最大ｐ＝２である。
【００４３】
【数１０】

【００４４】
数９は、温度変化が起きた場合の熱流束の変化を推定する連立方程式であり、各時間ステップにおいて、この数９を用いて両端の熱流束ｑを求める。まずは、前時間ステップでの熱電対位置での計算温度を初期Ｔ^*とし、数９によりｑを求める。このｑを、並行して計算している順問題熱伝導方程式モデルの境界条件として与え、温度分布を計算する。ここで求めた温度計算値を、次の温度参照値Ｔ^*として、ｑを再修正する（数９に代入してｑを再び求める）。この操作を、数５が一定残差以下になる（収束）まで、ｑとＴ^*の修正を繰り返し、各時間ステップにおける両端の熱流束（最終的なｑ）を求めていく。この計算手続きを繰り返すことにより、両端の熱流束ｑの変化を、２つ同時に推定することが可能となる。
【００４５】
数１０は、一種の感度行列を表しており、端的に言うと、境界端点での熱流束ｑの単位変化に対する熱電対位置での計算温度Ｔの変化の大きさの比率を示している。数１０は、逆解析と同時に計算している順問題計算によって、各時間ステップにおいて、単位時間ステップあたりの値の計算が可能である。
【００４６】
以下、１次元の逆問題解析を例にして、より望ましい解法について説明する。上述のように、２つの端面（反応容器の内表面及び外表面）の熱流束を未知の境界条件とした１次元逆問題を構成（定式化）しても、原理上は解を求めることができる。
【００４７】
ただし、熱電対の数や材料の熱物性条件等によって多解となる場合があり、計算が不安定となる可能性がある。その理由の一つは、「未知両端面の熱流束差」の組み合わせを適当に選ぶことができれば、離散的な温度測定点の温度変化を表現する熱流束の組み合わせは無数に存在する可能性があるためである。特に、熱伝導度の低い物質の場合、表面温度が極端に大きくなったり、小さくなったりしてしまうような境界条件を推定してしまう場合でも、離散測定点の温度の変化だけを再現すれば、一つの解として認識してしまうことも起こり得る。これは、現実の現象としてはあり得ないことであるばかりでなく、逆問題計算を非常に不安定なものとする。
【００４８】
また、実際の問題として、逆問題解析を開始する時の熱電対の温度（離散測定点の温度）は既知として与えられるが、その他の解析領域での温度分布の初期条件は不明であることが一般的である。このため、任意に与えた仮初期温度分布から計算を始め、計算ステップを進める中で、実際の温度分布を探索・推定し、妥当な温度分布へと徐々に修正しながら、安定的に計算を進めていけるような計算ロジックにすることが求められる（ここで言う温度分布とは、例えば、逆問題解析の計算手続きの中で、上記数９の解を修正するために並行して計算している順問題熱伝導方程式モデルの計算値である）。このように、初期温度分布が不確定であることも、逆問題計算を不安定なものとする大きな要因の一つとなる。
【００４９】
以上のことは、逆問題を安定化するためには、逆問題解析の過程で、ある程度の表面温度の目安（拘束条件）を与える必要性があることを示しているといえる。この考え方に基づき、拘束条件を適当に与える手法を、図３のフローチャートを参照して説明する。
【００５０】
まず、反応容器の内表面及び外表面のいずれか片側、ここでは外表面における熱流束として仮の熱流束ｑを与える。この仮の熱流束ｑの与え方として、熱伝達率ｈと参照温度Ｔ_bとを用いて、
ｑ＝ｈ（Ｔ_surf−Ｔ_b）
として与える（ステップＳ２０１）。
【００５１】
Ｔ_surfは未知境界、ここでは反応容器の外表面における温度を示している。この表面温度Ｔ_surfは、逆問題解析の過程で熱流束の値を修正するために、通常は順問題解析も同時に行うが、この順問題解析で求めた表面温度に相当する。
【００５２】
また、参照温度Ｔ_bは反応容器の内部及び内外表面以外での温度である。本実施の形態では、反応容器の冷却条件、例えば、水冷ならば水温等に基づいて定めるようにしている。
【００５３】
結果として、上式の左辺である熱流束ｑをあたかも既知の熱流束情報として与えることができる。このように仮の熱流束情報を与えることで、熱伝達率ｈと参照温度Ｔ_bという２つの拘束条件を与えることとなり、任意の熱流束を与えるのに比べて物理的な妥当性を確保して、極端な温度分布が生じることを防ぐことが可能となる。
【００５４】
次に、反応容器の外表面における仮の熱流束ｑ（＝ｈ（Ｔ_surf−Ｔ_b））を与えて、上記数２、又は、数５に示した温度Ｔ、Ｙの差の二乗の和が最小となる反応容器の内表面における熱流束を、反応容器の内表面における仮の熱流束として算出する（ステップＳ２０２）。このステップは、逆問題解析のメインの計算手続きであり、具体的な解法の一つとして、数４から数９に示した定式化と計算手続きが、そのまま適用できる。この場合では、数９を解く際に、反応容器の外表面における仮の熱流束ｑ（＝ｈ（Ｔ_surf−Ｔ_b））は既知として与え、反応容器の内表面における仮の熱流束を未知として解くことを意味する。
【００５５】
ここで、上記のように片側（反応容器の外表面）の熱流束情報を与えて、逆問題解析により求めた反対側（反応容器の内表面）の熱流束は、一つの解の可能性を示しているに過ぎない。また、既知と仮定した熱伝達率ｈや参照温度Ｔ_bも概算値であり、本来ならば未知の値である。
【００５６】
そこで、熱伝達率ｈ及び外部参照温度Ｔ_bの両方或いはいずれかを数点変化させて、すなわち、反応容器の外表面における仮の熱流束ｑの値を数点（Ｋ点）振って、反応容器の外表面における仮の熱流束ｑと、各仮の熱流束情報ｑを与えたとき温度Ｔ、Ｙの差の二乗の和が最小となる反応容器の内表面における熱流束との組み合わせをＫ個得る（ステップＳ２０３）。
【００５７】
そして、下記の数１１に示すように、反応容器の外表面における仮の熱流束ｑと、各仮の熱流束情報ｑに対応して得られた反応容器の内表面における熱流束とのＫ個の組み合わせのうち、温度Ｔ、Ｙの差の二乗の値が最も小さくなる組み合わせを選び出し、その組み合わせを反応容器の内表面及び外表面における熱流束とする（ステップＳ２０４）。
【００５８】
【数１１】

【００５９】
上式の大括弧の中は、片側の熱流束を既知として逆問題解析した１ケースの計算結果を示し、その計算をＫケース計算した中から更に最小二乗差の最も小さな結果を選び出すことを意味する。
【００６０】
この手続を、各時間ステップにおいて繰り返し行うことにより、反応容器の内表面及び外表面における熱流束経時変化を逐次同時計算していくことができる。
【００６１】
以上述べたように、反応容器の内表面及び外表面における熱流束変化を同時に求めるような１次元逆問題解析を安定して実行することができる。そして、反応容器の内表面及び外表面における温度変化や熱流束変化を同時推定できれば、例えば、ある温度測定点における温度変動が、反応容器の内表面における熱流束変化によるものなのか、反応容器外に設置された冷却装置の接触不良等によって引き起こされるような反応容器の外表面における熱流束変化によるものかを区別するようなことが可能となる。
【００６２】
上記手法は１次元逆問題解析に適用すると簡便であり、実際問題として有効である場合が多い。その理由は、一般的には、反応容器の上端と下端とは断熱条件（対称）とする場合が多く、実用的にも問題ないからである。
【００６３】
したがって、図２の破線で区切られた範囲での厚み方向１次元を仮定して逆問題解析し、その結果を上下方向に組み合わせることで、２次元化することも可能である。
【００６４】
より厳密に図２の上下方向の熱流れも考慮したい場合には、２次元逆問題解析が必要である。このような２次元解析は、図２の左右両端部の熱流束分割を上方向に細かくして、これらの熱電対位置での温度を最小二乗的に最小な熱流束分布を求めることと等価であり、上述した特開２００２−２０６９５８号公報に開示した逆問題定式化と同様の手法に従って本発明を適用すればよいこととなる。
【００６５】
この場合に、図２の上端下端の熱流束に関しては、未知としても、既知としても構わないが、計算の安定性を考慮すると、物理的な考察から適当な熱流束（例えば、断熱等）を与えて既知とした方が望ましい。
【００６６】
同様の考えに基づいて、３次元解析への拡張も容易に行うことができる。
【００６７】
ここで、例えば上述した逆問題解析により反応容器の表面における熱流束変化の「現在」の変化を推定するためには、できるだけ近い過去の反応容器内の状況を推定する必要がある。
【００６８】
そのためには、逆問題計算の単位時間ステップを短くすることが考えられる。逆問題計算の単位時間ステップが長いと、少なくともその時間ステップ分、大きく過去に遡った変化を推定してしまうことになる。また、逆問題計算は、計算時間ステップで時間平均した計算結果であるので、復元した熱流束変化も、鈍った変化を捉えることになり、時間ステップ以下の急激な変化は捉えられない。
【００６９】
ところが、逆問題計算の単位時間ステップを短くするということは、その短い時間の間に、熱流束の変化位置（例えば、溶銑と接触する高炉の内表面）から熱電対位置まで伝わってくる小さな温度変化を捉えなければならないことを意味する。特に、熱電対位置と熱流束の変化位置とが離れていて、その間の物質の熱伝導率が小さい場合、短時間での温度の動きは非常に小さなものとなる。
【００７０】
このような小さな温度変化から、熱流束変化による温度変化を理論的に抽出して、その熱流束変化を復元する（短時間の逆問題解析）のためには、温度測定方法や温度データ処理方法に工夫が必要である。
【００７１】
そこで、本願発明者らは、熱電対で測定される温度データを活用しつつ、短い時間ステップで、できるだけ近い過去の熱流束分布を推定できるようにすべく、鋭意検討を重ねた。
【００７２】
高炉等の操業の場合、温度の変化だけを見ていることもあって、大局的な温度の流れを見るに際しては、小数点以下の温度データを切り捨てた値で評価することが一般的である。
【００７３】
また、現在の温度の時間変化を細かく見る場合でも、熱電対の接点電圧変化から温度に変換する温度測定装置のメーカ保証範囲がほとんどの場合小数点１桁までであることもあり、小数点１桁までの温度データを使用することが多い。すなわち、採取して、記録している温度データは、小数点０桁か１桁といった粗いものである。
【００７４】
ところが、非常に微妙な温度変化を引き起こす原因となる、遠く離れた位置での境界条件変化（熱流束変化）を推定する（逆問題解析）ためには、温度変化の小数点以下の挙動を推定することが非常に重要となる。
【００７５】
すなわち、逆問題解析においては、温度そのものの精度よりは、温度の時間推移・温度変化の流れが重要である。温度そのものの精度だけを追求して、階段状に温度が変化してしまうと、逆問題解析を不安定なものとしてしまう。特に、短い時間ステップでの逆問題解析を指向すると、その時間ステップでは考えられないような階段状の温度変化が生じてしまうこともある。そのため、この階段状の温度変化で逆問題解析して熱流束解が得られた場合でも、物理的にありえない解に至る可能性を秘めている。実際、そのほとんどの場合、解が発散してしまい、逆問題解析の続行が不可能となる。
【００７６】
合理的な方法で、小数点以下の桁数の少ない温度データから、人工的に小数点以下の桁数を増やした温度データを生成することは、階段状の時間変化を滑らかにする上で有効な手法である。
【００７７】
その一つの手法として、演算部１０２において、逆問題解析での時間ステップより短いサンプリング時間で小数点以下０桁や１桁の温度データを採取し、それらを時間平均して、逆問題解析での計算時間ステップで用いる温度データの代表値とする。例えば、逆問題解析での時間ステップ（１時間）より短い５分間隔でサンプリングした小桁数の温度を、１時間範囲で単純平均して、１時間ステップの温度データの代表値とする。これにより、見かけ上、１時間温度データの小数点数を増やすことが可能となる。この手法を用いると、熱電対で測定される精度が保証されている小数点以下０桁、小数点以下１桁の温度データ（熱電対で測定されたもの）に基づいて、逆問題計算に用いる温度データの小数点以下の桁数を増やすことが可能である。
【００７８】
また、別の手法により、逆問題計算に用いる温度データの小数点以下の桁数を増やすようにしてもよい。例えば、熱電対において、温度データそのものの保証範囲は小数点１桁である場合でも、実際には熱電対接点での電圧から温度へと換算する変換式に基づいて計算するか、キャリブレーションを取った電圧−温度換算表から温度に変換するかの方法が一般的である。
【００７９】
前者の方法では、多数桁の温度計算値を四捨五入等して小数点以下１桁にしているわけであり、実際には計算機の能力に応じた小数点桁数まで求めることが可能である。例えば、Ｋタイプの熱電対での電圧から温度に変換する変換式として、下式に示すようなものがある。
Ｔ＝ａ＋ｂ・Ｘ＋ｃ・Ｘ²＋ｄ・Ｘ³＋ｅ・Ｘ⁴＋ｆ・Ｘ⁵＋ｇ・Ｘ⁶＋ｈ・Ｘ⁷＋ｉ・Ｘ⁸
Ｔ：温度 [℃]
Ｘ：接点電圧 [Ｖ]
ａ＝０．２２６５８４６０２
ｂ＝２４１５２．１０９００
ｃ＝６７２３３．４２４８
ｄ＝２２１０３４０．６８２
ｅ＝−８６０９６３９．９
ｆ＝４．８３５０６×１０¹⁰
ｇ＝−１．１８４５２×１０¹²
ｈ＝１．３８６９０×１０¹³
ｉ＝−６．３３７０８×１０¹³
上記のような変換式を用いて、接点電圧から温度に直接換算するようにすれば、簡単に温度データの小数点以下の数値を求めることが可能である。なお、ここでいう接点電圧とは、冷接点温度補償を経た真の接点電圧であることはいうまでもない。
【００８０】
（第２の実施の形態）
上記第１の実施の形態で説明したように小数点以下の桁数を増やす以外にも、図４に示すように、入力部１０１にフィルタ１０１ａを設けておき、温度データ前処理にフィルタ処理を施すことが挙げられる。すなわち、フィルタ１０１ａを用いて、熱電対で測定された温度データを修正し、その修正後の温度データを逆問題解析に用いる。フィルタを用いると、元の温度データが小数点以下の桁数の少ないものであっても、見かけ上、温度データに小数点以下の桁を生成することができるので、温度の変化が滑らかになって逆問題解析の計算が非常に安定する。
【００８１】
理論的に最も有効と思われるのは、ローパスフィルタを用いることである。ローパスフィルタは、等時間間隔でサンプリングした実測データに混入している高周波ノイズによる変動の影響を除去するために、実測データを補正するものである。
【００８２】
逆問題解析においてローパスフィルタを用いて温度データを修正することが有効な理由は、熱電対から遠く離れた位置での熱流束の変化は、耐火物内を熱拡散現象によって伝わってくる間に鈍化して、熱電対位置まで伝わった時には、温度変化が非常になまっている（低周波の変化が伝達される）と考えられるためである。
【００８３】
これは、遠く離れた炉内内部の変化によって引き起こされた熱電対位置での温度変化は、時間的に低周波な変化であることを意味する。ところが、実際の熱電対の接点電圧変化（温度変換前の元信号）には、他の原因による高周波のノイズが重畳されている。そこで、実際の測定温度変化から高周波ノイズを取り除く処理を施した温度データを逆問題解析で用いることは、物理的に非常に意味のある手法となる。したがって、この手法により、逆問題解析の計算安定性を向上させることができるのみならず、熱流束の推定精度を向上させることもできる。
【００８４】
ローパスフィルタ処理は、先に述べた時間平均をとったり、変換式を使ったりして桁数を増やした温度データは勿論のこと、その他なんらかの手法で小数点以下の桁数を増やした温度データや、元の小数点以下の桁数の少ない測定温度データに対して付与しても非常に有効である。すなわち、第１の実施の形態で説明した小数点以下の桁数を増やすことと、ローパスフィルタ処理を施すことを併用してもよく、それにより相乗的に温度変化の平滑化を図ることができる。
【００８５】
フィルタについては多種多様な種類が考えられるが、例えば、R.W.HammingのDigital filters,Dover publications Inc.1998等を参照するとよい。その他、参考資料は多数存在する。
【００８６】
ローパスフィルタの古典的な手法としては、等時間間隔でサンプリングした測定温度データＹを用いて中間点の温度を修正するものが知られている。その修正のための補正式として、下記の数１２に示す式（１）〜（３）の式が導出される。下付き文字のｉは時間のサンプリング回数を示しており、現在時刻ｉに対して、前後の時間の測定温度により修正を掛けることを意味する。
【００８７】
【数１２】

【００８８】
式（１）〜（３）はそれぞれ３項法、５項法、７項法の式を示しており、修正に用いる測定温度データ数が増えるほど、ローパスフィルタとしての性能が向上する。すなわち、より低周波の信号のみを抽出することができる。基本的には、各項の係数で表現できるので、ウインドウズ表示で表現すると、それぞれ式（４）〜（６）のようになる。
【００８９】
一般的に、項数が増えれば増えるほど、ローパスフィルタとしての性能は向上し、フィルタ設計の自由度が増える。つまり、多項フィルタにより、高周波の外乱が広範囲に除去することができるので、逆問題解析の計算は極めて安定する。その反面、項数が増えると、ある時刻点の修正温度値を得るために、その前後の複数時間点の測定温度を使用しなければならず、これは実質的に過去に遡っていることになるため、「逆問題解析で現在の熱流束を推定する」という場合の「現在」そのものが、ローパスフィルタによって引き戻されることを意味する。したがって、ローパスフィルタによる逆問題解析の安定性を確保することと、「現在」の推定時間をできるだけ今に近づけることの関係はトレードオフなので、両者の関係を鑑みながら、適当な組み合わせを選ぶ必要がある。
【００９０】
こういう観点からいえば、上記第１の実施の形態で説明した手法により、できるだけ小数点以下の桁数を増やした温度データを用い、階段状の温度変化を抑えた上で、項数の少ないローパスフィルタを用いることが、「現在」の熱流束を推定する上で有効となる。
【００９１】
（実施例）
図５には、高炉炉底の底盤に設置されているカーボン煉瓦内に埋め込まれた熱電対を用いて、１次元の非定常熱伝導モデルを仮定して、逆問題解析を試みた例を模式的に示す。熱電対は、冷却面側に偏って２本埋め込まれており、ＴＣ１は高温側熱電対、ＴＣ２は低温側熱電対を示す。これら熱電対ＴＣ１、ＴＣ２で測定される温度変化から、上述した逆問題解析により、高温熱流束面での非定常熱流束ｑ１と、冷却面での非定常熱流束ｑ２とを同時に推定する。高温熱流束面位置は冷却面から４．０ｍ奥の定点である。仮定したカーボン煉瓦の熱物性値は、比熱Ｃｐ＝７１２Ｊ／（ｋｇ・Ｋ）、密度ρ＝２３００ｋｇ／ｍ³、熱伝導度ｋ＝２１．２Ｗ／（ｍ・Ｋ）である。
【００９２】
図６（ａ）に、本実施例で得られた高温側熱流束ｑ１の推定結果を示す。図６の横軸は日付を表し、９月１目から１０月１目までの結果を示している。実際には、９月１目以前より計算を始めて、計算終了は１０月１日である。なお、逆問題解析の時間ステップは、８時間ステップである。
【００９３】
ケース１は、熱電対により測定される温度データの小数点１桁の数値を四捨五入して、小数点以下０桁としたデータに対して、フィルタ処理を実行した場合の計算結果である。ケース２は、温度データの小数点２桁の数値を四捨五入して、小数点以下１桁としたデータに対して、フィルタ処理を実行した場合の計算結果である。
【００９４】
図６（ｂ）、（ｃ）には、高温側熱電対ＴＣ１及び低温側熱電対ＴＣ２のフィルタ処理後の温度データ（逆問題解析に用いる温度データ）をプロットしたグラフを示す。また、図６（ｄ）には、熱電対ＴＣ１、ＴＣ２それぞれの位置でのケース１とケース２の温度差（ケース１−ケース２）を示す。
【００９５】
フィルタについては、過去の温度データには極力多くの項数のローパスフィルタをかけて計算を安定化すべく、最大２１項のスペンサーの式（式（１０））を用い、現在に近づくにつれて徐々に項数を減らしていく手法を用いている。具体的には、ウインドウズ形式の表記で、式（４）〜（１０）を用いている。即ち、１０月１日に近づくにつれて項数が減少し、最終的には式（４）を用いており、最後の温度データにはフィルタ処理を行っていない。
【００９６】
ケース１、ケース２のいずれの場合でも、過去の解析に関してはフィルタ効果で非常に安定化しているが、現在に近づくにつれて、徐々に不安定化して、最終的な推定熱流束は、現実的には考えられないほどの大きな振動を起こしている。これは、ローパスフィルタの効果が徐々に低下して、計算が不安定になっていることを示すと考えられる。
【００９７】
ケース１とケース２を比較すると、ケース２の方が最終的な振動の大きさを小さく抑えることができていることが分かる。これは、小数点以下の桁数が多い温度データの方が、項数の少ないローパスフィルタでも、比較的安定した計算が可能であることを示している。
【００９８】
更に、いずれの場合も、最後の推定結果で振動を起こしていることから、同じ桁数の温度データで比較した場合、ローパスフィルタの項数を増やして、温度データを極力平滑化することが、安定した解を得る上で有効な手段であることも示している。特に、式（４）〜（１０）の形式のローパスフィルタを用いた場合、式（４）を用いても、最後のデータにはフィルタがかからないので、元の温度データの小数点以下の値は多い方が、見かけ上平滑化されやすいので、逆問題解析解が安定する。
【００９９】
ケース１及びケース２において、温度データにフィルタを全くかけないで逆問題解析を試みたが、解が発散してしまい、計算を進めることができなかった。
【０１００】
図７には、熱電対の設置深さの異なるケースについて、図５と同様の模式図を示す。図８（ａ）に、逆問題解析の計算時間ステップを８時間とした場合と、６時間とした場合とについて、高温側熱流束ｑ１の推定結果を示す。これらの計算結果ついては、もっと長い期間の計算結果の中での９月１目から１０月１日までを抽出して示す。
【０１０１】
計算に用いた温度データは小数点以下１桁のものを用い、更にフィルタとして式（１０）を用いて平滑化を行っている。図８（ｂ）、（ｃ）には、その平滑化後の高温側熱電対ＴＣ１及び低温側熱電対ＴＣ２の温度データをプロットしたグラフを示す。また、図８（ｄ）は、実質的には、定常の熱流束に対応する、高温側熱電対温度と低温側熱電対温度との差（ＴＣ１−ＴＣ２）を示す。
【０１０２】
図８（ｄ）の疑似定常熱流束の推移と、図８（ａ）の逆問題解析により推定した非定常熱流束のピーク同志の対応関係を比較して見ると、直感的ながら、逆問題解析が遅れ時間を補正して、非定常熱流束の方が早めにピークが表れていることが分かる（図８（ａ）中の破線で示す）。
【０１０３】
図８（ｄ）のグラフの疑似定常熱流束は、熱電対位置周辺を通過している平均的な熱流束と解釈でき、４ｍ奥の高温面での値を推定した非定常熱流束の値と比較して必ずしもピーク同志が対応しているとは限らないが、この場合は、概略２〜３日のピーク遅れが観察される。換言すると、逆問題解析による非定常熱流束は、従来の定常法よりも、概略２〜３日早めに炉内の動きを検知していることになる。
【０１０４】
また、図８（ａ）のグラフでの８時間ステップ計算と６時間ステップ計算を比較すると、６時間ステップの計算の方がより輪郭がはっきりした動きを示している。時間ステップを長くすると、その時間内の平均化した動きしか捉えることができないので、極力時間ステップを短くとった方が、実際の内部の動きを表現できる可能性が高い。
【０１０５】
（他の実施の形態）
以上説明した逆問題解析装置は、コンピュータのＣＰＵ或いはＭＰＵ、ＲＡＭ、ＲＯＭ、ＲＡＭ等で構成されるものであり、上述のようにＲＡＭやＲＯＭ等に記憶されたプログラムが動作することによって実現される。
【０１０６】
したがって、プログラム自体が上述した実施の形態の機能を実現することになり、本発明を構成する。プログラムの伝送媒体としては、プログラム情報を搬送波として伝搬させて供給するためのコンピュータネットワーク（ＬＡＮ、インターネット等のＷＡＮ、無線通信ネットワーク等）システムにおける通信媒体（光ファイバ等の有線回線や無線回線等）を用いることができる。
【０１０７】
さらに、上記プログラムをコンピュータに供給するための手段、例えばかかるプログラムを格納した記憶媒体は本発明を構成する。かかる記憶媒体としては、例えばフレキシブルディスク、ハードディスク、光ディスク、光磁気ディスク、ＣＤ−ＲＯＭ、磁気テープ、不揮発性のメモリカード、ＲＯＭ等を用いることができる。
【０１０８】
なお、上記実施の形態において示した各部の形状及び構造は、何れも本発明を実施するにあたっての具体化のほんの一例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその精神、又はその主要な特徴から逸脱することなく、様々な形で実施することができる。例えば、本発明をネットワーク環境で利用すべく、一部のプログラムが他のコンピュータで実行されるようになっていてもかまわない。
【０１０９】
【発明の効果】
以上述べたように本発明によれば、反応容器の壁内部に配置された温度測定点で測定される温度データに基づいて、非定常熱伝導方程式を用いた逆問題解析を行うことにより、上記反応容器の表面における温度或いは熱流束を求めるに際して、例えば、できるだけ近い過去の温度或いは熱流束の変化を推定するために逆問題解析での時間ステップを短くするような場合でも、温度データの小数点以下の桁数を増やしたり、フィルタ処理を施したりすることにより、小さな温度変化を捉えることが可能となり、逆問題解析を安定化させることができる。
【図面の簡単な説明】
【図１】第１の実施での形態の逆問題解析装置の構成を示す図である。
【図２】複数の熱電対が埋め込まれた反応容器（高炉）の炉壁近くの２次元断面を示す図である。
【図３】逆問題解析における処理を説明するためのフローチャートである。
【図４】第２の実施での形態の逆問題解析装置の構成を示す図である。
【図５】実施例における熱電対ＴＣ１、ＴＣ２の配置関係を示す図である。
【図６】実施例における解析結果を説明するための図である。
【図７】他の実施例における熱電対ＴＣ１、ＴＣ２の配置関係を示す図である。
【図８】他の実施例における解析結果を説明するための図である。
【符号の説明】
１０１ 入力部
１０１ａ フィルタ
１０２ 演算部
１０３ 逆問題解析部
１０４ 出力部[0001]
BACKGROUND OF THE INVENTION
The present invention relates to an inverse problem analysis method, apparatus, and computer suitable for estimating the temperature or heat flux at the surface of a reaction vessel with a temperature change reaction such as a blast furnace, a steel heating furnace by combustion, and a coal gasification reactor. The present invention relates to a program and a computer-readable storage medium.
[0002]
[Prior art]
For example, as a technique for monitoring the molten metal flow at the bottom of a blast furnace furnace, a thermocouple is embedded in a brick at the bottom of the furnace, and a determination is made from the temperature measured by the thermocouple.
[0003]
However, this temperature change does not reflect the current hot water flow because it has passed a delay time after the change of the hot water flow in the furnace. This is because temperature abnormalities in the blast furnace are transmitted to the surface of the blast furnace as heat flux changes, and then heat accumulates inside the bricks of the blast furnace, while heat is gradually transferred by the heat conduction effect, and finally the temperature changes to the thermocouple In principle, the heat conduction phenomenon in a solid having a heat capacity has a time delay (unsteadiness).
[0004]
On the other hand, as disclosed in Patent Document 1, for example, the heat conduction phenomenon inside the wall of a reaction vessel such as a blast furnace is considered as an unsteady heat conduction inverse problem, and the reaction is caused by a temperature change at a thermocouple. A method for estimating the heat flux change on the inner surface of the container has been proposed.
[0005]
[Patent Document 1]
JP 2001-234217 A
[0006]
[Problems to be solved by the invention]
By the way, when estimating the heat flux change on the surface of a reaction vessel such as a blast furnace by inverse problem analysis, it is necessary to stabilize the inverse problem analysis.
[0007]
The present invention has been made in view of the above points, and performs inverse problem analysis using an unsteady heat conduction equation based on temperature data measured by a thermocouple or the like embedded in the wall of a reaction vessel. Thus, the object of the present invention is to stabilize the inverse problem analysis when obtaining the temperature or heat flux on the surface of the reaction vessel.
[0008]
[Means for Solving the Problems]
  In the present invention, the reaction vessel is analyzed by performing inverse problem analysis using an unsteady heat conduction equation based on temperature data measured at a temperature measurement point arranged inside the wall of the reaction vessel with a temperature change reaction. When calculating the temperature or heat flux at the surface of the surface, increase the number of digits after the decimal point of the temperature data measured at the temperature measurement point, and use the temperature data with the increased number of digits after the decimal point for the inverse problem analysis. Has characteristics. In this case, to increase the number of digits after the decimal point in the temperature data,Temperature data is collected at a sampling time shorter than the calculation time step in the inverse problem analysis, and is averaged to obtain a representative value of the temperature data used in the calculation time step in the inverse problem analysis.
[0010]
Further, in the present invention, by performing an inverse problem analysis using an unsteady heat conduction equation based on temperature data measured at a temperature measurement point arranged inside a reaction vessel wall with a temperature change reaction, When determining the temperature or heat flux at the surface of the reaction vessel, increase the number of digits after the decimal point of the temperature data measured at the temperature measurement point, filter the temperature data with the increased number of digits after the decimal point, It is characterized in that the temperature data after filtering is used for the inverse problem analysis. In this case, for example, low-pass filter processing is used as the filter processing.
[0011]
DETAILED DESCRIPTION OF THE INVENTION
DESCRIPTION OF EMBODIMENTS Hereinafter, preferred embodiments of an inverse problem analysis method, apparatus, computer program, and computer-readable storage medium according to the present invention will be described with reference to the drawings.
[0012]
(First embodiment)
FIG. 1 shows the configuration of the inverse problem analysis apparatus in the present embodiment. In the figure, reference numeral 101 denotes an input unit for inputting temperature data measured by a thermocouple (see FIG. 2) embedded in the reaction vessel wall.
[0013]
Reference numeral 102 denotes a calculation unit, which performs a calculation to increase the number of digits after the decimal point using the temperature data input to the input unit 101 as will be described later.
[0014]
Reference numeral 103 denotes an inverse problem analysis unit, which performs an inverse problem analysis using an unsteady heat conduction equation based on the temperature data obtained by increasing the number of digits after the decimal point by the calculation unit 102, whereby the temperature on the surface of the reaction vessel or Obtain heat flux.
[0015]
An output unit 104 displays, for example, a change in temperature or heat flux on the surface of the reaction vessel calculated by the inverse problem analysis unit 103 on a display (not shown).
[0016]
First, the inverse problem analysis used in this embodiment will be described. As described in Patent Document 1 described above, the heat conduction phenomenon inside the reaction vessel wall is considered as an unsteady one-dimensional inverse heat conduction problem, and a single thermocouple temperature change or a plurality of one-dimensional directions are arranged. A method has been proposed for estimating the heat flux change on the inner surface of the reaction vessel from the thermocouple temperature change.
[0017]
FIG. 2 shows a two-dimensional cross section near the furnace wall of a reaction vessel (blast furnace) in which a plurality of thermocouples “x” are embedded. Although the boundary is indicated by a broken line in the furnace wall, one-dimensional means that only the heat flow in the direction along the broken line is considered. That is, for example, when heat conduction in the directions of 1a → 1b → 1c and 1d → 1e is assumed, the heat flux on the furnace inner surface is estimated. At this time, assuming that the cooling condition of the outer surface of the furnace is known, the heat flux on the unknown inner surface of the furnace is generally obtained. Of course, it is also possible to reverse the known and unknown boundary conditions.
[0018]
However, the original unsteady one-dimensional inverse heat conduction problem is to simultaneously estimate the boundary conditions on the inner and outer surfaces of the furnace, and in the inverse solution method assuming that one of the boundary conditions is known, We can only get an approximate answer to the boundary condition. For example, a reaction that is caused by a temperature fluctuation measured by a certain thermocouple due to a change in the heat flux on the inner surface of the reaction vessel as described above or due to poor contact of a cooling device installed outside the reaction vessel. It cannot be distinguished whether it is due to a heat flux change on the outer surface of the container.
[0019]
In order to evaluate more strictly, the heat conduction phenomenon should also occur in the upward direction across the broken line shown in FIG. 2, and it is necessary to solve the inverse heat conduction problem in two dimensions. In this case, even when the upper and lower boundaries in FIG. 2 are assumed to be adiabatic, it is necessary to construct a two-dimensional inverse problem for estimating a fine heat flux distribution at the left and right boundaries.
[0020]
Therefore, the applicant of the present application has proposed inverse problem analysis for enabling simultaneous estimation of changes in heat flux and temperature on the inner and outer surfaces of the reaction vessel. Hereinafter, the inverse problem analysis proposed by the present applicant will be described in detail. The unsteady heat conduction equation used for the inverse problem analysis is expressed as shown in the following Equation 1.
[0021]
[Expression 1]

[0022]
In Equation 1, ρ is the density of the material of the reaction vessel, C_pIs the specific heat of the material in the reaction vessel, T is the calculated value of the temperature inside the reaction vessel, t is the time, and k is the thermal conductivity of the material in the reaction vessel.
[0023]
Inverse heat conduction problem analysis refers to estimating boundary conditions such as temperature and heat flux at the boundary of the region based on the unsteady heat conduction equation that governs the calculation region, assuming the temperature inside the region as known. . On the other hand, the heat conduction order problem analysis is to estimate the temperature inside the region from the known boundary conditions such as the temperature at the region boundary and the heat flux.
[0024]
As an example of the method of the two-dimensional inverse problem analysis, for example, the one disclosed by the applicant of the present application in Japanese Patent Application Laid-Open No. 2002-206958 can be applied to the one-dimensional inverse problem analysis as it is. As an example of one-dimensional inverse problem analysis, an analysis method proposed by Beck et al. Is known (Beck et al., Inverse Heat Conduction, 1985, Wiley, New York).
[0025]
In addition, as a recent method of inverse problem analysis, it is conceivable to apply a stochastic estimation method such as Kalman filter theory or projection filter theory. At present, this method has been studied for application to the steady-state heat conduction equation (observation equation) with the left side of Equation 1 set to zero, but if the observation matrix can be appropriately configured including the unsteady term, A similar inverse problem analysis may be possible. An example of the application of probability estimation to this stationary differential equation is detailed in Tosaka et al., “Mathematical and Inverse Problems and Inverse Analysis of Partial Differential Equations” (The University of Tokyo Press (1999)).
[0026]
In the present embodiment, the idea disclosed in the above Japanese Patent Laid-Open No. 2002-206958 is used as the inverse problem analysis method. That is, as shown in the following equation 2, the temperature Y measured by each thermocouple arranged in a certain one-dimensional direction (1a → 1b → 1c, 1d → 1e, etc. shown in FIG. 2), and the inside of the reaction vessel The assumed values that minimize the sum of the squares of the differences from the temperature T at each thermocouple position calculated from the assumed values of the heat fluxes on the surface and the outer surface by the unsteady heat conduction equation are the inner and outer surfaces of the reaction vessel. Obtained as heat flux at J represents the number of thermocouples.
[0027]
[Expression 2]

[0028]
In this way, instead of obtaining a solution (heat flux on the inner and outer surfaces of the reaction vessel) that perfectly matches the temperatures T and Y at a plurality of thermocouple positions, a solution that satisfies the least squares is obtained. By obtaining this, it is possible to estimate a realistic heat flux change. The reason is that since the measurement temperature includes various measurement error factors, it can be said that it is practically meaningless to make them completely match.
[0029]
In order to stabilize the calculation, a regularization term may be added. Equation 3 below shows an example of a zero-order regularization term. p is the number of divisions of the estimated heat flux and α₀Is a regularization parameter obtained from experience values.
[0030]
[Equation 3]

[0031]
More specifically, an example of a formulation and calculation procedure for estimating the heat flux on the inner surface and the outer surface of the reaction vessel with the temperatures Y at a plurality of thermocouple positions as known will be described below.
[0032]
Number 4 S below_mRepresents the overall objective function, and Equation 5 below represents the objective function representing the deviation between the measured temperature Y and the calculated temperature T. The following Equation 6 is an objective function added to stabilize the calculation, and has a function of suppressing a rapid change in the value in the space division direction. Α in Equation 6₀And α₁Is a regularization parameter obtained from a certain empirical value.
[0033]
[Expression 4]

[0034]
[Equation 5]

[0035]
[Formula 6]

[0036]
In the above formula 5, the objective function is set so that the square of the difference between the temperature Y measured by a certain thermocouple and the temperature T calculated by the heat conduction equation model from the assumed value of the heat flux is minimized. In the above equation 6, an objective function for regularizing the spatial direction is set so that the solution is stable even if there is a temperature measurement error. Then, using Equation 4 as an overall objective function, a minimum point is searched for an unknown heat flux division region as shown in Equation 7 below.
[0037]
[Expression 7]

[0038]
Here, as shown in Equation 8, for the purpose of stabilizing the solution, it is assumed that the heat flux value at each time step remains unchanged until a certain future time. The time step varies depending on the thermophysical properties and shape of the target material. In Equation 8, q represents the heat flux, and the heat flux q in the m time step._mFrom the heat flux q in the future time m + r-1 time step._{m + r-1}Is assumed to be constant.
[0039]
[Equation 8]

[0040]
Then, if the minimization of Expression 7 is expanded using the assumption of Expression 8, it can be expanded into a matrix form as shown in Expression 9.
[0041]
[Equation 9]

[0042]
X of number 9^TX is derived from the first term on the right side of Equation 4, and X^T2 terms following X (α₀H₀ ^TH₀+ Α₁H₁ ^TH₁) Is derived from the second term on the right side of Equation 4 (the superscript T represents a transposed matrix). The configuration of X is as follows:_{j, i, k}As shown. Here, i indicating the number of divisions in the time direction changes up to maximum M time steps, j indicating the number of thermocouples changes up to J, and k indicating the number of divisions of the heat flux distribution is maximum. It changes to p. Note that the superscript * in Equation 9 indicates a reference value in the repeated convergence calculation, and T^*Is the temperature reference value, q^*Is the heat flux reference value. In the case of the one-dimensional case, the boundary condition at both ends is estimated, so the division number k of the heat flux distribution is p = 2 at the maximum.
[0043]
[Expression 10]

[0044]
Equation 9 is a simultaneous equation for estimating a change in heat flux when a temperature change occurs. In each time step, Equation 9 is used to obtain the heat flux q at both ends. First, the calculated temperature at the thermocouple position in the previous time step is set to the initial T^*And q is obtained by Equation (9). This q is given as a boundary condition of the forward problem heat conduction equation model calculated in parallel, and the temperature distribution is calculated. The calculated temperature value obtained here is the following temperature reference value T^*Then, q is corrected again (substituting into Equation 9 to obtain q again). This operation is repeated until q 5 becomes equal to or less than a certain residual (convergence).^*Are repeated, and the heat flux (final q) at both ends in each time step is obtained. By repeating this calculation procedure, two changes in the heat flux q at both ends can be estimated simultaneously.
[0045]
Equation 10 represents a kind of sensitivity matrix. In short, it represents the ratio of the magnitude of the change in the calculated temperature T at the thermocouple position to the unit change in the heat flux q at the boundary end point. Equation 10 can calculate the value per unit time step at each time step by the forward problem calculation calculated simultaneously with the inverse analysis.
[0046]
Hereinafter, a more preferable solution will be described by taking a one-dimensional inverse problem analysis as an example. As described above, even if a one-dimensional inverse problem is formulated (formulated) using the heat fluxes of the two end faces (inner surface and outer surface of the reaction vessel) as unknown boundary conditions, a solution can be obtained in principle. it can.
[0047]
However, there may be many solutions depending on the number of thermocouples, the thermophysical condition of the material, etc., and the calculation may become unstable. One reason for this is that if a combination of “difference in heat flux at unknown end faces” can be selected appropriately, there may be an infinite number of heat flux combinations that represent temperature changes at discrete temperature measurement points. Because there is. In particular, in the case of a substance with low thermal conductivity, even if the boundary condition that the surface temperature becomes extremely large or small is estimated, if only the temperature change at the discrete measurement point is reproduced. It can happen that it is recognized as one solution. This is not only impossible as a real phenomenon, but also makes the inverse problem calculation very unstable.
[0048]
Moreover, as an actual problem, the temperature of the thermocouple when starting the inverse problem analysis (the temperature of the discrete measurement points) is given as known, but the initial conditions of the temperature distribution in other analysis areas may be unknown. It is common. For this reason, the calculation is started from the provisional initial temperature distribution given arbitrarily, and while proceeding with the calculation step, the actual temperature distribution is searched and estimated, and the calculation is stably performed while gradually correcting to an appropriate temperature distribution. It is required to have a calculation logic that can be advanced (the temperature distribution here is calculated in parallel in order to correct the solution of Equation 9 in the calculation procedure of the inverse problem analysis, for example) Is the calculated value of the forward problem heat conduction equation model). Thus, the fact that the initial temperature distribution is indefinite is one of the major factors that make the inverse problem calculation unstable.
[0049]
The above indicates that in order to stabilize the inverse problem, it is necessary to provide a standard (constraint condition) of a certain surface temperature in the process of inverse problem analysis. Based on this concept, a method for appropriately giving the constraint condition will be described with reference to the flowchart of FIG.
[0050]
First, a temporary heat flux q is given as a heat flux on one side of the inner surface and the outer surface of the reaction vessel, here the outer surface. As a method of giving this temporary heat flux q, heat transfer coefficient h and reference temperature T_bAnd
q = h (T_surf-T_b)
(Step S201).
[0051]
T_surfIndicates the temperature at the unknown boundary, here the outer surface of the reaction vessel. This surface temperature T_surfIn order to correct the heat flux value during the inverse problem analysis, the forward problem analysis is usually performed at the same time, which corresponds to the surface temperature obtained by the forward problem analysis.
[0052]
Reference temperature T_bIs the temperature outside the inside and outside surfaces of the reaction vessel. In the present embodiment, the cooling condition of the reaction vessel is determined based on the water temperature or the like for water cooling, for example.
[0053]
As a result, the heat flux q, which is the left side of the above equation, can be given as known heat flux information. Thus, by giving temporary heat flux information, the heat transfer coefficient h and the reference temperature T_bTherefore, it is possible to secure physical validity as compared with the case where an arbitrary heat flux is applied, and to prevent an extreme temperature distribution from occurring.
[0054]
Next, the provisional heat flux q (= h (T_surf-T_b)), The heat flux on the inner surface of the reaction vessel that minimizes the sum of the squares of the differences between the temperatures T and Y shown in Equation 2 or Equation 5 above is expressed as the temporary heat flow on the inner surface of the reaction vessel. It is calculated as a bundle (step S202). This step is the main calculation procedure of the inverse problem analysis, and the formulation and calculation procedure shown in Equations 4 to 9 can be applied as they are as one specific solution. In this case, when solving Equation 9, the provisional heat flux q (= h (T_surf-T_b)) Is given as known and means that the temporary heat flux on the inner surface of the reaction vessel is solved as unknown.
[0055]
Here, the heat flux information on one side (outer surface of the reaction vessel) is given as described above, and the heat flux on the opposite side (inner surface of the reaction vessel) obtained by inverse problem analysis indicates the possibility of one solution. It only shows. Further, the heat transfer coefficient h and the reference temperature T assumed to be known._bIs also an approximate value, which is originally unknown.
[0056]
Therefore, heat transfer coefficient h and external reference temperature T_bAre changed by several points, that is, the temporary heat flux q on the outer surface of the reaction vessel is shaken by several points (K points), and the temporary heat flux q on the outer surface of the reaction vessel is changed, When each provisional heat flux information q is given, K combinations of heat fluxes on the inner surface of the reaction vessel that minimize the sum of the squares of the differences between the temperatures T and Y are obtained (step S203).
[0057]
Then, as shown in Equation 11 below, K temporary heat fluxes q on the outer surface of the reaction vessel and heat fluxes on the inner surface of the reaction vessel obtained corresponding to each temporary heat flux information q Of these combinations, the combination with the smallest square value of the difference between the temperatures T and Y is selected, and the combination is used as the heat flux on the inner and outer surfaces of the reaction vessel (step S204).
[0058]
## EQU11 ##

[0059]
The brackets in the above formula show the calculation results of one case where the inverse heat problem analysis was performed with the heat flux on one side known, and it means that the result with the smallest least squares difference is selected from the K case calculations. To do.
[0060]
By repeating this procedure at each time step, it is possible to sequentially and simultaneously calculate the heat flux change over time on the inner and outer surfaces of the reaction vessel.
[0061]
As described above, it is possible to stably execute a one-dimensional inverse problem analysis that simultaneously determines the heat flux change on the inner surface and the outer surface of the reaction vessel. If the temperature change and heat flux change on the inner and outer surfaces of the reaction vessel can be estimated simultaneously, for example, whether the temperature fluctuation at a certain temperature measurement point is due to the heat flux change on the inner surface of the reaction vessel. It is possible to distinguish whether it is due to a change in the heat flux on the outer surface of the reaction vessel, which is caused by a poor contact of the cooling device installed in the reactor.
[0062]
The above method is simple when applied to one-dimensional inverse problem analysis, and is often effective as an actual problem. The reason is that, generally, the upper and lower ends of the reaction vessel are often adiabatic (symmetric), and there is no practical problem.
[0063]
Therefore, an inverse problem analysis is performed assuming a one-dimensional thickness direction in a range delimited by a broken line in FIG. 2, and the result can be combined in the vertical direction to be two-dimensional.
[0064]
To more strictly consider the heat flow in the vertical direction in FIG. 2, a two-dimensional inverse problem analysis is necessary. Such a two-dimensional analysis is equivalent to finely dividing the heat flux at the left and right ends of FIG. 2 upward and obtaining a heat flux distribution that minimizes the temperature at these thermocouple positions in a least-square manner. In other words, the present invention may be applied according to a method similar to the inverse problem formulation disclosed in the above-mentioned JP-A-2002-206958.
[0065]
In this case, the heat flux at the upper and lower ends of FIG. 2 may be unknown or known, but considering the stability of calculation, an appropriate heat flux (for example, heat insulation) is determined from physical considerations. It is better to give it and make it known.
[0066]
Based on the same idea, the extension to the three-dimensional analysis can be easily performed.
[0067]
Here, for example, in order to estimate the “current” change of the heat flux change on the surface of the reaction vessel by the inverse problem analysis described above, it is necessary to estimate the past situation in the reaction vessel as close as possible.
[0068]
To that end, it is conceivable to shorten the unit time step of the inverse problem calculation. If the unit time step of the inverse problem calculation is long, a change going back to the past is estimated at least by the time step. Further, since the inverse problem calculation is a calculation result obtained by time averaging at the calculation time step, the restored heat flux change also captures a dull change, and a sudden change below the time step cannot be captured.
[0069]
However, shortening the unit time step of the inverse problem calculation means that during that short time, a small temperature transmitted from the heat flux change position (for example, the inner surface of the blast furnace in contact with the hot metal) to the thermocouple position. It means that change must be captured. In particular, when the thermocouple position and the heat flux change position are separated from each other and the thermal conductivity of the substance therebetween is small, the temperature movement in a short time becomes very small.
[0070]
In order to theoretically extract the temperature change due to the heat flux change from such a small temperature change and restore the heat flux change (short-time inverse problem analysis), the temperature measurement method and the temperature data processing method Ingenuity is necessary.
[0071]
Therefore, the inventors of the present application have made extensive studies to make it possible to estimate the past heat flux distribution as close as possible in a short time step while utilizing temperature data measured by a thermocouple.
[0072]
In the case of operation of a blast furnace or the like, only the change in temperature is observed, and when looking at the general temperature flow, it is common to evaluate the temperature data after the decimal point by rounding down.
[0073]
Even if you look closely at the current temperature change over time, the manufacturer-guaranteed range of temperature measurement devices that convert thermocouple contact voltage changes to temperature may be up to one decimal place in most cases. Often use temperature data. That is, the temperature data collected and recorded is coarse, such as 0 or 1 decimal point.
[0074]
However, in order to estimate the boundary condition change (heat flux change) at a distant position, which causes a very subtle temperature change (inverse problem analysis), the behavior of the temperature change after the decimal point is estimated. It becomes very important.
[0075]
That is, in the inverse problem analysis, the time transition of temperature and the flow of temperature change are more important than the accuracy of temperature itself. Pursuing only the accuracy of the temperature itself, if the temperature changes stepwise, the inverse problem analysis becomes unstable. In particular, when an inverse problem analysis is performed at a short time step, a step-like temperature change that cannot be considered at that time step may occur. For this reason, even if a heat flux solution is obtained by performing an inverse problem analysis with this stepwise temperature change, there is a possibility of reaching a physically impossible solution. In fact, in most cases, the solution diverges, making it impossible to continue the inverse problem analysis.
[0076]
Generating temperature data that artificially increases the number of decimal places from temperature data with a small number of decimal places in a reasonable way is an effective method for smoothing the step-like time change. It is.
[0077]
As one of the methods, the calculation unit 102 collects temperature data of 0 or 1 digit after the decimal point with a sampling time shorter than the time step in the inverse problem analysis, averages them, and calculates in the inverse problem analysis. This is the representative value of the temperature data used in the time step. For example, the temperature of a small number of digits sampled at intervals of 5 minutes shorter than the time step (1 hour) in the inverse problem analysis is simply averaged over a 1 hour range to obtain a representative value of the temperature data of the 1 hour step. This makes it possible to increase the number of decimal points of the temperature data for one hour. Using this method, temperature data used for inverse problem calculation based on temperature data (measured with a thermocouple) with zero digits after the decimal point and one digit after the decimal point for which the accuracy measured by the thermocouple is guaranteed It is possible to increase the number of digits after the decimal point.
[0078]
Further, the number of digits after the decimal point of the temperature data used for the inverse problem calculation may be increased by another method. For example, in a thermocouple, even if the guaranteed range of the temperature data itself is a single decimal point, it is actually calculated based on a conversion formula that converts the voltage at the thermocouple contact to temperature, or calibration was performed. A method of converting from a voltage-temperature conversion table to a temperature is common.
[0079]
In the former method, the calculated temperature value of a large number of digits is rounded off to one digit after the decimal point, and it is actually possible to obtain the number of decimal digits according to the ability of the computer. For example, as a conversion formula for converting from voltage to temperature in a K-type thermocouple, there is a formula shown in the following formula.
T = a + b · X + c · X²+ D · X^Three+ E · X^Four+ F · X^Five+ G · X⁶+ H · X⁷+ I · X⁸
T: Temperature [° C]
X: Contact voltage [V]
a = 0.2526584602
b = 24152.10900
c = 672333.4248
d = 2210340.682
e = −86096339.9
f = 4.83506 × 10^Ten
g = −1.18452 × 10¹²
h = 1.38690 × 10¹³
i = −6.333708 × 10¹³
By directly converting the contact voltage to the temperature using the conversion formula as described above, it is possible to easily obtain a numerical value after the decimal point of the temperature data. In addition, it cannot be overemphasized that the contact voltage here is the true contact voltage which passed through cold junction temperature compensation.
[0080]
(Second Embodiment)
In addition to increasing the number of digits after the decimal point as described in the first embodiment, a filter 101a is provided in the input unit 101 as shown in FIG. Can be mentioned. That is, the temperature data measured by the thermocouple is corrected using the filter 101a, and the corrected temperature data is used for the inverse problem analysis. By using a filter, even if the original temperature data has few digits after the decimal point, it is possible to generate digits after the decimal point in the temperature data. The problem analysis calculation is very stable.
[0081]
The theoretically most effective method is to use a low-pass filter. The low-pass filter corrects the actual measurement data in order to remove the influence of the fluctuation caused by the high frequency noise mixed in the actual measurement data sampled at equal time intervals.
[0082]
The reason why it is effective to correct the temperature data using a low-pass filter in the inverse problem analysis is that the change in the heat flux at a position far away from the thermocouple is slowed down while being transmitted through the refractory by the thermal diffusion phenomenon. This is because when the temperature is transmitted to the thermocouple position, the temperature change is considered to be very low (low frequency change is transmitted).
[0083]
This means that the temperature change at the thermocouple position caused by a change in the furnace interior far away is a low frequency change in time. However, high-frequency noise due to other causes is superimposed on the actual contact voltage change (original signal before temperature conversion) of the thermocouple. Therefore, using temperature data that has been processed to remove high-frequency noise from actual measured temperature changes in an inverse problem analysis is a physically very meaningful technique. Therefore, this method can not only improve the calculation stability of the inverse problem analysis but also improve the estimation accuracy of the heat flux.
[0084]
Low-pass filter processing is not only temperature data with the number of digits increased by using the time average mentioned above or using a conversion formula, but also temperature data with the number of decimal places increased by some method, Even if it is applied to measured temperature data with a small number of digits after the decimal point, it is very effective. That is, increasing the number of digits after the decimal point described in the first embodiment and performing low-pass filter processing may be used in combination, thereby synergistically smoothing the temperature change.
[0085]
Various types of filters can be considered. For example, R.W.Hamming's Digital filters, Dover publications Inc. 1998 may be referred to. There are many other reference materials.
[0086]
As a classic method of the low-pass filter, a method of correcting the temperature at the intermediate point using measured temperature data Y sampled at equal time intervals is known. As correction equations for the correction, the following equations (1) to (3) are derived. The subscript i indicates the number of times of sampling, and means that the current time i is corrected by the measured temperatures of the preceding and succeeding times.
[0087]
[Expression 12]

[0088]
Equations (1) to (3) show the equations of the 3-term method, 5-term method, and 7-term method, respectively, and the performance as a low-pass filter improves as the number of measured temperature data used for correction increases. That is, only lower frequency signals can be extracted. Basically, since it can be expressed by the coefficient of each term, when expressed in the Windows display, equations (4) to (6) are obtained, respectively.
[0089]
Generally, as the number of terms increases, the performance as a low-pass filter improves and the degree of freedom in filter design increases. That is, the polynomial filter can remove high-frequency disturbances over a wide range, so that the calculation of the inverse problem analysis is extremely stable. On the other hand, when the number of terms increases, in order to obtain a corrected temperature value at a certain time point, the measured temperature at multiple time points before and after that must be used, which is substantially retroactive. Therefore, it means that the “current” in the case of “estimating the current heat flux by inverse problem analysis” is pulled back by the low-pass filter. Therefore, there is a trade-off between ensuring the stability of inverse problem analysis using a low-pass filter and bringing the "current" estimation time as close as possible to the present, so it is necessary to select an appropriate combination in consideration of the relationship between the two. is there.
[0090]
From this point of view, the method described in the first embodiment uses the temperature data with the number of digits after the decimal point as much as possible to suppress a step-like temperature change and reduce the number of terms. Is effective in estimating the “current” heat flux.
[0091]
(Example)
FIG. 5 is a schematic example of an inverse problem analysis using a thermocouple embedded in a carbon brick installed at the bottom of the blast furnace bottom, assuming a one-dimensional unsteady heat conduction model. Indicate. Two thermocouples are embedded on the cooling surface side in a biased manner, TC1 indicates a high temperature side thermocouple, and TC2 indicates a low temperature side thermocouple. From the temperature changes measured by these thermocouples TC1 and TC2, the unsteady heat flux q1 on the high temperature heat flux surface and the unsteady heat flux q2 on the cooling surface are simultaneously estimated by the inverse problem analysis described above. The high-temperature heat flux surface position is a fixed point at a depth of 4.0 m from the cooling surface. The assumed thermal properties of the carbon brick are: specific heat Cp = 712 J / (kg · K), density ρ = 2300 kg / m^ThreeThermal conductivity k = 21.2 W / (m · K).
[0092]
FIG. 6A shows the estimation result of the high temperature side heat flux q1 obtained in this example. The horizontal axis in FIG. 6 represents the date, and shows the results from September 1 to October 1. Actually, the calculation is started before September 1 and the calculation ends on October 1st. In addition, the time step of the inverse problem analysis is an 8 hour step.
[0093]
Case 1 is a calculation result when a filter process is performed on data obtained by rounding off a numerical value of one decimal point of temperature data measured by a thermocouple to zero decimal places. Case 2 is a calculation result when the filtering process is executed on the data obtained by rounding the numerical value of the decimal point of the temperature data to one digit after the decimal point.
[0094]
FIGS. 6B and 6C are graphs plotting temperature data (temperature data used for inverse problem analysis) after filtering of the high temperature side thermocouple TC1 and the low temperature side thermocouple TC2. FIG. 6D shows the temperature difference between case 1 and case 2 at each position of thermocouples TC1 and TC2 (case 1 to case 2).
[0095]
As for the filter, the past temperature data is subjected to a low-pass filter with as many terms as possible to stabilize the calculation, and the Spencer equation (Equation (10)) with a maximum of 21 terms is used. A method of reducing the number is used. Specifically, the expressions (4) to (10) are used in a Windows format notation. In other words, the number of terms decreases as it approaches October 1, and finally equation (4) is used, and the last temperature data is not filtered.
[0096]
In both cases 1 and 2, the past analysis is very stable due to the filter effect, but as it approaches the present, it gradually becomes unstable and the final estimated heat flux is realistic. Is causing an unexpectedly large vibration. This is considered to indicate that the effect of the low-pass filter gradually decreases, and the calculation becomes unstable.
[0097]
Comparing Case 1 and Case 2, it can be seen that Case 2 is able to keep the magnitude of the final vibration smaller. This indicates that temperature data with a larger number of digits after the decimal point can perform relatively stable calculation even with a low-pass filter with a smaller number of terms.
[0098]
Furthermore, in any case, since vibration has occurred in the last estimation result, when comparing with temperature data of the same number of digits, increasing the number of terms of the low-pass filter and smoothing the temperature data as much as possible, It is also shown that this is an effective means for obtaining a stable solution. In particular, when a low-pass filter of the form of equations (4) to (10) is used, even if equation (4) is used, since the last data is not filtered, there are many values after the decimal point of the original temperature data. However, the inverse problem analysis solution is more stable because it is easier to smooth.
[0099]
In Case 1 and Case 2, the inverse problem analysis was attempted without applying any filter to the temperature data, but the solution diverged and the calculation could not proceed.
[0100]
FIG. 7 shows a schematic diagram similar to FIG. 5 for cases where the installation depth of the thermocouple is different. FIG. 8A shows the estimation result of the high temperature side heat flux q1 when the calculation time step of the inverse problem analysis is 8 hours and when it is 6 hours. These calculation results are shown by extracting from September 1 to October 1 in the calculation results of a longer period.
[0101]
The temperature data used for the calculation is one digit after the decimal point, and smoothing is performed using the equation (10) as a filter. FIGS. 8B and 8C are graphs plotting temperature data of the high temperature side thermocouple TC1 and the low temperature side thermocouple TC2 after the smoothing. FIG. 8D substantially shows a difference (TC1−TC2) between the high temperature side thermocouple temperature and the low temperature side thermocouple temperature corresponding to the steady heat flux.
[0102]
Comparing the transition of the quasi-steady heat flux in FIG. 8 (d) with the peaks of the unsteady heat flux estimated by the inverse problem analysis in FIG. 8 (a), it is intuitive that the inverse problem analysis is performed. It can be seen that the delay time is corrected and the peak appears earlier in the unsteady heat flux (indicated by the broken line in FIG. 8A).
[0103]
The quasi-steady heat flux in the graph of FIG. 8D can be interpreted as an average heat flux passing around the thermocouple position, and the value of the unsteady heat flux estimated on the high temperature surface at the depth of 4 m. In comparison, the peaks do not always correspond to each other, but in this case, a peak delay of approximately 2 to 3 days is observed. In other words, the unsteady heat flux by the inverse problem analysis detects the movement in the furnace approximately 2-3 days earlier than the conventional steady method.
[0104]
Further, when comparing the 8-hour step calculation and the 6-hour step calculation in the graph of FIG. 8A, the 6-hour step calculation shows a clearer movement. If the time step is lengthened, only the averaged motion within that time can be captured. Therefore, it is more likely that the actual internal motion can be expressed if the time step is shortened as much as possible.
[0105]
(Other embodiments)
The inverse problem analysis apparatus described above is composed of a CPU or MPU of a computer, RAM, ROM, RAM, and the like, and is realized by operating a program stored in the RAM, ROM, etc. as described above. .
[0106]
Therefore, the program itself realizes the functions of the above-described embodiment, and constitutes the present invention. As a program transmission medium, a communication medium (wired line or wireless line such as an optical fiber) in a computer network (LAN, WAN such as the Internet, wireless communication network, etc.) system for propagating and supplying program information as a carrier wave Can be used.
[0107]
Furthermore, means for supplying the above program to a computer, for example, a storage medium storing such a program constitutes the present invention. As such a storage medium, for example, a flexible disk, a hard disk, an optical disk, a magneto-optical disk, a CD-ROM, a magnetic tape, a nonvolatile memory card, a ROM, or the like can be used.
[0108]
It should be noted that the shapes and structures of the respective parts shown in the above embodiments are merely examples of implementation in carrying out the present invention, and these limit the technical scope of the present invention. It should not be interpreted. That is, the present invention can be implemented in various forms without departing from the spirit or the main features thereof. For example, in order to use the present invention in a network environment, some programs may be executed on another computer.
[0109]
【The invention's effect】
As described above, according to the present invention, the inverse problem analysis using the unsteady heat conduction equation is performed based on the temperature data measured at the temperature measurement point arranged inside the reaction vessel wall. When calculating the temperature or heat flux at the surface of the reaction vessel, for example, even if the time step in inverse problem analysis is shortened to estimate the past temperature or heat flux change as close as possible, the decimal point of the temperature data By increasing the number of digits or by applying a filter process, it becomes possible to capture a small temperature change and stabilize the inverse problem analysis.
[Brief description of the drawings]
FIG. 1 is a diagram illustrating a configuration of an inverse problem analysis apparatus according to a first embodiment.
FIG. 2 is a diagram showing a two-dimensional cross section near the furnace wall of a reaction vessel (blast furnace) in which a plurality of thermocouples are embedded.
FIG. 3 is a flowchart for explaining processing in inverse problem analysis;
FIG. 4 is a diagram illustrating a configuration of an inverse problem analysis apparatus according to a second embodiment.
FIG. 5 is a diagram illustrating an arrangement relationship of thermocouples TC1 and TC2 in the embodiment.
FIG. 6 is a diagram for explaining an analysis result in the example.
FIG. 7 is a diagram showing an arrangement relationship of thermocouples TC1 and TC2 in another embodiment.
FIG. 8 is a diagram for explaining an analysis result in another embodiment.
[Explanation of symbols]
101 Input section
101a filter
102 Calculation unit
103 Inverse problem analysis
104 Output section

Claims (10)

Based on the temperature data measured at the temperature measurement point located inside the reaction vessel wall with temperature change reaction, the inverse temperature analysis using the unsteady heat conduction equation is performed, so that the temperature at the surface of the reaction vessel is Or an inverse problem analysis method for obtaining heat flux,
The temperature measurement is performed by collecting temperature data with a sampling time shorter than the calculation time step in the inverse problem analysis, and averaging them to obtain a representative value of the temperature data used in the calculation time step in the inverse problem analysis. Has a procedure to increase the number of decimal places of temperature data measured at points,
An inverse problem analysis method characterized by using temperature data with an increased number of digits after the decimal point for the inverse problem analysis. In addition to the procedure to increase the number of decimal places the temperature data measured above SL temperature measurement point,
The method further includes a step of filtering the temperature data with the number of digits after the decimal point increased,
Increase the number of digits following the decimal point, and the inverse problem analysis method according to claim 1 which comprises using the temperature data after the filtering process to analyze the inverse problem. The inverse problem analysis method according to claim 2 , wherein the filter processing is low-pass filter processing. A plurality of temperature measurement points are arranged at least in the thickness direction inside the reaction vessel wall, and the temperature or heat flux at the inner and outer surfaces of the reaction vessel is obtained by the inverse problem analysis. Item 4. The inverse problem analysis method according to any one of Items 1 to 3 . The temperature measured at each temperature measurement point and the temperature at each temperature measurement point position calculated by the unsteady heat conduction equation from the assumed temperature or heat flux at the inner and outer surfaces of the reaction vessel. 5. The inverse problem analysis method according to claim 4 , wherein the assumed value that minimizes the sum of the squares of the differences is obtained as temperature or heat flux on the inner surface and outer surface of the reaction vessel. Based on the temperature data measured at the temperature measurement point located inside the reaction vessel wall with temperature change reaction, the inverse temperature analysis using the unsteady heat conduction equation is performed, so that the temperature at the surface of the reaction vessel is Or an inverse problem analysis device for obtaining heat flux,
The temperature measurement is performed by collecting temperature data with a sampling time shorter than the calculation time step in the inverse problem analysis, and averaging them to obtain a representative value of the temperature data used in the calculation time step in the inverse problem analysis. A means to increase the number of digits after the decimal point of temperature data measured at a point,
An inverse problem analysis apparatus characterized in that temperature data with an increased number of digits after the decimal point is used for the inverse problem analysis. In addition to the means for increasing the upper Symbol temperatures number of decimal places the temperature data measured at the measurement point, further comprising a filter means for performing filtering process on the temperature data with an increased number of places the decimal point,
Increase the number of digits following the decimal point, and, inverse analysis apparatus according to the temperature data after the filtering process to claim 6, characterized by using the analysis the inverse problem. Based on the temperature data measured at the temperature measurement point located inside the reaction vessel wall with temperature change reaction, the inverse temperature analysis using the unsteady heat conduction equation is performed, so that the temperature at the surface of the reaction vessel is Alternatively, a computer program for causing a computer to execute a process for obtaining heat flux,
The temperature measurement is performed by collecting temperature data with a sampling time shorter than the calculation time step in the inverse problem analysis, and averaging them to obtain a representative value of the temperature data used in the calculation time step in the inverse problem analysis. Execute the process to increase the number of digits after the decimal point of the temperature data measured at the point,
A computer program characterized by using temperature data with an increased number of digits after the decimal point for the inverse problem analysis. In addition to the above SL temperature measurement points increase the number of decimal places the temperature data measured in the process, to further execute processing performing filter processing on the temperature data with an increased number of places the decimal point,
Increase the number of digits following the decimal point, and a computer program according to the temperature data after the filtering process to claim 8, characterized by using the analysis the inverse problem. A computer-readable storage medium storing the computer program according to claim 8 or 9 .

Images