计算流体驱动方腔程序

驱动方腔（粘性不可压流），方腔示意图如左下：

物理背景：如左图所示,装满黏性不可压流体的方腔上底面以速度U 运动,其它壁面固定,其内部流体将作类似湍流的运动.流体运动复杂度将依赖于以下要素:

(1)方腔形状 ；(2)U 的大小；(3)流体的黏性,即雷诺数Re 的大小。

本题只讨论理想流体在二维正方形空腔中的流动，如左图，作无量纲化处理后,方腔长度和上底面运动速度均取为1,通过差分数值模拟计算出运动的流函数，涡量函数，画出Re=100和200时的流线图和涡量图。并标出窝心

位置。

采用流函数-涡方法：

对于二维不可压缩的理想流体,流速 ),(u v u = 应满足如下方程组:：

（

1

）

连

续

方

程

：

0u =⋅∇

（1）

(2) N-S 方程：0)(2222=∂∂-∂∂+∂∂=∂∂+∂∂+∂∂x p

y

u x u u y u v x u u t u ρ （2）

0)(2222=∂∂-∂∂+∂∂=∂∂+∂∂+∂∂y

p

y u x u u y u v x v u t v ρ

（3）

引入流函数：),,(t y x ψ和),,(t y x Ω：

y u ∂∂=

ψ，x

v ∂∂-=ψ

，y u x v ∂∂-∂∂=Ω （4）

（4）代入（1）的流函数的方程：

）

＋－（2222y

x ∂∂∂∂=Ωψ

ψ； （5） 把（5）代入)2()3(x y ∂-∂，的涡量方程：

)(Re 12222y

x y v x u t ∂Ω∂+∂Ω∂=∂Ω∂+∂Ω∂+∂Ω∂ （6） 以上方法称之为流函数－涡方法。

数值方法：将区域分为N ×N 个方格单元，只有（N+1）×（N+1）格网格点要计算，令h=1/N 为空间步长，τ为时间步长。 1：流函数与速度的边界条件：

;

0.........;0,0,1,1....,........2,1,0,,,x 1,y x 0,,,,,,i ==========≤≤j i N i N i N i N i N i j v u else v u y if N j i jh y ih ψψ，

2涡量的边界条件： 固壁

上

的

涡

量

：

N j i h

h

h

H j

N j N j j j i N i ..........0,,2;2;2;/22

,1,2

,1,02

,10,,=-

=Ω-

=Ω-

=Ω=Ω-ψψψ

计算：1，方腔内的流函数：

采用中心五点差分格式对（5）进行离散并简化得：

1...........1,;4/)(,21,1,,1,1,-=+++-=-+-+N j i h j i j i j i j i j i j i ψψψψψψ

算出j i ,ψ后，对（4）采用中心差分离散可以求得u,v. 2,方腔内得涡量：

采用迎风格式对涡量方程（6）进行离散：

);22(Re

12

1

2/1,1,11,2

1,11,1,112/1,12/1,1

,12

/11

,2/11,,,h

h h

v h u n j i n j i n j i n j

i n j

i n j

i n j i n j i n j

i n i n j

i n j

i n

j

i n j i +-++++-++++-++++-+++Ω+Ω-Ω+Ω

+Ω-Ω=

Ω-Ω+Ω-Ω+Ω-Ωτ

程序实现过程如下：

clear all clc

jishi=cputime; timestep=2000;

[x,y]=meshgrid(0:1/100:1); Nx=length(x);Ny=length(y); dh=1/100;dt=0.6*dh; renold=200;

fluidfun=zeros(Nx,Ny); vortex=zeros(Nx,Ny); velocityx=zeros(Nx,Ny); velocityy=zeros(Nx,Ny); %%%初始条件 velocityx(:,Ny)=1; vortex(:,Ny)=-2/dh; %fluidfun(Ny,:)=0; %velocityy(Ny,:)=0; for t=1:timestep for i=2:Nx-1 for j=2:Ny-1

fluidfun(i,j)=(fluidfun(i+1,j)+fluidfun(i-1,j)+fluidfun(i,j+1)+fluidfun(i,j-1)+dh^2*vortex(i,j))/4;

end

end

for i=2:Nx-1

for j=2:Ny-1

velocityx(i,j)=(fluidfun(i,j+1)-fluidfun(i,j-1))/(2*dh);

velocityy(i,j)=(fluidfun(i-1,j)-fluidfun(i+1,j))/(2*dh);

end

end

for i=2:Nx-1

for j=2:Ny-1

a=renold*(abs(velocityx(i,j))+abs(velocityy(i,j)))/dh+4/dh^2+renold/(dt);

b=renold*abs(velocityx(i,j))/dh+1/dh^2;

c=renold*abs(velocityy(i,j))/dh+1/dh^2;

if (velocityx(i,j)>=0)&&(velocityy(i,j)>=0)

vortex(i,j)=((vortex(i+1,j)+vortex(i,j+1))/dh^2+b*vortex(i-1,j)+c*vortex(i,j-1)+renold*vortex(i,j)/ (dt))/a;

elseif (velocityx(i,j)>=0)&&(velocityy(i,j)<0)

vortex(i,j)=((vortex(i+1,j)+vortex(i,j-1))/dh^2+b*vortex(i-1,j)+c*vortex(i,j+1)+renold*vortex(i,j)/ (dt))/a;

elseif (velocityx(i,j)<0)&&(velocityy(i,j)>=0)

vortex(i,j)=((vortex(i-1,j)+vortex(i,j+1))/dh^2+b*vortex(i+1,j)+c*vortex(i,j-1)+renold*vortex(i,j)/ (dt))/a;

elseif (velocityx(i,j)<0)&&(velocityy(i,j)<0)

vortex(i,j)=((vortex(i-1,j)+vortex(i,j-1))/dh^2+b*vortex(i+1,j)+c*vortex(i,j+1)+renold*vortex(i,j)/ (dt))/a;

end

end

end

end

c1=max(abs(fluidfun));

b1=find(c1==max(c1));%%%在第几列

d1=find(abs(fluidfun(:,b1))==max(c1));%%%在第几行

figure

contour(x,y,fluidfun,10)

%text(0,0.2,['re=200时驱动方腔流线图']);

hold on

plot(d1*0.01,b1*0.01,'ks')

title('re=200时驱动方腔流线图');