上海交通大学硕士研究生课程《传热流动的数值分析大作业》

“传热流动的数值分析”2015年大作业

1. 2维条件下的无粘、不可压缩流体通过出口和入口流过箱体，具体情况如图所示，求该箱体内的流线情况，所有单位为厘米。

（1）流线方程为：22220x y

ψψ

∂∂+=∂∂

使用Gauss-Seidel 线迭代，0.1x y ∆=∆=，误差0.005ξ=，结果输出中，包括在y=0,1,2,3,4,5 处的所有X 处对应的流函数值。

（2）设出口处纵向速度V =0，试采用PSOR 方法，0.1x y ∆=∆=，计算在x=0,1,2,3,4,5 处的所有Y 处对应的流函数值，以及不同的松弛系数和迭代次数的关系曲线（至少三个系数）。 答：（1）该问题为稳态问题，流线方程为椭圆型方程，在求解方程时，首先对方程在计算域内进行离散。计算域为：{}(,)05,05x y x y Ω=≤≤≤≤，在离散时，0.1x y h ∆=∆=∆=，因此可以得到流线方程的差分方程为：

1,,1,,1,,1

22

220i j i j i j i j i j i j h h ψψψψψψ+-+--+-++=∆∆ （1）

整理后可得：

1,1,,1,1

,4

i j i j i j i j i j ψψψψψ+-+-+++=

（2）

在本题中，采用Gauss-Seidel 线迭代方法进行求解，扫描方向选为自左向右，此时有

111,1,11,1,1,4

n n n n i j i j i j i j

n i j

ψψψψψ

++++--+++++=

（3）

由于是线推进，因此在当前线方向求解时，之前的扫描线上的参数已经得到更新，所以

方程可改写为：

11,1,11

1,1,,111444n n

i j i j n n n i j i j i j ψψψψψ+-++++-++-+-=

（4） 其中1

1,n i j ψ+-看着当前迭代层中的已知变量。

从方程（4）中可以明显看出，在扫描线方向上，离散方程组的系数矩阵为三对角矩阵，因此该方程组可以采用追赶法进行求解。

从左向右如此推进，在做完了全区内各列的求解后就完成了一轮迭代，具体的程序代码如附录1所示。执行程序请参见文件包。

对于该方法，进出口边界的流函数假设为ax by ψ=+，对于进口，由于0,00,0.75120, 0ψψ==，因此()1600.75y ψ=-。而对于出口边界，由于5,55,4.250, 120ψψ==，

且进口流量相等，因此可以获得()1605y ψ=-。

计算输出在y=0,1,2,3,4,5 处的所有X 处对应的流函数值，具体数值请参见文件夹PSOR 中执行程序的输出文件RESULTS.TXT ，为了方便观察计算结果，本文采用Tecplot 作出云图，通过云图可以清晰的观察到计算域内的流线，如图1所示。图中每条网格线的交点即为流函数值。

此方法的迭代次数为645次。

图1 Gauss-Seidel 线迭代方法求解结果

（2）本题采用PSOR 方法进行计算。在该方法中，需要对,k i j ψ进行预计算，本文采用Gauss-Seidel 点迭代对该值进行预计算，即：

11

1111/L M k n

n

n k

kl l kl l k kk l l k a a b a ψψψ⨯--==+⎛⎫=++ ⎪⎝⎭∑∑ （5）

在本文中，Gauss-Seidel 点迭代实施步骤方式为先对同一条X 线逐点计算，然后再从左向右进行更新。在程序计算时，方案实施非常简便，在同一轮迭代中，将计算域内所有的流函数值采用同一个矩阵进行存储，当计算更新某个点时，采用矩阵赋值的方法及时进行更新。

PSOR 方法的计算方程为：

()11n n n k k k ψαψαψ-=-+ （6）

和题（1）相似，边界条件也采用相同的，此外，出口纵向速度V=0进一步确定()

ψ=-。

y

1605

在x=0,1,2,3,4,5 处的所有Y处对应的流函数值通过云图进行了展示，如图2所示。此外详细的数值结果也可以在该题的程序文件夹PSOR中的RESULTS.TXT中获得，在此不做详细列出。具体程序代码如附录2所示。

图2 PSOR方法求解结果

此外，本文采用了7个不同的松弛系数进行了计算，分布是1.2,1.4,1.5,1.6,1.8,1.9,2.0,其中当松弛系数为2.0时，计算结果不再收敛，其余计算结果如图3所示。随着松弛系数的增大，迭代次数减小。

图3 不同松弛系数下的迭代次数

2. 试使用一阶波动方程，计算波在一个封闭管道内从t=0s 到t=0.15s 的传播过程。假设声

速V=200m/s ，在t=0s 的初始波形为如下图所示的三角波，请分别用以下方法和步长求解。

10

20

30

40

50

60

70

0510152025

303540u (x )

x

（1） 一阶上风法 （2） Lax-Wendroff （3） BTCS 隐式法 使用以下的不同的步长： （a ） 1.0, 0.005x t ∆=∆= （b ） 1.0, 0.0025x t ∆=∆= （c ） 1.0, 0.00125x t ∆=∆=

请以0.025s 为间隔图示0s 到0.15s 的波传递情况。

答：一阶波动方程为：u u

V t x

∂∂=-∂∂ （7）

按照题目要求，对上述波动方程进行离散可以获得相应的离散方程并进行求解。

（1） 一阶上风法

111n n n

i i i V t V t u u u x x +-∆∆⎛⎫=-+

⎪∆∆⎝⎭ （8） （2） Lax-Wendroff ()()221

111

12222n n n n

n n n i

i

i i i i i V t V t u

u u u u u u x x

++-+-∆∆=--+-+∆∆ （9） （3） BTCS 隐式法 111

1122n n n n i i i i V t V t u u u u x x

+++-+∆∆-

++=∆∆ （10）

从上面的方程中可以看出，方程（8）和（9）两种格式为显示格式，因此可以通过直接逐点求解，而方程（10）为隐式格式，需要对所有点进行联立求解，明显可以看出，该方程的系数矩阵为三对角矩阵，因此可以采用追赶法进行求解。

初值条件：

()0

0 5.0210 5.015.0,0

25015.025.00

25.070.0x x x u x x x x ≤≤⎧⎪-≤≤⎪

=⎨-+≤≤⎪⎪≤≤⎩ ()0,00

(70,0)0

u u ==