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

“传热流动的数值分析”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 ==
声速V=200.0m/s
在本文中，将三种方法整理在同一个程序中，因此可以在相同程序中实现所有的功能，具体的程序代码如附录3所示。

通过计算，可以获得三种方法，三种不同步长的相关计算结果，一共9张结果，如图4-12所示。

从图中可以看出，对于不同的时间步长选择，三种离散格式都存在假扩散现象，即流动扩散。

对于一阶上风法，固定空间离散大小，减小时间步长会增强流动扩散现象，对于Lax-Wendroff 格式也是如此。

而对于BTCS 隐式法，减小时间步长会削弱流动扩散现象。

u (x )
x
图4 一阶上风法0.005T s ∆=的计算结果
u (x )
x
图5 一阶上风法0.0025T s ∆=的计算结果
u (x )
x
图6 一阶上风法0.00125T s ∆=的计算结果
u (x )
x
图7 Lax-Wendroff 格式0.005T s
∆=的计算结果
u (x )
x
图8 Lax-Wendroff 格式0.0025T s ∆=的计算结果
u (x )
x
图9 Lax-Wendroff 格式0.00125T s ∆=的计算结果
u (x )
x
图10 BTCS 隐式法0.005T s ∆=的计算结果
u (x )
x
图11 BTCS 隐式法0.0025T s ∆=的计算结果
u (x )
x
图12 BTCS 隐式法0.00125T s ∆=的计算结果
附录1 Gauss-Seidel线迭代程序
C--------------------------------------------------------------------------
C 主程序
C--------------------------------------------------------------------------
Program GS
USE MATH
implicit none
real*8 A(1:51,1:51), B(1:51), X(1:51,1:51), Y(1:51)
integer I,J,K,ITER
REAL*8 M(49,4),N(49)
REAL*8 E1,E2,EPS
open(6,file="RESULTS.TXT")
open(7,file="tecplot.DAT")
A=0.0d0
B=0.0d0
X=0.0d0
Y=0.0D0
M=0.0d0
N=0.0D0
E1=0.0D0
E2=1.0D0
EPS=0.005D0
ITER=0
C---------------------------------------------------------
C 边界条件
C---------------------------------------------------------
DO I=1,51
X(I,1)=120.0D0
X(I,51)=0.0D0
ENDDO
DO J=1,51
X(1,J)=0.0D0
X(51,J)=120.0D0
ENDDO
X(1,1)=120.0D0
X(51,51)=0.0D0
DO J=2,8
X(1,J)=0.1D0*(9-J)*160.0D0 !进口边界，此处Y=a*x+(120-0)/0.75*y m/s ENDDO
DO J=50,44,-1
X(51,J)=0.1D0*(51-J)*160.0D0 !出口边界，此处假设a=0.0d0m/s
ENDDO
A=X
DO WHILE(E2>EPS)
E2=0.0D0
DO I=2,50 !沿X方向线推进计算
DO J=2,50
IF(J.EQ.2)THEN
M(J-1,1)=0.0D0
M(J-1,2)=1.0D0
M(J-1,3)=-0.25D0
M(J-1,4)=0.25*(X(I-1,J)+X(I,J-1)+X(I+1,J))
ELSEIF(J.EQ.50)THEN
M(J-1,1)=-0.25D0
M(J-1,2)=1.0D0
M(J-1,3)=0.0D0
M(J-1,4)=0.25*(X(I-1,J)+X(I+1,J)+X(I,J+1))
ELSE
M(J-1,1)=-0.25D0
M(J-1,2)=1.0D0
M(J-1,3)=-0.25D0
M(J-1,4)=0.25*(X(I-1,J)+X(I+1,J))
ENDIF
ENDDO
CALL ZG(49,M,N)
C DO K=1,49
C WRITE(6,"(5F)")M(K,1),M(K,2),M(K,3),M(K,4),N(K) C ENDDO
C PAUSE
C------------------------------------------------------
C 更新X方向的值，这样就可以线推进
C------------------------------------------------------
DO J=2,50
X(I,J)=N(J-1)
ENDDO
ENDDO
DO I=1,51
DO J=1,51
E1=ABS(X(I,J)-A(I,J))
IF(E1>E2)THEN
E2=E1
ELSE
E2=E2
ENDIF
ENDDO
ENDDO
A=X
ITER=ITER+1
WRITE(*,*)E2,ITER
ENDDO
WRITE(6,"(51F)")X
WRITE(7,*)'V ARIABLES="X","Y","F"'
WRITE(7,*)"ZONE I= 51 ,J= 6 ,F=POINT"
DO J=1,51
IF(J.EQ.1 .OR. J.EQ.11 .OR. J.EQ.21 .OR.
1 J.EQ.31 .OR.J.EQ.41 .OR. J.EQ.51) THEN
DO I=1,51
WRITE(7,"(3F9.4)")(I-1)*0.1D0,(J-1)*0.1D0,X(I,J)
ENDDO
ENDIF
ENDDO
CLOSE(7)
CLOSE(6)
END PROGRAM
C------------------------------------------------------------------------------
C 追赶法
C------------------------------------------------------------------------------ MODULE MATH
CONTAINS
SUBROUTINE ZG(N1,M,X)
implicit real*8(a-h,o-z)
INTEGER::N1
REAL*8 MXLX(N1,4),X(N1),M(N1,4)
MXLX=M
DO i=1,N1
IF(i==1)then
MXLX(i,2)=MXLX(i,2)
MXLX(i,4)=MXLX(i,4)/MXLX(i,2)
ELSE
MXLX(i,2)=MXLX(i,2)-MXLX(i,1)*MXLX(i-1,3)
MXLX(i,4)=(MXLX(i,4)-MXLX(i,1)*MXLX(i-1,4))/MXLX(i,2) ENDIF
MXLX(i,3)=MXLX(i,3)/MXLX(i,2)
ENDDO
DO i=N1,1,-1
IF(i==N1)then
X(i)=MXLX(i,4)
ELSE
X(i)=MXLX(i,4)-MXLX(i,3)*X(i+1) ENDIF
ENDDO
END SUBROUTINE
END MODULE
附录2 PSOR方法程序
C---------------------------------------------------------------------
C 主程序
C---------------------------------------------------------------------
Program PSOR
implicit none
real*8 A(1:51,1:51), B(1:51), X(1:51,1:51), Y(1:51),X1(1:51,1:51)
integer I,J,K,ITER
REAL*8 M(49,4),N(49)
REAL*8 E1,E2,EPS
REAL*8 AF,BF
open(6,file="RESULTS.TXT")
open(7,file="tecplot.DAT")
A=0.0d0
B=0.0d0
X=0.0d0
Y=0.0D0
M=0.0d0
N=0.0D0
E1=0.0D0
E2=1.0D0
EPS=0.005D0
ITER=0
AF=1.9
C---------------------------------------------------------
C 边界条件
C---------------------------------------------------------
DO I=1,51
X(I,1)=120.0D0
X(I,51)=0.0D0
ENDDO
DO J=1,51
X(1,J)=0.0D0
X(51,J)=120.0D0
ENDDO
X(1,1)=120.0D0
X(51,51)=0.0D0
DO J=2,8
X(1,J)=0.1D0*(9-J)*160.0D0 !进口边界，此处Y=a*x+(120-0)/0.75*y m/s ENDDO
DO J=50,44,-1
X(51,J)=0.1D0*(51-J)*160.0D0 !出口边界，由于V=0，所以a=0.0
ENDDO
A=X
X1=X
DO WHILE(E2>EPS)
E2=0.0D0
DO I=2,50 !采用高斯赛德尔迭代高斯赛德尔迭代作为预计算值
DO J=2,50
BF=X(I,J)
X(I,J)=0.25*(A(I,J+1)+A(I,J-1)+A(I-1,J)+A(I+1,J)) !此处用A代替X，这样可以通过在后续的设定A=X体现已考虑及时更新X的计算值
X(I,J)=(1-AF)*BF+AF*X(I,J) !AF为松弛因子
A=X !及时更新所计算的值，体现为高斯赛德尔迭代
ENDDO
ENDDO
DO I=1,51
DO J=1,51
E1=ABS(X(I,J)-X1(I,J))
IF(E1>E2)THEN
E2=E1
ELSE
E2=E2
ENDIF
ENDDO
ENDDO
X1=X
A=X
ITER=ITER+1
WRITE(*,*)E2,ITER
ENDDO
WRITE(6,"(51F)")X
WRITE(7,*)'V ARIABLES="X","Y","F"'
WRITE(7,*)"ZONE I= 6 ,J= 51 ,F=POINT"
DO J=1,51
DO I=1,51
IF(I.EQ.1 .OR. I.EQ.11 .OR. I.EQ.21 .OR.
1 I.EQ.31 .OR.I.EQ.41 .OR. I.EQ.51) THEN
WRITE(7,"(3F9.4)")(I-1)*0.1D0,(J-1)*0.1D0,X(I,J)
ENDIF
ENDDO
ENDDO
CLOSE(7)
CLOSE(6)
END PROGRAM
附录3 波动方程求解程序
C----------------------------------------------------------------
C 主程序
C----------------------------------------------------------------
PROGRAM MAIN
USE MATH
IMPLICIT NONE
REAL*8 DT,DX,V !时间步长，空间步长，声速
REAL*8 NP !方法选择，1为一阶上风法，2为Lax-Wendroff，3为BTCS隐式法
REAL*8 UPRE(0:70),UNEX(0:70) !当前时层和下一时层的速度值
INTEGER I,J
REAL*8 TIME
REAL*8 X
REAL*8 TI
REAL*8 A(69,4),Y(69)
TIME=0.0D0
C---------------------------------------------------
C 初值确定
C---------------------------------------------------
DO I=0,70
X=DBLE(I)
IF(X.LE.5.0D0)THEN
UPRE(I)=0.0D0
ELSEIF(X.GE.5.0D0 .AND. X.LE.15.0D0)THEN
UPRE(I)=2.0D0*X-10.0D0
ELSEIF(X.GE.15.0D0 .AND. X.LE.25.0D0)THEN
UPRE(I)=-2.0D0*X+50.0D0
ELSE
UPRE(I)=0.0D0
ENDIF
ENDDO
UNEX=UPRE
V=200.0D0
c DT=0.005
c DT=0.0025
DT=0.00125
DX=1.0
NP=3
OPEN(7,FILE="INITIAL.DAT")
OPEN(8,FILE="RESULT.DAT")
DO I=0,70
WRITE(7,"(1F8.5,4X,1F8.5)")DBLE(I),UPRE(I)
ENDDO
CLOSE(7)
IF(NP.EQ.1)THEN
WRITE(*,*)"数值格式为一阶上风法"
WRITE(*,*)"DT=",DT,"DX=",DX
ELSEIF(NP.EQ.2)THEN
WRITE(*,*)"数值格式为Lax-Wendroff"
WRITE(*,*)"DT=",DT,"DX=",DX
ELSEIF(NP.EQ.3)THEN
WRITE(*,*)"数值格式为BTCS隐式法"
WRITE(*,*)"DT=",DT,"DX=",DX
ENDIF
C-----------------------------------------------------
TIME=TIME+DT !此处设置可以避免多算一步
DO WHILE(TIME.LE.0.15D0)
IF(NP.EQ.1)THEN !一阶上风法
DO I=1,69
UNEX(I)=(1.0D0-V*DT/DX)*UPRE(I)+V*DT/DX*UPRE(I-1) ENDDO
ELSEIF(NP.EQ.2)THEN
DO I=1,69
UNEX(I)=UPRE(I)-V/2.0D0*DT/DX*(UPRE(I+1)-UPRE(I-1))+
2 (V*DT/DX)**2/2*(UPRE(I+1)-2.0D0*UPRE(I)+UPRE(I-1))
ENDDO
ELSEIF(NP.EQ.3)THEN
DO I=1,69
IF(I.EQ.1)THEN
A(I,1)=0.0D0
A(I,2)=1.0D0
A(I,3)=V/2.0D0*DT/DX
A(I,4)=V/2.0D0*DT/DX*UPRE(I-1)+UPRE(I)
ELSEIF(I.EQ.69)THEN
A(I,1)=-V/2.0D0*DT/DX
A(I,2)=1.0D0
A(I,3)=0.0D0
A(I,4)=-V/2.0D0*DT/DX*UPRE(I+1)+UPRE(I)
ELSE
A(I,1)=-V/2.0D0*DT/DX
A(I,2)=1.0D0
A(I,3)=V/2.0D0*DT/DX
A(I,4)=UPRE(I)
ENDIF
ENDDO
CALL ZG(69,A,Y)
DO I=1,69
UNEX(I)=Y(I)
ENDDO
ENDIF
UPRE=UNEX
TI=ABS(CEILING(TIME/0.025D0)-TIME/0.025D0)
IF(TI.LE.1.0D-5)THEN
DO I=0,70
WRITE(8,"(1F8.5,4X,1F8.5,4X,1F8.5)")TIME, DBLE(I),UNEX(I)
ENDDO
ENDIF
TIME=TIME+DT
ENDDO
CLOSE(8)
ENDPROGRAM
其中，子程序ZG()见附录1。