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

岳俊宏,李 明,牛瑞萍

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

?

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

岳俊宏,李 明,牛瑞萍

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

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

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

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

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

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

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

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

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

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

(1)

(2)

(3)

式中：ui(x,t),i=1,2,3分別是區(qū)域Di,i=1,2,3內(nèi)的溫度分布；T是模型運行的最大時間；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)問題滿足如下的交界面條件:

(4)

(5)

(6)

(7)

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

(8)

(9)

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

在t=0時刻的初值條件為：

(10)

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

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

(11)

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

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

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

2.1 廣義邊界控制法

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

(12)

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

(13)

(14)

式中：

B1F1+B2fn=B1F1+B2a.

(15)

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

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

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

KU=F.

式中：

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

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

式中：

則有:

式中：

因此,上式展開為U=B1F1+B2a,即:

提取第1行得:

顯然,兩個單元時,B2是恒為1的列向量,故可求得a。

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

那么,n=k時,

3 數(shù)值實例

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

(16)

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

(17)

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

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

(18)

(19)

(20)

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

下面將用3個數(shù)值實例來研究該問題。

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

(21)

(22)

(23)

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

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

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

圖3 不同噪聲水平δ下的移動邊界的數(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ù)值移動邊界與精確移動邊界Fig.4 Numerical results using different initial condition

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

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

(24)

(25)

(26)

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

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

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

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

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

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

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

4) 初值條件:

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

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

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

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

4 結(jié)論

本文研究了使用基于有限元的廣義邊界控制法重構(gòu)多層區(qū)域上一維熱傳導(dǎo)問題的移動邊界,并證明了該方法的可行性。大量的數(shù)值方法說明了該方法的穩(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] 王明清.偏微分方程中兩個不適定問題數(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

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

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

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

O241

A

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