利用matlab程序解决热传导问题

课程作业报告

作业名称：传热学大作业——利用matlab程序解决热传导问题

                       院系：能源与环境学院

 专业：建筑环境与设备工程

                       学号：*******

                       ******

                           2015年6月8日

一、题目及要求

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)/4

ti=(100+e+m+j)/4;  tj=(f+i+k+n)/4;  tk=(g+j+l+o)/4;  tl=(2*k+h+q)/4

tm=(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)>=eps

    x0=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.5639

  108.1040  108.61  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:10

t(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'

end

contour(t',50);

ans =

  100.0000  200.0000  200.0000  200.0000  200.0000

  100.0000  136.05  146.9674  149.8587  150.7444

  100.0000  102.3012  103.2880  103.8632  104.3496

  100.0000   70.62   61.9465   59.8018   59.6008

  100.0000   19.0033   14.03   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)>=eps

    x0=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  =97

Z =

  139.6088  150.3312  153.0517  153.5639

  108.1040  108.61  108.3119  108.1523

   84.1429   67.9096   63.3793   62.4214

   20.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℃,=10℃，所以，

热流量Φ1=λ*[++++]

=1*（100/2+(200-139.6088)+(200-150.3312)+(200-153.0517)+(200-153.5639)/2)

=230.22W

下边界

热流量Φ2=|λ*[+++]-

h*(++++)|

=|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)| = |-4.925|W =4.25W

六、温度等值线

Gauss:

Yacobi:

七、计算小结

导热问题进行有限差分数值计算的基本思想是把在时间、空间上连续的温度场用有限个离散点温度的集合来代替，即有限点代替无限点，通过求解根据傅里叶定律和能量守恒两大法则建立关于控制面内这些节点温度值的代数方程，获得各个离散点上的温度值。

要先划分查分网格，在建立差分代数方程组，用MATLAB或者其他软件编程求解。

高斯-赛德尔迭代法和雅克比迭代法区别在于使用新植和旧值进行下一次迭代，而采用新值迭代的高斯-赛德尔迭代收敛的更快些，但其求解代数方程是不一定得到收敛的解，其原因可能由于迭代方式不合适造成。

在计算热流量过程中，主要是正确利用傅里叶定律和牛顿冷却公式，本题中需要特别注意的一点是后边界是绝热的，因而左右方向上几乎不存在热量的传递，所以看似是二维稳态问题实际上是一维稳态的问题。求解也比较简单。程序运行出来的等温线结果也很好的说明了这一点，温度总体是从上向下递减，热量传递方向是自上而下。    

课后练习——例题4-2（ P173 ）：

“有一各向同性材料的方形物体，其导热系数为常量。已知各边界的温度如图所示，试用高斯-赛德尔迭代求其内部网格节点1、2、3、和4的温度。”

请同学们课后编写计算程序对本题进行迭代计算。编写语言自由。课程结束前提交。作为平时成绩的参考。