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

作业条件和要求 设计参数 进气总压 101325 0 Pa 0 P 进气总温 288 15 K 0 T 质量流量G 550 kg s 总增压比 1 5 k 绝热效率 0 90 k 气体常数R 287 06 J kg K 定压比热 1004 7 J kg KCp 比热比 k 1 4 确定 说明理由 转速 流道几何 各排叶片参数的展向分布 最少根 中 尖三个截面 叶 片尾缘半径 D 因子 压比 效率 轮缘功 子午速度 相对速度 绝对速度 相对气流角 绝对气流角 静压 静温 总压 总温 密度 各叶片的基元几 何 最少根 中 尖三个截面 叶轮机高等气动力学 大作业 风扇设计风扇设计 院 系 名称 专业名称 学号 学生姓名 一 控制方程 1 连续性方程 连续性方程 积分形式的连续性方程 对展向计算站而言 tip GGx hub GVAKF W d A 2 cos G sin FR 2 运动方程 运动方程 以焓熵形式描述的展向平衡方程 12 x xx W FWFW 1 lnsin sincoscos cos cos x m DWtg F rdx 2 2 cos uu vv ris FT r 3 能量方程 能量方程 212 21 1 uu iiv rv r 4 状态方程 状态方程 pRT 5 熵增关系式 熵增关系式 1 21 2 21 ln1 ln pRR s T ssc T ssR 6 沿流线斜率曲率 沿流线斜率曲率 2 2 11 3 2 2 1 tantan 1 s s m s d f dfxdrD dx dxdxrdm df dx 二 数值过程二 数值过程 1 离散方式 离散方式 i 为计算站标号 j 为流线标号 1 11 111 1 2 jj iiii ijjiiii j FF FFFFFDF dxxxxx 2 运动方程 运动方程 流线平均参数 1 1 2 jjj FFF 1 1121 1 2 2 xjxj xjxjjj xjxj WW WWFF WW 3 连续方程 连续方程 1 11 2 xjxj jjGiGjjj WW GGKF 1 0 JN G GG 三 求解流场流程图 2m S 四 压气机设计 1 压气机级数的确定 压气机级数的确定 压气机设计流量为 550kg s 压比 1 5 可以看出 属于民用风扇范畴 因此初步选定为单级轴流压气机 一级转子与一级静子的 形式 2 压气机流道形式的选择 压气机流道形式的选择 对于级数较少的压气机而言 流道采用等 外径设计规律 可以充分利用压气机较高的叶尖切线速度 实现较大的加功量 因此本压气机选择等外径气流通道设计 3 压气机轮毂比的确定 压气机轮毂比的确定 转子轮毂比的选择需同时考虑气动性能和结构 强度等方面的因素 较小的进口轮毂比对气动性能有好处 可以提高风扇的效 率 但不利于转子的结构强度 且过小的轮毂比会导致根部的加功量足 对于 发动机进口机风扇而言 其轮毂比多在 0 32 0 4 之间 因此选定压气机轮毂比 为 0 35 并定义转子进口机匣与轮毂半径分别为 且满足 t R h R 0 35 ht RR 4 4 压气机进口参数的确定 压气机进口参数的确定 压气机进口轴向马赫数大小与压气机性能有密切的关系 当进口轴向马赫数 增加 压气机效率和裕度都要下降 根据以往经验 高压压气机进口马赫数通 常为 0 48 0 52 之间 对于风扇 进口马赫数通常在 0 6 0 68 之间 因此本设计选定进口马赫数 Ma 0 6 并认为进口流场为均匀场 来流方向轴向 无预旋 进口静温 由 可求得 268 80 K 2 1 1 2 TK Ma T 1 T 当地音速 c 328 67m s m s kRT 1 0 6197 20 a Ca 进口静压 由 计算得79439 20Pa 2 1 1 1 2 k k Pk Ma P 1 P 进口空气密度 1 030 1 1 1 P RT 3 kg m 空气质量流量 G 550 kg s 由连续方程 风扇进口面积 2 1 11 A2 708m a G C 由轮毂比和面积公式 可得 0 991m 0 347m 22 1 th ARR 1 R 1 r 5 压气机出口参数的确定 压气机出口参数的确定 风扇出口 Ma 数范围为 0 48 0 5 选定出口 Ma 为 0 48 压气机级出口总压101325 1 5 151987 50Pa 由压气 30 PP 机绝热效率 可求得压气机出口总温为327 47K 1 31 31 1 1 k k PP TT 3 T 313 04K 129816 75Pa 1 45 3 T 1 3 33 3 k k T PP T 3 3 3 P RT 3 kg m 由压气机出口截面积以及出口外径 可得压气 3 170 25 a vm s 3 A 3 0 991 Rm 机出口轮毂半径为 3 0 5224rm 6 转子出口截面参数确定 转子出口截面参数确定 确定静子总压恢复系数 转子 23 0 98 出口总压 转子压比 1 531 由转子效率公式 2 P155089 29Pa R 可得转子效率 0 948 因此由轮缘功公式可得 1 2 R 1 11 k k R T T 39518 53J kg 由单级风扇功分配的特点 应适当 K 1 K 1 1 1 Lu k k RT k k 减小根部和尖部的加功量 对于叶根处 加功量应尽可能小 以保证根部气流 转折角不至于过大 对于尖部 同样应减小加功量 来减小尖部气流的分离 由轮缘功的公式可知 可通过给定压比和效率沿展向的分布来保证不同展高位 置的功分配 表 1 给出了不同截面处 转子叶片轮缘功压比和效率的分配情况 图 1 给出了不同截面处压比分配曲线 表 1 转子叶片不同展高位置参数分布 百分比展高压比效率轮缘功 J kg 01 40 932458 3 5 1 4250 934253 6 10 1 450 936026 5 20 1 50 939507 8 30 1 540 942233 7 50 1 59310 945775 1 70 1 6050 946557 2 90 1 570 91543519 7 1001 450 9235243 3 平均值 1 5430 905839508 4 R 0 00 0 10 0 20 0 30 0 40 0 50 0 60 0 70 0 80 0 90 1 00 1 3001 4001 5001 6001 700 压比 相对展高 转子出口环量展向分布 0 0 10 0 20 0 30 0 40 0 50 0 60 0 70 0 80 0 90 0 100 0 0 050 0100 0150 0 环量 相对展高 图 1 转子压比沿展向分布 图 2 转子出口环量展向分布 7 叶尖切线速度和转速的确定 叶尖切线速度和转速的确定 根据轮缘功公式 可知 切线速度越大轮缘功越大 有利于压比的提高 当轮 uuu wwu rL 缘功不变时 切线速度越大 扭速越小 即叶片的弯度越小 当然 叶尖切线 速度过大时 尖部相对马赫数过高 激波损失增强 从而导致效率的下降 这 对于现代压气机设计来说已经成为了不可避免的问题 对于相同的角速度 为 得到更大的加功量 根部需要更高的扭速 即叶片的弯角更大 由公式 知 切线速度与转速密切关联 本压气机选取叶尖切线速度为 30 u nr r 380m s 其中 r 为转子外径 所以转速 n 3661 69rpm 角速度 所以对应的根中尖三个截面处的负荷系数分别为 389 733 rad s 2 23 0 601 和 0 274 由压气机出口气流速度 所以对应根中 3 170 25 a vm s 尖的流量系数分别为 1 280 0 864 0 448 图 2 给出了本压气机尖部流量系数 负荷系数的分配在 Smith 图上的位置 图中 可以看到 本压气机 2 2 Lu U 中部和尖部的流量系数与负荷系数的分配比较合理 坐标分别为 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 因子的大小 因此对 112 21wwwwD u 转子根中尖稠度分别取 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 avg Zrb 表 2 所示为根据以上计算所得压气机各相关参数 表 2 压气机各相关参数 动叶根中尖 展弦比1 822 5 弦长 b m 0 3800 3420 274 稠度 b t2 31 71 2 栅距 t m 0 1650 2010 228 D 因子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 232 D 因子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 静压 Pa 534 18 0 00 155 37 220 06163 96 1 40 0 90 1 27 0 43 320 44 296 34 141834 2107877 2 608 82 14 41 151 66 208 88178 241 42 0 90 1 31 0 42 322 07 300 36 144143 5112900 8 674 79 27 15 149 45 201 82196 35 1 44 0 90 1 34 0 41 323 65 303 38 146405 5116748 3 734 10 38 61 148 08 200 46211 491 49 0 90 1 37 0 41 326 63 306 63 150735 120830 7 788 35 49 08 147 22 198 21 228 341 52 0 90 1 40 0 40 328 73 309 18 153859 9124142 1 838 59 58 79 146 62 195 91245 501 54 0 90 1 42 0 39 330 42 311 31 156383 2 126964 885 38 67 82 146 29 195 78260 011 58 0 90 1 45 0 39 332 84 313 76 160069 4130195 8 929 45 76 33 146 11 193 62276 891 60 0 90 1 46 0 38 333 96 315 30 161794 7132306 2 971 34 84 42 145 94 190 58 294 681 60 0 90 1 47 0 36 334 33 316 26 162377 8 133672 1011 50 92 18 145 70 185 24 315 511 59 0 91 1 48 0 34 333 03 315 95 160923 8133848 3 1052 00 100 00 144 70 168 31 354 881 45 0 92 1 42 0 26 323 23 309 13 146919 3125686 7 表 4 静子出口各参数沿展向分布 半径 mm 相对展高 子午速度 m s 绝对速度 m s 相对速度 m s 密度D 因子总温 K 静温 K 总压 Pa 静压 Pa 601 970 00 155 70155 70155 701 390 03320 44308 38 140415 9122767 662 8213 52 155 32155 32155 321 400 03322 07310 07 142702 1124935 4 718 3025 85 155 14155 14155 141 420 02323 65311 67 144941 4127019 5 769 0137 12 155 18155 18155 181 450 03326 63314 64 149227 7130927 6 816 0347 56 155 34155 34155 341 470 03328 73316 72 152321 3133718 1 860 0657 35 155 50155 50155 501 490 02330 42318 38 154819 4135967 2 901 3466 52 155 73155 73155 731 510 03332 84320 77 158468 7139250 6 940 5375 23 155 97155 97155 971 520 02333 96321 85 160176 8140758 4 978 0883 58 156 06156 06156 061 530 01334 33322 21160754141263 4 1014 4991 66 155 97155 97155 971 520 00333 03320 92 159314 5139948 5 1052 00 100 00 155 01155 01155 011 430 00323 23311 27 145450 1127471 3 图 5 相对马赫数等值线 图 6 静压分布等值线 六 叶片造型 转子采用圆弧中弧线 由于转子叶片跨音 所以对尖部以下区域采用双圆 弧叶型 适应来流马赫数为 0 8 1 2 的跨音流动 其最大厚度相对位置和最 大挠度相对位置均为弦长的 50 采用重心积叠的方式 但由于本压气机转子 叶尖相对马赫数最大值超过了 1 2 所以对转子叶尖采用多圆弧叶型 最大厚度 相对位置和最大挠度相对位置均为弦长的 50 转子进出口绝对气流角和相对 气流角沿展向分布如下表 5 所示 表 5 转子进出口绝对 相对气流角 Y 进口绝对气流角进口相对气流角出口绝对气流角出口相对气流角 0 00 0 00 37 9145 09 18 62 49 08 0 00 60 0342 03 49 85 100 00 0 00 68 8830 07 65 94 叶片根部和尖部给定攻角 2 中部给定攻角 2 由叶型弯角公式 1 iX t m b 其中 2 2 0 92 0 0020 18ma 可计算得出根中尖位置处的叶型弯角 并由此画出中弧线 其中 根部基 元 叶中基元 尖部基元 计算出根部31 62 hub 15 07 mid 11 76 tip 叶型安装角为 26 26 中部叶型安装角 56 9 尖部叶型安装角为 65 4 因 此可画出转子根中尖三个截面的基元 图 7 给出了转子叶片根中尖三个基元截 面的叶型 图 7 转子叶片根中尖截面叶型 对静子叶片采用 NACA65 010 叶型 采用尾缘积叠的方式 并采用圆弧 中弧线 同理 可计算得静子叶片的叶型弯角 根中尖分别为 58 9 hub 根中尖截面的安装角分别为 14 93 11 39 59 8 mid 43 6 tip 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 3RG 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