中科大计算热物理大作业4.3

E − Boundary
bi = 4
p = −ci
i
ai pi−1 +bi
q = di−aiqi−1
i
ai pi−1 +bi
θi = piθi+1 + qi θN = qN
N = 11
取收敛判据为
| θP(n+1)−θP(n)
θP(n)
|max
≤
10−4，编写程序计算得到稳态温度分布为
(17) (18) (19)
图 2: ADI 离散格式下稳态无量纲化温度分布
4.3 控制容积法隐式格式
同理，按照显式计算思路，可以写出控制容积法隐式格式边界条件
aBθB = aEθE + aSθS + aN θN + b
W − Boundary :
aE
=
∆Y ∆X
aS
=
aN
=
∆X 2∆Y
a0P
=
∆X ∆Y 2∆τ
(1)
2 控制方程离散
2.1 控制容积法
下面利用控制容积法将方程 (1) 离散成显式格式。选择控制体 P，在时间间隔 [τ, τ + ∆τ ] 内将方程对时间
和空间积分，有
∫ n ∫ e ∫ τ +∆τ
s wτ
∂θ ∂τ
dτ dX dY
∫ τ +∆τ ∫ n ∫ e
=
τ
sw
∂2θ ∂X2
dX
dY
dτ


gS1/−2N
=
∆X∆Y /∆τ −4 sin2(kX ∆X/2) ∆X∆Y /∆τ +4 sin2(kY ∆Y /2)

gW 1/2−E
=
∆X∆Y /∆τ −4 sin2(kY ∆Y /2) ∆X∆Y /∆τ +4 sin2(kX ∆X/2)
(13)
要保证 ADI 交替前进一次完整时间步中 |g| ≤ 1，即 |gS1/−2N · gW 1/2−E| ≤ 1，计算可知该不等式恒成立，所以 ADI 格式为恒稳定格式。
∆X aS = aN = ∆Y
a0P
=
∆X ∆Y ∆τ
− aW
− aE − aS − aN
(3)
aP = aW + aE + aS + aN + a0P b = a0P θP0
显见，全隐式格式为
aP θP = aW θW + aEθE + aSθS + aN θN + b
∆Y aW = aE = ∆X
其中 ρ、c、λ 分别为导热体的密度、比热、和导热系数，且 ρ = 7820 kg/m3，c = 460 J/(kg · K)，λ = 15 W /(m · K)。将定义域各方向均匀分成 10 个区块，取点中心网格分割方式，各边界 11 个节点；或者取块中 心网格分割，各边界共十个节点。按要求数值求解该问题。
条件的方程写为
θP = 0.1(θW 0 + θE0 + θS0 + θN0 ) + 0.6θP0

θP = 0.2θE0 + 0.1(θS0 + θN0 ) + 0.6θP0 + 0.02 W − Boundary θP = 0.2θW 0 + 0.1(θS0 + θN0 ) + 0.6θP0 + 0.02 E − Boundray
∫ τ +∆τ ∫ e ∫ n
+
τ
ws
∂2θ ∂X2dY源自dX dτ(2)假定非稳态项中温度随时间为阶梯分布，扩散项中温度随空间为分段线性分布，随时间为阶梯分布，取显式格
式，积分可以得到
aP θP = aW θW 0 + aE θE0 + aS θS0 + aN θN0 + b
∆Y aW = aE = ∆X



W − E direction


−θS(n+1/2) + 4θP(n+1/2) − θN(n+1/2) = 2θEn + 0.2 −θS(n+1/2) + 4θP(n+1/2) − θN(n+1/2) = θW n + θEn −θS(n+1/2) + 4θP(n+1/2) − θN(n+1/2) = 2θW n + 0.2
(20)
S − Boundary : θ(X, 0, τ ) = 1
N − Boundary : θ(X, 1, τ ) = 0
4
4.3.1 Jaboci 简单点迭代
由上可知，控制容积法隐式格式具有恒稳定性。此处取 ∆τ = 0.01，则隐式格式的 Jaboci 点迭代形式可写
为

θP(n+1) = 0.2(θW n + θEn + θSn + θNn + θpn)
(9)
由于题中所选取 ∆X 和 ∆Y 相等，则可以代入系数得到放大因子：
g
=
1
−
4∆τ ∆X ∆Y
(sin2(kX ∆X/2)
+
sin2(kY
∆Y
/2))
(10)
要保证 |g| ≤ 1，可得
∆X ∆Y ∆τ ≤
(11)
4
则控制容积法显式格式为有条件稳定格式。
同理，可以得到控制容积法隐式格式的方法因子
(15)
取收敛判据为
| θP(n+1)−θP(n)
θP(n)
|max
≤
10−4，编写程序计算得到稳态温度分布为
图 1: 控制容积法显式格式下稳态无量纲化温度分布
4.2 ADI 格式
ADI 离散格式的方程为式 (7)，下面采用点中心网格，引入点中心网格的边界条件，取 ∆τ = 0.01，可得
到如下形式
∆X∆Y /∆τ g = ∆X∆Y /∆τ + 4(sin2(kX ∆X/2) + sin2(kY ∆Y /2))
(12)
可见，|g| ≤ 1 恒成立，该格式恒稳定。
下面对方程的 ADI 离散格式进行分析。将 (8) 分别代入到 (7) 的两个方程中，进一步代入系数，可以得到
ADI 格式两步的放大因子
根据 von Neumann 稳定性分析方法，将 (1) 中方程的解写为 Fourier 级数形式，而不同级数下的解满足叠
加原理，故可取一项进行分析，取
θ(X, Y, τ ) = gnei(kX X+kY Y )
(8)
将上式代入到 (3) 中，可以得到
aP gn+1 = aW gne−ikX ∆X + aE gneikX ∆X + aS gne−ikY ∆Y + aN gneikY ∆Y + a0P
aN bN
θN
dN
W −E :




a1 = 0 ai =−1
ci = −1 cN = 0
bi = 4

W − Boundary
di = 2θEn + 0.2 dN = 2θEn + 1.2
M iddle
di = θW n + θEn dN = θW n + θEn + 1
1 无量纲化处理
将题中所给无量纲化参数直接带入方程和边界条件中，即可得到：


∂θ ∂τ
=
(
∂2θ ∂X2
+
∂2θ ∂Y 2
)
qw
=
−
∂θ ∂X
|X
=0
=
qe
=
∂θ ∂X
|X =1
=
1

θ(X, 0, τ ) = θ1 = 1 θ(X, 1, τ ) = θ0 = 0 θ(X, Y, 0) = 1 − Y
=
θin++11,j/2 − 2θin,j+1/2 + θin−+11,/j2 ∆X 2
+
θin,j++11 − 2θin,j+1 + θin,j+−11 ∆Y 2
(6)
可以整理为


−aS θSn+1/2 + aP θPn+1/2 − aN θNn+1/2 = aW θW n + aE θEn + b −aW θW n+1 + aP θPn+1 − aE θEn+1 = aS θSn+1/2 + aN θNn+1/2 + b
∆Y ∆X
aS
=
aN
=
∆X 2∆Y
a0P
=
∆X ∆Y 2∆τ
− aW
− aS − aN
aB = aW + aS + aN + a0P b = a0P θP0 + qe∆Y
(14)
S − Boundary : θ(X, 0, τ ) = 1
N − Boundary : θ(X, 1, τ ) = 0
根据前面的分析，显式格式的稳定性条件 (11) 为 ∆τ ≤ 0.0025，在编程计算中取 ∆ = 0.001，可以将引入边界
W − Boundary E − Boundary

S − N direction

θ(X, 0, τ ) = 1 S − Boundary −θW (n+1) + 4θP(n+1) − θE(n+1) = θW (n+1/2) + θE(n+1/2) θ(X, 1, τ ) = 0 N − Boundary