流函数-涡量法的二维方腔流数值模拟

流函数- 涡量法的二维方腔流数值模拟

基本方程：

22v t

ξ

γξψξ∂=∇-∇∂-∇=⎧⎨⎩

t ———时间步长, s;

γ———流体的运动粘度, m2 /s;

ξ———涡量, s- 1;

ψ———流函数, m2 /s;

v ———速度矢量, m/s;

差分格式：

采用FTCS 格式，对于ξ有：

1,,1,1,,1,1

1,,1,,1,,1

,2

2

22(

)22()

()

n n n n

n n

n n n

n n n

i j i j

i j i j

i j i j i j i j i j

i j i j i j n

n

i j

i j

u

v

t

x

y

x y ξξξξξξξξξξξξυ++-+-+-+-+----+-+=--++

∆∆∆∆∆

采用FTCS 格式，对于ξ有：

1,,1,,1,,1,,1

,2

2

22)()()n n

n n n n n n i j i j

i j i j i j i j i j i j i j t

x y ψψψψψψψψξ++-+---+-+=

+

+∆∆∆

边界条件：

速度的边界条件：固体壁面处设置无滑移边界条件, 即u=0, v=0; 流函数的边界条件：在固体壁面及平板驱动处的流函数ξ=0; 涡量的边界条件： 平板驱动处：2u

x

ξ=

∆ 左右壁面上:;

,1,,2

2()()

i j i j i j y ξξξ+-=

∆

网格划分：

采用等距结构化网格划分（40*40）

编程计算：

本算例采用MATLAB进行编译，其主要优势是语言简单，可以方便地描绘出方腔环流的等值线图等。

主要语句：

while norm(c1)>1e-4|norm(c2)>1e-4

n=n+1;

O1=O;

E1=E;

t=t+dt;

for i=2:I

for j=2:J

O(i,j)=O(i,j)-dt*(u(i,j)*(O(i+1,j)-O(i-1,j))/(2*dx)+v(i,j)*(O(i,j+1)-O(i,j-1))/ (2*dy))+dt/re*((O(i+1,j)-2*O(i,j)+O(i-1,j))/dx^2+(O(i,j+1)-2*O(i,j)+O(i,j-1))/d y^2);

end

end

for i=2:I

for j=2:J

E(i,j)=E(i,j)+0.25*(E(i+1,j)+E(i-1,j)+E(i,j+1)+E(i,j-1)-4*E(i,j)+dx^2*O(i,j)); end

end

% 0为涡量，E流函数%

计算效果图：

分别设置雷诺数为200，500,1000

雷诺数Re=200时的流函数图及速度矢量图

雷诺数Re=500时的流函数图及速度矢量图

雷诺数Re=1000时的流函数图及速度矢量图