计算传热学数值模拟

1、Jacobi 迭代
在Jacobi 迭代法中任一点上未知值的更新是用上一轮迭代中所获得的各邻 点之值来计算的，即
kk k k
l l n l k n k a b T a T /)(1)1()(+=∑≠=- k=1,2,...,L 1×M 1
这里带括号的上角标表示迭代轮数。

所谓一轮是指把求解区域中每一节点之值都更新一次的运算环节。

显然，采用Jacobi 迭代式，迭代前进的方向（又称扫描方向）并不影响迭代收敛速度。

这种迭代法收敛速度很慢，一般较少采用。

但对强烈的非线性问题，如果两个层次的迭代之间未知量的变化过大，容易引起非线性问题迭代的发散。

在规定每一层次计算的迭代轮次数的情况下，有利于Jacobi 迭代有利于非线性问题迭代的收敛。

2、Gauss-Seidel 迭代
在这种迭代法中，每一种计算总是取邻点的最新值来进行。

如果每一轮迭代按T 的下角标由小到大的方式进行，则可表示为：
kk k M L k l n l
kl k
l l n l
kl n k
a b T
a T a T
/)(1
11
)
1(1
1)
()(++
=∑∑⨯+=--≠=
此时迭代计算进行的方向（即扫描方向）会影响到收敛速度，这是与边界条件的影响传入到区域内部的快慢有关的。

3、例题：
一矩形薄板几何尺寸如图所示，薄板左侧的边界温度T L =100K ，右侧温度T R =300K ，上侧温度T T =200K ，下侧温度T B =200K ，其余各面绝热，求板上个节点的温度。

要求节点数目可以变化，写出程序。

解析：
⑴列出描述问题的微分方程和定解条件。

22
220t t x y
∂∂+=∂∂；对于离散化的问题，其微分方程根据热平衡原理得到：
1,1,,1
,1
0i j
i j
i j i j y y x x t t t t x x y y λλλ
λ
-++-∆+∆+∆+∆=⎛⎫
⎛⎫∂∂∂∂⎛⎫⎛⎫ ⎪ ⎪
⎪ ⎪∂∂∂∂⎝⎭
⎝⎭
⎝⎭
⎝⎭
定解条件(边界条件):
T L =100K ，T R =300K ，T T =200K ，T B =200K 。

⑵网格划分示意图：
如下图所示，将薄板划分成m n ⨯(m=n)个网格，求m n ⨯个节点的温度分布。

⑶内部节点的离散化代数方程：
1,,1,,,1
,,1
,0
i j
i j
i j
i j
i j i j
i j i j
y y x y x
x
y
x
t t
t t
t t
t t
λ
λ
λ
λ
-++-----∆+∆+∆+∆=∆∆∆∆即
1,1,,1,1
,40i j i j
i j i j i j
t
t
t
t
t -+-++++-=
边界节点的的离散化代数方程即各节点的温度等于对应边界的温度，不做赘
述。

⑷源程序：
① 采用高斯-赛德尔迭代的程序，如下： m=input('h'); n=input('l'); t=zeros(m,n);
t0=zeros(m,n);
dteps=0.01;
for i=1:m
t(i,1)=200;
t(i,n)=200;
end
for j=1:n
t(1,j)=100;
t(m,j)=300;
end
for k=1:1000
for i=2:m-1
for j=2:n-1
t(i,j)=(t(i-1,j)+t(i+1,j)+t(i,j-1)+t(i,j+1))/4;
end
end
dtmax=0;
for i=2:m-1
for j=2:n-1
dtmax=max(abs(t(i,j)-t0(i,j)),dtmax);
end
end
dtmax
k
t0=t;
contour(t',40);
pause;
if dtmax<dteps break; end
end
②采用雅克比迭代的程序，如下：
m=input('h');
n=input('l');
t=zeros(m,n);
t0=zeros(m,n);
dteps=0.01;
for i=1:m
t(i,1)=200;
t(i,n)=200;
end
for j=1:n
t(1,j)=100;
t(m,j)=300;
end
t0=t;
for k=1:1000
for i=2:m-1
for j=2:n-1
t(i,j)=(t0(i-1,j)+t0(i+1,j)+t0(i,j-1)+t0(i,j+1))/4;
end
end
dtmax=0;
for i=2:m-1
for j=2:n-1
dtmax=max(abs(t(i,j)-t0(i,j)),dtmax);
end
end
dtmax
k
t0=t;
contour(t',40);
pause;
if dtmax<dteps break; end
end
⑸两种方法的收敛速度对比
下面是在相同的条件（m=n=20）下利用高斯-赛德尔迭代和雅克比迭代的得到的最终结果：
高斯-赛德尔迭代（只给出最后部分）
……
dtmax =
0.0099
k =
250
>>
雅克比迭代
……
dtmax =
0.0100
k =
444
>>
由此可以看出，高斯-赛德尔迭代的收敛速度要比雅克比迭代的收敛速度快，因此高斯-赛德尔迭代更加优越。

⑹不同节点数对收敛速度的影响
我们利用高斯-赛德尔迭代法，在m=n=20和m=n=30两种不同的条件下计算节点的温度，结果如下：（只给出m=n=30的结果）
……
dtmax =
0.0099
k =
509
>>
由结果可见迭代后一种情况迭代次数是前一种情况的两倍。

收敛速度明显比前者慢。

画出等温线图如下：（m=n=20的情况下利用高斯-赛德尔迭代的结果）
m=n=30的情况下利用雅克比迭代的结果
计算小结
数值计算是传热学比较重要的研究方法之一。

利用数值计算可以将复杂的解微分方程的问题转化为解代数方程的问题，而解代数方程的问题相对比较简单，完全可以在计算机上实现。

将微分方程转化为代数方程，我们利用网格划分的方法将所研究的物理现象发生的区域离散化，将求所有点参数的问题，转化为求有限节点的问题，这样就可以使问题简单化。

对于上述上述问题我们可以用行立式解代数方程，对于节点数目较少的情况，这种方法比较方便，但节点数目较多时，行立式很难列出来，因此此法就行不通了，迭代法就相对方便的多了。

迭代法包括高斯-赛德尔迭代和雅克比迭代。

前者在计算,i j t 时，1,i j t +、1,i j t -、,1i j t -、,1i j t +的值全部为新值，
即刚刚被迭代得到的值，而后者则利用的是1,i j t +、1,i j t -、,1i j t -、,1i j t +上一次迭代得到的值。

比较而言，同等条件下高斯-赛德尔迭代的收敛速度更快，因此，也根据有优越性，因此我们往往都用这种迭代法进行数值计算分析。

当节点的数目变化时，收敛的速度也随之而变，节点数目越多，收敛的速度越慢，这是显而易见的。

总之，数值计算是传热学非常重要的研究方法，特别是在导热问题的讨论中尤为适用。

研究稳态导热问题，我们常利用高斯-赛德尔迭代法和雅克比迭代法解代数方程，高斯-赛德尔迭代的收敛速度更快，比较常用。