二维热传导方程有限容积法的MATLAB实现

2012，48（24）1引言热传导广泛存在于各种材料成型过程中，多晶硅铸坯的凝固过程就可以看成是具有均匀厚度的无限大板的传热过程，另外轴对称金属铸件的定向凝固过程可以被简化作为二维截面数值模型从而实现传热模拟。要掌握这些成型过程的最优工艺参数，关键是正确求解关于传热问题的偏微分方程[1]，典型的数值方法有有限元[2]、有限差分[3]和有限容积法[4]。MATLAB是一款具有精确数值计算功能和丰富图形处理函数的软件[5]，传热偏微分方程的数值解可直观地以二维、三维图形方式显示在屏幕，之前有一些文章介绍用MATLAB实现有限差分方法求解传热问题[6]，本文主要讨论一、二维热传导方程有限容积法的MATLAB实现。2传热数值建模与推导2.1有限容积法简介有限容积法是一种通过代数方程组代替原有难以直接求解的传热偏微分方程的数值计算方法[7]，与有限差分法、有限元方法非常相似，都是求解离散单元网格节点值从而获得传热问题的解。如图1所示，对于二维问题，“有限体积单元”是指以节点P为中二维热传导方程有限容积法的MATLAB实现薛琼1，肖小峰2XUEQiong1,XIAOXiaofeng21.武汉理工大学理学院，武汉4300702.武汉纺织大学机械工程与自动化学院，武汉4300731.SchoolofScience，WuhanUniversityofTechnology，Wuhan430070，China2.CollegeofMechanicalEngineeringandAutomation，WuhanTextileUniversity，Wuhan430073，ChinaXUEQiong,XIAOXiaofeng.TwodimensionalheatconductionequationfinitevolumemethodtoachieveonMATLAB.ComputerEngineeringandApplications,2012,48（24）：197-200.Abstract：Twodimensionaldiffusionlessheatconductionphenomenonhasbeendescribedonpartialdifferentialequation.Basedonfinitevolumemethod,discretizedalgebraicequationofpartialdifferentialequationhasbeende-duced.Differentcoefficientsandsourcetermshavebeendiscussedunderdifferentboundaryconditions,whichin-cludeprescribedheatflux,prescribedtemperature,convectionandinsulation.TransientheatconductionanalysisofinfiniteplatewithuniformthicknessandtwodimensionalrectangleregionisrealizedbyprogrammingusingMAT-LAB.Itisusefultomaketheheatconductionequationmoreunderstandablebyitssolutionwithgraphicalexpres-sion.Feasibilityandstabilityofnumericalmethodhavebeendemonstratedbyrunningresult.Keywords：heatconductionequation;finitevolumemethod;MATLAB摘要：通过偏微分方程描述了二维无扩散热传导现象。基于有限容积法推导了该方程的离散代数方程组，针对恒定热流强度、恒定温度、对流换热和绝热这四种不同边界条件，分别讨论了热传导代数方程的离散系数和源项。通过MATLAB编程，分析了一维具有均匀厚度无限大板和二维矩形区域的瞬态传热现象。采用图形显示方式使得偏微分方程求解更为直观和容易理解，计算结果证明了有限容积求解方法是可行、稳定的。关键词：热传导方程；有限容积法；MATLAB文章编号：1002-8331（2012）24-0197-04文献标识码：A中图分类号：O551基金项目：国家自然科学基金（No.10901067）；中央高校基本科研业务费专项资金（No.2011-1a-023）。作者简介：薛琼（1980—），女，博士，讲师，主要研究领域为微分几何及其应用；肖小峰（1979—），男，讲师。E-mail：[e-mail]收稿日期：2011-12-08修回日期：2012-01-30CNKI出版日期：2012-05-21DOI：10.3778/j.issn.1002-8331.2012.24.044计算机工程与应用197ComputerEngineeringandApplications计算机工程与应用2012，48（24）心，由东（e）、西（w）、南（s）、北（n）方向界面所形成的封闭面积（图1中的虚线框），相邻的体积单元通过各自界面传递热流量（图1中flux），某个有限体积单元的流入热量和与之相邻有限体积单元的流出热量相等，所以它是一种能量守恒的计算方法，因此有限容积法被广泛应用于热量、动量传递问题的数值求解[8]。2.2数学模型与离散方程只考虑扩散二维瞬态的传热问题的偏微分方程如下式：ρc¶T¶τ=¶¶xæèöøλ¶T¶x+¶¶yæèçöø÷λ¶T¶y+S（1）x、y分别表示水平、垂直方向的空间坐标；ρ、c、τ、T、λ、S分别表示密度、热容量、时间、温度、热传导系数、能源项。根据传热偏微分方程式（1），有限容积单元积分如下式：ττ+Dτwesnρc¶T¶τdxdydτ=ττ+Dτwesnæèçöø÷¶¶xæèöøλ¶T¶x+¶¶yæèçöø÷λ¶T¶y+Sdxdydτ（2）对式（2）积分推导可得二维传热偏微分方程的等效离散代数方程如下式：ìíîïïïïïïïïïïïïaPTP=aETE+aWTW+aNTN+aSTS+SuDV+a0PT0PaP=aE+aW+aN+aS+a0P-SPDVa0P=ρcDVDτaE=λeDyδxeaW=λwDyδxwaN=λnDxδynaS=λsDxδysDV=DxDy（3）2.3边界条件图2所示有限容积单元的w界面不和其他有限容积单元相连，处于物理介质的边界上，所以规定式（3）中对应界面温度系数aW为0，但这并不意味w边界对整个传热系统没有影响，这需要根据不同边界类型定义不同源项S，并将其代入到方程组的迭代求解中，从而在数学物理模型上体现不同方式向物理介质内的有限容积单元传递热量，下面分三种情况进行讨论：（1）给定热流密度q边界条件热流密度q=-λT2j-T1jδxw，边界温度系数aW=aW2j=0，又有源项Su2j=qDyDV；则可计算出边界单元温度：T1j=qδxwλ+T2j。（2）对流换热的边界条件热流密度q=-λT2j-T1jδxw=h(Tf-T1j)，其中h、Tf分别是对流换热系数、流体温度。温度系数aW=aW2j=0，又有源项：Su2j=Tf(1/h+δxw/λ)DyDVSP2j=-1(1/h+δxw/λ)DyDV则可计算出边界单元温度：T1j=hTf+(λ/δxw)×T2jh+λ/δxw。（3）绝热边界条件温度系数aW=aW2j=0，则计算出边界单元温度：T1j=T2j。对于东边界上的有限容积单元，边界条件的讨论及推导过程原理完全一样，只需将索引1用N代替、2用N-1代替，以此类推，置0的温度系数aW用aE代替。3MATLAB数值算例3.1均匀厚度无限大板传热分析（1）数学模型均匀厚度无限大板传热问题可以被描述为只考虑扩散、一维瞬态热传导偏微分方程，该方程很容易NfluxδynδysδxeδxwDxSfluxfluxDyEensfluxwPW图1有限容积法的网格单元与节点图2西边界上的有限容积单元WestBoundaryDxδxw(1j)(2j)Dy1982012，48（24）通过式（1）和式（3）推导出，具体如下式：ìíîïïïïïïïïïïïïïïïïρc¶T¶τ=¶¶xæèöøλ¶T¶x+S|||¶T¶xx=0=hλ()T1-Tf|||¶T¶xx=0.06=0aPTP=aETE+aWTW+SuDx+a0PT0PaP=aE+aW+a0P-SPDxa0P=ρcDVDτaE=λeδxeaW=λwδxw（4）数值计算中采用的物理量及参数如表1所示，因为板厚方向温度分布是关于中间平面对称的，所以厚度方向空间坐标X的定义域为板厚的一半，具体为表1。（2）关键MATLAB程序代码%initializationfork=1：ni；t（k）=Tinit；told（k）=Tinit；ap（k）=density*cp*dx/dt；ae（k）=conductivity/dx；aw（k）=conductivity/dx；if（k==2）；aw（k）=0；end；if（k==ni－1）；ae（k）=0；end；d（k）=ap（k）－aw（k）－ae（k）；end；%algebraicequationdd=d（k）*told（k）；dd=d（k）*told（k）+（Tfluid－told（k））/（1/alfa+dx/（2*con-ductivity））；t（k）=（ae（k）*told（k+1）+aw（k）*told（k－1）+dd）/ap（k）；%Boundaryt（ni）=t（ni－1）；t（1）=（alfa*Tfluid+2*conductivity/dx*t（2））/（alfa+2*con-ductivity/dx）；data（j，1：32）=told；（3）计算结果与讨论取时间等于2s、10s、50s、100s、160s，运行程序，板厚方向的温度分布如图3所示。另取时间等于100s，运行下段程序：%Analyticalsolutionastime=100sJT（ni－kk+1）=（1.2102*exp（－（1.1925）2*0.277778）*cos（1.1925*（（kk－1）*dx－dx/2）/L））*（Tinit－Tfluid）+Tfluid；plot（x，data（1001，1：32），‘*’）；plot（x，JT，‘－’）；holdon；有限容积法的数值解和精确解析解如图4所示，发现二者结果非常接近。继续运行下段程序：RuncodesbellowonMATLAB：%Numericalsolutiononthreedimensionalcoordinatew（1：32）=（p－1）*2；q=（p－1）*10*2+1；plot3（w，x，data（q，：））；板厚方向的瞬态温度分布通过三维坐标直观显示在屏幕上，这使得偏微分方程的解非常容易理解。3.2二维矩形区域的传热分析（1）数学模型参数类型无限大板厚度L厚度方向坐标x密度ρ比热c大板初始温度Ti流体温度Tf热传导系数λ换热系数h数值0.120~0.0680005008020402000单位mmkg/m3J/(kg·K)℃℃W/(K×m)W/(K×m2)表1数值计算中采用的物理量85807570656055504540353000.010.020.030.040.050.06X坐标/m温度/(℃)t=160t=100t=50t=2t=10图3典型时间的数值计算结果706560555045403500.010.020.030.040.050.06X坐标/m温度/(℃)numericalsolutionanalyti