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

扬州大学 理想流体的平面圆柱绕流 计算水力学题目：用fortran语言编写程序解决理想流体的平面圆柱绕流问题，如下图所示。由于流动的对称性，可以只研究其中的四分之一区域，如图中abcde所示。在理想流体的平面运动中，流函数和势函数均满足拉氏方程：，其边界条件如下表所示。流函数势函数在下边界上在出口边界上在上边界上在进口边界上说明： 是切向流速， 是法向流速。下面就流函数进行讨论，为便于分析，把边界条件写成： 在上 其中：为具有本质b、c的边界 在上 为具有自然b、c的边界解题步骤：(1)写出积分表达式通过分部积分，可得：(2)区域剖分横向剖分数为9，纵向剖分数为10，其中圆弧段剖分数为5。利用作业三中的程序实现(由于网格内要画流速矢量图，故单元编号未写出)，另外，还需要建立本质b.c表。(3)确定单元基函数设网格划分后任意三角形单元的三个结点的坐标值别为，函数值分别为，根据基函数的构造思想，单元内近似函数可表示为式：。 在单元内作线性插值函数如下：根据基函数的插值条件，得到系数：。则基函数为：，。(4)单元分析将代入积分表达式：得单元有限元方程组为：(i=1,2,3)由(i=1,2,3)，可得：于是：，自然b.c处理：由于自然边界条件，则。(5)总体合成单元矩阵总体矩阵；单元矩阵行号整体矩阵行号，单元矩阵列号整体矩阵列号。(6)本质b.c处理即为了满足本质b.c，要对总体系数矩阵进行处理，具体处理方法见作业二。(7)解总体方程组，求出有关物理量解方程组的方法见作业一，由 及 得：三结点三角形单元，线性插值函数，每个单元只有一个流速，与单元内坐标无关，可理解为单元平均流速，位于单元中心(三中线交点)。源程序如下：说明：程序的部分说明作业一、二、三中已有，这里不再赘述；其中绘图子程序在作业三中也已有，这里略去。program yzrlimplicit noneinterfacesubroutine linear_equation_bc(n,a1,b,x_result)integer:i,j,k,imax integer,intent(in):nreal:max,creal,dimension(:,:),intent(in):a1real,dimension(:),intent(in):b real,dimension(:,:),allocatable:a,m real,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,jdz !定义单元、节点编号integer:m,n,m0integer,dimension(:,:),allocatable:zt !定义单元结点整体编号数组integer:a,b,c,jj,kkreal,parameter:pi=3.1415926536 !定义pi为常量，值为圆周率real, dimension(:),allocatable:x,y !定义整体结点坐标数组real:x1,y1,r !定义网格划分区域real:bcx1,bcx2,bcy1,bcx,bcy,bcyh !定义x，y方向及圆弧段计算步长real:x0,y0 !定义原点坐标real:ux !定义来流速度integer:z !定义本质b.c点个数integer, dimension(:),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:f !定义结点流函数值real, 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，圆弧等分数m0read*,m,n,m0print*,请输入网格划分部分尺寸x1，y1，rread*,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)*(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(1,file=单元中整体结点编号.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+1jd=(n+1)*(i-1)+jx(jd)=(i-1)*bcx1+(i-1)*(j-1)*(bcx2-bcx1)/ny(jd)=y1-(j-1)*bcy1end doend dobcyh=pi/(2.0*m0) !定义圆弧段的等分角大小do i=m-m0+2,m+1jd1=(n+1)*(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-1)*bcxy(jd)=y1+(j-1)*bcyend doend doopen(2,file=整体结点坐标.txt) !将整体结点坐标保存在txt文档中write(2,(a,i3,a,f8.4,3x,f8.4),(第,i,点坐标为:,x(i),y(i),i=1,jdz)close(2)print*,请输入图形名称read(*,(a)namelength=len_trim(name) !用内部函数求出name去掉尾部空格后的长度file=name(:length)/.dxf !连接字符串name与.dxfopen(3,file=file)call staplot !调用绘制dxf文件的开头子程序ly=wglength=len_trim(ly)x0=0 !给定绘图原点y0=0do i=1,dyza=zt(i,1)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*,请输入来流速度uxread*,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)thenbcjd(i)=ibcz(i)=bcy1*(n+1-i)*ux else if(i=m+n+1)thenbcjd(i)=(i-n)*(n+1)bcz(i)=0.0 elsebcjd(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(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，编号均为逆时针，故面积均为正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(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,*),(dya(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=0do k=1,z ! 处理本质b.ccs(bcjd(k)=bcz(k)do j=1,jdzif(j=bcjd(k)thenzta(j,j)=1elsezta(bcjd(k),j)=0end ifend doend doopen(7,file=本质b.c处理后的整体系数矩阵.txt)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 linear_equation_bc(jdz,zta,cs,f) !调用子程序解线性方程组open(9,file=结点流函数值.txt)do k=1,jdzwrite(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,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(11,(a,i3,a,f10.6),第,i,个单元流速,v(i) end do close(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(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) !画流线jtx1(i)=zxx(i)-sin(-jx(i)+8*pi/18.0)*v(i)*0.03 !箭头终点x坐标jty1(i)=zyy(i)-cos(-jx(i)+8*pi/18.0)*v(i)*0.03 !箭头终点y坐标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 !箭头终点x坐标jty2(i)=zyy(i)-sin(jx(i)-pi/18.0)*v(i)*0.03 !箭头终点y坐标call line(zxx(i),zyy(i),jtx2(i),jty2(i),ly(:length) !画箭头end doclose(12)call endplotclose(3)end program yzrlsubroutine linear_equation_bc