哈尔滨工程大学传热学大作业数值计算matlab程序内容

1. 表1 用三次样条方法插值计算0-90 度内整数点的sin 值和0-75 度内整数点的正切值，然后用5 次多项式拟合方法计算相同的函数。

a（度）0 15 30 45 60 75 90Sin（a）0 0.2588 0.5000 0.7071 0.8660 0.9659 1.0000tan（a）0 0.2679 0.5774 1.0000 1.7320 3.732解：分别对应的程序如下：正弦函数：x = pi*(0:90)/180;y = sin(x);xx = pi*(0:.25:90)/180;yy = spline(x,y,xx);plot(x,y,'o',xx,yy)正切函数：x = pi*(0:75)/180;y = tan(x);xx = pi*(0:.25:75)/180;yy = spline(x,y,xx);plot(x,y,'o',xx,yy)正弦拟合：figurex=pi*(0:15:90)/180;y=[0,0.2588,0.5,0.7071,0.866,0.9659,1.0]; xx=pi*(1:0.05:90)/180;p2=polyfit(x,y,5);yy=polyval(p2,xx);plot(x,y,'-ro',xx,yy);正切拟合：figurex=pi*(0:15:75)/180;y=[0,0.2679,0.5774,1,1.732,3.732];xx=pi*(1:0.05:75)/180;p2=polyfit(x,y,5);yy=polyval(p2,xx);plot(x,y,'-ro',xx,yy);legend('描点显示','五次拟合')2. 采用最近点法、线性法和3 次样条法插值计算1-100 整数间平方根n 1 4 9 16 25 36 49 64 81 100Sqtr（n）1 2 3 4 5 6 7 8 9 10解：程序如下：x=[1,4,9,16,25,36,49,64,81,100];y=[1,2,3,4,5,6,7,8,9,10];xx=1:100;yy=interp1(x,y,xx)subplot(2,2,1)plot(x,y,'-ro',xx,yy,'dr');title('线性法');subplot(2,2,2);y2=interp1(x,y,xx,'nearest');plot(x,y,'-ro',xx,y2,'dr');title('最近点法')subplot(2,2,3);y3=interp1(x,y,xx,'spline');plot(x,y,'-ro',xx,y3,'dr');title('3次样条法')仿真的结果：3. 已知p(x)=2x^4-3x^3+5x+13，求p（x）的全部根，由方程p（x）=0 的根构造一个多项式f（x），并和p（x）比较。

东南大学能源与环境学院课程作业报告课程名称：传热学作业名称：传热学大作业——利用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)区域离散化：如下图所示，用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。

数值计算软件作业1、利用mathematica 求下列函数的极限。

⑴求极限12)3131(lim +-+∞→-x x x 。

⑵求极限)ln (lim x x x x ++∞→。

⑶求极限x xx cot ln lim0→。

2、判别函数42()sin(1)f x x x x =+-在区间[-2, 2]上的单调性, 并给出单调区间。

3、求函数2369128)(xx x x f --+=的最大值、最小值、并画出函数的图形。

4、求曲线104--=x x y 与2062++=x x y 所围成图形的面积，并画出图形。

5、求解下列微分方程的数值解，并画出解函数的图形。

0)1(22=+'-+''y y y y ，1)0(,0)0(='=y y6、计算矩阵4124120235200117⎛⎫⎪⎪⎪⎪⎝⎭的行列式值、逆矩阵、特征值、特征向量、特征多项式。

7、以向量),1,1,1,1(1--=α),1,1,1,0(2-=α),1,1,0,0(3-=α)1,0,0,0(4=α为基，求向量)1,1,1,1(=β的坐标表达式。

8、绘制曲面222221, , 1x y z z x y x y z ++==+++=相交的空间图形。

9、咖啡馆配制两种饮料，甲种饮料每杯含奶粉9克、咖啡4克、糖3克，乙种饮料每杯含奶粉4克、咖啡5克、糖10克．已知每天原料的使用限额为奶粉3600克、咖啡2000克、糖3000克．如果甲种饮料每杯能获利0.7元，乙种饮料每杯能获利1.2元，每天在原料的使用限额内饮料能全部售出，每天应配制两种饮料各多少杯能获利最大？10、圆盘上有如下图所示的二十个数，请找出哪四个相邻数之和为最大，并指出它们的起始位置及最大和的值。

11、编写程序，求出能同时被2、3、5、7整除后余1的正整数，在1到10000以内的整数中有多少个？12、一副扑克牌有各种花色的牌各13张（不包括大小王），假设4个人玩牌，试编写程序，实现发牌的过程，使每家手里都有13张牌。

matlab传热计算程序
传热计算在工程学和科学领域中是一个重要的应用。

Matlab是一个功能强大的工程计算软件，可以用于传热计算。

在Matlab中，你可以使用各种方法来进行传热计算，比如有限元法、差分法、有限体积法等。

以下是一些常见的传热计算程序的示例：
1. 热传导方程求解，你可以编写一个Matlab程序来求解热传导方程，根据给定的边界条件和初始条件，使用差分法或有限元法来离散方程，并进行时间步进求解，得到温度场的分布。

2. 对流换热计算，对于流体内部的对流换热问题，你可以编写一个Matlab程序来求解Navier-Stokes方程和能量方程，结合有限体积法来进行流场和温度场的耦合求解。

3. 辐射换热计算，针对辐射换热问题，你可以编写一个Matlab程序来计算辐射传热，比如使用辐射传热方程和辐射传热模型，结合离散方法进行求解。

4. 传热系统优化，除了单一的传热计算，你还可以使用Matlab进行传热系统的优化设计，比如通过建立传热模型和耦合其
他工程模型，使用优化算法来寻找最优的传热系统设计参数。

总之，Matlab提供了丰富的工具和函数，可以用于传热计算的各个方面。

通过编写程序，你可以灵活地进行传热计算，并且可以根据具体的问题需求进行定制化的计算和分析。

希望这些信息对你有所帮助。

matlab大作业学号姓名：年级：专业：1、产生一个10 10的随机矩阵A，要求A中元素均为整数，范围[1,50]。

1)求出A中所有元素之和S，平均值M。

2)找到所有小于平均值，且能被3整除的元素。

3)绘制出A的二维纵向柱状图，横坐标为[8 5 9 1 2 3 4 7 10 13],条形宽度为0.7的“stacked”样式。

代码如下：clc,clear all,close allA=round(rand(10,10)*50);disp(A)S=sum(sum(A));P=mean(mean(A));disp(S)disp(P)disp('所有小于平均数且能被三整除的元素')XPS=H((mod(H,3)==0)&(H<P));disp(XPS')subplot(1,1,1),bar(A,0.7,'stacked'),title('ygh');set(gca,'XTickLabel',{'8','5','9','1','2','3','4','7','10','13'})2、产生一个随机四位密码。

用户用“input”进行输入对比。

猜错提示“WRONG”，正确提示“RIGHT”同时退出程序，最多五次机会。

代码如下：clc,clear all,close alldisp('请输入密码')A=round(8999*rand(1,1))+1000;m=1;while m<=5N=input('请输入一个四位数：');if A==N;disp('RIGHT');breakelsedisp('WRONG');endm=m+1;enddisp('密码是：')disp(A)disp('输入结束')3、按照脚本文件的编程风格，用for和while循环嵌套输出如下的乘法口诀表。

数值计算大作业一、用数值方法求解尺度为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 画图。

传热学二维稳态导热问题的数值解法杨达文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时，比较=数值解-解析解。

MATLAB 练习题1、已知 x=[2 3 5 6] ，y=-1:2 ，计算 z = x.^y ，解析 z 的值是如何计算的？程序：x=[2,3,5,6];y=-1:2;z=x.^y运行结果：z =0.5000 1.0000 5.0000 36.00002、设 1.2a =， 4.6b =-，8.0c =， 4.0e =-，d=3.2，计算22arctan e a bc t d ππ⎛⎫+ ⎪= ⎪ ⎪⎝⎭。

程序a=1.2;b=-4.6 ;c=8.0 ;e=-4.0;d=3.2;t=atan((2*pi*a+e/(2*pi*b*c))/d);t运行结果：t =1.17023、设矩阵311212123A ⎛⎫ ⎪= ⎪ ⎪⎝⎭，111210111B -⎛⎫ ⎪=- ⎪ ⎪-⎝⎭，求：（1）2A B +；（2）2243A B -；（3）AB ；（4）BA ；（5）AB BA -（1）程序A=[3,1,1;2,1,2;1,2,3];B=[1,1,-1;2,-1,0;1,-1,1];2*A+B运行结果：ans =7 3 16 1 43 3 7（2）程序：A=[3,1,1;2,1,2;1,2,3];B=[1,1,-1;2,-1,0;1,-1,1];4*A^2-3*B^2运行结果：ans =42 21 3840 19 4640 33 56 （3）程序：A=[3,1,1;2,1,2;1,2,3];B=[1,1,-1;2,-1,0;1,-1,1];A*B运行结果：ans =6 1 -26 -1 08 -4 2 （4）程序：A=[3,1,1;2,1,2;1,2,3];B=[1,1,-1;2,-1,0;1,-1,1];B*A运行结果：ans =4 0 04 1 02 2 2 （5）程序：A=[3,1,1;2,1,2;1,2,3];B=[1,1,-1;2,-1,0;1,-1,1];A*B-B*A运行结果：ans =2 1 -22 -2 06 -6 04、设矩阵21102041100A m ⎛⎫ ⎪- ⎪= ⎪ ⎪- ⎪⎝⎭，其中m 是你的学号的后四位，求1A -，A 的特征值和特征向量。

核科学与技术学院《传热学》二维稳态导热问题的数值解法作业姓名：罗晓学号：2014151214班级：20141512任课教师：李磊，张智刚哈尔滨工程大学核科学与技术学院2016年11月28日问题重述：第一题：如图所示，一个无限长矩形柱体，其横截面的边长分别为L 1和L 2，常物性。

该问题可视为二维稳态导热问题，边界条件如图中所示，其中L 1=0.6m ，L 2=0.4m ，T w 1=60℃，T w 2=20℃，λ=200W/(mK)。

1）编写程序求解二维导热方程。

2）绘制x =L 1/2和y=L 2/2处的温度场，并与解析解进行比较。

已知矩形内的温度场的解析解为：()()()()1211w2w1sh sh sin ,L L L y L x t t y x t πππ+=。

第二题将第一题中2y L =处的边界条件变为2w t t =，其他条件不变。

1)编写程序求解二维导热方程并计算从y =0处导入的热量2Φ。

2)当21L L 时，该二维导热问题可简化为一维导热问题。

在一维的近似下，试计算从y =0处导入的热量1Φ，并比较不同L 2/L 1下21ΦΦ的比值。

由该问题的解析解可知：L 2/L 10.0070.010.050.080.121ΦΦ0.99870.99120.9560.930.912解：（第一题第一问）对于此问题，由于可以视为二维稳态导热问题，由二维稳态热传导1,1,,,1,1,22220m n m n m nm n m n m nt t t t t t xy+-+-+-+-+=∆∆基本方程：22220t tx y∂∂+=∂∂用数值法对该区域进行节点划分(如下图所示）：x 方向上一共划分M 个节点，y 方向上一共划分N 个节点。

可以将以上方程改为：如果我们取x y ∆=∆则有：1,1,1,1,,4m n m n m n m nm n t t t t t +-+-+++=考虑到节点位置的特殊性，我们在此将节点的种类进行如下划分，并给出节点的离散方程。

5-2 matlab程序及结果分析说明：设边界条件中的，，共进行了15等分，产生了16个节点，要绘制的曲线为与的变化图线，程序中用y表示，用x来表示，用p来表示。

p=input('输入p=')x=0:1/15:1;Pe=15*p;y=(exp(Pe*x)-1)/(exp(Pe)-1);plot(x,y,'-*k') %精确解hold ony(1)=0,y(16)=1;for i=2:15y(i)=((1+0.5*p)*y(i-1)+(1-0.5*p)*y(i+1))/2;endplot(x,y(1:16),'-or') %中心差分hold onfor i=2:15y(i)=((1+p)*y(i-1)+y(i+1))/(2+p);endplot(x,y(1:16),'->g') %一阶迎风hold onfor i=2:15if p==1y(i)=((1+0.5*p)*y(i-1)+(1-0.5*p)*y(i+1))/2;elsey(i)=y(i-1)endendplot(x,y(1:16),'-+y') %混合格式hold onfor i=2:15if i==2y(i)=y(i+1)/(2+p)+(1+p)*y(i-1)/(2+p)+(p/(2+p))*(6*y(i)-3*y(i-1)-3*y(i+1))/8elsey(i)=y(i+1)/(2+p)+(1+p)*y(i-1)/(2+p)+(p/(2+p))*(5*y(i)-y(i-1)-y(i-2)-3*y(i+1))/8endendplot(x,y(1:16),'-<b')hold onlegend('精确解','中心差分','一阶迎风','混合格式','QUICK格式')运行结果见背面，由结果来看，数较小时，各种方式都与精确解比较接近，几种情况相比下，一阶迎风格式的准确度要差一些；较大时，中心差分方式会引起解的振荡，与准确解相差过大，几种格式相比，混合格式与准确解吻合较好。

传热学数值计算大作业航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)、由于各个节点都写成了差分显示表达，可用高斯—赛德尔迭代法求解。

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℃，试用数值计算方法确定内、外墙壁面温度在一天内的变化。

二维导热物体温度场的数值模拟姓名小明学号 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语言的使用。

这次大作业对于我们以后的学习和可能的研究来说是一个很好的锻炼机会。

传热学大作业（第四章）姓名：张宝琪学号：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次方时，跳出循环得出结果。

传热学数值计算大作业传热学是研究物体内部和之间热量传递的科学，其应用范围广泛，例如在工程领域中，传热学的数值计算被广泛用于优化热传递过程，提高能源利用效率。

本文将介绍传热学数值计算的大作业，主要内容包括问题陈述、计算方法和结果分析等。

问题陈述：本次大作业的问题是研究一个热管的热传递特性。

具体来说，热管由内外两个半圆形的金属管组成，内管壁与外管壁之间是一种导热的传热介质。

问题要求计算热管内外壁的温度分布，并分析传热过程的效率和优化热管的设计。

计算方法：计算热传递过程需要运用一些热传导定律和传热方程。

首先，根据Fourier 热传导定律，可得到内外壁的温度梯度。

然后，使用热传导方程来描述热传递过程，其中包括热扩散项和传热源项。

在计算热传导时需要注意材料的热导率、导热介质的热传导性质等参数。

在计算中，可以使用一些数值方法来离散化热传导方程，例如有限差分法、有限元法等。

其中，有限差分法是一种常见的数值方法。

通过将热传导方程中的导数用差分表达式替代，可以将偏微分方程转化为代数方程。

然后，可以使用迭代方法求解代数方程，得到温度分布的数值解。

结果分析：通过数值计算，可以得到热管内外壁的温度分布。

根据温度分布，可以分析热传递过程中的热流分布和传热效果。

例如，可以计算内外壁之间的热传导率，评估热管的热传递效率。

同时，可以对热管的设计进行优化。

例如，可以通过改变热导率高低、加大导热介质的厚度等方式，来提高热传递效果。

此外，对于热管的材料选择和导热介质的设计，还可以进行参数敏感性分析。

通过改变各个参数的数值，可以研究其对热传递过程的影响程度。

这有助于优化热管的设计，并提供一些实际应用方面的建议。

总结：传热学的数值计算是研究热传递现象的重要工具，可以帮助我们深入了解传热过程，优化传热装置的设计。

通过本次大作业，我们可以学习和练习传热学数值计算的方法和技巧，提升对传热现象的理解和分析能力。

希望通过这次大作业，能够更好地应用所学知识，解决实际问题。

实验一 MATLAB 基本操作及运算一、 实验目的1、 理解Matlab 数据对象的特点；2、 掌握基本Matlab 运算规则；3、 掌握Matlab 帮助的使用方法；二、 实验的设备及条件计算机一台（带有MATLAB7.0以上的软件环境）。

三、 实验内容要求建立一个名为experiment01.m 的，把与实验内容1-7相关的实验命令都放入该文件中，题与题之间用相应注释分割。

注意对实验中出现的相关函数或变量，请使用help 或doc 查询相关帮助文档，学习函数的用法。

1、 建立以下标量：1） a=102） b=2.5×10233） c=2+3i ，（i 为虚数单位）4） d=3/2πj e ,（j 为虚数单位，这里要用到exp ，pi ）2、 建立以下向量：1） aVec=[3.14 15 9 26]2） bVec=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡18228871.2 3） cVec=[5 4.8 … -4.8 -5 ] （向量中的数值从5到-5，步长为-0.2）4） dVec=[100 100.01 … 100.99 101] （产生1到10之间的等对数间隔向量，参考logspace ，注意向量的长度）3、 建立以下矩阵：1）⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=2222 aMat aMat 一个9×9的矩阵，其元素全为2；（参考ones 或zeros ）2）⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1000005000001 bMat bMat 是一个9×9的矩阵，除主对角上的元素为[1 2 3 4 5 4 3 2 1]外，其余元素均为0。

（参考diag ）。

3）10020109212291111=cMatcMat 为一个10×10的矩阵，可有1:100的向量来产生（参考reshape ）4）⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=NaN NaN NaN NaN NaN NaN NaN NaNNaN NaN NaN NaNdMatdMat 为3×4的NaN 矩阵，（参考nan ）5）⎥⎦⎤⎢⎣⎡---=8710225113eMat 6）产生一个5×3随机整数矩阵fMat ，其值的范围在-3到3之间。

MATLAB 大作业请各位同学，自己完成matlab 大作业的内容。

禁止相互抄袭，如有雷同，零分计算。

大作业的格式按照实验报告的格式书写，务必标明题号，作业完成后，将生成的报告打印出来提交。

正文的字号以宋体五号字，1.5倍行距的格式打印。

请与18周五前将大作业报告由班级负责人统一收齐交给我，过期不收。

1、 试编写名为test01.m 的MATLAB 函数，用以计算下述的值：⎪⎩⎪⎨⎧-<->=t t n t t t n t f 的对所有其他情况的对所有)4/sin()(si )4/sin()sin()4/sin()(si )4/sin()(ππππ要绘制t 关于函数f (t )的图形，其中t 的取值范围为ππ66≤≤-t ，间距为10/π。

（提示：注意要产生一系列的点，这里可考虑t 的输入是向量形式，可以利用find 函数找出所需限定值的元素的位置，对其按需要赋值后，再进行绘图；其次，另外一种思路，也可考虑使用循环的形式来实现）2、 编写函数，在同一窗口的4个子图中利用plot 等语句绘制y=at 2图像，其中a=[1 2 5 10],t 范围[-2，5]。

3、 求函数32)(3-+=x x x f 在区间[-5,5]上的最大值和最小值。

4、 求解函数⎰102dx e x 的数值积分和符号积分，并比较结果。

5、 求解微分方程3|;1|2)1(002='='=''+==x x y y y x y x 的精确解和解析解，并绘制图形。

假设求解区间为[0,10] 。

6、 说说你对MATLAB 及应用这门课程学习后的体会，另外请说明在所学章节中哪一章的内容你最感兴趣，为什么？哪一章的内容你认为是没有必要学习的，为什么？如果可以选择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函数，得出曲线图，结果如下。