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

作业条件和要求

设计参数 进气总压 进气总温 质量流量 总增压比 绝热效率 气体常数 定压比热

*=101325.0 Pa P0*=288.15 K T0G=550 kg/s

*=1.5 ?k*=0.90 ?kR= 287.06 J/kg/K

Cp= 1004.7 J/kg/K

比热比 k = 1.4 确定（说明理由）

转速，流道几何，各排叶片参数的展向分布（最少根、中、尖三个截面）：叶片

尾缘半径，D因子，压比，效率，轮缘功，子午速度，相对速度，绝对速度，相对气流角，绝对气流角，静压，静温，总压，总温，密度，各叶片的基元几何（最少根、中、尖三个截面）

叶轮机高等气动力学

大作业

——风扇设计

院（系）名称 专业名称 学

号

学生姓名

一、控制方程

（1）连续性方程

积分形式的连续性方程（对展向计算站而言）

G??VA?KG?FG?2?R??tip?hubFG?Wxd?

sin(???)

cos?（2）运动方程

以焓熵形式描述的展向平衡方程

?Wx?FW1x?F2/Wx ??F1??DlnWxsin??tg??sin?cos??cos(???)cos?

rmcos???dx?vu?(vur)?i*?s? F2?cos?????T???????r??2（3）能量方程

**i2?i1??(vu2r2?vu1r1)

（4）状态方程

p??RT

（5）熵增关系式

*?*?*T1s2?s1?cpln??R??1??R?*?T2? ?s2?s1?Rln?s*（6）沿流线斜率曲率

??tan?1?1D??dr??1?dfs?x?? ?tan ???????3rmdm?dx??dx???dfs?2?2?1????dx??????d2fsdx2二、数值过程

（1）离散方式（i为计算站标号，j为流线标号）

Fj?Fj?1??F??DF?1?Fi?1?FiFi?Fi?1?? ?? ???????????dx2x?xx?x??i??jjj?1ii?1??i?1i（2）运动方程

流线平均参数 Fj?1?Fj?F?j1? 2?Wxj?Wxj?1?2Wxj?Wxj?1??F1?F2???j??j?1?

2Wxj?Wxj?1????（3）连续方程

?Gj??Gj?1?KGi?FGj?Wxj?Wxj?12??j??j?1?

??G1?0 ???GJN?G三、求解S2m流场流程图

四、压气机设计

（1）压气机级数的确定：压气机设计流量为550kg/s，压比1.5，可以看出

属于民用风扇范畴，因此初步选定为单级轴流压气机。一级转子与一级静子的形式。

（2）压气机流道形式的选择：对于级数较少的压气机而言，流道采用等外

径设计规律，可以充分利用压气机较高的叶尖切线速度，实现较大的加功量。因此本压气机选择等外径气流通道设计。

（3）压气机轮毂比的确定：转子轮毂比的选择需同时考虑气动性能和结构

强度等方面的因素．较小的进口轮毂比对气动性能有好处，可以提高风扇的效率，但不利于转子的结构强度，且过小的轮毂比会导致根部的加功量足。对于发动机进口机风扇而言，其轮毂比多在0.32-0.4之间，因此选定压气机轮毂比为0.35。并定义转子进口机匣与轮毂半径分别为Rt,Rh，且满足Rh/Rt?0.35。

（4）压气机进口参数的确定：

压气机进口轴向马赫数大小与压气机性能有密切的关系，当进口轴向马赫数

增加，压气机效率和裕度都要下降。根据以往经验，高压压气机进口马赫数通常为0.48-0.52之间，对于风扇，进口马赫数通常在0.6-0.68之间。

因此本设计选定进口马赫数Ma=0.6，并认为进口流场为均匀场，来流方向轴向，无预旋。

T*K?1?1?Ma2，可求得T1= 268.80 K 进口静温：由T2当地音速：c=kRT=328.67m/s，?C1a?0.6a?197.20m/s。

*进口静压：由P?(1?k?1Ma2)k?1 计算得P1?79439.20Pa

P2k进口空气密度：?1?P1?1.030kg/m3 RT1空气质量流量：G=550 kg/s 由连续方程，风扇进口面积A1?G?2.708m2 ?1C1a由轮毂比和面积公式A1??(Rt2?Rh2)，可得R1=0.991m，r1=0.347m。

（5）压气机出口参数的确定：风扇出口Ma数范围为0.48-0.5，选定出口

***Ma为0.48。压气机级出口总压P3?P0*??101325*1.5=151987.50Pa，由压气机

k?1k绝热效率?*?(P/P)?1，可求得压气机出口总温为T3*?327.47K。**(T3/T1)?1kk?1*3*1?T??T3?313.04K，P3??3*??T3?P3*?129816.75Pa，?3?P3?1.45kg/m3，RT3v3a?170.25m/s。由压气机出口截面积A3以及出口外径R3=0.991m，可得压气机出口轮毂半径为r3?0.5224m。

（6）转子出口截面参数确定：确定静子总压恢复系数?23=0.98，?转子出

**口总压P2子=1.531，由转子效率公式?155089，.2转9Pa压比?Rk?1??*k?R???R?1????T2*??*?1?可得转子效率0.948，因此由轮缘功公式可得?T1?k?1k?*?Lu?RT1??Kk?1??K=39518.53J/kg。由单级风扇功分配的特点，应适当k?1??减小根部和尖部的加功量：对于叶根处，加功量应尽可能小，以保证根部气流转折角不至于过大；对于尖部，同样应减小加功量，来减小尖部气流的分离。由轮缘功的公式可知，可通过给定压比和效率沿展向的分布来保证不同展高位置的功分配。表1给出了不同截面处，转子叶片轮缘功压比和效率的分配情况，图1给出了不同截面处压比分配曲线。

表1 转子叶片不同展高位置参数分布

百分比展高 0 5% 10% 20% 30% 50% 70% 90% 100 平均值

压比 1.4 1.425 1.45 1.5 1.54 1.5931 1.605 1.57 1.45 1.543 效率 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.915 0.92 0.9058 轮缘功J/kg 32458.3 34253.6 36026.5 39507.8 42233.7 45775.1 46557.2 43519.7 35243.3 39508.4 R1.000.900.800.700.600.500.400.300.200.100.001.300转子出口环量展向分布1.4001.500压比1.6001.700100.0?.0?.0p.0`.0P.0@.00.0 .0.0%0.0%0.0相对展高相对展高50.0100.0环量150.0

图1转子压比沿展向分布 图2 转子出口环量展向分布

（7）叶尖切线速度和转速的确定：根据轮缘功公式：Lu?u?wu?r??wu可

知，切线速度越大轮缘功越大，有利于压比的提高，当轮缘功不变时，切线速度越大，扭速越小，即叶片的弯度越小。当然，叶尖切线速度过大时，尖部相对马赫数过高，激波损失增强，从而导致效率的下降，这对于现代压气机设计来说已经成为了不可避免的问题。对于相同的角速度，为得到更大的加功量，根部需要

r?n更高的扭速，即叶片的弯角更大。由公式：u?r??知，切线速度与转速密

30切关联，本压气机选取叶尖切线速度为380m/s，其中r为转子外径，所以转速n=3661.69rpm，角速度?=389.733rad/s，所以对应的根中尖三个截面处的负荷系数分?别为2.23、0.601和0.274。由压气机出口气流速度v3a?170.25m/s，所以对应根中尖的流量系数?分别为：1.280,0.864,0.448。图2给出了本压气机尖部流量系数负荷系数的分配在Smith图上的位置（图中?=2Lu/U2）。可以看到，本压气机中部和尖部的流量系数与负荷系数的分配比较合理（坐标分别为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?1?w2w1??wu?2w1??，可知稠度的取值能够影响到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。

根据叶片数公式Z?2?ravg?/b可求得转子叶片数为23片，静子为24片。 表2所示为根据以上计算所得压气机各相关参数。

表2 压气机各相关参数

动叶 展弦比 弦长b(m) 稠度τ=b/t 栅距t(m) D因子 轴向长度Δz(m) 叶片数 进口速度V1（m/s） 进口总压P1*（Pa） 进口静压P1（Pa) 进口总温T1*(K) 进口静温T1(K) 出口总压P2*(Pa) 出口静压P2(Pa) 出口总温T2*(K) 出口静温T2(K) 出口轴向速度V2a（m/s) 压比π* 效率η* 展弦比 弦长b(m) 稠度τ=b/t 栅距t(m) D因子 轴向长度Δz(m) 叶片数 出口速度V3（m/s） 出口总压P3*（Pa) 出口静压P3（Pa) 出口总温T3*(K) 出口静温T3(K) 170.25 151987.5 137494.26 327.47 313.05 166.05 101325 87538.94 288.15 274.43 140841.75 105120.95 319.74 284.18 150.13 1.39 0.90 静叶 2 0.251 1.9 0.132 0.363 0.243 2 0.251 1.3 0.193 0.242 0.303 24 170.25 151987.5 137494.26 327.47 313.05 170.25 151987.5 137494.26 327.47 313.05 1.8 0.278 1.2 0.232 0.021 根 1.8 0.380 2.3 0.165 0.449 0.331 中 2 0.342 1.7 0.201 0.576 0.332 23 166.05 101325 87538.94 288.15 274.43 162120 136899.32 334.16 309.06 150.13 1.6 0.90 166.05 101325 87538.94 288.15 274.43 146921.25 131791.70 324.01 308.95 150.13 1.45 0.90 尖 2.5 0.274 1.2 0.228 0.279

（9）计算站给定及流道几何设计：根据以上数据，初步设计压气机流道

几何如图3所示。并在转子叶片之前给定9个计算站，转静之间给定两个计算站，静子后给定4个计算站，没有给定叶片内部计算站。

图4 压气机流道图

五、通流计算结果

转子出口截面各参数展向分布如表3所示；静子出口截面个参数展向分布如表4所示。图5、图6分别给出了通流计算全流场的相对马赫数等值线和静压等值线分布。

表3 转子出口截面各参数沿展向分布 半径 mm 相对展高 子午速度 绝对速度 相对速度 压比 效率 密度 D因子 总温 m/s m/s m/s m/s K 静温 K 总压 Pa 静压 Pa 534.18 0.00% 608.82 14.41% 674.79 27.15% 734.10 38.61% 788.35 49.08% 838.59 58.79% 885.38 67.82% 929.45 76.33% 971.34 84.42% 1011.50 92.18% 155.37 220.06 163.96 1.40 0.90 1.27 0.43 320.44 296.34 141834.2 107877.2 151.66 208.88 178.24 1.42 0.90 1.31 0.42 322.07 300.36 144143.5 112900.8 149.45 201.82 196.35 1.44 0.90 1.34 0.41 323.65 303.38 146405.5 116748.3 148.08 200.46 211.49 1.49 0.90 1.37 0.41 326.63 306.63 150735 120830.7 147.22 198.21 228.34 1.52 0.90 1.40 0.40 328.73 309.18 153859.9 124142.1 146.62 195.91 245.50 1.54 0.90 1.42 0.39 330.42 311.31 156383.2 126964 146.29 195.78 260.01 1.58 0.90 1.45 0.39 332.84 313.76 160069.4 130195.8 146.11 193.62 276.89 1.60 0.90 1.46 0.38 333.96 315.30 161794.7 132306.2 145.94 190.58 294.68 1.60 0.90 1.47 0.36 334.33 316.26 162377.8 133672 145.70 185.24 315.51 1.59 0.91 1.48 0.34 333.03 315.95 160923.8 133848.3 1052.00 100.00% 144.70 168.31 354.88 1.45 0.92 1.42 0.26 323.23 309.13 146919.3 125686.7

表4 静子出口各参数沿展向分布

半径 相对展高 子午速度 绝对速度 相对速度 密度 D因子 mm m/s m/s m/s 601.97 0.00% 155.70 155.70 155.32 155.14 155.18 155.34 155.50 155.73 155.97 156.06 155.97 155.01 155.70 155.32 155.14 155.18 155.34 155.50 155.73 155.97 156.06 155.97 155.01 1.39 1.40 1.42 1.45 1.47 1.49 1.51 1.52 1.53 1.52 1.43 0.03 0.03 0.02 0.03 0.03 0.02 0.03 0.02 0.01 0.00 0.00 662.82 13.52% 155.32 718.30 25.85% 155.14 769.01 37.12% 155.18 816.03 47.56% 155.34 860.06 57.35% 155.50 901.34 66.52% 155.73 940.53 75.23% 155.97 978.08 83.58% 156.06 1014.49 91.66% 155.97 1052.00 100.00% 155.01 总温 静温 K K 总压 Pa 静压 Pa 320.44 308.38 140415.9 122767 322.07 310.07 142702.1 124935.4 323.65 311.67 144941.4 127019.5 326.63 314.64 149227.7 130927.6 328.73 316.72 152321.3 133718.1 330.42 318.38 154819.4 135967.2 332.84 320.77 158468.7 139250.6 333.96 321.85 160176.8 140758.4 334.33 322.21 160754 141263.4 333.03 320.92 159314.5 139948.5 323.23 311.27 145450.1 127471.3

图5 相对马赫数等值线 图6 静压分布等值线

六、叶片造型

转子采用圆弧中弧线，由于转子叶片跨音，所以对尖部以下区域采用双圆弧叶型，适应来流马赫数为0.8——1.2的跨音流动，其最大厚度相对位置和最大挠度相对位置均为弦长的50%，采用重心积叠的方式。但由于本压气机转子叶尖相对马赫数最大值超过了1.2，所以对转子叶尖采用多圆弧叶型。最大厚度相对位置和最大挠度相对位置均为弦长的50%。转子进出口绝对气流角和相对气流角沿展向分布如下表5所示。

表5 转子进出口绝对/相对气流角 Y% 0.00% 49.08% 100.00% 进口绝对气流角 进口相对气流角 出口绝对气流角 出口相对气流角 0.00 0.00 0.00 -37.91 -60.03 -68.88 45.09 42.03 30.07 -18.62 -49.85 -65.94 叶片根部和尖部给定攻角-2°，中部给定攻角2°，由叶型弯角公式： ??*?i*??X ??t1?mb*其中 m?0.92(a)2?0.002?2?0.18

可计算得出根中尖位置处的叶型弯角，并由此画出中弧线。其中，根部基元

?hub?31.62°，叶中基元?mid?15.07°，尖部基元?tip?11.76°计算出根部叶型

安装角为26.26°，中部叶型安装角56.9°，尖部叶型安装角为65.4°，因此可画出转子根中尖三个截面的基元。图7给出了转子叶片根中尖三个基元截面的叶型。

图7转子叶片根中尖截面叶型

对静子叶片采用NACA65—010叶型，采用尾缘积叠的方式，并采用圆弧中弧线。同理，可计算得静子叶片的叶型弯角，根中尖分别为?hub??58.9°，

?mid??59.8°，?tip??43.6°；根中尖截面的安装角分别为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 W=0.0 ELSE 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

END

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

WRITE(*,*) 'CONVERGENCE NOT REACHED...'

ELSEIF(IPASS.EQ.IPASSMAX) THEN ENDIF

ENDDO

SUBROUTINE OUTPUT IMPLICIT NONE INCLUDE 'COMM.FOR' INTEGER I,J

OPEN(30,FILE='FIELD.plt') WRITE(30,*) 'TITLE=\

WRITE(30,'(A300)') 'VARIABLES=\

!,\

WRITE(30,*) 'ZONE T=\

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)

$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 RECOVERY COEFFI...

$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

COMMON /STATE/ DENS,ST,SP,ENTRA,ENTHO,EFF_RS,EFF,TT,TP,FG DOUBLE

PRECISION

ROTER-IFFICIENCY,FOR

STATOR

TP

!SECTION OF L.E. COMMON

DOUBLE PRECISION DENS(IMAX,JMAX),ST(IMAX,JMAX),SP(IMAX,JMAX),

ENDWALL BLOCKAGE

/GEOM/

R_ORI,X_ORI,R_CUR,X_CUR,RM,BR,SLD,SOLE,SOTE,WBLK,SLDT,YITA,YITA_NEW

VX(IMAX,JMAX),WX_OLD(IMAX,JMAX),VR(IMAX,JMAX),VM(IMAX,JMAX),V_REL(IMAX,JMAX),

!DELTA_G DOUBLE

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

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

READ(30,*)

PRECISION

SIGMA(IMAX,JMAX),THETA(IMAX,JMAX),ARPHA(IMAX,JMAX),BETA(IMAX,JMAX)

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_O

ENDDO DO J=1,JMAX

X_CUR(I,J)=X_ORI(I,J) R_CUR(I,J)=R_ORI(I,J)

RI(I,1))

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=\

WRITE(30,'(A300)') 'VARIABLES=\ WRITE(30,*) 'ZONE T=\

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 DO J=1,JMAX

THETA(I,J)=ATAN(TEMP(I))+PI ENDDO

ELSE

ENDIF

!SOLE(1)

ENDDO DO I=1,IMAX

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 DO I=1,NOBR

DO J=1,JMAX

IF(R_CUR(SOTE(I),J).LT.BR(1,I)) THEN

QNT_NEW(SOTE(I),J)=(QNT(2,I)-QNT(1,I))

IF(FLG.EQ.2) THEN

!FLG=1--L.E.;FLG=2--T.E.

$*(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

IF(R_CUR(SOLE(I),J).LT.BR(1,I)) THEN

QNT_NEW(SOLE(I),J)=(QNT(2,I)-QNT(1,I))

ELSEIF(FLG.EQ.1) THEN

$*(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

END

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

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 W=0.0 ELSE 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) 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

DO I=SOTE(2*N)+1,IMAX

ENTHO(I,J)=ENTHO(SOTE(2*N),J)

ENDDO

$+(1-EFF(SOTE(2*N-1),J))*TT(SOLE(2*N-1),J)/TT(SOTE(2*N-1),J))

ELSE

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 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))/ 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

$(X_CUR(I+1,J)-X_CUR(I,J)))

DO I=1,IMAX

DO K=1,IPASSL

YITA_NEW(M,J)=YITA(M,J) ENDDO

AREA(M)=PI*(R_CUR(M,JMAX)**2.0-R_CUR(M,1)**2.0) DO M=1,IMAX

YITA(M,1)=0.0 YITA_NEW(M,1)=0.0 DO J=2,JMAX

SUBROUTINE SL_ADJ IMPLICIT NONE INCLUDE 'COMM.FOR'

INTEGER I,J,M,IPASS,K,IPASSL PARAMETER (IPASSL=50) DOUBLE

PARAMETER (ERROR=1E-4)

PRECISION

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

AREA(IMAX),WX_TEMP(IMAX),ERROR,ADENS(IMAX),YITA_TEMP(IMAX,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)

ENDDO

1))

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-

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 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

WX_TEMP(I)=WX(I,1)-0.9*(DETG(I,JMAX)-MF*1E6)/(AREA(I)*ADENS(I)) WX(I,1)=WX_TEMP(I) ELSE

ENDIF

IF((K.EQ.IPASSL).AND.(ABS(DETG(I,JMAX)/MF-1).GT.(0.05*ERROR))) THEN ENDIF ENDDO DO J=2,JMAX-1

ENDDO

ENDDO

OPEN(30,FILE='MESS_ADJ.plt') WRITE(30,*) 'TITLE=\

WRITE(30,'(A300)') 'VARIABLES=\ WRITE(30,*) 'ZONE T=\

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) END

SUBROUTINE SLC_CAL

IMPLICIT NONE INCLUDE 'COMM.FOR' INTEGER IPASS,I,J

DOUBLE PRECISION F1_TEMP(IMAX,JMAX),F2_TEMP(IMAX,JMAX),TEMP(IMAX,JMAX) DOUBLE PRECISION F1(IMAX,JMAX),F2(IMAX,JMAX) DOUBLE PRECISION WX_TEMP(IMAX,JMAX) CALL CURVANTURE CALL CAL_STATE DO I=1,IMAX

DO J=1,JMAX

WX_OLD(I,J)=WX(I,J) ENDDO

ENDDO DO J=1,JMAX

DO I=1,IMAX-1

TEMP(I,J)=(LOG(WX(I+1,J))-LOG(WX(I,J)))/(X_CUR(I+1,J)-X_CUR(I,J)) ENDDO

TEMP(IMAX,J)=TEMP(IMAX-1,J)

ENDDO DO J=1,JMAX-1

DO I=1,IMAX

IF(I.NE.1) THEN

F1_TEMP(I,J)=-SIN(THETA(I,J))*RM(I,J)/COS(SIGMA(I,J))-

$(TAN(SIGMA(I,J+1))-TAN(SIGMA(I,J)))/(YITA_NEW(I,J+1)-YITA_NEW(I,J))* $SIN(SIGMA(I,J))*COS(SIGMA(I,J))+0.5*(TEMP(I,J)+TEMP(I-1,J))* $COS(THETA(I,J)-SIGMA(I,J))*COS(SIGMA(I,J))

F2_TEMP(I,J)=(-VU(I,J)/R_CUR(I,J)*((VUR(I,J+1)-VUR(I,J))

$/(YITA_NEW(I,J+1)-YITA_NEW(I,J)))+(ENTRA(I,J+1)-ENTRA(I,J)) $/(YITA_NEW(I,J+1)-YITA_NEW(I,J))-(ENTHO(I,J+1)-ENTHO(I,J)) $/(YITA_NEW(I,J+1)-YITA_NEW(I,J)))*COS(SIGMA(I,J))**2.0

ELSE

F1_TEMP(I,J)=F1_TEMP(2,J) F2_TEMP(I,J)=F2_TEMP(2,J) ENDIF

ENDDO

ENDDO DO I=1,IMAX

DO J=2,JMAX-1

F1(I,J+1)=0.5*(F1_TEMP(I,J)+F1_TEMP(I,J-1)) F2(I,J+1)=0.5*(F2_TEMP(I,J)+F2_TEMP(I,J-1))

ENDDO

F1(I,2)=F1_TEMP(I,1) F2(I,2)=F2_TEMP(I,1)

END

DO J=2,JMAX

WX_TEMP(I,J)=WX(I,J-1)+(0.5*F1(I,J)*(WX(I,J)+WX(I,J-1))+

!-FR(I,J)*SIN(THETA(I,J))-FX(I,J)

WX(I,J)=WX_TEMP(I,J)

$2.0*F2(I,J)/(WX(I,J)+WX(I,J-1)))

ENDDO

ENDDO

本文来源：https://www.bwwdw.com/article/mkxr.html