哈工程传热学数值计算大作业

传热学大作业班级:２０１21515 学号:201215１531 姓名：张永宽第一题:如图所示,一个无限长矩形柱体，其横截面的边长分别为L 1和L 2,常物性。

该问题可视为二维稳态导热问题,边界条件如图中所示,其中Ｌ１=０.６m,L 2=0.4m,T w １=60℃，T w ２=2０℃,λ=200W ／(m·K ）。

（１) 编写程序求解二维导热方程。

(2） 绘制x =Ｌ１／２和y =L 2/2处的温度场,并与解析解进行比较。

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

（１）根据课本164页公式（b ）Tm ，ｎ=(Tm ＋1，ｎ+T ｍ—1,n ＋Tm,ｎ+1+Ｔm ，n-1）/4;取步长为1cm 。

编出以下程序迭代求解内部个点温度。

a=ｚeros(41，６1）; ％生成４1＊60的矩阵. k=0:6０； ａ（4１，：）=２0＊s ｉn(pi.＊k/６０）；％矩形上边温度满足Tw2=sin （pi*x ／L １）。

a ＝ａ+６0; ％使四周都为给定的边界条件。

ｆor x ＝1：１00０0%迭代10000次（估计能满足要求精度）。

for i=2：４0for j ＝２：６0 a(i ，j)=（a （i —1,ｊ)+a （ｉ，ｊ－1）+a （ｉ+1，ｊ)＋a （i,j+1))/4； ％内部每一个点都为周围四个点温度和的四分之一。

en ｄ e ｎd ｅndmes ｈ（a ）t ｉtle(’第一题（张永宽作请勿抄袭）’,’Ｆonts ｉze ＇,18） x ｌabel('x 轴张永宽作请勿抄袭，单位cm ＇，＇Fontsi ｚｅ＇，１４) ylabe ｌ（’y 轴，单位ｃm ’,'F ｏntsize'，1４） zlabel(＇t 轴，单位℃'，＇Fon ｔsize ’,1４） 迭代一万次后个点温度数据：迭代法温度分布图：x 轴张永宽作请勿抄袭,单位cmy 轴,单位cmt 轴,单位℃（2)Y=L2/２时的温度曲线即把第一问中第21行数据画出图即可.x 轴,单位cmt 轴,单位℃Y ＝L2／2处的温度误差,即用第一问中行列式第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张牌。

传热学大作业班级：20121515 学号：2012151531 ：永宽第一题：如图所示，一个无限长矩形柱体，其横截面的边长分别为L 1和L 2，常物性。

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

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

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

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

（1）根据课本164页公式（b ）Tm ，n=（Tm+1,n+Tm-1,n+Tm,n+1+Tm,n-1）/4;取步长为1cm 。

编出以下程序迭代求解部个点温度。

a=zeros(41,61); %生成41*60的矩阵。

k=0:60;a(41,:)=20*sin(pi.*k/60);%矩形上边温度满足Tw2=sin(pi*x/L1). a=a+60; %使四周都为给定的边界条件。

for x=1:10000%迭代10000次（估计能满足要求精度）。

for i=2:40 for j=2:60a(i,j)=(a(i-1,j)+a(i,j-1)+a(i+1,j)+a(i,j+1))/4; %部每一个点都为周围四个点温度和的四分之一。

end end end mesh(a)title('第一题（永宽作请勿抄袭）','Fontsize',18) xlabel('x 轴永宽作请勿抄袭,单位cm','Fontsize',14) ylabel('y 轴,单位cm','Fontsize',14) zlabel('t 轴,单位℃','Fontsize',14) 迭代一万次后个点温度数据：迭代法温度分布图：（2）Y=L2/2时的温度曲线即把第一问中第21行数据画出图即可。

计算传热学作业1、 一块厚度为2h=200mm 的钢板，放入T f =1000℃的炉子中加热，两表面换热系数h=174W/(m 2.℃),钢板的导热系数k=34.8 W/(m. ℃),热扩散率a=5.55×10-6m 2/s,初始温度T i =20℃. 求温度场的数值解；分别用显示、C-N 、隐式 解： 1、数学模型该问题属于典型的一维非稳态导热问题。

由于钢板两面对称受热，板内温度分布必以其中心截面为对称面。

因此，只要研究厚度为δ的一半钢板即可。

将x 轴的原点置于板的中心截面上。

这一半钢板的非稳态导热的数学描述为2、计算区域离散化：该一维非稳态导热问题可当做二维问题处理，有时间坐标τ和空间坐标x 。

采用区域离散方法A ，将空间区域等分为m 个子区域，得到m+1个节点。

如下图所示，纵坐标为时间，从一个时到另一个时层的间隔即时间步长为∆t ,每个时层都会对下一时层产生影响。

空间与时间网格交点(i ，k )，代表了时空区域的一个节点，其温度为，离散方法如下图。

综合考虑计算效率同时保证数值计算格式的稳定性，本文取空间步长∆x =0.01m ，时间步长∆t =5s ，对半平板空间的离散共得到11个节点。

x TaT 22∂∂=∂∂τ==τT T 00==∂∂x xT δλ=-=∂∂-x T T h xT f )(图 时间-空间区域离散化3、离散方程组对于一维非稳态方程，扩散项采用中心差分，非稳态项取时间向前差分。

扩散项根据时层采用不同的处理方法，得到了三种格式的离散方程组，即显式、隐式、C-N 格式，等式左右分属不同的时层。

（1） 显示差分格式： 内部节点：()]][[]][1[]][[2]][1[]1][[2j i T j i T j i T j i T xt a j i T +-+*-+∆∆*=+左边界：]][0[21]][1[2]1][0[22j T x t a j T xt a j T ⎪⎭⎫⎝⎛∆∆**-+∆∆**-=+ 右边界：()f T j T x k t a h j T x t a j T xt a j T -∆*∆***+⎪⎭⎫ ⎝⎛∆∆**-+∆∆**-=+]][10[2]][10[21]][9[2]1][10[22（2） 隐式差分格式： 内部节点：]][[]1][1[]1][[21]][1[222j i T j i T x t a j i T x t a j i T x t a -=⎪⎪⎭⎫ ⎝⎛+-∆∆*++⎪⎭⎫⎝⎛∆∆**+-+∆∆* 左边界：]][0[]1][0[)21(]1][1[222j T j T xt a j T xt a -=+∆∆**+-+∆∆**右边界：]][10[2]1][9[)2]1][10[)21(2j T xk t h a j T xt a j T xk t h a +∆*∆***=+∆∆**++∆*∆***+（3）C-N 差分格式：内部节点：()]][1[]][[2]][1[2]][[]1][1[]1][[21]1][1[22222j i T j i T j i T x t a j i T j i T x t a j i T x t a j i T x t a -+-+∆*∆*--=⎪⎪⎭⎫ ⎝⎛+-∆∆*++⎪⎭⎫⎝⎛∆∆**+-++∆*∆*左边界：]][1[]][0[)1(]1][1[)]1][0[)1(222j T j T xt a j T xt a j T xt a -∆∆*--=+∆∆*++∆∆*--右边界：fT xk t h a j T xt a j T xt a xk t h a j T xt a j T xt a xk t h a ∆*∆***-∆∆*-∆∆*+∆*∆**--=+∆∆*++∆∆*-∆*∆**--2]][9[]][10[)1(]1][9[)]1][10[)1(22224、计算结果源程序代码： 显式：#include<stdio.h>#include<time.h> #include<cstdlib> #include<math.h> #include<stdlib.h> #include <process.h> double T[11][5000]; main()int i,j;double k;/*µ¼ÈÈϵÊý*/double h;/*»»ÈÈϵÊý*/double a;/*ÈÈÀ©É¢ÂÊ*/double x1,t1;/*x1±íʾλÖò½³¤£¬ti±íʾʱ¼ä²½³¤*/ double T0;/*T0±íʾ³õʼζÈ*/double Tf;/*Tf±íʾ¯ÎÂ*/double p,q;h=174;k=34.8;a=0.00000555;T0=20;Tf=1000;x1=0.01;t1=5;/*T[199][j]=(T[198][j]+h*x1*Tf/k)/(1+h*x1/k);*/for(i=0;i<=10;i++) T[i][0]=T0;for(j=0;j<4999;j++){ T[0][j+1]=2*a*t1*(T[1][j]-T[0][j])/(x1*x1)+T[0][j];for(i=1;i<10;i++){p=a*(T[i+1][j]-2*T[i][j]+T[i-1][j])/(x1*x1);/*q=(T[i][j+1]-T[i][j])/t1;q=p;*/T[i][j+1]=p*t1+T[i][j];}T[10][j+1]=2*h*a*t1*(Tf-T[10][j])/(x1*k)+2*a*t1*(T[9][j]-T[10][j])/(x1*x1)+T[10][j];}for(i=0;i<=10;i++){printf("%f",T[i][4999]);/*´òÓ¡Êä³ö*/printf("\n");}system("pause");}隐式：#include<stdio.h>#include<time.h>#include<cstdlib>#include<math.h>#include<stdlib.h>#include <process.h>double T[11][5000];main(){int i,j;double k;/*µ¼ÈÈϵÊý*/double h;/*»»ÈÈϵÊý*/double a;/*ÈÈÀ©É¢ÂÊ*/double x1,t1;/*x1±íʾλÖò½³¤£¬t1±íʾʱ¼ä²½³¤*/ double T0;/*T0±íʾ³õʼζÈ*/double Tf;/*Tf±íʾ¯ÎÂ*/double A[11],B[11],C[11],D[11],P[11],Q[11];h=174;k=34.8;a=0.00000555;T0=20;Tf=1000;x1=0.01;t1=5;for(i=0;i<=10;i++)T[i][0]=T0;for(j=1;j<=4999;j++){for(i=1;i<=9;i++) A[i]=a*t1/(x1*x1);A[0]=0;A[10]=2*a*t1/(x1*x1);for(i=0;i<=9;i++)B[i]=-(1+2*a*t1/(x1*x1));B[0]=-(1+2*a*t1/(x1*x1));B[10]=-(1+2*a*t1*h/(k*x1))-2*a*t1/(x1*x1);for(i=1;i<=9;i++)C[i]=a*t1/(x1*x1);C[0]=2*a*t1/(x1*x1);C[10]=0;for(i=0;i<=9;i++)D[i]=-T[i][j-1];D[10]=-2*a*t1*h*Tf/(k*x1)-T[10][j-1];for(i=1;i<=10;i++){A[i] = A[i] / B[i-1];B[i] = B[i] - C[i-1] * A[i];D[i] = D[i] - A[i] * D[i-1];}T[10][j] = D[10] / B[10];for(i=9;i>=0;i--)T[i][j] = (D[i] - C[i] * T[i+1][j]) / B[i];}for(i=0;i<=9;i++){printf("%f",T[i][4999]);/*´òÓ¡Êä³ö*/printf("\n");}system("pause");}C-N：#include<stdio.h>#include<time.h>#include<cstdlib>#include<math.h>#include<stdlib.h>#include <process.h>double T[11][5000];main(){int i,j;double k;/*µ¼ÈÈϵÊý*/double h;/*»»ÈÈϵÊý*/double a;/*ÈÈÀ©É¢ÂÊ*/double x1,t1;/*x1±íʾλÖò½³¤£¬t1±íʾʱ¼ä²½³¤*/double T0;/*T0±íʾ³õʼζÈ*/double Tf;/*Tf±íʾ¯ÎÂ*/double A[11],B[11],C[11],D[11],P[11],Q[11];h=174;k=34.8;a=0.00000555;T0=20;Tf=1000;x1=0.01;t1=5;for(i=0;i<=10;i++)T[i][0]=T0;for(j=1;j<=4999;j++){for(i=1;i<=9;i++) A[i]=a*t1/(2*x1*x1);A[0]=0;A[10]=a*t1/(x1*x1);for(i=0;i<=9;i++)B[i]=-(1+a*t1/(x1*x1));B[0]=-(1+a*t1/(x1*x1));B[10]=-(1+a*t1*h/(k*x1))-a*t1/(x1*x1);for(i=1;i<=9;i++)C[i]=a*t1/(2*x1*x1);C[0]=a*t1/(x1*x1);C[10]=0;for(i=1;i<=9;i++)D[i]=-T[i][j-1]-(a*t1/(2*x1*x1))*(T[i+1][j-1]-2*T[i][j-1]+T[i-1][j-1]);D[0]=(-1+a*t1/(x1*x1))*T[0][j-1]-(a*t1/(x1*x1))*T[1][j-1];D[10]=(-a*t1*h/(k*x1)-a*t1*h/(k*x1))*Tf+(-1+a*t1*h/(k*x1)+a*t1/(x1*x1))*T[10][j-1]-a*t1*T[9][j-1]/(x1*x1);for(i=1;i<=10;i++){A[i] = A[i] / B[i-1];B[i] = B[i] - C[i-1] * A[i];D[i] = D[i] - A[i] * D[i-1];}T[10][j] = D[10] / B[10];for(i=9;i>=0;i--)T[i][j] = (D[i] - C[i] * T[i+1][j]) / B[i];}for(i=0;i<=9;i++){printf("%f",T[i][4999]);/*´òÓ¡Êä³ö*/printf("\n");}system("pause");}。

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

传热学数值计算大作业传热学数值计算大作业一选题《传热学》第四版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时，比较=数值解-解析解。

不同居住建筑户间传热温差的探讨1110811005 崔晓蒙摘要: 本文采用非稳态传热数值求解法, 分析传统混住宅、钢筋混凝土住宅，外墙内保温住宅, 采用分室调节最低值时, 房间室温和与邻室的温差及围护结构温度分布的变化规律, 为住宅按户计算供暖热负荷,提供了计算依据, 也为按户计量收费和建筑保温，能源节约研究提供参考。

1.前言：随着我国能源结构的调整和供暖收费改革制度, 采用燃气、电能、太阳能、地热能、江河湖海等清洁能源的多种能源供暖方式代替传统燃煤锅炉房供暖已成为必然趋势, 这也是我国节约能源、保护环境, 实施可持续发展战略的重要措施之一。

无论采用哪种热源集中供暖, 按用热量收费必然促使用户按需分室调节。

在进行住宅按户计量供暖设计时的热负荷确定中, 用户分室调节导致的户间传热计算, 成为工程设计与研究部门关注的热点问题, 其中, 分室调节时室温和户间传热温差及围护结构温度分布的变化规律, 是户间传热计算的关键依据。

目前, 我国户间传热计算的关键性参数的确定尚处于探讨之中, 迫于工程设计急需, 采取先行界定、在过程中逐步完善的策略, 因而, 在一些地方的设计规程中对这些关键参数采取/暂定的方法。

目前我国户间传热的研究主要基于稳态传热分析, 而用户分室调节导致的户间传热, 是通过围护结构的非稳态传热。

本文分别以目前较为普遍的传统砖混住宅、钢筋混凝土住宅和外墙内保温结构的三种居住建筑为对象, 采用非稳态传热分析方法, 数值求解, 分别得出了三种住宅标准层和顶层的典型房间( 分别有一面外墙和有两面外墙) , 采用分室调节最低值时( 即停止采暖时) , 典型房间室温和与邻室的户间传热温差及围护结构温度分布的变化规律, 并进行分析比较。

2.房间非稳态传热的数学模型描述典型房间室温及围护结构温度变化规律的不稳态传热的数学模型, 由围护结构的非稳态导热微分方程式及其定解条件和房间热平衡方程组成。

2. 1 非稳态导热微分方程式围护结构的不稳态导热, 可视为平壁一维不稳态导热问题。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

计算传热学大作业报告戴平0708180209选题：题目1题目3题目1已知：一块厚度为0.1mm的无限大平板，具有均匀内热源，q＝50×103W/m3，，导热系数K＝10W/m.℃，一侧边界给定温度为75℃，另一侧对流换热，T f＝25℃，，h=50W/m2.℃,求解稳态分布。

（边界条件用差分代替微分和能量平衡法），画图。

(内,外节点)解：由题目分析可得此情况是有内热源的一维稳态导热问题。

采用均匀网格，外节点法，网格间距取为0.01mm，将无限大的平板沿其厚度方向均匀分为10份。

如图：1 2 3 4 5 6 7 8 9 10 11 网格的具体分布：当I=2,3,4，…9，10边界条件：当I=1时，当I=11时：控制微分方程为： 0d dTkq dx dx⎛⎫+= ⎪⎝⎭ 边界条件为：一边为第一类边界条件0075x T C==另一边为第三类边界条件T f ＝25℃，，h=50W/m 2.℃ 方程的离散化：对控制微分方程进行积分0e wdT dT k k xq dx dx ⎛⎫⎛⎫-+∆= ⎪ ⎪⎝⎭⎝⎭ 设相邻网格之间的温度是线性分布的，而导热系数k 是常数，所以得()()0P WE P e w T T T T k k xq x x δδ⎛⎫⎛⎫---+∆= ⎪ ⎪ ⎪⎪⎝⎭⎝⎭T1=75℃Tf=25℃h经整理得：1120.5I I I T T T -+=++ I=2,3，…，9,10）边界条件：左边是第一类边界条件，得：23275.5T T =+右边是第三类边界条件，在节点11的半个控制容积内对控制微分方程进行积分：10.511102q q q x -+∆=而由条件得：()1111f q h T T =- 且 101110.5T T q kxδ-=所以得：()10111112f T T q x kh T T xδ-∆+=- 整理得：11101.05 1.5T T =+从而可得三对角矩阵，利用TDMA 解法解之得到温度的稳态分布。

传热学作业数值计算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行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差。

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