计算热物理

目录

题目 (2)

1.方程和定解条件的无量纲化 (4)

2.控制方程的离散 (4)

3.稳定性分析 (5)

4.温度场的计算求解 (6)

5.温度场的等值线图 (9)

6.分析比较 (12)

7.个人心得体会 (14)

附录 (16)

题目：定义在正方形区域（如图所示）H*H=3m*3m 的二维非稳态导热问题的控制方程及其定解条件为

其中ρ，c，λ分别为导热体的密度、比热和导热系数，且ρ=7820kg/m3，c=460J/(kg•K)，λ=15W/(m•K)。将定义域各方向均匀分成10个区块，取点中心网格分割方式，各边界共11个节点; 或者取块中心网格分割，各边界共10个节点，（两种网格，至少做一种），按以下要求，数值求解该问题。

（1）取（X，Y）=（x，y）

H ，θ=T−T0

T1−T0

，τ=t

ρcH2/λ

，将方程和定解条

件无量纲化，并在无量纲化条件下求解。

（2）采用以下方法对控制方程离散：

（a）用控制容积法离散成显式、全隐式格式；

（b）用ADI格式离散。

（3）用Fourier分析方法对以上3种格式的稳定性进行分析。（4）嵌入初、边条件，分别对以上3种格式构成的定解问题编制计算程序，作数值求解，直到达到稳态温度分布——能直接求解显式格式和ADI格式构成的代数方程）直接求解，不能直接求解的（隐式格式构成的代数方程）采用以下方法迭代求解：

（a）Jacobi简单点迭代；

（b）Guess-Seidel点迭代；

（c）带松弛的Guess-Seidel点迭代；

（d）带松弛的Guess-Seidel线迭代；

（e）基于Guess-Seidel的交替方向线迭代。

并对收敛结果做出比较。

（5）选择全隐式格式在5个不同时刻（包括1个到达稳态的时刻）的计算结果，将无量纲计算量还原成有量纲量，分别画出温度T在平面空间的等值线图。从5个不同时刻的等值线图，分析时间发展中T的变化特性。

（6）根据数值计算情况，（精度高低，CPU时间消耗，编制程序难易等）对各种不同的离散方法及代数解法进行分析比较，并做出结论。

（7）综合上述内容和作业情况，写出作业报告，并对此类作业提出

建议和要求。

解：

（1）方程和定解条件的无量纲化：

将（X ，Y ）=（x ，y ）

H

，θ=

T−T 0T 1−T 0

，τ=

t

ρcH 2/λ

代入原方程组，无

量纲化，得到：

（2） 控制方程的离散： A ． 控制容积法

假定非稳态项中温度随空间为阶梯分布，扩散项中温度随时间线性分布、随空间阶梯分布 显式格式:

(θP −

θP 0)∆X∆Y

=(θE 0−θP

(

δX )e

−

θP 0−θW

(δX )w

)∆Y∆τ+(θN 0−θP 0

(δY )n

−

θP 0−θS

(δY )s

)∆X∆τ

写成：a P θP =a E θE 0+a W θW 0+a N θN 0+a S θS 0

+b 的形式，

则：a E =∆Y

(δX )e ，a W =∆Y

(

δX )w

，a N =∆X

(

δY )n

，a S =∆X

(

δY )s

a P =

∆X∆Y

∆τ

，a P 0=

∆X∆Y ∆τ

−a E −a W −a N −a S ，b =a P 0θP 0

全隐式格式：

(θP −θP 0)∆X∆Y =(θE

−θP (

δX )e

−

θP −θW (δX )w

)∆Y∆τ+(θN

−θ

P

(δY )

n

−θP −θS (δY )s

)∆X∆τ

写成：a P θP =a E θE +a W θW +a N θN +a S θS +b 的形式，

则：a E=∆Y

(δX)e ，a W=∆Y

(δX)w

，a N=∆X

(δY)n

，a S=∆X

(δY)s

a P=∆X∆Y

∆τ+a E+a W+a N+a S，a P0=∆X∆Y

∆τ

，b=a P0θP0

B．ADI格式：

θi,j∗−θi,j n ∆τ/2=

θi+1,j

∗−2θ

i,j

∗+θ

i−1,j

∗

∆X2

+

θi,j+1

n−2θ

i,j

n+θ

i,j−1

n

∆Y2

θi,j n+1−θi,j∗∆τ/2=

θi+1,j

∗−2θ

i,j

∗+θ

i−1,j

∗

∆X2

+

θi,j+1

n+1−2θ

i,j

n+1+θ

i,j−1

n+1

∆Y2

（3）稳定性分析:

A．显示格式：

a PθP=a EθE0+a WθW0+a NθN0+a SθS0+b

a P g=a E e ik X∆X+a W e−ik X∆X+a N e ik Y∆Y+a S e−ik Y∆Y+a P0(*)

本题中定义域各方向均匀分成10个区块，则

(δX)e=(δX)w=(δY)n=(δY)s=∆X=∆Y=0.1

a E=a W=a N=a S=1，a P=∆X∆Y

∆τ

则(*)式可化为：g=2cos(k X∆X)+2cos(k Y∆Y)+(∆X∆Y/∆τ−4)

∆X∆Y/∆τ

假设k X=k Y，则g=1−8∆τ

∆X2sin2(k∆X

2

)

有|g|≤1恒成立，则∆τ

∆X2sin2(k∆X

2

)≤1

4

∆τ∆X2≤1

4

，∆τ≤1

4

∆X2。

即∆τ≤1

4

∆X2时，该显式格式收敛。

本题中，取∆X=0.1，故∆τ≤0.0025。

B．全隐式格式：

a PθP=a EθE+a WθW+a NθN+a SθS+b

a P g=g(a E e ik X∆X+a W e−ik X∆X+a N e ik Y∆Y+a S e−ik Y∆Y)+a P0(**)