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