传热学数值计算大作业

合集下载

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

传热与流体数值计算大作业

大作业一、假设0,1x y≤≤的方腔内充满不可压缩流体，左、右、下壁面固定，上壁面以()22161u x x=--运动。

试求腔内的定常解。

（流体的物性取20℃的水。

同时，可以使用20℃的甘油作为对比）二、求解二维圆柱坐标中的Poisson-Nernst-Plack（PNP）方程，PNP方程来描述纳米孔内带电离子在浓度梯度及电场作用下的迁移行为和离子浓度分布。

具体方程如下所示：其中i=+/-，分别代表阴阳离子。

以及连续性方程:其中Φ是局域的电动势，c i表示i种离子的浓度，左侧边界上c+=10,c-=10，右侧边界上c+=1,c-=1。

j i表示离子流，D i为离子的扩散系数2×10-9，z i为离子的带电量,z i=1，T为溶液的温度,T=300。

e是电子电量1.602×10-19，ε0×εr=80，k B为波尔兹曼常数,k B=1.38×10-23。

边界上的电势Φ由高斯定律决定：对于带电的纳米孔壁（图中红色实线所示），有σs=σ（σ为纳米孔的表面电荷密度，数值为0.05）；对于其余区域有σs=0。

离子流j i在边界上的法向分量为零，即，求解φ、浓度c i以及ij的场。

（备注：求解区域为一圆柱形区域，长度为1200，直径为d=10。

建议步骤：可首先猜想浓度场c+和c-，并求解电动势场φ，通过连续性方程修正离子流场ij）大作业要求：1-3人为一组，完成以上任选一题目。

最终截止时间为12月26日。

在最终截止时间之前可以提交1次，若不满意得分可以继续修改。

大作业以报告形式提交，内容至少包括计算域的网格划分、方程的离散化、边界条件的处理、计算收敛的判据、计算的结果、结果的图形化显式、结果分析等。

源代码作为附录附在报告的最后。

数值传热_数值传热学大作业3gg

数值传热学 2009-2010 学年第一学期大作业 3
阶梯形标量场（长宽均为 1 ）的纯对流传递（ θ = 340 ），控制方程为 u ∂φ + v ∂φ = 0 ，上游的边界条件都是第一类的，即给定了φ 的分布，下游按开
∂x ∂y 口边界的方式处理。
图 1 示意图
（2）从图中可以得知：数值计算在剧烈变化区域（y=0.5 处），采用 CD、 SUD、QUICK 格式时，产生越界现象。
（3）计算过程中采用 STOIC 格式，计算效果最好。（4）采用不同的格式时，其表达式在 CBC 线内，则会出现稳定解，否则会出现解得越界现象，越界现象与对流稳定性不同。（5）编程过程中，采用 SOR 低松弛迭代方法，所得结果比较理想，计算速度远远 Gauss—Seidel 方法。在本次作业中，松弛系数取为 0.8。（6）编程过程中，要将边值点带入循环进行计算，否则会导致计算结果出错。进而也证明了，对流现象是具有方向性，其扰动只能沿下游传播。
要求：
1、分别利用 FUD,CD,SUD,QUICK,CLAM,EULER,MINMOD,MUSCL,OSHER,SMART, STOIC 来离散对流项，观察它们的计算结果有何不同。 2、用延迟修正进行求解。 3、写出详细的离散过程和求解方法。 4、用 TECPLOT 软件画出整个流场中φ 的分布，并画出 y = 0.5 时，φ 随 x 的分布。 5、编程采用 C/C++或 FORTRAN 语言。 6、将源程序附于作业之后，程序中要有详细的注释，以反映出思路。 7、源程序电子版和打印版上交时间：截止 2009 年 12 月 28 日。 8、每组交一份作业，给出每位同学的贡献度。(总和为 100%)

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

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

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

近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");}。

传热学数值计算大作业

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

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.8137 W/m²，下边界平均热流密度为1000 W/m²。

物体平均温度为78.4058℃。

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

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

物体平均温度为77.0089℃。

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

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

物体平均温度为73.1784℃。

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

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

物体平均温度为72.7267℃。

使用100×100网格数，温度分布图和热流密度分布图如下：
温度云图三维图为：
此时，下边界平均热流密度为408.3915 W/m²，右边界平均热流密度为-623.0208W/m²，左边界平均热流密度为-605.3508 W/m²，下边界平均热流密度为1000 W/m²。

物体平均温度为73.1784℃。

结论：采用多数的网格数可以增加温度分布图的精确度，也使温度和热流密度计算得更加精准，使云图更加的准确，但是同样会使程序运算的时间加长。

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²。

物体平均温度为77.0639℃。

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

当h=10时，下边界平均热流密度为263.9854 W/m²，右边界平均热流密度为-507.5412W/m²，左边界平均热流密度为-84.8137 W/m²，下边界平均热流密度为1000 W/m²。

物体平均温度为78.4058℃。

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

当h=100时，下边界平均热流密度为304.8856 W/m²，右边界平均热流密度为-598.6677W/m²，左边界平均热流密度为-371.0580 W/m²，下边界平均热流密度为1000 W/m²。

物体平均温度为77.0089℃。

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

当h=1000时，下边界平均热流密度为408.3915 W/m²，右边界平均热流密度为-623.0208W/m²，左边界平均热流密度为-605.3508 W/m²，下边界平均热流密度为1000 W/m²。

物体平均温度为73.1784℃。

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

传热学数值计算大作业

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

传热与流体数值计算大作业

数值传热_数值传热学大作业3gg

最新清华大学传热大作业-无限大平壁简单数值计算报告

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

计算传热学大作业

传热学数值计算大作业

传热学数值计算大作业