非線性逆熱傳導問題之理論發展與實驗驗證(I)

行政院國家科學委員會專題研究計畫 期中進度報告

非線性逆熱傳導問題之理論發展與實驗驗證(1/2)

計畫類別： 個別型計畫 計畫編號： NSC92-2212-E-151-007- 執行期間： 92 年 08 月 01 日至 93 年 07 月 31 日 執行單位： 國立高雄應用科技大學模具工程系 計畫主持人： 楊慶煜 報告類型： 精簡報告 處理方式： 本計畫可公開查詢

中 華 民 國 93 年 5 月 31 日

行政院國家科學委員會專題研究計畫成果期中報告 非線性逆熱傳導問題之理論發展與實驗驗證(1/2) 計畫編號: NSC92-2212-E-151-00 7-執行期限: 92 年 8 月 1 日至 93 年 7 月 31 日 主持人: 楊慶煜 國立高雄應用科技大學模具工程系 中文摘要 逆運算熱傳導問題， 一般是應用 熱電耦或紅外線光學法量得固體內部 點 或 表 面 的 溫 度 值 ， 以 反 求 邊 界 條 件、熱傳係數、表面熱傳量、內部熱 源等。逆解熱傳導問題方法已被廣泛 的應用在許多設計與製造的問題，尤 其是當物體表面之狀態無法直接量測 時。例如量測機械切削時刀具界面的 溫度或熱傳量、量測燃燒室內壁溫度 或內外壁之熱傳遞係數、高速飛行器 之表面溫度等等。在過去的逆運算領 域的研究中，多數學者發展的方法僅 限於具線性熱力性質的系統，亦即線 性的熱傳系統。然而多數的工程實際 問題中，熱力性質大多與溫度相依， 熱力系統呈現非線性的性質，影響系 統的溫度分佈、熱傳現象與熱穩定的 問題甚巨，也因此造成非線性逆熱傳 問題邊界估算上的困難。原因在於非 線性熱力性質隨溫度變化的特性，影 響熱系統溫度輸出的變化甚大。因此 發展一精確、穩定與快速的非線性逆 熱傳邊界估算方法，為本計畫預計達 成之目標。本計畫以兩年時間進行二 維暫態非線性逆熱傳邊界的熱行為預 測的理論架構與實驗驗證。第一年的 研究重點結合數值運算與未來時間的 概 念 之 循 序 方 法 ， 同 時 結 合 修 正 Newton-Raphson 法於每一時間區間求 解邊界條件。在第二年的執行期中， 將進行實驗設備的架構與實驗流程的 設計，並驗證第一年的理論結果。 關鍵詞: 逆運算、非線性熱傳導、修正 Newton-Raphson 法 Abstract

The determination of the boundary conditions from the measured temperature profiles is the inverse boundary estimation problems. The boundary estimation problems have been widely applied in many design and manufacturing problems especially when the direct measurements of the surface conditions are not possible such as the measurement of the heat flux or temperature at the inner surface of a heated pipe, at the inside of a combustion chamber, at the outer surface of a reentry vehicle, or at the tool-work interface of a machine cutting. In the past researches, most of the inverse techniques are confined to the problems with temperature independent thermal properties, which is the linear heat conduction equation. However, most of the practical engineering problems, the thermal properties are temperature dependent and lead the heat conduction equation to a nonlinear form. The determination of the thermal system with temperature dependent properties is more difficult than that with the temperature independent properties because the magnitude of the thermal properties has a significant influence on the analysis of temperature distribution, heat flow rate, and thermal instability problems. Therefore, it is necessary to develop an inverse algorithm in order to estimate the boundary condition in the field of the nonlinear inverse heat conduction problem.

Keywords: nonlinear heat conduction, modified Newton-Raphson method

緣由與目的

逆 運 算 熱 傳 導 問 題 (Inverse heat conduction problem) 簡稱為 IHCP， 一 般是應用熱電耦或紅外線光學法量得 固體內部點或表面的溫度值，以反求 邊界條件、熱傳係數、表面熱傳量、 內部熱源等。逆解熱傳導問題方法已 被廣泛的應用在許多設計與製造的問 題，尤其是當物體表面之狀態無法直 接量測時。例如量測機械切削時刀具 界面的溫度或熱傳量、量測燃燒室內 壁溫度或內外壁之熱傳遞係數、高速 飛行器之表面溫度等等。這些問題即 可 歸 納 為 逆 運 算 問 題 (Inverse problem)。 在過去的逆運算領域的研究中， 多數學者發展的方法僅限於具線性熱 力性質的系統，亦即線性的熱傳系 統。然而多數的工程實際問題中，熱 力性質大多與溫度相依，熱力系統呈 現非線性的性質，影響系統的溫度分 佈、熱傳現象與熱穩定的問題甚巨， 也因此造成非線性逆熱傳問題邊界估 算上的困難。 一般而言逆熱傳估算結果的精確 度，是非常的不穩定，即使小小的量 測誤差也會造成無法接受的估算結 果。因此可採用正則化方法 (regularization method) [1]或未來時間觀 念(concept of the future time) [2]解決估 算不穩定的問題。正則化方法將待定 參數的變化量轉化成非線性均方根的 懲罰函數，而未來時間觀念是假設未 來幾個時間的待定參數量是已知的。 使用此二方法的相關研究，可參考文 獻 3-9。至於其他相關的逆運算方法如 疊代正規法(iterative regularization method)[5]、動態規劃法( dynamic programming) [10]、the mollification method [11]、共軛梯度法( adjoint equation approach coupled to the

conjugate gradient method )[12]、基因演 算法(genetic algorithms) [13]、符號循序 法( symbolic sequential method )[14]以及

數值循序法( numerical sequential method) [15]。然而上列方法多數僅能 用於線性熱傳系統，對於非線性熱傳 系統的估算僅有少數的研究曾探討過 [16-17]。 估算非線性熱傳問題的邊界，較 線性熱傳導系統困難。原因在於非線 性熱力性質隨溫度變化的特性，影響 熱系統溫度輸出的變化甚大。例如 Gregory and Daniel [17]使用牛頓法 (Newton method )估算線性與非線性逆 熱傳導邊界問題。對於具線性熱傳性 質的問題，在數次的疊代後即可找到 精確的結果(例題三，參考文獻 17)；但 是對於非線性熱傳性質的問題，牛頓 法或 BFGS–Newton 法 (Broyton-Fletcher-Goldfarb-Shanno combined with Netwon method)均無法找 到收斂的結果 (例題四，參考文獻 17) 。也因此發展一精確、穩定與快速 的非線性逆熱傳邊界估算方法，為本 計畫預計達成之目標。 本計畫擬以兩年時間進行二維暫 態非線性逆熱傳邊界的熱行為預測的 理論架構與實驗驗證。第一年的研究 重點擬結合數值運算與未來時間[9]的 概 念 之 循 序 方 法 ， 同 時 結 合 修 正 Newton-Raphson 法[18]於每一時間區間 求解邊界條件。研究中所用之方法可 避免以非線性最小均方根法(Nonlinear least-squares error) 架構估 算 熱邊界 問 題，避免增加問題的非線性。求解步 驟中以混合有限元素與有限差分法結 合未來時間，循序解出待估測之物理 量，再予以穩定化。除了上述之優點 外，本方法能將未知條件假設可以不 假設任何的形式直接由計算機算出對 應數值。這對於處理較複雜形式的未 知條件，或是未知條件的形式難以預 測時，均可準確的描述其結果。在第 二年的執行期中，將進行實驗設備的 架構與實驗流程的設計。本實驗的目 的是藉由試件邊界溫度量測估算熱邊 界。因此試件必須裝置熱電偶，將熱

電偶連接至數據擷取系統測量溫度， 以估算邊界值。同時需要電源控制器 透過加熱片在邊界處提供熱通量。為 了確定所輸入熱通量之大小，需由瓦 特計來量測並監控，同時以 AD/DA 卡 來調整熱通量的輸入。為了能控制各 個儀器的運作與量測，將自行撰寫控 制 程 式 。 透 過 實 驗 以 真 實 的 溫 度 分 佈，估算二維非線性逆熱傳導問題邊 界，並驗證第一年的理論結果。 總之本研究以新的方法，估算二維非線 性逆熱傳導問題邊界，並以實驗驗證理 論的可行性。透過理論與實際的結合， 能使逆運算法，推廣應用於實際工程問 題領域。 問題描述 考慮二維物體邊界條件有(1) 溫 度 T = T_b on Γ_T, (2) 熱傳量 q_n = q₀ on Γ_q (3) 熱 對 流 q_n = h(T − T_f) on .此物體內部定義為V,物體邊界定 義為Γ = Γc Γ_T∪ ,其數學表示 式可寫為: Γq ∪ Γc t T T C y T T k y x T T k x x y ∂ ∂ = ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ ) ( ) ) ( ( ) ) ( ( ρ (xi,yi)∈V (1) ) , , ( ) , , (x y t T x y t T _{i i} = _{b i i} (2) T i i y x, )∈Γ ( ) , , ( ) , , (x y t q₀ x y t q_{n i i} = _{i i} (3) q i i y x, )∈Γ ( )] ( ) , , ( [ ) , , (x y t hT x y t T t qn i i = i i − f c i i y x, )∈Γ ( (4) ) , ( ) 0 , , (x_i y_i T₀ x_i y_i T = V y x_i, _i)∈Γ∪ ( (5) T表 示 溫 度 分 佈T (x, y, t) 。 和 表示熱傳導係數與 ) (T k_x ) (T k_y ρC(T)是單 位容積之熱容量係數。本研究中以方 塊邊界之量測溫度來反求未知邊界大 小。 理論分析 計 畫 中 擬 於 每 一 時 間 疊 代 過 程 中，估算同一時間點的未知邊界。因 此求解過程中，包括了直接問題、敏 感度問題、非線性方程組問題與中止 疊代的問題。現分述如下。 計畫中以循序方式估算未知邊界，因 此逆運算問題被侷限於一個時間點。 所以方程式(1-5)被限制在某一特定時 間t= 。因此方程式可改寫為 t_m t T T C y T T k y x T T k x m m m m y m m x _∂ ∂ = ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ ) ( ) ) ( ( ) ) ( ( ρ V y x_i, _i)∈ ( at t = (6) t_m T i T m m i i y t x T( , , )=φ , T i i y x, )∈Γ ( (7) q i q m m i i n x y t q ( , , )=φ , q i i y x, )∈Γ ( (8) ] ) , , ( [ ) , , ( c,ic m m i i m i i n x y t hT x y t q = −φ c i i y x, )∈Γ ( (9) 1 1) , , (xi yi tm₋ =Tm₋ T V y xi, i)∈Γ∪ ( (10 在此 Tm =T(xi,yi,tm) ),, (0 , mi i iqmqqtyx =φ ) , , ( , m i i b i T m T x y t T = φ ) , , ( 0 , m i i i q m q x y t q = φ ) ( , m f i c m T t c = φ T i T m , φ , qiq, 為未知的邊界 m , φ cic m , φ 逆 運 算 問 題 是 屬 於 非 適 定 問 題 （ill-posed），尤其對非線性問題，非適 定現象更為明顯。因此需要引用未來 時間的觀念來穩定化估算結果。例如 當t =t_m時 ， 邊 界 條 件 在 t=t₁ 與 1 − =tm t 間已循序求出。估算 時的 邊 界 需 要 做 如 下 的 假 設 （ concept of future time）: m t t = T T T T T Ti m i T r m i T r m i T m i T m , , 1 , 2 , 2 , 1 φ φ φ φ φ ₊ = ₊ =_L= ₊₋ = _{+ −} = q q q q q qi m i q r m i q r m i q m i q m , , 1 , 2 , 2 , 1 φ φ φ φ φ ₊ = ₊ =_L= ₊₋ = ₊₋ = c c c c c ci m i c r m i c r m i c m i c m , , 1 , 2 , 2 , 1 φ φ φ φ φ ₊ = ₊ =_L= ₊₋ = ₊₋ =

(11) 在此 r 未知時間數 當邊界條件假設已知，直接解可 代入方程式(6-10)求得物體內之溫度分 佈。求解的方法是將方程式(6-10)，以 有限元素法予以離散化，直接求解。 由於式(6)為非線代方程組，因此以 Newton-Raphson 法疊代求解。在本部份 所得到的解答再代入敏感度問題中， 進行敏感度分析。 敏感度問題 以疊代的方式求解問題時，需要 透 過 ， 求 得 每 次 疊 代 過 程 中 的 修 正 量。為了求得必要之修正量需先將方 程式對待定量求微分。因此將(6-10)式 取 * *,i m ∂φ ∂ 可得到: t X T C y X T k y x X T k x m m m m y m m x _∂ ∂ = ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ ) ( ) ) ( ( ) ) ( ( ρ (12) V y x_i, _i)∈ ( t =t_m * *, , ) , , ( _i m i T m m i i T t y x X φ φ ∂ ∂ = (13) T i i y x, )∈Γ ( * * *, , *, ) , , ( i m i q m i m m i i n q t y x q φ φ φ ∂ ∂ = ∂ ∂ (14) q i i y x, )∈Γ ( * * *, , *, ) , , ( i m i c m i m m i i n c h t y x q φ φ φ ∂ ∂ − = ∂ ∂ (15) c i i y x, )∈Γ ( 0 ) 1 , , (x_i y_i m− = X_m−1 = X (16) V y x_i, _i)∈Γ∪ ( 在此 * *, ) , , ( i m m i i m t y x T X ∂φ ∂ = 方程式 (12-16) 為求解敏感度係 數的方程組，當 和 已知時即可 順利求出。 * *, i m φ T_m 修正 Newton-Raphson 法 本研究中，非線性問題是由量測 溫 度Φ_meas( , )i j 與 計 算 溫 度Φ_c( , )i j 比

對。因此假設在i -grid 空間與 j -grid 時 間時之計算與量測溫度整理成式(17): Φ( , )i j =Φc( , )i j −Φmeas( , )i j =0 (17) where j =1, 2,3 ,..., N_t and N_t是時間 格點。在式(17)中 Nt個方程當僅有一個 熱電偶使用時。詳細步驟如下： 將時間指標 j 由 1 到 Nt代入式(17) ， Φ=

[

Φ( , ), ( , ), ( , ),..., ( ,i1 Φ i2 Φ i3 Φ i Nt)

]

T =

{ }

Φ$_u (18) 待定係數可表為如下： 式(18)中含有未知邊界 , , ，因此求解邊界可透過式(18)求 得 。 在 式 (18) 中 共 有 個 變 數，共有 T i T m , φ qiq m , φ c i c m , φ c q T n n n + + p r 個方程式( p 為感測器數 目， r 為未來時間數)。將式(18)中指標 從 代到 j m m+ r−1和指標 i 從1代到 p ，可得到詳細步驟如下： Φ = T r m i m i m i m i, ), (, 1), (, 2), , ( , 1)] ( [Φ Φ + Φ + LΦ + − )] 1 , ( , ), 2 , ( ), 1 , ( ), , ( | ) 1 , 2 ( , ), 2 , 2 ( ), 1 , 2 ( ), , 2 ( ) 1 , 1 ( , ), 2 , 1 ( ), 1 , 1 ( ), , 1 ( [ T r m p m p m p m p r m m m m r m m m m − + Φ + Φ + Φ Φ − + Φ + Φ + Φ Φ − + Φ + Φ + Φ Φ = L L L L L L L L L L L L L L L L L L L L L L L L ｜ =

{ }

Φ$u (19) 式中Φ$_u對x_v之導數可求得並表示為： Ψ_{u v} Φu v x , $ = ∂ ∂ (20) 敏感度矩陣Ψ可定義為： Ψ=

{ }

Ψ_{u v,} (21) u =1, 2, 3,..., Nt 與v =1, 2, 3,..., m +n ，Ψu v, 是Ψ在 u− th列v 行的 元素，初始假設值 = − th x0

[

x1,0, x2,0, x3,0, x4,0,..., xm + n ,0

]

T (22) 最後可得到: x_{λ +1}= x_λ + ∆_λ (23) Ψ(x_λ)∆_λ = −Φ(x_λ) (24) ∆_λ=−[ΨT (x_λ)Ψ(x_λ)]−1ΨT(x_λ)Φ(x_λ) (25) 上述的推導，僅以一點量測為

例，可以很容易推廣到多點量測。例

如若有p個空間量測點，則式(19)中將

有 p× Nt個元素，亦即 p×Nt個方程

式。

The stopping criterion

研究中的疊代過程（式 23-25）是 用來求得未知向量x。從 到x 的修 正量 是由式(23)求得。因此疊代過程 中 x 不斷的被新的∆ 修正，直到下 列條件滿足： x_{λ λ +1} ∆_λ λ +1 λ J(x_{λ +}1) <ε (26) J(x_{λ +}1) = ∆t Nmeas×Nt [ j=1 N_t

∑

i=1 Nmeas

∑

Φ_c( , )i j − Φmeas( , )]i j 2 (27)

ε 是 stopping criterion 以及Nmeas是 t 空間量測點 當未考慮量測誤差時，ε可假定 為一極小的值。但實驗過程中，不容 易控制在零誤差的情況因此差異性觀 念(discrepancy principle)，可用來設定停 止疊代的條件，在此溫度餘數(residuals) 可表為 Φmeas( , )i j −Φc( , )i j ≈ σ (28) 其中σ是量測標準差，在特定的實 驗條件下，為一常數。 所以 ε = ∆tσ2 (29) 例題 圖一：一中空長管之剖面 例題一：考慮一長管具有如上圖之幾 何形狀，其熱傳導係數與單位熱容量 係數分別為 C m W T T k_x =1.5+0.01 +0.002 2 / o ， ) ( 01 . 0 4 . 1 Sin T k_y = + ρC= 0.15 J m3 0C 。長管的外側為絕 熱。初始溫度為T₀ = 20 0C，突然有液 體溫度為 T_f(t) (h= 50 W m2 0C) 流 經長管內側。溫度函數 T_f(t) 表示如 下: T_f(t)= 20 +100t 當 0< t≤1.0秒 T_f(t)= 220 −100t 當 1.0< t ≤ 2.0秒 此例題具有非線性之材料性質，並以 上列之輸入溫度來估算溫度邊界，以 驗證擬發展方法的可行性。圖一結果 顯示當σ=0.5 時，r= 4 有較佳之估算結 果。 圖二：管內溫度估算結果 σ =0.5, 2 = r and 4. 圖三：例題二之示意圖 例題二：: 考慮一長方板底部與兩側均 為絕熱，上方有熱傳量傳入。其初始 溫度為T0 = 0 0 C。長方板之熱傳導係 數與單位熱容量係數分別為 2 0005 . 0 00125 . 0 6081 . 0 T T k_x= + + W/m0C 2 0004 . 0 00725 . 0 7 . 0 T T ky = + + W m C 0 /

and ρC =0.9+0.001T J/m3 0C. 本研究之目的是從底部之溫度量測估 算輸入之熱傳量。 量測之溫度值如圖四所示。當不考慮 量測誤差時可估得正確之熱傳量（σ= 0 and r= 1）。當考慮誤差值為σ =0.02 時若取r= 1 結果偏離正確值甚遠。若 取 則可估得合理之熱傳值（見圖 五）。 3 = r 圖二與圖五的結果顯示,本研究 可 估 算 飛 現 線 系 統 之 溫 度 與 熱 傳 值。，且在較差之量測環境上亦可準 確估算熱傳導係數與熱容量係數。 圖四：例題二中之量測溫度 圖五：例題二中之估算熱傳值當 02 . 0 = σ 與r =3。 計畫結果自評 如期完成第一年預定之結果，以 新發展的方法(修正 Newton-Raphson 法) 同時估算非線性逆熱傳導問題之溫度 與熱傳量。從已完成的結果，可知適 當的未來時間可增進估算結果之穩定 性。在第二年的執行期中，將進行實 驗設備的架構與實驗流程的設計驗證 第一年的理論結果。 參考文獻

1. Tikhonov, A. N. and Arsenin, V. Y.,

Solutions of Ill-posed Problems, 1st ed., Winston, Washington, DC, 1977.

f

2. Beck, J. V., Blackwell, B. and Clair, C. R. St., Inverse Heat Conduction-Ill Posed Problem, 1st ed., Wiley, New York, 1985.

3. Morozov, V. A. and Stessin, M.,

Regularization Methods or Ill-Posed Problems, 1st ed., CRC Press, Inc., 1993.

4. Alifanov, O. M., Inverse Heat Transfer Problems, 1st ed., Springer-Verlag, Berlin Heidelberg, 1994.

5. Liu, J., "A Stability Analysis on Beck's procedure for Inverse Heat Conduction Problems," Journal of Computational Physics, Vol. 123, 1995, pp. 65-73. 6. Reinhardt, H. J., "A Numerical method

for the Solution of the

Two-Dimensional Inverse Heat Conduction Problem," International Journal for Numerical Methods in Engineering, Vol. 32, 1991, pp. 363-383.

7. Alifanov, O. M. and Artyukhin, E. A., "Regularized Numerical Solution of Nonlinear Inverse Heat Conduction Problem," Journal of Engineering Physics, Vol. 29, No. 1, 1975, pp. 934-938.

8 Alifanov, O. M. and Mikhailov, V. V., "Solution of the Nonlinear Inverse Thermal Conductivity Problem by the Iteration Method," Journal of

1978, pp. 1501-1506.

9. Beck, J. V. and Murio, D. A., "Combined Function

Specification-regularization Procedure for Solution of Inverse Heat

Conduction Problem," AIAA Journal, 1986, Vol. 24, pp. 180-185.

10. Busby, H. R. and Trujillo, D. M., "Numerical Solution to A

Two-Dimensional Inverse Heat Conduction Problem," International Journal for Numerical Methods in Engineering, Vol. 21, 1985, pp. 349-359.

11. Murio, D. A., The Mollification Method and the Numerical Solution of

Ill-Posed Problems ,

Wiley-Interscience, New York, 1993. 12. Jarny, Y., Ozisik, M. N. and Bardon, J.

P., "A General Optimization Method Using Adjoint Equation for Solving Multidimensional Inverse Heat Conduction," International Journal of Heat and Mass Transfer, Vol. 34, No. 11, 1991, pp. 2911-2919.

13. Li., H. Y. and Yang, C. Y., "A Genetic Algorithm for Inverse Radiation Problems," International Journal of Heat and Mass Transfer, Vol. 40, No. 7, 1997, pp. 1545-1549.

14. Yang, C. Y., “A Sequential Method to Estimate The Strength of The Heat Source Based on Symbolic

Computation," International Journal of Heat and Mass Transfer, Vol. 41, No. 14, 1998, pp. 2245-2252.

15. Yang, C. Y., “Inverse Estimation of Mix-typed Boundary Conditions in Heat Conduction Problems," AIAA Journal of Thermophysics and Heat Transfer, Vol. 12, No. 4, 1998, pp. 552-561.

16. Beck J. V., “Nonlinear Estimation Applied to The Nonlinear Inverse Heat

Conduction Problem,＂ International Journal of Heat and Mass Transfer, Vol. 13, 1970, pp. 703-716.

17. Dorai, G. A. And Tortorelli, D. A., “Transient Inverse Heat Conduction Problem Solutions Via Newton＇s Method,＂ International Journal of Heat and Mass Transfer, Vol. 40, No. 17, 1997, pp. 4115-4127

18. Yang, C. Y., “Determination of Temperature Dependent

Thermophysical Properties from Temperature Responses Measured ad Medium＇s Boundaries," International Journal of Heat and Mass Transfer, Vol.43, 2000, pp. 1261-1270