计算传热学数值模拟

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 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;

相关文档
最新文档