“传热流动的数值分析”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,,122220i j i j i j i j i j i j h h ψψψψψψ+-+--+-++=∆∆ （1）整理后可得：1,1,,1,1,4i j i j i j i j i j ψψψψψ+-+-+++=（2）在本题中，采用Gauss-Seidel 线迭代方法进行求解，扫描方向选为自左向右，此时有111,1,11,1,1,4n n n n i j i j i j i jn i jψψψψψ++++--+++++=（3）由于是线推进，因此在当前线方向求解时，之前的扫描线上的参数已经得到更新，所以方程可改写为：11,1,111,1,,111444n ni j i j n n n i j i j i j ψψψψψ+-++++-++-+-=（4） 其中11,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 点迭代对该值进行预计算，即：111111/L M k nnn kkl 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进一步确定()ψ=-。

y1605在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 的初始波形为如下图所示的三角波，请分别用以下方法和步长求解。

102030405060700510152025303540u (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 uV t x∂∂=-∂∂ （7）按照题目要求，对上述波动方程进行离散可以获得相应的离散方程并进行求解。

（1） 一阶上风法111n n ni i i V t V t u u u x x +-∆∆⎛⎫=-+⎪∆∆⎝⎭ （8） （2） Lax-Wendroff ()()22111112222n n n nn n n iii i i i i V t V t uu u u u u u x x++-+-∆∆=--+-+∆∆ （9） （3） BTCS 隐式法 1111122n n n n i i i i V t V t u u u u x x+++-+∆∆-++=∆∆ （10）从上面的方程中可以看出，方程（8）和（9）两种格式为显示格式，因此可以通过直接逐点求解，而方程（10）为隐式格式，需要对所有点进行联立求解，明显可以看出，该方程的系数矩阵为三对角矩阵，因此可以采用追赶法进行求解。

初值条件：()00 5.0210 5.015.0,025015.025.0025.070.0x x x u x x x x ≤≤⎧⎪-≤≤⎪=⎨-+≤≤⎪⎪≤≤⎩ ()0,00(70,0)0u 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 GSUSE MATHimplicit nonereal*8 A(1:51,1:51), B(1:51), X(1:51,1:51), Y(1:51)integer I,J,K,ITERREAL*8 M(49,4),N(49)REAL*8 E1,E2,EPSopen(6,file="RESULTS.TXT")open(7,file="tecplot.DAT")A=0.0d0B=0.0d0X=0.0d0Y=0.0D0M=0.0d0N=0.0D0E1=0.0D0E2=1.0D0EPS=0.005D0ITER=0C---------------------------------------------------------C 边界条件C---------------------------------------------------------DO I=1,51X(I,1)=120.0D0X(I,51)=0.0D0ENDDODO J=1,51X(1,J)=0.0D0X(51,J)=120.0D0ENDDOX(1,1)=120.0D0X(51,51)=0.0D0DO J=2,8X(1,J)=0.1D0*(9-J)*160.0D0 !进口边界，此处Y=a*x+(120-0)/0.75*y m/s ENDDODO J=50,44,-1X(51,J)=0.1D0*(51-J)*160.0D0 !出口边界，此处假设a=0.0d0m/sENDDOA=XDO WHILE(E2>EPS)E2=0.0D0DO I=2,50 !沿X方向线推进计算DO J=2,50IF(J.EQ.2)THENM(J-1,1)=0.0D0M(J-1,2)=1.0D0M(J-1,3)=-0.25D0M(J-1,4)=0.25*(X(I-1,J)+X(I,J-1)+X(I+1,J))ELSEIF(J.EQ.50)THENM(J-1,1)=-0.25D0M(J-1,2)=1.0D0M(J-1,3)=0.0D0M(J-1,4)=0.25*(X(I-1,J)+X(I+1,J)+X(I,J+1))ELSEM(J-1,1)=-0.25D0M(J-1,2)=1.0D0M(J-1,3)=-0.25D0M(J-1,4)=0.25*(X(I-1,J)+X(I+1,J))ENDIFENDDOCALL ZG(49,M,N)C DO K=1,49C WRITE(6,"(5F)")M(K,1),M(K,2),M(K,3),M(K,4),N(K) C ENDDOC PAUSEC------------------------------------------------------C 更新X方向的值，这样就可以线推进C------------------------------------------------------DO J=2,50X(I,J)=N(J-1)ENDDOENDDODO I=1,51DO J=1,51E1=ABS(X(I,J)-A(I,J))IF(E1>E2)THENE2=E1ELSEE2=E2ENDIFENDDOENDDOA=XITER=ITER+1WRITE(*,*)E2,ITERENDDOWRITE(6,"(51F)")XWRITE(7,*)'V ARIABLES="X","Y","F"'WRITE(7,*)"ZONE I= 51 ,J= 6 ,F=POINT"DO J=1,51IF(J.EQ.1 .OR. J.EQ.11 .OR. J.EQ.21 .OR.1 J.EQ.31 .OR.J.EQ.41 .OR. J.EQ.51) THENDO I=1,51WRITE(7,"(3F9.4)")(I-1)*0.1D0,(J-1)*0.1D0,X(I,J)ENDDOENDIFENDDOCLOSE(7)CLOSE(6)END PROGRAMC------------------------------------------------------------------------------C 追赶法C------------------------------------------------------------------------------ MODULE MATHCONTAINSSUBROUTINE ZG(N1,M,X)implicit real*8(a-h,o-z)INTEGER::N1REAL*8 MXLX(N1,4),X(N1),M(N1,4)MXLX=MDO i=1,N1IF(i==1)thenMXLX(i,2)=MXLX(i,2)MXLX(i,4)=MXLX(i,4)/MXLX(i,2)ELSEMXLX(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) ENDIFMXLX(i,3)=MXLX(i,3)/MXLX(i,2)ENDDODO i=N1,1,-1IF(i==N1)thenX(i)=MXLX(i,4)ELSEX(i)=MXLX(i,4)-MXLX(i,3)*X(i+1) ENDIFENDDOEND SUBROUTINEEND MODULE附录2 PSOR方法程序C---------------------------------------------------------------------C 主程序C---------------------------------------------------------------------Program PSORimplicit nonereal*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,ITERREAL*8 M(49,4),N(49)REAL*8 E1,E2,EPSREAL*8 AF,BFopen(6,file="RESULTS.TXT")open(7,file="tecplot.DAT")A=0.0d0B=0.0d0X=0.0d0Y=0.0D0M=0.0d0N=0.0D0E1=0.0D0E2=1.0D0EPS=0.005D0ITER=0AF=1.9C---------------------------------------------------------C 边界条件C---------------------------------------------------------DO I=1,51X(I,1)=120.0D0X(I,51)=0.0D0ENDDODO J=1,51X(1,J)=0.0D0X(51,J)=120.0D0ENDDOX(1,1)=120.0D0X(51,51)=0.0D0DO J=2,8X(1,J)=0.1D0*(9-J)*160.0D0 !进口边界，此处Y=a*x+(120-0)/0.75*y m/s ENDDODO J=50,44,-1X(51,J)=0.1D0*(51-J)*160.0D0 !出口边界，由于V=0，所以a=0.0ENDDOA=XX1=XDO WHILE(E2>EPS)E2=0.0D0DO I=2,50 !采用高斯赛德尔迭代高斯赛德尔迭代作为预计算值DO J=2,50BF=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 !及时更新所计算的值，体现为高斯赛德尔迭代ENDDOENDDODO I=1,51DO J=1,51E1=ABS(X(I,J)-X1(I,J))IF(E1>E2)THENE2=E1ELSEE2=E2ENDIFENDDOENDDOX1=XA=XITER=ITER+1WRITE(*,*)E2,ITERENDDOWRITE(6,"(51F)")XWRITE(7,*)'V ARIABLES="X","Y","F"'WRITE(7,*)"ZONE I= 6 ,J= 51 ,F=POINT"DO J=1,51DO I=1,51IF(I.EQ.1 .OR. I.EQ.11 .OR. I.EQ.21 .OR.1 I.EQ.31 .OR.I.EQ.41 .OR. I.EQ.51) THENWRITE(7,"(3F9.4)")(I-1)*0.1D0,(J-1)*0.1D0,X(I,J)ENDIFENDDOENDDOCLOSE(7)CLOSE(6)END PROGRAM附录3 波动方程求解程序C----------------------------------------------------------------C 主程序C----------------------------------------------------------------PROGRAM MAINUSE MATHIMPLICIT NONEREAL*8 DT,DX,V !时间步长，空间步长，声速REAL*8 NP !方法选择，1为一阶上风法，2为Lax-Wendroff，3为BTCS隐式法REAL*8 UPRE(0:70),UNEX(0:70) !当前时层和下一时层的速度值INTEGER I,JREAL*8 TIMEREAL*8 XREAL*8 TIREAL*8 A(69,4),Y(69)TIME=0.0D0C---------------------------------------------------C 初值确定C---------------------------------------------------DO I=0,70X=DBLE(I)IF(X.LE.5.0D0)THENUPRE(I)=0.0D0ELSEIF(X.GE.5.0D0 .AND. X.LE.15.0D0)THENUPRE(I)=2.0D0*X-10.0D0ELSEIF(X.GE.15.0D0 .AND. X.LE.25.0D0)THENUPRE(I)=-2.0D0*X+50.0D0ELSEUPRE(I)=0.0D0ENDIFENDDOUNEX=UPREV=200.0D0c DT=0.005c DT=0.0025DT=0.00125DX=1.0NP=3OPEN(7,FILE="INITIAL.DAT")OPEN(8,FILE="RESULT.DAT")DO I=0,70WRITE(7,"(1F8.5,4X,1F8.5)")DBLE(I),UPRE(I)ENDDOCLOSE(7)IF(NP.EQ.1)THENWRITE(*,*)"数值格式为一阶上风法"WRITE(*,*)"DT=",DT,"DX=",DXELSEIF(NP.EQ.2)THENWRITE(*,*)"数值格式为Lax-Wendroff"WRITE(*,*)"DT=",DT,"DX=",DXELSEIF(NP.EQ.3)THENWRITE(*,*)"数值格式为BTCS隐式法"WRITE(*,*)"DT=",DT,"DX=",DXENDIFC-----------------------------------------------------TIME=TIME+DT !此处设置可以避免多算一步DO WHILE(TIME.LE.0.15D0)IF(NP.EQ.1)THEN !一阶上风法DO I=1,69UNEX(I)=(1.0D0-V*DT/DX)*UPRE(I)+V*DT/DX*UPRE(I-1) ENDDOELSEIF(NP.EQ.2)THENDO I=1,69UNEX(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))ENDDOELSEIF(NP.EQ.3)THENDO I=1,69IF(I.EQ.1)THENA(I,1)=0.0D0A(I,2)=1.0D0A(I,3)=V/2.0D0*DT/DXA(I,4)=V/2.0D0*DT/DX*UPRE(I-1)+UPRE(I)ELSEIF(I.EQ.69)THENA(I,1)=-V/2.0D0*DT/DXA(I,2)=1.0D0A(I,3)=0.0D0A(I,4)=-V/2.0D0*DT/DX*UPRE(I+1)+UPRE(I)ELSEA(I,1)=-V/2.0D0*DT/DXA(I,2)=1.0D0A(I,3)=V/2.0D0*DT/DXA(I,4)=UPRE(I)ENDIFENDDOCALL ZG(69,A,Y)DO I=1,69UNEX(I)=Y(I)ENDDOENDIFUPRE=UNEXTI=ABS(CEILING(TIME/0.025D0)-TIME/0.025D0)IF(TI.LE.1.0D-5)THENDO I=0,70WRITE(8,"(1F8.5,4X,1F8.5,4X,1F8.5)")TIME, DBLE(I),UNEX(I)ENDDOENDIFTIME=TIME+DTENDDOCLOSE(8)ENDPROGRAM其中，子程序ZG()见附录1。