王婷婷，王發(fā)杰，張耀明

(山東理工大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院， 山東 淄博 255049)

數(shù)學(xué)物理反問(wèn)題源于物理、生物、醫(yī)學(xué)、地質(zhì)等眾多科學(xué)領(lǐng)域中的實(shí)際問(wèn)題。反問(wèn)題是相對(duì)于正問(wèn)題而言，主要包括Chauchy問(wèn)題、源項(xiàng)識(shí)別問(wèn)題、參數(shù)識(shí)別問(wèn)題、初值反問(wèn)題、幾何形狀反問(wèn)題和自由邊界反問(wèn)題等。在實(shí)際工程應(yīng)用中，由于受到工況條件或設(shè)備條件的限制，僅能測(cè)量部分結(jié)構(gòu)的邊界溫度和熱流，其余邊界上的邊界條件無(wú)法獲取，這類反問(wèn)題就是所謂的Chauchy熱傳導(dǎo)反問(wèn)題[1]。Chauchy熱傳導(dǎo)反問(wèn)題由于可利用邊界信息較少，傳統(tǒng)數(shù)值方法大多對(duì)測(cè)量邊界誤差非常敏感，通常會(huì)產(chǎn)生不適定問(wèn)題，從而難以得到可接受的數(shù)值解。近年來(lái)，諸多研究者開(kāi)始關(guān)注和探究這類反問(wèn)題的數(shù)值方法。

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

本文是平均源邊界點(diǎn)法在模擬平面熱傳導(dǎo)Cauchy反問(wèn)題的第一次嘗試。眾所周知，數(shù)值模擬方法所建立的不適定系統(tǒng)往往高度病態(tài)，微小的邊界數(shù)據(jù)擾動(dòng)可能會(huì)導(dǎo)致極大的計(jì)算誤差。無(wú)網(wǎng)格法已廣泛應(yīng)用于Helmholtz方程問(wèn)題[7]、含源項(xiàng)熱傳導(dǎo)方程[8]等方面。本文利用新的無(wú)網(wǎng)格方法—平均源邊界點(diǎn)法求解平面二維熱傳導(dǎo)Cauchy反問(wèn)題，采用截?cái)嗥娈愔捣纸?TSVD)和Tikhonov正則化技術(shù)[9-10]求解病態(tài)線性方程組，通過(guò)廣義交叉校驗(yàn)準(zhǔn)則(GCV)來(lái)確定正則化參數(shù)[11-12]。所提出的方法無(wú)需網(wǎng)格劃分，避免了奇異積分計(jì)算，是一種簡(jiǎn)單精確的無(wú)網(wǎng)格方法。本文為平面熱傳導(dǎo)反問(wèn)題的研究開(kāi)辟了新的途徑，拓展了平均源邊界點(diǎn)法的應(yīng)用領(lǐng)域。數(shù)值試驗(yàn)結(jié)果表明了所提方法的有效性、穩(wěn)定性和精確性。

1 平面熱傳導(dǎo)反問(wèn)題的平均源邊界點(diǎn)法

1.1 控制方程

本文假定問(wèn)題的區(qū)域?yàn)橛薪鐓^(qū)域Ω?R2，其邊界為Γ=?Ω?？偸羌俣ㄟ吔缬蓛刹糠纸M成，Γ=Γ1∪Γ2，這里Γ1,Γ2≠?，Γ1∩Γ2=?，函數(shù)u(x)滿足Laplace控制方程：

(1)

邊界條件為：

(2)

(x1,x2)∈Γ1

(3)

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

(4)

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

1.2 平均源邊界點(diǎn)法

為了避免直接計(jì)算強(qiáng)、弱奇異積分，本文方法基于文獻(xiàn)[3-4]的間接平均源邊界積分方程：

(5)

其中：Gij、Hij為系數(shù)矩陣；N為總的邊界節(jié)點(diǎn)數(shù)；uj、qj是節(jié)點(diǎn)j處的溫度和法向通量。Gij、Hij表示為：

(6)

(7)

遠(yuǎn)離邊界的內(nèi)部點(diǎn)y的溫度和溫度梯度可表示為：

(8)

(9)

(10)

其中rk=xk-yk。

2 正則化方法

2.1 截?cái)嗥娈愔捣纸?TSVD)

考慮如下線性方程組：

Ax=b

(11)

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

(12)

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

截?cái)嗥娈愔捣纸?TSVD)是一種常用的正則化方法，其基本思想[9]是用K階矩陣AK來(lái)逼近M×N階矩陣A。其中，AK可以表示為

(13)

這里K為正則化參數(shù)，與之對(duì)應(yīng)的截?cái)嗥娈愔禐?/p>

(14)

2.2 Tikhonov正則化方法

Tikhonov方法是一種應(yīng)用比較廣泛且能有效地解決病態(tài)問(wèn)題的正則化方法，其基本思想[10]如下：

把正則化泛函

(15)

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

(16)

2.3 廣義交叉校驗(yàn)準(zhǔn)則

廣義交叉校驗(yàn)準(zhǔn)則(GCV)由Golub[12]提出，該方法以正則化參數(shù)K為參變量，尋求GCV函數(shù)的最小值。當(dāng)求得GCV的極小值時(shí)，對(duì)應(yīng)的K值為最優(yōu)值。G(K)的計(jì)算公式為

(17)

其中AI滿足：

可用來(lái)產(chǎn)生正則化參數(shù)。

3 數(shù)值算例

在實(shí)際的工程計(jì)算中，被測(cè)量的邊界值會(huì)不可避免地受擾動(dòng)的影響。本文的數(shù)值算例中采用以下公式[13]來(lái)增加邊界數(shù)據(jù)擾動(dòng)：

(18)

其中：bi是精確解；rand是一個(gè)-1～1的隨機(jī)數(shù)；δ表示擾動(dòng)幅度。為了評(píng)估數(shù)值解的有效性，引入下列平均相對(duì)誤差[13]：

Average relative error=

(19)

3.1 算例1：齒輪型區(qū)域

首先，考慮齒輪型區(qū)域中的穩(wěn)定溫度場(chǎng)問(wèn)題。該區(qū)域由兩部分組成：內(nèi)部邊界Γ2為半徑r=0.5的圓域；外部邊界Γ1的參數(shù)方程為

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

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

u(x1,x2)=ex1cosx2

圖1 齒輪型區(qū)域

計(jì)算時(shí)，內(nèi)、外邊界分別劃分為60、120個(gè)節(jié)點(diǎn)。分別使用TSVD和Tikhonov正則化方法對(duì)問(wèn)題進(jìn)行求解，并用GCV法選取正則化參數(shù)。如圖2、3所示，兩種方法分別在k=126和λ=0.01處取得最優(yōu)正則化參數(shù)。

圖2 GCV法參數(shù)選取

圖3 未知邊界上數(shù)值解與解析解對(duì)比

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

3.2 算例2：三圓環(huán)區(qū)域

考慮三圓環(huán)區(qū)域的穩(wěn)定溫度場(chǎng)問(wèn)題，多連通區(qū)域如圖4所示。內(nèi)邊界Γ2和Γ3的溫度和通量都不可獲得，但外邊界Γ1的溫度和通量已知且溫度場(chǎng)的精確解為

u=r2cos(2θ)+rsinθ

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

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

圖5 未知邊界上溫度和通量在無(wú)擾動(dòng)及添加3%擾動(dòng)下數(shù)值解與解析解對(duì)比

正則化方法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：四圓環(huán)區(qū)域

考慮四圓環(huán)區(qū)域上的穩(wěn)態(tài)溫度場(chǎng)問(wèn)題，如圖6所示。外邊界Γ1的邊界條件已知且溫度場(chǎng)的精確解為

其中內(nèi)邊界的溫度和通量未知。假設(shè)邊界節(jié)點(diǎn)在這4部分邊界上的分布為Γ1∶Γ2∶Γ3∶Γ4=7∶1∶1∶1。

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

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

表2 不同節(jié)點(diǎn)數(shù)下正則化參數(shù)的選取

圖7 不同節(jié)點(diǎn)數(shù)下未知邊界溫度和通量的平均相對(duì)誤差

4 結(jié)束語(yǔ)

本文首次采用平均源邊界點(diǎn)法(ASBNM)求解二維Cauchy熱傳導(dǎo)反問(wèn)題。對(duì)求解中產(chǎn)生的不適定線性系統(tǒng)，采用TSVD和Tikhonov技術(shù)結(jié)合GCV的正則化方法來(lái)求解。與其他現(xiàn)有的求解病態(tài)的Cauchy熱傳導(dǎo)反問(wèn)題的方法相比，該方法無(wú)積分計(jì)算、無(wú)網(wǎng)格，具有計(jì)算精度高、收斂速度快、程序?qū)崿F(xiàn)簡(jiǎn)單、適合于高維問(wèn)題等優(yōu)點(diǎn)，是一種真正意義上的無(wú)網(wǎng)格數(shù)值方法。

通過(guò)3個(gè)經(jīng)典的數(shù)值算例測(cè)試了復(fù)雜多聯(lián)通區(qū)域上的熱傳導(dǎo)反問(wèn)題，數(shù)值結(jié)果表明了方法的精度和有效性。本文為Cauchy熱傳導(dǎo)反問(wèn)題求解提供了新的求解思路，同時(shí)拓寬了平均源邊界點(diǎn)法的應(yīng)用范圍。

[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] 張耀明,屈文鎮(zhèn),陳正宗.三維位勢(shì)問(wèn)題新的規(guī)則化邊界元法[J].中國(guó)科學(xué)(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] 師晉紅,傅卓佳,陳文.用邊界粒子法求解柯西反問(wèn)題的數(shù)值計(jì)算軟件包[J].計(jì)算機(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.