fortran理想流体的平面圆柱绕流程序

1、扬州大学理想流体的平面圆柱绕流计算水力学题目用Fortran语言编写程序解决理想流体的平面圆柱绕流问题，如下图所示。由于流动的对称性，可以只研究其中的四分之一区域，如图中abcde所示。2 2 2 2水研2010级第17页共15页水利学院其边界条件如下表所示流函数势函数在下边界abc上0 0 n在出口边界Cd上0nCon St在上边界de上2 0 n在进口边界ea上y yUX1n22 0， 2 2 0XyXy说明： 是切向流速， 一是法向流速。nn下面就流函数进行讨论，为便于分析，把边界条件写成: 在1上 其中：1为具有本质B、C的边界-0 在2上2为具有自然B、C的边界n解题步骤：写出 aQ

2、 puh积分表达式2 22 一2 dxdy 0 X y通过分部积分，可得：dxdy g ds y y2(2)区域剖分横向剖分数为9,纵向剖分数为10,其中圆弧段剖分数为5。利用作业三中的程序实现(由于网格内要画流速矢量图,故单元编号未写出)，另外，还需要建立本质B.C表。(3)确定单元基函数ie设网格划分后任意三角形单元的三个结点的坐标值别为(Xie, y)(i1,2,3),函数(e)值分别为ie(i 1,2,3),根据基函数的构造思想，单元内近似函数可表示为式:在单元内作线性插值函数如下：毎 b?x qy;3e3 b3XC3ye1a b1eqy; 2 根据基函数的插值条件，得到系数：a,b,

3、G(i 1,2,3)。则基函数为：eiibx Ciy,i 1,2,3。(4)单元分析ee将iei 代入 aepuH积分表达式：eeeeedxdyge dsXXy y2得单元有限元方程组为：Aje jeeFi (i=1,2,3)eiie(i 1,2,3)O由ieai于是:ejbj，eiXey-cj,-Xbb1C1C1b1b2C?b1b3c1c3eAb2b1c2c1b2b2C2C2b2b3c2c3,b3b1c3c1b3b2C3C2b3b3c3c3biX cy(i=1,2,3),可得:AjeFieiyCieg i ds自然B.C处理:由于自然边界条件 一nO ,则 Fie(5)总体合成单元矩阵Aj

4、e总体矩阵Anm ；单元矩阵行号整体矩阵行号，单元矩阵列号号,整体矩阵列号。(6)本质B.C处理即为了满足本质B.C ,要对总体系数矩阵进行处理,具体处理方法见作业二。(7)解总体方程组，求出有关物理量解方程组的方法见作业一，由 iaibi XCiyi i得:VXiCi i , vyy ybi i；(i1,2,3)三结点三角形单元，线性插值函数,每个单元只有一个流速,与单元内坐标无关,可理解为单元平均流速，位于单元中心(三中线交点)。源程序如下:说明：程序的部分说明作业一、二、三中已有，这里不再赘述；其中绘图子程序 在作业三中也已有，这里略去。PrOgram yzrlimplicit none

5、in terfaceSubrOutine Iinear_equation_bc(n,a1,b,x_result)in teger:i,j,k,imaxin teger,i nten t(i n): nreal:max,c real,dime nsio n(:,:),i nte nt(i n):a1real,dime nsio n(:),i nten t(i n):breal,dime nsio n(:,:),allocatable:a,mreal,dime nsio n(:),i nten t(i no ut):x_resultend SUbrOUt ine Iin ear_equati on

6、_bcend in terfaceCharaCter*12 file, name,ly*8in teger*2 Ien gthin teger:i,j,k,dy,dyz,jd,jd1,jd2,jdzin teger:m, n,min teger,dime nsio n(:,:),allocatable:Ztin teger:a,b,c,jj,kkreal,parameter:pi=3.1415926536 real, dime nsio n(:),allocatable:x,y real:x1,y1,r!定义单元、节点编号!定义单元结点整体编号数组!定义Pi为常量，值为圆周率!定义整体结点坐标

7、数组!定义网格划分区域real:bcx1,bcx2,bcy1,bcx,bcy,bcyh !定义 X，y 方向及圆弧段计算步长 real:xO,yO!定义原点坐标real:UX!定义来流速度integer:Z!定义本质B.C点个数in teger, dime nsio n( :),allocatable:bcjd!定义本质结点编号real, dimension(:),allocatable:bcz!定义本质 B.C 值real,dime nsio n(3,3):dya!定义单元系数矩阵real,dime nsio n(:,:),allocatable:Zta!整体系数矩阵real,dime ns

8、io n(:,:),allocatable:bb,cc !定义基函数系数矩阵 real:dym,d!定义单元三角形面积!定义结点流函数值real, dime nsion (:),allocatable:CS!定义常数项数组 real, dime nsio n(:),allocatable:freal, dimension(:),allocatable:vx,vy,v!定义单元 X、y 方向流速及合成流速real, dime nsion (:),allocatable:zx,zy,zxx,zyy!定义单元中心及流速终点坐标real, dimension(:),allocatable:jx!定义合

9、成流速与 X 方向夹角real, dime nsio n(:),allocatable:jtx1,jty1,jtx2,jty2!箭头坐标Print*,'请输入X方向等分数m, y方向等分数n,圆弧等分数m'read*,m, n,mPrint*,'请输入网格划分部分尺寸x1，y1，r'read*,x1,y1,rdyz=m*n *2!计算单元总数jdz=(m+1)*( n+1)!计算结点总数allocate(x(jdz),y(jdz),zt(dyz,3)do i=1,m!计算局部节点编号与整体节点编号的关系数组do j=1,2* ndy=(i-1)*2* n+jif

10、(mod(j,2)=0)the n!计算偶数单元的节点对应关系数组zt(dy,1)=j2+n+( n+1)*(i-1)+1zt(dy,2)=j2+( n+1)*(i-1)+1zt(dy,3)=j2+n+( n+1)*(i-1)+2else!计算奇数单元的节点对应关系数组zt(dy,1)=j2+( n+1)*(i-1)+1 zt(dy,2)=j2+( n+1)*(i-1)+2 zt(dy,3)=j2+n+( n+1)*(i-1)+2end ifend do end do open（1,file='单元中整体结点编号.txt'）Write(1,'(a,i4,a,3i5)&#

11、39;),('第',i,'个单元:',(zt(i,j),j=1,3),i=1,dyz)close(1)bcx1=x1m!定义x1边的等分步长bcx2=(x1-r)(m-m0)!定义x2边的等分步长bey仁y1n !定义y1边的等分步长do i=1,m-m+1do j=1, n+1jd=( n+1)*(i-1)+jx(jd)=(i-1)*bcx1+(i-1)*(j-1)*(bcx2-bcx1)/ny(jd)=y1-(j-1)*bey1end doend dobcyh=pi(2.0*m0)!定义圆弧段的等分角大小do i=m-m0+2,m+1jd 1=(n+1)*(

12、i-1)+1jd2=( n+1)*ix(jd1)=(i-1)*bcx1!计算与圆弧段对应的x1边上的等分点的X坐标y(jd1)=y1!计算与圆弧段对应的x1边上的等分点的y坐标x(jd2)=x1-r*cos(i-(m-m0+1)*bcyh)!计算圆弧段 m0等分点的各个X坐标y(jd2)=r*sin(i-(m-m0+1)*bcyh)!计算圆弧段 m0等分点的各个y坐标bcx=(x(jd2)-x(jd1)/n!计算该计算区域的X方向的步长bcy=(y(jd2)-y(jd1)n!计算该计算区域的y方向的步长do j=2 ,njd=( n+1)*(i-1)+jx(jd)=(i-1)*bcx1+(j-

13、1)*bcxy(jd)=y1+(j-1)*bcyend doend doopen(2,file='整体结点坐标.txt')!将整体结点坐标保存在txt文档中Write(2,'(a,i3,a,f84,3x,f84)'),('第',i,'点坐标为:',x(i),y(i),i=1,jdz)close (2)Print*,'请输入图形名称read(*,'(a)' )n ameIength=len_trim(name)!用内部函数求出name去掉尾部空格后的长度file=name(:Iength)'.dxf!

14、连接字符串"name"与".dxf"OPe n(3,file=file)call staplot !调用绘制dxf文件的开头子程序ly='wg'Ien gth=le n_trim(Iy)x0=0!给定绘图原点y0=0do i=1,dyza=zt(i,1)b=zt(i,2)c=zt(i,3)call li ne(x(a),y(a),x(b),y(b),ly(:LENGTH)!调用绘图子程序画线call li ne(x(b),y(b),x(c),y(c),ly(:LENGTH)call li ne(x(c),y(c),x(a),y(a),ly

15、(:LENGTH)end docall textc(x(jdz-n )+0.1,y(jdz-n)+0.05,0.03,file,ly(:Ie ngth),0.0)Print*,'请输入来流速度ux'read*,uxz=2*(m+1)+n-1!计算具有本质B.C节点个数allocate(bcjd(z),bcz(z)do i=1,z!定义本质B.C,找到具有本质B.C的节点和其对应的值if(i<=n+1)the nbcjd(i)=ibcz(i)=bcy1*( n+1-i)*uxelse if(i<=m+n+1)the nbcjd(i)=(i-n)*( n+1)bcz(i

16、)=0.0elsebcjd(i)=( n+1)*(z-i+1)+1bcz(i)=y1*uxend ifend doOPen(4,file='本质 B.C 表.txt') !建立本质 B.C 表Write(4,'(a,3x,a,3x,a)'),序号','对应序号整体编号','本质B.C值do i=1,zWrite(4,'(i3,10x,i3,12x,f6.3)'),i,bcjd(i),bcz(i)end doclose(4)allocate(zta(jd z,jdz),bb(dyz,3),cc(dyz,3)zta=0

17、!确定单元基函数do i=1,dyzd=x(zt(i,2)*y(zt(i,3)+x(zt(i,3)*y(zt(i,1)+x(zt(i,1)*y(zt(i,2)-x(zt(i,2)*y(zt(i,1)-x(zt(i,1)*y(zt(i,3)-x(zt(i,3)*y(zt(i,2)dym=d2.0!计算单元三角形面积dym,编号均为逆时针，故面积均为正bb(i,1)=(y(zt(i,2)-y(zt(i,3)d!计算基函数系数bb(i,2)=(y(zt(i,3)-y(zt(i,1)dbb(i,3)=(y(zt(i,1)-y(zt(i,2)dcc(i,1)=(x(zt(i,3)-x(zt(i,2)dc

18、c(i,2)=(x(zt(i,1)-x(zt(i,3)dcc(i,3)=(x(zt(i,2)-x(zt(i,1)d!计算各单元系数矩阵，由于自然B.C为0,故Fi(e)为0open(5,file='单元系数矩阵.txt')Write(5,'(a,i3)'),'单元',ido j=1,3!总体系数矩阵的合成jj=zt(i,j)do k=1,3kk=zt(i,k)dya(j,k)=(bb(i,j)*bb(i,k)+cc(i,j)*cc(i,k)*dymzta(jj,kk)=zta(jj,kk)+dya(j,k)end doWrite(5,*),(dy

19、a(j,k),k=1,3)!输出单元系数矩阵end doend doclose(5)OPen(6,file='整体系数矩阵.txt')Write(6,*),(zta(i,j),j=1,jdz),i=1,jdz)!输出整体系数矩阵close (6)allocate(cs(jdz),f(jdz)CS=Odo k=1,z!处理本质B.Ccs(bcjd(k)=bcz(k)do j=1,jdzif(j=bcjd(k)the nzta(j,j)=1elsezta(bcjd(k),j)=0end ifend doend doopen(7,file='本质B.C处理后的整体系数矩阵.t

20、xt')Write(7,*),(zta(k,j),j=1,jdz),k=1,jdz)!输出本质 B.C 处理后的整体系数矩阵close(7)open(8,file='方程组右侧矢量项.txt')do k=1,jdz!输出方程组右侧矢量项Write(8,'(a,i3,a,f6.3)'),'第',k,'行矢量项',cs(k)end doclose(8)call li near_equation_bc(jdz,zta,cs,f)!调用子程序解线性方程组open(9,file='结点流函数值.txt')do k=1

21、,jdzWrite(9,'(a,i3,a,f106)'),'第',k,'点函数值',f(k)end doclose(9)allocate(vx(dyz),vy(dyz),v(dyz)OPen(10,file='单元 X、y 向流速.txt')OPen(11,file='单元流速.txt')do i=1,dyz !计算单元流速 vx(i)=cc(i,1)*f(zt(i,1)+cc(i,2)*f(zt(i,2)+cc(i,3)*f(zt(i,3) vy(i)=-(bb(i,1)*f(zt(i,1)+bb(i,2)*f(

22、zt(i,2)+bb(i,3)*f(zt(i,3) v(i)=sqrt(vx(i)*2+vy(i)*2)Write(10,'(a,i3,a,f10.6,2x,a,f10.6)'),第',i,'单元 X 向流速',vx(i),'y 向流速',vy(i)Write(11,'(a,i3,a,f10.6)'),'第',i,'个单元流速',v(i)end doclose(10)close(11)allocate(zx(dyz),zy(dyz),zxx(dyz),zyy(dyz),jx(dyz)allo

23、cate (jtx1(dyz),jty1(dyz),jtx2(dyz),jty2(dyz)open(12,file='夹角.txt')!绘制流场矢量图do i=1,dyz!找到各个单元的中心zx(i)=(x(zt(i,1)+x(zt(i,2)+x(zt(i,3)3zy(i)=(y(zt(i,1)+y(zt(i,2)+y(zt(i,3)3if(vx(i)=0)then!计算流速与X轴的夹角jx(i)=pi2.0elsejx(i)=ata n(vy(i)vx(i)end ifWrite(12,'(a,i3,a,f10.6)'),'第',i,'

24、单元流速与 X 轴夹角',jx(i) zxx(i)=zx(i)+cos(jx(i)*v(i)*0.1!流线终点 X 坐标zyy(i)=zy(i)+sin(jx(i)*v(i)*0.1!流线终点 y 坐标call li ne(zx(i),zy(i),zxx(i),zyy(i),ly(:le ngth)!画流线jtx1(i)=zxx(i)-si n(-jx(i)+8*pi18.0)*v(i)*0.03 jty1(i)=zyy(i)-cos(-jx(i)+8*pi18.0)*v(i)*0.03CaIlli ne(zxx(i),zyy(i),jtx1(i),jty1(i),ly(:Ie ngt

25、h) jtx2(i)=zxx(i)-cos(jx(i)-pi18.0)*v(i)*0.03 jty2(i)=zyy(i)-s in (jx(i)-pi18.0)*v(i)*0.03call li ne(zxx(i),zyy(i),jtx2(i),jty2(i),ly(:le ngth) end doclose(12)call en dplotclose (3)end PrOgram yzrlSUbrOUtine Iinear_equation_bc(n,a1,b,x_result) implicit nonein teger:i,j,k,imaxin teger,i nten t(i n):

26、nreal:max,creal,dime nsio n(:,:),i nte nt(i n):a1real,dime nsio n(:),i nten t(i n):breal,dime nsio n(:,:),allocatable:a,mreal,dime nsio n(:),i nten t(i no ut):x_resultallocate(a( n,n+1),m( n,n+1)do i=1, ndo j=1, n+1if(j<=n )the na(i,j)=a1(i,j)elsea(i,j)=b(i)end ifend do!箭头终点X坐标!箭头终点y坐标!画箭头!箭头终点X坐标!箭头终点y坐标!画箭头!解方程组子程序end dodo k =1, n-1max=abs(a(k,k)imax=kdo i=k+1, n+1if (abs(a(i,