二维方腔环流计算资料下载.pdf

二维方腔环流计算资料下载.pdf

《二维方腔环流计算资料下载.pdf》由会员分享，可在线阅读，更多相关《二维方腔环流计算资料下载.pdf（19页珍藏版）》请在冰点文库上搜索。

2.差分方程：

+1=14（,2+1,+1+,1+1+1,+,+1）+（11）,+1=14+1,+,+1+1,+1+,1+1+4（,+1,1+1）（+1,1,+1）（+1,1,+1）（,+1,1+1）+（11）,+1=,+22（1+2）+1,+1,+1+2,+2,1+1（）2（）,2（12）,（）,=2（+1,2,+1,（）2）（,+12,+,1（）2）（+1,+1+1,11,+1+1,14）2其中k为迭代次数，1为亚松弛因子,2为超松弛因子。

一般亚松弛因子1为011,超松弛因子2为122。

当松弛因子选择恰当，会大大加快迭代收敛的速度。

3.边界条件流函数边界条件：

=0，在四边涡量边界条件：

（1）左边（AD）：

+1=2（+1,）2

（2）右边（BC）：

+1=2（1,）2（3）底边（CD）：

+1=2（,+1,）2（4）顶边（AB）：

+1=2（,1,+）2压强边界条件：

=

（2）cos（,）+

（2）cos（,）4.用fortran编程（附件1）进行计算，求出不同雷诺数（400，,1000,3200等）下流场、2涡量场和压力场的数值解（调整控制收敛速度）；

再用matlab程序进行绘图（附件3），进行分析和比较。

（三）使用ADI方法和亚松弛迭代法对非定常问题进行求解1.概述：

每个时间步长内，使用ADI方法对涡量进行求解，再用亚松弛迭代法对流函数进行求解；

2.差分方程

（1）涡量（D-R格式）：

预报：

（,+1,1）4221,+1+（22+）,+1+（,+1,1）422+1,+1=,+42（+1,1,）（,+1,1）+2（,+12,+,1）校正：

（+1,1,）422,1+1+（22+）,+1+（+1,1,）422,+1+1=,42（,+1,1）（+1,1,）+2（+1,2,+1,）

（2）流函数（FTCS格式）：

（,+1）+1=14（,+12+（1,+1）+1+（,1+1）+1+（+1,+1）+（,+1+1）+（11）（,+1）其中n为时间步数，k为迭代次数，1为亚松弛因子，。

一般亚松弛因子1为01err.and.counterE1）E1=abs（U（i,j）-U0（i,j）if（abs（Z（i,j）-Z0（i,j）E1）E1=abs（Z（i,j）-Z0（i,j）enddoif（E1E2）E2=E1enddowrite（*,*）counter,E2,E1enddo!

endwhile流函数涡量计算完成!

计算速度场VX（:

N）=111doi=2,N-1doj=2,N-1VX（i,j）=（U（i,j+1）-U（i,j-1）/（2.0*h）VY（i,j）=-（U（i+1,j）-U（i-1,j）/（2.0*h）enddoenddo!

计算SPSP=0doi=2,N-1doj=2,N-1SP（i,j）=2.0*rhow*（U（i+1,j）-2.0*U（i,j）+U（i-1,j）*（U（i,j+1）-2.0*U（i,j）+U（i,j-1）/h*4.0-（U（i+1,j+1）-U（i+1,j-1）-U（i-1,j+1）+U（i-1,j-1）/4.0/h*2.0）*2.0）enddoenddo100p=0PN=0!

计算压强counter=0E2=1dowhile（E20.00001.and.counterE1）E1=abs（P（i,j）-P0（i,j）enddoif（E1E2）E2=E1enddowrite（*,*）counter,E2,E1enddowrite（*,*）计算完成open（unit=8,file=outputSOR.txt）doi=1,Ndoj=1,Nwrite（8,*）h*REAL（I）-h,h*REAL（J）-h,U（i,j）,Z（i,j）,P（i,j）enddoenddoendprogram

（2）ADI.f90文件!

使用说明：

更改参数dt（时间步长），Re（雷诺数），rho（松弛因子）和err（允许误差）!

运行成功后在程序目录下可以找到输出文件”output.txt”programmainimplicitnonereal:

a,b,cREAL,PARAMETER:

DT=0.05real,parameter:

Re=5000INTEGER,PARAMETER:

TN=500000!

1000-2100013real,parameter:

err=0.0001!

允许误差real,parameter:

rho=0.05!

松弛因子REAL,PARAMETER:

l=1.0REAL,PARAMETER:

h=0.01INTEGER,PARAMETER:

N=L/H+1INTEGER,parameter:

Tw=5000real:

E1,E2real:

U（N,N）,U0（N,N）,Z（N,N）,Z0（N,N）,Z1（N,N）,W（N,N）,G（N,N）,F（N,N）,Vx（N,N）,Vy（N,N）!

U记录流场，Z记录涡量场，Vx是水平速度，Vy是竖直速度。

integer:

I,J,t,counteropen（unit=8,file=outputADI.txt）open（unit=9,file=outputUV.txt）U=0Z=0!

迭代求解dot=1,TNG=0W=0F=0U0=UZ0=ZVx=0Vx（:

N）=1Vy=0!

第一步-预估!

系数a=1b=-2c=1!

计算涡量边界值doi=1,NZ1（i,N）=-2.0*（U（i,N-1）-U（i,N）+h）/H/HZ1（i,1）=-2.0*（U（i,2）-U（i,1）/H/HZ1（1,i）=-2.0*（U（2,i）-U（1,i）/H/HZ1（N,i）=-2.0*（U（N-1,I）-U（N,i）/H/Henddo!

计算涡量内点!

涡量追DOJ=2,N-1Doi=2,N-1a=（-Re*dt*（U（i,j+1）-U（i,j-1）/4.0-dt）/h/h14b=Re+2.0*dt/h/hc=（Re*dt*（U（i,j+1）-U（i,j-1）/4.0-dt）/h/hif（i=2）thenF（i,J）=Re*Z（i,j）+Re*（U（i+1,j）-U（i-1,j）*（Z（i,j+1）-Z（i,j-1）*dt/4.0/h/h+（Z（i,j+1）-2.0*Z（i,j）+Z（i,j-1）*dt/h/h-a*Z1（i-1,j）G（i,j）=F（i,j）/bW（i,j）=c/belseif（i=N-1）thenF（i,J）=Re*Z（i,j）+Re*（U（i+1,j）-U（i-1,j）*（Z（i,j+1）-Z（i,j-1）*dt/4.0/h/h+（Z（i,j+1）-2.0*Z（i,j）+Z（i,j-1）*dt/h/h-c*Z1（i+1,j）G（i,j）=（F（i,j）-G（i-1,j）*a）/（b-a*W（i-1,j）W（i,j）=C/（b-a*W（i-1,j）elseF（i,J）=Re*Z（i,j）+Re*（U（i+1,j）-U（i-1,j）*（Z（i,j+1）-Z（i,j-1）*dt/4.0/h/h+（Z（i,j+1）-2.0*Z（i,j）+Z（i,j-1）*dt/h/hG（i,j）=（F（i,j）-G（i-1,j）*a）/（b-a*W（i-1,j）W（i,j）=C/（b-a*W（i-1,j）endifENDDO!

i完成涡量追!

涡量赶DOi=N-1,2,-1if（i=N-1）thenZ1（i,j）=G（i,j）elseZ1（i,j）=G（i,j）-W（i,j）*Z1（i+1,j）endifenddo!

i完成涡量赶enddo!

j完成涡量第一步!

完成第一步!

第二步-校正!

计算涡量边界值doi=1,NZ（i,N）=-2.0*（U（i,N-1）-U（i,N）+h）/H/HZ（i,1）=-2.0*（U（i,2）-U（i,1）/H/HZ（1,i）=-2.0*（U（2,i）-U（1,i）/H/HZ（N,i）=-2.0*（U（N-1,I）-U（N,i）/H/Henddo!

涡量追DOi=2,N-1Doj=2,N-1a=（Re*dt*（U（i+1,j）-U（i-1,j）/4.0-dt）/h/hb=Re+2.0*dt/h/hc=（-Re*dt*（U（i+1,j）-U（i-1,j）/4.0-dt）/h/h15if（j=2）thenF（i,j）=Re*Z1（i,j）-（U（i,j+1）-U（i,j-1）*（Z1（i+1,j）-Z1（i-1,j）*Re*dt/4.0/h/h+（Z1（i+1,j）-2.0*Z1（i,j）+Z1（i-1,j）*dt/h/h-a*Z（i,j-1）G（i,j）=F（i,j）/bW（i,j）=c/belseif（j=N-1）thenF（i,J）=Re*Z1（i,j）-（U（i,j+1）-U（i,j-1）*（Z1（i+1,j）-Z1（i-1,j）*Re*dt/4.0/h/h+（Z1（i+1,j）-2.0*Z1（i,j）+Z1（i-1,j）*dt/h/h-c*Z（i,j+1）G（i,j）=（F（i,j）-G（i,j-1）*a）/（b-a*W（i,j-1）W（i,j）=C/（b-a*W（i,j-1）elseF（i,J）=Re*Z1（i,j）-（U（i,j+1）-U（i,j-1）*（Z1（i+1,j）-Z1（i-1,j）*Re*dt/4.0/h/h+（Z1（i+1,j）-2.0*Z1（i,j）+Z1（i-1,j）*dt/h/hG（i,j）=（F（i,j）-G（i,j-1）*a）/（b-a*W（i,j-1）W（i,j）=C/（b-a*W（i,j-1）endifENDDO!

涡量赶DOj=N-1,2,-1if（j=N-1）thenZ（i,j）=G（i,j）elseZ（i,j）=G（i,j）-W（i,j）*Z（i,j+1）endifenddo!

完成涡量求解!

迭代法求解流量counter=0E2=1dowhile（E2err.and.counterE1）E1=abs（U（i,j）-U0（i,j）enddoif（E1E2）E2=E1enddowrite（*,*）t,counter,E2,E116enddo!

enddowhile!

输出计算结果if（mod（t,Tw）=0）then!

计算速度doi=2,N-1doj=2,N-1VX（i,j）=（U（i,j+1）-U（i,j-1）/（2.0*h）VY（i,j）=-（U（i+1,j）-U（i-1,j）/（2.0*h）enddoenddoDOI=1,NDOJ=1,Nwrite（8,*）h*real（I）-h,h*real（J）-h,U（i,j）,Z（i,j）write（9,*）h*real（I）-h,h*real（J）-h,Vx（i,j）,Vy（i,j）ENDDOENDDOendifenddoend（3）Draw.m文件%说明：

该文件是用于绘制等流线、等涡线、等压线图的文件。

clearall;

clc;

num=importdata（outputSOR.txt）;

x=0:

0.01:

1;

y=0:

x1=num（:

1）;

y1=num（:

2）;

z=num（:

3）;

w=num（:

4）;

%p=num（:

5）;

forn=1:

1temp=z（length（x）*length（y）*（n-1）+1）:

（length（x）*length（y）*n）;

temp2=w（length（x）*length（y）*（n-1）+1）:

fori=1:

length（x）forj=1:

length（y）z1（i,j）=temp（i-1）*length（y）+j）;

w1（i,j）=temp2（i-1）*length（y）+j）;

%p1（i,j）=p（i-1）*length（y）+j）;

17endendfigure（3）,contour（y,x,w1,1000）;

colorbar;

title（Re=1000，涡量场）;

figure

（2）,contour（y,x,z1,100）;

title（Re=1000，流场）;

%figure

（1）,contour（y,x,p1,1000）;

%colorbar;

%title（Re=400，压力场）;

%title（num2str（n）%title（Re=10000,第,num2str（n*10000）,步）drawnow;

end（4）DrawUV.m文件%说明：

该文件是用于绘制速度矢量图文件clearall;

num=importdata（outputUV.txt）;

u=num（:

v=num（:

y,x=meshgrid（0:

1,0:

1）;

1temp=u（length（x）*length（y）*（n-1）+1）:

temp2=v（length（x）*length（y）*（n-1）+1）:

length（y）u1（i,j）=temp（i-1）*length（y）+j）;

v1（i,j）=temp2（i-1）*length（y）+j）;

UZ=sqrt（u1（i,j）2+v1（i,j）2）;

u1（i,j）=u1（i,j）/UZ;

v1（i,j）=v1（i,j）/UZ;

endendquiver（x,y,u1,v1）;

axis（0101）;

drawnow;

end18