热传导问题的MATLAB数值计算

《热传导问题的MATLAB数值计算》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、第30卷第9期华中科技大学学报（自然科学版）Vol.30No.92002年9月J.HuazhongUniV.ofSci.&Tech（.NatureScienceEdition）Sep.2002热传导问题的MATLAB数值计算李灿高彦栋黄素逸（湖南冶金职业技术学院）（华中科技大学）冶金系能源与动力工程学院摘要：分析了应用MATLAB中PDE工具箱解热传导问题的方法和步骤，编制了三个难以用解析方法求解的算例.采用有限元法求解导热偏微分方程，应用PDE工具箱得到数值解.对适合圆柱坐标描述的问题，通过公式变化将

2、其转换为能用PDE工具箱求解的形式.算例表明，用MATLAB对复杂形状和有内热源的非稳态导热问题进行数值计算和图形处理是方便高效的.关键词：热传导；非稳态导热；MATLAB；数值计算中图分类号：TK124文献标识码：A文章编号：1671-4512（2002）09-0091-03""许多工程问题需要确定物体内部的温度场或#-·（c）+a=f；""确定其内部温度到达某一限定值所需要的时间，d（!／!t）-·（c）+a=f；"""（1）因此研究导热问题特别是非稳态导热问题十分重d（!2／!t2）-·（c）+a

3、=f；""要.目前非稳态导热问题的描述方程为多维非线$-·（c）+a=!d，性的偏微分方程，这些方程只在几何形状与边界式中，为域"上的求解变量；!为特征值；d，c，条件都较简单的情况下才能求得理论解，而对于a，f为常数或变量；t为时间变量.前3个方程分几何形状和边界条件复杂的情况多用数值解法，别称为椭圆方程、抛物线方程和双曲线方程，第4需借助于计算机将时间和空间坐标划分成数量巨个方程称为特征值方程.大的网格才能得到较精确的数值解.本文应用导热问题的通用微分方程可写成［2］""MATLAB中PDE工具箱，

4、求解复杂边界条件的#c1（!／!t）=（$）+ZV，（2）热传导问题.式中，为求解变量，此处表示被求解物体内的温度；$为导热系数；Z为热源的发热率密度；为V#!求解方法密度；c为定压比热容.1可以看出，式（1）和式（2）中的抛物线方程求解方法是基于数值解法中的有限元法［1］，有着类似的形式.其中，求解变量为区域的温度，其基本原理是把计算区域划分成一系列的三角形d与#c，c与$，f-a与Z可以一一对应.1V单元，每个单元上取一个节点，选定一个形状函数MATLAB中的PDE工具箱定义了两类边界（抛物线形或双

5、曲线形），并通过单元中节点上的条件被求解变量值表示该函数.通过对控制方程作积h=T；｛"分来获得离散方程.有限元法的最大优点是对不，（3）!·（c）+Z=g规则区域适应性好，故用MATLAB方法求解的式中，!为垂直于边界的单位矢量；h，T，Z和g结果在边界上也较精确.对于适合圆柱坐标和球为常量或与有关的变量.方程（3）中的第1个方坐标描述的问题，通过对其热传导方程的变换，也程称为狄利克雷（Dirichlet）边界条件，第2个方程能在MATLAB中求解.称为纽曼（Neumann）边界条件.可以看出，导热问

6、应用MATLAB的PDE（PartialDifferential题中的第一类边界与狄利克雷边界条件对应，第!ECuation）工具箱可以解如下四类偏微分方程二类和第三类边界条件与纽曼边界条件对应.这些对应关系可以使用MATLAB中的PDE工具箱收稿日期：2002-05-09.作者简介：李灿（1968-），女，副教授；株州，湖南冶金职业技术学院冶金系（412000）.!PartialdifferentialeCuationtoolboxuser/sguide.TheMathWorks，Inc.，2000.

7、92华中科技大学学报（自然科学版）第30卷来求解.如果导热体物性系数!为温度的函数，只对一个导热问题的计算可以按图l的步骤要写出!（）的函数关系式，就可以得到解进行.!!有内热源的圆柱体非稳态导热问题有一半径为02m，长为3m的圆柱形核电站用燃烧棒置于f=l00C的水中，由于链式反应，棒内有恒定的产热率密度=20000W／m3，V计算l0h后燃烧棒内的温度分布已知，燃烧棒的密度／m3，=500W·s／（kg·C），#=7800kg!=40W／（m·C），0=0C，燃烧棒右端恒温图lMATLAB计算流程图

8、=l00C，左端有一恒热流=5000W／m2，燃rl烧棒外表面与外界水的对流换热系数=50!算例W／（m2·C）此问题宜采用圆柱坐标，由于燃烧棒内温度!."三维非稳态无内热源的导热问题沿半径对称分布，因此可以转换为（，）坐标的边长为0.5m，0.7m和l.0m的长方形钢二维问题将圆柱坐标内的热传导方程改写为锭，置于炉温f=l200C的加热炉内，计算5h!!!!!#（!）（!）V，后钢锭的温度已知钢锭的!=405W／（m·C），!!!!!"=