廣義邊界控制法在多層熱傳導邊界識別問題中的應用

岳俊宏,李 明,牛瑞萍

(太原理工大學 數學學院,太原 030024)

?

廣義邊界控制法在多層熱傳導邊界識別問題中的應用

岳俊宏,李 明,牛瑞萍

(太原理工大學 數學學院,太原 030024)

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

有限元;廣義邊界控制法;多層熱傳導;邊界識別

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

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

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

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

1 多層區域上的熱傳導邊界識別問題

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

在上述3層區域上滿足如下熱傳導控制方程:

(1)

(2)

(3)

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

該熱傳導問題滿足如下的交界面條件:

(4)

(5)

(6)

(7)

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

(8)

(9)

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

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

(10)

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

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

(11)

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

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

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

2.1 廣義邊界控制法

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

(12)

式中：K是剛度矩陣；Ut+Δt是由t+Δt時刻的節點溫度組成的向量。為了后續討論的方便,將Ut+Δt簡記作U。因此,整體的有限元系統方程變為：

(13)

(14)

式中：

B1F1+B2fn=B1F1+B2a.

(15)

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

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

考慮兩個單元(圖2),推導基于FEM的廣義邊界控制法。考慮該網格下的整體有限元系統方程:

KU=F.

式中：

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

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

式中：

則有:

式中：

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

提取第1行得:

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

假設n=k-1時，有:

那么,n=k時,

3 數值實例

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

(16)

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

(17)

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

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

(18)

(19)

(20)

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

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

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

(21)

(22)

(23)

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

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

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

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

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

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

(24)

(25)

(26)

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

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

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

實例3:假設移動邊界函數為:

且滿足u3(s(t),t)=0，Neumann邊界條件q0(t)=-1/3,在x=0處的Dirichlet邊界條件能夠通過相應的正問題獲得,正問題式(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給出了不同噪聲水平δ下移動邊界的數值結果。表4給出了不同噪聲水平下T=1時刻溫度值u和不同時刻右邊界s的均方根誤差。這表明反問題得到的結果能很好地吻合精確解。這些結果也再次證實了上面實例得到的結論。

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

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

4 結論

本文研究了使用基于有限元的廣義邊界控制法重構多層區域上一維熱傳導問題的移動邊界,并證明了該方法的可行性。大量的數值方法說明了該方法的穩定性和有效性。

[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] 王明清.偏微分方程中兩個不適定問題數值解法的研究[D].太原:太原理工大學,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

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

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

李明,教授,博導,主要從事計算數學,(E-mail)liming01@tyut.edu.cn

O241

A

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