第1页共16页一维非稳态热传导热源反问题研究摘要本文是关于热传导的正反问题的研究,即利用偏微分方程中典型热传导方程求解含有内热源的金属细杆t时刻温度分布与热源位置。本文从解偏微分方程出发,由已知条件最终得出温度分布函数及热源位置函数并建立了两个数学模型。模型一:利用偏微分方程及初始温度分布函数建立了一段时间后的温度分布与热源强度、位置之间的数学模型,最终解出一段时间后长杆上的温度分布。模型二:通过一类抛物型偏微分方程模型,解决已知初始温度分布函数、一段时候后的温度分布函数及热源强度的确定热源位置和中间任意时刻的温度分布函数。根据模型一建立偏微分方程组,用分离变量法求解(,)uxt,即t时刻的温度分布函数,并通过Matlab中的PDE(偏微分方程)工具箱求解偏微分方程组,且使解可视化。根据模型二依然建立偏微分方程组,通过测得(,)uxT,结合抛物型方程,运用离散正则法,确定热源位置,并通过论证说明问题的唯一性和确定性,给出反问题的数值解法。最后再简单介绍差分法解决热传导在非稳态导热问题中的应用。最后是结论部分,主要总结本文的结果并提出一些尚待进一步研究的问题,以及研究该反问题的应用前景。相同t不同x的温度变化曲线相同x不同t的温度变化曲线第2页共16页一维非稳态热传导热源反问题研究一、问题的提出在金属细秆的传热过程中,温度差是导致其发生必要条件,有无热源决定传导效率的高低。从一维非稳态传导问题的数学模型和初始条件出发,经过对有内热源问题的进一步分析,在初始温度分布已知的情况下,对分布函数的处理显得很关键。对热源反问题的处理中,我们的问题是如何寻找某种合理的附件条件,通过已知方程来解决方程右端的热源的具体位置并使其具有唯一性。本文利用微分方程并建立了满足温度分布的数学物理模型,从理论上导出了温度分布函数和热源位置的求解,并借助计算机软件画出了温度分布图。二、问题的分析对于热传导问题,为了使函数解决起来更容易,对于细秆的初始温度分布()gx我们可以设它在区间[0,L]连续,那么()gx可以展成正弦或余弦级数,对于有内热源的处理,由于细秆边界条件是齐次的,我们采用叠加原理把一根金属细秆的导热问题分解为有热源的具有其次边界条件的稳态导热问题和一个非稳态其次问题,则原问题的解为(,)1(,)2()uxtuxtux。对于源反问题的解决有如下3个问题:1、反问题的唯一性:附加条件给得是否合理,也就是说,这个附加条件是否可以唯一确定热源的具体位置。2、反问题的稳定性:反演所得到的热源的具体位置,该热源是否是连续地依赖于测量数据()ht?3、反问题的数值解法:如何用可行的数值方法反演该热源的具体位置。用离散正则法将温度分布离散化,由已知初始温度分布再利用计算机软件得出热源位置三、模型假设1、金属细杆边界与外界无热量交换,即与外界绝缘第3页共16页2、热源强度在整个时间段里始终保持常量。3、在求解源反问题时,热源分布相对细杆长度来说,可假设为点热源。四、符号说明(,)uxt:温度分布函数,即x处在t时刻的温度(,)fxt:热源强度a2:热扩散系数,单位为m2/sL:细杆长度()gx:初始温度分布函数,即t=0时的杆在x处温度其他运算过程中使用符号在步骤中说明,再不赘述。五、模型的建立与求解5.1建立热传导微分方程,并求出温度分布函数(,)uxt由题意可的模型如下:222()uuafxtx0xLt0(1)(0,)(,)0utuLt0x(2)(,0)()uxgx0xL(3)根据可以把有内热源的非稳态导热问题分解为有热源的具有齐次边界条件的稳态导热问题和一个非稳态齐次问题,即原方程的解为:(,)1(,)2()uxtuxtux(4)式中,1(,)uxt是如下非稳态的解:(方程组A)22211uuatx(A-1)1(0,)1(,)0utult(A-2)第4页共16页1(,0)()2()()uxgxuxpx(A-3)2()ux是如下稳态问题的解:(方程组B)2222()0uafxx0xL2(0)2()0uul对方程组A:利用分离变量法,假定偏微分方程的解是两个独立变数的乘积,即设1(,)()()uxtXxYt代入微分方程(A-1)中,可得'2''XYaYX或'''2XYXaY为常数(A-4)(若为+,则推导后所得出的解,其结果将对取任何值都不能满足边界条件),再由边值条件,有u1(0,t)=X(0)Y(t)=0u1(L,t)=X(L)Y(t)=0必有(0)()0xxL。由(A-4),有''0XX;X(0)=X(L)=0(A-5);'20YaY(A-6);至此,通过分离变量,我们把微分方程的边值问题转化成为常微分方程的边值问题。要是(A-5)有非零解,只要取=n=222nL(n=1,2,3···)对应的非零解为sinnxL,类似方程(A-6)的解2()()exp(())nnaYtYttL于是得到方程组A的一组线性无关的解:2(,)sinexp(())nnxnauxttLL(n=1,2,3···)第5页共16页通过适当取值,设p(x)是函数sinnxL的有限线性组合即1()sinNnnxpxcL最终得到21(,)sinexp(())NnnxnauxtctLL其中,02()sinLnnxcpxdxLL2121(,)()sinsinexp(())NnxnxnauxtpxdxtLLLL显然,(,)uxt满足边值条件(,0)()uxpx对方程组B:即考虑有内热源的情况下,利用傅里叶变换求解2222()0uafxx0xL2(0)2()0uul的解,得出()2()()FLuxFxxL式中,()Fx是21()fxa的二次积分。最后,利用叠加定理将方程组A与B合并,即为所建模型温度分布函数的解:21(,)1(,)2()2()()sinsinexp(())()NuxtuxtuxnxnxnaFLpxdxtFxxLLLLL5.2已知T时刻的温度分布和初始温度分布,求解热源位置。对于求源问题,由已知初始温度分布和一段时间后细秆的温度分布,那即为(,)()uxTfx来求热源位置。a、离散正则化一维问题的一般提法为:第6页共16页(,)(,)(,)0,(0,)(,0)()0(0,)()(0,)(0,)(,)0(0,)txxxuxtuxtfxxtTuxgxxutttTutulttT(A)式中,边界条件(0,)()(0,)xutttT当(,)fxt,()gx,()t已知的情况下,原方程变为一维热传导方程的正问题,问题5.1就是属于一类一维热传导问题,对此借助参考文献[1]得出另一解为:000(,)2(,)()()(()())ttuxtfdkxtgkxtkxtd(B)其中,21(,)exp()44xkxttt对于已知(0,)(,)0utult,将其换为附加条件(0,)()0utttT(C)此时根据已知条件,原问题就构成一维热传导方程的寻源反问题,把C式代入B式中可得到:000()2(0,)()()(()())ttfdtktdgktktd(D)以下为使用离散正则法解决源反问题:式(D)可写成第一类积分算子方程的形式:,,AfFfZFU(E)其中,A为第一类积分算子,Z为解空间,U为数据空间。由于原问题是不适定的,故将积分算子离散化后得到的线性方程组常常是病态的,而且随着维数的增加其病态状况更加严重。这时,直接求解离散方程显然是不行的,必须将其正则化后求解才能得到稳定的近似解。离散正则化的过程可归结如第7页共16页下:1)将积分方程离散化得,,hhhhhhhAfFfWFL其中,hA为积分算子的离散矩阵;hf为方程的离散解;hF为近似数据;h为离散步长;2)形成相应的法方程ThThhhhAAfAF3)法方程的正则化'()()ThhhhaaAACfCF其中,C为单位矩阵;4)利用偏差原理决定正则参数()a,使得((ha,其中为实际数据与测量数据之间的误差,((hhhhaAfF;5)求得离散正则解1()()hTThahhhfAACAFb、对离散正则模型进行数值实验(注:本文实验模拟特殊函数借助文献2)设T时刻细秆的温度分布:T=1,2()0.00001(1),01Tgxxxx3410.0002(14(exp(1))),01324ttttttt4.5256(),01315ttt则根据拟最优准则4得345()325TtttFt由式E得22(1)Tftt。在实际应用中,得到的只是一些离散的数据,并且存在着一定的误差故给()TFt一个扰动作为式E的右端项()fFt,应用离散正则化方法来求解()fFt作为右端项的第一类积分算子方程,求解结果为()fft。将Tf与()fft进第8页共16页行比较,从而验证离散正则化方法用来求解热传导方程寻反问题的正确性。C、寻求源的函数表达式1)离散第一类积分算子应用矩形求积公式的和来代替积分的方法对积分算子进行离散,可得到求式E的近似解的Euler方法10(),iifijhfFti=1,2,···,N其中,if是()fift的近似值,等价的,定义12((),(),~())NTffffNFFtFtFt121(,,~)NTNNffffR得到下面的矩阵方程NNNfAfF其中,NA是一个NN的下三角矩阵,其元素依赖于步长h。111NA2)由于NA为病态矩阵,并且随着维数的增加,病态特性逐渐增加,应用正则化方法求解该线性方程组。取T=1,当300.01,,10,100N时,模拟图形见图1;当300.01,,10,200N时,模拟图形见图2;图像见附录。D、求解热源位置随机设初始时刻温度分布函数3()0.0001()Tgxxx函数图像:第9页共16页最大值取的X位置约为0.58处。五、结果分析与检验(MATLAB模拟与算法)对问题5.1的偏微分方程组的matlab解法(PDEtools使用):(1)用GUI求解。打开PDEToolboxGUI窗口,在应用模型中选HeatTransfer。(2)设置坐标的显示范围。单击OptionsGridSpacing,在对话框中选择Auto;打开AxesLimits对话框,设置X,Y范围为[0,10](3)绘制区域:设置矩形长8,宽0.5,热源长1,宽0.5,位置任意。并满足默认边界条件的Dirichlet边界条件。再在PDEmode模式下,分别双击两个区域设置抛物型(Parabolic)Mode参数:矩形区域:密度rho为7.8,比热C为0.465,导热系数k为54(以钢介质为例)热源:rho为1,比热C为0.1,导热系数k为2.热源Q为4,取h=Text=0。打开SolveParameters,输入时间向量time为logspace(-2,-1,10)(4)模拟实现。a、在plotselection对话框选择color,height(3-dplot)和animation,单击plot,可显示解的动画过程。b、在plotselection对话框选择color,contour和arrows,单击plot。附截图见附录对问题5.2有解的结果图,见附录。流程图:第10页共16页六、模型的优缺点及改进方向本文采用两个模型,对问题的讨论由正问题到反问题,都给出了详细的推导并进行了严格的数学分析与推导,为使模型建立在坚定数学基础上,模型的建立采用导热微分方程及初始条件得出温度分布和热原位置,对于用离散正则法求源