fortran理想流體的平面圓柱繞流程序

1、題目：用fortran語言編寫程序解決理想流體的平面圓柱繞流問題，如下圖所示。由于流動的對稱性，可以只研究其中的四分之一區(qū)域，如圖中abcde所示。在理想流體的平面運動中，流函數(shù)和勢函數(shù)均滿足拉氏方程：，其邊界條件如下表所示。流函數(shù)勢函數(shù)在下邊界上在出口邊界上在上邊界上在進(jìn)口邊界上說明： 是切向流速， 是法向流速。下面就流函數(shù)進(jìn)行討論，為便于分析，把邊界條件寫成： 在上 其中：為具有本質(zhì)b、c的邊界 在上 為具有自然b、c的邊界解題步驟：(1)寫出積分表達(dá)式通過分部積分，可得：(2)區(qū)域剖分橫向剖分?jǐn)?shù)為9，縱向剖分?jǐn)?shù)為10，其中圓弧段剖分?jǐn)?shù)為5。利用作業(yè)三中的程序?qū)崿F(xiàn)(由于網(wǎng)格內(nèi)要畫流速矢量圖

2、，故單元編號未寫出)，另外，還需要建立本質(zhì)b.c表。(3)確定單元基函數(shù)設(shè)網(wǎng)格劃分后任意三角形單元的三個結(jié)點的坐標(biāo)值別為，函數(shù)值分別為，根據(jù)基函數(shù)的構(gòu)造思想，單元內(nèi)近似函數(shù)可表示為式：。 在單元內(nèi)作線性插值函數(shù)如下：根據(jù)基函數(shù)的插值條件，得到系數(shù)：。則基函數(shù)為：，。(4)單元分析將代入積分表達(dá)式：得單元有限元方程組為：(i=1,2,3)由(i=1,2,3)，可得：于是：，自然b.c處理：由于自然邊界條件，則。(5)總體合成單元矩陣總體矩陣；單元矩陣行號整體矩陣行號，單元矩陣列號整體矩陣列號。(6)本質(zhì)b.c處理即為了滿足本質(zhì)b.c，要對總體系數(shù)矩陣進(jìn)行處理，具體處理方法見作業(yè)二。(7)解總體方

3、程組，求出有關(guān)物理量解方程組的方法見作業(yè)一，由 及 得：三結(jié)點三角形單元，線性插值函數(shù)，每個單元只有一個流速，與單元內(nèi)坐標(biāo)無關(guān)，可理解為單元平均流速，位于單元中心(三中線交點)。源程序如下：說明：程序的部分說明作業(yè)一、二、三中已有，這里不再贅述；其中繪圖子程序在作業(yè)三中也已有，這里略去。program yzrlimplicit noneinterfacesubroutine linear_equation_bc(n,a1,b,x_result)integer:i,j,k,imax integer,intent(in):nreal:max,creal,dimension(:,:),intent(

4、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 !定義單元、節(jié)點編號integer:m,n,m0integer,dimension(:,:),allocatable:zt

5、!定義單元結(jié)點整體編號數(shù)組integer:a,b,c,jj,kkreal,parameter:pi=3.1415926536 !定義pi為常量，值為圓周率real, dimension(:),allocatable:x,y !定義整體結(jié)點坐標(biāo)數(shù)組real:x1,y1,r !定義網(wǎng)格劃分區(qū)域real:bcx1,bcx2,bcy1,bcx,bcy,bcyh !定義x，y方向及圓弧段計算步長real:x0,y0 !定義原點坐標(biāo)real:ux !定義來流速度integer:z !定義本質(zhì)b.c點個數(shù)integer, dimension(:),allocatable:bcjd !定義本質(zhì)結(jié)點編號real

6、, dimension(:),allocatable:bcz !定義本質(zhì)b.c值real,dimension(3,3):dya !定義單元系數(shù)矩陣real,dimension(:,:),allocatable:zta !整體系數(shù)矩陣real,dimension(:,:),allocatable:bb,cc !定義基函數(shù)系數(shù)矩陣real:dym,d !定義單元三角形面積real, dimension(:),allocatable:cs !定義常數(shù)項數(shù)組real, dimension(:),allocatable:f !定義結(jié)點流函數(shù)值real, dimension(:),allocatable:

7、vx,vy,v !定義單元x、y方向流速及合成流速real, dimension(:),allocatable:zx,zy,zxx,zyy !定義單元中心及流速終點坐標(biāo)real, dimension(:),allocatable:jx !定義合成流速與x方向夾角real, dimension(:),allocatable:jtx1,jty1,jtx2,jty2 !箭頭坐標(biāo)print*,'請輸入x方向等分?jǐn)?shù)m，y方向等分?jǐn)?shù)n，圓弧等分?jǐn)?shù)m0'read*,m,n,m0print*,'請輸入網(wǎng)格劃分部分尺寸x1，y1，r'read*,x1,y1,rdyz=m*n *2

8、 !計算單元總數(shù)jdz=(m+1)*(n+1) !計算結(jié)點總數(shù)allocate(x(jdz),y(jdz),zt(dyz,3)do i=1,m !計算局部節(jié)點編號與整體節(jié)點編號的關(guān)系數(shù)組do j=1,2*ndy=(i-1)*2*n+jif(mod(j,2)=0)then !計算偶數(shù)單元的節(jié)點對應(yīng)關(guān)系數(shù)組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 !計算奇數(shù)單元的節(jié)點對應(yīng)關(guān)系數(shù)組zt(dy,1)=j/2+(n+1)*(i-1)+1zt(dy,2)=j/2+(n+1)*

9、(i-1)+2zt(dy,3)=j/2+n+(n+1)*(i-1)+2end ifend doend doopen(1,file='單元中整體結(jié)點編號.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)=

10、(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 !計算與圓弧段對應(yīng)的x1邊上的等分點的x坐標(biāo)y(jd1)=y1 !計算與圓弧段對應(yīng)的x1邊上的等分點的y坐標(biāo)x(jd2)=x1-r*cos(i-(m-m0+1)*bcyh) !計算圓弧段m0等分點的各個x坐標(biāo)y(jd2)=r*sin(i-(m-m0+1)*bcyh) !計算圓弧段m0等

11、分點的各個y坐標(biāo)bcx=(x(jd2)-x(jd1)/n !計算該計算區(qū)域的x方向的步長bcy=(y(jd2)-y(jd1)/n !計算該計算區(qū)域的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='整體結(jié)點坐標(biāo).txt') !將整體結(jié)點坐標(biāo)保存在txt文檔中write(2,'(a,i3,a,f8.4,3x,f8.4)'),('第',i,'點坐標(biāo)為:',x(i),y(i),i=1,jdz

12、)close(2)print*,'請輸入圖形名稱'read(*,'(a)')namelength=len_trim(name) !用內(nèi)部函數(shù)求出name去掉尾部空格后的長度file=name(:length)/'.dxf' !連接字符串"name"與".dxf"open(3,file=file)call staplot !調(diào)用繪制dxf文件的開頭子程序ly='wg'length=len_trim(ly)x0=0 !給定繪圖原點y0=0do i=1,dyza=zt(i,1)b=zt(i,2)c

13、=zt(i,3)call line(x(a),y(a),x(b),y(b),ly(:length) !調(diào)用繪圖子程序畫線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 !計算具有本質(zhì)b.c節(jié)點個數(shù)allocate(bcjd(z),bcz(z)do

14、i=1,z !定義本質(zhì)b.c,找到具有本質(zhì)b.c的節(jié)點和其對應(yīng)的值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='本質(zhì)b.c表.txt') !建立本質(zhì)b.c表write(4,'(a,3x,a,3x,a)'),'序號','對應(yīng)序號整體編號','

15、本質(zhì)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 !確定單元基函數(shù)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，編

16、號均為逆時針，故面積均為正bb(i,1)=(y(zt(i,2)-y(zt(i,3)/d !計算基函數(shù)系數(shù)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!計算各單元系數(shù)矩陣，由于自然b.c為0，故fi(e)為0open(5,file='單元系數(shù)矩陣.txt')write(5,'(a,i3)'),'單元

17、9;,ido j=1,3 !總體系數(shù)矩陣的合成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) !輸出單元系數(shù)矩陣end doend doclose(5)open(6,file='整體系數(shù)矩陣.txt')write(6,*),(zta(i,j),j=1,jdz),i=1,jdz) !輸出整體系數(shù)矩陣close(6)allocate(cs(jdz),f(jd

18、z) cs=0do k=1,z ! 處理本質(zhì)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='本質(zhì)b.c處理后的整體系數(shù)矩陣.txt')write(7,*),(zta(k,j),j=1,jdz),k=1,jdz) !輸出本質(zhì)b.c處理后的整體系數(shù)矩陣close(7)open(8,file='方程組右側(cè)矢量項.txt')do k=1,jdz !輸出方程組右側(cè)矢量項write(8,'(a,i3,a

19、,f6.3)'),'第',k,'行矢量項',cs(k)end doclose(8)call linear_equation_bc(jdz,zta,cs,f) !調(diào)用子程序解線性方程組open(9,file='結(jié)點流函數(shù)值.txt')do k=1,jdzwrite(9,'(a,i3,a,f10.6)'),'第',k,'點函數(shù)值',f(k)end doclose(9)allocate(vx(dyz),vy(dyz),v(dyz)open(10,file='單元x、y向流速.txt'

20、;)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(

21、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)

22、=(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坐標(biāo)zyy(i)=zy(i)+sin(jx(i)*v(i)*0.1 !流線終點y坐標(biāo)call line(zx(i),zy(i),zxx(i),zyy(i),ly(:le

23、ngth) !畫流線jtx1(i)=zxx(i)-sin(-jx(i)+8*pi/18.0)*v(i)*0.03 !箭頭終點x坐標(biāo)jty1(i)=zyy(i)-cos(-jx(i)+8*pi/18.0)*v(i)*0.03 !箭頭終點y坐標(biāo)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坐標(biāo)jty2(i)=zyy(i)-sin(jx(i)-pi/18.0)*v(i)*0.03 !箭頭終點y坐標(biāo)call 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 none