一種求解平面熱傳導反問題的新型無網格方法

王婷婷，王發杰，張耀明

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

數學物理反問題源于物理、生物、醫學、地質等眾多科學領域中的實際問題。反問題是相對于正問題而言，主要包括Chauchy問題、源項識別問題、參數識別問題、初值反問題、幾何形狀反問題和自由邊界反問題等。在實際工程應用中，由于受到工況條件或設備條件的限制，僅能測量部分結構的邊界溫度和熱流，其余邊界上的邊界條件無法獲取，這類反問題就是所謂的Chauchy熱傳導反問題[1]。Chauchy熱傳導反問題由于可利用邊界信息較少，傳統數值方法大多對測量邊界誤差非常敏感，通常會產生不適定問題，從而難以得到可接受的數值解。近年來，諸多研究者開始關注和探究這類反問題的數值方法。

有限元法(FEM)和有限差分法(FDM)長久以來是實際工程應用中的主流數值計算方法，已被廣泛應用于反問題的研究中。但是，這些基于網格類的方法需要對整個計算區域進行離散，過程繁瑣且耗時，特別是對于無限域問題和高維的不規則幾何形狀問題。此外，反問題已知部分邊界條件時，有限元模擬相當困難且計算精度不高。邊界元法(BEM)將問題維數降低一維，并具有邊界離散和半解析性質，是目前較有潛質的一種數值算法[2]，在熱傳導反問題的求解中，比有限元法更有效、更準確。然而，邊界元法不可避免地涉及奇異積分和幾乎奇異積分的計算，這就為解決實際問題帶來了很大困擾。針對此情況，筆者[3-4]提出了一種簡單精確的無網格方法—平均源邊界點法。該方法基于規則化邊界積分方程[5-6]，通過加減去奇異和平均積分的思想，消除了基本解的源點奇異性，具有無網格、無積分、僅需邊界離散、半解析的特性，目前已成功應用于位勢問題的研究。

本文是平均源邊界點法在模擬平面熱傳導Cauchy反問題的第一次嘗試。眾所周知，數值模擬方法所建立的不適定系統往往高度病態，微小的邊界數據擾動可能會導致極大的計算誤差。無網格法已廣泛應用于Helmholtz方程問題[7]、含源項熱傳導方程[8]等方面。本文利用新的無網格方法—平均源邊界點法求解平面二維熱傳導Cauchy反問題，采用截斷奇異值分解(TSVD)和Tikhonov正則化技術[9-10]求解病態線性方程組，通過廣義交叉校驗準則(GCV)來確定正則化參數[11-12]。所提出的方法無需網格劃分，避免了奇異積分計算，是一種簡單精確的無網格方法。本文為平面熱傳導反問題的研究開辟了新的途徑，拓展了平均源邊界點法的應用領域。數值試驗結果表明了所提方法的有效性、穩定性和精確性。

1 平面熱傳導反問題的平均源邊界點法

1.1 控制方程

本文假定問題的區域為有界區域Ω?R2，其邊界為Γ=?Ω?？偸羌俣ㄟ吔缬蓛刹糠纸M成，Γ=Γ1∪Γ2，這里Γ1,Γ2≠?，Γ1∩Γ2=?，函數u(x)滿足Laplace控制方程：

(1)

邊界條件為：

(2)

(x1,x2)∈Γ1

(3)

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

(4)

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

1.2 平均源邊界點法

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

(5)

其中：Gij、Hij為系數矩陣；N為總的邊界節點數；uj、qj是節點j處的溫度和法向通量。Gij、Hij表示為：

(6)

(7)

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

(8)

(9)

(10)

其中rk=xk-yk。

2 正則化方法

2.1 截斷奇異值分解(TSVD)

考慮如下線性方程組：

Ax=b

(11)

其中：A∈Rm×n；x∈Rn；b∈Rm，且m≥n。對矩陣A進行奇異值分解，有

(12)

其中：U=(u1,…,uM)和V=(v1,…,vN)分別滿足UTU=IM，VTV=IN；∑=diag(σ1,…,σN)表示非負對角矩陣。

截斷奇異值分解(TSVD)是一種常用的正則化方法，其基本思想[9]是用K階矩陣AK來逼近M×N階矩陣A。其中，AK可以表示為

(13)

這里K為正則化參數，與之對應的截斷奇異值為

(14)

2.2 Tikhonov正則化方法

Tikhonov方法是一種應用比較廣泛且能有效地解決病態問題的正則化方法，其基本思想[10]如下：

把正則化泛函

(15)

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

(16)

2.3 廣義交叉校驗準則

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

(17)

其中AI滿足：

可用來產生正則化參數。

3 數值算例

在實際的工程計算中，被測量的邊界值會不可避免地受擾動的影響。本文的數值算例中采用以下公式[13]來增加邊界數據擾動：

(18)

其中：bi是精確解；rand是一個-1～1的隨機數；δ表示擾動幅度。為了評估數值解的有效性，引入下列平均相對誤差[13]：

Average relative error=

(19)

3.1 算例1：齒輪型區域

首先，考慮齒輪型區域中的穩定溫度場問題。該區域由兩部分組成：內部邊界Γ2為半徑r=0.5的圓域；外部邊界Γ1的參數方程為

2(n+1)cos(nφ-nπ/2)], 0≤φ≤2π}

其中n=6。邊界Γ2上的溫度和通量未知，但邊界Γ1上的溫度和通量均可獲得，問題的解析解為

u(x1,x2)=ex1cosx2

圖1 齒輪型區域

計算時，內、外邊界分別劃分為60、120個節點。分別使用TSVD和Tikhonov正則化方法對問題進行求解，并用GCV法選取正則化參數。如圖2、3所示，兩種方法分別在k=126和λ=0.01處取得最優正則化參數。

圖2 GCV法參數選取

圖3 未知邊界上數值解與解析解對比

圖4描述了兩種正則化方法下求得的未知邊界的溫度和通量的數值解與解析解的比較，從圖中可以看出，所求結果與解析解十分吻合。

3.2 算例2：三圓環區域

考慮三圓環區域的穩定溫度場問題，多連通區域如圖4所示。內邊界Γ2和Γ3的溫度和通量都不可獲得，但外邊界Γ1的溫度和通量已知且溫度場的精確解為

u=r2cos(2θ)+rsinθ

圖4 多連通區域(r1=2,r2=0.5,r3=0.25,a=1.0)

圖5描繪了兩種正則化方法的溫度、通量在無擾動和3%擾動程度下的精確解與本文計算結果的比較。表1給出了不同擾動下GCV法選取的正則化參數，發現在k=130和λ=1.00×10-3處可獲得最優參數。從圖中可以看出，在已知條件存在隨機擾動的情況下，邊界溫度和通量與精確解能較好地吻合，兩種方法均能求得較好的結果。

圖5 未知邊界上溫度和通量在無擾動及添加3%擾動下數值解與解析解對比

正則化方法0%1%3%5%TSVD法130130130130Tikhonov法4.040 3×10-54.040 3×10-53.418 2×10-53.418 2×10-5

3.3 算例3：四圓環區域

考慮四圓環區域上的穩態溫度場問題，如圖6所示。外邊界Γ1的邊界條件已知且溫度場的精確解為

其中內邊界的溫度和通量未知。假設邊界節點在這4部分邊界上的分布為Γ1∶Γ2∶Γ3∶Γ4=7∶1∶1∶1。

圖6 四圓型區域(r1=2,r2=r3=r4=0.25,a=1.0)

圖7給出了未知邊界溫度和通量的平均相對誤差曲線。分別將邊界劃分成100、200、300、400、500、600個節點，表2給出了不同節點數時 GCV法和TR法選擇的最優正則化參數。從圖7可以看出：隨著單元數的增加，邊界溫度、通量的平均相對誤差呈現明顯下降的趨勢，表明該方法具有較好的收斂性；僅使用100個邊界節點所求得的溫度及通量的相對誤差在10-2～10-1，說明該方法對位勢反問題可進行有效的計算求解。

表2 不同節點數下正則化參數的選取

圖7 不同節點數下未知邊界溫度和通量的平均相對誤差

4 結束語

本文首次采用平均源邊界點法(ASBNM)求解二維Cauchy熱傳導反問題。對求解中產生的不適定線性系統，采用TSVD和Tikhonov技術結合GCV的正則化方法來求解。與其他現有的求解病態的Cauchy熱傳導反問題的方法相比，該方法無積分計算、無網格，具有計算精度高、收斂速度快、程序實現簡單、適合于高維問題等優點，是一種真正意義上的無網格數值方法。

通過3個經典的數值算例測試了復雜多聯通區域上的熱傳導反問題，數值結果表明了方法的精度和有效性。本文為Cauchy熱傳導反問題求解提供了新的求解思路，同時拓寬了平均源邊界點法的應用范圍。

[1] WANG Fajie,CHEN Wen,QU Wenzhen,et al.A BEM formulation in conjunction with parametric equation approach for three-dimensional Cauchy problems of steady heat conduction [J].Engineering Analysis with Boundary Elements,2016,63(114):1-14.

[2] CHENG A H D,CHENG D T.Heritage and early history of the boundary element method [J].Engeering Analysis with Boundary Elements,2005,29(3):268-302.

[3] ZHANG Yaoming,SUN Fangling,YOUNG Derliang,et al.Average source boundary node method for potential problems [J].Engeering Analysis with Boundary Elements,2016,70:114-125.

[4] SUN Fangling,ZHANG Yaoming,YOUNG Derliang,et al.A new boundary meshfree method for potential problems [J].Advances in Engineering Software,2016,100(C):32-42.

[5] ZHANG Yaoming,SLADEK V,SLADEK J,et al.A new boundary integral equation formulation for plane orthotropic elastic media [J].Applied Mathematical Modelling,2012,36(10):4862-4875.

[6] 張耀明,屈文鎮,陳正宗.三維位勢問題新的規則化邊界元法[J].中國科學(G輯),2013,43(3):297-308.

[7] CHEN Wen,FU Zhuojia.Boundary Particle Method for Inverse Cauchy Problems of Inhomogeneous Helmholtz Equations [J].Journal of Marine Science and Technology-Taiwan,2009,17(3):157-163.

[8] FU Zhuojia,CHEN Wen,ZHANG Chuanzeng.Boundary Particle Method for Inverse Cauchy inhomogeneous potential problems [J].Inverse Problems in Science and Engineering,2012,20(2):189-207.

[9] HANSEN P C,SEKII T,SHIBAHASHI H.The modified truncated SVD method for regularization in general form [J].Society for Industrial and Applied Mathematics,1992,13(5):1142-1150.

[10] SCHERZER O.The use of Morozov’s discrepancy principle for Tikhonov regularization for solving nonlinear ill-posed problems [J].Computing,1993,51(1):45-60.

[11] 師晉紅,傅卓佳,陳文.用邊界粒子法求解柯西反問題的數值計算軟件包[J].計算機輔助工程,2013,22(1):64-70.

[12] GOLUB G H,HEATH M,WAHBA G.Generalized cross-validation function for ill-conditioned least squares problems [J].Technometrics,1979,21(2):215-223.

[13] GU Yan,CHEN Wen,FU Zhuojia.Singular boundary method for inverse heat conduction problems in general anisotropic media [J].Inverse Problems in Science & Engineering,2014,22(6):889-909.