北叶轮机高等气动力学大作业

作业条件和要求设计参数进气总压 =101325.0 Pa进气总温 =288.15 K质量流量 G=550 kg/s总增压比 =1.5绝热效率 =0.90气体常数 R= 287.06 J/kg/K定压比热 = 1004.7 J/kg/K比热比 k = 1.4确定（说明理由）转速，流道几何，各排叶片参数的展向分布（最少根、中、尖三个截面）：叶片尾缘半径，D因子，压比，效率，轮缘功，子午速度，相对速度，绝对速度，相对气流角，绝对气流角，静压，静温，总压，总温，密度，各叶片的基元几何（最少根、中、尖三个截面） 叶轮机高等气动力学大作业——风扇设计院（系）名称 专业名称学号 学生姓名一、控制方程（1）连续性方程 积分形式的连续性方程（对展向计算站而言）（2）运动方程以焓熵形式描述的展向平衡方程（3）能量方程（4）状态方程（5）熵增关系式（6）沿流线斜率曲率二、数值过程（1）离散方式（i为计算站标号，j为流线标号）（2）运动方程流线平均参数 （3）连续方程三、求解流场流程图四、压气机设计（1）压气机级数的确定：压气机设计流量为550kg/s，压比1.5，可以看出属于民用风扇范畴，因此初步选定为单级轴流压气机。

一级转子与一级静子的形式2）压气机流道形式的选择：对于级数较少的压气机而言，流道采用等外径设计规律，可以充分利用压气机较高的叶尖切线速度，实现较大的加功量因此本压气机选择等外径气流通道设计3）压气机轮毂比的确定：转子轮毂比的选择需同时考虑气动性能和结构强度等方面的因素．较小的进口轮毂比对气动性能有好处，可以提高风扇的效率，但不利于转子的结构强度，且过小的轮毂比会导致根部的加功量足对于发动机进口机风扇而言，其轮毂比多在0.32-0.4之间，因此选定压气机轮毂比为0.35并定义转子进口机匣与轮毂半径分别为,，且满足4）压气机进口参数的确定： 压气机进口轴向马赫数大小与压气机性能有密切的关系，当进口轴向马赫数增加，压气机效率和裕度都要下降根据以往经验，高压压气机进口马赫数通常为0.48-0.52之间，对于风扇，进口马赫数通常在0.6-0.68之间因此本设计选定进口马赫数Ma=0.6，并认为进口流场为均匀场，来流方向轴向，无预旋进口静温：由，可求得= 268.80 K当地音速：c==328.67m/s，m/s进口静压：由 计算得79439.20Pa进口空气密度：1.030空气质量流量：G=550 kg/s由连续方程，风扇进口面积 由轮毂比和面积公式，可得=0.991m，=0.347m。

5）压气机出口参数的确定：风扇出口Ma数范围为0.48-0.5，选定出口Ma为0.48压气机级出口总压101325*1.5=151987.50Pa，由压气机绝热效率，可求得压气机出口总温为327.47K313.04K，129816.75Pa，1.45，由压气机出口截面积以及出口外径，可得压气机出口轮毂半径为6）转子出口截面参数确定：确定静子总压恢复系数，转子出口总压，转子压比=1.531，由转子效率公式可得转子效率0.948，因此由轮缘功公式可得=39518.53J/kg由单级风扇功分配的特点，应适当减小根部和尖部的加功量：对于叶根处，加功量应尽可能小，以保证根部气流转折角不至于过大；对于尖部，同样应减小加功量，来减小尖部气流的分离由轮缘功的公式可知，可通过给定压比和效率沿展向的分布来保证不同展高位置的功分配表1给出了不同截面处，转子叶片轮缘功压比和效率的分配情况，图1给出了不同截面处压比分配曲线表1 转子叶片不同展高位置参数分布百分比展高压比效率轮缘功J/kg01.40.932458.35%1.4250.934253.610%1.450.936026.520%1.50.939507.830%1.540.942233.750%1.59310.945775.170%1.6050.946557.290%1.570.91543519.71001.450.9235243.3平均值1.5430.905839508.4 图1转子压比沿展向分布 图2 转子出口环量展向分布（7）叶尖切线速度和转速的确定：根据轮缘功公式：可知，切线速度越大轮缘功越大，有利于压比的提高，当轮缘功不变时，切线速度越大，扭速越小，即叶片的弯度越小。

当然，叶尖切线速度过大时，尖部相对马赫数过高，激波损失增强，从而导致效率的下降，这对于现代压气机设计来说已经成为了不可避免的问题对于相同的角速度，为得到更大的加功量，根部需要更高的扭速，即叶片的弯角更大由公式：知，切线速度与转速密切关联，本压气机选取叶尖切线速度为380m/s，其中r为转子外径，所以转速n=3661.69rpm，角速度，所以对应的根中尖三个截面处的负荷系数分别为2.23、0.601和0.274由压气机出口气流速度，所以对应根中尖的流量系数分别为：1.280,0.864,0.448图2给出了本压气机尖部流量系数负荷系数的分配在Smith图上的位置（图中）可以看到，本压气机中部和尖部的流量系数与负荷系数的分配比较合理（坐标分别为0.864,1.202；0.448,0.548 ），且在图中效率较高位置，然而根部的负荷过大，这有可能导致根部分离较严重，效率不高图3 流量系数与负荷系数对应关系（8）叶片弦长及稠度的确定：叶片展弦比与发动机使用成本、叶尖速度、压气机效率、裕度有直接关系在压气机设计时首先选定压气机的平均展弦比低展弦比对气动性能的影响主要表现在效率和抗失速能力上，低的展弦比就意味着有更高的弦长雷诺数，子午面内激波更斜，这些因素都有利于提高压气机效率同时有利于抗失速能力；另外，宽弦叶片可以省掉阻尼台，这毫无疑问有利于效率的提升。

宽弦叶片边界层较厚这不利于压气机效率的提高综上所述：由于本压气机根部区域负荷过高，气流在根部区域有较大的转折角，因此需要对根部采用较小的展弦比来获得更大的弦长，从而防止根部叶型弯角过大而对于尖部区域则采用较大的展弦比来获得较小的弦长，从而提高尖部效率因此本压气机转子平均展弦比为2，根部展弦比为1.8，中部为2.0，尖部为2.5由转子进口前缘叶高L=0.684m可得转子根中尖弦长分别为0.38m,0.342m,0.274m同理，给定静子根中尖展弦比分别为2, 2, 1.8，由静子叶高L=0.4215m可得静子弦长分别为0.25,0.25,0.278稠度的确定：稠度的选取与气流的转角密切相关，同时也反映了叶片间的通道面积，稠度过大会导致通道面积减小，容易发生堵塞根据D因子表达式：，可知稠度的取值能够影响到D因子的大小因此对转子根中尖稠度分别取2.3,1.7,1.2；静子根中尖稠度分别取1.9,1.3,1.2根据D因子表达式得转子根中尖D因子分别为0.448,0.576,0.279；静子根中尖D因子分别为0.363,0.242,0.021根据叶片数公式可求得转子叶片数为23片，静子为24片表2所示为根据以上计算所得压气机各相关参数。

表2 压气机各相关参数动叶根中尖展弦比1.822.5弦长b(m)0.3800.3420.274稠度τ=b/t2.31.71.2栅距t(m)0.1650.2010.228D因子0.4490.5760.279轴向长度Δz(m)0.3310.332叶片数23进口速度V1（m/s）166.05166.05166.05进口总压P1*（Pa）101325101325101325进口静压P1（Pa)87538.9487538.9487538.94进口总温T1*(K)288.15288.15288.15进口静温T1(K)274.43274.43274.43出口总压P2*(Pa)140841.75162120146921.25出口静压P2(Pa)105120.95136899.32131791.70出口总温T2*(K)319.74334.16324.01出口静温T2(K)284.18309.06308.95出口轴向速度V2a（m/s)150.13150.13150.13压比π*1.391.61.45效率η*0.900.900.90静叶展弦比221.8弦长b(m)0.2510.2510.278稠度τ=b/t1.91.31.2栅距t(m)0.1320.1930.232D因子0.3630.2420.021轴向长度Δz(m)0.2430.303叶片数24出口速度V3（m/s）170.25170.25170.25出口总压P3*（Pa)151987.5151987.5151987.5出口静压P3（Pa)137494.26137494.26137494.26出口总温T3*(K)327.47327.47327.47出口静温T3(K)313.05313.05313.05（9）计算站给定及流道几何设计：根据以上数据，初步设计压气机流道几何如图3所示。

并在转子叶片之前给定9个计算站，转静之间给定两个计算站，静子后给定4个计算站，没有给定叶片内部计算站图4 压气机流道图五、通流计算结果转子出口截面各参数展向分布如表3所示；静子出口截面个参数展向分布如表4所示图5、图6分别给出了通流计算全流场的相对马赫数等值线和静压等值线分布表3 转子出口截面各参数沿展向分布半径mm相对展高m/s子午速度m/s绝对速度m/s相对速度m/s压比效率密度D因子总温K静温K总压Pa静压Pa534.18 0.00%155.37 220.06163.96 1.40 0.90 1.27 0.43 320.44 296.34141834.2107877.2608.82 14.41%151.66 208.88178.241.42 0.90 1.31 0.42 322.07 300.36144143.5112900.8674.79 27.15%149.45 201.82196.35 1.44 0.90 1.34 0.41 323.65 303.38146405.5116748.3734.10 38.61%148.08 200.46211.491.49 0.90 1.37 0.41 326.63 306.63150735120830.7788.35 49.08%147.22 198.21 228.341.52 0.90 1.40 0.40 328.73 309.18153859.9124142.1838.59 58.79%146.62 195.91245.501.54 0.90 1.42 0.39 330.42 311.31156383.2126964885.38 67.82%146.29 195.78260.011.58 0.90 1.45 0.39 332.84 313.76160069.4130195.8929.45 76.33%146.11 193.62276.891.60 0.90 1.46 0.38 333.96 315.30161794.7132306.2971.34 84.42%145.94 190.58 294.681.60 0.90 1.47 0.36 334.33 316.26162377.81336721011.50 92.18%145.70 185.24 315.511.59 0.91 1.48 0.34 333.03 315.95160923.8133848.31052.00 100.00%144.70 168.31 354.881.45 0.92 1.42 0.26 323.23 309.13146919.3125686.7表4 静子出口各参数沿展向分布半径mm相对展高子午速度m/s绝对速度m/s相对速度m/s密度D因子总温K静温K总压Pa静压Pa601.970.00%155.70155.70155.701.390.03320.44308.38140415.9122767662.8213.52%155.32155.32155.321.400.03322.07310.07142702.1124935.4718.3025.85%155.14155.14155.141.420.02323.65311.67144941.4127019.5769.0137.12%155.18155.18155.181.450.03326.63314.64149227.7130927.6816.0347.56%155.34155.34155.341.470.03328.73316.72152321.3133718.1860.0657.35%155.50155.50155.501.490.02330.42318.38154819.4135967.2901.3466.52%155.73155.73155.731.510.03332.84320.77158468.7139250.6940.5375.23%155.97155.97155.971.520.02333.96321.85160176.8140758.4978.0883.58%156.06156.06156.061.530.01334.33322.21160754141263.41014.4991.66%155.97155.97155.971.520.00333.03320.92159314.5139948.51052.00100.00%155.01155.01155.011.430.00323.23311.27145450.1127471.3 图5 相对马赫数等值线 图6 静压分布等值线六、叶片造型转子采用圆弧中弧线，由于转子叶片跨音，所以对尖部以下区域采用双圆弧叶型，适应来流马赫数为0.8——1.2的跨音流动，其最大厚度相对位置和最大挠度相对位置均为弦长的50%，采用重心积叠的方式。

但由于本压气机转子叶尖相对马赫数最大值超过了1.2，所以对转子叶尖采用多圆弧叶型最大厚度相对位置和最大挠度相对位置均为弦长的50%转子进出口绝对气流角和相对气流角沿展向分布如下表5所示表5 转子进出口绝对/相对气流角Y%进口绝对气流角进口相对气流角出口绝对气流角出口相对气流角0.00%0.00 -37.9145.09 -18.6249.08%0.00 -60.0342.03 -49.85100.00%0.00 -68.8830.07 -65.94叶片根部和尖部给定攻角-2°，中部给定攻角2°，由叶型弯角公式：其中 可计算得出根中尖位置处的叶型弯角，并由此画出中弧线其中，根部基元°，叶中基元°，尖部基元°计算出根部叶型安装角为26.26°，中部叶型安装角56.9°，尖部叶型安装角为65.4°，因此可画出转子根中尖三个截面的基元图7给出了转子叶片根中尖三个基元截面的叶型图7转子叶片根中尖截面叶型对静子叶片采用NACA65—010叶型，采用尾缘积叠的方式，并采用圆弧中弧线同理，可计算得静子叶片的叶型弯角，根中尖分别为°，°，°；根中尖截面的安装角分别为14.93°，11.39°，8.48°。

附：源程序 PROGRAM MAIN IMPLICIT NONE INCLUDE 'COMM.FOR' INTEGER I,J,K,N,IPASS DOUBLE PRECISION RES_WX(IMAX,JMAX),RES_MF(IMAX,JMAX),A,W,RES1,RES2,ERROR PARAMETER (ERROR=3E-3) RADPS=RPM转速*PI/30.0 CALL INPUT CALL INITIAL DO IPASS=1,IPASSMAX CALL SL_ADJ DO I=1,IMAX DO J=1,JMAX RES_WX(I,J)=ABS(WX(I,J)/WX_OLD(I,J)-1) RES_MF(I,J)=ABS(DETG(I,J)/1E6-MF*(J-1)/(JMAX-1)) ENDDO ENDDO RES1=MAXVAL(RES_MF) RES2=MAXVAL(RES_WX) WRITE(*,*) IPASS,RES1,RES2 IF((RES1.LT.ERROR).AND.(RES2.LT.ERROR)) THEN CALL CAL_STATE CALL INTP(SLD,SLDT,2) DO I=1,IMAX DO J=1,JMAX DO N=1,NOBR/2 IF((I.GE.SOLE(2*N-1)).AND.(I.LE.SOTE(2*N-1))) THEN W=RADPS ELSE W=0.0 ENDIF ENDDO WU(I,J)=VU(I,J)-W*R_CUR(I,J)/1000.0 V_REL(I,J)=SQRT(VM(I,J)**2.0+WU(I,J)**2.0) A=SQRT(KG*RG*ST(I,J)) MA(I,J)=V(I,J)/A REL_MA(I,J)=V_REL(I,J)/A ARPHA(I,J)=ATAN(VU(I,J)/VM(I,J)) BETA(I,J)=ATAN(WU(I,J)/VM(I,J)) ENDDO ENDDO DO J=1,JMAX DO N=1,NOBR DF(SOTE(N),J)=1-V_REL(SOTE(N),J)/V_REL(SOLE(N),J) $+(WU(SOTE(N),J)-WU(SOLE(N),J))/(2.0*SLDT(SOTE(N),J)*V_REL(SOLE(N),J)) ENDDO ENDDO CALL OUTPUT EXIT ELSEIF(IPASS.EQ.IPASSMAX) THEN WRITE(*,*) 'CONVERGENCE NOT REACHED...' ENDIF ENDDO END SUBROUTINE OUTPUT IMPLICIT NONE INCLUDE 'COMM.FOR' INTEGER I,J OPEN(30,FILE='FIELD.plt') WRITE(30,*) 'TITLE="FIELD"' WRITE(30,'(A300)') 'VARIABLES="CoordinateX","CoordinateR","WX","VR","VM","VU","V","WU", $"V_REL","MA","MA_REL","TT","ST","TP","SP","RHO","DF","ARPHA","BETA"' !,"J","I" WRITE(30,*) 'ZONE T="Zone',1,'"' WRITE(30,*) 'I=',IMAX,'J=',JMAX,',K=',1,',ZONETYPE=Ordered' WRITE(30,*) 'DATAPACKING=BLOCK' WRITE(30,*) WRITE(30,'(E16.9)') ((X_CUR(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((R_CUR(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((WX(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((VR(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((VM(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((VU(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((V(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((WU(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((V_REL(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((MA(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((REL_MA(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((TT(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((ST(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((TP(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((SP(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((DENS(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((DF(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') (((ARPHA(I,J)*180.0/PI),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') (((BETA(I,J)*180.0/PI),I=1,IMAX),J=1,JMAX) WRITE(30,*) CLOSE(30) END INTEGER IMAX,JMAX,NOBR,IPASSMAX,JN !NUMBER OF BLADE ROWS PARAMETER (IMAX=19,JMAX=11,NOBR=2,IPASSMAX=5E3,JN=9) DOUBLE PRECISION RPM,PI,RG,EFFI,TT0,TP0,CP,KG,MF,RADPS,TPRC,LAMDA PARAMETER (RPM=3900,PI=3.14159265358,RG=287.06) PARAMETER (TT0=288.15,TP0=101325,CP=1004.7,KG=1.4,MF=550) PARAMETER (EFFI=0.9,TPRC=0.98,LAMDA=0.52) DOUBLE PRECISION R_ORI(IMAX,JMAX),X_ORI(IMAX,JMAX),R_CUR(IMAX,JMAX), $X_CUR(IMAX,JMAX),RM(IMAX,JMAX),BR(JN,NOBR),SLD(JN,NOBR),WBLK(IMAX), $SOLE(NOBR),SOTE(NOBR),SLDT(IMAX,JMAX),YITA(IMAX,JMAX),YITA_NEW(IMAX,JMAX) !SECTION OF L.E. ENDWALL BLOCKAGE COMMON /GEOM/ R_ORI,X_ORI,R_CUR,X_CUR,RM,BR,SLD,SOLE,SOTE,WBLK,SLDT,YITA,YITA_NEW DOUBLE PRECISION DENS(IMAX,JMAX),ST(IMAX,JMAX),SP(IMAX,JMAX), $ENTRA(IMAX,JMAX),ENTHO(IMAX,JMAX),EFF_RS(JN,NOBR),EFF(IMAX,JMAX),TT(IMAX,JMAX), $TP(IMAX,JMAX),FG(IMAX,JMAX) !EFF_RS:FOR ROTER-IFFICIENCY,FOR STATOR TP RECOVERY COEFFI... COMMON /STATE/ DENS,ST,SP,ENTRA,ENTHO,EFF_RS,EFF,TT,TP,FG DOUBLE PRECISION VX(IMAX,JMAX),WX_OLD(IMAX,JMAX),VR(IMAX,JMAX),VM(IMAX,JMAX),V_REL(IMAX,JMAX), $DETG(IMAX,JMAX),VU(IMAX,JMAX),WX(IMAX,JMAX),VUR(IMAX,JMAX),VURI(JN,NOBR),V(IMAX,JMAX),WU(IMAX,JMAX) COMMON /VELOCITY/ VX,VR,VM,V_REL,VU,WX,WX_OLD,VUR,VURI,V,DETG,WU !DELTA_G DOUBLE PRECISION SIGMA(IMAX,JMAX),THETA(IMAX,JMAX),ARPHA(IMAX,JMAX),BETA(IMAX,JMAX) COMMON /ANGLE/ SIGMA,THETA,ARPHA,BETA DOUBLE PRECISION FU(IMAX,JMAX),FR(IMAX,JMAX),FX(IMAX,JMAX) COMMON /BLADE_FORCE/ FU,FR,FX DOUBLE PRECISION MA(IMAX,JMAX),REL_MA(IMAX,JMAX),DF(IMAX,JMAX) COMMON /DESIGN/ MA,REL_MA,DF SUBROUTINE INPUT IMPLICIT NONE INCLUDE 'COMM.FOR' INTEGER I,J,K DOUBLE PRECISION TEMP(JMAX) OPEN(30,FILE='BOUNDARY.DAT') DO J=1,NOBR READ(30,*) SOLE(J),SOTE(J) ENDDO READ(30,*) DO I=1,IMAX READ(30,*) X_ORI(I,1),R_ORI(I,1),X_ORI(I,JMAX),R_ORI(I,JMAX),WBLK(I) ENDDO CLOSE(30) OPEN(01,FILE='BLADE.DAT') DO I=1,NOBR READ(01,*) DO J=1,JN READ(01,*) BR(J,I),VURI(J,I),SLD(J,I),EFF_RS(J,I) ENDDO ENDDO CLOSE(01) END SUBROUTINE INITIAL IMPLICIT NONE INCLUDE 'COMM.FOR' INTEGER I,J DOUBLE PRECISION AREA(IMAX),ACR,TEMP(IMAX),A WRITE(*,*) 'INITIALLING...' A=PI*R_ORI(1,JMAX)**2.0 DO I=1,IMAX AREA(I)=PI*(R_ORI(I,JMAX)**2.0-R_ORI(I,1)**2.0) DO J=2,JMAX-1 R_ORI(I,J)=SQRT((J-1)*AREA(I)/((JMAX-1)*PI)+R_ORI(I,1)**2.0) X_ORI(I,J)=X_ORI(I,1)+(X_ORI(I,JMAX)-X_ORI(I,1))*(R_ORI(I,J)-R_ORI(I,1))/(R_ORI(I,JMAX)-R_ORI(I,1)) ENDDO DO J=1,JMAX X_CUR(I,J)=X_ORI(I,J) R_CUR(I,J)=R_ORI(I,J) ENDDO ENDDO DO I=1,IMAX DO J=1,JMAX DF(I,J)=0.0 ENDDO ENDDO OPEN(30,FILE='MESS_INI.plt') WRITE(30,*) 'TITLE="MESS_INI"' WRITE(30,'(A300)') 'VARIABLES="CoordinateX","CoordinateR"' WRITE(30,*) 'ZONE T="Zone',1,'"' WRITE(30,*) 'I=',IMAX,'J=',JMAX,',K=',1,',ZONETYPE=Ordered' WRITE(30,*) 'DATAPACKING=BLOCK' WRITE(30,*) WRITE(30,'(E16.9)') ((X_CUR(I,J),I=1,IMAX),J=1,JMAX) WRITE(30,*) WRITE(30,'(E16.9)') ((R_CUR(I,J),I=1,IMAX),J=1,JMAX) CLOSE(30) DO I=1,IMAX TEMP(I)=(R_CUR(I,JMAX)-R_CUR(I,1))/(X_CUR(I,JMAX)-X_CUR(I,1)) IF(TEMP(I).GE.0) THEN DO J=1,JMAX THETA(I,J)=ATAN(TEMP(I)) ENDDO ELSE DO J=1,JMAX THETA(I,J)=ATAN(TEMP(I))+PI ENDDO ENDIF ENDDO DO I=1,IMAX !SOLE(1) DO J=1,JMAX TT(I,J)=TT0 TP(I,J)=TP0 ENTRA(I,J)=CP*TT0 ENTHO(I,J)=0.0 ENDDO ENDDO ACR=SQRT((2.0*KG/(KG+1))*RG*TT0) DO I=1,IMAX DO J=1,JMAX WX(I,J)=LAMDA*ACR*AREA(1)/A ENDDO ENDDO DO I=1,IMAX DO J=1,JMAX VUR(I,J)=0.0 EFF(I,J)=0 ENDDO ENDDO WRITE(*,*) 'INITIALLING DONE' END SUBROUTINE INTP(QNT,QNT_NEW,FLG) !,NB,FLG) IMPLICIT NONE INCLUDE 'COMM.FOR' DOUBLE PRECISION QNT(JN,NOBR),QNT_NEW(IMAX,JMAX),QNT_TEMP(IMAX,JMAX) INTEGER I,J,K,NB,FLG !FLG=1--L.E.;FLG=2--T.E. DO I=1,NOBR DO J=1,JMAX IF(FLG.EQ.2) THEN IF(R_CUR(SOTE(I),J).LT.BR(1,I)) THEN QNT_NEW(SOTE(I),J)=(QNT(2,I)-QNT(1,I)) $*(R_CUR(SOTE(I),J)-BR(1,I))/(BR(2,I)-BR(1,I))+QNT(1,I) ELSEIF(R_CUR(SOTE(I),J).GE.BR(JN,I)) THEN QNT_NEW(SOTE(I),J)=(QNT(JN,I)-QNT(JN-1,I)) $*(R_CUR(SOTE(I),J)-BR(JN,I))/(BR(JN,I)-BR(JN-1,I))+QNT(JN,I) ELSE DO K=1,JN-1 IF((BR(K,I).LE.R_CUR(SOTE(I),J)).AND.(BR(K+1,I).GT.R_CUR(SOTE(I),J))) THEN QNT_NEW(SOTE(I),J)=(QNT(K+1,I)-QNT(K,I)) $*(R_CUR(SOTE(I),J)-BR(K,I))/(BR(K+1,I)-BR(K,I))+QNT(K,I) ENDIF ENDDO ENDIF ELSEIF(FLG.EQ.1) THEN IF(R_CUR(SOLE(I),J).LT.BR(1,I)) THEN QNT_NEW(SOLE(I),J)=(QNT(2,I)-QNT(1,I)) $*(R_CUR(SOLE(I),J)-BR(1,I))/(BR(2,I)-BR(1,I))+QNT(1,I) ELSEIF(R_CUR(SOLE(I),J).GE.BR(JMAX,I)) THEN QNT_NEW(SOLE(I),J)=(QNT(JN,I)-QNT(JN-1,I)) $*(R_CUR(SOLE(I),J)-BR(JN,I))/(BR(JN,I)-BR(JN-1,I))+QNT(JN,I) ELSE DO K=1,JN-1 IF((BR(K,I).LE.R_CUR(SOLE(I),J)).AND.(BR(K+1,I).GT.R_CUR(SOLE(I),J))) THEN QNT_NEW(SOLE(I),J)=(QNT(K+1,I)-QNT(K,I)) $*(R_CUR(SOLE(I),J)-BR(K,I))/(BR(K+1,I)-BR(K,I))+QNT(K,I) ENDIF ENDDO ENDIF ENDIF ENDDO ENDDO END SUBROUTINE CAL_STATE IMPLICIT NONE INCLUDE 'COMM.FOR' INTEGER I,J,K,N DOUBLE PRECISION W RADPS=2*PI*RPM/60.0 CALL INTP(VURI,VUR,2) CALL INTP(EFF_RS,EFF,2) DO I=1,NOBR-1 DO K=SOTE(I)+1,SOLE(I+1) DO J=1,JMAX VUR(K,J)=VUR(SOTE(I),J) ENDDO ENDDO ENDDO DO I=1,IMAX DO J=1,JMAX VU(I,J)=VUR(I,J)*1000/R_CUR(I,J) ENDDO ENDDO DO N=1,NOBR/2 DO J=1,JMAX DO I=2,IMAX IF((I.GE.SOLE(2*N-1)).AND.(I.LE.SOTE(2*N-1))) THEN W=RADPS ELSE W=0.0 ENDIF ENTRA(I,J)=ENTRA(I-1,J)+W*(VUR(I,J)-VUR(I-1,J)) TT(I,J)=ENTRA(I,J)/CP ENDDO ENTHO(SOTE(2*N-1),J)=ENTHO(SOLE(2*N-1),J)-CP*LOG(EFF(SOTE(2*N-1),J) $+(1-EFF(SOTE(2*N-1),J))*TT(SOLE(2*N-1),J)/TT(SOTE(2*N-1),J)) DO I=SOTE(2*N-1)+1,SOLE(2*N) ENTHO(I,J)=ENTHO(SOTE(2*N-1),J) ENDDO ENTHO(SOTE(2*N),J)=ENTHO(SOLE(2*N),J)-RG*LOG(EFF(SOTE(2*N),J)) IF(N.NE.(NOBR/2)) THEN DO I=SOTE(2*N)+1,SOLE(2*N+1) ENTHO(I,J)=ENTHO(SOTE(2*N),J) ENDDO ELSE DO I=SOTE(2*N)+1,IMAX ENTHO(I,J)=ENTHO(SOTE(2*N),J) ENDDO ENDIF ENDDO ENDDO DO I=SOLE(1)+1,IMAX DO J=1,JMAX TP(I,J)=TP(I-1,J)*EXP((CP*LOG(TT(I,J)/TT(I-1,J))-(ENTHO(I,J)-ENTHO(I-1,J)))/RG) ENDDO ENDDO DO I=1,IMAX DO J=1,JMAX VM(I,J)=WX(I,J)/COS(SIGMA(I,J)) VR(I,J)=WX(I,J)*TAN(SIGMA(I,J)) V(I,J)=SQRT(VM(I,J)**2.0+VU(I,J)**2.0) ST(I,J)=TT(I,J)-(VM(I,J)**2.0+VU(I,J)**2.0)/(2.0*CP) SP(I,J)=TP(I,J)*(ST(I,J)/TT(I,J))**(KG/(KG-1)) DENS(I,J)=SP(I,J)/(RG*ST(I,J)) FG(I,J)=2*PI*R_CUR(I,J)*DENS(I,J)*SIN(THETA(I,J)-SIGMA(I,J))/COS(SIGMA(I,J)) ENDDO ENDDO END SUBROUTINE CURVANTURE IMPLICIT NONE INCLUDE 'COMM.FOR' INTEGER I,J DOUBLE PRECISION THETA_TEMP(IMAX,JMAX),TEMP(IMAX,JMAX) DO J=1,JMAX DO I=1,IMAX-1 THETA_TEMP(I,J)=ATAN((R_CUR(I+1,J)-R_CUR(I,J))/ $(X_CUR(I+1,J)-X_CUR(I,J))) TEMP(I,J)=SQRT((R_CUR(I+1,J)-R_CUR(I,J))**2.0+(X_CUR(I+1,J)- $X_CUR(I,J))**2.0) ENDDO THETA_TEMP(IMAX,J)=THETA_TEMP(IMAX-1,J) TEMP(IMAX,J)=TEMP(IMAX-1,J) ENDDO DO J=1,JMAX DO I=2,IMAX-1 RM(I,J)=-2.0*(THETA_TEMP(I,J)-THETA_TEMP(I-1,J))/(TEMP(I,J)+TEMP(I-1,J)) SIGMA(I,J)=ATAN(0.5*((R_CUR(I+1,J)-R_CUR(I,J))/(X_CUR(I+1,J)- $X_CUR(I,J))+(R_CUR(I,J)-R_CUR(I-1,J))/(X_CUR(I,J)-X_CUR(I-1,J)))) ENDDO RM(1,J)=RM(2,J) RM(IMAX,J)=0.0 SIGMA(1,J)=ATAN((R_CUR(2,J)-R_CUR(1,J))/(X_CUR(2,J)-X_CUR(1,J))) SIGMA(IMAX,J)=ATAN((R_CUR(IMAX,J)-R_CUR(IMAX-1,J))/ $(X_CUR(IMAX,J)-X_CUR(IMAX-1,J))) ENDDO END SUBROUTINE SL_ADJ IMPLICIT NONE INCLUDE 'COMM.FOR' INTEGER I,J,M,IPASS,K,IPASSL PARAMETER (IPASSL=50) DOUBLE PRECISION AREA(IMAX),WX_TEMP(IMAX),ERROR,ADENS(IMAX),YITA_TEMP(IMAX,JMAX) PARAMETER (ERROR=1E-4) DO I=1,IMAX DO K=1,IPASSL DO M=1,IMAX YITA(M,1)=0.0 YITA_NEW(M,1)=0.0 DO J=2,JMAX YITA(M,J)=SQRT((R_CUR(M,J)-R_CUR(M,1))**2.0+(X_CUR(M,J)-X_CUR(M,1))**2.0) YITA_NEW(M,J)=YITA(M,J) ENDDO AREA(M)=PI*(R_CUR(M,JMAX)**2.0-R_CUR(M,1)**2.0) ENDDO CALL SLC_CAL DETG(I,1)=0.0 DO J=2,JMAX DETG(I,J)=DETG(I,J-1)+0.5*WBLK(I)*FG(I,J)*(WX(I,J)+WX(I,J-1))*(YITA_NEW(I,J)-YITA_NEW(I,J-1)) ENDDO ADENS(I)=0.0 DO J=1,JMAX ADENS(I)=ADENS(I)+DENS(I,J) ENDDO ADENS(I)=ADENS(I)/JMAX IF((ABS(DETG(I,JMAX)/(MF*1E6)-1).LT.(0.05*ERROR))) THEN EXIT ELSE WX_TEMP(I)=WX(I,1)-0.9*(DETG(I,JMAX)-MF*1E6)/(AREA(I)*ADENS(I)) WX(I,1)=WX_TEMP(I) ENDIF IF((K.EQ.IPASSL).AND.(ABS(DETG(I,JMAX)/MF-1).GT.(0.05*ERROR))) THEN ENDIF ENDDO DO J=2,JMAX-1 YITA_TEMP(I,J)=2.0*DETG(I,JMAX)/((JMAX-1)*FG(I,J)*WBLK(I)*(WX(I,J)+WX(I,J-1)))+YITA(I,J-1) YITA_NEW(I,J)=YITA(I,J)+0.5*(YITA_TEMP(I,J)-YITA(I,J)) X_ORI(I,J)=X_CUR(I,J) R_ORI(I,J)=R_CUR(I,J) X_CUR(I,J)=YITA_NEW(I,J)*COS(THETA(I,J))+X_ORI(I,1) R_CUR(I,J)=YITA_NEW(I,J)*SIN(THETA(I,J))+R_ORI(I,1) ENDDO ENDDO OPEN(30,FILE='MESS_ADJ.plt') WRITE(30,*) 'TITLE="MESS_。