广义边界控制法在多层热传导边界识别问题中的应用

岳俊宏,李 明,牛瑞萍

(太原理工大学 数学学院,太原 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