岳俊宏,李 明,牛瑞萍

(太原理工大學(xué) 數(shù)學(xué)學(xué)院,太原 030024)

?

廣義邊界控制法在多層熱傳導(dǎo)邊界識(shí)別問(wèn)題中的應(yīng)用

岳俊宏,李 明,牛瑞萍

(太原理工大學(xué) 數(shù)學(xué)學(xué)院,太原 030024)

基于有限元(FEM)的廣義邊界控制法可應(yīng)用于求解固體力學(xué)的反問(wèn)題。帶有柯西數(shù)據(jù)(部分邊界的溫度值與熱流值)的多層熱傳導(dǎo)邊界識(shí)別問(wèn)題是一類反向熱傳導(dǎo)問(wèn)題。研究用該方法求解帶有柯西數(shù)據(jù)的一維多層熱傳導(dǎo)邊界識(shí)別問(wèn)題,并證明了該方法的可行性。數(shù)值實(shí)例證實(shí)基于有限元的廣義邊界控制法對(duì)多層熱傳導(dǎo)邊界識(shí)別問(wèn)題是有效且穩(wěn)定的。

有限元;廣義邊界控制法;多層熱傳導(dǎo);邊界識(shí)別

熱傳導(dǎo)方程的邊界識(shí)別問(wèn)題是通過(guò)測(cè)量部分邊界或某些內(nèi)部位置上的溫度值和熱流值來(lái)得到未知區(qū)域信息的問(wèn)題,并可根據(jù)所得信息和移動(dòng)邊界所滿足的Dirichlet邊界條件確定該移動(dòng)邊界[1]。這類問(wèn)題有廣泛的物理背景,如在冰山勘探過(guò)程中,希望通過(guò)露在水面上冰山的信息來(lái)探測(cè)出冰山在水中邊界形狀及大小;在高爐煉鐵過(guò)程中,由于高爐壁與融化的鐵水、各類雜質(zhì)長(zhǎng)時(shí)間相互作用使得高爐壁的內(nèi)層厚度出現(xiàn)變化進(jìn)而造成破裂,因此,進(jìn)行實(shí)時(shí)監(jiān)控高爐壁內(nèi)層的厚度變化是必要的。然而,高爐壁(由多層耐高溫、絕熱、高強(qiáng)度的材料組成)內(nèi)層的情況不易獲得,這就需要通過(guò)外部邊界的溫度和熱流值反向求得內(nèi)層邊界的情況,這類問(wèn)題就是反向熱傳導(dǎo)問(wèn)題。和其它反問(wèn)題一樣,它們通常是病態(tài)的,即輸入數(shù)據(jù)的任何小的改變,都會(huì)造成結(jié)果的誤差很大。為了克服這個(gè)問(wèn)題,許多研究人員提出了不同的解決辦法,其中正則化方法是一種使解穩(wěn)定的有效方法[2]。

目前,在穩(wěn)態(tài)熱傳導(dǎo)方程的邊界識(shí)別和重構(gòu)問(wèn)題上已有大量的研究[3-4]。求解這類熱傳導(dǎo)方程的正問(wèn)題已經(jīng)很成熟且有很多種數(shù)值方法,如有限元[5],光滑有限元[6],無(wú)網(wǎng)格法[7]等。然而求解相應(yīng)反問(wèn)題的方法主要是無(wú)網(wǎng)格法,如:基本解法(Method of Fundamental Solution,MFS)[8]、Kansa方法[9]。但無(wú)網(wǎng)格法求解速度很慢,不利于應(yīng)用在工程中。眾所周知,有限元的計(jì)算效率比無(wú)網(wǎng)格法高的多。因此,人們提出了基于有限元的廣義邊界控制法(GBCM)[10],它已經(jīng)應(yīng)用于懸臂梁的反問(wèn)題中。廣義邊界控制法的思想是將反問(wèn)題轉(zhuǎn)化為正問(wèn)題,然后通過(guò)求解相應(yīng)的正問(wèn)題得到反問(wèn)題的解。為了克服數(shù)據(jù)測(cè)量中隨機(jī)擾動(dòng)帶來(lái)的不適定性,該方法使用Tikhonov正則化,并通過(guò)L曲線法選取正則化參數(shù)[2]。事實(shí)上,對(duì)于一維問(wèn)題的廣義邊界控制法并不需要采用正則化。

目前,關(guān)于多層材料瞬態(tài)熱傳導(dǎo)方程的邊界識(shí)別問(wèn)題的研究并不多[11-12],但是這類問(wèn)題的唯一性已經(jīng)被證明[13]。在這些文獻(xiàn)中主要使用基本解法和Kansa方法。但這些方法計(jì)算效率不高,而且結(jié)果受參數(shù)的影響很大,而有限元能比較高效、經(jīng)濟(jì)地模擬各類問(wèn)題,且結(jié)果較穩(wěn)定。因此,我們用新提出的基于有限元的廣義邊界控制法來(lái)求解這類問(wèn)題。假設(shè)多層材料區(qū)域在y方向是無(wú)限長(zhǎng)的,即熱源和任何溫度只與x有關(guān),與y無(wú)關(guān)。在這種情況下,該問(wèn)題在空間上變?yōu)橐痪S瞬態(tài)熱傳導(dǎo)問(wèn)題。

本文首先介紹了一維瞬態(tài)熱傳導(dǎo)方程在多層區(qū)域上的邊界識(shí)別問(wèn)題。其次,給出了基于有限元的廣義邊界控制法的公式,用加權(quán)殘值法將偏微分方程轉(zhuǎn)化為弱形式下的有限元系統(tǒng)方程,進(jìn)而導(dǎo)出一維熱傳導(dǎo)問(wèn)題的廣義邊界控制法的公式,并對(duì)該方法的可行性進(jìn)行了分析。最后,3個(gè)數(shù)值實(shí)例表明該方法對(duì)多層區(qū)域上的一維瞬態(tài)熱傳導(dǎo)邊界識(shí)別問(wèn)題穩(wěn)定且有效。

1 多層區(qū)域上的熱傳導(dǎo)邊界識(shí)別問(wèn)題

本節(jié)考慮具有移動(dòng)邊界s(t)的多層區(qū)域上的一維熱傳導(dǎo)邊界識(shí)別問(wèn)題。為了簡(jiǎn)化問(wèn)題,下面以3層區(qū)域?yàn)槔?見(jiàn)圖1)。對(duì)于更復(fù)雜的情況也可使用類似的方法。

在上述3層區(qū)域上滿足如下熱傳導(dǎo)控制方程:

(1)

(2)

(3)

式中：ui(x,t),i=1,2,3分別是區(qū)域Di,i=1,2,3內(nèi)的溫度分布；T是模型運(yùn)行的最大時(shí)間；ki,ρi,ci,i=1,2,3分別是區(qū)域Di,i=1,2,3的熱傳導(dǎo)系數(shù),密度和比熱。D1=[l1,l2)×[0,T],D2=[l2,l3)×[0,T],D3=[l3,s(t)]×[0,T]分別是控制方程(式1-式3)解所在的區(qū)域。

該熱傳導(dǎo)問(wèn)題滿足如下的交界面條件:

(4)

(5)

(6)

(7)

在固定端x=l1處的邊值條件為：

(8)

(9)

其中，式(8)、式(9)分別是控制方程式(1)-(3)的Dirichlet邊界條件和Newman邊界條件。通常將式(8)、式(9)稱為控制方程式(1)-(3)的柯西條件,并稱具有柯西條件的熱傳導(dǎo)問(wèn)題為熱傳導(dǎo)方程的柯西問(wèn)題,即邊界識(shí)別問(wèn)題。

在t=0時(shí)刻的初值條件為：

(10)

式(10)為該問(wèn)題的初值條件。

熱傳導(dǎo)邊界識(shí)別問(wèn)題可通過(guò)如下的Dirichlet邊界條件確定邊界移動(dòng)函數(shù)s(t):

(11)

式中，us(t)是一個(gè)給定的函數(shù)。它通常表示介質(zhì)的融點(diǎn),在我們的數(shù)值實(shí)例中取us(t)≡0。

2 基于FEM的廣義邊界控制法

考慮上述熱傳導(dǎo)方程的柯西問(wèn)題,首先用加權(quán)殘值法將偏微分方程轉(zhuǎn)化為弱形式下的有限元系統(tǒng)方程,然后用廣義邊界控制法將已知的Dirichlet邊界條件轉(zhuǎn)化為待定邊界的條件,并進(jìn)行求解。

2.1 廣義邊界控制法

由于FEM對(duì)于多層區(qū)域交界面上的Dirichlet和Newman邊界條件是自然滿足的,所以單層區(qū)域上的一維熱傳導(dǎo)問(wèn)題的公式可以應(yīng)用到多層區(qū)域。由加權(quán)殘值法很容易得到t+Δt時(shí)刻節(jié)點(diǎn)溫度滿足的整體有限元系統(tǒng)方程：

(12)

式中：K是剛度矩陣；Ut+Δt是由t+Δt時(shí)刻的節(jié)點(diǎn)溫度組成的向量。為了后續(xù)討論的方便,將Ut+Δt簡(jiǎn)記作U。因此,整體的有限元系統(tǒng)方程變?yōu)椋?/p>

(13)

(14)

式中：

B1F1+B2fn=B1F1+B2a.

(15)

從上式中解出a,并將fn=a代入(14)式即可求得U。根據(jù)求得的數(shù)值解U和移動(dòng)邊界滿足的Dirichlet邊界條件u(s(t+Δt),t+Δt)=0得到t+Δt時(shí)刻溫度值為零對(duì)應(yīng)的位置s(t+Δt)。

2.2 基于FEM的廣義邊界控制法的可行性分析

考慮兩個(gè)單元(圖2),推導(dǎo)基于FEM的廣義邊界控制法?？紤]該網(wǎng)格下的整體有限元系統(tǒng)方程:

KU=F.

式中：

圖2 具有兩個(gè)單元的一維離散區(qū)域Fig.2 One dimensional discrete domain with two element

按2.1節(jié)的討論,施加邊界之后可將整體的有限元系統(tǒng)方程寫(xiě)成如下的分塊矩陣:

式中：

則有:

式中：

因此,上式展開(kāi)為U=B1F1+B2a,即:

提取第1行得:

顯然,兩個(gè)單元時(shí),B2是恒為1的列向量,故可求得a。

假設(shè)n=k-1時(shí)，有:

那么,n=k時(shí),

3 數(shù)值實(shí)例

在數(shù)值試驗(yàn)中,假設(shè)模型的最大運(yùn)行時(shí)間T=1,為了檢驗(yàn)數(shù)值結(jié)果的計(jì)算誤差,定義T=1時(shí)刻區(qū)域上的溫度值u的均方根誤差R(u)為

(16)

式中,n是有限元中節(jié)點(diǎn)個(gè)數(shù),在下面計(jì)算中,n=11。不同時(shí)刻右邊界s的均方根誤差R(s)為

(17)

式中,m是在區(qū)間[0,T]內(nèi)均勻分布的時(shí)間點(diǎn)的總數(shù)。在計(jì)算中,m=21,l1=0,l2=0.2,l3=0.4,s(0)=1。

由于實(shí)際測(cè)量得到的邊值條件(柯西條件)和初值條件存在誤差,為了檢驗(yàn)該方法的有效性,使用如下的噪聲數(shù)據(jù):

(18)

(19)

(20)

式中:u0(ti)和q0(ti)分別表示ti時(shí)刻的問(wèn)題域的左邊界的精確溫度和熱流值；φj(x)表示t=0時(shí)刻問(wèn)題域的精確溫度；γ為一個(gè)均勻分布在[-1,1]之間的隨機(jī)數(shù)；δ表示相對(duì)噪聲水平。

下面將用3個(gè)數(shù)值實(shí)例來(lái)研究該問(wèn)題。

實(shí)例1:多層介質(zhì)的熱傳導(dǎo)問(wèn)題(式1-式3)的精確解為:

(21)

(22)

(23)

表1 不同噪聲水平下T=1時(shí)刻溫度值u和不同時(shí)刻右邊界s的均方根誤差

圖3給出了不同噪聲水平δ下移動(dòng)邊界的數(shù)值結(jié)果。表1給出了不同噪聲水平下T=1時(shí)刻溫度值u和不同時(shí)刻右邊界s的均方根誤差。由此可得:1) 基于噪聲數(shù)據(jù)得到結(jié)果的誤差數(shù)量級(jí)與噪聲水平δ的數(shù)量級(jí)基本相同或更小,說(shuō)明該方法較穩(wěn)定;2) 由表1的運(yùn)行時(shí)間可得,該方法具有有限元的高效性,它比基于RBF的Kansa法[9]與基本解法[13]的效率明顯高;3)t>0.1且噪聲水平δ為10%時(shí),基于FEM的廣義邊界控制法也有很好的效果；但是t≤0.1時(shí)誤差較大。這些誤差是由于初值和左邊界的擾動(dòng)(僅考慮溫度值的擾動(dòng))引起的,那么下面分析一下誤差較大的主要原因。

圖4給出了當(dāng)左邊界值為精確解時(shí),不同初值下的移動(dòng)邊界。表2給出了T=1時(shí)刻溫度值u在不同初值下的均方誤差。注意到當(dāng)初值為精確值時(shí),移動(dòng)邊界的數(shù)值結(jié)果和精確解吻合的相當(dāng)好(圖4(a))。圖4(b)-(d)分別是初值恒為0,5,10時(shí)移動(dòng)邊界的數(shù)值結(jié)果和精確解。它們均是在t≤0.3時(shí)誤差很大,t>0.3時(shí)吻合的比較好。這說(shuō)明基于FEM的廣義邊界控制法得到的T=1時(shí)刻溫度值u的均方誤差對(duì)初值沒(méi)有依賴。因此,該問(wèn)題在T=1時(shí)刻溫度值u的均方誤差主要是由左邊界的溫度值引起的。

圖3 不同噪聲水平δ下的移動(dòng)邊界的數(shù)值結(jié)果Fig.3 Numerical result s(t) using different noisy level δ

(a)exact initial condition;(b)the initial value is 0;(c)the initial value is 5;(d)the initial value is 10圖4 不同初值下的數(shù)值移動(dòng)邊界與精確移動(dòng)邊界Fig.4 Numerical results using different initial condition

表2 T=1時(shí)刻溫度值u在不同初值下的均方誤差

實(shí)例2:多層介質(zhì)的熱傳導(dǎo)問(wèn)題式(1)-(3)的精確解為:

(24)

(25)

(26)

圖5給出了不同噪聲水平δ下的移動(dòng)邊界的數(shù)值結(jié)果。表3給出了不同噪聲水平下T=1時(shí)刻溫度值u和不同時(shí)刻右邊界s的均方根誤差。由此可得：1) 基于噪聲數(shù)據(jù)得到結(jié)果的誤差數(shù)量級(jí)與噪聲水平δ的數(shù)量級(jí)基本相同或更小,再次說(shuō)明該方法較穩(wěn)定；2) 本例的移動(dòng)邊界s的均方根誤差比實(shí)例1稍微大些,這是因?yàn)橐苿?dòng)邊界s是隨時(shí)間的遞增函數(shù),需要使用外插求解s,而外插的誤差相對(duì)于內(nèi)插稍大些；3) 表3中運(yùn)行時(shí)間也再次說(shuō)明了該方法的高效性；4)t>0.2且噪聲水平δ=10%時(shí),基于FEM的廣義邊界控制法有很好的效果，但t≤0.2時(shí)誤差較大。同例1一樣,這些誤差仍主要是由初值的擾動(dòng)引起的。故本例也說(shuō)明使用基于FEM的廣義邊界控制法得到的T=1時(shí)刻溫度值u的均方誤差對(duì)初值沒(méi)有依賴。

表3 不同噪聲水平下T=1時(shí)刻溫度值u和不同時(shí)刻右邊界s的均方根誤差

圖5 不同噪聲水平δ下的移動(dòng)邊界的數(shù)值結(jié)果Fig.5 Numerical result s(t) using different noisy level δ

實(shí)例3:假設(shè)移動(dòng)邊界函數(shù)為:

且滿足u3(s(t),t)=0，Neumann邊界條件q0(t)=-1/3,在x=0處的Dirichlet邊界條件能夠通過(guò)相應(yīng)的正問(wèn)題獲得,正問(wèn)題式(1)-式(3)滿足的條件:

2) Dirichlet邊界條件u3=(s(t),t)=0.

3) 交界面的條件(式(4)-(7)).

4) 初值條件:

該正問(wèn)題可用FEM求解,并取ρ1c1=2/3，ρ2c2=3/4，ρ3c3=2，k1=6，k2=3，k1=2.

圖6給出了不同噪聲水平δ下移動(dòng)邊界的數(shù)值結(jié)果。表4給出了不同噪聲水平下T=1時(shí)刻溫度值u和不同時(shí)刻右邊界s的均方根誤差。這表明反問(wèn)題得到的結(jié)果能很好地吻合精確解。這些結(jié)果也再次證實(shí)了上面實(shí)例得到的結(jié)論。

表4 不同噪聲水平下T=1時(shí)刻溫度值u和不同時(shí)刻右邊界s的均方根誤差

圖6 不同噪聲水平δ下的移動(dòng)邊界的數(shù)值結(jié)果Fig.6 Numerical result s(t) using different noisy level δ

4 結(jié)論

本文研究了使用基于有限元的廣義邊界控制法重構(gòu)多層區(qū)域上一維熱傳導(dǎo)問(wèn)題的移動(dòng)邊界,并證明了該方法的可行性。大量的數(shù)值方法說(shuō)明了該方法的穩(wěn)定性和有效性。

[1]HONYC,WEIT.Afundamentalsolutionmethodforinverseheatconductionproblem[J].EngineerAnalysiswithBoundaryElements,2004,28:489-495.

[2]ENGLHW,HANKEM,NEUBAUERA.Regularizationofinverseproblems[M].Dordrecht,USA:KluwerAcademicPublishersGroup,1996.

[3]ALESSANDRINIG.Examplesofinstabilityininverseboundary-valueproblems[J].InverseProblems,1997,13(4):8874-897.

[4]ALESSANDRINIG,MORASSIA,ROSSETE.Detectingcavitiesbyelectrostaticboundarymeasurements[J].InverseProblems,2002,18(5):1333-1353.

[5]LIUGR,QUEKSS.Finiteelementmethod:apracticalcourse[M].BurlingtonMA:BH, 2003.

[6]LIUGR,NGUYEN-THOIT.Smoothedfiniteelementmethods[M].BocaRaton,USA:CRCPress,2010.

[7]LIUGR.Meshfreemethods:movingbeyondthefiniteelementmethod[M].2ndEd.BocaRaton,USA:CRCPress,2009.

[8]WEIT,LIYS.Aninverseboundaryproblemforone-dimensionalheatequationwithamultilayerdomain[J].EngineerAnalysiswithBoundaryElements,2009,33:225-232.

[9]NIURP,LIUGR,LIM.Reconstructionofdynamicallychangingboundaryofmultilayerheatconductioncompositewalls[J].EngineeringAnalysiswithBoundaryElements,2013,42:92-98.

[10] 王明清.偏微分方程中兩個(gè)不適定問(wèn)題數(shù)值解法的研究[D].太原:太原理工大學(xué),2014.

[11]BADIAAEI,MOUTAZAIMF.Aone-phaseinverseStefanproblem[J].InverseProblems,1999,15(6):1507-1522.

[12]FREDMANTP.Aboundaryidentificationmethodforaninverseheatconductionproblemwithanapplicationinironmaking[J].HeatMassTransfer,2004,41:95-103.

[13]LIYS,WEIT.Identificationofamovingboundaryforaheatconductionprobleminamultilayermedium[J].HeatMassTransfer,2010，46:779-789.

(編輯：李文娟)

A General Boundary Control Method Based on FEM for Boundary Identification in a Multi-Layer Heat Conduction

YUE Junhong,LI Ming,NIU Ruiping

(College of Mathematics,Taiyuan University of Technology,Taiyuan 030024,China)

The general boundary control method (GBCM) based on finite element method (FEM) has been applied to solve inverse solid mechanics problems. The identification of a moving boundary for multi-layer heat conduction with Cauchy boundary data is a typical inverse heat conduction problem. In this paper, the boundary identification problem of one-dimensional multi-layer heat conduction is solved using the GBCM based on FEM. Meanwhile, the feasibility of this method is proved. Numerical experiments are presented to demonstrate that the method is effective and stable.

finite element method;general boundary control method;multi-layer heat conduction;boundary identification

1007-9432(2016)04-0545-07

2015-03-27

國(guó)家自然科學(xué)基金資助項(xiàng)目:可商業(yè)化光滑有限元分析軟件的研發(fā)(11472184);國(guó)家青年科學(xué)基金資助項(xiàng)目:牙種植技術(shù)中的多參數(shù)識(shí)別問(wèn)題的計(jì)算方法(11401423)

岳俊宏(1989-),男,山西晉中人,博士生,主要從事計(jì)算數(shù)學(xué),(E-mail)woyuejunhong@163.com

李明,教授,博導(dǎo),主要從事計(jì)算數(shù)學(xué),(E-mail)liming01@tyut.edu.cn

O241

A

10.16355/j.cnki.issn1007-9432tyut.2016.04.022