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

1、题目:用Fortran语言编写程序解决理想流体的平面圆柱绕流问题，如下图所示。由于流动的对称性，可以只研究其中的四分之一区域，如图中 abcde所示。222222- 0， 22 0xyxy其边界条件如下表所小。流函数势函数在下边界abc上0 0 n在出口边界cd上 0 nconst在上边界de上2 0 n在进口边界ea上y y- Ux 1 n说明： 是切向流速， 一是法向流速。nn下面就流函数进行讨论，为便于分析，把边界条件写成： 在1上 其中：1为具有本质B、C的边界0 在2上2为具有自然B、C的边界n解题步骤:(1)写出r a jie pkhh积分表达式通过分部积分，可得:x222dxdy

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

3、)单元分析ei 代入r ajiepKnh积分表达式:dxdy y ye ds得单元有限元方程组为:AjeeFi (i=1,2,3)由ieaibiX qy(i=1,2,3),可得:bj，Ci于是:b1blGGbb*2b1b3c1c3b2blc2clb2b2c2c2b2 b3aqb3 blc3clb3b2c3c2b3 b3AjAeeFie2 g i ds自然B.C处理:由于自然边界条件 n0。(5)总体合成总体矩阵Anm ；单元矩阵行号整体矩阵行号，单元矩阵列号整体矩阵列号。(6)本质B.C处理即为了满足本质B.C,要对总体系数矩阵进行处理,具体处理方法见作业(7)解总体方程组，求出有关物理量解方

4、程组的方法见作业一，由iaibi x ci yi ifVxi ' Ci i , vyy ybi i；(i1,2,3)三结点三角形单元，线性插值函数,每个单元只有一个流速,与单元内坐标无关,可理解为单元平均流速，位于单元中心(三中线交点)。源程序如下:其中绘图子程序说明：程序的部分说明作业一、二、三中已有，这里不再赘述;在作业三中也已有，这里略去。program yzrl implicit none interface subroutine linear_equation_bc(n,a1,b,x_result)integer二i,j,k,imaxinteger,intent(in):nr

5、eal:max,c real,dimension(:,:),intent(in):a1real,dimension(:),intent(in):breal,dimension(:,:),allocatable:a,mreal,dimension(:),intent(inout):x_resultend subroutine linear_equation_bcend interfacecharacter*12 file,name,ly*8integer*2 lengthinteger二i,j,k,dy,dyz,jd,jd1,jd2,jdzinteger:m,n,m0integer,dimens

6、ion(:,:),allocatable:zt integer二a,b,c,jj,kkreal,parameter二pi=3.1415926536real, dimension(:),allocatable:x,yreal:x1,y1,r!定义单元、节点编号!定义单元结点整体编号数组!定义pi为常量，值为圆周率!定义整体结点坐标数组!定义网格划分区域real:bcx1,bcx2,bcy1,bcx,bcy,bcyh!定义 x, y 方向及圆弧段计算步长real:x0,y0!定义原点坐标real:ux!定义来流速度integer:z!定义本质B.C点个数integer, dimension(:),

7、allocatable:bcjd!定义本质结点编号real, dimension(:),allocatable:bcz!定义本质 B.C 值real,dimension(3,3):dya!定义单元系数矩阵real,dimension(:,:),allocatable:zta!整体系数矩阵real,dimension(:,:),allocatable:bb,cc!定义基函数系数矩阵real:dym,d!定义单元三角形面积!定义结点流函数值real, dimension(:),allocatable:cs!定义常数项数组 real, dimension(:),allocatable:freal,

8、dimension(:),allocatable:vx,vy,v!定义单元 x、y 方向流速及合成流速real, dimension(:),allocatable:zx,zy,zxx,zyy!定义单元中心及流速终点坐标real, dimension(:),allocatable:jx!定义合成流速与 x 方向夹角real, dimension(:),allocatable:jtx1,jty1,jtx2,jty2 !箭头坐标print*,'请输入x方向等分数m, y方向等分数n,圆弧等分数m0'read*,m,n,m0print*,'请输入网格划分部分尺寸x1 , y1,

9、 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(mod(j,2)=0)then!计算偶数单元的节点对应关系数组zt(dy,1)=j/2+n+(n+1)*(i-1)+1zt(dy,2)=j/2+(n+1)*(i-1)+1zt(dy,3)=j/2+n+(n+1)*(i-1)+2else!计算奇数单元的节点对应关系数组zt(dy,1)=j/2+(n+1)*(

10、i-1)+1zt(dy,2)=j/2+(n+1)*(i-1)+2zt(dy,3)=j/2+n+(n+1)*(i-1)+2end ifend doend doopen(1file='单元中整体结点编号.txt')write(1,'(a,i4,a,3i5)'),('第',i,'个单元:',(zt(i,j),j=1,3),i=1,dyz)close(1)bcx1=x1/m!定义x1边的等分步长bcx2=(x1-r)/(m-m0)!定义x2边的等分步长bcy1=y1/n !定义y1边的等分步长do i=1,m-m0+1do j=1,n+1

11、jd=(n+1)*(i-1)+jx(jd)=(i-1)*bcx1+(i-1)*(j-1)*(bcx2-bcx1)/ny(jd)=y1-(j-1)*bcy1end doend do bcyh=pi/(2.0*m0)!定义圆弧段的等分角大小do i=m-m0+2,m+1 jd1=(n+1)*(i-1)+1 jd2=(n+1)*i x(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-(

12、m-m0+1)*bcyh)!计算圆弧段 m0等分点的各个y坐标bcx=(x(jd2)-x(jd1)/n!计算该计算区域的x方向的步长bcy=(y(jd2)-y(jd1)/n!计算该计算区域的y方向的步长do j=2,n jd=(n+1)*(i-1)+j x(jd)=(i-1)*bcx1+(j-1)*bcx y(jd)=y1+(j-1)*bcy end doend doopen(2file='整体结点坐标.txt')!将整体结点坐标保存在txt文档中write(2,'(a,i3,a,f8.4,3x,f8.4)'),('第',i,'点坐标为:

13、',x(i),y(i),i=1,jdz)close(2)print*,'请输入图形名称read(*,'(a)')namelength=len_trim(name)!用内部函数求出name去掉尾部空格后的长度file=name(:length)/'.dxf!连接字符串"name"与".dxf"open(3,file=file)call staplot !调用绘制dxf文件的开头子程序ly='wg'length=len_trim(ly)x0=0!给定绘图原点y0=0do i=1,dyza=zt(i,1)

14、b=zt(i,2)c=zt(i,3)call line(x(a),y(a),x(b),y(b),ly(:LENGTH)!调用绘图子程序画线call line(x(b),y(b),x(c),y(c),ly(:LENGTH)call line(x(c),y(c),x(a),y(a),ly(:LENGTH)end docall textc(x(jdz-n)+0.1,y(jdz-n)+0.05,0.03,file,ly(:length),0.0)print*,'请输入来流速度ux'read*,uxz=2*(m+1)+n-1!计算具有本质B.C节点个数allocate(bcjd(z),b

15、cz(z)do i=1,z!定义本质B.C,找到具有本质B.C的节点和其对应的值if(i<=n+1)thenbcjd(i)=ibcz(i)=bcy1*(n+1-i)*uxelse if(i<=m+n+1)thenbcjd(i)=(i-n)*(n+1)bcz(i)=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)'),序号','对应序号整体编号,'本质

16、B.C值'do i=1,zwrite(4,'(i3,10x,i3,12x,f6.3)'),i,bcjd(i),bcz(i)end doclose(4)allocate(zta(jdz,jdz),bb(dyz,3),cc(dyz,3)zta=0!确定单元基函数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=d/2.0!计算单元三角形面积dym,编号均为逆

17、时针，故面积均为正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)/dcc(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(5file='单元系数矩阵.txt')write(5,'(a,i3)'),'单元',id

18、o 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,*),(dya(j,k),k=1,3)!输出单元系数矩阵end do end do close(5) open(6file='整体系数矩阵.txt')write(6,*),(zta(i,j),j=1,jdz),i=1,jdz)!输出整体系数矩阵close(6)allocate(cs(jdz),f(jdz) cs=0

19、do k=1,z!处理本质B.Ccs(bcjd(k)=bcz(k)do j=1,jdzif(j=bcjd(k)thenzta(j,j)=1 elsezta(bcjd(k),j)=0end ifend do end doopen(7,file='本质B.C处理后的整体系数矩阵.txt')write(7,*),(zta(k,j),j=1,jdz),k=1,jdz)!输出本质 B.C 处理后的整体系数矩阵close(7)open(8file='方程组右侧矢量项.txt')do k=1,jdz!输出方程组右侧矢量项write(8,'(a,i3,a,f6.3)&#

20、39;),'第',k,'行矢量项',cs(k)end doclose(8)call linear_equation_bc(jdz,zta,cs,f)!调用子程序解线性方程组open(9file='结点流函数值.txt')do k=1,jdzwhte(9,'(a,i3,a,f10.6)'),'第',k,'点函数值',f(k)end doclose(9)allocate(vx(dyz),vy(dyz),v(dyz)open(10,file='单元 x、y 向流速.txt')open(11

21、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(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(1

22、1,'(a,i3,a,f10.6)'),'第',i,'个单元流速',v(i)end doclose(10)close(11)allocate(zx(dyz),zy(dyz),zxx(dyz),zyy(dyz),jx(dyz)allocate (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

23、(zt(i,2)+y(zt(i,3)/3if(vx(i)=0)then!计算流速与x轴的夹角jx(i)=pi/2.0elsejx(i)=atan(vy(i)/vx(i)end ifwrite(12,'(a,i3,a,f10.6)'),'第',i,'单元流速与 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 line(zx(i),zy(i),zxx(i),zyy(i),ly(:length) !画流线

24、jtx1(i)=zxx(i)-sin(-jx(i)+8*pi/18.0)*v(i)*0.03 jty1(i)=zyy(i)-cos(-jx(i)+8*pi/18.0)*v(i)*0.03 call line(zxx(i),zyy(i),jtx1(i),jty1(i),ly(:length) jtx2(i)=zxx(i)-cos(jx(i)-pi/18.0)*v(i)*0.03 jty2(i)=zyy(i)-sin(jx(i)-pi/18.0)*v(i)*0.03call line(zxx(i),zyy(i),jtx2(i),jty2(i),ly(:length) end doclose(12)call endplotclose(3)end program yzrlsubroutine linear_equation_bc(n,a1,b,x_result)implicit noneinteger二i,j,k,imaxinteger,intent(in):nreal:max,creal,dimension(:,:),intent(in):a1 real,dimension(:),intent(in):b real,dimension(:,:),allocatable:a,m r