ラプラス変換を用いた熱伝導逆問題解析手法の鋼材冷却への適用

Abstract

When we evaluate complicated cooling heat transfer characteristics of high-temperature steel with water, a solution of the inverse heat conduction problem (IHCP) is essential to specify the surface heat flux and surface temperature change with time from the measurement of the inside solid temperatures. In the present study, the analytical solution of the IHCP with the Laplace transform technique proposed by Monde is applied to steel cooling processes. This technique deals uniform and constant thermal properties of the solid, and the measured temperature change at two depths is approximated with half-power polynomials of time. The steel cooling processes show a very high cooling rate and a large temperature drop after quenching. Therefore, the temperature dependence of thermophysical properties could not be negligible, and much better accuracy of the approximate functions is also required for the complicated measured temperature history. We applied modifications to the analysis, such as a simple and explicit treatment of the temperature-depending properties and superimposing approximate functions. The accuracy of the modified technique was assessed for quenching temperature histories in a steel specimen. From the sensibility of approximation parameters constructing the superimposed functions, the optimum parameters were specified for the measured temperature of quenching experiments. The explicit treatment of the thermal properties ensures accurate estimation results for a higher cooling rate, in which the product of the temperature gradient near the surface and the temperature gradient of the thermal conductivity is small.

1. 緒言

鋼材を冷却する連続鋳造や熱処理設備ではラミナー,ジェット,スプレーなどの水冷却設備がコストと冷却能力の高さから広く用いられる¹⁾。鋼材温度が高いことから水冷時,膜沸騰から遷移沸騰域を経て核沸騰に至る複雑な非定常沸騰伝熱特性が支配する。この伝熱特性を解明するには,例えばビオ数:Biが0.06未満では集中熱定数系近似で計算されるが²⁾,ビオ数の大きい鋼材の水冷では内部温度を一様と見做せないので固体内部測定点の温度変化から固体表面上の熱流束を逆算する必要がある。

熱伝導逆算法としては例えばBeck³⁾が提案した推定法に代表される探索的な手法がある^3–6)。これらは,未知の表面熱流束を変化させながら,固体内部測定点の温度推定誤差を表す評価関数(残差)が最小になるように反復計算する数値解法である。本研究で提案する逆算法は,Mondeら⁷⁾が提案したラプラス変換法を用いた解析的解法に基づく厳密解法のため,反復計算なしに温度・熱流束が数式で得られ,解の安定性の点でも優れている。一方,本手法は解析解が得られる定物性の平板や円筒などの単純形状に対象が限定される。また,クエンチを伴う不規則な“実測”の温度変化に対して測定点の温度変化を“数式”で与える必要がある。これに対し従来研究では,(1)小時間区分で分割近似する改良版時間区分法⁸⁾にて,(2)温度変化を低次の時間の半値多項式(t^1/2の多項式)で表し^7–9)高い精度の近似式を解析に適用している。

ところで鉄鋼材料は扱う温度と変動幅が大きく熱物性の温度依存性が顕在化するため,熱伝導方程式の熱物性の非線形項の影響が顕著となる。本手法は比較的小さい冷却速度の高温鋼材の冷却の解析にも用いられているが^10,11),解析的解法を用いた従来研究^7–9)では熱物性の温度依存項を無視しているため,高温域の強冷却時の伝熱評価に適用するに当たり,温度依存性が推定解に及ぼす影響とクエンチ温度履歴に対する近似精度の検証が課題となっていた。ただし,実際の鋼材冷却では境界条件となる沸騰伝熱特性が未知なため,実験データによる逆算結果の精度検証は一般に困難となる。

そこで,本報では定物性或いは変物性の非定常熱伝導順問題解析の境界条件として,想定した3つの鋼材冷却の沸騰曲線に対する鋼材内部温度履歴をそれぞれ数値計算で求めた。実測値の代わりに既知の境界条件に対する内部温度履歴から推定した表面温度,表面熱流束,沸騰曲線に及ぼす種々の解析パラメータの影響を評価し,鋼材のクエンチ実験に特化したパラメータの最適値を定めた。さらに,定物性の逆解析解で熱物性の温度依存性を簡易的に取り扱う手法を提案し,変物性条件での表面熱流束の推定誤差が10%以内に収まる計測条件を定めて実際のクエンチ実験評価での便宜を図ることを目的とする。

2. 逆問題の定式化

2・1　定熱物性熱伝導方程式の解法

Fig.1に示す熱物性が均質一定な平板内部のx₁,x₂の2点の測定温度履歴から表面x=0の境界条件を推定する汎用的な一次元非定常熱伝導逆問題解について,ラプラス変換を用いて定式化する。熱伝導方程式(1a)を一様温度T₀の初期条件式(1b)と2点の内部温度履歴の境界条件式(1c),(1d)の下で解く。式(1a)の両辺にラプラス変換を施すと2階常微分方程式(2)が得られる。本稿では原関数にチルダ(~)を付加してラプラス変換後の像関数を表記する。像関数は,ラプラス演算子sと空間座標xを独立変数とする複素関数である。式(2)の一般解は式(3)で与えられる。

  

$\frac{\partial T}{\partial t}=a\frac{{\partial }^{2}T}{\partial x^2};\left(0\le x\le L_x\right)$(1a)

  

$$t=0;T=T_0;\left(0\le x\le L_x\right)$$(1b)

  

$$x=x_1;T=T_1\left(t\right)$$(1c)

  

$$x=x_2;T=T_2\left(t\right)$$(1d)

  

$\frac{{\partial }^2\tilde{T}}{\partial x^2}=\frac{s}{a}\tilde{T}\left(s,x\right)-\frac{1}{a}{T}_0$(2)

  

$\tilde{T}=\frac{1}{s}{T}_0+A{\text{e}}^{x\sqrt{s/a}}+B{\text{e}}^{-x\sqrt{s/a}}$(3)

Fig. 1.

Calculation domain of one-dimensional inverse heat conduction problem and temperature measuring positions in a finite solid.

式(1c),(1d)の各計測位置x₁, x₂の温度履歴T_k(t)は,初期温度0を基準温度とする温度履歴f_k(t)を用いて式(4)で表した後,ラプラス変換を施すと式(5)を得る。f_k(t)は測定温度履歴を近似する時刻の半値多項式(t^1/2の多項式)で,詳細は後述する。式(5)を式(3)に代入してA, Bを消去し,sinh x=(e^x+e^-x)/2を用いて整理すると,温度分布の逆問題解式(6)を得る。なお,熱伝導方程式と境界条件の線形性からT₀の値の相違は温度基準の相違だけで解の導出過程に影響しないので,以後の展開は,T₀=0として記述する。

  

$$T_k\left(t\right)=f_k\left(t\right)+T_0\left(k=1,2\right)$$(4)

  

${\tilde{T}}_k\left(s\right)={\tilde{f}}_k\left(s\right)+\frac{1}{s}{T}_0$(5)

  

$\tilde{T}\left(s,x\right)={\tilde{f}}_1\left(s\right)\frac{\text{sinh}\left\{\left(x_2-x\right)\sqrt{s/a}\right\}}{\text{sinh}\left\{\left(x_2-x_1\right)\sqrt{s/a}\right\}}-{\tilde{f}}_2\left(s\right)\frac{\text{sinh}\left\{\left(x_1-x\right)\sqrt{s/a}\right\}}{\text{sinh}\left\{\left(x_2-x_1\right)\sqrt{s/a}\right\}}$(6)

また,熱流束はFourier則に基づき式(6)をxで微分した温度勾配を適用すると式(7)で表すことができる。

  

$$\frac{\tilde{q}\left(s,x\right)}{\lambda }={\tilde{f}}_1\left(s\right)\frac{\sqrt{s/a}\text{ cosh}\left\{\left(x_2-x\right)\sqrt{s/a}\right\}}{\text{sinh}\left\{\left(x_2-x_1\right)\sqrt{s/a}\right\}}-{\tilde{f}}_2\left(s\right)\frac{\sqrt{s/a}\text{ cosh}\left\{\left(x_1-x\right)\sqrt{s/a}\right\}}{\text{sinh}\left\{\left(x_2-x_1\right)\sqrt{s/a}\right\}}$$(7)

複素平面上での厳密解式(6),(7)はf ~_kと双曲線関数の分数式(以後カーネル関数と呼ぶ)との積で与えられる。これらをラプラス逆変換して時間領域での温度T(t,x),熱流束q(t,x)の厳密解が得られるが,sのべき関数で与えられるf ̃_kとカーネル関数との積は複雑な関数形でラプラス変換表を用いた逆変換は不可能となる。Mondeら⁷⁾は厳密解を級数展開した後,項別に変換表を用いてラプラス逆変換を行い時間領域での厳密解を求める手法を提案している。本稿でもこの手法を採用する。

2・1・1　改良版時間区分法による温度近似式作成

ここでは,サンプリングされた離散値で与えられる測定温度履歴の近似式f_k(t)の扱いを説明する。

ある時刻t(>0)の解析領域の温度,熱流束の推定に用いるf_k(t)は,時間領域[0,t]の全体([t,∞]の未来データも一部含める)の測定温度履歴データに対する近似式とする。ただし,拡散型の熱伝導方程式の解は時間と伴に減衰する性質を有するので,tから過去に遡るある一定区間内[t₀, t]の測定温度に対して高い近似精度が保証できれば十分で,時刻tにおいて減衰してしまう遠い過去の温度履歴の近似精度は必要としない。この性質を活用した近似式作成法として時間区分法が提案されている^8,9)。

過去の時刻t₀近傍の温度履歴が現在時刻t(>t₀)に及ぼす影響の減衰は,無次元減衰時間である式(8)の減衰フーリエ数Fo_decayで与えられ,減衰に必要な時間は式(9)の減衰時間Δt_decayで与えられる。

  

$F{o}_{decay}\equiv \frac{a\left(t-t_0\right)}{L^2}$(8)

  

$\Delta t_{decay}\equiv t-t_0=F{o}_{decay}{L}^2/a$(9)

固体内の温度波振幅の減衰⁸⁾は,式(1a)の固有値β_n(次数n=0~∞)を用いてe−βⁿ² Fo^decayで規定される。Fo_decay,つまりΔt_decayの増加と伴に指数関数的に減衰が大きくなるので,t₀以前の温度履歴の影響は無視できる。従って,半値近似式を適用するデータを全範囲[0, t]からより狭い時間区間[t -Δt_decay, t]に限定できるため近似精度を容易に向上できる。例えば,第1種境界条件の平板に対する最小固有値β₀=π/2に対するFo_decayは3.733となり,これ以上に取れば区間[0, t₀]の温度履歴の影響を10^-4以下にできる。以後,本解析でのFo_decayはこの値を用いる。厚さL=2 mm,温度伝導率a=5.37×10^-6 m²/sの炭素鋼ではΔt_decay=2.8 sとなり,f_k(t)は時刻tから2.8 sだけ過去に遡る時間区間の測定データに対して作成すれば十分となる。

ところで,冷却中にクエンチを生じる鋼材冷却では不規則性が強い温度履歴を示すので,区間Δt_decay内の温度履歴に対する近似式を適用する時間区分法でも十分な精度が得られない。Woodfieldら⁸⁾は時間区間の温度履歴に対して複数の半値多項式(t^1/2の多項式)を重ね合わせて近似精度を向上させる改良時間区分法を提案している。Fig.2に示すように時間区間幅Δt_decayを重複するN_corr個の小区間:[(t_o,1)_n, t_n+Δt_min/2]…[(t_o,_j)_n, t_n+Δt_min/2]…[(t_o,_Ncorr)_n, t_n+Δt_min/2]に対するそれぞれの近似式を重ね合わせて精度を改善する。ただし,時間を過去に遡るほど温度履歴が現在に及ぼす影響が小さくなる熱伝導固有の性質に基づき,重複させるj番目の小区間の始点(t_o,_j)_nは,式(10)で表されるように,区間幅が公比r(<1)の単調減少等比数列に設定し,時刻が推定時刻t_n(= nΔt)に近づくにつれて,より多くの小区間を重ねて近似精度を高める。

  

${\left(t_{o,j}\right)}_n=t_n-\Delta t_{decay}{\left(1-r\right)}^{j-1}$(10)

Fig. 2.

Comparison of measured data at x=1 mm with superimposed approximate function for time range of [(t_o,1)_n¸ t_n + Δt_min/2] = [5.235 s, 8.080 s].

近似式を重ね合わせる小区間数N_corrは,区間幅がΔt_minを下回らない最大の整数値に設定する。Δt_minは任意に設定可能な最小時間区間幅で,式(11b)で示すN次の半値多項式の未定係数を決定にはN+1個以上のデータが必要なため,Δt_min≧(N+1)Δt(Δtは温度履歴サンプリング周期)の下限値制約がある。

以上より,改良時間区分法ではn番目の温度サンプリング時刻:t_n=nΔtの測温点番号k(=1, 2)におけるΔt_decay区間内の温度近似式(11a)は,小区間に対するN次の半値多項近似式(11b)を重ね合わせる。式(11b)の関数形は,熱伝導方程式の基本解がt^1/2で表現されるので少ない次数でも高い近似精度が得られる。また分母のガンマ関数はラプラス変換後の温度f ̃_k(s)の表記を簡単にするため,予め割られている。

  

$f_{k,n}\left(t\right)={\displaystyle \sum _{j=1}^{N_{corr}}{f}_{j,k,n}}\left(t\right)$(11a)

  

$$f_{j,k,n}\left(t\right)=\{\begin{array}{cc}{\displaystyle \sum }_{i=0}^N\frac{b_{i,j,k,n}{\left\{t-{\left(t_{o,j}\right)}_n\right\}}^{\frac{i}{2}}}{\Gamma \left(i/2+1\right)}& ;\left(t\geqq {\left(t_{o,j}\right)}_n\right)\\ 0& ;\left(t<{\left(t_{o,j}\right)}_n\right)\end{array}$$(11b)

近似式(11b)の係数b_i,j,k,nは,測温点番号k毎にj=1の小区間より順々に求めるが,その手順を以下に述べる。1)j=1の小区間[(t_o,1)_n, t_n+Δt_min/2]の測定温度に対して最小二乗法で近似式f_1,k,,n(t)の係数b_i,1,k,nを求める。2)j=2の小区間[(t_o,2)_n, t_n+Δt_min/2]での前ステップの近似式の残差R_k,2(t)= f_k,n(t)- f_1,k,,n(t)に対して最小二乗法でf_2,k,,n(t)の係数b_i,2,k,nを求める。3)以後,j=3~N_corrまで2)と同様の計算を繰り返す。

Fig.2にクエンチ温度履歴と上述のアルゴリズムで構成された近似式(11a)との比較の一例を示す。温度履歴(破線)は,後述するFig.3の急冷境界条件(Type1)に対して数値計算で評価されたx=1 mmでの温度変化,近似式(実線)は,最大の近似誤差を記録した時刻t_n=8.055 sの時間区間に対する近似式をそれぞれ示す。 近似式決定に必要なパラメータは,最高次数N=3, 最小時間幅Δt_min=0.05 s,減衰フーリエ数Fo_decay=3.733,N_corr=9となる。Fig.2よりt_nから時間を遡るに従い分割した小区間の開始時刻付近での近似式の偏倚を生じるが,時刻t=t_nに近づくにつれて近似式と測定温度履歴は非常良く一致している。

Fig. 3.

Boundary condition (Postulated boiling curves).

2・1・2　時間領域での温度・熱流束の厳密解

Mondeら⁷⁾の解析的手法に従い,複素平面上での温度および熱流束の逆問題解式(6),(7)をラプラス演算子sのべき級数展開したのち,項別にsのべき乗のラプラス逆変換を施すと時間領域での厳密解式(12),(13)が得られる。係数C_m,j,k,,n, D_m,j,k,,nは近似式の係数b_i,j,k,nとカーネル関数をsのべき級数展開した項の係数との積から導かれる定数である。

  

$T_n\left(t_n,x\right)={\displaystyle \sum }_{j=1}^{N_{corr}}\left\{{\displaystyle \sum }_{m=-1}^N\frac{C_{m,j,1,n}{\left\{t_n-{\left(t_{o,j}\right)}_n\right\}}^{m/2}}{\Gamma \left(m/2+1\right)}-{\displaystyle \sum }_{m=-1}^N\frac{C_{m,j,2,n}{\left\{t_n-{\left(t_{o,j}\right)}_n\right\}}^{m/2}}{\Gamma \left(m/2+1\right)}\right\}$(12)

  

$\frac{q_n\left(t_n,x\right)}{\lambda }={\displaystyle \sum }_{j=1}^{N_{corr}}\left\{{\displaystyle \sum }_{m=-1}^N\frac{D_{m,j,1,n}{\left\{t_n-{\left(t_{o,j}\right)}_n\right\}}^{m/2}}{\Gamma \left(m/2+1\right)}-{\displaystyle \sum }_{m=-1}^N\frac{D_{m,j,2,n}{\left\{t_n-{\left(t_{o,j}\right)}_n\right\}}^{m/2}}{\Gamma \left(m/2+1\right)}\right\}$(13)

温度・熱流束の厳密式(12),(13)は,x=0を代入して表面温度・表面熱流束を評価できるだけでなく,任意のx∈[0, L_x]を代入して内部の温度・熱流束変化も評価できる。ただし,測温点から位置が離れるに従い逆問題解の推定精度は一般に劣化する。本手法では,逆解析の時間ステップ毎に測定温度履歴の近似関数を決定する手間が掛かるが,代数計算式で与えられる厳密解は,反復計算を必要としないので解析時間は比較的短い。

なお,一般に高温鋼材は表面に数十~数百μmのスケールが生成しているので,スケールの熱抵抗層内に温度差を生じる¹⁰⁾。均質熱物性に対する本解析法を用いて実測定データから推定した表面温度・表面熱流束の解は,スケールと母材金属との界面温度・界面熱流束であることを指摘しておく 。

2・2　変物性の取り扱い

温度変動幅が大きくなる高温冷却条件では,材料の熱物性値ρ,c,λの温度依存性が無視できなくなる。変物性における熱伝導方程式(14)は式(15)で表すことが出来るが,非線形方程式となり,2・1節で示した解析的解法を適用できない。そこで定物性値に対する厳密解に対して,時間ステップごとに温度依存の均質な物性値を更新する簡易的な変物性の取り扱いを提案し,数値計算で得られた変物性の順問題解を用いて逆問題解の推定精度を検証する。

  

$\rho \left(T\right)c\left(T\right)\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left(\lambda \left(T\right)\frac{\partial T}{\partial x}\right)$(14)

  

$\frac{1}{a\left(T\right)}\frac{\partial T}{\partial t}=\frac{{\partial }^{2}T}{\partial x^2}+\frac{1}{\lambda \left(T\right)}\frac{\text{d}\lambda }{\text{d}T}{\left(\frac{\partial T}{\partial x}\right)}^2$(15)

式(12),(13)の係数C_m,j,k,n, D_m,j,k,nに含まれる温度伝導率aや式(14)のρcとλは,時間ステップごとに変化する固体の代表温度で評価する。熱物性値を評価する代表温度として,逆問題解析領域[0, x₁]の平均温度など考えられるが,いずれにせよ固体内の温度分布に基づく熱物性分布を一定値と見做すので,ここでは反復計算が不要な冷却面直下の測温点x=x₁での測定温度で評価した熱物性値を与える。なお,式(15)の右辺第2項の非線形項を無視すると均質一定熱物性の逆問題解式(12),(13)がそのまま適用できる。

この評価では,表面近傍の逆解析領域[0, x₁]の最大温度で評価した熱物性値に基づく解析となる。従って,表面付近の温度勾配が小さくなるほど,或いはx₁, x₂の位置が表面に近づくにつれ,熱物性値はより正確な解析領域の熱物性値を与える。簡易的な変物性の取り扱いによる逆問題解の精度は,後述する数値実験で評価する。

3. 数値実験による逆問題解析パラメータの検討

逆問題解析手法の推定精度は,既知の境界条件での検証が必要である。しかし,境界条件の真値が不明な実験データでは,精度検証が不可能なので,規定した境界条件に対する非定常熱伝導解析で得られた内部温度履歴を用いた数値実験により,逆問題解析の推定精度を検証する。

3・1　解析条件

一次元変物性の非定常熱伝導方程式(14)をPatanker¹²⁾のコントロールボリューム法で離散化した差分方程式を解いて内部温度を求めた。冷却開始の初期温度T₀=950°C,冷却水温度T_l=20°C,板厚L_x=10 mmとし,Al脱酸の0.08%C炭素鋼(以後,単に炭素鋼と記す)とSUS304の2鋼種を対象とし,報告されている高温物性値¹³⁾を用いた。変態潜熱は比熱として扱われる。境界条件は,スケールが存在しないx=0の冷却面上でFig.3に示す沸騰曲線を与え,裏面x=10 mmは断熱条件とした。Fig.3の沸騰冷却特性は,膜沸騰域,遷移沸騰域,限界熱流束域,核沸騰域,単相熱伝達域の表面熱伝達率を表面温度の折れ線関数で表現し,各領域の熱伝達率および温度域を変えて,強冷却,中冷却,緩冷却の条件を模擬したType1, Type2, Type3の3つの冷却条件で検証する。本稿での検証は,主に最も逆解析条件が厳しくなるType1での結果を代表例として提示する。ここでは,まず初期温度T=T₀での定物性条件での検証結果を示し,変物性条件での検証は4章で述べる。

3・2　最小時間幅が近似式品質と逆問題解析結果に与える影響

冷却速度が最速となるType1条件で求めたx=1, 2 mmでの炭素鋼の温度履歴T_trueとその近似式(11a)にΔt=5 msの時間ステップ毎に時刻t_nを代入して温度履歴を再現した近似温度変化T_appとの比較をFig.4(a)に示す。x=1, 2 mmでの近似式の絶対誤差ΔT_err(=T_app-T_true)をそれぞれFig.4(b)とFig.4(c)に示す。図中の各近似式は,最小時間幅Δt_min=0.05, 0.1, 0.5 sの条件で決定したもので,それ以外のパラメータ,最高次数N, 減衰フーリエ数Fo_delayは図中に示す推奨値とした。Fig.4(a)において,Δt_minに関わらずT_appとT_trueは,よく一致しているように見えるが,Fig.4(b),(c)の誤差の拡大スケールより,t=8 s付近の急冷開始点で高々数Kの誤差を生じている。温度変化が大きいx=1 mmの方がより大きな誤差を示す。また,同じ測温位置ではΔt_minを小さくするに従い,近似誤差を小さくできることが分かる。

Fig. 4.

Comparisons of measured data T_true and approximation curve T_app, and approximation error ΔT_err at two depthes. (effect of minimum time window size Δt_min).

次にType1炭素鋼のx₁=1 mm, x₂=2 mmでの温度変化を用いて,最小時間幅Δt_minが逆解析結果に与える影響を検証する。Fig.5(a)に表面温度T_wの推定値と真値との比較,T_wの推定値から真値を差し引いた絶対誤差ΔT_errを示す。Fig.5(b)に表面熱流束q_wの推定値と真値との比較を示す。

Fig. 5.

Effect of minimum time window size Δt_min on inverse solution for each Δt_min of (a) Surface temperature T_w and estimated error ΔT_err and (b) Surface heat flux q_w.

Fig.5(a)よりT_wの誤差ΔT_errは,Δt_min=0.1 s以下で最大6 K程度だが,0.5 sでは53 Kと非常に大きい。一方,Fig.5(b)でのq_wの推定精度も同様で,Δt_min=0.1 s以下ではほぼ真値と一致する結果が得られている。なお,Δt_minを 0.1 sから0.05 sに小さくしてもT_wの推定精度の改善に差は見られないが,q_wでは精度の向上が見られる。境界条件Type1においてクエンチ発生後の熱流束の急上昇に要する時間は0.15 sだが,その1/2程度のΔt_minを設定すれば,概ね表面上の非定常現象を逆解析で再現できることが分かった。さらにΔt_min=0.025 sでの解析結果を見ると0.05 sでの解析精度からの改善は見られなかった。

Fig.5(b)のピーク熱流束発生時刻(t=7.85 s近傍)に着目すると,Δt_minをいくら小さくしても熱流束ピーク検出位置に0.015 sの遅れが発生することが分かった。この遅れは3・4節で述べるが表面から測温点までの温度波の伝播遅れが原因で,これより短い周期の温度振動が表面で発生しても,測温点での検出が不可能なためである。従って,最小時間幅Δt_minは対象となる表面温度変動の速さに応じた時間分解能と測温位置から決まる最小予測時間の大きいほうを目安に設定すれば良い。つまり,高速現象を逆解析で捉える場合は,その速さに応じて表面からの測温点の深さ位置とその温度履歴を近似式で再現できるサンプリング周期を適切に選定する必要がある。

3・3　半値多項式の最高次数Nが近似式品質と逆問題解析結果に与える影響

Type1,炭素鋼のx₁=1 mm,x₂=2 mmでの温度変化を用いて,近似式の最高次数N=1~4が逆解析結果に及ぼす影響をFig.6に示す。Fig.6(a)は,T_wの推定値と真値との比較,T_wの推定値の絶対誤差ΔT_errを示し,Fig.6(b)はq_wの推定値と真値との比較を示す。

Fig. 6.

Effect of maximum order of polynomials N on inverse solution for each N of (a) Surface temperature T_w and estimated error ΔT_err and (b) Surface heat flux q_w.

Nに関わらずT_w,q_wの推定値は真値と比較的良い一致を示すが,低次の近似式ではT_wの推定結果に時間遅れが顕著となる。Nの増加により時間遅れは解消されるが,N=3以上で改善効果は飽和し,ΔT_errは最大6 Kとなる。ただし,真値に対する推定結果の最大偏倚を示す時刻7.85 sでの冷却速度は1419 K/sと非常に速いため,5 msの時間遅れで7 Kの温度誤差を生じることを指摘しておく。一方,N=3, 4でのq_wは5~10 ms程度の時間遅れで真値に追従しているが,最大熱流束付近でオーバーシュートを生じている。図示は省略するが,緩冷却条件のType3ではN=1, 2の低次近似式でも高次式と同じ精度が得られた。

以上より,温度変化が非常に速いクエンチ現象を捉える場合は,少なくともN=3で十分であることが分かった。なお,N >> 4の高次の近似式を用いると冷却条件によっては数値的不安定を生じて推定結果のハンチングが現れる場合もあった。このような場合は低次に留めておく必要であるが,経験的にN=2~4程度で十分な精度が得られる。なお,半値多項式の各項の振る舞いを考えると1次の項(t^1/2)は上に凸の関数,3, 4次の項(t^3/2, t^4/2)はそれぞれ下に凸と上に凸の関数となるため,3次程度を取れば区間内で上に凸・下に凸の不規則温度履歴のいずれに対しても近似式で表現可能となる。以後,特に断りがない限りN=4での結果を提示する。

3・4　測温位置が逆解析結果に与える影響

Type1,炭素鋼の条件で測温点間距離x₂-x₁=1.0 mmに固定してx₁を変化させたときに逆解析結果が受ける影響をFig.7に示す。Fig.7(a)は,T_wの推定値と真値との比較,T_wの絶対誤差ΔT_errを示し,Fig.7(b)はq_wの推定値と真値との比較を示す。

Fig. 7.

Effect of sensor depths x₁, x₂ of (a) Surface temperature T_w and estimated error ΔT_err and (b) Surface heat flux q_w.

Fig.7よりx₁が表面から離れるに従い,真値に対する推定値T_wとq_wの時間遅れが顕著となることが分かる。時間遅れの程度を評価するため,限界熱流束(ピーク熱流束)到達時刻の遅れをΔt_delayと定義し,Δt_delayをx₁²で整理した結果をFig.8に示す。図中の直線は,式(16)を満足する関係式である。式(16)は表面温度がt=0で突然変化するときの半無限固体の非定常熱伝導の厳密解で評価されるx₁の位置で初期温度から1%(= 0.01)の温度変化が検出可能となる最小予測時間Δt_{delay, min}の関係を示す。

  

$$\text{erfc}\left(x_1/2\sqrt{a\Delta t_{delay,min}}\right)=0.01$$(16)

Fig. 8.

CHF detection delay in an inverse solution.

Fig.8より逆解析での極大値検出時間遅れΔt_delayは,Δt_{delay, min}とほぼ一致することが分かる。つまり,表面温度変化が測定位置の深さx₁で検出可能となるΔt_dalayと温度伝導率aに応じて式(16)に基づき適切に選定する必要がある。

Fig.9にx₁=1.0 mmに固定し,x₂の位置のみ変化させた時の測温点間距離が逆解析結果に及ぼす影響を示す。x₂が4.0 mm以上の深さ,つまり測温点間隔が3 mmではT_wにハンチングが生じ,x₂=5.0 mmではq_wが負にアンダーシュートするなど推定精度の劣化が顕著となる。これは,x₂の測温点が表面およびx₁から離れるに従い,表面からの温度波の伝達時間遅れとx₁での温度変動に対する位相差が増大するためと考えられる。Fig.10に推定結果がハンチングしたx₁=1.0 mm,x₂=5.0 mmの測温位置での近似式の最高次数をN=2, 4としたときの推定結果を比較した結果を示す。前述の通りN=2の低次式を用いるとハンチングの大幅な抑制効果があるが,測温位置で決まる推定結果の時間遅れは改善できないことが分かる。

Fig. 9.

Inverse solution for each x₂ of (a) Surface temperature T_w and (b) Surface heat flux q_w.

Fig. 10.

Inverse solution for N=2 and 4 of (a) Surface temperature T_w and (b) Surface heat flux q_w.

以上より,測温位置x₁, x₂は表面温度変動が伝搬するまでの最小予測時間(式(16))を決定づける重要なパラメータであり,最も厳しい急冷条件となるType1での評価では,x₁=1 mmおよびx₁とx₂の測温点間隔を3 mm以下に選択することが推奨される。

4. 逆問題解析の精度検証

4・1　解析条件

3章と同様にType1~Type3の3つの境界条件に対する順問題解析結果に基づき逆解析の推定精度を検証する。ただし,対象鋼種を0.08%C炭素鋼とSUS304の2種類とし,定物性と変物性の各条件で解析を行った。近似関数決定に必要なパラメータは3章の結果を反映させて Δt =0.005 s,Δt_min=0.05 s,Fo_decay=3.733,N=4,x₁=1.0 mm,x₂=2.0 mmとした。

4・2　解析結果

Fig.11, 12に0.08%C炭素鋼とSUS304の逆解析結果を示す。それぞれ境界条件Type1~Type3における定物性と変物性条件での逆解析で再現された沸騰曲線を真値と比較した結果を示す。定物性の逆解析結果は鋼種に関わらず最大熱流束付近でのオーバーシュートを除いて真値とほぼ一致する。一方,変物性の結果は,定物性の場合と同様にオーバーシュートを示すが,炭素鋼では真値に対して過大評価,SUS304は過小評価側の推定結果を与える。これらの推定精度の低下は,支配方程式での変物性の影響を簡易的に取り扱っているためである。また,変物性条件での予測値の偏倚は,最大の冷却速度となるType1で最も大きく,冷却が緩慢となるType3では小さくなることが分かる。

Fig. 11.

Comparison of the inverse solution for 0.08%C steel.

Fig. 12.

Comparison of the inverse solution for SUS304.

精度が最も低下する限界熱流束の推定値の相対誤差を炭素鋼に対して評価すると,Type1の急冷条件で最大の26%を示し,Type2では10%,Type3では3%まで縮小する。なお,非定常沸騰冷却実験では,再現性や測定の不確かさによって誤差が数10%に達することも少なくない。従って,Type2条件で達成される限界熱流束3.45 MW/m²以下となる遅い冷却速度条件では,簡易的な変物性評価を用いた逆算法を用いても実用上許容できる程度の推定精度が得られることが分かった。

4・3　支配方程式の変物性の非線形項の影響

前節で変物性条件での沸騰曲線の推定誤差は,冷却速度の影響を強く受けることが分かった。温度変動幅が大きく変物性の影響が無視できない鉄鋼冷却実験に本解析法を適用する際には,変物性の影響による誤差解析に基づく適用条件を明らかにすることが不可欠である。

物性値は,逆解析の対象領域である表面x=0から表面に直近の測温点x=x₁の間の温度分布に応じて変化する。そこで,変物性の熱伝導方程式(15)の両辺に熱伝導率λを乗じた式(17)を用いて,変物性の影響を評価する。

  

$\rho c\frac{\partial T}{\partial t}=\lambda \frac{{\partial }^{2}T}{\partial x^2}+\frac{\text{d}\lambda }{\text{d}T}{\left(\frac{\partial T}{\partial x}\right)}^2$(17)

式(17)の左辺の物性値ρcに注目すると,x=0, x₁での代表温度に対するρcの物性値差に基づき,左辺の誤差Err_leftを式(18)で評価する。Err_leftより物性値差の影響は,冷却速度で増幅されることが分かる。

  

$$Er{r}_{left}=\left[\rho \left(T\left(0,t\right)\right)c\left(T\left(0,t\right)\right)-\rho \left(T\left(x_1,t\right)\right)c\left(T\left(x_1,t\right)\right)\right]\frac{\partial T}{\partial t}$$(18)

一方,式(17)右辺の誤差は,本解析で無視している非線形の第2項の式(19)で評価する。

  

$Er{r}_{right}=\frac{\text{d}\lambda }{\text{d}T}{\left(\frac{\partial T}{\partial x}\right)}^2$(19)

式(18),(19)の値を炭素鋼の境界条件Type1,核沸騰領域である350°C付近でのx=0, x₁(=1 mm)での温度に基づく直線温度勾配(=(T(x₁,t)- T(0,t))/x₁)と冷却速度で評価すると,Err_left=0.43×10⁸ W/m³, Err_right=3.22×10⁸ W/m³となり,式(17)の右辺の熱伝導率の温度依存項の誤差が支配的となることが分かる。

Fig.13にx₁=1.0 mmと表面での温度差(∝温度勾配)でq_wの推定誤差を整理した結果を示す。推定誤差は,Fig.3における限界熱流束域の中間温度(Type1: 350°C,Type2: 271°C,Type3: 193°C)での厳密解熱流束q_w,n.b.-exactと逆解析熱流束q_w,n.b.の比で評価した。Fig.13より逆解析誤差はx=x₁と0での温度差の増加とともに真値と一致する1から偏倚し,炭素鋼の誤差は温度差に対して正の勾配,SUSが負の勾配を示すことが分かる。この理由は,Fig.14に示す各鋼種のλの温度依存性から分かるようにdλ/dTが負を示す炭素鋼では表面熱流束が過大評価,正となるSUSでは熱流束が過小評価されるためである。Fig.13よりq_wの推定誤差が10%以内となるx₁(=1 mm)と表面との間の最大許容温度差は,炭素鋼で56 K,SUSで95 K,つまり温度勾配ΔT/Δxがそれぞれ56 K/mm, 95 K/mmを超えるとq_wの推定誤差は10%を超える。

Fig. 13.

Relationship between estimated error and temperature difference.

Fig. 14.

Thermal conductivity.

次にFig.15に熱流束逆解析推定の絶対誤差Δq_w,n.b.=q_w,n.b.-q_w,n.b.-exactを解析では無視されている式(19)の右辺第2項で整理した結果えお示す。なお,第2項の微分勾配は表面(x=0)とx₁との間の差分近似式-(Δλ/ΔT)(ΔT/Δx)²で評価した。Fig.15より誤差Δq_w,n.bは,右辺第2項と直線関係が示され,第2項が誤差を支配する主要要因であることが分かる。Type1の様にクエンチ時に1400 K/sを超える極めて冷却速度の大きい実験で逆解析を行う場合,熱伝導方程式の非線形項を考慮した解析が不可欠であることが分かった。非線形項のより厳密な取り扱いについては,稿を改めて報告する。

Fig. 15.

Relationship between estimated error and the nonlinear term of the heat conduction equation.

5. 結論

本研究では,ラプラス変換を用いた逆問題解析手法を簡易に熱物性の温度依存性を解析できるように拡張し,鋼材冷却に適用した際の近似精度を評価した。その結果,以下の結論を得た。

(1)近似関数の最小時間幅Δt_minは,逆解析で評価対象の変動周期の1/2程度に設定すれば十分な応答で解析が可能である。

(2)半値多項式の最高次数Nは3~4程度で十分高い近似精度が得られる。一方,冷却速度が1400 K/sを超える急速冷却条件で2点の温度測定位置が3 mm以上離れた条件では,推定値の発散抑制のため,低次(N=2)の設定を推奨する。

(3)温度測定位置が表面に近いほど,表面伝熱挙動を高い時間応答性で解析することが出来る。逆問題解析の応答遅れは式(16)で規定されるように,表面側測定位置x₁の2乗に比例して増加する。

(4)支配方程式の非線形項を無視した簡易的な温度依存性の取り扱い法を提案し,その手法を用いてq_wを10%の精度で推定可能となる条件として,表面温度と表面直下の測温位置での温度差を鋼種ごとに定めた。また,推定誤差は2つ鋼種に対して-Δλ/ΔT(ΔT/Δx)²で評価できることを示した。

記 号

t: 時刻[s]

s: ラプラス演算子

x: 表面からの距離[m]

T: 温度[K]

a: 温度伝導率[m/s²]

q: 熱流束[W/m²]

λ: 熱伝導率[W/mK]

f: 測温点の温度実測値,推定値[K]

L_x: 試験片厚さ[m]

L: 代表長さ(=x₂)[m]

Fo_decay: 減衰フーリエ数[-]

Δt_decay: 減衰時間[s]

N_corr: 改良版時間区分法の小区間数

N: 半値多項式の最高次数

r: 時間区分の等比数列分割の公比

t_o: 近似多項式定義開始時刻[s]

t_n: 解析対象時刻[s]

Δt: 温度サンプリング周期[s]

Δt_min: 最小時間区分幅[s]

A, B: 非定常熱伝導方程式の一般解の係数(式(3))

C, D: 熱伝導逆問題解の係数(式(12),(13))

b: 温度近似式の係数(式(11b))

<添字>

i: 小区間温度履歴を近似する半値多項式の次数

j: 改良版時間区分法の重ね合わせ近似式の次数

k: 測温点番号

l: カーネル関数の級数展開式の次数

m: 温度・熱流束の厳密解を表す級数の次数

n: 測温データのサンプル番号

w: 表面(x=0)

0: 初期値(t=0)

1, 2: 測温点番号

文献

• 1)  M.Mitsutsuka: Tetsu-to-Hagané, 79(1993), N405 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.79.6_N405
• 2)  S.Nishio and M.Uemura: Trans. Jpn. Soc. Mech. Eng. B, 52(1986), 1811 (in Japanese). https://doi.org/10.1299/kikaib.52.1811
• 3)  J.V. Beck: Int. J. Heat Mass Transf., 13(1970),703. https://doi.org/10.1016/0017-9310(70)90044-X
• 4)  M.Mitsutsuka: Tetsu-to-Hagané, 54(1968), 1457 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.54.14_1457
• 5)  K.Fukuda and T.Ariyoshi: Tetsu-to-Hagané, 70(1984), 869 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.70.8_869
• 6)  T.Otsuka and S.Ito: Tetsu-to-Hagané, 102(2016), 29 (in Japanese). https://doi.org/10.2355/tetsutohagane.TETSU-2015-067
• 7)  M.Monde, Y.Mitsutake, T.Shiki and T.Nishimoto:Trans. Jpn. Soc. Mech. Eng. B,66(2000), 519 (in Japanese). https://doi.org/10.1299/kikaib.66.519
• 8)  P.L. Woodfield, M. Monde and Y. Mitsutake: Int. J. Heat Mass Transf., 49(2006), 2864. https://doi.org/10.1016/j.ijheatmasstransfer.2006.01.050
• 9)  M.Monde, W.Liu, H.Arima, Y.Mitsutake and J.A.Hammad:Trans. Jpn. Soc. Mech. Eng. B,68(2002), 2093 (in Japanese). https://doi.org/10.1299/kikaib.68.2093
• 10)  H.Fukuda, N.Nakata, H, Kijima, T.Kuroki, A.Fujibayashi, Y.Takata and S.Hidaka: Tetsu-to-Hagané, 100(2014), 1514 (in Japanese). https://doi.org/10.2355/tetsutohagane.100.1514
• 11)  Y.Umehara, T.Okawa and K.Enoki: Tetsu-to-Hagané, 105(2019), 1050 (in Japanese). https://doi.org/10.2355/tetsutohagane.TETSU-2019-035
• 12)  S.V.Patanker: Computation of Conduction and Duct Flow Heat Transfer, CRC Press, Boca Raton, (1991).
• 13)  Heat transfer experiment and calculation algorithm in a continuous reheating furnace, ISIJ, Tokyo, (1970), (in Japanese).