Regular Article
An Identification Method of Multi-point Heat Transfer Coefficient
2016 Volume 102 Issue 1 Pages 29-33
Details
2016 Volume 102 Issue 1 Pages 29-33
A heat transfer coefficient identification model has been developed which considers multiple boundary with different conditions using influence coefficient matrix derived from a mathematically and numerically obtained Jacobian matrix. The model was applied to a top surface and bottom surface heat transfer estimation problem, i.e. two unknown quantities problem, for an infinite flat steel plate which is subject to different boundary conditions. The method yields thus fast and accurate heat transfer coefficient results as well as the temperature distribution which is in excellent agreement with measured one. Further application of the model was two-dimensional problem of a work roll temperature and boundary condition estimation while rolling a hot strip. The two-dimensional temperature distribution of a cross section of the roll was calculated by using heat transfer coefficient as input values which are inversely identified by the already calculated temperature distribution. The identified heat transfer coefficient directly confronted with the input value exhibiting the method’s capacity for multiple-point identification of heat transfer coefficient.
鋼材の焼入れ過程において，表面の熱的境界条件である熱伝達係数は，材質を決定する重要な因子である冷却速度を左右するため，非常に重要な要素である1,2,3)。炭素鋼の焼入れを例に取ると，冷却速度が比較的速い場合には，マルテンサイトという硬質な相ができる一方，冷速が遅い場合には，フェライトやパーライトといった相ができる。マルテンサイト相は，通常非常に硬質であり，かつ，脆い性質を持っている一方で，フェライト相は軟質で，延性に富む。これらの相を目的の性質となるように混合するためには，冷却速度を高精度に制御する必要があり，熱処理シミュレーションにおいても重要な因子となっている4)。
このような相分率と同時に，鋼材の寸法精度や残留応力も大きな品質を決定づける要素である。近年，焼入れ後の寸法変化や残留応力を精度よくコンピュータシミュレーションによって予測する技術が開発されており5,6,7,8,9)，最適な焼入れ条件の探索に利用されている。
精度のよいシミュレーションを実現するためにも，熱伝達係数の正確なインプットが必須となるため，銀棒を用いた熱伝達係数の推定方法が開発されている10)。本方法では，銀の熱伝導率の高さを利用しており，銀棒内の温度分布が無視できるために集中熱容量法などの比較的簡易な方法で熱伝達率を同定することが可能となっている。しかし，水焼入れの場合には，焼入れ中の銀棒の表面と鋼の表面性状は著しく異なっており，実際の熱伝達率が正確に測定できないという問題がある。
そこで，実際の焼入れを行う部材と同じ材料を用いた熱伝達率の推定方法が考案されている。本方法では，無限遠平板を仮定した1次元熱伝導方程式を用いた逆問題的手法を用いることが一般的である11)。従って，表裏で異なる冷却を受けるため，無限遠平板の仮定が成り立たない場合や，複雑な形状を持つ部材の場合など，複数の熱伝達境界を同定対象とする場合には適用できなかった。
そこで本研究では，複数の熱伝達境界を持つ冷却について，その熱伝達係数を同時に，かつ，高速に同定する方法を示し，その適用例として，表裏面から異なる冷却を受けた場合の表裏熱伝達係数の同定結果を実験結果と共に示す。さらに，本方法を2次元に拡張し，圧延中のロールの温度変化のシミュレーション結果から，ロール周方向の複数点の熱伝達係数を同定し，シミュレーションに用いたインプット条件と比較し，その有効性を議論する。
本モデルでは，n個の温度測定点を持ちn個の熱伝達係数を同定する場合を想定している。ただし，nは1以上の整数である。
このとき，n=1，すなわち，同定対象となる熱伝達係数が1個の場合は，従来通りの方法で，熱伝導方程式から逆問題的に境界条件を算出することが可能である。つまり，未知数である熱伝達係数hを適宜変更しながら，熱伝達係数hを仮定したときの温度測定点における計算温度と，測定温度が一致するまで繰返し計算を行えばよい。
ところで，この同定対象(=測定点)が2点以上の場合には，同定に要する計算回数が指数関数的に増加するため，通常の方法では困難である。そこでまず，同定対象となる位置iにおける熱伝達係数と温度との関係を線形化する。なお，本論文では下付き添字の総和規約を採用する。
| (1) |
式(1)では，測定点iにおける温度変化に全ての位置における熱伝達係数の影響が考慮されている。ここで，以下ではヤコビアンAij=∂Ti/∂hjを影響係数マトリクスと呼ぶことにする。
初期条件として，n個の熱伝達係数の初期値hi(i=1, 2, …, n)を仮決定し，現在の時刻tにおける温度分布が既知であるとして，時間変化Δt後の温度分布を，数値計算などを用いて求める。
次に，i=1の位置における熱伝達係数をdh1変化させ(hi=hi+dhi)，温度分布を同様に計算する。以下，dhiを試行熱伝達率変化と呼ぶ。このとき，上記の熱伝達係数の初期値hiを用いた場合の温度分布との差dTiを算出する。この操作で，影響係数マトリクスの一部である，∂Ti/∂h1(i=1, 2, …, n)が得られる。この操作を測定点の個数nだけ繰り返すことで，全ての影響係数マトリクスが得られることになる。ただし，i≠1の場合には，前回の変更を元に戻す操作(hi−1=hi−1−dhi−1)を行うこととする。
ここで，熱伝達係数の初期値hiを用いた場合の計算結果と測定結果との差のベクトルΔTiを予め求めておく。このとき，影響係数マトリクスは既に得られているので，測定値に計算値を近づけるための熱伝達係数の修正量Δhiは，影響係数マトリクスの逆行列を用いて，以下のように表わすことができる。
| (2) |
このようにして修正を繰り返すことによって，高速に複数点の熱伝達係数を同定することが可能となる。本計算シーケンスをFig.1に示す。計算開始時の初期温度は，均一温度等を仮定し，ある時間ステップ後の温度を計算する。一方，次ステップの温度計算には，前ステップにおける温度計算結果が初期温度として利用されることとなる。
Flow of heat transfer identification procedure.
鉄鋼材料の製造過程において，加熱炉や圧延ラインにおける板材の板厚方向温度分布を求めることは，荷重予測や材質制御の観点から非常に重要である。例えば，スラブの加熱過程においては，スキッドの影響などによってスラブの表裏面の熱伝達係数が異なることが想定され，さらに圧延ラインにおいては，デスケーリングやテーブル搬送における接触熱伝導などによって，表裏面の熱伝達係数が変化することが想定される。
そこで，ここでは平板の板厚方向1次元問題として，板の表裏面は異なる熱的境界条件である場合を考える。
3・1 実験条件供試材は，板厚20 mm，板幅100 mm，板長200 mmのS45C(JIS)鋼を用いる。本鋼板にFig.2に示すように板幅方向・板長方向の中心部における表面から1.0 mm位置および裏面から2.5 mm位置にシース熱電対を溶接し，板厚方向2か所の温度を測定する。本鋼板を加熱炉で1100°Cまで加熱し，5 min保持する。温度保持後に，炉外へ抽出し，空冷を開始する。
Schematic illustration of the sample and thermo-couples.
本実験に用いる実験設備をFig.3に示す。空冷開始後60 sec後に，板を，水冷ノズルを配備したチャンバー内を往復させ，板の表面のみ水冷を行う。すなわち，往復の間に板の表面は2回水冷を受けることになる。このとき，裏面は空冷および漏れ水の影響のみの冷却となる。表面のみ水冷した後，引き続き表裏共空冷を行い，温度を測定する。
Schematic illustration and photograph of the experimental apparatus.
境界条件である熱伝達率の同定には，密度，熱伝導率および比熱の熱物性値が温度毎の関数として与える必要がある。S45Cにおけるこれらのデータは，材料データベースMATEQ12)のデータを用いる。
3・3 計算条件板厚方向の1次元問題を仮定し，両端が境界条件となる100節点(99要素)のFEM温度計算を行う。節点の配置は，両端付近が細かいメッシュとなるように等比級数で決定し，境界付近の大きな熱流束にも対応できるように設定する。このとき，等比級数の公比は1.01とする。
加熱炉から試験片を抽出した時間を冷却開始0 secと定め，このときの表面付近の熱電対の測定温度が板内で均一であると仮定した後，温度計算を開始する。ここで，各点における熱伝達率の初期値を仮定する必要があるが，この値は正解に近い方が収束が速くなる。そこで，計算開始時の熱伝達率の初期値は，表裏ともhi=10 W/(m2K)for ∀iとし，次ステップからの初期値は，前ステップにおける収束後の各点の熱伝達率を用いる。また，試行熱伝達率変化は，計算開始時はdhi=1000 W/(m2K)for ∀iであり，次ステップの初期値は前ステップにおける収束後の値を用いるのは上記と同様であるが，各iterationにおいても各点における計算温度と実測温度との差dT1およびその点の熱伝達率を変化させた前後の計算温度同士の差dT2を用いて，下記のようにダイナミックに変化させる。
| (3) |
ここで，係数の0.1は計算を安定化させる値であり，各ステップの時間変化や温度変化などによって最適値が変化するが，本研究では簡単のため全計算過程で同一の値を用いる。また，dhiの上限は1000 W/(m2K)とする。
3・4 実験結果および計算結果実験結果として得られた表裏面付近に着けた熱電対の測定温度の結果をFig.4に示す。また，本測定温度から，表裏の熱伝達係数を以上で述べた方法を用いて同定した結果をFig.5に示す。Fig.5(b)は，(a)の縦軸を拡大したものである。ここで，同定した熱伝達係数を用いて，熱電対測定位置の温度を再度計算した結果をFig.4に同時に示してある。Fig.4から明らかなように，2回の水冷によって表層のみ温度が降下しており，逆に裏面の温度は水冷によって大きく変化していない。ここで，実測の温度と計算温度は完全に一致しており，Fig.5に示す，推定した熱伝達係数が確からしい値であることが示される。さらに，2回の水冷によって，熱伝達率が2回大きく上昇しており，水冷による冷却効果が再現できていることが分かる。このとき，裏面の熱伝達係数は，水冷初期において漏れ水の影響と思われる軽微な上昇が認められるが，ほぼ空冷の状態と変わらない状態を示していることから，表裏の熱伝達係数が同時に，また精緻に導かれていることが分かる。
Comparison of measured and calculated temperatures at near top and bottom surface positions.
Identified heat transfer coefficient of top and bottom surface.
次に熱間圧延中のワークロール表層における熱伝達係数の同定について述べる。熱間圧延中のワークロールは，圧延材からの入熱と，空冷およびロール冷却装置による水冷を繰返し受けるため，一般的にワークロール表層の境界条件は非常に複雑で，これまで，水冷や空冷などのゾーン毎に熱伝達係数を設定し，実際のロールの膨張(サーマルクラウン)量と計算値が合うデータを用いてきた。ところが，本方法では，試行錯誤的に複数回の熱計算を行う必要があることや，ゾーンの分け方によっては，正しい境界条件が得られない可能性があった。そこで，多点熱伝達係数の同定法を適用し，圧延ロール表面の熱伝達係数の同定を試みる。ここでは，インプット条件として予め定めた熱伝達係数によってロール温度計算を行い，次に計算結果として得られたロールの温度分布から，表面の熱伝達係数を逆解析によって求め，インプット値と比較する。
4・1 圧延中のワークロール温度変化計算圧延中のワークロール温度変化は，ロール径80 mmにおける断面の2次元問題として取り扱う。2次元熱伝導方程式の離散化にはFEMを用いることとし，このときの断面のメッシュ分割をFig.6に示す。メッシュは，表層部の急激な温度変化が表現できるように，表層部に近づくに従い細かくなっている。ここで，総節点数は49245，総要素数は48220である。
Finite element division of work roll cross section.
対象とするロールは，直径80 mmの中実円であり，熱物性値は密度7850 Kg/m3，比熱640 J/(Kg K)，熱伝導率23 W/(m K)とする。熱伝達境界条件はFig.7に示すように，角度によって圧延，水冷，空冷を考慮し決定する。ただし，Twはロール表面が接触する物体の温度である。初期温度20°Cであるロールを角速度3.33 rad/sec(約8 mpm)で回転させた場合のロール温度を計算する。ただし，ワークロールの回転は反時計回りであるとする。計算時間は0.000225 sec/stepである。本条件で計算した結果として，ロール120周後のワークロール断面の温度分布をFig.8に示す。
Input heat transfer coefficient as a boundary condition for 2-dimensional temperature calculation.
Calculated temperature distribution of the work roll during hot rolling.
Fig.8から，圧延材との接触によってロール表層に大きな熱流束が発生し，急激に温度上昇を生じるが，ロールの回転と共に冷却されており，徐々に表層と内部との温度差が緩和している様子が分かる。この熱は，徐々に内部に浸透していき，ワークロールのサーマルクラウンに帰結する。
4・2 ロール表層の熱伝達係数推定次に，計算で得られた時系列温度分布のうち，同定に用いる温度データとして，周方向に均等に10分割した各周方向位置(Fig.7に示す0.314 radの位置から0.628 radずつ回転させた位置)における表層から1 mm中心部へ入った部分の計算温度を用い，表面の熱伝達率を逆算することを考える。すなわち，フローチャートFig.1におけるTjmeasは上記位置における4・1節での計算結果であり，同定計算によってTj0を順次計算し，各時間・各点におけるhiを求める。熱伝達率の同定計算は，ワークロールが初期温度20°Cの状態からスタートし，ロールが120周までの同定計算を上記同様に0.000225 sec/stepで行う。同定する熱伝達係数は周方向に10点であるため，熱伝達の影響係数マトリクスは10×10のマトリクスとなるが，その逆行列算出には，ガウスの消去法を用いる。初期値であるhiと試行熱伝達率変化dhiは3・3節に示す板の問題と同様の手法で決定するが，その値はhi=115 W/(m2K)for ∀i，dhi=100 W/(m2K)for ∀iとする。
ここで，圧延によって入熱される熱量は極めて大きく，かつ，瞬間的であるため，圧延による入熱も同時に同定する場合，Tjmeasとして使用する温度はより表層付近である必要があり，実用上問題がある。このため，本研究では，圧延によって入熱されるロール周方向位置および当該位置における熱伝達率200000 W/(m2K)は既知であるとして，残りのロール周方向位置の水冷や空冷による熱伝達率を同定する。
このとき，入力として与えたワークロール120周後の熱伝達係数と逆問題手法で得られた熱伝達係数との比較をFig.9に示す。Fig.9から，圧延材による加熱以外の部分である空冷や水冷部においては，インプットした熱伝達係数とほぼ同様の値が同定できており，2次元問題において，多点の熱伝達係数を同時に同定可能であるということが示された。
Comparison of heat transfer coefficient between input value and inversely estimated value.
鉄鋼材料の熱処理シミュレーションにおいて重要な境界条件である熱伝達係数について，複数点に及ぶ熱伝達係数を同時に同定する方法として，各点における熱伝達の影響係数マトリクスを作成し，その逆行列を用いる方法を示した。
次に，本方法を表裏異なる平板の冷却実験に適用し，表裏の熱伝達係数を同定した。当該境界条件によって計算された温度と測定温度を比較することで，精度よく熱伝達係数が同定されていることを確認した。
さらに，本方法を2次元問題へ拡張した。圧延中のワークロールの表層における入熱と冷却の境界条件を予めインプットし，ロール断面の2次元温度計算を行い，表層付近の温度を算出した。当該計算温度から，ロール1周を10分割したデータとして熱伝達係数を逆問題的に推定し，インプットした熱伝達率とほぼ同等の値となることを確認した。
