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