传热学MATLAB温度分布大作业完整版
传热大作业 第4版4-23
东南大学能源与环境学院课程作业报告课程名称：传热学作业名称：传热学大作业——利用matlab程序解决热传导问题院（系）：能源与环境学院专业：热能与动力工程姓名：姜学号：完成时间：2012 年11 月8日评定成绩：审阅教师：目录一．题目及要求 (3)二．各节点离散化的代数方程..............................3&13 三．源程序......................................................5&16 四．不同初值时的温度分布情况...........................7&18 五．冷量损失的计算.......................................12&24 六．计算小结 (27)传热大作业——利用matlab 程序解决复杂热传导问题姓名:姜 学号: 班级:成绩：____________________一、题目及要求计算要求：一个长方形截面的冷空气通道的尺寸如附图所示。
假设在垂直于纸面的方向上冷空气及通道墙壁的温度变化很小，可以忽略。
试用数值方法计算下列两种情况下通道壁面中的温度分布及每米长度上通过壁面的冷量损失：（1） 内、外壁面分别维持在10℃及30℃；（2） 内、外壁面与流体发生对流传热，且有110f t C =︒、2120/()h W m K =⋅，230f t C =︒、224/()h W m K =⋅。
(取管道导热系数为0.025/()W m K λ=⋅)二、各节点的离散化的代数方程1、基本思想：将导热问题的温度场，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。
2、基本步骤：(1)建立控制方程以及定解条件：对于(1)问有：2.2m3m 2m1.2m h 1、t f1h 1、t f2导热微分方程22220t t x y ∂∂+=∂∂定解条件为第一类边界条件对(2)问有： 导热微分方程22220t t x y ∂∂+=∂∂定解条件为第三类边界条件(2)区域离散化：如下图所示，用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。
(完整word版)matlab绘制温度场
通过在室内的某些位置布置适当的节点，采集回来室内的温湿度以及空气质量等实际参数。
首先对室内空间建模，用一个无限细化的三维矩阵来模拟出室内的温度分布情况，针对采集回来的数据，采用插值法和适当次数的拟合函数的拟合，得出三维矩阵的实际值的分布，最后结合matlab软件绘制出计算出的温度场的三维图像。
一．数据的采集与处理因为影响人的舒适感的温度层只是室内的某一高度范围内的温度，而温度传感器虽然是布置在一个平面内，但是采用插值法和拟合函数法是可以大致再现出影响人的舒适感的温度层的温度变化的。
同时，在构建出的三维模型中，用第三维表示传感器层面的温度。
在传感器层面，传感器分布矩阵如下：X=【7.5 36.5 65.5】（模型内单位为cm）Y=【5.5 32.5 59.5】Z=【z1 z2 z3；z4 z5 z6；z7 z8 z9；】（传感器采集到的实时参数）采用meshgrid（xi，yi，zi，…）产生网格矩阵；首先按照人的最小温度分辨值，将室内的分布矩阵按照同样的比例细化，均分，使取值点在坐标一定程度上也是接近于连续变化的，从而才能最大程度上使处理数据得来的分布值按最小分辨值连续变化！根据人体散热量计算公式：C=hc（tb-Ta）其中hc为对流交换系数；结合Gagge教授提出的TSENS热感觉指标可以计算出不同环境下人的对环境温度变化时人体温度感知分辨率，作为插值法的一个参考量，能使绘制出的温度场更加的符合人体的温度变化模式。
例如按照10cm的均差产生网格矩阵（实际上人对温度的分辨率是远远10cm大于这个值的，但是那样产生的网格矩阵也是异常庞大的，例如以0.5cm为例，那么就可以获得116*108=12528个元素，为方便说明现已10cm为例）：[xi yi]=meshgrid(7.5:10:65.5,5.5:10:59.5)xi =7.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.5000yi =5.5000 5.5000 5.5000 5.5000 5.5000 5.500015.5000 15.5000 15.5000 15.5000 15.5000 15.500025.5000 25.5000 25.5000 25.5000 25.5000 25.500035.5000 35.5000 35.5000 35.5000 35.5000 35.500045.5000 45.5000 45.5000 45.5000 45.5000 45.500055.5000 55.5000 55.5000 55.5000 55.5000 55.5000产生网格矩阵之后，就可以在测得的实时数据的基础上，通过相关的温度场的专业的估算函数，以及相关的数值处理函数来估计整个分布面（有最小的分辨率）上的温度了。
报告
大作业报告【9-16】题目略 分析过程：● 第一步，建立节点温度差分方程。
有题目可以看出，本题有7类不同的节点，需要7个不同的方程。
由于y x ∆=∆在原始方程中将它们消去后有： 1、底边的中间5个点满足的方程：x q t t t t w j i j i j i j i ∆+++=-+-1,,1,1,222λλλλ2、右边对流换热的中间3个点满足的方程：∞∆+--∆+++=+t B t t t t B i j i j i j i j i i 22)42(1,1,,1,3、左边绝热部分的中间3个点满足的方程：1,1,,1,222+-+++=j i j i j i j i t t t t λλλλ4、上边恒温点满足的方程：w j i t t =.5、中间5×3的网格里的点满足的方程：1,1,,1,1,4+-+-+++=j i j i j i j i j i t t t t t6、左下角的既属于绝热层又接受加热的点：0)(2)(22,1,,,1=-+-+∆⋅++j i j i j i j i w t t t t x q λλ 7、右下角的既被加热又与流体有对流换热的点：0)(2)(2)(22,,1,,,1=-⋅∆⋅+-+-+∆⋅+-j i f j i j i j i j i w t t y h t t t t x q λλ● 第二步，采用高斯-塞德尔迭代法进行迭代：建立矩阵，左下点为原点(1,1)，i 值向右增加，j 值向上增加。
按照上面7种类型的分组依次进行计算，每次计算直接用计算结果代替原值，完成一轮迭代。
程序设定200次迭代。
根据最终结果判断，100次迭代就已经结果收敛。
● 结果分析： Matlab 运算结果：表1 节点的温度分布（i为横行，j为竖列）5 400 400 400 400 400 400 4004 397.1890 396.8375 395.7156 393.5919 389.9486 383.6080 371.52583 394.2443 394.4454 392.4332 388.7034 382.5945 372.9575 357.94582 394.2443 393.4299 390.8685 386.1939 378.7685 367.6818 351.89201 395.0363 394.1616 391.4168 386.4354 378.6037 367.1090 351.1356 节点/˚C 123456 7做出三维分布图像：结果分析：由数值计算结果，可以看出：1、最顶部的点温度恒定，满足题设条件。
2013年matlab应用作业
程序为：
>> x=0:2:24;
>> y=[12 9 9 10 18 24 28 27 25 20 18 15 13];
>> x1=13;y1=interp1(x,y,x1,'spline')
2）本金 以每年 次，每次 的增值率（ 与 的乘积为每年增值额的百分比）增加，当增加到 时所花费的时间（单位：年）为
用MATLAB表达式写出该公式并用下列数据计算： 。
程序为：
>> r=2;p=0.5;n=12;
>> T=log(r)/n/log(1+0.01*p)
甲（百箱）
乙（百箱）
现有
原料（kg）
6
5
60
工人（名）
10
20
150
获利（万元）
10
9
程序为：
>> c=[-10 -9];
>> a=[10 20;6 5];
>> b=[150;60];
>> aeq=[];
>> beq=[];
>> vlb=[0 0];
>> vub=[8 12];
>> [x,f]=linprog(c,a,b,aeq,beq,vlb,vub)
>> p=polyfit(t,y,2)
p =
-0.0445 1.0711 4.3252
>> x=30;
>> f(x)=-0.0445*x^2+1.0711*x+4.3252
利用matlab求解圆柱内稳定温度分布
数学物理方法论文（圆柱体齐次边界条件以及matlab的可视化）学院: 信息院年级: 2011级班级: 通信一班姓名: ***一、应用背景:现在由于公寓楼采用集体供暖的措施,需要在地下铺设圆柱形长管道,而对于管道的半径选取需要进行一定的计算,才能使得热水在传送过程中不会因温度过低而达不到用户的要求。
对于以上的问题，我们可以简化到应用数学中来解决。
二、简化例题：若一供暖公司采用的运输管道为标准圆柱体，其半径为ρ。
，长为L ，设管入口有均匀分布的强度为q 。
的热流流入，出口有相同的热流流出，管道侧面保持温度为0℃。
求解管道内的稳恒温度。
三、例题解答：解： 因为上下底非齐次边界的非齐次项是常数，故可较易化成齐次边界。
这样本征值问题就变成傅里叶级数本征问题，而不是贝塞尔函数本征问题，同时系数的求解也较为简单。
令 z 01kq =νv kq u +=0则v 的定解问题为：⎪⎪⎪⎩⎪⎪⎪⎨⎧=-===∆===z k q v L z z z z 00200,00ρρννν 分离变数得到的本征值问题为：⎪⎩⎪⎨⎧===+==0',0'0"02L z z Z Z Z h Z 解得 ()......2,1,0,==n Ln h n πz Ln A Z n n πcos =，问题在柱内的有限解为：()⎪⎪⎭⎫⎝⎛∞=⋅=∑ρππρνL n n n I z L n A z 00cos ， 由初始条件，得z kq IA n z L n L n 00ncos -=∑∞=⋅⎪⎪⎭⎫⎝⎛ππρ将上式右端展为傅里叶余弦级数，则有()dz L zn z k q L n LI n LA ππρcos 2000⋅-=⎰L L z n z n L L z n n L Ln L k I q 02000si n c o s 2-⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛⋅+=πππππρ()()[]112-2000--⎪⎭⎫ ⎝⎛=nn L L n L k Iq ππρ ()()()⎪⎩⎪⎨⎧+===12,402000220m n L n I n k Lq m n πρπ但不等于，()()()00,2-01000000=-==⎰I kL q zdz k q LI A L ()∑∞=+⋅⎪⎭⎫ ⎝⎛++⎪⎭⎫⎝⎛++=∴0002020012cos 12121242-m z L m L m I m L m I k L q k L q ππρπρπνν+=∴kqu四、以上问题的Matlab 的可视化：建立圆柱体：t=linspace(-pi,pi,200); y=linspace(0,4,200); [T,Y]=meshgrid(t,y); X=sin(T); Z=cos(T); mesh(X,Y,Z); axis equal由于本题中需要四维的数学建模（三维立体图形，再加一维的温度坐标），对于专业知识要求较高，以现有的知识水平还无法对四维进行分析求解，因此采用三维坐标（二维管道的横截面，再加一维的平面温度分布），进行求解。
传热学数值计算大作业
数值计算大作业一、用数值方法求解尺度为100mm×100mm 的二维矩形物体的稳态导热问题。
物体的导热系数λ为1.0w/m·K。
边界条件分别为： 1、上壁恒热流q=1000w/m2; 2、下壁温度t1=100℃； 3、右侧壁温度t2=0℃; 4、左侧壁与流体对流换热，流体温度tf=0℃，表面传热系数 h 分别为1w/m2·K、10 w/m2·K、100w/m2·K 和1000 w/m2·K;要求：1、写出问题的数学描述；2、写出内部节点和边界节点的差分方程；3、给出求解方法；4、编写计算程序(自选程序语言)；5、画出4个工况下的温度分布图及左、右、下三个边界的热流密度分布图；6、就一个工况下（自选）对不同网格数下的计算结果进行讨论；7、就一个工况下（自选）分别采用高斯迭代、高斯——赛德尔迭代及松弛法（亚松弛和超松弛）求解的收敛性（cpu 时间，迭代次数）进行讨论；8、对4个不同表面传热系数的计算结果进行分析和讨论。
9、自选一种商业软件（fluent 、ansys 等）对问题进行分析，并与自己编程计算结果进行比较验证（一个工况）。
（自选项）1、写出问题的数学描述 设H=0.1m微分方程 22220t tx y∂∂+=∂∂x=0，0<y<H ：()f th t t xλ∂-=-∂ 定解条件 x=H ，0<y<H ：t=t 2 y=0，0<x<H ：t=t1t 1t 2h ；t fq=1000 w/m 2y=H ，0<x<H ：tq yλ∂-=∂ 2、写出内部节点和边界节点的差分方程 内部节点：()()1,,1,,1,,122220m n m n m nm n m n m n t t t t t t x y -+-+-+-++=∆∆左边界： (),1,,1,1,,,022m n m n m n m nm n m n f m n t t t t t t x x h y t t y y y xλλλ-++---∆∆∆-+++∆=∆∆∆右边界： t m,n =t 2上边界： 1,,1,,,1,022m n m n m n m nm n m n t t t t t t y y q x x x x yλλλ-+----∆∆∆+++∆=∆∆∆ 下边界： t m,n =t 13、求解过程利用matlab 编写程序进行求解，先在matlab 中列出各物理量，然后列出内部节点和边界节点的差分方程，用高斯-赛德尔迭代法计算之后用matlab 画图。
大作业题(2014版)——中科院传热学(雁栖湖校区)
（注意：如发现作业中互相拷贝及抄袭者，此次作业成绩全将记为 一矩形平板 0 x a , 0 y b ，内有均匀恒 定热源 g 0 ， 在 x 0 及 y 0 处绝热， 在x a及y b 处保持温度 T1 ， 初始时刻温度为 T0 ， 如右图 1 所示： 1、求 t 0 时，矩形区域内的温度分布 T x, y, t 的解 析表达式； 2、若 a 18m ， b 12m ， g 0 1W m3 ， T1 600 K ， T0 200 K ，热传导系数 热扩散系数 0.8 m 2 s 。 请根据 1 中所求温度分布用 MATLAB k 1.0W m K ， 软件绘出下列结果，并加以详细物理分析： (a) 300s 内，在同一图中画出点 (0,4) 、 (0,8) 、 6,0 、 (12,0) 、 (9,6) （单位：m） 温度随时间的变化 (b) 200s 内，画出点 (18,4) 、 (18,8) 、 6,12 、 (12,12) 、 (9,6) （单位：m）处，分 别沿 x、y 方向热流密度值随时间的变化 (c) 画出 t 50s、 75s、 100s、 125s、 150s 时刻区域内的等温线 (d) 300s 内，在同一图中画出点 9,0 ， 9, 6 （单位：m）在其他参数不变时 g 0 分 别等于 1W m 3 ， 2W m3 ， 3W m3 情况下的温度、热流密度 (e) 600s 内，在同一图中画出点(9,6)（单位：m）在其它参数不变情况时热导率 分别为 0.5W m K 、 1.0W m K 和 1.5W m K 时的温度、热流密度 (f) 600s 内，在同一图中画出点(9,6) （单位：m）在其它参数不变情况时热扩散 系数分别为 0.4 m 2 s 、 0.8 m 2 s 和 1.2 m 2 s 时的温度、热流密度 3、运用有限差分法计算 2 中(b)、(d)和(e)，并与解析解结果进行比较，且需将数 值解与解析解的相对误差减小到万分之一以下； 4、附上源！ ）
传热学数值计算大作业
传热学数值计算大作业传热学数值计算大作业一选题《传热学》第四版P179页例题 4-3二相关数据及计算方法1．厚2δ=0.06m的无限大平板受对称冷却，故按一半厚度作为模型进行计算2. δ=0.03m，初始温度t0=100℃,流体温度t∞=0℃；λ=40W/(m.K),h=1000W/(m2.K),Bi=h*△x/λ=0.25；3.设定Fo=0.25和Fo=1两种情况通过C语言编程（源程序文件见附件）进行数值分析计算；当Fo=0.25时，Fo<1/(2*(1+Bi)),理论上出现正确的计算结果；当Fo=1时，Fo>1/(2*(1+Bi))，Fo>0.5,理论上温度分布出现振荡，与实际情况不符。
三网格划分将无限大平面的一半划分为6个控制体，共7个节点。
△x=0.03/N=0.03/6=0.005，即空间步长为0.005m四节点离散方程绝热边界节点即i=1时，tij+1=2Fo△ti+1j+(1-2Fo△)tij 内部节点即0tij+1=tij（1-2Fo△Bo△-2Fo△）+2Fo△ti-1j+2Fo△Bo△tf五温度分布线图（origin）六结果分析1 空间步长，时间步长对温度分布的影响空间步长和时间步长决定了Bo和Fo，两者越小计算结果越精确，但同时计算所需的时间就越长。
2 Fo数的大小对计算结果的影响编程时对Fo=1及0.25的情况分别进行了计算，发现当Fo=1时，各点温度随时间发生振荡，某点的温度高反而会使下一时刻的温度变低，违反了热力学第二定律，因此在计算中对Fo的选取有限制。
为了保证各项前的系数均为正值，对于内节点，Fo>0.5；对于对流边界节点，Fo<1/(2*(1+Bi))。
3 备注在Fo=0.25时，为了反映较长时间后温度的分布，取T=600,并选取了其中部分时刻的温度输出进行画图。
图像显示，随着时间的增长，各点温度趋向一致。
而当Fo=1时由于结果会出现振荡，只取T=6观察即可。
哈工程传热学数值计算大作业
传热学二维稳态导热问题的数值解法杨达文2011151419赵树明2011151427杨文晓2011151421吴鸿毅2011151416第一题：a=linspace(0,0.6,121);t1=[60+20*sin(pi*a/0.6)];t2=repmat(60,[80 121]);s=[t1;t2]; %构造矩阵for k=1:10000000 %理论最大迭代次数，想多大就设置多大S=s;for j=2:120for i=2:80S(i,j)=0.25*(S(i-1,j)+S(i+1,j)+S(i,j-1)+S(i,j+1));endendif norm(S-s)<0.0001break; %如果符合精度要求，提前结束迭代elses=S;endendS %输出数值解数值解数据量太大，这里就不打印出来，只画出温度分布。
画出温度分布：figure(1)xx=linspace(0,0.6,121);yy=linspace(0.4,0,81);[x,y]=meshgrid(xx,yy);surf(x,y,S)axis([0 0.6 0 0.4 60 80])grid onxlabel('L1')ylabel('L2')zlabel('t（温度）').60.66666777778L 1L 2t （温度）A0=[S(:,61)];for k=1:81B1(k)=A0(81-k+1);endB1 %x=L1/2时y方向的温度A1=[S(41,:)] %y=L2/2时x方向的温度x=0:0.005:0.6;y=0:0.005:0.4;A2=60+20*sin(pi*x/0.6)*((exp(pi*0.2/0.6)-exp(-pi*0.2/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6) )/2) %计算y=L2/2时x方向的解析温度B2=60+20*sin(pi*0.3/0.6)*((exp(pi*y/0.6)-exp(-pi*y/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6))/ 2) %计算x=L1/2时y方向的解析温度figure(2)subplot(2,2,1);plot(x,A1,'g-.',x,A2,'k:x'); %画出x=L1/2时y方向的温度场、画出x=L1/2时y方向的解析温度场曲线xlabel('L1');ylabel('t温度');title('y=L2/2');legend('数值解','解析解');subplot(2,2,2);plot(x,A1-A2); %画出具体温度场与解析温度场的差值曲线xlabel('L1');ylabel('差值');title('y=L2/2时，比较=数值解-解析解');subplot(2,2,3);plot(y,B1,'g-.',y,B2,'k:x'); %画出y=L2/2时x方向的温度场、画出y=L2/2时x方向的解析温度场曲线xlabel('L2');ylabel('t温度');title('x=L1/2');legend('数值解','解析解');subplot(2,2,4);plot(y,B1-B2); %画出具体温度场与解析温度场的差值曲线xlabel('L2');ylabel('差值');title('x=L1/2时，比较=数值解-解析解');y=L2/2时x方向的温度：60 60.1635347276130 60.3269574318083 60.4901561107239 60.653018915996160.8154342294146 60.9772907394204 61.1384775173935 61.298884093677961.4584005332920 61.6169175112734 61.7743263876045 61.930519281669662.0853891461909 62.2388298405943 62.3907362037523 62.541004126057762.6895306207746 62.8362138946214 62.9809534175351 63.123649991570263.2642058188844 63.4025245687647 63.5385114436490 63.672073244095163.8031184326565 63.9315571966177 64.0573015095482 64.180265191631864.3003639687311 64.4175155301449 64.5316395850212 64.642657917384664.7504944397430 64.8550752452343 64.9563286582797 65.054185283707565.1485780543131 65.2394422768254 65.3267156762441 65.410338438521565.4902532515567 65.5664053444751 65.6387425251668 65.707215216057165.7717764880854 65.8323820928694 65.8889904930310 65.941562890665265.9900632539310 66.0344583417471 66.0747177265744 66.110813815270166.1427218680003 66.1704200151959 66.1938892725421 66.213113553990066.2280796827826 66.2387774004857 66.2451993740203 66.247341200688866.2452014111934 66.2387814706441 66.2280857775556 66.213121660833566.1938993747528 66.1704320919304 66.1427358942990 66.110829762085766.0747355608048 66.0344780262737 65.9900847476605 65.941586148577365.8890154662295 65.8324087286383 65.7718047299493 65.707245003846265.6387737950858 65.5664380291767 65.4902872802189 65.410373736929465.3267521668755 65.2394798789402 65.1486166840471 65.054224854168964.9563690796505 64.8551164248743 64.7505362822981 64.642700324897664.5316824570463 64.4175587638655 64.3004074590802 64.180308831415964.0573451895733 63.9316008058186 63.8031618582281 63.672116371626463.5385541572596 63.4025667512431 63.2642473518283 63.123690755529062.9809932921539 62.8362527587866 62.6895683527611 62.541040603677462.3907713045038 62.2388634418130 62.0854211252013 61.930549515936761.7743547548873 61.6169438897778 61.4584248018242 61.298906131798361.1384972055701 60.9773079591820 60.8154488635041 60.653030848523060.4901652273162 60.3269636197632 60.1635378760476 60x=L1/2时y方向的温度：60 60.1308958471008 60.2618814819943 60.3930468323419 60.524481948785060.6562770664196 60.7885226663977 60.9213095376979 61.054728839108661.1888721614654 61.3238315901874 61.4596997681540 61.596569958966661.7345361106384 61.8736929197574 62.0141358961654 62.155961428198162.2992668485325 62.4441505006859 62.5907118062120 62.739051332642462.8892708622179 63.0414734614594 63.1957635516239 63.352246980097063.5110310927684 63.6722248074423 63.8359386883315 64.002285021688564.1713778926236 64.3433332631650 64.5182690516120 64.696305213238964.8775638224022 65.0621691561100 65.2502477791090 65.441928630549065.6373431122839 65.8366251788694 66.0399114293203 66.247341200688866.4590566635297 66.6752029193167 66.8959280998773 67.121383468913967.3517235256817 67.5871061108928 67.8276925149213 68.073647588380968.3251398551535 68.5823416279436 68.8454291264398 69.114582598162569.3899864420822 69.6718293350911 69.9603043614169 70.255609145064670.5579459853794 70.8675219958221 71.1845492460516 71.509244907413471.8418314019312 72.1825365549057 72.5315937512233 72.889242095483173.2557265760494 73.6312982331452 74.0162143310978 74.410738534857774.8151410909089 75.2296990126956 75.6546962706925 76.090423987246276.5371806363247 76.9952722483076 77.4650126199600 77.946723529732178.4407349585321 78.9473853161230 79.4670216732992 8066666666L 1t 温度y =L 2/2--1.--0.-3L 1差值y =L 2/2时，比较=数值解-解析解66778L 2t 温度x =L 1/200.050.10.150.20.250.30.350.4--1.--0.-3L 2差值x =L 1/2时，比较=数值解-解析解。
传热学数值计算大作业
传热学数值计算大作业航14 艾迪2011011537 如图所示，有一个正方形截面的无限长的水泥柱，热导率为，密度为，比热容为。
水泥柱的边长为。
水泥柱的左侧靠墙，可以认为保持温度为。
水泥柱被包围在温度为°的热空气中。
三个面上均只考虑对流换热，并且对流换热系数分别为，，。
请编写程序数值求解该稳态导热问题（可使用Fortran 或C 或Matlab 语言）。
作业要求提交源代码和报告，报告内容包括：（1） 给出该导热问题的数学描述； （2） 描述所采用的差分格式和求解过程；（3） 验证求解结果的准确性，给出网格无关性验证； （4） 给出求解结果（温度云图、边界热流、平均温度等）； （5） （选做）讨论对流换热系数、热导率等参数对求解结果的影响。
解：(1)、因为无内热源，温度分布：222201230(0,0)(x,0)t(0,y)t ,((x,0))(,y)(x,)((,y)),((x,H))f f f t tx H y H x ydt h t t dx dt H dt H h t H t h t t dx dxλλλ∂∂+=<<<<∂∂⎧=-=-⎪⎪⎨⎪-=--=-⎪⎩（2）、采用热平衡法建立内节点和边界节点的离散方程，x 、y 方向各取n 个节点，即()()11n n -⨯- 个网格，且x y ∆=∆ 。
对于任意内节点（i ，j ），有:,1,1,,1,1t (t t t t )/4i j i j i j i j i j -+-+=+++D边界三边界一边界节点：边界1、 1,0(1j )j t t n =≤≤边界2、11,1,21,11,1h 2h 2(2)t 2t t t (1k n)k k k k f xxt λλ-+∆∆+=+++<<边界3、22,k n 1,k n,k 1,k 1h 2h 2(2)t 2t t t (1k n)n n f xxt λλ--+∆∆+=+++<<边界4、33k,n k,n 11,n k 1,n h 2h 2(2)t 2t t t (1k n)k fxxt λλ--+∆∆+=+++<<C 点、2121n,1n 1,1n,2(h h )(h )(2)t t t f xh xt λλ-+∆+∆+=++D 点、2323n,n ,n 11,n (h h )(h )(2)t t t n n f xh xt λλ--+∆+∆+=++(3)、由于各个节点都写成了差分显示表达，可用高斯—赛德尔迭代法求解。
铜芯电缆温度分布MATLAB计算模型
1.题目如图1所示铜芯电缆，电流为5000A ，内径为10mm ，外包材料聚氯乙烯的厚度为2mm ，导热系数为0.15+0.00013{t}K)W/(m ⋅。
电缆左半边为绝热边界条件，右半边为第三类边界条件，空气温度为20℃，绝缘层表面与环境间的复合表面传热系数为10K)W/(m 2⋅。
铜的电阻率为()[]2010t-a R R t +=，m 1075.180⋅Ω⨯=-R ，C /004.0︒=a ，t 的单位为摄氏度。
试通过数值方法求解温度分布。
图 12.编程计算2.1 控制方程根据题意，本题为二维稳态导热问题，其控制方程为：011=+⎪⎭⎫⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂S T r r r T r r r θλθλ 边界条件：πππθ2~23,2~0007.0==r ：()W f B T T h q -=23~2007.0ππθ==r ：0=B q其中：C 20︒=f T 。
2.2 方程离散为建立通用方程，考虑非稳态项的控制方程为：011=+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂S T r r r T r r r t T cθλθλρ 采用全隐格式，在t ∆时间内，对控制容积积分，整理后可得：b T a T a T a T a T a S S N N W W E E P P ++++=其中：()ee e E r r a λδθ/∆=，()ww w W r r a λδθ/∆=，()n n n N r r a λδθ/∆=，()ss s S r r a λδθ/∆=，V S a a a a a a P P S N W E P ∆-++++=0()tV c a P ∆∆=ρ0，00P P CT a V Sb +∆=，()θ∆∆+=∆r r r V s n 5.0采用通用表达式，各表达式如下表：表1 坐标及系数表达式2.3 边界条件处理对于北边界，采用附加源项法处理。
由于北边界(πππθ2~23,2~0007.0==r )为第三类边界条件，则最靠近边界的控制容积加入以下附加源项：()Bn f adC x h T V AS λδ//1,+∆=()Bn ad P x h V A S λδ//11,+∆-=其中：C 20︒=f T将附加源项加到相应控制容积后，再令相应的0N =a 。
计算传热学导热作业
题目：如图所示，A,B 为两种不同材料的物体，厚度相同，左边壁面温度为Tw1，右边为Tw2。
两物体的导热系数都是温度的函数()111T a T b λ=+，a,b可以自己设定，但是要满足12λλ ，图中壁面温度和长度值可以自己确定，上下面绝热。
问题解决在此设定左边壁面温度Tw1为100℃，右边壁面温度Tw2为0℃.设壁面长度为1。
物体导热系数的函数()111T a Tb λ=+中待定系数10.00001a =，10b =，220.000001,0a b ==。
上下面绝热。
本题中采用Matlab 进行程序编写 ，并进行画图处理。
在处理中利用一维矩阵将计算区域划分成100个节点。
并通过式子：()(),p p e e w w e ew we w ewp e w p p c p a T a T a T b A A a a x x a a a S A x b S A xλλδδ=++===+-∆=∆进行迭代运算。
在本题中没有内热源所以0pc SS ==，则0b =，p e w a a a =+。
而可以假设面积1ew p A A A ===，上式可简化为：()(),p p e e w w ewe w ewp e wa T a T a T a a x x a a a λλδδ=+===+界面上的当量导热系数采用算术平均法计算：()()()()eee p e e e x x x x δδλλλδδ+-⎡⎤⎡⎤⎢⎥⎢⎥=+⎢⎥⎢⎥⎣⎦⎣⎦网格独立性的判定：可以通过输入不同的节点数来进行网格独立性试验，在这里分别选网格数为10、50、100。
计算的图如下。
10个网格时：50个网格时：100个网格时：可见当划分到100个网格时，网格已经具有独立性。
在这个问题中具有两种导热材料，在材料分界面处，把它当作离散区域的分界面来处理。
由于本题中的边界条件是第一类边界条件，所以不需要特殊的边界处理。
最终得到的温度分布为：。
1_数值传热-传热学作业-matlab
L P Rn5. 一厚度为 250mm 无限大平壁，其导热系数λ=43+0.08t w/(m ·k)，平壁一侧温度为 250℃，另一侧温度为 46℃，试用数值方法确定平壁内的温度分布，并确定通过该平壁的热流密度。
解：（1） 建立模型本题属于非常物性，无热源的一维稳态导热过程将平壁沿厚度方向（x 方向）划分为 N 个均匀相等的间距。
节点布置如图所示12 3 4 N N+1本题给出 N=10。
（2） 通过热流法建立离散方程a) 内节点离散方程对节点 P(i)所代表的微元体，在 x 方向上与节点 P 相邻的节点分别为 L(i-1)和 R(i+1)。
由于节点之间的间距很小，可以认为相邻节点间的温度分布是线性的。
节点 P 所代表的网格单元与它周围各网格单元之间的导热量可根据傅里叶定律直接写为:ΦLP = Z LP其中t i -1 - t i ⨯1∆xΦRP = Z RPt i +1 - t i ⨯1∆xZ RPZ LP= 43 + 0.08⨯ (t i + t i +1 )2 = 43 + 0.08⨯ (t i -1 + t i )2对节点 P(i)所代表的微元写热平衡式，即可得节点 P(i)温度的离散方程ΦLP + ΦRP = 0Z LPt i -1 - t i ⨯1+ Z ∆x RPt i +1 - t i ⨯1 = 0 ∆xt = Z LP t i -1 + Z LP t i +1 Z LP + Z RPb) 边界节点离散方程由于本题的壁面温度属于第一类边界条件，因此t 1 = 250 t n +1 = 46c) 热流密度的计算Z∑Z RP(t i +1 - t i )q = ∆t =i =166i46, t i = 次数 m, t k=k+1k>m不收 停机（3） 编写流程图计算机程序中使用的变量标识符： i 节点的坐标变量 t(i) 节点 i 的温度tt 前一次算出的节点温度 k 迭代次数计算机程序中输入数据： n 沿 x 方向的网格划分数 e 控制迭代过程终止的误差 m 允许的最大迭代次数开始输入 n, e, k, t 1, t n+1, t iN=10, t 1=250, t n+1= 1, e=0.01, k=10000迭代 k=0tt i =t iNoYesNoYes打印i , q 打印“敛”（4）编写程序本题使用matlab 软件，所编写的程序如下：clear;clc;t=ones(1,11); % 设定各项初始值q=ones(1,11);t(1)=250;t(11)=46;e=0.01;k=1;while 1 % 迭代程序tt=t;for i=2:10a=43+0.08*(t(i-1)+t(i))/2;b=43+0.08*(t(i)+t(i+1))/2;t(i)=(a*t(i-1)+b*t(i+1))/(a+b);endk=k+1;ttt=abs(t(5)-tt(5));if(ttt<e)break;endendfor i=2:11 % 计算热流密度a=43+0.08*(t(i-1)+t(i))/2;q(i)=q(i-1)+a*(t(i)-t(i-1));endktq=q(i)/0.25w=0:25:250;v=w/25+1;y=t(v); %绘制温度分布图plot(w,y),xlabel('位置（mm）'),ylabel('温度（℃）')legend('平壁内的温度分布',0),grid（5）计算结果迭代次数k=75平壁温度分布见下表热流密度q= -4.47E+04 w/m2（6）温度分布图∂⎬ ⎝ ⎭ i i -1i 10. 一砖墙厚 240mm ，内、外表面的表面传热系数分别为 6.0 w/(m 2·k)和 15 w/(m 2·k)，墙体材料的导热系数λ=0.43 w/(m ·k)，密度ρ=1668kg/m 3，比热容 c=0.75kJ/(kg ·K)，室内空气温度保持不变为 20℃，室外空气温度周期性变化，中午 12 点温度最高为 3℃，晚上 12 点温度最低为-9℃，试用数值计算方法确定内、外墙壁面温度在一天内的变化。
传热学实验—-墙角matlab导热问题
二维导热物体温度场的数值模拟姓名小明学号 111111班级能动学院能动一、问题描述有一墙角模型，尺寸如图1所示，导热系数0.53W/(m·K)，墙角内外壁为第一类边界条件。
求解该模型的温度分布及导热量。
图1q=0二、计算原理根据热平衡法列出节点方程，各方向导入单元体的热量之和为零。
内节点和绝热边界点（图1点划线上的点）的方程形式不同。
图2 Array图2所示的内节点和绝热边界节点方程如下：内节点：)()()()(1,,1,,1,1,,1,=⎥⎦⎤⎢⎣⎡-+-+-+-••=+++-+-+x y t t x y t t y x t t y x t t j i j i j i j i j i j i j i j i W E S N ∆∆∆∆∆∆∆∆ΦΦΦΦλ绝热边界点：)(02)(2)(1,,1,1,,1,=⎥⎦⎤⎢⎣⎡-++-+-••=+++--+x y t t y xt t y x t t j i j i j i j i j i j i W E S N ∆∆∆∆∆∆ΦΦΦΦλ三、计算过程用Matlab7.1语言编写计算程序，初取网格步长m y x 1.0=∆=∆运行结果：图1：各个点的温度数值图2：分层设色等温线分布图3：等温线分布（每两条线间隔为三度）四、小结本次数值模拟是运用matlab程序用于数值计算。
小组成员共同讨论并复习了热传导问题的数学描述和热平衡法；从模拟过程中练习了不同节点迭代方程的建立；并简单学习了matlab语言的使用。
这次大作业对于我们以后的学习和可能的研究来说是一个很好的锻炼机会。
传热学MATLAB温度分布大作业完整版
传热学大作业（第四章）姓名：张宝琪学号：03110608一、题目及要求1.各节点的离散化的代数方程2.源程序3.不同初值时的收敛快慢4.上下边界的热流量（λ=1W/(m℃））5.计算结果的等温线图6.计算小结题目：已知条件如下图所示：二、方程及程序(1)各温度节点的代数方程ta=(300+b+e)/4 ; tb=(200+a+c+f)/4; tc=(200+b+d+g)/4; td=(2*c+200+h)/4 te=(100+a+f+i)/4; tf=(b+e+g+j)/4; tg=(c+f+h+k)/4 ; th=(2*g+d+l)/4ti=(100+e+m+j)/4; tj=(f+i+k+n)/4; tk=(g+j+l+o)/4; tl=(2*k+h+q)/4tm=(2*i+300+n)/24; tn=(2*j+m+p+200)/24; to=(2*k+p+n+200)/24; tp=(l+o+100)/12 (2)源程序【G-S迭代程序】【方法一】函数文件为：function [y,n]=gauseidel(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=G*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0;0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0;0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0;0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0;0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]';[x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6) xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.1523 84.1429 67.9096 63.3793 62.4214 20.1557 15.4521 14.8744 14.7746 【方法2】>> t=zeros(5,5);t(1,1)=100;t(1,2)=100;t(1,3)=100;t(1,4)=100;t(1,5)=100;t(2,1)=200;t(3,1)=200;t(4,1)=200;t(5,1)=200;for i=1:10t(2,2)=(300+t(3,2)+t(2,3))/4 ;t(3,2)=(200+t(2,2)+t(4,2)+t(3,3))/4;t(4,2)=(200+t(3,2)+t(5,2)+t(4,3))/4;t(5,2)=(2*t(4,2)+200+t(5,3))/4;t(2,3)=(100+t(2,2)+t(3,3)+t(2,4))/4;t(3,3)=(t(3,2)+t(2,3)+t(4,3)+t(3,4))/4; t(4,3)=(t(4,2)+t(3,3)+t(5,3)+t(4,4))/4; t(5,3)=(2*t(4,3)+t(5,2)+t(5,4))/4;t(2,4)=(100+t(2,3)+t(2,5)+t(3,4))/4;t(3,4)=(t(3,3)+t(2,4)+t(4,4)+t(3,5))/4;t(4,4)=(t(4,3)+t(4,5)+t(3,4)+t(5,4))/4;t(5,4)=(2*t(4,4)+t(5,3)+t(5,5))/4;t(2,5)=(2*t(2,4)+300+t(3,5))/24;t(3,5)=(2*t(3,4)+t(2,5)+t(4,5)+200)/24;t(4,5)=(2*t(4,4)+t(3,5)+t(5,5)+200)/24;t(5,5)=(t(5,4)+t(4,5)+100)/12;t'endcontour(t',50);ans =100.0000 200.0000 200.0000 200.0000 200.0000 100.0000 136.8905 146.9674 149.8587 150.7444 100.0000 102.3012 103.2880 103.8632 104.3496 100.0000 70.6264 61.9465 59.8018 59.6008 100.0000 19.0033 14.8903 14.5393 14.5117【Jacobi迭代程序】函数文件为：function [y,n]=jacobi(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=B*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0; 0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0; 0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0; 0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0; 0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]'; [x,n]=jacobi(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6); xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)n =97Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.152384.1429 67.9096 63.3793 62.421420.1557 15.4521 14.8744 14.7746三、不同初值时的收敛快慢1、[方法1]在Gauss 迭代和Jacobi 迭代中，本程序应用的收敛条件均为norm(y-x0)>=eps ，即使前后所求误差达到e 的-6次方时，跳出循环得出结果。
哈尔滨工程大学传热学大作业数值计算matlab程序内容
传热学作业数值计算1 程序内容： 数值计算matlab程序内容：>> tw1=10; % 赋初值赋初值 tw2=20; c=1.5; p2=20; p1=c*p2; L2=40; L1=c*L2; deltaX=L2/p2; a=p2+1; b=p1+1; =ones(a,b)*5; m1=ones(a,b); m1(a,2:b-1)=zeros(1,b-2); m1(2:a,1)=zeros(a-1,1); m1(2:a,b)=zeros(a-1,1); m1(1,:)=ones(1,b)*2; k=0; max1=1.0; tn= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (i,j); switch m case 0 tn(i,j)=tw1; case 1 tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j)); case 2 tn(i,j)=tw1+tw2*sin(pi*(j-1)/(b-1)); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k max1 t2=ones(a,b); %求解析温度场求解析温度场for i=a:-1:1 for j=1:1:b y=deltaX*(a-i); x=deltaX*(j-1); t2(i,j)=tw1+tw2*sin(pi*x/L1)*(sinh(pi*y/L1))/(sinh(pi*L2/L1)); end end t2 迭代次数k =706 数值解温度场 数值解每次迭代的最大误差max1 =9.8531e-07 解析温度场 t2 解析温度场取第11行的解析解和数值解的点行的解析解和数值解的点行的解析解的直线，散点为其数值解的点 曲线为第11行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差数值计算matlab 程序内容：程序内容： >> tw1=10; tw2=20; c=1.5; p2=20; p1=c*p2; L2=20; deltaX=L2/p2; L1=c*L2; a=p2+1; b=p1+1; =ones(a,b)*5; m1=ones(a,b); m1(a,2:b-1)=zeros(1,b-2); m1(2:a,1)=zeros(a-1,1); m1(2:a,b)=zeros(a-1,1); m1(1,:)=ones(1,b)*2; k=0; max1=1.0; tn= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (i,j); switch m case 0 tn(i,j)=tw1; case 1 tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j)); case 2 tn(i,j)=tw2; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k max1 tx=ones(a,b); for i=1:1:a for j=1:1:b y=(a-i)*deltaX; x=(j-1)*deltaX; m=sym('m'); g=(((-1)^(m+1)+1)/m)*sin(m*pi*x/L1)*sinh(m*pi*y/L1)/sinh(m*pi*L2/L1); h=symsum(g,m,1,100); tx(i,j)=2*h*(tw2-tw1)/pi+tw1; end end tx 迭代次数k = 695 数值解温度场 数值解每次迭代的最大误差max1 =9.8243e-07 解析温度场 tx = 解析温度场行的解析解和数值解的点取第11行的解析解和数值解的点行的解析解的直线，散点为其数值解的点 曲线为第11行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像： 数值温度场 图像图像 ：解析温度场tx图像：数值解与解析解的误差数值解与解析解的误差程序内容： 数值计算matlab程序内容：>> t0=90; =10; L=10; c=0.25; p2=20; p1=p2/c; B=c*L; d=0.5*B; h=10; a=p2+1; b=p1+1; deltaX=B/p2; lambda=160; Bi=h*deltaX/lambda; =ones(a,b)*10; m1=ones(a,b)*3; m1(2:a-1,1)=zeros(a-2,1); m1(a,2:b-1)=ones(1,b-2); m1(1,2:b-1)=ones(1,b-2)*6; m1(2:a-1,b)=ones(a-2,1)*2; m1(1,b)=ones(1,1)*4; m1(a,b)=ones(1,1)*5; m1(1,1)=7; m1(a,1)=8; tn= ; max1=1.0; k=0; while ( max1>1e-6) k=k+1; max1=0; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=t0; case 1 tn(i,j)=(2*tn(i-1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 2 tn(i,j)=(2*tn(i,j-1)+tn(i-1,j)+tn(i+1,j)-4* )/(4+2*Bi)+ ; case 3 tn(i,j)=0.25*(tn(i,j-1)+tn(i,j+1)+tn(i-1,j)+tn(i+1,j)); case 4 tn(i,j)=(tn(i,j-1)+tn(i+1,j)-2* )/(2*Bi+2)+ ; case 5 tn(i,j)=(tn(i,j-1)+tn(i-1,j)-2* )/(2*Bi+2)+ ; case 6 tn(i,j)=(2*tn(i+1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 7 tn(i,j)=t0; case 8 tn(i,j)=t0; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k ta=ones(a,b); Bi1=h*d/lambda; sbi=sqrt(Bi1); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end x=deltaX*(j-1); ta(i,j)=(cosh(sbi*(L-x)/d)+sbi*sinh(sbi*(L-x)/d))*(t0- )/(cosh(sbi*L/d)+sbi*sinh(sbi*L/d))+ ; end end ta 迭代次数k =1461 数值解温度场 解析温度场 ta 解析温度场行的解析解和数值解的点取第11行的解析解和数值解的点曲线为第11行的解析解的直线，散点为其数值解的点行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像如下图像如下数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差程序内容：数值计算matlab程序内容：>> tw=10; L2=15; c=0.75; L1=L2/c; p2=24 ; p1=p2/c; deltaX=2*L2/p2; a=p2+1; b=p1+1; lambda=16; qv0=24; =ones(a,b)*5; m1=ones(a,b); m1(1,:)=zeros(1,b); m1(2:a,b)=zeros(a-1,1); m1(2:a,1)=zeros(a-1,1); m1(a,2:b-1)=zeros(1,b-2); tn= ; max1=1.0; k=0; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=tw; case 1 tn(i,j)=0.25*(tn(i-1,j)+tn(i+1,j)+tn(i,j-1)+tn(i,j+1)+qv0*(deltaX^2)/lambda); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k; tx=ones(a,b); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end if j>(b+1)/2 x=(j-(b+1)/2)*deltaX; else x=-((b+1)/2-j)*deltaX; end m=sym('m'); xi=(2*m-1)*pi/2; g=((-1)^m)/(xi^3)*(cosh(xi*y/L1)/cosh(xi*L2/L1))*cos(xi*x/L1); h=symsum(g,m,1,100); tx(i,j)=2*qv0*L1^2/lambda*h+qv0*(L1^2-x^2)/(2*lambda)+tw; end end tx 数值温度场 解析温度场tx 取第13行的解析解和数值解的点行的解析解和数值解的点行的解析解的直线，散点为其数值解的点曲线为第13行的解析解的直线，散点为其数值解的点解析解 第13行的误差=[数值解(13行) –解析解(13行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差。
平壁温度及热流计算-传热学大作业
传热大作业一一、网格的划分空间：网格在沿壁厚方向共划分为10份，编号为1~11共11个节点；时间：每60s划分一个网格，取最大时间为90000s，编号为1~1501共1501个节点。
故共划分出1500×10个网格。
二、节点方程组k为时间节点变量，i为空间节点变量1.内部节点：t(k+1,i)=Fo*(t(k,i-1)+t(k,i+1))+(1-2*Fo)*t(k,i)；2.左侧（高温侧）边界节点：t(k+1,1)=2*Fo*(t(k,2)+Bi1*tf1)+(1-2*Bi1*Fo-2*Fo)*t(k,1)3.右侧（低温侧）边界节点：t(k+1,n)=2*Fo*(t(k,n-1)+Bi2*tf2)+(1-2*Bi2*Fo-2*Fo)*t(k,n)三、计算框图见下页。
开始输入n,np,tmax,dt,dx,t0,tf1,tf2,a,b,h1,h2计算Fo，Bi1，Bi2Fo>1/(2*Bi1+2)or输出“不稳定”Fo>1/(2*Bi2+2)结束k*dt<=tmax计算内部节点和边界节点的温度，k=k+1四、程序代码（注：用matlab编译）n=11;%沿壁厚方向划分网格np=6;%打印节点温度的时间间隔数tmax=90000;%终止计算的时间dt=60;%时间间隔dx=0.01;%沿壁厚划分网格的间隔k=tmax/dt;%时间变量t0=5;%初始温度tf1=50;%左侧流体温度tf2=5;%右侧流体温度a=0.3437*10^-6;%热扩散率b=0.43;%导热系数h1=11;%左侧对流换热系数h2=23;%右侧对流换热系数Bi1=h1*dx/b;Bi2=h2*dx/b;Fo=a*dt/(dx^2);if Fo>1/(2*Bi1+2) %判断是否稳定disp (‘不稳定’)elseif Fo>1/(2*Bi2+2) %判断是否稳定disp( ‘不稳定’)elset=ones(k,n);%定义一个k行n列的二维矩阵存放温度计算数据t1=ones(k/np,n);%定义一个矩阵用于存储输出的各节点温度q=ones(k/np,2);%定义一个矩阵用于存储输出的左右侧面热流密度j=1;%计数变量，用于控制打印各量k=1;%初始化kt(j,:)=t0*t(k,:);%给温度矩阵第一行赋值初始温度t1(j,:)=t(1,:);%给输出温度矩阵第一行赋值初始温度q(j,1)=h1*(tf1-t(k,1));%计算左侧面初始热流密度q(j,2)=h2*(t(k,n)-tf2);%计算右侧面初始热流密度j=j+1;while k*dt<=tmax %判断是否超出终止计算时间for i=2:(n-1)t(k+1,i)=Fo*(t(k,i-1)+t(k,i+1))+(1-2*Fo)*t(k,i); %计算内部节点温度endt(k+1,1)=2*Fo*(t(k,2)+Bi1*tf1)+(1-2*Bi1*Fo-2*Fo)*t(k,1);%计算左侧边界节点温度t(k+1,n)=2*Fo*(t(k,n-1)+Bi2*tf2)+(1-2*Bi2*Fo-2*Fo)*t(k,n);%计算右侧边界节点温度if mod(k,np)==0t1(j,:)=t(k,:);%输出各节点温度q(j,1)=h1*(tf1-t(k,1));%输出左侧面热流密度q(j,2)=h2*(t(k,n)-tf2);%输出右侧面热流密度j=j+1;%下一个输出点endk=k+1;%进入下一个时间节点endend五、计算结果分别调用mesh和plot函数，得出曲线图，结果如下。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
东南大学能源与环境学院课程作业报告作业名称：传热学大作业——利用matlab程序解决热传导问题院系：能源与环境学院专业：建筑环境与设备工程学号：姓名：2014年11月9日一、题目及要求1.原始题目及要求2.各节点的离散化的代数方程3.源程序4.不同初值时的收敛快慢5.上下边界的热流量（λ=1W/(m℃））6.计算结果的等温线图7.计算小结题目：已知条件如下图所示：二、各节点的离散化的代数方程各温度节点的代数方程ta=(300+b+e)/4 ; tb=(200+a+c+f)/4; tc=(200+b+d+g)/4; td=(2*c+200+h)/4 te=(100+a+f+i)/4; tf=(b+e+g+j)/4; tg=(c+f+h+k)/4 ; th=(2*g+d+l)/4ti=(100+e+m+j)/4; tj=(f+i+k+n)/4; tk=(g+j+l+o)/4; tl=(2*k+h+q)/4tm=(2*i+300+n)/24; tn=(2*j+m+p+200)/24; to=(2*k+p+n+200)/24; tp=(l+o+100)/12 三、源程序【G-S迭代程序】【方法一】函数文件为：function [y,n]=gauseidel(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=G*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0;0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0;0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0;0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0;0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]';[x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6) xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.1523 84.1429 67.9096 63.3793 62.4214 20.1557 15.4521 14.8744 14.7746 【方法2】>> t=zeros(5,5);t(1,1)=100;t(1,2)=100;t(1,3)=100;t(1,4)=100;t(1,5)=100;t(2,1)=200;t(3,1)=200;t(4,1)=200;t(5,1)=200;for i=1:10t(2,2)=(300+t(3,2)+t(2,3))/4 ;t(3,2)=(200+t(2,2)+t(4,2)+t(3,3))/4;t(4,2)=(200+t(3,2)+t(5,2)+t(4,3))/4;t(5,2)=(2*t(4,2)+200+t(5,3))/4;t(2,3)=(100+t(2,2)+t(3,3)+t(2,4))/4;t(3,3)=(t(3,2)+t(2,3)+t(4,3)+t(3,4))/4; t(4,3)=(t(4,2)+t(3,3)+t(5,3)+t(4,4))/4; t(5,3)=(2*t(4,3)+t(5,2)+t(5,4))/4;t(2,4)=(100+t(2,3)+t(2,5)+t(3,4))/4;t(3,4)=(t(3,3)+t(2,4)+t(4,4)+t(3,5))/4;t(4,4)=(t(4,3)+t(4,5)+t(3,4)+t(5,4))/4;t(5,4)=(2*t(4,4)+t(5,3)+t(5,5))/4;t(2,5)=(2*t(2,4)+300+t(3,5))/24;t(3,5)=(2*t(3,4)+t(2,5)+t(4,5)+200)/24;t(4,5)=(2*t(4,4)+t(3,5)+t(5,5)+200)/24;t(5,5)=(t(5,4)+t(4,5)+100)/12;t'endcontour(t',50);ans =100.0000 200.0000 200.0000 200.0000 200.0000 100.0000 136.8905 146.9674 149.8587 150.7444 100.0000 102.3012 103.2880 103.8632 104.3496 100.0000 70.6264 61.9465 59.8018 59.6008 100.0000 19.0033 14.8903 14.5393 14.5117【Jacobi迭代程序】函数文件为：function [y,n]=jacobi(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=B*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0; 0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0; 0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0; 0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0; 0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]'; [x,n]=jacobi(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6); xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)n =97Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.152384.1429 67.9096 63.3793 62.421420.1557 15.4521 14.8744 14.7746四、不同初值时的收敛快慢1、[方法1]在Gauss 迭代和Jacobi 迭代中，本程序应用的收敛条件均为norm(y-x0)>=eps ，即使前后所求误差达到e 的-6次方时，跳出循环得出结果。
将误差改为0.01时，只需迭代25次，如下[x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',0.01)运行结果为 将误差改为0.1时，需迭代20次，可见随着迭代次数增加，误差减小，变化速度也在减小。
[方法2]通过 i=1:10判断收敛，为迭代10次，若改为1:20，则迭代20次。
2、在同样的误差要求下，误差控制在e 的-6次方内，Gauss 迭代用了49次达到要求，而Jacobi 迭代用了97次，可见，在迭代中尽量采用最新值，可以大幅度的减少迭代次数，迭代过程收敛快一些。
在Gauss 中，初值为100，迭代46次达到精确度1.0e-6,初值为50时，迭代47次，初值为0时，迭代49次，初值为200时迭代50次，可见存在一个最佳初始值，是迭代最快。
这一点在jacobi 迭代中表现的尤为明显。
五、上下边界的热流量：上边界t=200℃,∞t =10℃，所以，热流量Φ1=λ*[2*100-200x y ∆∆+x y a∆∆t -200+x y ∆∆b t -200+x y ∆∆c t -200+2*t -200d x y ∆∆] =1*（100/2+(200-139.6088)+(200-150.3312)+(200-153.0517)+(200-153.5639)/2) =230.2264W 下边界热流量Φ2=|λ*[x y ∆∆mi t -t +x y∆∆o j t -t +x y ∆∆p k t -t +2*t -t q l x y ∆∆]- h*(2*10-100x y ∆∆+x *t -t n ∆∆∞y +x *t -t o ∆∆∞y +x *t-t m ∆∆∞y +2*t -t p x y ∆∆∞)|=|1*((84.1429-20.1557)+(67.9096-15.4521)+(63.3793-14.8744)+(62.4214- 14.7746)/2)-10*(90/2+(20.1557-10)+(15.4521-10)+(14.8744-10)+(14.7746-10)/2)| = |-489.925|W =489.25W六、温度等值线Gauss:Yacobi:七、计算小结导热问题进行有限差分数值计算的基本思想是把在时间、空间上连续的温度场用有限个离散点温度的集合来代替，即有限点代替无限点，通过求解根据傅里叶定律和能量守恒两大法则建立关于控制面内这些节点温度值的代数方程，获得各个离散点上的温度值。