传热学数值计算大作业2014011673

合集下载

传热学大作业

传热学大作业——二维物体热传导问题的数值解法1.二维热传导问题的物理描述：本次需要解决的问题是结合给定的边界条件，通过二维导热物体的数值解法，求解出某建筑物墙角稳态下的温度分布t以及单位长度壁面上的热流量φ。

1.1关于边界条件和研究对象选取的物理描述：如图所示为本次作业需要求解的建筑物墙壁的截面。

尺寸如图中所标注。

1.2由于墙角的对称性，A-A，B-B截面都是绝热面，并且由于对称性，我们只需要研究墙角的1/4即可(图中阴影部分)。

假设在垂直纸面方向上不存在热量的传递，我们只需要对墙角进行二维问题的研究即可。

1.3 关于导热量计算截面的物理描述：本次大作业需要解决对流边界条件和等温边界条件下两类边界条件的问题。

由于对称性，我们只需研究1/4墙角外表面和内表面的导热量再乘4，即是墙壁的总导热量。

2.二维热传导问题的数学描写：本次实验的墙角满足二维，稳态无内热源的条件，因此：壁面内满足导热微分方程：∂2t ∂x2+∂2t∂x2=0。

在绝热面处，满足边界条件：−λ(∂t∂n)=0。

在对流边界处满足边界条件：−λ(∂t∂n )x=x(x x−x x)3.二维热传导问题离散方程的建立：本次作业中墙角的温度场是一个稳态的连续的场。

本次作业中将1/4墙角的温度场离散化，划分成若干小的网格，每个网格的节点看成以它为中心的一个小区域的代表。

通过这些节点，采用“热平衡法”，建立起相应的离散方程，通过高斯-赛德尔迭代法，得到最终收敛的温度场，从而完成对墙角温度场的数值解。

对1/4墙角的网格划分如下：选取步长Δx=Δy=0.1m，为了方便研究，对导热物体的网格节点进行编码，编码规则如下：x,y坐标轴的方向如图所示，x,y轴的单位长度为步长Δx, 取左下角点为(1,1)点，其他点的标号为其在x,y轴上的坐标。

以此进行编码，进行离散方程的建立。

建立离散方程，要对导热物体中的节点根据其边界条件进行分类（特殊节点用阴影标出）：首先以对流边界条件下的墙角为例1.外壁面上，平直边界节点：建立离散方程：λΔy x x+1,x−x x,xΔx+λΔx2x x,x+1−x x,xΔy+λΔx2x x,x−1−x x,xΔy+hoΔx(x xx−x x,x)=0以(i,j)为中心节点，进一步整理得：x x,x=x2·(x x,x−1+x x,x+1)+x·x x+1,x+xx·Δx·x xx2x+xx·Δx2.外部角点：建立离散方程：ho·Δx(x xx−x x,x)+λΔy2x x,x+1−x x,xΔx+xΔx2·x x,x−1−x x,xΔx=0以(i,j)为中心节点，进一步整理得：x x,x=x2·(x x+1,x+x x,x−1)+xx·Δx·x xxx+xx·Δx3.绝热+对流边界角点：建立离散方程：ho·Δx2·(x xx−x x,x)+xΔx2·x x,x+1−x x,xΔx+xΔx2·x x+1,x−x x,xΔx=0以(i,j)为中心节点，进一步整理得：x x,x=x2·(x x,x+1+x x+1,x)+xx·Δx2·x xxx+xx·Δx24.内部角点：建立离散方程：hi·Δx·(x xx−x x,x)+x·Δx·x x,x+1−x x,xΔx+xΔx·x x−1,x−x x,xΔx+xΔx2·x x+1,x−x x,xΔx+xΔx2·x x,x−1−x x,xΔx=0以(i,j)为中心节点，进一步整理得：x x,x=x2·(x x+1,x+x x,x−1)+x(x x,x+1+x x−1,x)+xx·Δx·x xx3x+xx·Δx5.绝热平直边界节点：建立离散方程：x Δx2·x x,x+1−x x,xΔx+xΔx2·x x,x−1−x x,xΔx+xΔx·x x−1,x−x x,xΔx =0以(i,j)为中心节点，进一步整理得：x x,x=x2·(x x,x−1+x x,x+1)+x·x x−1,x2x6.对于普通内部节点：建立离散方程：xΔx·x x,x+1−x x,xΔx+xΔx·x x,x−1−x x,xΔx+xΔx·x x−1,x−x x,xΔx +xΔxx x+1,x−x x,xΔx=0以(i,j)为中心节点，进一步整理得：x x,x=x·(x x,x−1+x x,x+1+x x−1,x+x x+1,x)4x等温边界条件下：等温边界下内部节点和绝热边界下的节点离散方程与上述5，6式形式相同，在等温壁面处，节点方程只需写成x x,x=x x即可4.方程的求解：由上图可知，本题中有16*12=192个节点，相应地，就会有192个待求解的离散方程。

传热学习题答案(2014年10月15日)

习题答案（2014年10月15日）1. 解：采用热平衡法列出节点2、3、4的离散方程。

节点2：整理得同理节点3：节点4：肋端绝热，肋端对流，式中/3x H ∆=将题中已知条件代入可得方程组对于肋端绝热：联立求解得 t 2=92.2℃ t 3=87.7℃ t 4=86.2℃ 对于肋端对流：联立求解得 t 2=91.1℃ t 3=85.4℃ t 4=82.6℃2. 解：首先检验是否可用集总参数法。

为此计算B i V 数：3224/(4)(/)330.02524/()30.006060.033333/()V Rh R R hh V A Bi mW m K W m K ππλλλ⨯===⋅⨯==⋅＜可以用集总参数法。

22413324/()4(0.025)7.74107753/480/()(0.025)hA W m K m s cV kg m J kg K m πρ--⋅⨯⨯==⨯⨯⋅⨯ 据式（3-9）有4030030exp(7.7410)30t t t t τ-∞∞--==-⨯--℃℃450℃℃32122))2()0t t x h t t xxλδλδ∞--++∆⋅-=∆∆t （（t 1322(2)2()0t t t x h t t xλδ∞+-+∆⋅-=∆2433(2)2()0t t t x h t t xλδ∞+-+∆⋅-=∆344344()()0()()()0t t h x t t xt t h x t t xλδλδδ∞∞-+⋅∆-=∆-+⋅∆+-=∆322343432234342.045100.902.0450.901.02250.4502.045100.902.0450.901.03750.80t t t t t t t t t t t t t t -+=⎧⎪-++=⎨⎪-+=⎩-+=⎧⎪-++=⎨⎪-+=⎩由此解得τ=570 s=0.158 h3. 解：先判断本题能否利用集总参数法。

传热大作业-数值解法-清华-传热学

一维非稳态导热的数值解法一、导热问题数值解法的认识（一）背景所谓求解导热问题，就是对导热微分方程在规定的定解条件下的积分求解。

这样获得的解称为分析解。

近100年来，对大量几何形状及边界条件比较简单的问题获得了分析解。

但是，对于工程技术中遇到的许多几何形状或边界条件复杂的导热问题，由于数学上的困难目前还无法得出其分析解。

另一方面，在近几十年中，随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展十分迅速，并得到日益广泛的应用。

这些数值方法包括有限差分法、有限元法及边界元法等。

其中，有限差分法物理概念明确，实施方法简便，本次大作业即采用有限差分法。

（二）基本思想把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场，用有限个离散点上的值的集合来代替，将连续物理量场的求解问题转化为各离散点物理量的求解问题，将微分方程的求解问题转化为离散点被求物理量的代数方程的求解问题。

（三）基本步骤（1）建立控制方程及定解条件。

根据具体的物理模型，建立符合条件的导热微分方程和边界条件。

（2）区域离散化。

用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。

每一个节点都可以看成是以它为中心的一个小区域的代表，将小区域称之为元体。

（3）建立节点物理量的代数方程。

建立方法主要包括泰勒级数展开法和热平衡法。

（4）设立迭代初场。

（5）求解代数方程组。

（6）解的分析。

对于数值计算所获得的温度场及所需的一些其他物理量应作仔细分析，以获得定性或定量上的一些结论。

对于不符合实际情况的应作修正。

二、问题及求解（一）题目一厚度为0.1m 的无限大平壁，两侧均为对流换热边界条件，初始时两侧流体温度与壁内温度一致，1205f f t t t ===℃；已知两侧对流换热系数分别为h 1=11 W/m 2K 、h 2=23W/m 2K ，壁的导热系数λ=0.43W/mK ，导温系数a=0.3437×10-6 m 2/s 。

计算传热学大作业

计算传热学作业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");}。

大作业题(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观察即可。

传热学数值计算

传热学数值计算作业数值解程序：tw1=40 %三边温度tw2=100 %一边温度正弦变化幅度l1=40 %板长L1:40厘米l2=20 %板宽L2:20厘米m=41 %分划成40*20的网格n=21k=2dx=l1/(m-1)c=ones(n,m)for i=1:ma2(i)=tw1+tw2*sin(pi*dx*(i-1)/l1)c(1,i)=tw1 ,c(n,i)=a2(i)endfor j=1:nc(j,1)=tw1c(j,m)=tw1endwhile (abs(c(j,i)-k)>0.0001)k=c(j,i)for i=2:m-1for j=2:n-1c(j,i)=0.25*(c(j,i-1)+c(j,i+1)+c(j-1,i)+c(j+1,i)) endendend数值解中各网格点的温度值：数值二维温度分布图像：解析解程序： tw1=40 tw2=100 l1=40 l2=20 p=40 q=20 x(1)=0 for i=1:px(i+1)=x(i)+1 end y(1)=0 for j=1:qy(j+1)=y(j)+1 endfor i=1:p+1 for j=1:q+1n(j,i)=tw1+tw2*sinh(pi*y(j)/l1)*sin(pi*x(i)/l1)/sinh(pi*l2/l1) end end各网格点用解析式得到的温度值：50L1/cmnumerical calculation 2D temperature distributionL2/cmt e m p e r a t u r e /c e l s i u s d e g r e e解析二维温度分布图像：误差分析：取x=21，即位于板长一半处，温度随y （宽度）的变化曲线。

c1(:,1) 取自于数值解， c1(:,2) 取自于解析解 c1(:,1) c1(:,2) 40.0000 40.0000 43.3106 43.4164 46.6465 46.8538 50.0313 50.3335 53.4889 53.8771 57.0430 57.5062 60.7178 61.2434 64.5376 65.1117 68.5273 69.1350 72.7122 73.3381 77.1187 77.7470 81.7736 82.3888 86.7050 87.2922 91.9423 92.4875 97.5162 98.0068 103.4592 103.8840 109.8058 110.1555 116.5925 116.8600 123.8586 124.0388 131.6461 131.7363 140.0000 140.000050L1/cmanalytical method 2D temperature distributionL2/cmt e m p e r a t u r e /c e l s i u s d e g r e e误差曲线：由相对误差公式：d1= (c1(:,2) -c1(:,1))./ c1(:,2) 可得： d1 = 0 0.0024 0.0044 0.00600.00720.0081 0.0086 0.0088 0.0088 0.0085 0.0081 0.0075 0.0067 0.0059 0.0050 0.0041 0.0032 0.0023 0.0015 0.0007 0结论：数值解与解析解吻合很好。

传热学数值计算大作业

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

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

数值计算大作业
一、用数值方法求解尺度为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 t
x y
∂∂+=∂∂
x=0，0<y<H ：()f t
h t t x
λ
∂-=-∂ 定解条件 x=H ，0<y<H ：t=t 2 y=0，0<x<H ：t=t
1 t 1
t 2
h ；t f
q=1000 w/m 2
y=H ，0<x<H ：t
q y
λ
∂-=∂ 2、写出内部节点和边界节点的差分方程内部节点：
()
()
1,,1,,1,,1
2
2
220m n m n m n
m n m n m n t t t t t t x y -+-+-+-++
=∆∆
左边界：(),1,,1,1,,,022m n m n m n m n
m 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 n
m n m n t t t t t t y y q x x x x y
λ
λλ-+----∆∆∆+++∆=∆∆∆ 下边界： t m,n =t 1
3、求解过程
利用matlab 编写程序进行求解，先在matlab 中列出各物理量，然后列出内部节点和边界节点的差分方程，用高斯-赛德尔迭代法计算之后用matlab 画图。

代码（附件于文本档案中）使用100×100的网格数。

4、编写计算程序(自选程序语言) matlab 代码附附件文本档案中
5、画出4个工况下的温度分布图及左、右、下三个边界的热流密度分布图；当h=1时，温度分布图和热流密度分布图如下：
温度云图三维图：
此时，下边界平均热流密度为266.8284 W/m²，右边界平均热流密度为-462.8706 W/m²，左边界平均热流密度为-10.2150 W/m²，下边界平均热流密度为1000 W/m²。

物体平均温度为77.0639℃。

程序计算时间长度为41.8960s。

当h=10时，温度分布图和热流密度分布图如下：
温度云图三维图：
此时，下边界平均热流密度为263.9854 W/m²，右边界平均热流密度为-507.5412W/m²，左边界平均热流密度为-84.8137W/m²，下边界平均热流密度为1000 W/m²。

物体平均温度为78.4058℃。

程序计算时间长度为21.3020s。

当h=100时，温度分布图和热流密度分布图如下：
温度云图三维图为：
此时，下边界平均热流密度为304.8856 W/m²，右边界平均热流密度为-598.6677W/m²，左边界平均热流密度为-371.0580W/m²，下边界平均热流密度为1000 W/m²。

物体平均温度为77.0089℃。

程序计算时间长度为15.2290s。

当h=1000时，温度分布图和热流密度分布图如下：
温度云图三维图为：
此时，下边界平均热流密度为408.3915 W/m²，右边界平均热流密度为-623.0208W/m²，左边界平均热流密度为-605.3508W/m²，下边界平均热流密度为1000 W/m²。

物体平均温度为73.1784℃。

程序计算时间长度为10.8260s。

6、就一个工况下（自选）对不同网格数下的计算结果进行讨论；
当h=1000时：
使用50×50网格数，温度分布图和热流密度分布图如下：
温度云图三维图为
此时，下边界平均热流密度为340.6835 W/m²，右边界平均热流密度为-578.3104W/m²，左边界平均热流密度为-636.4214W/m²，下边界平均热流密度为
程序计算时间长度为2.7400s。

使用100×100网格数，温度分布图和热流密度分布图如下：
温度云图三维图为：
此时，下边界平均热流密度为408.3915 W/m²，右边界平均热流密度为-623.0208W/m²，左边界平均热流密度为-605.3508W/m²，下边界平均热流密度为
程序计算时间长度为10.8260s
结论：采用多数的网格数可以增加温度分布图的精确度，也使温度和热流密度计算得更加精准，使云图更加的准确，但是同样会使程序运算的时间加长。

7、就一个工况下（自选）分别采用高斯迭代、高斯——赛德尔迭代及松弛法（亚松弛和超松弛）求解的收敛性（cpu时间，迭代次数）进行讨论；
在h=1000工况下：
采用高斯迭代法，cpu时间为20.4491，迭代次数为1762次。

采用高斯——赛德尔迭代法，cpu时间为10.8260s，迭代次数为974次。

采用松弛法，cpu时间为7.9332s，迭代次数为761次。

斯迭代、高斯——赛德尔迭代及松弛法（亚松弛和超松弛）三种迭代法对应的迭代次数是逐渐减少的，也就是说三种方法在迭代的速度上是逐个增加的；
三种迭代法计算得到的解与严格计算方程组后的精确解在结果所示精度下是相同的，说明三种迭代法的求解精度是不低的。

8、对4个不同表面传热系数的计算结果进行分析和讨论。

当h=1时，下边界平均热流密度为266.8284 W/m²，右边界平均热流密度为-462.8706 W/m²，左边界平均热流密度为-10.2150 W/m²，下边界平均热流密度为1000 W/m²。