應用數值循序逆運算法估算二維多熱源問題與逆運算法之數學性質分析(II)

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

應用數值循序逆運算法估算二維多熱源問題與

逆運算法之數學性質分析 (2/2)

計畫編號: NSC89-2212-E151-028 執行期限: 89 年 8 月 1 日至 91 年 7 月 31 日 主持人: 楊慶煜 國立高雄應用科技大學模具系 中文摘要 在過去，近四十年中已有許多的 方法被發展出來，且成功的應用在熱 傳領域。但目前所使用之方法需要反 覆之疊代運算以及在非線性之領域中 求解。造成估算上時間的花費與準確 度之失真。本研究所討論之逆運算模 式包含逆向矩陣法，符號運算法與符 號循序方法與數值循序方法，可避免 疊代過程與避免在非線性領域求解逆 運算之問題。此四方法之特色在於當 連續之領域經離散化後並不需對題目 做任何簡化即可求解。透過所推導之 逆運算數學模式，求解逆運算問題之 過程可轉換為一組線性方程式來估算 未知之物理量。其特點為可直接求解 無需疊代且僅需求解一組線性方程式 就可估算未知之物理量。 逆向矩陣法與符號運算法屬於全 時域性的逆運算法，符號循序方法與 數值循序方法屬於循序性的逆運算 法。兩者的差別在於全域性的方法將 所有待定的參數在時間區間內匯集， 一次求解。而循序法是在每個時間求 出當時的待定參數。此種方法各有其 優缺點。本研究將針對逆向矩陣法、 符號運算法、符號運算之循序法與數 值循序法的計算法作一比較與分析， 同時提供熱傳逆問題與工程軟體結合 之方法。 關鍵詞: 逆運算，全時域方法，逆向 矩陣法，符號運算法，循序方法，未 來時間 Abstract

In the traditional researches, most of the inverse methods adopt the nonlinear least-squares formulation in solving the inverse problems. There are few disadvantages in the past approaches. One is that the iterative process in the computation not be avoided and the other is that the inverse problem can only be solved in a nonlinear domain.

To solve the above problems, two whole domain methods and two sequential methods, the reverse matrix method, the symbolic method, the symbolic combined with sequential method, and the numerical combined with sequential method, have been developed. Either the whole domain method or the sequential method is with the properties of both sides of advantage and disadvantage. In this research, the computational algorithms and the characteristic of the methods are discussed. For instance, the whole domain method needs more computational efforts but less effort to represent the undetermined parameters. As well, the sequential method is more efficiency but it needs more computer memory to distinct the parameters. In this study, an example, the estimation of the inner temperature of a fin is used to demonstrate the method to combine the inverse algorithm and a CAE software.

Keywords: whole domain method,

reverse matrix method, symbolic method, sequential method, concept of future time

緣由與目的 在熱傳導的問題中，若邊界條件、 初始條件與材料性質已知，可直接解 出內部溫度場，亦即以邊界條件求取 內部溫度分佈，其解稱為直接解。若 其邊界條件、初始條件或材料性質未 知，而已知內部溫度場之分佈，即由 量測之溫度反求邊界條件、初始條件 或材料性質，其解稱為逆解。 傳統上，求解逆運算問題包括兩步 驟：分析過程與最佳化過程。在分析 過程中，先假設未知的狀態。而後， 以有限差分法或有限元素法或其他數 值方法直接解出分析的結果。將上述 的結果與測量資料結合，產生一組非 線性平方項而對此平方項進行最小化 的過程 [1-11] 。 上述所使用之傳統方法需要反覆之 疊代運算以及在非線性之領域中求解 (含直接解、敏感度分析與最佳化求 解)。造成估算上時間的花費與準確度 之失真。在過去兩年中，已有四種新 方法可解決傳統方法之缺點---逆向矩 陣法 [12]，符號運算法 [13] ，與結合 符號運算之循序方法[14] 與結合數 值運算之循序方法[15]。此四方法之 特色在於當連續之領域經離散化後， 並不需對題目做任何簡化即可求解。 透過逆運算數學模式，求解逆運算問 題之過程可轉換為一組線性方程式來 估算未知之物理量。可直接求解無需 疊代且僅需求解一組線性方程式來估 算未知之物理量。 然而各種方法均有其存在的優缺 點。例如全域性的方法只要少數的幾 個參數即可描繪出大部分的未定函 數，但是卻需要較多的記憶空間與計 算時間。在全時域方法方面，當使用 逆向矩陣法處理多維暫態問題時如果 時間區間數目取的太多，將造成系統 矩陣過大，無法有效率的求出反矩， 導至使用上之不便。例如二維問題在 空間座標上共 100 個離散點以及時間 座標上取 100 個區間，則系統矩陣為 10000*10000。以目前之電腦運算速度 需耗費相當長的時間來計算，不切實 際。循序方法可節省記憶空間與計算 時間但是每個不同的空間與時間點都 需要一個獨立的變數。若以符號運算 之循序方法處理多維暫態問題時，則 面臨低效率之運算[16]。 本研究所討論的方法，是直接解 出未知條件，不須經由多次的疊代逼 近精確值。提高暫態熱傳逆問題的運 算效率。研究中以求解二維逆熱運算 問題為例，說明方法的使用程序。 問題描述 考 慮 二 維 物 體 內 有 許 多 熱 源 為 例。假設物體內部定義為 V ,物體邊界 定 義 為  有 熱 源 的 部 份 定 義 為 S  SQi i1 n_Q ,其數學表示式可寫為:  x k T x   _y kT_y      Qi(t) i1 nQ



(Vi S Qi₎CT t (x, y)V (1) T (x, y, 0) T₀ (x, y)  V t  0 (2) q_n(x, y, t) 0 (x, y) t  0 (3) SQi表示有熱源作用的位置，nq表示熱 源作用的位置數目，Q_i(t) 是作用在S_Qi 處 的 熱 源 大 小 ，(V_i SQ_i ) 是 Dirac delta 函數，qn是在邊界的熱傳量，T表 示溫度分佈T (x, y, t)，k 表示熱傳導係 數與C 是單位容積之熱容量係數。 逆向矩陣法 方程式(1)-(3)經由數值化簡後可得到 AT  b (4) 其中 A 為熱傳矩陣 T 為溫度場 b 為邊 界條件、初始條件與熱源所構成的向

量 AT  S (5) 因此可得到溫度 T  A1S =C  (6) 再將有量測點的溫度場進一步分離得 到 Tmeas Tnull    _= _DD null    _ (7) Tmeas =D 據此得到估算的物理量   ( DT D)1DTTmeas (8) 符號運算法 求解步驟中以一個符號(a_m )來表示目 前待估測之物理量，採用 ADI 之分析 方 法 求 出 以 符 號 表 示 之 溫 度T_{i, j,k} (Function of a_m )，再與量測溫度(T_measi , j,k) 比 較 ， 以 構 成 一 線 性 方 程 式 (T_{i, j,k}=T_measi , j,k)，循序解出待估測之物理 量，再予以穩定化。本方法可直接解 出未知條件，不須經由多次的疊代逼 近精確值。除了上述之優點外， 本方 法能將未知條件假設可以不假設任何 的形式直接由計算機算出對應數值。 這對於處理較複雜形式的未知條件， 或是未知條件的形式難以預測時，較 能準確的描述其結果。同時本方法一 次僅處理一個符號，可避免大矩陣之 運算以及減少大量電腦記憶容量，進 而快速求得逆運算之解。求解步驟可 圖示如下:

Symbolic ADI Solver

f (t_m) a_m Let Ti, j,k Tmeas i, j ,k m m 1 a_m Solve Input boundary condition at t t_m 圖.: 逆運算之求解步驟 符號循序法與數值循序法 研究中使用線性有限元素將空間 座標離散化，並用有限差分法對時間 座標 離散 化。假設 空間有n_p個 離散 點，當t  t_j時，式(1-3)可被表為下列 矩陣方程：

 

B

   

Tj  Sj 

 

A

 

Tj (9) [ A]是具np維的熱矩陣 [B]是n_p維的暫態矩陣 {S_j}是熱源向量 {T_j} 是溫度向量

 

 

1 1 1 } { _     j j j T t T t T (10) 在此t 是時間增量 將(10)式代入(9)式中可得到下式 [K]{T_j} 1 t[B]{Tj1}{Sj} (11) 當t t_j時，由(11)式可推得溫度分佈 {Tj}如下： {T_j} 1 t[K] 1 [B]{T_j1}[K]1{S_j}  [C]{Tj1} [D]{Sj} 在此[C] 1 t[K] 1 [B]與[D] [K]1 在 t  t_m,t_m1,...,t_{m r 1}時溫度分佈可表 示如下：

{Tm} [C]{Tm1} [D]{Sm} {T_m1} [C]{T_m} [D]{S_m1} =[C]2{Tm1}[C][D]{Sm}[D]{Sm1} {T_{mr 1}} [C]{T_{mr 2}} [D]{S_{m r 1}}  [C]r {Tm1} [C] r1 [D]{Sm} [C]r2[D]{S_m1}... [C][D]{S_{mr 2}} [D]{S_{m r 1}} (12) 因此可得到溫度向量在(m + n )個時間 離散點。 {Tmn} [C] n1 {Tm1} [C] l l0 n



[D]{Sm nl} (13) n 為整數且n =0, 1, 2,..., r1 在一單位列向量

 

ui 乘到(13)式之左 右兩邊後，在第i 個空間離散點的溫度 可表示如下： Tm n i _{ u}i  [C]n1{Tm1} u i  [C]l l 0 n



[D]{Sm n l}  ui

 

[C]n1{T_m1} ui

 

[C]l[D] uiQ j

 

m n l iQ j j1 nQ



l 0 n



(14) 在此_{S mn l} u iQ j

 

_miQ_{ nl} j j1 nQ



在此

 

ui 是一單位列向量且在i 位置有 一單位的數值，而 i 即表示量測點。 {uiQ j }是一單位行向量，在iQ j 處有值。同 時_miQ_{ nl} j 表 示 在 i_Qj 處 的 熱 源 強 度 大 小。n_Q表示未知熱源的數量。因此溫 度在m + n 的時間離散點與i 位置離散 點可表示為： Tm n i _ m n i _  mn,m l i ,i_Qj m nl i_Qj j1 n_Q



l 0 n



(15) 在此 m n i _{ u}i

 

[C]n1{Tm1} _mi,i_{ n,m l}Q j  ui

 

[C]l[D] uiQ j

 

將n 由 0 遞增到r1可得： T_mi _mi  _{m, m}i ,iQ j m i_Qj j1 n_Q



Tm1 i _ m1 i _{ } m1, m i,iQ j m1 iQ j m1,m1 i ,iQ j m iQ j





j1 n_Q



(16) 加入未來時間後式(16)可改寫成： T_mi _mi  e₀i ,iQ j m i_Qj j1 n_Q



T_mi_1 _mi_1 e₀i,iQ j  e₁i, i_Qj





j1 n_Q



m i_Qj T_mi_{ 2}_mi_{ 2} (e₀i,iQ j  e₁i ,i_Qj  e₂i,i_Qj ) j1 n_Q



m i_Qj ... T_mi r1m r1 i _ (e₀i, iQ j  e1 i,iQ j  e2 i, iQ j  er1 i ,iQ j ) j1 nQ



m iQ j (17) 當量測溫度Y_ji (t  t_j 與 x  x_i)代入 時， 中之元素可由線性最小均方根 求得如下結果： ?  (T_)1_T _ (18) 符號循序運算法與數值循序運算法所 選用的計算邏輯幾乎是一樣的，因此 在此一並說明。使用符號運算法時需 要較大的記憶體但是由於是符號表示 因此許多算式可由電腦推衍。數值運 算法時可節省記憶空間有較快的執行 速度。兩者都使用未來時間的觀念來 穩定化估算結果。。 工程軟體耦合 本計畫選用 ANSYS 為有限元素 法分析之核心，透過 APDL 可程式語 言與本文所發展之逆運算法結合。其 過程與步驟如下：1.ANSYS 建立幾何 模組 2.建立有限元數模組 3.建立物理 性質模組 4.以 APDL 語言撰寫逆運算 邏輯 5.呼叫 ANSYS 所建立的模組 6. 求解。今從一散熱鰭片之外部估算內 未知熱傳量與溫度為例。

Figure 2: The boundary conditions, the locations of the sensor, and the

Figure 3: A mesh configuration with 868-node and 682-element for a finned heat sink

Figure 4: The estimation of the heat flux in example when  0 and r 1

Figure 5: The estimation of the heat flux in example when  0.5 and

4 

r and  1 and r 4

Figure 6: The estimation of

temperature at the base of the sink in example one when r 1 and

, 3 , 1 , 0   and 5 計畫結果自評 本研究所討論之逆運算模式包含 逆向矩陣法，符號運算法與符號循序 方法與數值循序方法。計畫中整理四 種方法的運算邏輯並對其運算性質提 出說明。並以散熱鰭片之外部估算內 未知熱傳量與溫度為例，說明逆運算 法與工程軟體的耦合。 參考文獻

1. Stolz, G. Tr.,"Numerical Solution to An Inverse Problem of Heat Condition for Simple Shapes", ASME Journal of Heat Transfer, Vol.82, pp.20-26, Feb. 1960

2. Beck, J.V., "Surface Heat Flux Determination Using an Integral Method," Nucl. Eng. Des., Vol.7, pp.170-178, 1968.

3. Beck, J.V., "Nonlinear Estimation Applied to the Nonlinear Inverse Heat Conduction Problem," Intermational Journal of Heat and Mass Transfer, Vol.13, pp.713-71, 2970.

4. Beck, J.V., Litkouhi, B., and St. Clair, C.R., "Efficient Sequential Solution of Nonlinear Inverse Heat Conduction Problem," Numerical Heat Transfer, Vol.5, pp.275-286, 1982. 0 2000 4000 6000 8000 1 104 0 800 1600 2400 3200 Exact  = 0.5  = 1 Heat Fl u x ( W/ m 2) Temporal-coordinate (seconds) 0 2000 4000 6000 8000 1 104 0 800 1600 2400 3200 Exact  = 0 Heat Fl u x ( W/ m 2) Temporal-coordinate (seconds) 20 40 60 80 100 120 140 160 0 800 1600 2400 3200  = 0  = 1  = 3  = 5 Esti mat e d Te mper at ure ( 0C) Temporal-coordinate (seconds)

5. Zienkiewicz, O.C., "The Finite Element Method in Engineering Science", McGraw-Hill, London, 1971.

6. Tohio Yoshimura, and Kazushige Ikuta, "Inverse Heat-Conduction Problem by Finite-Element Formulation", Int.J. System Sci., Vol.16, pp.1365-1376, 1985.

7. Zabaras, N., and Jashua C. Liu, "An Analysis of Two-Dimensional Linear Inverse Heat Transfer Problem Using an Intrgral Method", Numerical Heat Transfer, Vol.13,pp.527-533,1988.

8. Bass, B.R., "Application of the Finite Element Method to the Nonlinear Inverse Heat Conduction Problem Using Beck's Second Method," ASME Journal of Engineering for Industry, Vol.102, pp.168-176, 1980.

9. O.M. Alifanov, "Solution of an Inverse Problem of Heat Conduction by Iteration Method", Journal of Engineering Physics, Vol.25, pp.471-476, 1974.

10. Tikhonov, A.N. and Arsenin, V.Y., Solutions of Ill-Posed Problems, V.H. Winston & Sons, Washington,D.C.,1977.

11. Marquardt, D.W., Generalized Inverses, Ridge Regression, Biased Linear Estimated, and Nonlinear Estimation, Technometrics, 12, pp591-612, 1970.

12. Yang, Ching-yu, "Noniterative Solution of Inverse Heat Conduction Problems in One Dimension," Communications in Numerical Methods in Engineering, Vol. 13, pp. 419-427, 1997.

13. Yang, Ching-yu, A Symbolic Approach to Estimate Two-Sided Boundary Conditions in Two- Dimensional Conduction Problems, AIAA Journal of Thermophysics and Heat Transfer, Vol. 11, No. 3, at press, 1997.

14. Yang, Ching-yu., 1998, “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, pp. 2245-2252.

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

16. Noor, A. K. and Andersen, C. M., "Computerized Symbolic Manipu-lation in Structural Mechanics of Progress and Potential," Computers and Structures, Vol 10, 1979, pp 95-118.