基于时间推进的通流计算方法：现状及展望

杨金广， 王春雪， 王大磊， 邵伏永， 杨晨， 吴虎

1.大连理工大学 能源与动力学院， 大连 116024 2.中国航天科工集团 北京动力机械研究所， 北京 100074 3.西北工业大学 动力与能源学院， 西安 710072

基于时间推进的通流计算方法：现状及展望

杨金广1,*， 王春雪2， 王大磊2， 邵伏永2， 杨晨3， 吴虎3

1.大连理工大学 能源与动力学院， 大连 116024 2.中国航天科工集团 北京动力机械研究所， 北京 100074 3.西北工业大学 动力与能源学院， 西安 710072

为探究基于时间推进的通流方法在叶轮机械设计分析中的应用潜力，报告了国内外相关研究现状，详细阐述了该方法在应用过程中的若干关键问题，包括：叶片堵塞的几何描述、叶片力的模拟、叶片前后缘间断问题的处理和激波的捕获等，并对不同模型进行了比较。总结可知，相比传统通流方法，该方法具有捕捉激波，适应堵塞，流场模拟更为精细以及计算扩展性良好等诸多优点，在先进航空发动机或燃气轮机的部件设计和分析、整机稳态、过渡态性能模拟中，都具有相当的优势与应用潜力；同时尽管该方法的研究取得了很大进展，但在上述若干关键问题上仍有提高空间。通过进一步的发展与完善，该方法有望成为叶轮机械设计和分析的标准工具。

叶轮机械； 通流； 时间推进； 叶片堵塞； 叶片力； 叶片前后缘间断； 激波

吴仲华[1]于1952年提出三元流动理论，将叶轮机械内部全三维非定常黏性流动近似为两类流面内的准三维定常流动，即S1流面(叶片截面)和S2流面(子午面)，并在两类流面间迭代求解。三元流动理论提出后被广泛地应用于叶轮机械的设计与性能预测中，极大地促进了相关行业和学科的发展。而其中又以在S2流面上无叶和有叶区域内都进行计算的通流计算方法应用最为广泛，在风扇、压气机以及涡轮等叶轮机械的设计体系中占据着重要地位。尽管全三维定常和非定常黏性计算迅速发展并且开始应用于叶轮机的设计中，发展先进的通流计算方法仍然非常有必要。实际上，以20世纪90年代以来美、英等国的发动机公司的风扇、压气机设计体系为例，设计过程大致可以分为5个密切相关的步骤：初步设计、通流计算、叶片造型(二元气动优化造型)、叶片造型(三元气动优化造型)和放大尺寸的实验研究。这5个步骤环环相扣，每个阶段采用不同层次的数学、物理模型和经验数据，相互补充，相互交叉检验，最终将设计风险降到最小[2]。相较于三维的黏性计算，通流计算所需要的时间要降低几个数量级，尤其在多级计算中这一优势更为明显。因此设计中可以将通流计算和三维计算作为互补的工具，后者可以在前者的计算结果基础上进一步验证；同时通流方法与三维计算相耦合，提供更加灵活的精度和计算周期选择。此外通流方法在航空发动机/燃气轮机的整机模拟中具有很大的潜力，2015年Petrovic和Wiedermann[3]采用基于流函数的通流方法，结合简化的燃烧室模型，对整个燃气轮机在S2流面上成功地进行了性能模拟。

流行的通流计算方法有2种，Novak[4]提出的求解完全径向平衡方程的流线曲率法，以及Marsh[5]提出的引入流函数的矩阵通流方法。这些方法具有快速稳定的优点，并且发展时间长，理论较为完备，在工程实际中得到了广泛的应用。然而这2种方法均存在一个显著的缺点，即由于激波的存在以及堵塞工况的出现，在跨、超声速流动中的应用受到限制。因此在跨、超声速风扇、压气机和涡轮等叶轮机械设计中，需要一种实用的能够处理高速流动的通流方法，基于时间推进求解Euler或Navier-Stokes方程的通流方法是就是一种非常合适的选择。相比于流函数法或者流线曲率法，基于时间推进的通流方法——计算流体力学(Computational Fluid Dynamics，CFD)通流方法具有以下优点：① 质量流量由流场计算本身得出，因而可以预测叶片排的堵塞工况；② 可以自动捕获特定形式的激波，从而进行跨、超声速流场计算；③ 可以应用CFD领域发展成果，如有限体积方法、有限元方法以及高精度数值格式等，同时时间推进的引入使得CFD通流方法可以自然地应用于过渡态的计算。

本文对国内外基于时间推进的通流计算方法进行总结，对该方法在具体实施中的关键问题处理方法进行讨论，最后对其未来的发展方向进行展望。

1 国内外研究进展

在过去的几十年里，基于时间推进的通流计算方法有了一定程度的发展。对于国外而言，Spurr[6]最早于1980年在通流计算中采用时间推进的方法对周向平均的Euler方程进行了求解。英国剑桥大学Whittle实验室的Dawes[7]对三维Navier-Stokes方程求解器进行了修改，使其能够进行Euler通流模型的求解，并以多级压气机为例，使用CFD通流方法提供多级环境，仅对单排叶片进行三维计算，对CFD通流方法与三维黏性数值计算的结合进行了探索性研究。比利时布鲁塞尔自由大学的Yao和Hirsch[8]基于三维求解器修改实现了CFD通流方法计算，并细致分析了叶片区弦向气流角和堵塞分布的影响。美国雪城大学的Damle等[9-10]发展了求解Euler方程的完全通流求解器，并实现了对跨声速和超声速叶片的设计，展示了其在跨、超声速设计问题中的潜力。瑞典查尔姆斯理工大学的Baralon等[11-13]对CFD通流方法进行了深入的研究，提出了一种基于伪时间步长的叶片力计算模型，并针对叶片前缘的间断问题提出了2种处理方法，使得CFD通流方法的求解稳定性和实用性得到了提升；Baralon等同时对Euler通流模型的激波捕获机理进行了研究，提出了一种计算流面法向堵塞因子的新方法，并在Euler通流模型中引入了叶型损失、落后角和端壁损失，此外还引入了径向掺混模型。Sturmayr和Hirsch[14]对通道激波的捕获机理进行了研究，并对设计问题和分析问题进行了对比。进入21世纪，CFD通流方法研究进入高潮，比利时列日大学的Simon等[15-17]发展了求解周向平均Navier-Stokes方程的通流求解器，采用2个亚声速算例进行了验证，并与Euler通流模型进行了对比，此外还对叶片力项、周向应力、确定应力等在周向平均通流模型中的相对重要性进行了详细研究。意大利米兰理工大学的Persico和布雷西亚大学的Rebay[18]共同提出了一种基于罚函数的叶片力模型。都灵理工的Taddei和Larocca[19]对CFD通流方法的设计问题进行了持续的研究，并且提出了2种处理叶片前后缘间断问题的新方法[20-21]。佛罗伦萨大学的Pacciani等[22]提出了一种自适应流面修改方法，以解决前后缘间断问题。值得注意的是俄罗斯中央航空发动机研究院(Central Institute of Aviation Motors，CIAM)的Nigmatullin和Ivanov[23]基于高阶精度的Godunov格式，采用贴体坐标系，完全独立地建立了一整套CFD通流计算理论模型，包括叶片力的模拟和叶片前后缘间断问题处理模型等，并且已经成功应用于航空发动机整机流场计算。

相比于国外CFD通流方法的长期发展，国内相关研究起步较晚。袁宁等[24]对流函数法、流线曲率法以及求解Euler方程的CFD通流方法3种S2流面计算程序进行了比较。季路成等[25]在国内较早地将时间推进思想应用于通流计算，并设计了一轴向流动跨声速的激波转子进行验证。施发树和刘兴洲[26-27]与CIAM联合开发了涡轮发动机整机通流计算程序，并对一台小型双涵道涡扇发动机进行了整机S2计算。清华的于龙江[28]、北航的曹志鹏[29]和金东海[30]等使用CFD通流方法进行了整机流场仿真计算。李德英等[31]基于CFD通流方法开展了多级涡轮的气动分析，并得到了与三维计算较为一致的结果。万科等[32-33]基于周向平均的Navier-Stokes方程通流计算方法考察了周向平均方法在风扇/增压级分析中的准确性，并对周向不均匀项对于周向平均模型的影响进行了研究。

尽管CFD通流方法已有几十年的发展时间，理论上相比传统的通流方法也更具优势，但要在工程实践中成为一种成熟的设计手段仍然面临诸多困难。本文将就各国研究者在CFD通流方法研究方面所做的工作进行介绍，同时针对一些关键问题展开讨论，分析影响计算准确性的关键问题，提出可能的解决方案，并展望未来的发展趋势。

2 CFD通流方法基本计算模型

通流计算的基本思想在于认为叶轮机中的流动在无叶区满足轴对称假设，在叶片通道内则沿S2流面流动。因此，CFD通流计算模型一般通过对三维的Euler/Navier-Stokes方程在相应方向进行平均得到，从而实现从全三维的流动问题(图1(a))向平均叶片S2流面求解(图1(b))的降维简化[18]。如果此时进一步在流场中施加适当形式的体积力来模拟叶片的作用，迫使当地流动速度方向满足与平均叶片流面相切，则问题可以进一步简化到子午平面内的流场计算问题(图1(c))。

图1 从全三维模型向通流模型的简化[18]Fig.1 Reduction from fully 3D to throughflow models[18]

不同的叶片力模型会决定计算模型是否需要对叶片区周向动量方程进行求解。因此，根据物理计算域的简化程度以及是否在有叶区求解周向动量方程，目前的CFD通流方法可以分为3类：

第1类，在叶片区不需求解周向动量方程(即下文式(1)中第4个方程)，周向动量方程用来推导体积力形式的周向叶片力。该类计算模型相当于在叶片区求解了二维的流场问题。该方法可以同时进行分析问题(正问题)和设计问题(反问题)的求解。求解正问题时，在叶片区每一时间步中通过子午速度计算满足流面相切条件的周向速度，并对叶片力进行更新，计算得到新的子午速度，直到得到流场的收敛解。在反问题中，则首先规定载荷或者环量分布，通过流面相切条件计算得到新的平均叶片流面几何；作为一个初值问题，流面几何的时间推进计算需要设定边界条件，通常通过确定前缘倾斜形状或者积叠规律实现；根据新的流面几何用来更新叶片力，直到得到收敛解。Dawes[7]、Damle[9]和Taddei[20]等均采用了这种计算模型，所不同的是Dawes的通流求解器由三维Navier-Stokes方程求解器修改而来，通过在周向仅布置一个网格单元结合周期边界条件实现CFD通流方法求解；同时Dawes只进行了通流正问题的求解，主要用于为单叶片排三维计算提供多级环境，而Damle和Taddei均同时进行了正、反问题的求解。

第2类，在叶片区求解周向动量方程，通过添加适当形式的体积力，模拟无黏叶片力的强迫作用，使得在最终收敛的流场计算结果中，流体流动方向满足流面相切条件。该种计算模型在子午平面内进行求解，不同于第1类模型。在正问题求解中，该类模型的周向速度完全通过求解得出，可以认为具有准三维的特征；同时叶片力模型的数学性质更有利于计算迭代的稳定和收敛，因此得到了广泛使用。Sturmayr[14]、Baralon[11]、Simon[15]及Persico[18]等均采用了此种CFD通流方法基本计算模型。其中Sturmayr和Hirsch并没有求解真实的叶型几何，而是通过幂函数的形式模拟了中弧线以及堵塞因子的分布，使其符合真实叶型的主要特征，以此探索激波捕捉机理以及不同的参数分布对正问题求解的影响。Baralon[11]和Simon[34]等均尝试了采用法向堵塞因子，希望借此提高流场的计算精度。由于叶片力模型的不同，该类方法在反问题中的使用与第1类模型有所不同。Sturmayr和Hirsch[14]对正问题和反问题进行了比较，在反问题中首先指定叶片尾缘的周向速度，然后根据来流周向速度，在叶片内从前缘到尾缘按幂函数形式分布作为目标周向速度，当计算收敛时，计算得到的周向速度与目标周向速度趋同，此时的流面即设计目标。Simon[34]采用了相同的方法进行反问题的求解。

前两类CFD通流方法均在子午平面内进行流场计算。以周向平均后的轴对称Euler方程为例，一般的控制方程表达形式为

(1)

式中：

其中：U为守恒量；F和G为对流(无黏)通量；S为源项；ρ为密度；p为压力；z、r和φ分别为轴向、径向和周向坐标；u、v和w分别为轴向、径向和周向速度分量；e为总能；ε(T)为温度T时的单位质量内能；b(r,z)为堵塞因子，模拟叶片厚度带来的堵塞效应；FB为无黏叶片力项；fBz、fBr和fBφ分别为FB在3个方向的分量；FL为黏性力项；fLz、fLr和fLφ分别为FL在3个方向的分量；在叶片区以外FB为0；t为时间；ω为角速度。

第3类，直接在三维平均叶片流面上求解。目前的公开文献中，仅有Nigmatullin和Ivanov[23]通过对相对圆柱坐标系下的三维控制方程进行空间坐标平均，建立了基于贴体坐标系(ξ,η,ζ)的CFD通流方法基本计算模型。该计算模型在叶片区和无叶区采用了不同的控制方程，在无叶区满足轴对称条件，在叶片区对每个流面网格单元求解流面法向平均的控制方程。相比于前两类计算模型在子午平面进行通量的求解，该计算模型通量的求解无疑带有更多的三维特征，同时该模型的叶片力模拟也更为方便，通过理论推导即可得出。该叶片力模型能够与流场的迭代同步进行。该模型对于叶片前后缘间断问题采用了任意间断分解方法进行处理，较为完整地考虑了流动中可能出现的各种情况，同时鲁棒性较好，是非常具有应用潜力的CFD通流方法。

3 CFD通流方法中的若干关键问题

CFD通流方法在应用中面临诸多问题，如叶片堵塞的描述、叶片力的模拟、叶片前后缘间断问题的处理以及激波的捕获等，只有处理好这些问题，CFD通流方法才可能作为一种快速且准确的设计分析手段应用于叶轮机械领域。下文将对这些问题以及目前常见的处理方法进行详述。

3.1 叶片堵塞的几何描述

通过对三维控制方程进行平均，CFD通流方法实现了流场计算的降维，但与此同时，在叶片区需要使用堵塞因子对真实叶片厚度造成的流道堵塞效应进行模拟。堵塞因子的定义一般分为2种，一种为周向堵塞因子，另一种为法向堵塞因子。周向堵塞因子通过对三维Navier-Stokes方程进行周向平均后直接导出，定义为

(2)

式中：φp-φs为叶片在周向以弧度计量的厚度；Nb为叶片数；π为圆周率。显然在无叶区堵塞因子为1。Sturmayr和Hirsch[14]采用了周向堵塞因子进行计算，但没有使用真实的叶型几何，而是近似计算叶型法向厚度最大点对应的周向厚度，即

dmax≈dnmax/cosβnmax

(3)

式中：dnmax为叶型法向最大厚度；βnmax为法向厚度最大位置点对应的叶片角。以dmax为基准，沿流向按幂函数分布确定叶片周向厚度，并在叶片前后缘附近对堵塞因子分布进行光滑。通过试取不同的幂指数和调整dmax的弦向位置实现不同的堵塞因子分布。基准算例幂指数为2，dmax位于0.5cc处。他们研究了流场计算结果对于堵塞因子分布的敏感度，如表1[14],其中:c为弦长;Case 1 为基准算例。由表1的计算结果发现，不同堵塞因子分布对堵塞流量以及对应效率的计算结果有明显影响。

为了更好地估计实际叶片厚度对流动的影响，Baralon等[11]认为采用法向堵塞因子比周向堵塞因子更为合适，并给出了相应的法向堵塞因子计算方法，表达式为

(4)

表1 堵塞流量对堵塞因子分布的敏感度[14]

Table 1 Sensitivity of blocking mass flow to blockage

factor distribution[14]

CaseModificationMassflow/(kg·s-1)PressureratioEfficiency/%1dmaxat0.5c30.81.53782dmaxat0.75c32.71.52863dmaxat0.25c27.21.53664Nosmoothing28.91.53725dmax-20%32.21.54826dmax+20%29.31.51757dmaxexponent1.531.81.53818dmaxexponent329.11.5275

式中：S为栅距；β为平均流面角；β1和β2为叶片角；t1和t2为相应的叶片厚度。该表达式仅限于计算点A和点B之间的堵塞因子，如图2[11]所示，从上方叶片前缘点向平均流面(在S1面上为平均流线)作垂线，垂线与平均流面交点即为A点；同理，点B由下方叶片尾缘点向平均流面作垂线得到；点P为点A和点B之间的平均流线上任一点。对于前缘到点A和尾缘到点B的部分，堵塞因子定义为从bnA和bnB线性变化到1。可以发现，该方法定义的法向堵塞因子与周向堵塞因子差别较为明显，对于某一展向叶片截面，计算得到的流道面积大小和喉道轴向位置均不相同，这些差异会对流场计算结果造成较为明显的影响，尤其在激波存在的跨、超声速流动中。该种堵塞因子定义方法存在实际操作复杂、前后缘部分分布过于简单、精确度较差的缺点。

Simon[34]借鉴Baralon的思路，提出了一种简化的法向堵塞因子计算模型

bn=bcosβ

(5)

式中：b为周向堵塞因子，cosβ为对应叶片角余弦值。Simon使用Rotor67转子作为算例，采用2种堵塞因子分别进行了计算，并对得到的结果进行了比较，发现采用法向堵塞因子的堵点流量计算结果与实验值符合得更好，并且激波更为靠前。

图2 法向堵塞因子[11]Fig.2 Normal blockage factor[11]

Nigmatullin和Ivanov[23]的计算模型将周向堵塞因子与贴体坐标系和圆柱坐标系之间的变换进行了关联，定义b=1/ζφ，应用于度量系数的计算当中，从而模拟叶片堵塞效应。考虑到ζ=const对应于平均S2流面，因此该定义可以作为法向堵塞因子的一种标准定义。

从流动的物理特征上理解，基于周向平均的堵塞因子定义不能真切反应叶片真实存在造成的物理堵塞效应，而在这方面法向堵塞因子更具优势，但其计算方法仍有待改进。另外堵塞因子的定义要与控制方程相匹配，如果采用周向平均堵塞因子，则可以采用周向堵塞因子；如果沿平均流面法向进行平均，则采用法向堵塞因子更好，这也是很多文献中忽略的问题。总之，目前对于堵塞因子的研究仍然不够系统和全面，堵塞因子的定义以及堵塞效应的数学描述需要从理论上更进一步完善。

3.2 叶片力

通流计算中，除了需要考虑实际叶片厚度带来的影响，另一个主要问题便是如何在叶片区模拟叶片对流体的作用力。Marble[35]最早提出采用分布的体积力模拟叶轮机流场，这种理念也被拓展应用于通流计算中，用分布的体积力模拟叶片力[36-37]。通常将叶片力分为无黏叶片力和黏性力两部分进行模拟，从而达到简化的目的。对于无黏叶片力，其作用方向通常定义为叶片平均流面法向，即

FBW=uρfBz+vρfBr+wρfBφ=0

(6)

式中：FB为无黏叶片力；W为相对速度矢量；定义平均流面α=φ-g(r,z)，g为r和z的函数，流面单位法向量

单位法向量

则FB可表示为

FB=B(r,z)·ng

(7)

叶片力的轴向和径向分量为

由此FB的方向已经确定，B(r,z)为叶片力。

一般存在2种思路去模拟无黏叶片力，一种通过周向动量方程直接计算得到无黏叶片力周向分量ρfBφ，如式(8)所示，通过投影即可得到B(r,z)及各分量的大小，Dawes[7]和Damle[9]等均采用了此种方法。

(8)

式中：m为子午方向；qm为质量流量；Vφ为周向速度。

另一种思路为将叶片力看做时间相关的未知量，通过求解时间相关的微分方程对无黏叶片力进行计算，相比于求解式(6)的刚性约束，该方法数学性质上更易于稳定求解。该方法基本思想在于周向叶片力的表达必须满足流面相切条件，基于此Baralon等[11]提出的一种采用伪时间步长的叶片力计算模型，时间相关的微分方程定义为

(9)

式中：λ为典型响应时间平方的倒数，s-2。ρfBφ的响应时间应该和流场计算中的当地网格时间步长处于同一量级，因此λ可以表达为

式中：Δx为典型网格尺寸；Δt为当地网格时间步长。该方程可以理解为一种压力梯度，一旦流面相切条件不符合就会产生，该方法相当于对周向叶片力在时间上的一种离散。Simon和Léonard[15]采用了类似的叶片力模拟方法，只不过采用叶片力的模作为求解变量。万科等[33]指出该种叶片力模型在实际使用过程中，λ的取值对迭代收敛的影响较大，取值不当会影响收敛速度，严重时甚至会导致迭代不收敛。

(10)

式中：C为松弛因子。

Persico和Rebay[18]同样基于流面相切条件，提出了一种基于罚函数的叶片力模型为

B(r,z)=-KW·ng

(11)

式中：K为惩罚因子。该模型可以理解为将流面相切条件施加于有弹性的叶片上，叶片的刚度即惩罚因子K。在实践中发现当K≥104时，最终的计算结果中流面相切条件的偏差不大于0.1°。同时采用上述模型的好处在于不需要求解额外的微分方程。此外可以在计算进程中调节刚度K，在开始迭代阶段K取小值，从而容许流面相切条件有较大的偏差，防止求解的发散，随着计算的进行，刚度K不断增大，从而收敛到满意的解。该模型的根本目的在于将式(6)的刚性约束转化为柔性约束，利于求解。

不同于上述2种思路，Nigmatullin和Ivanov[23]采用了另外一种显式的叶片力计算方法。他们首先计算出轴向、径向和周向动量方程中不考虑无黏叶片力项的残差，对贴体坐标系下离散后的动量方程进行理论推导，并利用流面相切条件，可将无黏叶片力通过计算出的残差项进行表达，这样随着流场计算的收敛，叶片力的分布也一同收敛。叶片力具体表达式为

(12)

式中：J为坐标系变换的雅克比行列式；ζz、ζr和ζφ为与ζ=const网格面对应的轴向、径向和周向的度量系数；S2、S3和S4分别为轴向、径向和周向动量方程中不考虑无黏叶片力项的残差。该方法简单通用，不需要引入伪时间项，亦避免了因此带来的常系数取值问题，同时具有较好的收敛性。

对于黏性力，一般可以写成如式(13)所示的表达式。

(13)

式中：黏性力方向与当地速度方向相反。L(r,z)为一正值，代表黏性损失的大小。对于黏性损失的计算，Bosman和Marsh[38]提出通过给定流场的熵分布进行计算，该损失模型简单且使用广泛。由此L(r,z)可以和熵分布联系起来，即

(14)

式中：T为当地静温；s为熵。

3.3 叶片前、后缘间断问题

叶片前缘冲角以及尾缘落后角在CFD通流计算中是需要被重点关注的问题。由于采用了轴对称假设，CFD通流计算方法无法在无叶区感受到叶片存在而提前产生旋流速度；而一旦进入叶片区，流动则必须满足流面相切条件，周向速度会被平均叶片流面角度所约束，因此在叶片前缘网格界面会产生速度的间断。如果采用叶片力源项模拟前缘的速度转折，就会产生体积力的突变，即便不断加密网格也无济于事。在实践中发现，只要前缘来流冲角大于2°，就会产生比较明显的数值损失[11]，而这种损失所对应的熵增是非物理的。对于叶片尾缘，同样会产生类似的间断问题。

上述问题在很长一段时间内没有得到足够的重视，随着CFD通流方法的发展，Baralon等[11]较早关注和深入研究了该问题，并给出了2种处理办法。第1种通过在流场计算过程中人为修改叶片前缘部分流面，使平均叶片流面角度与来流角度相吻合，冲角变为零，这种方法相当于引入了长度尺度，从而消除了间断问题。对于尾缘可以采用类似的方法，根据输入的落后角大小对尾缘部分流面进行修改。这种方法在过去的十几年得到了较为广泛的应用，包括Simon[34]和Persico[18]等学者均采用了此种前后缘处理方法。图3给出了修改叶片前缘流面和不进行处理计算所产生的熵增对比[34]，其中l.e.(leading edge)指叶片前缘。

尽管简单明了且应用广泛，此种方法仍存在明显的缺点：首先，在冲角较大时，对平均叶片流面的修改程度较大，对于流场计算结果的影响会比较大；其次，在超声速流场计算时对激波的捕获会产生明显影响。由于前缘流面的修改会导致叶片流道面积的变化，对激波是否产生以及激波位置均有决定性影响，进一步对流场的速度分布以及损失预估的精确度产生重要影响。为此，Baralon等[11]同时给出了第2种方法：动量间断处理。这种方法通过在前后缘两侧建立质量守恒、能量守恒关系，同时在动量方程中引入另外的目标间断力项。所谓“间断力”用来控制前缘下游气流角的间断以及所产生的熵增。该方法主要分为2个求解步骤：第1步求解目标间断力，通过引入流面相切、前缘展向连续以及给定前缘总压损失3个条件，与上述守恒关系联立求解得到；第2步即在间断力已知的情况下求解两侧流动参数，并对前缘网格界面通量进行修改。该方法通过间断关系式对前、后缘问题进行了解决，但在实际应用中仍然存在问题，即流动参数的迭代求解在数学上存在困难，尤其在马赫数接近1时，因此无法应用于跨声速流动的情况，如跨声速转子中。此外该方法还需要人为输入额外的参数，即总压损失。

基于对前后缘流面进行修改的思路，Taddei和Larocca[20]提出了采用反问题修正前缘流面的方法。相比于人为地直接对平均叶片流面进行修改，该方法通过给定来流到前缘下游部分的旋流速度分布，求解出新的前缘部分流面。由于直接对速度分布进行操作，该方法对于前缘流场的控制更好，弥补了Baralon第1种方法的一些缺点，但仍然带有一定的任意性，并且在来流冲角较大且叶片扭曲程度较强时容易计算发散。为此Taddei和Larocca[21]提出了一种基于激盘AD

图3 叶片前缘的修改[34]Fig.3 Adaptation of leading edge [34]

(Actuator Disk)模型的前后缘间断问题处理方法。通过在叶片前缘或尾缘采用如图4[21]所示的激盘模型，使得前缘或尾缘网格线i+1/2两侧保证质量守恒、能量守恒和熵守恒，同时集成了冲角或者落后角关系。其中：x为空间轴;t为时间轴;Vx为x方向速度；a为当地声速;i为网格标号。通过对激盘模型的求解，得到1、2两侧的质量、能量以及间断的动量通量，用于上下游流场的求解。由于激盘模型的建立包含了非堵塞工况的假设，因此该前后缘间断处理不能适用于堵塞工况。同时由于没有包含激波关系式，该模型同样不能对激波进行描述。

与Taddei和Larocca[20]的做法类似，Pacciani等[22]对叶片前后缘间断问题进行了研究。他们基于对平均叶片流面进行修改的思想，提出一种时间相关的适应性公式对流面进行更新，如式(15)所示。

(15)

图4 前后缘的激盘模型[21]Fig.4 Model of AD at leading or trailing edge[21]

Nigmatullin和Ivanov[23]同样较早关注这个问题，并对其进行了系统且深入的研究，提出了叶片前后缘任意间断分解的处理模型，较好地解决了该问题。如图5和图6所示[23]，由于在无叶区(轴对称区)和叶片区采用了不同的控制方程，因此前后缘网格界面的法线方向在轴对称侧和叶片侧并不相同。前后缘网格线在S1流面上可分解为CA、CB。图5和图6中：x为CA面的法向量;ξ为CB面法向量；N为平均流面的法向量，与ξ相互垂直；M为AB的法向量；l为AB的切向量。在前缘，2处和1处分别表示轴对称侧和叶片侧，在尾缘则相反。2处参数与L处参数，R处参数与1处参数均通过精确的波特征关系相联系；在前缘，L处与R处通过质量守恒、能量守恒以及熵守恒关系关联；在尾缘处则需要多补充一个方程，一般为AB向的动量方程。

图5 前缘任意间断模型[23]Fig.5 Arbitrary discontinuity model at leading edge[23]

图6 尾缘任意间断模型[23]Fig.6 Arbitrary discontinuity model at trailing edge[23]

该前后缘间断处理模型实际上相当于对一般黎曼问题中的中间波进行修改，这里的中间波不再是接触间断，而是在L和R之间考虑气流转折和流动图景的一种“中间”波，在正常工况下具有熵守恒的性质。同时Nigmatullin和Ivanov[23]还给出了前缘存在脱体激波、前缘堵塞、尾缘堵塞以及出现回流等各种特殊情况的处理方法，使得在任意流态下前后缘间断问题均能被较好地处理。

目前在CFD通流方法前后缘间断问题的处理中，Baralon等[11]的平均叶片流面修改方法虽然存在一些缺点，但由于较为简单和易于实现，得到了较为广泛的应用。Nigmatullin和Ivanov[23]提出的前后缘任意间断分解模型是对前后缘间断问题在物理流动图景下的一种解析解，可以适用于各种工况的计算。Taddei和Larocca[21]的激盘处理可以看做是任意间断分解方法中的正常工况处理。因此任意间断分解法是一种较为完善的前后缘处理模型，但由于较为复杂，采用较少。

3.4 激波的捕获

通常认为CFD通流方法相比于流线曲率法和流函数法最大的优势在于捕获激波的能力。然而，在正问题和反问题中，CFD通流方法表现出来的激波捕获能力以及捕获激波的性质并不相同。Sturmayr和Hirsch[14]采用CFD通流方法对Rotor67分别进行正问题和反问题求解，并对流场马赫数进行了对比，发现反设计得到的流场分布中没有出现激波。Baralon等[11]对通流计算中的激波问题进行了详细分析。结果发现，对于正问题，须给定连续光滑的流动角度分布，此时朗金-雨贡纽关系描述的是流向的正激波关系，因此当流向马赫数大于1时，CFD通流方法能够捕捉到如图7(a)所示的流向正激波。

对于设计问题，需要给定连续光滑的周向速度分布，此时朗金-雨贡纽关系式描述的是轴向正激波，因此CFD通流方法只有在轴向马赫数(考虑径向速度分量时为子午马赫数)大于1时才能捕捉到如图7(b)所示的轴对称正激波。由于在传统叶轮机中，子午马赫数很难达到1以上，因此可以认为在反问题计算中，CFD通流方法不具备自动捕获激波的能力，因此需要添加相应的激波损失模型。

图7 S1流面的流向正激波，轴对称正激波 Fig.7 Normal shock and axisymmetric shock on S1stream surface

3.5 周向应力

在通流计算中，为了得到更加精准的流场模拟结果，需要对流动的三维特征以及非定常效应进行模拟。传统方法均采用经验关系式，适用的普遍性较差。为实现CFD通流方法的精细化预测，减小对经验关系的依赖，必须进一步对CFD通流方法进行完善。Adamczyk[39]针对多级叶轮机非定常流动提出了通道平均理论，通过对非定常Navier-Stokes方程进行一系列平均，利用雷诺平均消除湍流脉动，用时间平均消除非定常影响，用通道平均消除非周期性的影响，依次产生了雷诺应力、确定应力和通道平均应力等附加项。在此基础上进一步进行周向平均，便会产生周向应力项，该项反映了流动周向非均匀性的平均效应。Baralon等[13]以Rotor67为例对周向应力项进行了研究，发现周向应力项对流场的影响大于黏性应力项，产生周向不均匀的主要原因包括间隙泄漏涡、通道激波、激波-边界层干涉以及轮毂角区失速等物理现象。Baralon等还对不同的应力项进行了具体的对比，并在通流计算中引入周向应力项，计算结果表明相比一般的CFD通流方法，考虑周向应力的CFD通流方法与三维结果吻合得更好。Simon等[17]以单级压气机为算例，对不同通流模型进行计算比较，包括只计入无黏叶片力、计入无黏叶片力和黏性应力、计入无黏叶片力和周向应力以及计入全部模型。研究结果表明，在CFD通流计算中无黏叶片力的作用占主要地位，周向应力作用大于黏性应力并且随着负荷的提高而增大。Thomas和Léonard[40-42]首先对一低速压气机级中间展向位置处的S1面流动进行了研究，证实周向应力项具有谐波特征。通过使用扰动场的傅里叶级数对周向应力进行模拟，并以附加项的形式代入到了一维计算中，在一维流场中重现了二维流动特征。Thomas和Léonard进一步将周向应力通过三维扰动场的傅里叶级数进行近似，并通过谐波对近似周向应力进行重构，与CFD通流方法进行耦合求解，成功地反映流动的周向非均匀性，验证了谐波重构的周向应力在CFD通流方法中对流动三维特征的再现能力。Thomas和Léonard还使用非线性谐波结合浸入边界方法对周向应力进行了建模。浸入边界方法用来构建连续的力场以取代真实叶片，解决叶片表面无滑移边界条件与傅里叶级数表达对叶片表面流动连续性要求的冲突。该模型被应用于无黏圆柱算例，并与解析解进行了对比，结果吻合较好。该方法模拟周向不均匀项时，谐波法的级数越高，获得的周向脉动动能与实际状态越接近。

4 结论及展望

CFD通流计算方法可以完成传统流线曲率通流方法的设计和分析功能；同时还具有更多优点：精确捕捉激波，良好的堵塞适应性，能够提供更为精细的流场结构，更加正确的几何描述，也可以更方便地扩展到带有双涵道、引气、冷却以及叶尖间隙等复杂结构的情况。

在正问题中，CFD通流方法对于堵塞因子以及流动角度的分布非常敏感；而对于反问题则没有此类问题；目前的CFD通流方法在传统的叶轮机正问题中可以捕捉流向正激波，在反问题中只能捕获轴向正激波。除此之外，无黏叶片力模型对于CFD通流方法求解得到的流场参数分布具有决定性影响，目前的叶片力模型仍然稍显粗糙，计算结果与三维真实流场相差较大。基于改型黎曼问题的任意间断分解方法为解决前后缘间断问题提供了很好的基础，但在理论上还需进一步验证和完善。上述问题是目前制约着CFD通流方法在工程实际中大规模应用的主要问题。

未来下列问题的研究与解决将是建立高精度CFD通流方法的关键：首先，需要建立法向堵塞因子的标准计算方法，排除人为调整因素，得到准确的流场几何描述，以期更加正确地反应真实几何对流场的物理堵塞效应。其次，通过添加边界层厚度实现对堵塞更为精确的描述，同时需要研究边界层厚度的添加对于叶片力计算的影响；其三，对跨、超声速流动中激波的三维图景进行建模描述，克服CFD通流方法正反问题中分别只能捕获流向、轴对称正激波损失的缺点；其四，将基于改型黎曼问题的前后缘任意间断分解模型与流面修正相结合，实现通流计算中对叶片流动分离工况的处理。此外，以通道平均理论为基础的高阶CFD通流方法也是值得关注的方向。

与传统的流线曲率法类似，影响CFD通流方法计算精度的最主要因素在于所采用的损失和落后角(滑移)模型。相关的研究较多，并且目前尚无统一的高精度模型存在，因此没有在本文中进行讨论。需要指出的是，高精度的通流CFD计算是本文所讨论的堵塞、激波、叶片力、前后缘间断等诸因素与合适的损失和落后角模型全局耦合匹配的结果。

CFD通流方法相对于传统流线曲率通流方法的优点决定了其不仅在传统的叶轮机部件设计和分析中能够占有一席之地，而且在下一代先进航空发动机或燃气轮机的部件设计和分析、整机稳态数值模拟，甚至是过渡态性能模拟中，都具有相当的优势与应用潜力。目前中国在这方面取得了一些进展，但主要还是对世界研究水平的跟随，CFD通流方法的潜力还没有被完全发挥出来。通过进一步的发展与完善，该方法有希望成为叶轮机械以及燃气轮机设计和分析的标准工具。

[1] WU Z H. A general theory of three-dimensional flow in subsonic and supersonic turbomachines of axial-, radial-, and mixed-flow types: NACA TN 2604[R]. Washington， D.C.: NASA, 1952.

[2] 蒋浩兴. 国外发展风扇/压气机设计体系的一些经验和启示[J]. 航空发动机, 2001(2): 45-51.

JIANG H X. The experiences and enlightenments of foreign fan/compressor design system development[J]. Aeroengine, 2001(2): 45-51 (in Chinese).

[3] PETROVIC M V, WIEDERMANN A. Fully coupled through-flow method for industrial gas turbine analysis: GT2015-42111[R]. New York: ASME, 2015.

[4] NOVAK R A. Streamline curvature computing procedures for fluid-flow problems[J]. Journal of Engineering for Power, 1967, 89(4): 478-490.

[5] MARSH H. A digital computer program for the through-flow fluid mechanics in an arbitrary turbomachine using a matrix method: RM 3509[R]. London: National Gas Turbine Establishment, 1968.

[6] SPURR A. The prediction of 3D transonic flow in turbomachinery using a combined throughflow and blade-to-blade time marching method[J]. International Journal of Heat and Fluid Flow, 1980, 2(4): 189-199.

[7] DAWES W N. Toward improved throughflow capability: The use of three-dimensional viscous flow solvers in a multistage environment[J]. Journal of Turbomachinery, 1992, 114(1): 8-17.

[8] YAO Z, HIRSCH C. Throughflow model using 3D Euler or Navier-Stokes solver[J]. VDI Berichte, 1995, 1185: 51.

[9] DAMLE S V, DANG T Q, REDDY D R. Throughflow method for turbomachines applicable for all flow regimes[J]. Journal of Turbomachinery, 1997, 119(2): 256-262.

[10] DAMLE S V. Throughflow method for turbomachines using Euler solvers: AIAA-1996-0010[R]. Reston, VA: AIAA, 1996.

[11] BARALON S, ERIKSSON L E, HÅLL U. Validation of a throughflow time-marching finite-volume solver for transonic compressors: 98-GT-47[R]. New York: ASME, 1998.

[12] BARALON S, HÅLL U, ERIKSSON L E. Viscous throughflow modelling of transonic compressors using a time-marching finite volume solver[C]//the 13th International Symposium on Air Breathing Engines. Reston: AIAA, 1997: 502-510.

[13] BARALON S, ERIKSSON L E, HÅLL U. Evaluation of higher-order terms in the throughflow approximation using 3D Navier-Stokes computations of a transonic compressor rotor: 99-GT-74[R]. New York: ASME, 1999.

[14] STURMAYR A, HIRSCH C. Throughflow model for design and analysis integrated in a three-dimensional Navier-Stokes solver[J]. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 1999, 213(4): 263-273.

[15] SIMON J F, LÉONARD O. A throughflow analysis tool based on the Navier-Stokes equations[C]//Proceedings of ETC 6th European Conference on Turbomachinery. Florence: European Turbomachinery Society, 2005: 7-11.

[16] SIMON J F, LÉONARD O. Modeling of 3-D losses and deviations in a throughflow analysis tool[J]. Journal of Thermal Science, 2007, 16(3): 208-214.

[17] SIMON J F, THOMAS J P, LÉONARD O. On the role of the deterministic and circumferential stresses in throughflow calculations[J]. Journal of Turbomachinery, 2009, 131(3): 031019.

[18] PERSICO G, REBAY S. A penalty formulation for the throughflow modeling of turbomachinery[J]. Computers & Fluids, 2012, 60(10): 86-98.

[19] TADDEI S R, LAROCCA F. Axisymmetric design of axial turbomachines: An inverse method introducing profile losses[J]. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 2008, 222(6): 613-621.

[20] TADDEI S R, LAROCCA F. CFD-based analysis of multistage throughflow surfaces with incidence[J]. Mechanics Research Communications, 2013, 47(47): 6-10.

[21] TADDEI S R, LAROCCA F. An actuator disk model of incidence and deviation for RANS-based throughflow analysis[J]. Journal of Turbomachinery, 2014, 136(2): 021001.

[22] PACCIANI R, RUBECHINI F, MARCONCINI M, et al. A CFD-based throughflow method with an explicit body force model and an adaptive formulation for the S2 streamsurface[J]. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 2016,230(1): 16-28.

[23] NIGMATULLIN R Z, IVANOV M J. The mathematical models of flow passage for gas turbine engines and their components: LS 198[R]. Neuillysurseine: AGARD, 1994.

[24] 袁宁, 张振家, 顾中华, 等. 涡喷发动机压气机三种S2流面计算程序的比较[J]. 推进技术, 1998, 19(1): 51-57.

YUAN N, ZHANG Z J, GU Z H, et al. A comparison of three kinds of calculation program of S2 stream surface in the compressor of aero-engine[J]. Journal of Propulsion Technology, 1998, 19(1): 51-57 (in Chinese).

[25] 季路成, 孟庆国, 周盛. 叶轮机通流计算的时间推进方法[J]. 航空动力学报, 1999, 14(1): 23-26.

JI L C, MENG Q G, ZHOU S. Time-marching method for through-flow computation of turbomachinery[J]. Journal of Aerospace Power, 1999, 14(1): 23-26 (in Chinese).

[26] 施发树, 刘兴洲. 多部件模型在全尺寸小型双函道涡扇发动机气流数值模拟中的应用[J]. 推进技术, 1998, 19(4): 22-26.

SHI F S, LIU X Z. Multicomponent models in application to numerical simulation of a small full-sized by-pass turbofan engine[J]. Journal of Propulsion Technology, 1998, 19(4): 22-26 (in Chinese).

[27] 施发树. 一体化弹用小涡扇发动机系统的气动热力数值模拟[D]. 南京: 南京航空航天大学, 1999.

SHI F S. Aero-thermodynamic numerical simulation of integrated small turbofan engine system[D]. Nanjing: Nanjing University of Aeronautics and Aerospace, 1999 (in Chinese).

[28] 于龙江, 陈美宁, 朴英. 航空发动机整机准三维流场仿真[J]. 航空动力学报, 2008, 23(6): 1008-1013.

YU L J, CHEN M N, PIAO Y. Quasi-3D simulation of aero engine full flow field[J]. Journal of Aerospace Power, 2008, 23(6): 1008-1013 (in Chinese).

[29] 曹志鹏, 刘大响, 桂幸民, 等. 某小型涡喷发动机二维数值仿真[J]. 航空动力学报, 2009, 24(2): 439-444.

CAO Z P, LIU D X, GUI X M, et al. Two dimensional numerical simulation of small turbojet engine[J]. Journal of Aerospace Power, 2009, 24(2): 439-444 (in Chinese).

[30] 金东海, 桂幸民. 某涡扇发动机考虑级间引气的二维数值模拟[J]. 航空动力学报, 2011, 26(6): 1346-1351.

JIN D H, GUI X M. Two dimensional numerical simulation of a turbofan engine with air bleeding in compressor[J]. Journal of Aerospace Power, 2011, 26(6): 1346-1351 (in Chinese).

[31] 李德英， 宋彦萍， 陈浮, 等. 任意曲线坐标系Euler方程S2流面的计算方法[J]. 西安交通大学学报, 2015, 49(7): 42-48.

LI D Y, SONG Y P, CHEN F, et al. Euler S2 stream surface calculation for arbitrary curvilinear coordinate system[J]. Journal of Xi’an Jiaotong University, 2015, 49(7): 42-48 (in Chinese).

[32] WAN K, JIN H, JIN D, et al. Influence of non-axisymmetric terms on circumferentially averaged method in fan/compressor[J]. Journal of Thermal Science, 2013, 22(1): 13-22.

[33] 万科, 朱芳, 金东海, 等. 周向平均方法在某风扇/增压级分析中的应用[J]. 航空学报, 2014, 35(1): 132-140.

WAN K, ZHU F, JIN D H, et al. Application of circumferentially averaged method in fan/booster[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(1): 132-140 (in Chinese).

[34] SIMON J F. Contribution to throughflow modelling for axial flow turbomachines[D]. Liege: University of Liege, 2007.

[35] MARBLE F. Three-dimensional flow in turbomachines[J]. High Speed Aerodynamics and Jet Propulsion, 1964, 10: 83-166.

[36] DENTON J D. Throughflow calculations for transonic axial flow turbines[J]. Journal of Engineering for Power, 1978, 100(2): 212-218.

[37] DANG T Q, WANG T. Design of multi-stage turbomachinery blading by the circulation method: actuator duct limit: 92-GT-286[R]. New York: ASME, 1992.

[38] BOSMAN C, MARSH H. An improved method for calculating the flow in turbo-machines, including a consistent loss model[J]. Journal of Mechanical Engineering Science, 1974, 16(1): 25-31.

[39] ADAMCZYK J J. Model equation for simulating flows in multistage turbomachinery[J]. Lecture Series—van Kareman Institute for Fluid Dynamics, 1996, 5: N1-N28.

[40] THOMAS J P, LÉONARD O. Investigating circumferential non-uniformities in throughflow calculations using a harmonic reconstruction: GT2008-50328[R]. New York: ASME, 2008.

[41] THOMAS J P, LÉONARD O. Towards a high order throughflow: Part Ⅰ—Investigating the effectiveness of a harmonic reconstruction for 3D flows: GT2010-22841[R]. New York: ASME, 2010.

[42] THOMAS J P, LÉONARD O. Toward a high order throughflow—Investigation of the nonlinear harmonic method coupled with an immersed boundary method for the modeling of the circumferential stresses[J]. Journal of Turbomachinery, 2012, 134(1): 011017.

(责任编辑： 鲍亚平， 张晗)

*Corresponding author. E-mail: jinguang_yang@dlut.edu.cn

Time marching based throughflow method: Current status and future development

YANG Jinguang1,*, WANG Chunxue2, WANG Dalei2, SHAO Fuyong2, YANG Chen3, WU Hu3

1.SchoolofEnergyandPower,DalianUniversityofTechnology,Dalian116024,China2.BeijingPowerMachineryResearchInstitute,ChinaAerospaceScienceandIndustryCorporation,Beijing100074,China3.SchoolofPowerandEnergy,NorthwesternPolytechnicalUniversity,Xi’an710072,China

Aiming to have an exploration on application potential of the throughflow computation method for turbomachines based on the time marching, the research progress home and abroad is concluded, and several key issues of application of this method are generalized, including geometric description of blade blockage, simulation of blade force, discontinuity problems in blade leading and trailing edges, and shock capture. Different computation models are also compared. It is summarized that the throughflow computation method based on the time marching has many advantages compared to the traditional throughflow method, and thus has the potential to be a competitive tool in advanced gas turbine engine components design and analysis, as well as in total engine steady and transient state simulation. This method is expected to have much more improvements in the aspects mentioned above, although it has achieved great progress. It is believed that this method will be adopted as a new standard tool in turbomachinery design and analysis with further developments and improvements.

turbomachine； throughflow； time marching； blade blockage； blade force； discontinuity at blade leading edge and trailing edge； shock

2016-11-28； Revised： 2016-12-19； Accepted： 2017-01-11； Published online： 2017-03-09 09:05

URL： www.cnki.net/kcms/detail/11.1929.V.20170309.0905.004.html

the Fundamental Research Funds for Central Universities (DUT15RC(3)035)

V231.3

A

1000-6893(2017)09-520996-13

2016-11-28； 退修日期： 2016-12-19； 录用日期： 2017-01-11； 网络出版时间： 2017-03-09 09:05

www.cnki.net/kcms/detail/11.1929.V.20170309.0905.004.html

中央高校基本科研业务费专项资金 (DUT15RC(3)035)

*通讯作者.E-mail: jinguang_yang@dlut.edu.cn

杨金广， 王春雪， 王大磊， 等． 基于时间推进的通流计算方法： 现状及展望 [J]． 航空学报, 2017, 38(9): 520996. YANG J G, WANG C X, WANG D L, et al. Time marching based throughflow method: Current status and future development[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(9): 520996.

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2017.620996