考虑端区黏性和径向掺混的流线曲率通流算法研究

刘志伟，季路成，蔡 军

（1.中国科学院工程热物理研究所，北京100190；2.北京理工大学宇航学院，北京100081）

考虑端区黏性和径向掺混的流线曲率通流算法研究

刘志伟1，季路成2，蔡 军1

（1.中国科学院工程热物理研究所，北京100190；2.北京理工大学宇航学院，北京100081）

高负荷多级压气机的技术难点之一在于如何快速、完整、准确地获得全工况下叶片前尾缘气动参数分布和整机特性结果，而通流计算凭借速度快、指向性明确的优点，一直是多级压气机设计的核心技术。为了提高通流代码的适用范围和预测精度，推导了统一适用于轴流、斜流、离心压气机的流线曲率方程，引入端区黏性和径向掺混模型，应用Fortran语言针对轴流压气机全新开发了其通流特性计算程序，对NASA Rotor37、Rotor67转叶和GY1-2J压气机进行了计算，计算结果表明：新模型的加入提高了程序的计算精度，初步验证了程序的有效性。并提出提高程序的通用性是今后通流计算的发展方向。

通流特性；端区黏性；径向掺混；流线曲率；压气机；航空发动机

0 引言

高负荷多级压气机是航空发动机关键部件之一，其气动设计重点在于如何快速、完整、准确地获得全工况下各排叶片前尾缘由叶根到叶尖气流参数和整机性能变化，并以此为据进行设计调整，其预测质量直接决定多级压气机的设计成败。近20年来，3维定常/非定常流动模拟技术日趋成熟，并在设计中得到充分应用，然而，除了物理模型仍然存在严重失真外，为实现上述目的而导致的海量计算、海量数据信息还使得设计周期过长、设计导向性很不明朗，而通流计算不但能够得到平均流面上详细的流场信息，而且计算时间比3维CFD的要少很多，适用于在设计阶段快速分析压气机的性能，评定压气机的设计效果。相比之下，通流特性计算虽有模型经验性强的缺点，但由于其计算速度快以及设计指向性明确，在经历CFD发展带来的冲击、反省后，始终作为多级压气机设计技术核心而发展[1-4]。通流计算方法基于吴仲华的S1/S2 2元流面理论[5-6]，在20世纪50～70年代获得长足发展，主要体现在有关损失、脱轨角、堵塞模型等大量经验积累和总结，并由Novak等首先将其引入工程应用[7]；1990年前后，Gallimore/Adkins等的研究考虑了展向掺混和端壁黏性；Casey/Robinson等则在流线曲率算法改进方面取得重要进展[8]；国内自1990年引进国外通流特性计算程序，很好地支撑了设计工作，但该程序适用负荷范围小、未考虑掺混和端壁黏性，制约了国内多级压气机技术发展。

本文面向高负荷组合压气机特性计算应用，推导了统一适用于轴流、斜流、离心压气机的流线曲率方程，引入端区黏性和径向掺混模型，针对轴流压气机全新开发了其通流特性计算程序，并以3个算例进行初步验证。

1 流线曲率通流特性计算方法

面向高负荷组合压气机特性计算应用，推导了统一适用于轴流、斜流、离心压气机的流线曲率方法，其最终形式为

其中

式中：m为流线方向；n为拟法线方向；s为熵；I为转焓；I*为总焓；Mam为子午马赫数；Maθ为切向马赫数。

式中：τ为切应力；ρ为气体密度。

在推导上述方程过程中，直接从可压缩黏性N-S方程入手，仅作了周向平均假设，未作等熵假设和无黏简化。该方程的精确解为

式中：Vm为子午速度；Δn为同一计算站上相邻流线间距。

为使方法统一适用于轴流、斜流、离心压气机，重新定义了计算站角度和流通面积。流线曲率法角度如图1所示，α表示计算站和流线方向的夹角。角度的正负作如下规定：γ为负，当位于r轴另一侧为正；φ为正，当位于z轴另一侧为负。流通面积如图2所示，其大小是同一个计算站上相邻网格点的连线绕z轴旋转所形成凸台的侧面积，如此避免了计算离心压气机时由于流线斜率过大而导致的面积计算错误。

图1 流线曲率法角度

图2 流通面积

在具体编程过程中，对每个计算站都要进行内、外2层循环迭代。在内层循环中从中间流线到两边流线迭代求解出合适的子午速度，在外面1层循环中还需要根据质量守恒，通过子午速度和每2条流线之间的流管面积插值求出新的流线位置。以上2层循环的收敛完成了1个计算站的通流计算。

2 经验模型

通流计算依赖大量的经验模型，经验模型的精确与否是通流计算成败的关键。需要用到经验关系的计算主要有损失计算、边界层拥塞和流动分离[11]。

2.1 攻角、落后角的计算方法

程序中攻角、落后角的经验模型借鉴NASA以及俄罗斯的研究成果。

选用NACA-SP36中的根据叶栅吹风试验结果求最小损失攻角和根据气流与叶背相切的条件计算最佳攻角，分别为

式中：（ki）sh、（ki）t为修正系数；（i0）10为相对厚度为10%零弯度叶型叶栅的最佳攻角；c和n是经验指数；θ为叶型弯度；a/b为最大弯度相对位置；Ma为进口相对马赫数。

程序中落后角计算为

式中：（kδ）sh、（kδ） 为修正系数；（δ0）10为相对厚度为10%零弯度叶型叶栅的最小损失落后角；m和σ为经验指数；θ为叶型弯度；a/b为最大弯度相对位置；b/t为稠度。由此计算得到的落后角还需进行攻角、子午速度、速度系数和径向间隙的修正。

2.2 端壁黏性的计算方法

在代码中通过黏性损失和黏性切应力来计算端壁黏性。

在整个子午流场，端壁边界层除了导致流场拥塞，边界层内的黏性效应也是压气机损失的1个重要因素。压气机中的损失通常用损失系数来表征，常用压力损失系数和基于熵变的损失系数，而二者可以通过数学变换相互转化[12]。在靠近壁面的网格处考虑黏性耗散导致的损失。定义黏性耗散系数Φ为

耗散系数表征了由黏性导致的机械能向热能的不可逆转化。在实际使用中，可以通过如下经验关系计算耗散系数[13]

式中：V为边界层内的平均速度。通过上面的计算可以得到边界层内黏性效应导致的熵增，进而计算出损失系数。

程序中的黏性切应力模型参考Denton的黏性边界层算法。对于方程中的黏性项使用Denton的边界层内切应力为

通过该模型可以计算出黏性切应力，对应通流方程中保留的黏性项，将该模型引入通流方程求解中，计及边界层内的黏性应力对压气机流场产生的影响。

2.3 径向掺混方法

以往的通流代码最大的缺陷之一就是没有考虑流动沿展向的掺混效应。在实际情况下，压气机内的流动具有强烈的3维特性。由于流面扭曲、二次流动、尾迹以及湍流扩散等的影响导致了平均流线的展向输运。在上文介绍的通流计算核心算法中，只能通过迭代求解出气体沿着流线流动的子午速度，无法计及气体沿着展向的掺混，因此，必须通过额外的经验模型来考虑展向掺混。在实际流动中，沿展向输运的参数有角动量、转焓和熵[11]，分别表征了动量、能量和损失的展向输运。本程序中采用Denton提出的掺混模型来处理这些参数的展向输运。径向掺混关系如图3所示。根据Denton的掺混模型[14]，掺混参数的关系为

图3 径向掺混

式中：P可以为角动量、转焓或熵。

先通过不考虑径向掺混的通流代码计算得到流场信息，然后按照式（13）对角动量、转焓和熵进行修正，得到新的流场信息作为下次迭代的起点。

该方法存在的问题是，当流线数目变化时，公式中的系数f也需修正。为了避免该问题，Denton建议取为0.5。

3 结果校验

选取Rotor67转子、Rotor37转子和GY1-2J压气机3个算例对程序的正确性进行验证，同时对比程序中加入掺混、端区黏性模型对计算结果的影响。

3.1 Rotor67转子算例

Rotor67转子是1个低展弦比、无阻尼环设计的某2级风扇的第1级转子[15]，设计转速为16043 r/s，叶尖速度为429 m/s，相对马赫数为1.38，转子有22个叶片，展弦比为1.56。Rotor67转子在100%转速下的压力和效率特性曲线如图4所示，对比了加入掺混、黏性模型、未加入新模型以及试验结果。在设计转速、30.5 kg/s流量下代码改进前、后子午流场的对比结果如图5所示。

图4 Rotor6转子在100%转速下的特性

图5 在100%转速、30.5 kg/s流量下的子午速度分布

根据上述结果可知，目前的代码对于跨声速转子的计算还不够精准，主要体现在压比计算结果比试验结果约有6%的误差，说明经验模型的适用范围还需进一步调研和优化。对于跨声速转子，模型中的经验系数还需通过大量计算来进一步标定。但是计算所得的压气机特性趋势与试验结果一致，模拟具有一定的精度，这对于后续将该程序延伸至风扇的特性计算具有参考意义。对比加入黏性损失模型和掺混模型前、后的流场结果，其区别主要体现在叶根处，新模型的加入使得该区域结果更加细腻，速度梯度更加明显，主要体现了黏性应力模型在通流代码中的作用。3.2 Rotor37转子算例

Roto37转子的设计转速为17188.7 r/s，叶尖速度为454.14 m/s，转子有36个叶片，设计压比为2.106，设计流量为20.19 kg/s。Rotor37转子在100%转速下的压力和效率特性曲线如（a）效率特性（b）压力特性图6所示，对比了加入掺混、黏性模型、未加入新模型以及试验结果。在设计转速、19.42 kg/s流量下代码改进前、后子午流场的对比结果如图7所示。

图6 Rotor37转子在100%转速下的特性

图7 在100%转速、19.42 kg/s流量下的子午速度分布

从Rotor37转子的流场对比结果较明显可见掺混模型的影响。展向掺混使得整个流场沿径向更趋向于均匀。特性预测方面，新模型的加入使得压力特性线更趋近于试验结果，而效率预测的精度相差无几。上述2个算例显示，在特性预测方面，黏性和掺混模型的加入对结果产生的影响均不明显，主要原因是对于单级转子，黏性效应对压比和效率的影响有限，黏性耗散产生的损失相对于叶型损失、二次损失等较小。展向掺混最主要的影响体现在对流场结构的改变，虽然掺混会改变熵量，但其改变也相对较小，对于单排转子表现并不明显。与试验结果对比，压比特性线的误差都比较明显，可能是因分离角的计算存在误差造成的。分离角经验模型中除了用到叶型几何参数，还与马赫数等流场信息相关，上述2个算例的出口马赫数有明显不同，这可能是导致1个算例的压比较小，而另1个算例较大的原因。

3.3 GY1-2J压气机算例

对GY1-2J压气机的模拟计算包括进口导叶和2级叶片。在80%转速、13 kg/s流量下算法改进前、后流场的计算结果对比如图8所示；在设计转速、17 kg/s流量下的流场计算结果如图9所示；该压气机特性线的计算结果对比如图10所示；算法改进前、后第2级转叶在70%转速、9 kg/s流量下进口气流角的对比如图11所示；第2级转叶在100%转速、17.84 kg/s流量下进口气流角的对比如图12所示。

图8 在80%转速、13.0 kg/s流量下的子午速度分布

图9 在100%转速、17.0 kg/s流量下的子午速度分布

图10 GY1-2J压气机特性计算结果

通过对子午流场的对比计算可知，黏性和掺混的加入使得计算结果变化比较明显，整个子午流场的速度分布比未加入黏性和掺混模型前更均匀。通过特性线的对比可知，在非设计工况下二者计算结果相差比较明显，而在设计工况下比较接近。原因之一是在非设计工况下流动的径向掺混效应比设计工况下的更为明显。从图11、12中也可以得到这一结论。在非设计工况下，新模型的加入使得进口气流角沿着展向变化显著，而在设计工况下则变化不明显。同时可见，由于展向掺混和端区黏性产生的额外损失并不明显，因此在效率和压比特性线上并无显著差别。从子午流场中还可见，黏性模型的加入改变了端区的速度分布，使得速度梯度在该区域更加明显。

图11 在70%转速、9 kg/s流量下的进口气流角对比

4 结束语

图12 在100%转速、17.84 kg/s流量下的进口气流角对比

本文回顾了流线曲率法的主要内容，推导了统一适用于轴流、斜流、离心压气机的流线曲率方程，全新开发了1套模块化、可维护的通流计算代码框架，并在框架内加入了展向掺混和端壁黏性模型。通过对经验模型改进前、后计算结果以及与试验结果的对比，初步校验代码框架的准确性，同时也显示出了模型的通用性还存在局限，如何提高代码的适用范围是后续研究方向。

[1]程荣辉.轴流压气机设计技术发展 [J].燃气涡轮试验与研究，2004，17（2）：1-8.CHENG Ronghui.Review on development of design technique of axial compressors[J].Gas Turbine Experiment and Re search，2004，17（2）：1-8.（in Chinese）

[2]季路成.高性能叶轮机全3维叶片技术趋势展望[J].航空发动机，2013，39（4）：9-18.JILucheng.Trend of full three-dimensionalblading techniques for high performance turbomachinery[J].Aeroengine，2013，39（4）：9-18.（in Chinese）

[3]周盛.涡轮风扇发动机气动热力学三项新概念探讨[C]//空天安全的重大基础问题学术讨论会文集.北京：香山会议，2001：53-57.ZHOU Sheng.Investigation on three new aerothermodynamic concepts of turbofan engine [C]//Symposium on Fundamental Problems of Space Security.Beijing：Xiangshan Conference,2001:53-57.（in Chinese）

[4]邵卫卫.风扇/轴流压气机最大负荷设计技术探索 [D].北京：中国科学院研究生院，2008.SHAO Weiwei.Exploration on maximum loading design tech niques for fan and axial-flow compressor[D].Beijing：Graduate University of Chinese Academy of Sciences，2008.（in Chinese）

[5]Wu C H.A general through flow theory of three dimensional flow in subsonic and supersonic turbomachines of axial,radial andmixed-flow types[R].NACA-TN-2604,1952.

[6]季路成，田勇，李伟伟，等.叶身/端壁融合技术研究[J].航空发动机，2012，38（6）：5-11.JI Lucheng，TIAN Yong，LI Weiwei，et al.Investigation on blended blade and endwall technique[J].Aeroengine，2012，38（6）：5-11.（in Chinese）

[7]Novak R A.Streamline curvature computing procedures for fluid flow problems[J].Journal of Engineering for Power，1967，89（4）：478-490.

[8]Denton JD.Throughflow calculations for axial flow turbines[J].Journal of Engineering for Power，1978，100（2）：212-218.

[9]Adkins G G，Smith L H.Spanwisemixing in axial-flow turbomachines[J].Journal of Engineering for Power，1982，104（1）：97-110.

[10]Gallimore S J.Spanwise mixing in multistage axial flow com pressors:Part II-throughflow calculations includingmixing[J].Journal of Turbomachinery，1986，108（1）：10-16.

[11]Casey M，Robinson C.A new stream line curvature through flow method for radial turbomachinery[J].Journal of Turbomachinery，2010，132（3）：21-31.

[12]Brown L E.Axial flow compressor and turbine loss coefficients:a comparison of several parameters[J].Journal of Engineering for Gas Turbines and Power，1972，94 （3）：193-201.

[13]Hall D K，Greitzer E M，Tan C S.Performance limits of axial compressor stages[R].ASME 2012-GT-69709.

[14]Denton JD，Hirsch C.Through flow calculations in axial turbomachines[R].AGARD-AR-175,1981.

[15]Strazisar A J.Laser anemometermeasurements in a transonic axial-flow fan rotor[R].NASA-TP-2879,1989.

Stream line Curvature Algorithm Considering Viscosity and RadialMixing on Tip Region

LIU Zhi-wei1,JILu-cheng2,CAIJun1
（1.Institute of Engineering Thermophysics,Chinese Academy of Sciences,Beijing 100190,China;2.Schoolof Aerospace Engineering,Beijing Institute of Technology,Beijing 100081,China）

One technical difficulty of high loadingmulti-stage compressor is how to get the distribution of aerodynamic parameters along blade leading and trailing edge underalloperation conditionsand thewhole engine characteristics quickly,completely and accurately.The through flow calculation is the core technology in multi-stage compressor design based on its quick calculating speed and explicit orientation.In order to improve the applicable scope and accuracy of the code,a uniform streamline curvature equation for both axial and centrifugal compressorwasderived.The through flow performance calculation codeofaxial compressorwasnewly developed by importing the tip region viscosity and radial nuxingmodel,and applying the Fortran language.The performance of NASA Rotor 37,Rotor 67 and GY1-2J compressorwas calculated,and the results show that calculation precision of the code is improved by the newmodel,and the effectiveness of the code is validated.To improve the universality of the code is the developmentdirection of through flow calculation in the future.

through flow calculation;tip region viscosity;radialmixing;stream line curvature;compressor;aeroengine

V211.6

A

10.13477/j.cnki.aeroengine.2014.03.007

2013-03-18 基金项目：国家自然科学基金（51176012、51006100）资助

刘志伟（1987），男，在读硕士研究生，研究方向为叶轮机械气动热力学和通流计算；E-mail：zhiwei19871987@163.com。

刘志伟，季路成，蔡军.考虑端区黏性和径向掺混的流线曲率通流算法研究[J].航空发动机，2014,40（3）：34-39.LIU Zhiwei,JILucheng,CAI Jun.Stream line curvaturealgorithm considering viscosityand radialmixingon tip region[J].Aeroengine,2014,40（3）：34-39.