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

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

近年來，許多學(xué)者對(duì)熱傳導(dǎo)反問題進(jìn)行了深入研究[1]，并廣泛應(yīng)用于冶煉、航空、化學(xué)工程等方面.

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

本文是平均源邊界節(jié)點(diǎn)法在模擬三維熱傳導(dǎo)反問題的一次嘗試.采用截?cái)嗥娈愔捣纸?truncated singular value decomposition TSVD)和Tikhonov正則化技術(shù)[7-8]求解病態(tài)線性方程組，利用廣義交叉校驗(yàn)準(zhǔn)則(generalized cross validation criterion GCV)來確定正則化參數(shù)[9].

1 三維熱傳導(dǎo)反問題的平均源邊界節(jié)點(diǎn)法

1.1 控制方程

本文假定問題的區(qū)域?yàn)橛薪鐓^(qū)域Ω?R3，其邊界為Γ=?Ω.邊界Γ由兩部分組成，即Γ=Γ1∪Γ2，這里Γ1∩Γ2=?，函數(shù)u(x)滿足Laplace控制方程：

(1)

邊界條件為

(x1,x2,x3)∈Γ1

(2)

(3)

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

(4)

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

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

為了避免直接計(jì)算強(qiáng)、弱奇異積分，本文基于文獻(xiàn)[4-5]的平均源邊界節(jié)點(diǎn)法：

(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)

k=1,2,3

(10)

式中rk=xk-yk.

2 正則化方法

2.1 截?cái)嗥娈愔捣纸?/h3>

考慮如下線性方程組

Ax=b

(11)

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

截?cái)嗥娈愔捣纸獾幕舅枷隱7]是用K階矩陣AK來逼近M×N階矩陣A.其中，AK可以表示為

(12)

這里U=(u1,…,uM)和V=(v1,…,vN)分別滿足UTU=IM;VTV=IN；Σ=diag(σ1,…σN)表示非負(fù)對(duì)角矩陣；K為正則化參數(shù)，與之對(duì)應(yīng)的截?cái)嗥娈愔禐?/p>

(13)

2.2 Tikhonov正則化方法

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

把正則化泛函

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

(14)

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

(15)

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

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

(16)

其中AI滿足

(17)

可用來產(chǎn)生正則化參數(shù).

3 數(shù)值算例

數(shù)值算例中采用以下公式[8]來施加邊界數(shù)據(jù)擾動(dòng)：

(18)

這里，bi是精確解，rand是一個(gè)-1到1之間的隨機(jī)數(shù)，δ表示擾動(dòng)幅度.為了對(duì)比不同方法數(shù)值解的有效性，引入下列平均相對(duì)誤差

Average relative error =

(19)

算例1球型區(qū)域

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

x1+x2+x3+2

(20)

圖1 球型區(qū)域Fig.1 Problem sketch

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

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

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

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

表1 不同單元數(shù)下正則化參數(shù)的選取

Tab.1 Optimal regularization parameters for different boundary nodes

節(jié)點(diǎn)數(shù)GCV法參數(shù)kTR法參數(shù)λ2001550.0076254820.0018406420.0029607340.0011 2009140.001

算例2空心球區(qū)域

考慮空心球區(qū)域的穩(wěn)態(tài)溫度場(chǎng)問題.內(nèi)邊界Γ2的溫度和通量都未知，外邊界Γ1的溫度和通量已知，且溫度場(chǎng)精確解為

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

(21)

(a) 溫度 (b) 通量圖5 不同邊界節(jié)點(diǎn)數(shù)下未知邊界溫度和通量的平均相對(duì)誤差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法選取的正則化參數(shù)，在k=924和l=1.00×10-3處獲得最優(yōu)參數(shù).利用兩種正則化方法及其相應(yīng)的最優(yōu)化參數(shù)，圖6給出了在3%擾動(dòng)程度下溫度和通量的精確解與本文計(jì)算結(jié)果的比較.從圖中可以看出，在已知條件存在隨機(jī)擾動(dòng)的情況下，邊界溫度和通量與精確解能夠較好地吻合，兩種方法均能求得較好的結(jié)果.

表2 不同擾動(dòng)下正則化參數(shù)的選取

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òng) (d) 溫度：TR-GCV-3%擾動(dòng) 圖6 未知邊界上溫度的數(shù)值解與解析解的對(duì)比圖Fig.6 Relative error surfaces of the temperature on the unavailable boundary compared with the analytical solution

4 結(jié)束語

本文采用平均源邊界節(jié)點(diǎn)法(ASBNM)求解三維Cauchy熱傳導(dǎo)反問題，對(duì)求解中產(chǎn)生的不適定線性系統(tǒng)，采用TSVD和Tikhonov技術(shù)求解結(jié)合GCV的正則化方法來求解.與其他現(xiàn)有的求解病態(tài)的Cauchy熱傳導(dǎo)反問題的方法相比，該方法無積分計(jì)算、無網(wǎng)格，具有計(jì)算精度高、收斂速度快、程序?qū)崿F(xiàn)簡(jiǎn)單，適合于高維問題等優(yōu)點(diǎn).為Cauchy熱傳導(dǎo)反問題提供了新的求解思路，且拓寬了平均源邊界節(jié)點(diǎn)法的應(yīng)用范圍.