微分方程生成网格技术以及程序.doc

《微分方程生成网格技术以及程序.doc》由会员分享，可在线阅读，更多相关《微分方程生成网格技术以及程序.doc（43页珍藏版）》请在三一文库上搜索。

1、喷管网格生成报告一、 PDE方法建立坐标转换假定物理平面为，计算平面为，物理平面和计算平面之间的Jacob矩阵 为，为了满足两平面之间一一对应的关系，应该是非平凡的。利用微分计算的链式法则 ，可以得到 将其代入Laplace方程中，并做和运算，容易得到如下逆变换方程， 其中， ， ，此为非线性椭圆形方程。二、 逆变换PDE的差分方程 根据差分知识，对于椭圆形方程，其离散方法宜采用中心型差分格式，因此差分方程是隐式的。x，y的各阶导数的中心差分格式为，将上述式子代入到逆变换方程的系数中，代入逆变换方程中，并令得逆变换PDE的差分方程， 三、 边界条件1. 几何边界收敛圆弧： 收敛錐： 喉部圆弧：

2、 扩散錐： 直线入口边界： 直线出口边界： 圆弧入口边界： 圆弧出口边界： 2. 四边界节点边值的确定(1) 直线进出口边界时，喷管壁及喷管轴上的逆变换为： 其中为喷管壁入口处的轴向坐标，为喷管壁出口处的轴向坐标，A为调节参数。进出口边界等距分割，即等参变换： (2) 圆弧进出口边界时，以喉部为界，设喷管壁在喉部上游最临近喉部的节点轴向坐标为，圆弧入口、出口边界的轴向坐标为和，则喉部上游的喷管轴上的逆转换为：喉部下游的喷管轴上的逆转换为：喷管壁上的逆变换为：左右进出口边界等弧分割：进口圆弧： 出口圆弧：四、计算方法1. 简单迭代法（Jacob法）2. Gauss-Seidel迭代法迭代时，列按

3、从左至右，行按从下到上的顺序。当计算时，总是启用前面最新计算出的，计算时也如此，3. 逐次松弛迭代（SOR法）SOR迭代格式是第k+1次迭代的结果看成是第k次迭代的结果加上一个校正值乘上一个松弛因子，使得改进后的迭代方法收敛速度较快。这里将Gauss-Seidel迭代值作为，再将与加权平均，得这里采用超松弛，取4. 线（块）迭代法每次用直接法（解三对角矩阵方程）求解一行或一列未 知数(1) Jacob法逆变换PDE差分离散方程经过移项可化为三对角矩阵方程如下（这里是一行一行的求解）：即 用TDMA（追赶法）求解上面三对角矩阵方程(2) Gauss-Seidel线迭代即 (3) SOR法将Gau

4、ss-Seidel线迭代值作为，再将与加权平均，得5. 交替方向隐式迭代（ADI法SOR法）第一步横扫然后采用的超松弛第二步竖扫然后采用的超松弛6. 迭代初值选择采用轴对称形式五、程序框图程序结构见下页，原程序见附件的F90文件。参数A，ibdy，iopt赋值 直线进出口边界时四边界节点边值计算ibdy=1 .TRUE.圆弧进出口边界时四边界节点边值计算ibdy=2 .FALSE. .TRUE. .FALSE.计算迭代初值iopt=1点迭代Jocobi法计算 .TRUE. .FALSE.iopt=2点迭代G-S法计算 .TRUE. .FALSE.ADI法计算iopt=7.TRUE. .FALS

5、E.写网格文件和输出信息文件END六、计算结果及结果分析讨论本次结果分析取Imax=41，Jmax=11，收敛阀值为0.0001。 E-005 E-005迭代次数占用CPU时间（s）E-001A0.5直线进出口边界点迭代Jocobi法3.3704087618579.9511358318428623.90625点迭代G-S法9.81609022832113.169459033864462.18750点迭代SOR法9.9900488815302.9920590484781800.93750线迭代Jocobi法9.8165017377972.5591313595802111.25000线迭代G-S法

6、9.9689350027911.5054647974912041.09375线迭代SOR法2.7994429318725.898914990165310.15625ADI法7.3489510961084.798758796509290.31250圆弧进出口边界点迭代Jocobi法2.9348747438929.9098331366237923.75000点迭代G-S法9.93864719421311.574220225354141.56250点迭代SOR法9.9391854221052.7491749819301650. 78125线迭代Jocobi法9.6628048530885.65421

7、73883842181.25000线迭代G-S法9.7427598296391.9283155215712121.09375线迭代SOR法4.9806108101957.700308198321300. 15625ADI法8.5634297638667.851372482292290. 46875A1.0直线进出口边界点迭代Jocobi法2.9985585890349.94784834427110065.00000点迭代G-S法9.8387874997579.7588501422945352.81250点迭代SOR法9.6710598854492.5123367277662150. 93750

8、线迭代Jocobi法9.5533554941652.8407889049122111.25000线迭代G-S法9.5116746823021.6691808693902031.09375线迭代SOR法7.9221099926645.560956054617320. 15625ADI法8.3215634050145.871829348924280.31250圆弧进出口边界点迭代Jocobi法2.8150440212899.9381689703919344.53125点迭代G-S法9.9650917774258.7940420395945012.34375点迭代SOR法9.444402555658

9、2.2810253874622001.09375线迭代Jocobi法9.5781496081726.0734195198802171.09375线迭代G-S法9.8369558812152.3411005322322101.25000线迭代SOR法7.4194045472345.418009688717320.15625ADI法6.2482344587477.339576387899290. 31250A1.5直线进出口边界点迭代Jocobi法2.9716871508239.93753876961110855.15625点迭代G-S法9.9837604241068.57229835840058

10、42.65625点迭代SOR法9.5254844870492.2669275827612351.25000线迭代Jocobi法9.9582058823983.1757688389692101.09375线迭代G-S法9.4885995980591.8134785935152020.93750线迭代SOR法8.7260806807594.673832856338330.15625ADI法9.7129710142587.487822192375270.31250圆弧进出口边界点迭代Jocobi法2.8379017942309.92680483373710134.84375点迭代G-S法9.8733

11、601750647.5870519882235512.50000点迭代SOR法9.8848285471382.1988384559362191.09375线迭代Jocobi法9.7759061226226.3702511132302161.25000线迭代G-S法9.8951507180392.5828436038372091.25000线迭代SOR法8.5320183584814.793943165193330.15625ADI法4.6795299315726.765606103886290.15625表1 计算结果比较在A0.5，1.0，1.5时，直线进出口边界和圆弧进出口边界情况下算出的

12、物理平面上网格如图1、图2所示，可看出A值越大，喷管喉部的网格线越密，出口入口处网格越稀疏。从表一可以看出：a 所以相对于Jocobi迭代法，G-S迭代法收敛速度加快一倍左右，占用CPU时间也少一半左右；b SOR法是G-S迭代格式的一种加速方法，比G-S迭代法收敛速度至少加快一倍，占用CPU时间也减少；c 线迭代每次用直接法求解一行或一列未知数，即解三对角矩阵方程，相当于显式求解，而点迭代是联立求解所有未知数，相当于隐式求解。所以线迭代比点迭代收敛速度要快，占用CPU时间要少；d ADI法是七种方法中所需迭代次数最少的，占用CPU时间也非常低；e 从整体效果看线迭代SOR法和ADI法方法均比

13、较优秀，但就编程的难易程度而言，采用线迭代SOR法是很好的选择。 图1 直边入口出口网格（Jocobi点迭代法） 图2 圆弧入口出口网格（Jocobi点迭代法）七、建议和要求程序：c*c * 喷管物理平面和计算平面之间的坐标转换 *c *入口和出口边界均为垂直于轴向X轴的直线 *c * 分别用点迭代，线迭代和ADI迭代进行计算 *c * program main integer I_M,J_Mparameter(A=1.5,I_M=41,J_M=11)integer i,j,opt,opt1,opt2c 计算平面上点的坐标值 integer x_c(1:I_M,1:J_M),y_c(1:I_M

14、,1:J_M)c 物理平面上点的坐标值 real x(1:I_M,1:J_M),y(1:I_M,1:J_M)c 给计算平面上的点赋初值do 10,i=1,I_M do 20,j=1,J_M x_c(i,j)=i-1 y_c(i,j)=j-1 20 continue 10continuec 计算物理平面上喷管壁及喷管轴的坐标值call bond_cal(x_c,x,y)c 给中间坐标赋迭代初值 call initial(x,y)c 开始进行迭代求解 write(*,*) 请选择您要的迭代类型:write(*,*) 点迭代请输入1，线迭代请输入2，ADI迭代请输入3read *,optif (op

15、t.eq.1) then write(*,*) 请选择您要的点迭代类型: write(*,*)JACOBI迭代请输入1，GUASS迭代请输入2，松弛迭代请输入3 read *,opt1 if (opt1.eq.1) then call dian_j(x,y) write(*,*) dd else if(opt1.eq.2) then call dian_g(x,y) else call dian_s(x,y) endif else if(opt.eq.2) then write(*,*) 请选择您要的线迭代类型: write(*,*)JACOBI迭代请输入1，GUASS迭代请输入2，松弛迭代请

16、输入3 read *,opt2 if (opt2.eq.1) then call xian_j(x,y) else if (opt2.eq.2) then call xian_g(x,y) else call xian_s(x,y) endif else call adi(x,y)endifc 输出数据 open(1,file=grid.dat) write(1,*)variables= x y write(1,*) ZONE T=GRIDS F=POINT I=,i_M,J=,j_M do j=1,j_Mdo i=1,i_Mwrite(1,*)x(i,j),y(i,j)enddo enddo

17、 close(1)endc *c * 计算喷管边界的子程序 *c *subroutine bond_cal(a,b,c)integer I_M,J_M parameter(AA=1.5,I_M=41,J_M=11) integer mreal ishd,ishx,biintegera(1:I_M,1:J_M)real b(1:I_M,1:J_M),c(1:I_M,1:J_M)c 计算壁面坐标ishd=log(AA*30.5)+sqrt(1+(AA*30.5)*(AA*30.5) ishx=log(AA*(-56.86)+sqrt(1+(AA*(-56.86)*(AA*(-56.86)do 30

18、,m=1,I_M bi=a(m,J_M)*1.0/a(I_M,J_M) b(m,J_M)=1/AA*sinh(bi*(ishd-ishx)+ishx)c write(*,*) m,b(m,J_M) if (b(m,J_M).lt.-42.49) then c(m,J_M)=sqrt(20.32*20.32-(b(m,J_M)+56.86)*(b(m,J_M)+56.86) c(m,J_M)=c(m,J_M)+43.18 endif if(b(m,J_M).ge.-42.49).and.(b(m,J_M).lt.-8.98) then c(m,J_M)=(-1)*b(m,J_M)+15.06 e

19、ndif if(b(m,J_M).ge.-8.98).and.(b(m,J_M).lt.3.287) then c(m,J_M)=33.02-sqrt(12.7*12.7-b(m,J_M)*b(m,J_M) endif if(b(m,J_M).ge.3.287).and.(b(m,J_M).le.30.5) then c(m,J_M)=0.268*b(m,J_M)+19.867 endifc 计算轴上坐标 b(m,1)=b(m,J_M) c(m,1)=0 30 continuec 计算出入口坐标 do 40,n=1,J_M b(1,n)=-56.86 c(1,n)=6.35*(n-1) b(I

20、_M,n)=30.5 c(I_M,n)=2.804*(n-1) 40 continue endc *c * 给中间坐标赋迭代初值 *c *subroutine initial(b,c)integer I_M,J_M parameter(AA=1.5,I_M=41,J_M=11) real b(1:I_M,1:J_M),c(1:I_M,1:J_M)real xs,xdinteger m,n do 70,m=2,I_M-1 do 80,n=2,J_M-1 xs=(n-1)*1.0/(J_M-1) xd=(J_M-n)*1.0/(J_M-1) b(m,n)=xd*b(m,1)+xs*b(m,J_M)

21、 c(m,n)=xd*c(m,1)+xs*c(m,J_M) 80 continue 70continue endc *c * 点迭代JACOBI求解 *c *subroutine dian_j(b,c)integer I_M,J_Mreal Eparameter(E=0.001,AA=0.5,I_M=41,J_M=11)real b(1:I_M,1:J_M),c(1:I_M,1:J_M)real b0(1:I_M,1:J_M),c0(1:I_M,1:J_M)integer m,n,num real errx,erry,alpha,beta,gamma,xi,xj,yi,yj,temp,temp

22、1real xia,xja,yia,yja errx=1.0erry=1.0num=0 do 90,while(errx.gt.E).or.(erry.gt.E) num=num+1 errx=0 erry=0 do 110,m=1,I_M do 120,n=1,J_M b0(m,n)=b(m,n) c0(m,n)=c(m,n) 120continue 110 continue do 130,m=2,I_M-1 do 140,n=2,J_M-1 xi=b0(m+1,n)-b0(m-1,n) xj=b0(m,n+1)-b0(m,n-1) yi=c0(m+1,n)-c0(m-1,n) yj=c0(

23、m,n+1)-c0(m,n-1) alpha=(xj*xj+yj*yj)/4.0 gamma=(xi*xi+yi*yi)/4.0 beta=(xi*xj+yi*yj)/4.0 xia=b0(m+1,n)+b0(m-1,n) xja=b0(m,n+1)+b0(m,n-1) yia=c0(m+1,n)+c0(m-1,n) yja=c0(m,n+1)+c0(m,n-1) temp=b0(m+1,n+1)-b0(m+1,n-1)+b0(m-1,n-1)-b0(m-1,n+1) temp1=c0(m+1,n+1)-c0(m+1,n-1)+c0(m-1,n-1)-c0(m-1,n+1) b(m,n)=(a

24、lpha*xia-0.5*beta*temp+gamma*xja)/(2.0*(alpha+gamma) c(m,n)=(alpha*yia-0.5*beta*temp1+gamma*yja)/(2.0*(alpha+gamma) 140 continue 130 continue do 150,m=2,I_M-1 do 160,n=2,J_M-1 errx=errx+abs(b(m,n)-b0(m,n)erry=erry+abs(c(m,n)-c0(m,n) 160 continue 150 continue 90 continuewrite(*,*) JACOBI点迭代次数为：,num e

25、ndc *c * 点迭代GUASS求解 *c * subroutine dian_g(b,c)integer I_M,J_Mreal Eparameter(E=0.001,AA=0.5,I_M=41,J_M=11)real b(1:I_M,1:J_M),c(1:I_M,1:J_M)real b0(1:I_M,1:J_M),c0(1:I_M,1:J_M)integer m,n,numreal errx,erry,alpha,beta,gamma,xi,xj,yi,yj,temp,temp1real xia,xja,yia,yjaerrx=1.0erry=1.0num=0do 170,while(

26、errx.gt.E).or.(erry.gt.E) num=num+1 errx=0 erry=0 do 180,m=1,I_M do 190,n=1,J_M b0(m,n)=b(m,n) c0(m,n)=c(m,n) 190 continue 180 continue do 200,m=2,I_M-1 do 210,n=2,J_M-1 xi=b0(m+1,n)-b(m-1,n) xj=b0(m,n+1)-b(m,n-1) yi=c0(m+1,n)-c(m-1,n) yj=c0(m,n+1)-c(m,n-1) alpha=(xj*xj+yj*yj)/4.0 gamma=(xi*xi+yi*yi

27、)/4.0 beta=(xi*xj+yi*yj)/4.0 xia=b0(m+1,n)+b(m-1,n) xja=b0(m,n+1)+b(m,n-1) yia=c0(m+1,n)+c(m-1,n) yja=c0(m,n+1)+c(m,n-1) temp=b0(m+1,n+1)-b0(m+1,n-1)+b(m-1,n-1)-b(m-1,n+1) temp1=c0(m+1,n+1)-c0(m+1,n-1)+c(m-1,n-1)-c(m-1,n+1) b(m,n)=(alpha*xia-0.5*beta*temp+gamma*xja)/(2.0*(alpha+gamma) c(m,n)=(alpha*

28、yia-0.5*beta*temp1+gamma*yja)/(2.0*(alpha+gamma) 210 continue 200 continue do 220,m=2,I_M-1 do 230,n=2,J_M-1 errx=errx+abs(b(m,n)-b0(m,n) erry=erry+abs(c(m,n)-c0(m,n) 230 continue 220 continue 170 continue write(*,*) GAUSS-SEIDEL点迭代次数为：,num endc *c * 点迭代松弛求解 *c * subroutine dian_s(b,c) integer I_M,J

29、_M real E,PIparameter(PI=3.14159,E=0.001,AA=0.5,I_M=41,J_M=11)real b(1:I_M,1:J_M),c(1:I_M,1:J_M)real b0(1:I_M,1:J_M),c0(1:I_M,1:J_M)integer m,n,numreal errx,erry,alpha,beta,gamma,xi,xj,yi,yj,temp,temp1real xia,xja,yia,yja,werrx=1.0erry=1.0num=0do 240,while(errx.gt.E).or.(erry.gt.E) num=num+1 errx=0

30、erry=0 do 250,m=1,I_M do 260,n=1,J_M b0(m,n)=b(m,n) c0(m,n)=c(m,n) 260 continue 250 continue do 270,m=2,I_M-1 do 280,n=2,J_M-1 xi=b0(m+1,n)-b(m-1,n) xj=b0(m,n+1)-b(m,n-1) yi=c0(m+1,n)-c(m-1,n) yj=c0(m,n+1)-c(m,n-1) alpha=(xj*xj+yj*yj)/4.0 gamma=(xi*xi+yi*yi)/4.0 beta=(xi*xj+yi*yj)/4.0 xia=b0(m+1,n)+

31、b(m-1,n) xja=b0(m,n+1)+b(m,n-1) yia=c0(m+1,n)+c(m-1,n) yja=c0(m,n+1)+c(m,n-1) temp=b0(m+1,n+1)-b0(m+1,n-1)+b(m-1,n-1)-b(m-1,n+1) temp1=c0(m+1,n+1)-c0(m+1,n-1)+c(m-1,n-1)-c(m-1,n+1) b(m,n)=(alpha*xia-0.5*beta*temp+gamma*xja)/(2.0*(alpha+gamma) c(m,n)=(alpha*yia-0.5*beta*temp1+gamma*yja)/(2.0*(alpha+g

32、amma)c no=(alpha*cos(PI/(I_M-1)+gamma*cos(PI/(J_M-1)/(alpha+gamma)c w=2.0/(1+sqrt(1-no*no) w=1.4999 b(m,n)=w*b(m,n)+(1-w)*b0(m,n) c(m,n)=w*c(m,n)+(1-w)*c0(m,n) 280 continue 270 continue do 290,m=2,I_M-1 do 300,n=2,J_M-1 errx=errx+abs(b(m,n)-b0(m,n) erry=erry+abs(c(m,n)-c0(m,n) 300 continue 290 conti

33、nue 240 continue write(*,*) 松弛点迭代次数为：,numendc *c * 线迭代JACOBI求解 *c *subroutine xian_j(b,c) integer I_M,J_M real Eparameter(E=0.001,AA=0.5,I_M=41,J_M=11) real b(1:I_M,1:J_M),c(1:I_M,1:J_M) real b0(1:I_M,1:J_M),c0(1:I_M,1:J_M)integer m,n,numreal oo(2:I_M-1),pp(2:I_M-1),qq(2:I_M-1),ffx(2:I_M-1),ffy(2:I_

34、M-1)real errx,erry,alpha,beta,gamma,xi,xj,yi,yj,temp,temp1 real xia,xja,yia,yjaerrx=1.0erry=1.0num=0do 330,while(errx.gt.E).or.(erry.gt.E) num=num+1 errx=0 erry=0 do 340,m=1,I_M do 350,n=1,J_M b0(m,n)=b(m,n) c0(m,n)=c(m,n) 350 continue 340 continue do 360,n=2,J_M-1 do 370,m=2,I_M-1 xi=b0(m+1,n)-b0(m

35、-1,n) xj=b0(m,n+1)-b0(m,n-1) yi=c0(m+1,n)-c0(m-1,n) yj=c0(m,n+1)-c0(m,n-1) alpha=(xj*xj+yj*yj)/4.0 gamma=(xi*xi+yi*yi)/4.0 beta=(xi*xj+yi*yj)/4.0 xia=b0(m+1,n)+b0(m-1,n) xja=b0(m,n+1)+b0(m,n-1) yia=c0(m+1,n)+c0(m-1,n) yja=c0(m,n+1)+c0(m,n-1) temp=b0(m+1,n+1)-b0(m+1,n-1)+b0(m-1,n-1)-b0(m-1,n+1) temp1

36、=c0(m+1,n+1)-c0(m+1,n-1)+c0(m-1,n-1)-c0(m-1,n+1) if(m.eq.2) then oo(m)=0 else oo(m)=(-1)*alpha endif qq(m)=(-1)*alpha pp(m)=2*(alpha+gamma) if(m.eq.2) then ffx(m)=(-0.5)*beta*temp+gamma*xja+alpha*b0(1,n) ffy(m)=(-0.5)*beta*temp1+gamma*yja+alpha*c0(1,n) else if(m.eq.I_M-1) then ffx(m)=(-0.5)*beta*temp+gamma*xja-qq(m)*b0(I_M,n) ffy(m)=(-0.5)*beta*temp1+gamma*yja-qq(m)*c