跨声速压气机流动失稳的建模与仿真

跨声速压气机流动失稳的建模与仿真

流失稳定是叶片机械面临的主要挑战之一。轴速损失是轴流压气机和压缩系统中普遍存在的流量不稳定型。当压气机偏离设计工况时,就有可能进入这种不稳定状态。压气机流动失稳边界的预测一直是设计人员十分关注的问题。自从旋转失速现象被发现以来,为了认识这种流动现象的物理机理,学术界产生了大量关于压气机不稳定现象的理论和实验研究工作。20世纪80年代以前,对于失速的研究主要集中在对现象的观察和物理描述。随后,一些学者在理论上提出:失速的发生是从线性扰动发展到非线性失稳的一个过程。McDougall等在一台低速压气机实验中观察到一些沿周向旋转的正弦波动,即所谓的模态波。之后Day在实验中首次检测到一种所谓的短尺度“Spike”类型失速扰动。20世纪90年代,研究人员开展了大量的实验工作(比如Camp与Day;Day等;Wilson与Freeman),研究失速演变过程中的各种扰动类型,Tan等对此作了综述性总结。除了以上两类失速先兆波外,Day等在4台高速压气机上开展实验研究,发现了两类新的物理现象:高频失速与固定位置失速。研究发现:在同一台压气机上还可以发现多种类型失速先兆波,或者失速过程中不同类型扰动波之间可以相互触发。进气畸变、雷诺数和级间不匹配也会影响失速类型。总之,不同压气机中存在多种样式的失速过程。即使是某一特定压气机,也会因工作状态不同而呈现复杂多样的失速演化方式。由于先兆波类型的多样性,目前尚未有任何统一理论能够解释旋转失速瞬态的细节演变过程。然而,各种类型的失速先兆波都是由小扰动发展为最终失速团的。尽管失速演化过程的研究具有重要理论意义,但工程实践中设计和使用人员更为关注的是风扇/压气机流动失稳起始点的位置。迄今为止,有关压气机旋转失速起始判定的研究工作主要有两类:解析模型和数值计算。几十年来,从不可压缩到可压缩,从二维到三维,解析模型不断发展完善。各种解析模型虽然在描述流动画面、提炼物理本质上各有侧重和不同,但其根本的出发点都是要建立判断压气机失稳的定性准则,分析影响压气机流动失稳的主要因素。这些模型在形式上都具有简洁明晰的特点,并在早期研究中产生了积极的影响。但是从另一方面讲,这类模型对压气机,特别是叶栅通道内的流动细节缺乏深入的了解和详细的刻画,叶片几何造型对稳定性的影响不能精确地包含在模型中,而这些因素常常对稳定性起到至关重要的影响。这类模型采用激盘或半激盘的形式代替复杂的真实叶片结构,同时需要引入压气机特性线、总压损失、叶栅出口气流角等经验关系加以补充。但是,通常情况下这些经验关系都不能精确地得到,特别是在压气机的实际设计阶段。模型预测的精度严重依赖所采用的经验关系与实际情况的符合程度。如果新的造型设计超出原有数据库中的经验数据范围,此类模型则往往不足以应对这种挑战,不能为设计人员预测失速边界提供有效工具。另一方面,数值模型[13,14,15,16,17,18,19]通常采用非定常计算流体力学(CFD)来判定失速边界。这类方法可以包含更多的真实物理信息,提供更丰富的关于失速过程的物理图画,增进研究人员对于失速过程的理解。但从工程实际的角度来看,对于失稳点的预测是最关注的问题,而常见的数值模型对于失稳点的判断却不如解析模型那样明确。更重要的是,非定常计算对CPU、内存资源的大量需求和较长的时间周期(动辄几万至几十万CPU小时)严重妨碍了其在当前工程设计中的广泛应用。解析模型和数值模型各自的特点引领我们思考:是否存在一种新的建立稳定性模型的方式,一方面像解析模型那样,能够为判断失速边界或者失稳点提供清晰的判断依据或准则;同时像数值模型那样,能够包含复杂几何因素对流场结构的影响,并且不像解析模型那样严重依赖经验关系。可以设想:这样一种在数值解法上可以快速实现、物理层面上可以包含更多因素、失稳点判断清晰明确的模型方法,可以作为稳定性设计工具而成为压气机设计阶段中的关键一环。目前,国际上绝大部分稳定性模型方法都是基于低速主流、忽略压缩性的假设。这个假设利于求解不可压的扰动解析解,便于构造压气机稳定性解析模型。对于低速压气机,相对马赫数小于0.3,流场的压缩性不强,可简化为不可压处理。然而,对于现代高速压气机,由于负荷高、叶片切线速度大,可压缩性就会明显。Tryfonidis等在8台高速压气机上应用相同的数据采集程序开展动态实验,发现了两类失速演化行为。其中几台压气机在失速开始之前出现近似50%转频的先兆波,这些扰动波长等于压气机圆周周长。这类扰动对应于Moore-Greitzer(M-G)模型所描述的不可压缩模态。而当这些压气机工作在较高转速时,实验中却难以在频谱分析中检测出这类扰动波。而对于一些工作在设计转速下的压气机,随着失速点的接近,实验中发现了以接近转频速度进行周向传播的扰动波。Weigl等证实了这些扰动波为可压缩模态。Bonnaure开展了多级压气机的二维可压缩性分析,结果表明扰动波的压缩性对高速压气机的稳定性行为有较大影响,然而该模型依然应用近似计算方法得到主流数据,并且预测精度依赖于损失和落后角的动力学经验公式。总之,压缩性对于失稳预测的影响仍不清晰,发展有效的跨声速压气机失速预测模型依然至关重要。针对压气机失速等流动稳定性问题,Sun等提出了一个关于叶轮机械流动失稳开始的通用理论,本文将应用该理论开展深入的研究工作。为了保证内容的完整性,本文首先简要介绍该理论,之后将其应用于压气机子午流面的可压缩和不可压缩流场,通过对比预测结果,揭示流动压缩性对于高速压气机失稳预测的重要影响。最后,通过一个单级跨声速压气机的稳定边界预测算例,验证本文发展的模型适用于跨声速多排叶片的失速起始点预测。1u3000navier-stkies方程的建立叶轮机械内部流动的不稳定现象属于流动不稳定性问题中的一种。同时,作为旋转机械,其变工况特征、复杂的几何形状和固壁边界使这一问题更加棘手。旋转失速和喘振是最主要的两类风扇/压气机内部流动失稳问题。从物理上讲,这两类流动问题同样属于流动在受到一定扰动后不再能够维持原有的流动状态而发生的失稳现象。如果只关注失速起始阶段的系统稳定特性,完全可以采用线性小扰动理论来预测失速的产生。另外,从工程实际的角度来看,对于失速起始点的判断是设计人员更加关注且至关重要的问题,并且能够满足实际工程的需要。本节将发展一个处理三维可压缩流动稳定性问题的通用特征值理论,其主控方程是NavierStokes方程组。式中:p为流场的压力;t为时间;ρ为流场的密度;γ为比热比;q为流场内热源热流量;V为速度向量;F为体积力;τ为斯托克斯假设下的黏性应力张量。假设黏性系数μ为常数,黏性项可以简化为由于本文关注流动失稳的起始阶段,因此假定流场由平均流和小扰动组成。这里上标“-”代表平均流场,上角标“′”代表扰动量,平均流动可以通过定常CFD计算求解,通常是计算三维定常雷诺时均Navier-Stokes(RANS)方程。将式(5)~式(8)代入方程组式(1)~式(4),得到线化形式的Navier-Stokes方程组。假设体积力为流场状态量的函数:式中:下标r、θ、z代表径向、周向和轴向方向;F为体积力的模;f为体积力扰动的模。对于特定的复杂流场,体积力的精确表达式仍较难以给出。理论上,可以建立体积力的简化模型来模化力源对流场的主要作用。需要说明的是,在建立通用流动稳定性问题的特征值方程的过程中,由复杂结构和固壁边界产生的体积力形式极其复杂,势必带来非常大的难度和挑战。另外,在特征值问题的建立过程中,需要将线化Navier-Stokes方程组从物理域转化为计算域。三维CFD计算中,通常贴体网格技术被用于处理复杂的固体几何。复杂的固壁边界将会增加方程组在两个区域转化时的难度。笔者认为:采用浸入式边界(IB)方法可以将空间扭曲的贴体网格方式替换掉,从而提供一种简单的方式来模拟复杂固壁边界。Peskin应用浸入式边界方法成功模拟了心脏瓣膜的动边界-血液相互作用的流场。这种方法是在动量方程的右端项处采用力源项来代表固壁几何对于流场的作用。数学上可以证明:如式(11)~式(13)所示,作用于流场上的物质边界等价于一个分布式力源。式中:u为速度矢量;x、r、s为位置坐标;χ、η′为体积力系数;S为表面积;εΩi(t)为体积元;δ为狄拉克函数。这样就不必生成复杂的贴体网格,用普通的正交网格就可以进行方程离散。在应用浸入式边界方法时,浸入式边界产生的作用可以用流体上的奇异力来反映,边界和流体的相互作用力采用狄拉克函数来刻画,而力的大小与流场速度等参数联系起来,将力源项写成各个流场参数相关的函数形式。当然,要把这样的力写成可以线化的形式仍然需要相当多的数学技巧,但理论上这种方式是存在的,并且可以根据具体实际问题加以简化。进一步,假定流场小扰动为谐波分解形式:式中:上标“~”代表流动扰动的振幅;m为周向模态数;ω为扰动复数频率;θ为周向弧度坐标。将式(14)~式(18)式代入线化Navier-Stokes方程组中,对应特定的流动采用适当的边界条件,构成封闭方程组,再经过坐标变换,借助于浸入式边界方法有可能将方程组从物理域转到相对简单的正交网格计算域上,最后得到一个二阶偏微分方程组如下:式中:A、B、C、D、Q、G、H为系数矩阵。则方程式(19)简写为式中:X(ω)为系数矩阵。因此,为了得到非零解Φ,可以表述为如下形式的特征值方程:在计算域上采用合适的数值离散方法,进而可以求解方程式(22),得到系统的特征复数频率为考虑到扰动波的表达式(14)~式(18),如果ωi小于零,则流动稳定;如果ωi大于零,则流动失稳。至此,本文已经详细论述了一个完整的流动稳定性通用理论。理论上,在整个模型建立的过程中,除了为研究流动失稳开始而做的线性假设外,没有做任何其他的假设,整个推导过程是严格成立的。笔者认为:这样的流动失稳理论有可能为以翼型失速、边界层转捩为代表的外流失稳问题,以及离心或轴流压气机失速为代表的内流失稳问题提供一种新的解决方案。接下来,本文将具体论述如何在这一通用理论基础上开展适用于叶轮机内部流动稳定性的研究。受湍流转捩问题研究中的稳定性理论启发,本节从最基本的Navier-Stokes方程组出发,应用线性化小扰动理论,结合浸入式边界方法和体积力方程,发展了一个通用流动稳定性模型。理论上,这样的模型可以用于研究任意边界情况下复杂流动失稳开始的问题。上文介绍了关于预测流动失稳开始的通用特征值模型理论的完整思路。然而,针对具体的带有复杂固体边界的三维可压缩复杂流动问题,如何构造体积力的表达式仍是难点;同时,最终得到的特征矩阵也会因为流动本身复杂性和数值求解精度要求等因素而变得非常庞大。一般来讲,主流的复杂程度直接决定了预测流动稳定性方法的难易水平。考虑到本文重点关注于叶轮机内部,特别是风扇/压气机流动稳定性问题,可以根据实际问题需要,将流场进行简化处理,得到通用稳定性理论的若干种简化形式。在马赫数较低的情况下,压缩性对于流动本身影响不大,可以将流动主控方程简化为三维不可压缩Navier-Stokes方程组;对于近似轴对称流动情况,可以重点关注子午平均流面;当径向方向流动可以忽略时,研究二维轴对称流动便能抓住问题的本质和关键。这些简化方式要根据具体实际问题的需要而适当选择。为了捕获影响流场失稳的主要特征,并且尽可能地与现有设计体系相适应,下文将首先研究包含体积力模型的子午平均流场的稳定性问题。实际上,下文所发展的体积力模型可以作为浸入式边界方法的一种简化后的特殊形式。2压气机失速预测简化模型2.1叶片流场动力学方程本文的控制方程是柱坐标系下包含体积力的三维Euler方程组,在绝对坐标下的形式如下:1)连续方程2)径向动量方程3)周向动量方程4)轴向动量方程5)能量方程式中:v和T为流场的速度及温度;下标r、θ、z代表径向、周向和轴向坐标方向;R和cv分别为热力学常数和定容比热。在叶片区域,下文中将把体积力进一步模化为流场状态量的响应函数。而在无叶片区域,体积力为零。2.2周向平均假设基础上的模型推导主控方程式(24)~式(28)中,流场状态量参数的平均量及其导数均来自定常CFD计算及数据后处理,通过求解RANS方程得到。根据模型简化程度的不同,也可以通过与模型相适应的其他方式得到,比如,可以通过求解通流计算程序得到一个子午面上的平均流场,这时模型方程中平均量的周向导数为零。下文主要讨论在周向平均假设基础上的模型推导。主流的计算采用发展较为成熟的三维CFD方法求解子午流面。理论上,这样可以使模型包含整个展向分布上的三维复杂叶片造型,特别是叶片前掠或后掠对主流流场的影响。值得一提的是,在Gordon与Gong的模型中,主流流场是通过求解不可压缩Euler方程组加体积力的方式来得到的。一方面,不可压缩前提简化了计算,但只能适用于低速风扇/压气机系统;另一方面,为了求得一条特性线上的主流流场,需要预先通过几个已知流场数据的工作点来得到体积力与主流流场状态量的参数表达关系式,也就是说,未知工作点的体积力需要通过体积力拟合曲线的内插或外推来得到。而据Denton与Dawes的工作表明:损失或落后角的拟合关系式适用范围非常有限,应该特别小心的处理才能得到较为合理的结果。对于旋转失速稳定性问题,主流的描述精度将严重制约最终稳定性的判定结果。因此有理由相信:通过较为成熟的CFD求解RANS方程来获得主流流场是更为可取且可行的办法。2.3叶片对流场的非定常扰动模型鉴于叶轮机流动稳定性的实际需要,为简化问题的计算复杂性,采用线化Euler方程组为基本方程,而黏性、叶片的作用及损失来源都通过体积力的方式来模化。这种模化方法曾被用于Longley、Chima、Gong的研究工作中,并得到了不断地完善。数学上,体积力模型并不是唯一的,只要恰当反映流场内的质量、动量和能量守恒关系,不同的建模方式均能在一定精度范围内刻画流动细节。相对于传统模型中所采用的激盘或半激盘方式,体积力方法可以包含真实复杂三维叶片造型的影响,而且相对于求解初边值问题分析稳定性的数值方法,计算量明显降低。本文将叶片对气流的影响分解为两部分:一部分垂直于叶片表面方向,为无黏体积力Fv,主要作用是引起气流偏转及做功;另一部分平行于叶片表面方向,主要是由叶片表面的黏性产生的,记为Ft。Xu的研究结果表明:相对于使气体偏转的部分,引起损失的黏性力部分相对较小。因此,本文忽略这一部分体积力,而将Fv直接标记为F。F的周向分量Fθ被分解为两部分:Fθ1使得气流转折,Fθ2产生流动损失。需要指出的是,之后的模型算例计算结果的精度将会验证这个简化的体积力模型以及接下来线化过程的合理性。采用三维定常CFD计算或者设计阶段的流线曲率方式获得风扇/压气机的子午流场数据后,根据周向动量方程,Fθ1可以由式(29)计算得到。与Gordon的分析类似,当地的叶片周向力会驱使气流相对速度转向叶片中弧线方向,即式中:β为叶片周向金属角;Ω为转子角频率;转折体积力系数kc(r,z)可由式(29)和式(30)计算得到。对于叶轮机内部的可压缩流场,定义绝对总压为转总压为由于叶片区域内的流动损失,转总压将沿流线减小。假设叶排内的流动损失正比于气流相对速度大小,两个点a和b位于叶片通道内的一条流线上,定义损失系数λ为另一方面,假定Fθ2正比于单位流线长度上的转总压变化,即式中:定义l为流线上坐标,当a与b两点很近时,因此,Fθ2正比于损失系数λ考虑到体积力与叶片中弧线的正交性,其轴向分量为假设叶片径向角γ′,体积力径向分量为本文考虑流道面积轴向变化不剧烈以及叶片的弯掠不明显的情况,叶片对流场的径向力较小,也可以忽略不计。接下来,把叶片对流场的非定常作用描述为对叶片区域当地非定常三维流动变量的响应函数,即体积力的非定常扰动被表述为带有延迟时间常数τ的一阶延迟方程。其中延迟时间常数τ表征了风扇/压气机的非定常特征,其量级一般为叶片通道内主流的对流时间,即叶栅通道长度除以叶栅内平均速度。本文经过算例验证表明:该模型对于这一常数并不敏感,因此下文算例中均采用了一倍通流时间作为延迟时间常数值。2.4扰动的频率和扰动频率在失速的早期发展过程中,先兆波的发展存在由小扰动到大尺度扰动的过程,因此在研究失速开始阶段,线性小扰动理论是符合流动物理特征的。另外,从工程实际的角度看,真实风扇/压气机不允许工作在旋转失速的失稳工作状况下,准确判断失速起始发生的位置已经足够。因此,本文应用小扰动理论线化主控方程组。在线性小扰动假设基础上,流场参数分解为平均量和扰动量之和的形式,如式(5)~式(8),假定速度、压力和密度的扰动量可以写成谐波形式(式(14)~式(18))。同样地,如果虚部为负ωi=Im(ω)<0,那么扰动是衰减的;如果虚部为正ωi=Im(ω)>0,那么扰动是增长的,此时系统失稳,实部ωr=Re(ω)为失稳的频率。实际上,对于叶轮机械内的失速开始问题,He证明:第一阶周向模态先兆波具有最不稳定的性质。Gordon的结果也表明:高阶模态通常都比低阶模态要趋于稳定,只有出现极其严重的级不匹配或者进气畸变时,高阶模态才会有可能比低阶模态提前失稳,而本文不研究这些极端情况。因此,在后面的算例校核中,本文也仅关注周向和径向最低阶模态的扰动频率。将式(5)~式(8)与式(14)~式(18)代入方程组式(24)~式(28),去除零阶量及高阶小量,假设进口无扰动,出口压力扰动衰减为零,固壁为滑移边界条件,得到最终线性化的扰动量将满足偏微分方程组:式中:偏导数。3资源值方程的解3.1特征方程的确定为提高离散精度,减少离散节点,进而减少计算量,本文采用基于Gauss-Lotatto节点的谱配置法(拟谱方法)对特征方程式(43)进行数值离散求解特征值。这种离散的格式广泛适用于各种边值问题,利用定义在区间-1≤ξi≤1,-1≤ηj≤1上的Chebyshev多项式TNr、TNz,Nr+1、Nz+1分别为径向和轴向总离散点数。其中一般地,在采用谱离散的时候,需要将物理域映射到适用于谱方法的正交计算域,例如:(r,z)→(ξ,η)。计算域中的一阶导数可以表示为式中:Erik和Ezkj为导数矩阵的元素。为了便于合理分配网格空间分布,在流场状态量梯度大的地方实现局部加密,本文采用分区谱方法,即将整个流场划分为进口区域、叶片区域和出口区域,在计算域分区界面应用扰动及其导数连续匹配条件。经过上述离散,并且结合适当的边界条件方程,在5×(Nr+1)×N(z+1)个节点上联立相应数目的封闭方程组,最终可以整理为这样的特征方程组形式:因此,特征频率ω由如下的特征方程来确定:为了体现失速先兆波的物理特征,将求得的特征频率实部和虚部分别进行无量纲化,即定义相对速度和衰减因子,分别代表扰动相对压气机的周向传播速度(RS)和扰动随时间的衰减程度(DF)。式中:f为转子转动频率;rT为转子前缘叶尖处半径;U0为转子进口远前方轴向来流速度。3.2基于实际应用由于离散点数依然较多,最终的矩阵珚A-iω珚B维数非常大,加上数值误差累计的影响,因此传统的求根方法(比如牛顿迭代法)都遇到了困难,使得直接求解方程式(57)几乎不可实现。于巍巍等曾经采用环绕积分方法求解稳定性特征值问题,取得了很好的结果。环绕积分法能够有效地处理复杂复数方程根问题,但由于这种方法需要求解积分区域非常多的行列式值,所以适用于矩阵维数比较小的问题或者计算函数值较为方便的问题。而对于本文模型,矩阵维数非常大,因此所需要的计算时间将非常长,且算法上难以通过并行计算来提高计算速度。Gordon曾经采用逆迭代方法求解特征值问题,Malik也曾用这种方法求解流动稳定性问题。相比于其他方法,这种方法虽然计算时间很短,但需要提供迭代初值,很多情况下这种方法的结果收敛性对初值非常敏感,容易收敛到错误的根,或者不收敛。另外,对于旋转失速问题,本身初始扰动的特征频率大致范围难以给定,初值也就难以给出。因此逆迭代方法也不适用于本文发展的风扇/压气机稳定性问题。理论上式(56)可以通过QZ分解方法一次性求得所有的特征根,然而由于矩阵珚A、珚B非常庞大,传统的QZ求特征值算法计算时间很长,并且精度不够,而且往往得到大量的伪根。鉴于上述方法的局限性,本文采用矩阵理论中的奇异值分解(SingularValueDecomposition,SVD)方法求特征频率。SVD可以将任意一个矩阵分解成为两个正交矩阵U和V,以及一个实对角矩阵D。其中D矩阵的对角元素即是矩阵A的奇异值。矩阵A的条件数被定义为它的最大与最小奇异值之比。严格来讲,一个矩阵是奇异的,即这个矩阵的条件数为无穷大。但是,考虑到实际计算机的存储能力,只要矩阵条件数大于计算机精度的倒数,就可以认为此时矩阵奇异。然而,由于实际数值计算时的残差及迭代误差积累,一个矩阵奇异很难真正达到非常大的条件数。因此,考虑实际情况,只要矩阵的条件数远大于矩阵元素的某种误差形式,即可以认为矩阵奇异,对于本文所应用的维数非常大的矩阵特征值问题来讲,对应的ω就是待求的特征频率。综上,本文在复平面离散网格点上逐点求解矩阵珚A-iω珚B的奇异值,实际计算结果中,选取复平面上条件数倒数与周围位置点相比显著减小处的极小点作为特征值。这种数值方法也在类似特征值求解工作中得到应用。同时,数值运算中也将检测对应条件数非常大的ω,矩阵的最大奇异值并非很大,并且矩阵珚A-iω珚B行列式的确很小。另外,由于所求的ω具有明确的物理意义,代表着扰动波的频率。因此通过变换不同网格节点数,真正有物理意义的ω值应该不会变化位置。依靠这一点可以进一步去除伪模态的影响,最终得到具有真正物理意义的特征频率。4rans方程的主要特征如果将特征值方程中的密度扰动强制为零,便可将本文发展的三维可压缩稳定性模型转化为不可压模型。本节将应用可压缩和不可压模型,借助于NASARotor37(R37)高低两个转速算例,对比研究扰动的可压缩性对稳定边界预测精度的影响。需要指出的是,主流依然是求解三维可压缩RANS方程而得到。R37是20世纪70年代由美国国家航空航天局(NASA)Glenn研究中心设计的4个高压压气机进口级之一。Moore与Reid对其进行了详细介绍。Suder与Celestina用激光测速仪和探针对R37流场进行了详细测量,获得了不同截面上的总压和总温分布。正是由于该转子排除了转静干涉的影响,有详尽的实验数据,同时又包含了激波、边界层分离、叶尖间隙流动、二次流、径向掺混和尾迹等复杂流动,因此在跨声速转子研究中非常具有代表性。1993年的美国机械工程学会国际气轮机分会(ASMEIGTI)将R37作为盲题算例,用来评价当时叶轮机械CFD程序对压气机总性能和流场细节的预测能力。对于该转子的主流计算,本文不再详述,特性线计算结果见文献,将主流计算结果代入模型中即可得到系统的特征频率。4.1叶尖失稳转速变化图1为R37在60%转速下应用不可压模型计算得到的系统特征频率。在所选择的求解区域中,每个工作点上共有5个特征模态频率解。随着截流过程流量减小,衰减因子大体趋势上逐渐增大,模态1最先在10.89kg/s处跨过失稳临界值,与实验失稳点流量相对误差为3.1%。模型计算的失稳点对应的相对速度为32.5%转频,并且失稳前的变化都比较平缓。R37在60%转速下叶尖最大相对马赫数约为0.85,属于高亚声速流动。由表1对比可知:扰动的压缩性使该系统趋于稳定,但相比于可压缩性模型而言,流量误差仅增加1.3%,然而失稳扰动波的传播速度计算结果却降低了21.9%。因此,对于这个高亚声速流动算例,扰动的压缩性对于失稳点流量计算影响不大,而对扰动频率影响较大。4.2扰动的压缩性图2为R37在设计转速下应用不可压模型计算得到的R37系统特征频率。在所选择的求解区域中,每个工作点上仅有一个特征模态频率解。因此,笔者将可压缩模型和不可压模型计算得到的最不稳定(衰减因子最小)的特征模态解放到一起比较。随着截流过程,可压缩模态在19.37kg/s处出现失稳。然而,不可压模态的衰减因子在整个定常收敛解的流量范围内均为出现由负变正,即模型预测该算例流场一直处于稳定状态,而这显然是错误的。随着流量减小,不可压模态的相对速度先略微增大,后略微减小,整个变化趋势非常平缓,范围是40%~50%转频,明显小于可压缩模型计算结果。需要说明的是,对于该单转子算例,定常流场计算耗时约18CPU小时,而失稳边界预测总共耗时约990CPU小时。显然,相对于非定常数值模型求解失速瞬态过程时动辄几十万CPU小时的计算量,本文模型是具有优势的。表2对比了可压缩模型和不可压模型对于R37设计转速失稳点预测的计算结果,其中不可压模型计算失稳点为衰减因子曲线近似延长后的推断结果。R37在设计转速下叶尖最大相对马赫数达到1.4,属于跨声速流动。由该表对比可知:扰动的压缩性使该跨声速流动系统趋于稳定,相比于可压缩性模型而言,失稳点流量误差增加超过9.1%,失稳扰动波的传播速度计算也低了近40%。因此,对于这个跨声速流动算例,扰动的压缩性对于失稳点流量和传播速度计算影响都非常大,失稳点预测出现严重偏差。这也证实:对于跨声速流动稳定性预测,扰动压缩性的影响非常重要。因此,下文算例采用的是可压缩模型。5单级跨声速压气机动力学模型上文应用本文发展的稳定性模型计算了低速或高速单转子的稳定边界,本节讨论模型在跨声速单级压气机上的预测结果精度,以一台典型的高负荷、高压比压气机进口级NASAStage35作为研究对象,计算85%转速下的失速边界。单级跨声速压气机Stage35是NASA为先进核心压气机研制的一系列低展弦比进口级之一。Stage35在设计转速为17189r/min时的转子叶片切线速度为455m/s,设计压比为1.82,设计流量为20.2kg/s,其几何参数及实验结果见文献。5.1网格特征和转速本算例进行有间隙定常CFD计算,网格总数为144万,最小正交性为20°,最大长宽比为985,最大扩张比为3.1。进口区域采用H型网格,网格点数为I×J×K=33×41×33;转子区域采用O型网格,I×J×K=33×41×145;静子区域采用O型网格,I×J×K=33×41×145;出口区域采用H型网格,I×J×K=33×41×33;叶尖间隙为0.4mm,此区域采用蝶形网格,径向布置9个网格点。采用Spalart-Allmaras湍流模型,空间格式为中心差分,时间格式为四阶龙格库塔法,应用多重网格和局部时间步技术加速数值收敛。固定转速为14610.4r/min,即85%转速。出口边界为径向平衡方式,逐渐增大出口机匣背压,计算得到Stage35的特性曲线。图3对比了CFD和实验的85%转速总静压升特性曲线。由图可知:计算结果和实验特性曲线趋势一致,总静压升在设计点和近失速点非常接近,堵点流量减小了约1.3%。数值发散点位于14.38kg/s,而实验测得的失速开始点为15.608kg/s,两者相对误差为-7.87%,失速裕度被高估了超过8%。5.2算例2:网络流场+固定流场图4所示为Stage35算例求解系统特征频率时所采用的物理域网格分