计算流体力学与传热学大作业

########学院

计算流体力学与传热学

学号：

专业：

学生姓名：

任课教师：教授

2013年12月

目录

第一章验证显式格式的稳定性 (4)

1.1 概述 (4)

1.2 数学推导 (4)

1.3 问题描述 (4)

1.4 数值模拟 (4)

1.5 结果及分析 (5)

第二章判断肋片可以按一维问题处理的主要依据 (6)

2.1 概述 (6)

2.2 问题描述及算法 (6)

2.3 数值模拟 (7)

2.4 结果及分析 (8)

第三章三层墙导热 (9)

3.1 概述 (9)

3.2 问题描述 (9)

3.3 TDMA算法 (9)

3.4 结果 (10)

第四章一维无源稳态对流扩散问题 (11)

4.1 公式及初值 (11)

4.2 情况一 (11)

4.3 情况二 (12)

4.4 情况三 (13)

第五章用ADI算法计算长方肋内的温度分布 (14)

5.1 问题描述 (14)

5.2 初始参数 (14)

5.3 情况一，一列列扫 (14)

5.4 情况二，一行行扫 (14)

5.5 情况三，采用ADI算法 (15)

5.6 结果分析 (15)

参考文献 (16)

第一章 验证显式格式的稳定性

1.1 概述

将一维非稳态热传导方程用显式格式差分化为代数方程，在求解的迭代过程中必须满足一定的条件，才能使方程收敛且结果正确。此处即验证β≤?。

1.2 数学推导

方程： 22T t T x α??=??

（1）

显式离散格式： 此处时间向前差分，空间中心差分

111

22n n n n n

i i i i i T T T T T t x α+-+--+=??

1112(2)n n n n n

i i i i i t T T T T T x

α

+-+?-=-+? 令β=2

t

x α

??则： 111(2)n n n n n

i i i i i T T T T T β+-+-=-+ （2）

误差也应该满足上式，故：

()()1()()()2()()i

i

i

i

i

Ikx Ikx Ik x x Ikx Ik x x n n n n n

T e T e T e T e T e ψψβψψψ----?--+?+??-=-+??

()()()1()12()()()i

i

i

i

Ikx Ikx Ik x x Ik x x n n n n

T e T e T e T e ψβψβψψ----?-+?+??=-++??

()()1()12()()i

i

i

Ikx Ikx Ikx

n n Ik x Ik x n T e T e e e T e ψβψβψ---+-??=-++

()()1()

121()

n Ik x Ik x n

T e e T ψββψ+-??=-++≤ 因此 β≤?。即当β≤? 时方程（2）才会有收敛的解。

1.3 问题描述

在验证过程中同时可模拟一个实际问题，即冬季里墙壁中的温度分布。此时室内壁温设为Tl=30.0℃，室外壁温Tr=-25.0℃，墙壁以11号楼为例，L=1m ，热扩散系数ɑ=alfa=1.33e-6m 2/s

然后分别取β=0.4，n=10和β=0.6，n=10两种情况，看最后的结果是否收敛和正确。

1.4 数值模拟

程序流程如下：

1.5 结果及分析

㈠当beta=0.4、n=10时，运行程序可收敛，得结果如下：

30.00 24.47 18.94 13.42 7.91 2.41 -3.09 -8.58 -14.06 -19.53 -25.00 30.00 24.47 18.95 13.43 7.91 2.41 -3.09 -8.57 -14.05 -19.53 -25.00 30.00 24.47 18.95 13.43 7.92 2.41 -3.08 -8.57 -14.05 -19.53 -25.00 30.00 24.47 18.95 13.43 7.92 2.42 -3.08 -8.57 -14.05 -19.53 -25.00 30.00 24.48 18.95 13.44 7.92 2.42 -3.08 -8.56 -14.05 -19.52 -25.00

可看出最后已经稳定。 用图表显示可得：

可看出最后温度在墙壁内基本成线性分布。

输入ɑ、β、L 、n 值 设定初值：T(0)=Tl T(i)=Tr(i=1，n)

迭代计算：

TN(i)=T(i)+β×(T(i+1)+T(i-1)-2×T(i))

求解输出 判断是否收结束

是 否

T(i)=TN(i)

开始

㈡当beta=0.6、n=10时，运行程序不收敛，得到如下结果：

30.00 -184.31 395.52 -498.83 562.01 -562.79 483.59 -411.54 253.44 -160.76 -25.00 30.00 292.17 -488.99 674.28 -749.37 739.92 -681.32 524.52 -394.07 169.21 -25.00 30.00 -333.83 677.67 -877.87 998.39 ******* 894.93 -750.14 495.06 -285.28 -25.00 30.00 491.37 -862.55 1181.21 ******* 1337.27 ******* 984.02 -720.26 339.09 -25.00 30.00 -597.80 1176.06 ******* 1777.14 ******* 1639.36 ******* 937.92 -514.98 -25.00 30.00 843.19 ******* 2082.30 ******* 2410.96 ******* 1820.11 ******* 650.75 -25.00 30.00 ******* 2060.31 ******* 3169.92 ******* 2985.10 ******* 1746.07 -935.83 -25.00 30.00 1467.32 ******* 3688.80 ******* 4341.70 ******* 3337.52 ******* 1219.81 -25.00 综上所述，只有当beta≦?时，结果才会稳定收敛。

第二章判断肋片可以按一维问题处理的主要依据2.1 概述

在实际问题中，经常将肋片按一维处理，即假设表面上的换热热阻1/h远远大于肋片中的导热热阻δ/λ（即Bi很小），因而在任一截面上肋片温度可认为是均匀的。这样假设的确便于数学处理，但是否合理，下面就进行验证。

2.2 问题描述及算法

粗而短的肋片的分析

节点划分

计算在下表所示的两种条件下肋片的效率，并与一维分析解的结果比较。

程序设计中各个点的温度为：

1 T(i,j)=0.25*(T(i-1,1)+T(i+1,1)+2*T(i,2))

2 T(i,j)=0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1)) 3a T(i,j)=(T(i-1,N)+T(i+1,N)+2*T(i,N-1))/(4+2*Bi) 3b T(i,j)=(T(M,j-1)+T(M,j+1)+2*T(M-1,j))/(4+2*Bi) 4a T(i,j)=(T(M-1,N)+T(M,N-1))/(2+2*Bi)

4b T(i,j)=(T(M-1,1)+T(M,2))/(2+Bi)

2.3 数值模拟

程序流程图为：

工况 T0/℃ Tf/℃ h/W/(m 2.℃) λ

/W/(m.K) δ/m H/m 1 100 20 50 100 0.02 0.04 2

100

20

400

8

0.02

0.08

以计算得到的T(i,j)计算临近节点的温度

否

是

结束 判断是否稳定：

开始 输入计算参数 h ，λ，δ,L ，M ，N 求Bi

设定温度初值: J=1,N T(1,j)=T0-Tf 其余的节点为T(i,j)=0.0 令c(i,j)=T(i,j) 根据初值计算节点温度(共6类) 1 T(i,j)=0.25*(T(i-1,1)+T(i+1,1)+2*T(i,2)) 2 T(i,j)=0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1)) 3a T(i,j)=(T(i-1,N)+T(i+1,N)+2*T(i,N-1))/(4+2*Bi) 3b T(i,j)=(T(M,j-1)+T(M,j+1)+2*T(M-1,j))/(4+2*Bi) 4a T(i,j)=(T(M-1,N)+T(M,N-1))/(2+2*Bi) 4b T(i,j)=(T(M-1,1)+T(M,2))/(2+Bi)

2.4 结果及分析

肋效率的数值计算结果

工况1（Bi=0.01）

工况2 （Bi=1）

工况 节点 Bi 二维数值计算的η值 一维数值计算的η值 相对偏差 1 9×5 0.01 0.970 0.971 0.1％ 2

17×5

1

0.185

0.200

8.1％

当Bi=0.01时，得出的肋效率与一维公式计算结果的差别完全可以忽略，而对于Bi=1的长肋片，它们之间差别较大。由此可见，判断肋片中导热可否按一维问题处理的标准应该是Bi 尽量小。

第三章 三层墙导热

3.1 概述

因为各层材料的导温系数不同，在分界面处的导温系数需用调和平均数计算。用有限体积法计算时还要注意边界条件的处理。

3.2 问题描述

一居民的砖墙结构如图所示：

各层参数如下：

3.3 TDMA 算法

一维稳态不含源项的控制方程为

0)(=dx dT k dx d

用有限体积法可化为：

0)()(=---w x TW

TP kw e x TP TE ke

δδ

参数 k （W/㎡.k ） L/㎜ N 塑料板 0.16 12 10 杨木 0.141 25 20 砖块

0.72

200

100

室内空气 TL=（25±5）℃ hL=7W/(㎡.k)

室外空气 Tr=5℃

hr=28W/(㎡.k)

塑料板

杨木

一般建筑用砖块

设为：a p *TP=a e *TE+a w *TW

其中 a e =

e x ke

)(δ ；a

w=

w x kw )(δ

a p=a e +a w

为使方程具一般性，可设为A i *T i =B i *T i-1+C i *T i+1+D i 对于此问题，墙壁两端均为对流换热：

当i=1时，A1=hL+ke/xe ；B1=0 ；C1=ke/xe ；D1=hL*TL ；P1=C1/A1；Q1=D1/A1 当i=2--N-1时，Ai=Bi+Ci ；Bi=kw/xw ；Ci=ke/xe ；Di=0

1*--=

i i i i P B A C Pi 11

**---+=i i i i i i

i P B A Q B D Q

当i=N 时，Ai=hr+kw/xw ；Bi=kw/xw ；Ci=0；Di=hr*Tr ； Pi=0；

11**---+=

i i i i i i i P B A Q B D Q

由TDMA 算法知：Ti=Qi

然后由公式 T i-1=P i-1*T i +Q i-1 反向回带可依次求出其余温度值。

3.4 结果

当温度TL 在（25±5℃）范围内变动时，温度在三层墙壁内的分布曲线图如下：

第四章 一维无源对流扩散问题

4.1 公式及初值 公式： 其中ρ=1,Φ0=1,ΦL =0,L=1,Γ=1

TN 为解析解，T 为数值解，采用混合格式hybrid 编程计算。

4.2 情况一：

u=0.1 ，δx=0.2 ，Pe=0.2

数据：

x TN T 0.00 1.0000 1.0000 0.20 0.8711 0.8714 0.40 0.7138 0.7141 0.60 0.5215 0.5220 0.80 0.2868 0.2871

1.00 0.0000 0.0000

图表：

)

()(dx d dx d u dx d φφρΓ=

Pe<2时，程序中采用中心差分格式，结果较理想。

4.3 情况二：

u=2.5 ，δx=0.1 ，Pe=2.5 数据：

x TN T

0.00 1.0000 1.0000

0.10 1.0000 1.0000

0.20 1.0000 1.0000

0.30 1.0000 1.0000

0.40 1.0000 1.0000

0.50 1.0000 1.0000

0.60 1.0000 1.0000

0.70 0.9994 1.0000

0.80 0.9933 1.0000

0.90 0.9179 1.0000

1.00 0.0000 0.0000

图表：

Pe>2时，采用迎风格式，结果满足物理规律。

4.4 情况三：

u=2 ，δx=0.05 ，Pe=1

数据：

x TN T

0.00 1.0000 1.0000 0.55 0.5733 0.5733 0.05 0.9702 0.9702 0.60 0.5215 0.5216 0.10 0.9388 0.9388 0.70 0.4100 0.4100 0.15 0.9058 0.9058 0.75 0.3499 0.3500 0.20 0.8711 0.8712 0.80 0.2868 0.2868 0.25 0.8347 0.8347 0.85 0.2204 0.2204 0.30 0.7964 0.7964 0.90 0.1505 0.1506 0.35 0.7561 0.7561 0.95 0.0772 0.0772 0.40 0.7138 0.7138 1.00 0.0000 0.0000 0.45 0.6693 0.6693

0.50 0.6225 0.6225

图表：

Pe<2，采用中心差分，结果较理想。

第五章

用ADI 算法计算长方肋内的温度分布

5.1 问题描述

运用三种算法计算长方肋内的温度分布，根据结果好坏和计算时间长短判断各算法的优劣。

5.2 初始参数

5.3 情况一，一列列扫

只一列列的扫，得温度分布如下：

5.4 情况二，一行行扫

只一行行的扫，得温度分布如下：

参数 h/W/(m 2.℃

) k/W/(m.K) δ/m

L/m TL/℃ Tf/℃ M(肋高方向节点数) N(肋厚方向

节点数) 数据 400.0 8.0 0.04

0.12 100.0

20.0

25 20

先一列列扫，再一行行扫，即采用ADI算法，得温度分布如下：

5.6 结果分析

可看出得到的温度分布基本相同，即三种方法均可以得出结果，但计算速度

却差别很大，此例中花费时间由小到大依次为：情况三<情况一<情况二。

参考文献

[1] Professor Wang Ge ,Computational Fluid Dynamics and Heat

Transfer ,Harbin Engineering University,2013

[2] 杨世铭、陶文铨编著，传热学（第四版），高等教育出版社，2006

[3] 江春波、张永良、丁则平编著，计算流体力学，中国电力出版社，2007

[4] 彭国伦编著，Fortran 95程序设计，中国电力出版社，2002

[5] 陶文铨编著，数值传热学（第二版），西安交通大学出版社，2001