三維熱傳導邊界條件識別反問題的平均源邊界節點法研究

劉曌陽，王婷婷，王發杰，張耀明

(山東理工大學 數學與統計學院， 山東 淄博 255049)

近年來，許多學者對熱傳導反問題進行了深入研究[1]，并廣泛應用于冶煉、航空、化學工程等方面.

邊界元法[2-3]將問題的維數降低一維，且只需對物體邊界劃分單元，既降低了工作難度，又提高了計算精度，是目前很有潛質的一種數值算法，已廣泛應用于各種工程問題.然而，對于復雜幾何邊界的問題，劃分單元仍然是非常耗費精力的.此外，邊界元法需要計算大量復雜邊界積分，這些積分甚至是奇異的和幾乎奇異的，這就給解決實際問題帶來了很大的困難.張耀明提出了一種新的邊界型無網格方法——平均源邊界節點法[4-5](Average Source Boundary Node Method, ASBNM)，該方法基于完全規則化邊界積分方程和平均源技術[6]，具有無網格、無積分、半解析的特征.

本文是平均源邊界節點法在模擬三維熱傳導反問題的一次嘗試.采用截斷奇異值分解(truncated singular value decomposition TSVD)和Tikhonov正則化技術[7-8]求解病態線性方程組，利用廣義交叉校驗準則(generalized cross validation criterion GCV)來確定正則化參數[9].

1 三維熱傳導反問題的平均源邊界節點法

1.1 控制方程

本文假定問題的區域為有界區域Ω?R3，其邊界為Γ=?Ω.邊界Γ由兩部分組成，即Γ=Γ1∪Γ2，這里Γ1∩Γ2=?，函數u(x)滿足Laplace控制方程：

(1)

邊界條件為

(x1,x2,x3)∈Γ1

(2)

(3)

控制方程(1)的基本解為

(4)

式中|x-y|表示兩點x和y之間的歐幾里得距離.

1.2 平均源邊界節點法

為了避免直接計算強、弱奇異積分，本文基于文獻[4-5]的平均源邊界節點法：

(5)

式中：Gij,Hij為系數矩陣；N為總的邊界節點數；uj,qj是節點j處的溫度和法向通量.

Gij,Hij如下：

(6)

(7)

遠離邊界的內部點y溫度和溫度梯度表示為

(8)

(9)

k=1,2,3

(10)

式中rk=xk-yk.

2 正則化方法

2.1 截斷奇異值分解

考慮如下線性方程組

Ax=b

(11)

式中A∈Rm×n，x∈Rn，b∈Rm，且m≥n.

截斷奇異值分解的基本思想[7]是用K階矩陣AK來逼近M×N階矩陣A.其中，AK可以表示為

(12)

這里U=(u1,…,uM)和V=(v1,…,vN)分別滿足UTU=IM;VTV=IN；Σ=diag(σ1,…σN)表示非負對角矩陣；K為正則化參數，與之對應的截斷奇異值為

(13)

2.2 Tikhonov正則化方法

Tikhonov法的基本思想[8]如下：

把正則化泛函

Jα(x)=‖Ax-b‖2+α2‖x‖2,α>0

(14)

的極小元xα作為方程Ax=b的解，可表示為如下形式

(15)

2.3 廣義交叉校驗準則

廣義交叉校驗準則是由Golub.G.H.[9]提出，該方法以正則化參數K為參變量，當求得GCV的極小值時，對應的K值為最優值.其計算公式為

(16)

其中AI滿足

(17)

可用來產生正則化參數.

3 數值算例

數值算例中采用以下公式[8]來施加邊界數據擾動：

(18)

這里，bi是精確解，rand是一個-1到1之間的隨機數，δ表示擾動幅度.為了對比不同方法數值解的有效性，引入下列平均相對誤差

Average relative error =

(19)

算例1球型區域

u(x1,x2,x3)=x1x2+x1x3+x2x3+

x1+x2+x3+2

(20)

圖1 球型區域Fig.1 Problem sketch

計算時，球域外殼上配置200個邊界節點.分別用TSVD和Tikhonov正則化方法對問題進行求解，并用GCV法選取正則化參數.如圖2所示，兩種方法分別在k=155和l=0.007處取得最優正則化參數.圖3和圖4分別描述了結合TR-GCV、TSVD-GCV法求解未知邊界的溫度和通量的誤差曲面圖.從圖中可以看出，兩種方法均能對三維Cauchy反問題進行有效求解.

(a)TSVD正則化參數 (b)Tikhonov正則化參數圖2 GCV法選取的正則化參數Fig.2 GCV Curvers for TSVD and TR

(a)溫度 (b) 通量圖3 未知邊界上溫度和通量的誤差曲面圖(結合TR-GCV法)Fig.3 Relative error surfaces of the temperature and the flux on the unavailable boundary obtained using the TR-GCV method

(a) 溫度 (b) 通量圖4 未知邊界上溫度和通量的誤差曲面圖(結合TSVD-GCV法)Fig.4 Relative error surfaces of the temperature and the flux on the unavailable boundary obtained using the TSVD-GCV method

為了驗證方法的收斂性，計算時分別在邊界配置200、625、840、960、1 200個邊界節點.表1給出了不同節點數下，結合TSVD和TR技術，GCV法選擇的最優正則化參數.圖5給出了未知邊界溫度和通量的平均相對誤差收斂曲線.從圖中可以看出，隨著邊界節點數的增加，邊界溫度和通量的平均相對誤差呈現明顯下降的趨勢，表明該方法具有很好的收斂性，從而說明此方法對三維Cauchy反問題可以進行有效地求解.

表1 不同單元數下正則化參數的選取

Tab.1 Optimal regularization parameters for different boundary nodes

節點數GCV法參數kTR法參數λ2001550.0076254820.0018406420.0029607340.0011 2009140.001

算例2空心球區域

考慮空心球區域的穩態溫度場問題.內邊界Γ2的溫度和通量都未知，外邊界Γ1的溫度和通量已知，且溫度場精確解為

u(x1,x2,x3)=x1x2x3+10x1+10x2+10x3

(21)

(a) 溫度 (b) 通量圖5 不同邊界節點數下未知邊界溫度和通量的平均相對誤差Fig.5 Average relative error of temperatures and fluxes on underspecified boundary of different boundary nodes using the TSVD technique and the TR technique in conjunction with the GCV criterion

表2給出了GCV法選取的正則化參數，在k=924和l=1.00×10-3處獲得最優參數.利用兩種正則化方法及其相應的最優化參數，圖6給出了在3%擾動程度下溫度和通量的精確解與本文計算結果的比較.從圖中可以看出，在已知條件存在隨機擾動的情況下，邊界溫度和通量與精確解能夠較好地吻合，兩種方法均能求得較好的結果.

表2 不同擾動下正則化參數的選取

Tab.2Optimalregularizationparametersforinvestigatedregularizationmethodatvariousnoiselevels

正則化方法0%1%3%5%TSVD法924924924924Tikhonov法1.00×10-31.00×10-31.00×10-31.00×10-3

(a) 溫度:TSVD-GCV (b) 溫度：TR-GCV

(c)溫度：TSVD-GCV-3%擾動 (d) 溫度：TR-GCV-3%擾動 圖6 未知邊界上溫度的數值解與解析解的對比圖Fig.6 Relative error surfaces of the temperature on the unavailable boundary compared with the analytical solution

4 結束語

本文采用平均源邊界節點法(ASBNM)求解三維Cauchy熱傳導反問題，對求解中產生的不適定線性系統，采用TSVD和Tikhonov技術求解結合GCV的正則化方法來求解.與其他現有的求解病態的Cauchy熱傳導反問題的方法相比，該方法無積分計算、無網格，具有計算精度高、收斂速度快、程序實現簡單，適合于高維問題等優點.為Cauchy熱傳導反問題提供了新的求解思路，且拓寬了平均源邊界節點法的應用范圍.