Fortran语言——涡量流函数法中心差分格式的二维方腔顶盖驱动计算解析

涡量流函数法中心差分格式的二维方腔顶盖驱动计算
本题是关于粘性流体方腔顶盖驱动的问题。

采用涡量流函数法，在均分网格下，用中心差分格式进行计算，结果与文献中所采用的其他方法和格式进行比较，认为中心差分格式符合计算的精度，但是其明显缺点是计算过程不稳定。

1、参数的无量纲化
令 X =x H; Y =y H; U =u u t ; ⁄⁄⁄V =v u t ; ⁄Ω=w w 0 ; ⁄Ψ=φφ0⁄ 其中，由涡量的定义式：u v y x
ω∂∂=
-∂∂ 将速度的导数项无量纲化，并将上述定义式代入，即
ω=∂u ∂y −∂v ∂x =u t ∂U H ∂Y −u t ∂V H ∂X =u t H (∂U ∂Y −∂V ∂X
)
令u
t H ⁄=w 0 则涡量的无量纲形式可以写为： Ω=
∂U ∂Y
−
∂V ∂X
=
ωω0
再由流函数的定义式：u =∂ψ/∂y 同理，令 ψ0=u t ∗H
可以得到流函数的无量纲表达式为：Ψ=ψ
ψ0
流函数边界条件
(X,Y )∈AB Ψ=0
(X,Y )∈BC Ψ=0 (X,Y )∈CD Ψ=0 (X,Y )∈AD Ψ=0
涡量的边界条件（Thom 公式）
(X,Y )∈AB Ω1，j =
2(Ψ2,j −Ψ1,j )
δX 2 (X,Y )∈BC Ωi ，1=
2(Ψi,2−Ψi,1)
δY (X,Y )∈CD ΩL1,j =
2(ΨL2,j −ΨL1,j )
δX 2 (X,Y )∈AD Ωi ，M1=
2(Ψi,M2−Ψi,M1)
δY 2
+ 2
δY
2、方程的无量纲处理过程
2.1、非守恒型方程的处理
（1）将以上假设的各式代入到非守恒型方程组：
22
22()u v x y x y ωωωωρρμ∂∂∂∂+=+∂∂∂∂
2222
0x y ψψ
ω∂∂+-=∂∂
得到以下无量纲形式的方程： ①涡量控制方程
∂2Ω2+∂2Ω2=R e (U ∂Ω+V ∂Ω
) ②流函数控制方程
∂2Ψ∂X 2+∂2Ψ
∂Y 2
−Ω=0 （2）采用有限差分法离散非守恒型涡量的无量纲控制方程：
R e (U P ΩE −ΩW 2ΔX +V P ΩN −ΩS 2ΔY )=ΩW −2ΩP +ΩE ∆X 2+ΩS −2ΩP +ΩN
∆Y 2
整理得到：(2
∆X 2+2
∆Y 2)ΩP =(1
∆X 2+R e U P
2ΔX
)ΩW +(1
∆X 2−
R e U P 2ΔX )ΩE
+(1
∆Y 2+R e V P 2ΔY
)ΩS +(1
∆Y 2−R e V P
2ΔY
)ΩN
简化为：4ΩP =(1+R e U P ΔX
2)ΩW +(1−R e U P ΔX
2)ΩE
+(1+
R e V P ΔY
2
)ΩS +(1−
R e V P ΔY
2
)ΩN
写成一般形式是: a p Ωp =a W ΩW +a E ΩE +a S ΩS +a N ΩN 其中 a p =4 ；a W =(1+R e U P ΔX
2) ；a E =(1−
R e U P ΔX
2)
a S =(1+
R e V P ΔY
2
) ； a N =(1−
R e V P ΔY
2
)
2.2、守恒型方程的处理
（1）将假设的各个无量纲量代入守恒型涡量方程 得到其无量纲形式为：
Re[∂(UΩ)∂X +∂(VΩ)∂Y ]=∂2Ω∂X 2+∂2Ω∂Y
2
（2）采用有限容积法离散守恒型涡量的无量纲控制方程：
对于扩散项，有
∫∫[ ∂e
w
n s (∂Ω)+∂(∂Ω
)]dXdY
=
∂Ω∂X |w e ∆Y +∂Ω∂Y |s
n
∆X =[ΩE −ΩP (δX )e −ΩP −Ωw (δX)w ]∆Y +[ΩN −ΩP (δY )n −ΩP −ΩS
(δY)s
]∆X
对于对流项，有
∫∫[∂(UΩ)∂X +∂(VΩ)
∂Y ]e
w
n s dXdY =[(UΩ)e −(UΩ)w ]∆Y +[(VΩ)n −(VΩ)s ] ∆X
对于界面上的取值，采用如下形式：
U e ≥0 时，Ωe =ΩP +(f e +−ΩP )； U e ≤0 时，Ωe =ΩE +(f e −
−ΩE ) U w ≥0时，Ωw =ΩW +(f e +−ΩW )； U w ≤0 时，Ωw =ΩP +(f e −−ΩP ) V n ≥0时，Ωn =ΩP +(f e +−ΩP )； V n ≤0 时，Ωn =ΩN +(f e −−ΩN ) V s ≥0时，Ωs =ΩS +(f e +−ΩS )； V s ≤0 时，Ωs =ΩP +(f e −−ΩP )
其中，界面上涡量的插值采用中心差分格式，令
f e +=f e −=ΩP +ΩE
2 ；f w +=f w −
=ΩP +Ωw
2；
f n +=f n −=
ΩP +ΩN
2
；f s +=f s −=ΩP +ΩS
2
所以，Ωe =
ΩP +ΩE
2
Ωw =ΩP +ΩW
2
Ωn =
ΩP +ΩN
2
Ωs =
ΩP +ΩS
2
关于界面上速度的插值，可以用上一层次的流函数来表示，形式如下： U e =
ΨN +ΨNE −ΨS −ΨSE
4∆Y U w =
ΨNW +ΨN −ΨS −ΨSW
4∆Y
V n =
ΨNW +ΨW −ΨE −ΨNE 4∆X
V s =
ΨW +ΨSW −ΨSE −ΨE
4∆X
由以上各式，整理得到守恒型涡量控制方程无量纲形式的离散方程为：
R e [
ΨN +ΨNE −ΨS −ΨSE 4∆Y
ΩP +ΩE
2
−
ΨNW +ΨN −ΨS −ΨSW 4∆Y
ΩP +ΩW
2
]∆Y
+R e [ΨNW +ΨW −ΨE −ΨNE 4∆X
ΩP +ΩN
2
−
ΨW +ΨSW −ΨSE −ΨE 4∆X
ΩP +ΩS
2
] ∆X
= [
ΩN −ΩP
(δY )n
−
ΩP −ΩS (δY)s
]∆X +[
ΩE −ΩP (δX )e
−
ΩP −Ωw (δX)w
]∆Y
综上，可整理得到（同位的均分网格）：
4ΩP = (1−
R e 8
(ΨN +ΨNE −ΨS −ΨSE ))ΩE +(1+
R e 8
(ΨNW +ΨN −ΨS −ΨSW ))ΩW
+(1−R e
8
(ΨNW+ΨW−ΨE−ΨNE))ΩN+(1+
R e
8
(ΨW+ΨSW−ΨSE−ΨE))ΩS
写成一般形式为：
b pΩp=b EΩE+b WΩW+b NΩN+b SΩS 其中，b p=4
b E=1−R e
8
(ΨN+ΨNE−ΨS−ΨSE)
b W=1+R e
8
(ΨNW+ΨN−ΨS−ΨSW)
b N=1−R e
8
(ΨNW+ΨW−ΨE−ΨNE)
b S=1+R e
8
(ΨW+ΨSW−ΨSE−ΨE)
2.3、对于流函数方程，离散可得到：
(2
2
+
2
2
)ΨP=
1
2
ΨW+
1
2
ΨE+
1
2
ΨS+
1
2
ΨN−ΩP
整理并写成一般形式是：c pΨp=c WΨW+c EΨE+c SΨS+c NΨN−d
其中，c p=4 ；c W=1 ；c E=1 ；c S=1 ；c N=1；d=ΩP∙∆X2
并且，采用有限容积法离散的结果与采用这里的方法得到的结果是相同的。

2.4、将压力的泊松方程进行无量纲化并离散
令 P0=2ρu t2
代入压力泊松方程，得到其无量纲形式为：
∂2P ∂X2+∂2P
∂Y2
=[∂2Ψ
∂X2
∙∂2Ψ
∂Y2
−(∂2Ψ
∂X∂Y
)
2
]
边界上的压力满足N-S方程，同样将其无量纲化，并与上式组成压力的封闭方程
组，X方向上：∂P
∂X =1
R e
(∂2U
∂X2
+∂2U
∂Y2
)−(U∂U
∂X
+V∂U
∂Y
)
Y方向上：∂P
∂Y =1
R e
(∂2V
∂X2
+∂2V
∂Y2
)−(U∂V
∂X
+V∂V
∂Y
)
将以上三式所组成的方程组离散，得到：
内点压力计算式为：(2
∆X2+2
∆Y2
)P P=1
∆X2
P W+1
∆X2
P E+1
∆Y2
P S+1
∆Y2
P N
−[ΨE−2ΨP+ΨW
∆X2
∙ΨN−2ΨP+ΨS
∆Y2
−(ΨNE−ΨSE−ΨNW+ΨSW
4∆X∆Y
)2]
左边界压力为：1
∆X P P=1
∆X
P E−1
R e
(U P−2U E+U EE
∆X2
+U N−2U P+U S
∆Y2
)
+(U P∙U E−U P
∆X
+V P∙
U N−U P
∆Y
)
右边界压力为：1
∆X P P=1
∆X
P W+1
R e
(U P−2U W+U WW
∆X2
+U N−2U P+U S
∆Y2
)−(U P∙
U P−U W
+V P∙
U N−U P
)
上边界压力为：1
∆Y P P=1
∆Y
P S+1
R e
(V E−2V P+V W
∆X2
+V P−2V S+V SS
∆Y2
)−(U P∙V E−V P
∆X
)
下边界压力为：1
∆Y P P=1
∆Y
P N−1
R e
(V E−2V P+V W
∆X
+V P−2V N+V NN
∆Y
) +(U P∙
V E−V P
∆X
+V P∙
V N−V P
∆Y
)
由此可得到压力的封闭方程组，通过求得的速度场来计算压力场。

3、
4、计算结果
（1）有限差分法的情况下： A 流函数等值线（图1a/b/c ）
N=20×20
N=40×
40
N=80×80
B 涡量等值线（图2a/b/c ）
N=20×20 N=40×40 N=80×80
C 压力等值线（图3a/b/c ）
N=20×20 N=40×40 N=80×80
（2）有限容积法的情况下：
A 流函数等值线（图5a/b/c）
N=20×20N=40×40N=80×80 B 涡量等值线（图6a/b/c）
N=20×20N=40×40N=80×80 C 压力等值线（图7a/b/c）
N=20×20N=40×40N=80×80
4.3、雷诺数为1000时各物理量场的分布图形
（1）有限差分法的情况下：
A 流函数等值线（图9a/b/c）
N=20×20N=40×40N=80×80 B 涡量等值线（图10a/b/c）
N=20×20N=40×40N=80×80 C 压力等值线（图11a/b/c）
N=20×20N=40×40N=80×80
D 中心线的速度分布（图12a/b/c）
（2）有限容积法的情况下：
A 流函数等值线（图13a/b/c）
N=20×20N=40×40N=80×80 B 涡量等值线（图14a/b/c）
N=20×20N=40×40N=80×80
C 压力等值线（图15a/b/c）
N=20×20N=40×40N=80×80 D 中心线的速度分布（图16a/b/c）
4.4、文献[1]上的结果
（1）雷诺数100时候的流函数和涡量图形（图17a/b）
（2）雷诺数1000时的流函数和涡量图形（图18a/b）
5、讨论与小结
在图4a/b、图8a/b、图12a/b/c、图16a/b/c中，将计算方腔的中心线速度分布与U.GHIA等人用CSI-MG方法在129×129网格数情况下得到的结果[1]相比较，可以看出本题的计算结果与基准解符合得比较好，说明计算方法具有合适的精度。

但是，本题在整个计算过程中，所遇到的关于稳定性及收敛性方面的问题较多，正如文献[9]所讲，在涡量流函数的计算中，由于涡量控制方程是对流占优扩散的，使用中心差分的离散方法不能得到很好的计算结果。

在低雷诺数RE=100时，计算很容易就得到符合基准解的结果。

从中心线上的速度分布可以看到，结果与基准解吻合。

但是，在较高雷诺数1000时，由于中心差分格式是条件稳定[10]，这造成了计算的不稳定，并且非常容易发散。

关于涡量的上边界条件，若按照文献[9]上的形式: Ω=2/Y 来做，则得到的速度会出现越界，但是这时速度的图形却更接近高雷诺数时候的基准解，且涡量和流函数也更加接近高雷诺数时候的基准解；若采用题目的边界条件，则得出的速度均有物理意义，且涡量及流函数的图形比较符合低雷诺数的解。

这说明涡量边界条件的确定对于计算的准确性有很大的影响。

（如图19a/b）
a. 采用文献上的涡量边界在
b. 采用文献上的涡量边界在
RE=100时候的涡量等值线图RE=100时候的流函数等值线图
本题在计算雷诺数RE=1000的时候，过程是明显不稳定的。

将松弛因子取得很小（0.1），造成收敛性变差，计算很费时。

原因之一是采用高斯赛德尔迭代，且内迭代次数取200次，而外迭代次数在结果收敛时（80×80网格）竟达到115832次之多。

本题程序只适合用于计算雷诺数较低且网格数较密的情况。

参考文献
[1] U Ghia, K N Ghia and C T Shin. High-Re solutions for incompressible flow using
the Navier-Stokes equations and a multigrid method. Journal of Computational Physics. 1982, 48: 387-411
[2] E.Erturk1.Numerical solutions of 2-D steady incompressible driven cavity flow at
high Reynolds numbers. Journal for Numerical Methods in Fluids. 2005,
48:747–774
[3] Yih-Ferng Peng,Yuo-Hsien Shiau,Robert R.Hwang.Transition in a 2-D lid-driven
cavity puters & Fluids.2003;32:337–352
[4] Patankar SV. Numerical heat transfer and fluid flow. New York: McGraw-Hill;
1980.
[5] Charles-Henri Bruneau , Mazen Saad. The 2D lid-driven cavity problem revisited.
Computers & Fluids.2006;35:326–348
[6] 李明秀,陶文铨,王秋旺. LATTICE-BOLTZMANN方法及其在顶盖驱动流数值模拟
中的应用.工程热物理[J].2001,3.22(2):223-225
[7] 周高领,陈斌,宇波.不同形状空腔顶盖驱动流的数值模拟.中国工程热物理学
会学术会议论文.
[8] 陈庆光.张永建.钱宝光.QUICK和乘方格式在顶盖驱动方腔流动数值计算中的
比较[J]-山东科技大学学报(自然科学版) 2002,3.21(1):8-11
[9] 吴晓冬,陈立亮.流函数-涡量法的二维方腔流数值模拟.中国铸造装备与技术
[J].2007,1:36-38
[10] 王贤钢,俞茂铮,陶文铨. 一种稳定性与精度协调一致的二阶差分格式.西安
交通大学学报.2005,11.39(11):1194-1198
源程序：
（1）有限差分法源程序：
! THIS PROGRAM IS WRITTEN TO SOLVE A 2-D DRIVEN FLOW PROBLEM IN A SQUARE ! CAVITY USING THE VORTICITY-STREAM FUNCTION METHOD.
! L1 X方向的节点数
! M1 y方向的节点数
! DX X方向的步长
! DY Y方向的步长
!PHAI 流函数
!OMEGA 涡量
!RE 雷诺数
!ITS 迭代次数
!ERR1，ERR2 前后迭代步各点涡量之差
PROGRAM MAIN
IMPLICIT NONE
REAL DX,DY
REAL ERR1,ERR2,ERR3
INTEGER,PARAMETER ::L1=81,M1=81
REAL ::RE=100.0
INTEGER i,j,K
INTEGER ITS !计算涡量时的跌代数
INTEGER ITSS !计算压力时的跌代数
INTEGER N !内迭代次数
INTEGER ::ITSMAX =5000 !最大跌代数
REAL ::EPS =1E-6 !迭代精度
REAL ::ALPHA1 =0.3 !松弛因子
REAL ::ALPHA2 =0.3
REAL U(L1,M1),V(L1,M1) !速度分量
REAL P(L1,M1) !压力分布
REAL XL(L1,M1),YM(L1,M1) !存放每个点的横坐标和纵坐标
REAL PHAI(L1,M1)
REAL OMEGA(L1,M1)
REAL AW(L1,M1),AE(L1,M1),AS(L1,M1),AN(L1,M1) !涡量方程中的各项系数
REAL OME_K(L1,M1) !存放每一迭代步得到的涡量值
REAL PHAI_K(L1,M1) !存放每一迭代步得到的流函数值
REAL P_K(L1,M1) !存放压力泊松方程迭代时上一步的值DOUBLE PRECISION TIME_START,TIME_END
! *...............................................................................
CALL CPU_TIME(TIME_START)
CALL MESHING(L1,M1,DX,DY,PHAI,OMEGA,U,V,p,XL,YM)
ERR1=1.0
ITS=0
DO WHILE(ERR1>=EPS.OR.ERR2>=EPS.OR.ITS<ITSMAX)
DO I=1,L1
DO J=1,M1
OME_K(I,J)=OMEGA(I,J) !用OME_K储存前一迭代步的计算结果
PHAI_K(I,J)=PHAI(I,J) !用PHAI_K储存前一迭代步的计算结果
END DO
END DO
DO I=2,L1-1 !涡量方程系数的更新
DO J=2,M1-1
AW(I,J)=1.0+RE*U(I,J)*DX/2.0
AE(I,J)=1.0-RE*U(I,J)*DX/2.0
AS(I,J)=1.0+RE*V(I,J)*DY/2.0
AN(I,J)=1.0-RE*V(I,J)*DY/2.0
END DO
END DO
DO N=1,200
DO J=M1-1,2,-1 !解涡量方程
DO I=2,L1-1
OMEGA(i,j)=(AW(i,j)*OMEGA(i-1,j)+AE(i,j)*OMEGA(i+1,j)&
+AS(i,j)*OMEGA(i,j-1)+AN(i,j)*OMEGA(i,j+1))/4.0
OMEGA(I,J)=OME_K(I,J)+ALPHA1*(OMEGA(I,J)-OME_K(I,J))
END DO
END DO
END DO
WRITE(*,*) 'OMEGA(3,11)=',OMEGA(3,11)
DO N=1,200
DO J=M1-1,2,-1 !解流函数方程
DO I=2,L1-1
PHAI(i,j)=(PHAI(i-1,j)+PHAI(i+1,j)&
+PHAI(i,j-1)+PHAI(i,j+1)-OMEGA(i,j)*DX**2)/4.0
PHAI(I,J)=PHAI_K(I,J)+ALPHA2*(PHAI(I,J)-PHAI_K(I,J))
END DO
END DO
END DO
WRITE(*,*) 'PHAI(3,11)=',PHAI(3,11)
DO I=2,L1-1 !由流函数计算速度分布
DO J=2,M1-1
U(I,J)=(PHAI(I,J+1)-PHAI(I,J))/DY+(2.0*PHAI_K(I,J)-PHAI_K(I,J+1)-PHAI_K(I,J-1))/(
2.0*DY)
V(I,J)=(PHAI(I-1,J)-PHAI(I,J))/DY+(2.0*PHAI_K(I,J)-PHAI_K(I+1,J)-PHAI_K(I-1,J))/(
2.0*DX)
END DO
DO j=1,M1 !左、右边界的涡量值
OMEGA(1,J)=2.0*(PHAI(2,J)-PHAI(1,J))/DX**2
OMEGA(L1,J)=2.0*(PHAI(L1-1,J)-PHAI(L1,J))/DX**2
END DO
DO i=2,L1-1 !上、下边界的涡量值
OMEGA(i,M1)=2.0/DY
OMEGA(i,1)=2.0*(PHAI(i,2)-PHAI(i,1))/DY**2
END DO
OMEGA(1,1)=(OMEGA(1,2)+OMEGA(2,1))/2.0
OMEGA(1,M1)=(OMEGA(1,M1-1)+OMEGA(2,M1))/2.0
OMEGA(L1,1)=(OMEGA(L1-1,1)+OMEGA(2,L1))/2.0
OMEGA(L1,M1)=(OMEGA(L1,M1-1)+OMEGA(L1-1,M1))/2.0
ITS=ITS+1
ERR1=MAXVAL(ABS(OME_K-OMEGA))
ERR2=MAXVAL(ABS(PHAI_K-PHAI))
WRITE(*,*) 'ERR1=',ERR1
WRITE(*,*) 'ERR2=',ERR2
WRITE(*,*) 'ITS=',ITS
END DO
! *......................................................................................................
ERR3=1.0
ITSS=0 !用计算得到的速度求解内点压力场分布DO WHILE(ERR3>=EPS.AND.ITSS<=ITSMAX)
DO I=1,L1
DO J=1,M1
P_K(I,J)=P(I,J)
END DO
END DO
DO I=2,L1-1
DO J=2,M1-1
P(I,J)=(P(I-1,J)+P(I+1,J)+P(I,J+1)+P(I,J-1)-0.25*((U(I+1,J)-U(I-1,J))*(V(I,J+1)-V(I,J-1))-& (V(I+1,J)-V(I-1,J))*(U(I,J+1)-U(I,J-1))))/4.0
P(I,J)=P_K(I,J)+0.3*(P(I,J)-P_K(I,J))
END DO
END DO
! *............................................................................................................
DO J=2,M1-1 !压力场的左、右边界值
P(1,J)=4.0/3.0*P(2,J)-1.0/3.0*P(3,J)-(2.0/(3.0*RE))*(-2.0*U(2,J)+U(3,J))/D
P(L1,J)=4.0/3.0*P(L1-1,J)-1.0/3.0*P(L1-2,J)+(2.0/(3.0*RE))*(-2.0*U(L1-1,J)+U
(L1-2,J))/DX
END DO
DO I=2,L1-1 !压力场的上、下边界值
P(I,M1)=(4.0/3.0)*P(I,M1-1)-(1.0/3.0)*P(I,M1-2)+(2.0/(3.0*RE*DY))*(-2.0*V(I,M1-
1)+V(I,M1-2))
P(I,1)=4.0/3.0*P(I,M1-1)-1.0/3.0*P(I,M1-2)-(2.0/(3.0*RE))*(-2.0*V(I,2)+V(I,3))/DX END DO
!角点上的压力进行加权平均
P(1,1)=0.5*P(1,2)+0.5*P(2,1)
P(1,M1)=0.5*P(1,M1-1)+0.5*P(2,M1)
P(L1,1)=0.5*P(L1,2)+0.5*P(L1-1,1)
P(L1,M1)=0.5*P(L1,M1-1)+0.5*P(L1-1,M1)ITSS=ITSS+1
ERR3=MAXVAL(ABS(P_K-P))
END DO
CALL CPU_TIME(TIME_END)
WRITE(*,*)'总共运行时间(S)',TIME_END-TIME_START
! *..........................................................................................................
OPEN(1,FILE='PHAI.DAT')
OPEN(2,FILE='OMEGA.DAT')
OPEN(3,FILE='V.DAT')
OPEN(4,FILE='U.DAT')
OPEN(5,FILE='PRESSURE.DAT')
WRITE(1,*) 'TITLE="PHAI"'
WRITE(1,*) 'VARIABLE="X","Y","PHAI"'
WRITE(1,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
WRITE(2,*) 'TITLE="OMEGA"'
WRITE(2,*) 'VARIABLE="X","Y","OMEGA"'
WRITE(2,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
WRITE(3,*) 'TITLE="V"'
WRITE(3,*) 'VARIABLE="X","Y","V"'
WRITE(3,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
WRITE(4,*) 'TITLE="U"'
WRITE(4,*) 'VARIABLE="X","Y","U"'
WRITE(4,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
WRITE(5,*) 'TITLE="PRESSURE"'
WRITE(5,*) 'VARIABLE="X","Y","PRESSURE"'
WRITE(5,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
DO I=1,L1
DO J=1,M1
WRITE(1,*) XL(I,J),YM(I,J), PHAI(I,J)
WRITE(2,*) XL(I,J),YM(I,J), OMEGA(I,J)
WRITE(3,*) XL(I,J),YM(I,J), V(I,J)
WRITE(4,*) XL(I,J),YM(I,J), U(I,J)
WRITE(5,*) XL(I,J),YM(I,J), P(I,J)
END DO
END DO
CLOSE (1)
CLOSE (2)
CLOSE (3)
CLOSE (4)
CLOSE (5)
OPEN(6,FILE='V_X.TXT')
OPEN(7,FILE='U_Y.TXT')
DO I=1,L1
WRITE(6,*) V(I,(M1+1)/2.0)
END DO
DO J=1,M1
WRITE(7,*) U((L1+1)/2.0,J)
END DO
CLOSE (6)
CLOSE (7)
! *............................................................................................................ CONTAINS
!网格划分并且设定初场
SUBROUTINE MESHING(L,M,X,Y,PHI,OME,U_X,V_Y,P0,XL,YM)
INTEGER L,M
REAL X,Y,PHI(L,M),OME(L,M),U_X(L,M),V_Y(L,M),p0(L,M),XL(L,M),YM(L,M) X=1.0/(L-1)
Y=1.0/(M-1)
xL(1,:)=0
xL(L,:)=1
YM(:,1)=0
YM(:,M)=1
DO I=1,L
DO J=1,M
XL(I,J)=XL(1,J)+(I-1)*X
YM(I,J)=YM(I,1)+(J-1)*Y
END DO
END DO
DO i=1,L !流函数、压力初场DO j=1,M
PHI(i,j)=0
P0(I,J)=1.0
END DO
END DO
DO i=1,L !涡量初场
DO j=1,M-1
OME(i,j)=0
END DO
END DO
DO i=1,L
OME(i,M)=2.0/Y
END DO
DO I=1,L !速度初场
U_X(I,1)=0
V_Y(I,1)=0
U_X(I,M)=1.0
V_Y(I,M)=0
END DO
DO J=2,M-1
DO I=1,L
U_X(I,J)=0
V_Y(I,J)=0
END DO
END DO
END SUBROUTINE
END PROGRAM
（2）有限容积法源程序：
! THIS PROGRAM IS WRITTEN TO SOLVE A 2-D DRIVEN FLOW PROBLEM IN A SQUARE ! CAVITY USING THE VORTICITY-STREAM FUNCTION METHOD.
! L1 X方向的节点数
! M1 y方向的节点数
! DX X方向的步长
! DY Y方向的步长
!PHAI 流函数
!OMEGA 涡量
!RE 雷诺数
!ITS 迭代次数
!ERR1，ERR2 前后迭代步各点涡量之差
PROGRAM MAIN
IMPLICIT NONE
REAL DX,DY
REAL ERR1,ERR2
INTEGER,PARAMETER ::L1=21,M1=21
INTEGER :: RE=100
INTEGER i,j
INTEGER ITS !跌代数
INTEGER ::ITSMAX =1000 !最大跌代数
REAL ::EPS =1E-6 !迭代精度
REAL N !内迭代次数
REAL ::ALPHA =0.1 !松弛因子
REAL PHAI(L1,M1)
REAL OMEGA(L1,M1)
REAL AP(L1,M1),AW(L1,M1),AE(L1,M1),AS(L1,M1),AN(L1,M1) !一般形式的涡量方程中的各项系数
REAL OME_K(L1,M1) !存放每一迭代步得到的涡量值
REAL PHAI_K(L1,M1) !存放每一迭代步得到的流函数值
REAL CP,CW,CE,CS,CN !一般形式的流函数方程中的各项系数DOUBLE PRECISION TIME_START,TIME_END
CALL MESHING(L1,M1,DX,DY,PHAI,OMEGA)
CP=2.0/DX**2+2.0/DY**2
CW=1.0/DX**2
CE=1.0/DX**2
CS=1.0/DY**2
CN=1.0/DY**2
ERR1=1
ERR2=1
ITS=0
CALL CPU_TIME(TIME_START)
CALL MESHING(L1,M1,DX,DY,PHAI,OMEGA,U,V,p,XL,YM)
ERR1=1.0
ERR2=1.0
ITS=0
DO WHILE(ERR1>=EPS.OR.ERR2>=EPS.AND.ITS<ITSMAX)
DO I=1,L1
DO J=1,M1
OME_K(I,J)=OMEGA(I,J) !用OME_K储存前一迭代步的计算结果
PHAI_K(I,J)=PHAI(I,J) !用PHAI_K储存前一迭代步的计算结果
END DO
END DO
DO I=2,L1-1 !涡量方程系数的更新
DO J=2,M1-1
AW(I,J)=1.0+RE*(U(i-1,j)+U(i,j))*DX/4.0
AE(I,J)=1.0-RE*(U(i+1,j)+U(i,j))*DX/4.0
AS(I,J)=1.0+RE*(V(i,j-1)+V(i,j))*DY/4.0
AN(I,J)=1.0-RE*(V(i,j+1)+V(i,j))*DY/4.0
END DO
END DO
DO N=1,200 !内迭代次数取200
DO i=2,L1-1 !解涡量方程
DO j=2,M1-1
OMEGA(i,j)=(AW(i,j)*OMEGA(i-1,j)+AE(i,j)*OMEGA(i+1,j)&
+AS(i,j)*OMEGA(i,j-1)+AN(i,j)*OMEGA(i,j+1))/AP(i,j) END DO
END DO
END DO
DO N=1,200
DO i=2,L1-1 !解流函数方程
DO j=2,M1-1
PHAI(i,j)=(CW*PHAI(i-1,j)+CE*PHAI(i+1,j)&
+CS*PHAI(i,j-1)+CN*PHAI(i,j+1)-OMEGA(i,j))/CP
END DO
END DO
END DO
ITS=ITS+1
ERR1=MAXVAL(ABS(OME_K-OMEGA))
ERR2=MAXVAL(ABS(PHAI_K-PHAI))
DO j=1,M1 !用亚松弛处理左、右边界的涡量值OMEGA(1,J)=2*(PHAI(2,J)-PHAI(1,J))/DX**2
OMEGA(L1,J)=2*(PHAI(L1-1,J)-PHAI(L1,J))/DX**2
OMEGA(1,J)=OME_K(1,j)+ALPHA*(OMEGA(1,J)-OME_K(1,J))
OMEGA(L1,J)=OME_K(L1,J)+ALPHA*(OMEGA(L1,J)-OME_K(L1,J))
END DO
DO i=2,L1-1 !用亚松弛处理上、下边界的涡量值OMEGA(i,1)=2*(PHAI(i,2)-PHAI(i,1))/DY**2
OMEGA(i,M1)=2*(PHAI(i,M1-1)-PHAI(i,M1))/DY**2
OMEGA(i,1)=OME_K(i,1)+ALPHA*(OMEGA(i,1)-OME_K(i,1))
OMEGA(i,M1)=OME_K(i,M1)+ALPHA*(OMEGA(i,M1)-OME_K(i,M1))
END DO
END DO
! *...............................................................................
ERR3=1.0
ITSS=0 !用计算得到的速度求解内点压力场分布DO WHILE(ERR3>=EPS.AND.ITSS<=ITSMAX)
DO I=1,L1
DO J=1,M1
P_K(I,J)=P(I,J)
END DO
END DO
DO I=2,L1-1
DO J=2,M1-1
P(I,J)=(P(I-1,J)+P(I+1,J)+P(I,J+1)+P(I,J-1)-0.25*((U(I+1,J)-U(I-1,J))*(V(I,J+1)-V(I,J-1))-& (V(I+1,J)-V(I-1,J))*(U(I,J+1)-U(I,J-1))))/4.0
P(I,J)=P_K(I,J)+0.3*(P(I,J)-P_K(I,J))
END DO
END DO
! *...............................................................................
DO J=2,M1-1 !压力场的左、右边界值
P(1,J)=4.0/3.0*P(2,J)-1.0/3.0*P(3,J)-(2.0/(3.0*RE))*(-2.0*U(2,J)+U(3,J))/D
P(L1,J)=4.0/3.0*P(L1-1,J)-1.0/3.0*P(L1-2,J)+(2.0/(3.0*RE))*(-2.0*U(L1-1,J)+U
(L1-2,J))/DX
END DO
DO I=2,L1-1 !压力场的上、下边界值
P(I,M1)=(4.0/3.0)*P(I,M1-1)-(1.0/3.0)*P(I,M1-2)+(2.0/(3.0*RE*DY))*(-2.0*V(I,M1-
1)+V(I,M1-2))
P(I,1)=4.0/3.0*P(I,M1-1)-1.0/3.0*P(I,M1-2)-(2.0/(3.0*RE))*(-2.0*V(I,2)+V(I,3))/DX END DO
!角点上的压力进行加权平均
P(1,1)=0.5*P(1,2)+0.5*P(2,1)
P(1,M1)=0.5*P(1,M1-1)+0.5*P(2,M1)
P(L1,1)=0.5*P(L1,2)+0.5*P(L1-1,1)
P(L1,M1)=0.5*P(L1,M1-1)+0.5*P(L1-1,M1)ITSS=ITSS+1
ERR3=MAXVAL(ABS(P_K-P))
END DO
CALL CPU_TIME(TIME_END)
WRITE(*,*)'总共运行时间(S)',TIME_END-TIME_START
! *..........................................................................................................
OPEN(1,FILE='PHAI.DAT')
OPEN(2,FILE='OMEGA.DAT')
OPEN(3,FILE='V.DAT')
OPEN(4,FILE='U.DAT')
OPEN(5,FILE='PRESSURE.DAT')
WRITE(1,*) 'TITLE="PHAI"'
WRITE(1,*) 'VARIABLE="X","Y","PHAI"'
WRITE(1,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
WRITE(2,*) 'TITLE="OMEGA"'
WRITE(2,*) 'VARIABLE="X","Y","OMEGA"'
WRITE(2,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
WRITE(3,*) 'TITLE="V"'
WRITE(3,*) 'VARIABLE="X","Y","V"'
WRITE(3,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
WRITE(4,*) 'TITLE="U"'
WRITE(4,*) 'VARIABLE="X","Y","U"'
WRITE(4,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
WRITE(5,*) 'TITLE="PRESSURE"'
WRITE(5,*) 'VARIABLE="X","Y","PRESSURE"'
WRITE(5,*) 'ZONE T="S"', 'I=',L1, 'J=',M1, 'C=BLACK'
DO I=1,L1
DO J=1,M1
WRITE(1,*) XL(I,J),YM(I,J), PHAI(I,J)
WRITE(2,*) XL(I,J),YM(I,J), OMEGA(I,J)
WRITE(3,*) XL(I,J),YM(I,J), V(I,J)
WRITE(4,*) XL(I,J),YM(I,J), U(I,J)
WRITE(5,*) XL(I,J),YM(I,J), P(I,J)
END DO
END DO
CLOSE (1)
CLOSE (2)
CLOSE (3)
CLOSE (4)
CLOSE (5)
OPEN(6,FILE='V_X.TXT')
OPEN(7,FILE='U_Y.TXT')
DO I=1,L1
WRITE(6,*) V(I,(M1+1)/2.0)
END DO
DO J=1,M1
WRITE(7,*) U((L1+1)/2.0,J)
END DO
CLOSE (6)
CLOSE (7)
!*…………………………………………………………………………………………………………CONTAINS
SUBROUTINE MESHING(L,M,X,Y,PHI,OME) !网格划分并且设定初场INTEGER L,M
REAL X,Y,PHI(L,M),OME(L,M)
X=1.0/(L-1)
Y=1.0/(M-1)
DO i=1,L
DO j=1,M
PHI(i,j)=0
END DO
END DO
DO i=1,L
DO j=1,M-1
OME(i,j)=0
END DO
END DO
DO i=1,L
OME(i,M)=2.0/Y
END DO
END SUBROUTINE
END PROGRAM。