风扇压气机进气畸变周向平均体积力模型三维非定常流

STYLEREF北航章目录北京航空航天大学硕士学位论文PAGEivPAGEiii摘要压气机进口的均匀流动条件是工程设计所期望的，但是当把发动机安装接到飞机上时，在飞行过程中，就不能保证在各种状态下压气机或风扇的进口流场都是均匀的了。压气机进口畸变带来的直接严重后果是压气机的稳定裕度下降和增压比、效率降低，因而加强进气畸变对压气机性能和稳定性影响的研究就显得非常重要。本文采用三维欧拉方程加体积力模型的方法，发展了一套三维非定常进气畸变快速数值模拟程序，藉此可以实现风扇/压气机在进气畸变条件下的全环面三维非定常数值模拟。该方法从物理基本原理出发，能够较好地反映工作于非均匀流场中的压缩系统以及叶片通道内部的非定常三维流动本质。程序用极小的计算机资源即可以实现对风扇全环面三维非定常畸变流场的模拟。利用该程序模拟计算了某航空发动机风扇进口稳态总压畸变条件下特性和流场结构，并与试验结果进行了对比，结果表明所能够以较高的精度预估气流畸变对压气机/风扇特性的影响。另外还计算了稳态总温畸变，动态总压振荡和温度脉冲对该风扇流动结构的影响，计算结果与目前畸变试验和理论研究结论一致，说明该程序可以很好地模拟稳态和动态总温、总压畸变对压气机流动结构和特性的影响。关键词：风扇/压气机；进气畸变；周向平均体积力模型；三维非定常流

AbstractFortheevaluationofstabilityofaero-engine,amajorlimitingfactorofcompressorinstabilityisinletdistortion.Longlengthscaleunsteadinesssuchasthatcausedbyloworderinletdistortioninallaeroengine,haveasignificantinfluenceinengineperformanceandstability.Thestallmargin,pressureratioandefficiencycanbeconsiderablyreducedinoperatingenvironmentsforwhichtheinletconditionsareken-uniform.So,itveryimportanttoemphasizetheresearchofeffectofcompressorperformanceandstabilityresponsetoinletdistortion.Thispaperpresentsacomputationalmodelfordescribingcomplexthree-dimensionalnon-lineardisturbancesassociatedwithinletdistortionoftheflowthroughthefullannulusincompressors.Thedifficultyproblemtosolveinletdistortionflowfieldistreatwiththeflowinductsandintra-blade-rowgapswhichisdescribedbythethree-dimensionalunsteadyEulerequationsformass,momentum,andenergyconservationwithsourcetermsofbodyforce.Smallstructuresofbladerowpassagescaleofflowpassingthroughthefilterwillbesimulatedviabodyforcemodels.Largescalesaboveacertainthresholdaredirectlycalculatedfromthegoverningequationsoftheflow.Adequatenumericallysimulationofthesephenomenabysuingthiscomputationalmodelonlyrequiressmallercomputerresourcesthancurrentthree-dimensionalCFDofNavier-Stokesequationwithappropriateturbulencemodels,theywillbeinthereachofthedesignerwishingtoperformsuchcalculationsonaroutinebasis.Inthispaper,ourcomputationalresultsofcompressorresponsetototalpressureinletdistortionshavebeenverifiedexperimentallyatawiderrangeofdifferentoperationspeedofafanrotor.Thepredictionsofthemodelformassflow,pressureratioandefficiencyofthisfaninoperatingenvironmentsofinlettotalpressuredistortionsareinagreementwithexperimentaldata.Inaddition,thismodelsomepredictionsofperformanceresponseandflowstructuretoinlettemperaturedistortions,inletpulsanttotalpressureandinlettemperaturepulse.Themodelpredictionssatisfycurrentempirical,numericaloranalyticalsources.Theserepresentaadvantageinthepredictionsoftheinletdistortions.Keywords:fan/compressor,inletdistortions,circumferentialaveragebodyforcemodel,three-dimensionalunsteadyflow目录第一章绪论…………….11.1选题意义…………………….11.2课题背景……………………21.3本文工作内容……………….8第二章压缩系统的建模………………102.1管道中的流动……………….102.2叶片排中的流动…………….10第三章三维非定常欧拉方程求解器……………….123.1网格生成…………………….123.2基本方程………………….123.3计算格式…….…………….133.4程序的流程图………….………………..14叶片通道周向平均体积力模型建立…….………………….184.1周向平均体积力模型原理………………..184.2叶片力的分解…………….………………194.3叶片对气流切向摩擦力的模化………194.4叶片对气流法向压力的模化……………..214.5堵塞判断……..…………..234.6落后角计算………..………………….….25第五章进气畸变条件下风扇定常与非定常流场数值模拟………….265.1算例风扇转子的主要气动参数…………..265.2计算网格…………………..265.3均匀进气条件下数值模拟………………..275.4稳态进气畸变条件下模拟计算……….…315.4.1稳态压力畸变………………..…..315.4.1.1给定进口畸变场…………………….….315.4.1.2计算结果与试验结果对比……………..325.4.1.3进口总压畸变流场流动结构分析……..335.4.1.4风扇转静子在进口总压畸变条件下流场数值模拟计算……………..355.4.1.5利用扩散因子进行稳定性预估的初步研究…………..365.4.2稳态总温畸变…………………...375.5动态进气畸变条件下模拟计算………………….…….405.5.1进口周期压力振荡………………405.5.2进口温度快速瞬变……………….425.6结论……………………….45研究总结与展望……………………….47参考文献……………….49攻读硕士期间发表的论文及所取得的研究成果………………..……51致谢…………………….52PAGE52PAGE51第一章绪论1.1选题意义压气机进口的均匀流动条件是工程设计所希望的，当压气机装在压气机试验台上进行试验时，试验台的进气装置保证了压气机进口流场的均匀性，当发动机在地面试车台上试车时，气流经过工艺进气道而流入发动机的压气机，这时压气机的进口流场也是相当均匀的。但是当把发动机安装到飞机上时，在飞行过程中就不能保证在各种状态下压气机或风扇的进口流场都是均匀的了。飞机爬升、大偏航角、侧风和进气道不起动工况等都会造成压气机进口处总压、总温、速度和流向不均，从而引起压气机性能的恶化。进气畸变对压气机工作的影响有两个方面，一是空气动力学方面的影响，另一方面是气体力学方面的影响。前者的影响使压气机的性能和稳定性变坏，后者的作用是叶片上非定常力引起机械振动。本文主要研究压气机对进气畸变的气动响应。在航空发动机性能评估中，一个影响其性能的重要的因子就是进口流场畸变，而进气畸变对发动机性能影响最直接、最重要的表现就是对压气机性能和稳定性的影响。压气机进口畸变带来的直接严重后果是压气机的增压比、效率降低和稳定裕度下降。因而加强进气畸变对压气机性能影响的研究就显得非常重要。现在高性能航空燃气涡轮发动机的研制和发展过程中，为了实现高性能的发动机设计指标，压气机设计者总是希望能够对畸变与发动机相互作用的物理本质进行深入的认识和理解，压气机设计必须同时进行性能和抗畸变能力的综合考虑，设计者必须掌握能够预测不同设计参数情况下压气机抗畸变能力的计算工具。研究进气畸变对压气机性能影响主要有两条途径，第一条途径是试验研究，在压气机运转过程中通过在压气机前方安装畸变发生器来模拟各种状态下压气机进口畸变场，然后对进口畸变场和压气机性能进行详细测量，获得各种畸变条件下压气机特性，但模拟整个使用条件范围内的全科目试验成本非常昂贵，而且试验研究本身迄今为止对加深物理问题认识的帮助甚少。在国际上基本上也是通过大量试验来获取压气机设计参数与抗畸变能力的工程关联，但是有关的试验数据无法在公开的文献中获得，另一方面有限的研究结果都属于敏感情报，有关各国都极力保密。近些年来计算流体力学的发展为畸变问题研究开辟了另一个途径——数值模拟，我们希望通过数值模拟来获得压气机在进气畸变条件下的特性，在设计阶段就对压气机的抗畸变能力有所了解，利用数值计算来指导设计，提高压气机设计的成功率。进口畸变与多级压气机相互作用问题是一个复杂的三维、非定常、非线性问题，数值模拟方法既要保证抓住畸变问题的物理本质，其计算量又要能够被工程应用所接受。在过去半个多世纪压气机发展的过程中，有许多研究者通过采用各种近似简化假设，发展可许多预测压气机畸变性能的近似模型，但因其计算体系广泛采用了一维、二维流动假设，无法反映真实流动的三维性，而压气机众多设计新理念是基于三维流动过程控制和设计，因此这些计算系统不能反映出这些新的设计思想的效果，无法评估发动机新的设计理念对发动机抗畸变能力的影响。接下来很容易想到的是用非定常的RANS来直接模拟，它几乎不依赖于经验，能够计算叶片槽道内、外很细微的非定常流动结构，但用非定常RANS数值模拟来研究畸变流场面临的一个主要困境就是非定常扰动尺度范围过大。压气机在畸变条件下各种不同尺度的非定常扰动都将会显示出来，一方面，进口气流畸变是以压气机直径为特征尺寸的大尺度非定常扰动，求解畸变流场就要求计算压气机整个环面所用通道中的非定常流动。另一方面，为了正确估算粘性作用，还需要求解比边界层尺度还小的扰动，对于高雷诺数小尺度扰动的长度可以达到几十个微米的量级，精细网格和全环面多通道计算使得计算网格数量惊人的巨大，它对于计算机资源的高要求就不止包含内存，也包含CPU计算时间。因此利用非定常RANS求解畸变流场需要极为苛刻的计算机硬件能力，其计算所费机时是目前工程设计工作无法接受的。本课题的目的是寻找一个即能真实地反映压气机内畸变流动的三维、非定常物理本质，又适合于工程应用的快速压气机畸变数值模拟方法，可以利用它来探索畸变对压气机性能影响的机理，研究不同设计参数对压气机抗畸变能力的影响，在压气机设计阶段就可以对其畸变性能作出预估，对压气机的性能和抗畸变能力做出综合考虑，降低设计风险，提高压气机设计成功率。1.2课题背景在进气畸变数值模拟研究方面，已经比较成熟的技术是利用一维非线性模型加平行压气机理论来研究有畸变进气的压缩系统，美国的Arnold工程发展中心和俄罗斯中央航空研究院都利用该方法发展了非常完善的计算软件。所谓一维模型，也称stage-by-stage模型。它将整个压气机按照级沿着轴线方向划分成为一个一个的计算单元（一维）。针对每一个计算单元，列出气流运动的控制方程，通过对所有单元控制方程的联立求解分析，获得压气机性能。而平行压气机理论，是将压气机沿周向分解成许多独立段，每一段构成一个子压气机，不同的子压气机之间相互并行工作，具有相同的出口边界条件。畸变进气的流动条件则表现为：在不同的子压气机进口具有不同的进气条件，所有并行的子压气机都按无畸变的压气机特性工作，他们之间通过出口边界条件联系起来[1]。平行压气机理论可以成功地解释一些物理现象，但也具有很大的局限性，该模型不能考虑到不同的进口流场区域在压气机通道内流动过程中的气流掺混，而且一维模型由于对非定常流动的模拟方式过于简单，无法反映压缩系统内部三维的流动本质，直接影响到预测结果的可靠性。图1平行压气机理论示意图早在1978年，Greitzer和Strand在研究施加在周向涡流上的畸变[2]时就发现畸变流场表现为三维特性，为了确定进气畸变对压气机的影响，需要深入认识进气畸变下压气机流场细节。但是由于计算技术和试验技术发展的限制，针对畸变进气条件下发动机三维流场的细节认识及其理论计算研究工作一直进展不大。近年来数值方法和计算机科学的发展使计算流体力学有可能对畸变问题进行模拟研究。似乎我们可以通过时间精确的三维RANS方程来求解压气机内的非定常畸变流场，但问题在于畸变问题是个扰动尺度相差几个数量级的非定常问题，大尺度的扰动可以达到米的量级，想要捕捉到这种扰动在压气机中的传播和衰减就要求能够解压气机的整个环形通道，而不是象我们现在通常所做的非定常计算时每个叶排只求解几个通道。另一方面所要求解的最小尺度比边界层的尺度还要小，这样才可以正确的估算粘性作用，因此计算网格尺度必须非常小，对于高雷诺数可以达到几十个微米的量级。精细的网格加上整台压气机全环面所有通道都要进行非定常计算，这样的计算量对计算机硬件能力和计算时间的需求是工程设计工作所无法接受的。例如，由Unno等人（2001）做的模拟风扇—吊架结构用了大约500万个网格在大型的并行（166个CPU）计算机上计算[3]。目前仅有美国NASA和英国帝国理工学院两家有此计算能力并在进行相关真实流动情况下全压气机/发动机流场的非定常流动计算，但即使使用高密度计算机群，其海量的计算使一个简单算例的计算也要进行几周甚至几个月。这样的工作周期对工程分析和设计而言是远水解不了近渴，而且计算结果是详细的非定常流场数据，如何从中提取有用的相关物理信息，尚需复杂的后处理分析工作[4]。针对低维模型无法反映流动的三维特性，而直接求解RANS方程又因计算机硬件和计算时间的限制而无法应用于工程研究，1985年Adamcyzk给出了他的通道平均多级流动模型[5]。这种方法直接求解压气机中的三维流动，但并不直接计算叶片对气流的作用力，而是通过体积力模型来模拟，因为不需要直接求解叶片力，计算量比RANS要小很多。通道平均多级流动模型主要思想是忽略叶片通道周向不均匀性，物理上可以理解为周向有无穷多个叶片，这样叶栅通道内的某些流动信息将丢失，但因为进气畸变这种非定常扰动的尺度远远大于叶栅栅距，我们关注的是叶片排对非定常三维非均匀流的响应，叶栅通道内部小尺度的非定常响应往往不是我们关心的重点，因此这种忽略小尺度扰动的通道平均模型应用于进气畸变模拟是可以接受的。叶片对气流的作用可以用激盘和半激盘模型或者通道平均体积力模型来实现。1997年Joo和Hynes基于圆盘激励模型[6]，1998年Gong等人基于通道平均体力模型[7]来模拟计算进气畸变，这两种方法都需要通过均匀的轴对称特性来推断在非均匀扰动作用下的叶片排性能。通道平均体力模型比激盘模型的优势为可以计算叶片排内部的流动，体积力源项可以沿径向和轴向有不同的分布，在数值解法上体积力模型摆脱了激励圆盘模型中难应付的人工边界。尽管如此体积力模型仍然存在一个问题，体力的源项在空间的分布基于均匀的轴对称特性，但在在非设计点比如接近失速时，叶片负荷将会经历大的弦向和径向的再分布，这就还需要一些特殊的处理[4]。1998年Gong等将通道平均体积力模型应用于畸变这种非定常、三维问题的研究中[7]，他的作法是不把叶片放入流场中参与计算，用无粘的三维非定常欧拉方程求解管道流动，叶片的作用以非定常体积力形式加入到方程源项中。这种方法的关键技术是如何准确地建立体积力源项与当地动力学量的关联，尤其是在非设计点。Gong分别采用两种方法来建立体积力源项，对于低速大展弦比压气机仍采用通过均匀的轴对称特性来获得体积力分布，即通过每排叶片各截面的压力升特性等参数换算出体积力源项，但是源项沿叶片弦向的分布则需要依赖经验来获得。而对于跨音级高速小展弦比的压气机他选择采用另一种方法，这种方法首先利用CFD计算或测试获得压气机的三维流场，对流场进行周向流动推力平均后得到轴对称流场，然后对轴对称流场应用动量方程求出体积力，经过一系列的不同状态下的CFD计算后，总结出体积力源项与当地动力学参数的函数关系，这种方法的优点是对经验的依赖较前一种方法要少。Gong利用这个模型模拟计算了NASA的单级跨音压气机Stage35的稳态总压、总温和动态总压畸变流场，成功地捕捉到了高速压气机独有的二种特性：阻塞和激波，验证了周向平均体积力模型能够对压气机畸变引起的流动再分配及畸变在压气机中的传递有很好的刻画。图2和图3为Gong计算总压和总温畸变条件下压气机内速度系数沿流程的变化，可以明显地看到畸变流场的三维流动特性。图2Gong利用体积力模型计算的稳态总压畸变条件下速度系数沿流程的变化图3Gong利用体积力模型计算的稳态总温畸变条件下速度系数沿流程的变化2001年Xu等人提出了粘性体力模型[4]，并利用该模型来模拟畸变流场。Xu的粘性体积力模型与Gong的周向平均体积力模型的不同是他把叶片放入流场中参与计算，无粘叶片力直接由方程计算出来的，而粘性影响才作为当地源项在粗网格上进行模拟的，被模化的只是叶片槽道中的粘性力作用，因此这种模型相对于粘性力和非粘性叶片力都需要模化的模型对经验的依赖程度要小，而且它有能力求解在叶栅通道内的详细流动，可以解由于叶片扫过而产生的扰动，并可以很精确的给出叶片的非定常力。但是因为方程求解时要计算叶片非粘性力作用，计算耗时要比体积力模型高一个量级。Xu发现一般情况下粘性力相对于总的叶片力是很小的，它的非定常扰动则更小，只需要一个很简单的阻力系数模型就可以很好地模拟时均和动态粘性作用，但对于有分离的情况，阻力系数的变化将会变的很大，这时就需要更详细的模拟。Xu在研究模拟粘性体积力的方法时发现如果直接计算无粘的叶片力，而模拟粘性力将会十分有利，因为无粘非定常力通常对扰动有很强的动力学响应，因此模拟无粘非定常力非常困难。还有，由于无粘叶片力是总的非定常叶片力的主要部分，模拟整体叶片力的任何误差的影响将会比单独模拟粘性影响重要得的多。换句话说，与周向平均的体积力模型相比粘性体力模型的模拟可以承受更大的误差[4]。图4为Xu利用粗网格粘性体积力模型计算的进气畸变压力场与精细网格RANS计算计算结果对比，两者计算静压和总压场十分相近。粘性体积力模型RANS图4Xu利用粗网格粘性体积力模型计算的进气畸变压力场与精细网格RANS计算结果对比比较Gong的周向平均体积力模型和Xu的粘性体力模型，周向平均体积力模型由于在管道内只去解通道平均的流动，没有把叶片几何形状考虑进去（全部组合进源项中），因此可以采用更规则的网格，在计算速度上比粘性体力模型高效一个量级，更适合于工程应用。然而，通过通道平均叶栅通道内的流动有关信息将会丢失，并且不可能去研究叶片对非定常扰动的响应。粘性体力模型是较周向平均体积力模型进行了改进的，由于其在叶栅通道内解非粘性流动，因此它有能力在叶栅通道内解非定常流动并可以很精确的给出叶片的非定常力。综上所述，利用数值方法来模拟非定常进气畸变流动的计算模型主要包括高度简化的低维模型如平行压气机理论模型和激盘模型还有周向平均体积力模型、粘性体积力模型及RANS直接数值模拟等。2001年Xu等人为了解决叶轮机械内部非定常流动数值模拟所面对的非定常尺度跨度大的问题，提出了一个对各种尺度的非定常流动都适用的通用框架[4]，框架示意图见图5。框架最上端是低维模型，它对计算机资源的需求最小，但对流动做的简化最多，计算中被忽略的流动结构需要依靠经验来建模，简化模型越多，对经验的依赖也就越大，经验总结的正确性决定了它计算的精度。框架最底端是直接数值模拟，它对计算机资源有极高的要求，所有的流动细节都通过控制守恒方程直接计算，几乎不依赖于经验。框架中间的方法则介于两者之间。这个框架更像是广义大涡模拟，在某一预设的临界值以上尺度的非定常流通过流动的控制守恒方程直接计算，低于临界值的小尺度扰动将会通过不同的模型——通常是由经验、数值、解析方法或者这些方法的结合来模拟。当我们要对某一个特殊应用选择一个合适的模型时，就要在经验和计算机资源之间做一个折衷。图5各种计算模型对计算机资源和经验的需求层次1.3本文工作内容本文工作的目的是建立一个适用于工程研究应用的畸变流场计算程序，要求能够真实地反映压气机内畸变流动的三维、非定常物理本质，并且计算速度快，在压气机设计阶段结可以对其畸变特性做出预估，为方案筛选提供依据。目前高性能压气机设计更加注重三维性，宽弦、掠型、端弯等设计先进的设计手段在设计中频繁应用，简单的低维模型模型不能模拟出这些三维因素对压气机抗畸变能力的影响。而对于畸变问题直接求解三维RANS方程、大涡模拟(LES)和直接数值模拟(DNS)则需要耗费大量计算机资源和时间，在目前的情况下还不能应用于工程设计阶段。鉴于这种情况，本文采用了周向平均的体积模型来模拟计算畸变流场，周向平均的体积模型居于非定常数值模拟框架的中间，它的优点是对经验和计算机资源的依赖都相对较少，可以采用更规则的网格，计算速度快，在模型建立合理的前提下，数值更准确，更适合于工程应用。而且，因为在整个压气机流场内求解的是三维非定常流体控制方程，因此能够真实地反映压气机中畸变流场的三维、非定常流动本质，为研究进气畸变和压气机相互作用的物理机理提供了可能。本文的主要工作是构造一个三维非定常欧拉方程加叶片通道周向平均体积力模型的进气畸变条件下压气机流场计算程序。其中主要包括三部分内容：第一部分是三维非定常欧拉方程求解器的调试。本文所使用的欧拉方程求解器具有三阶时间精度和两阶空间精度的欧拉方程求解器，可以计算双涵道多级压气机非定常流场，作者主要负责程序的阅读、修改和调试及验算工作。第二部分是建立周向平均体积力模型，这是本文工作的重点，体积力模型建立的好坏直接关系到程序计算正确性及精度。与Gong相同我们都采用了周向平均的体积力模型来模拟叶片排对气流的作用。但是体积力模型的构造方法与Gong不同，Gong的体积力模型需要预先输入各个叶片排的基元特性然后推算出体积力，或者根据一系列CFD计算结果归纳出体积力与当地动力学参数的函数关系，而本文的体积力的模化方法是将叶片对气体的作用力分解为粘性摩擦力和无粘压力两部分分别模拟，粘性摩擦力用简单的壁面摩擦阻力系数来模化，无粘压力则通过气流流过叶片后的周向动量变化，利用角动量守恒方程计算。同时因为计算时没有加入真实的叶片排，无法考虑叶片排对压气机通道的堵塞，在程序中加入了堵塞判断模型，当判断出叶片某一截面发生堵塞时，程序将大幅度提高摩擦阻力系数，使近堵点的计算更接近真实情况。本文的体积力构造方法主要基于叶型基本设计参数，所需参数与通流设计所要求的参数略同，因而特别适合于设计阶段的方案评估。第三部分为程序验算，选取某航空发动机单排风扇转子作为计算对象，计算了该转子均匀进气条件下相对换算转速=0.6～1.0转速下风扇特性，并与试验进行了对比；还分别计算了进口总压稳态畸变条件下风扇流场特性，也与试验进行了对比；计算了进口总温稳态畸变条件下风扇流场、进口动态温度脉冲条件下风扇流场、进口动态压力振荡条件下风扇流场，充分验证了程序计算压气机在各种畸变类型包括动态畸变流场的能力。

第二章压缩系统建模对于进口流场畸变这个三维、非定常问题本文选择用非定常三维欧拉方程加周向平均体积力模型来建立计算模型。该计算模型与Gong所采用的模型[7]相似：在整个压气机通道中利用三维无粘非定常欧拉方程来求解管道流动，不把叶片放入流场中参与计算，叶片作用以体积力的形式加入到方程源项中。2.1管道中的流动气流流经压气机通道时，受到叶片对它的法向压力和切向摩擦力、上下流道的摩擦力和气流本身的粘性力，在正常的流动情况下，粘性力主要发生在附面层中，主流区中粘性力相对很小，因此本文忽略主流区的粘性作用，同时压气机上下壁面对气流的摩擦力也忽略不计，在这种情况下压气机前后的无叶片段和叶片排间隙段都可以直接利用无粘的欧拉方程求解。在压气机进口管道、出口管道和叶片排间隙中采用三维非定常欧拉方程可以如此表示：（2.1）其中：ρ为密度；p为压力；e为内能；Vx、Vr、Vθ为轴向、径向、周向速度2.2叶片排中的流动在叶片排区域内叶片对气流的法向压力和切向摩擦力也不直接计算，而是用一定的方法进行模化，模化后的叶片力以体积力的形式被加入到欧拉方程的源项中，因此在叶片排区域求解的是带有体积力源项的欧拉方程，方程具体形式如下：（2.2）其中Fx、Fr、Fθ分别为体积力源项在轴向、径向、周向的分量，为体积力对气流做的功。

第三章三维非定常欧拉方程求解器3.1网格生成程序有自己的网格生成程序，在已知压气机流道尺寸的情况下可以自动生成3维H形网格，程序所采用的网格有三个特点：第一，因为畸变流场模拟要求计算压气机整个环形通道内的流场，因此其网格也必须沿整个压气机通道生成，而不是象通常三维计算时只生成一个或几个叶片通道区域的网格；第二，程序在计算叶片区作用时叶片排对气流的作用力不是直接计算而是通过在叶片区域内加入体积力源项，程序在生成网格时不需要考虑叶片的叶型，周向网格沿压气机轴线是轴对称的，叶片排前后缘轴向位置利用两条径向网格线给定（见图6）；第三：因为求解欧拉方程时不求解粘性力，近壁区网格没有必要加密，因此采用的计算网格比较均匀，有利于数值计算的收敛。同时由于稀疏网格在显式计算格式中对应于大的计算时间步长，大幅度加快了数值计算的速度。图6转静子附近的计算网格3.2基本方程程序采用的流体控制方程是带有源项的三维非定常欧拉方程在绝对柱坐标系下的积分形式，在通道和叶片轴向间隙区等没有叶片作用的区域源项为零。(3.1)其中：ρ为密度；p为压力；e为内能；H为总焓；Vx、Vr、Vθ为轴向、径向、周向分速度；U为叶排旋转周向速度；nx、nr、nθ分别为x、r、θ三个坐标方向的单位矢量；Fx、Fr、Fθ分别为体积力源项在轴向、径向、周向的分量。3.3计算格式求解欧拉方程所用的计算格式与Denton采用的Scree格式相同。在空间采用有限体积积分，变量储存在有限体积的八个角点上，而所求参数的变化量在体积的中心：(3.2)此处为位于有限体积中心的守恒量（密度、轴向动量、径向动量、周向动量矩和总能量）的变化量。FLUX为流经有限体积表面的通量（流量，三个方向的动量和能量）。在计算通量时，假设各流动参数在有限体积表面呈线性分布，因此此格式的空间离散精度为二阶精度。时间采用一侧后向差分，其精度视所用的项数，可以是二阶或三阶以上。在本程序中使用非标准的三阶格式，即：(3.3)计算时需要储存前两步时间的残差，由于本程序所要求的网格密度不大，因此造成的额外的内存要求不足为道。所求得的即为五个守恒变量之一在有限体中心的残差，或时间变化量，将其均分到有限体的八个角点上（线性假设），即可得到有关守恒变量随时间的变化。进口边界条件：进口给定总温、总压和气流流入角度，进口静压利用公式（3.4）三点外差计算。(3.4)其中系数α在0.05～0.2之间取值。图7进口静压外差示意图出口边界条件：出口给定静压，静压沿径向满足简单径向平衡条件。程序的流程图程序的主要流程包括生成网格、读入进、出口气动参数和叶片几何参数，建立初始流场，求解控制方程和最后计算结果输出。主程序流程图见图8。图8主程序流程图求解控制方程的SETFLUX子程序是程序的主体，主要功能包括：给定进、出口边界条件，构造叶片排内代替叶片力的体积力源项，分别求解连续方程、动量方程和能量方程，并对流场进行光顺。SETFLUX子程序的主要流程见图9。

图9SETFLUX子程序流程图

第四章叶片通道周向平均体积力模型建立体积力源项模化方法可以千差百异，如根据叶片基元级特性推算出体积力，这种方法在已知试验特性的情况下可能获得很好的结果，但在压气机方案设计阶段无法使用，而且获得叶片基元级特性需要在每一排叶片后沿叶高安装大量温度、压力测头，尤其在多级压气机试验中实现非常困难。另一种方法是利用一系列CFD计算归纳出体积力与当地动力学项的函数关系，这种方法需要进行大量的CFD计算，使用起来比较麻烦，将耗费大量的计算时间，而且总结出的函数关系在不同的压气机之间并不通用。本文采用物理上较为直观的模化方法，基于叶型的基本设计参数得到体积力，所需参数与通流设计所要求的参数相同，因而这种模化方法特别适合于设计阶段对各种方案抗畸变能力进行评估。4.1周向平均体积力模型原理周向平均体积力模型建立在无穷多叶片假设的基础上，所谓无穷多叶片假设是将叶片通道内流动假设为局部轴对称，即在每个叶片通道内叶片力沿周向均匀地作用在流体上，在无穷多叶片假设的基础上我们才能将叶片力用周向平均的体积力模拟[7]。通常情况下进气畸变的扰动都属于大波长扰动，其波长远大于叶片排栅距，因此研究畸变场时，忽略叶片通道内细小流动结构，着重强调叶片排对非定常三维非均匀流的响应，在此前提下引入无穷多叶片或无穷大稠度的假设是合理的。图10叶片力周向平均后转化为气流的体积力源项叶片通道周向平均体积力模型的基本思想是将叶片对气流的作用力在叶片槽道内沿周向平均，平均后的叶片力不再以叶片力的形式存在，而是以欧拉方程中体积力源项的形式存在，其大小就是沿槽道周向平均了的叶片力，见图10。在求解叶片排区域时并不把真实的叶片放入流场中，叶片对气流的作用通过周向均布的体积力源项来实现。所要构建的体积力源项应该像真实的叶片一样随着当地的流动参数的变化而变化。需要指出的是，周向平均的尺度只是限于叶片通道。因此虽然对应于叶片通道的周向范围，流动参数是均匀的；但是在大于叶片通道的尺度上，体积力仍允许有变化。叶片对周向流动畸变的响应就是由于这些变化而产生。4.2叶片力的分解首先将叶片对气流的作用力进行分解，叶片对气流中的力可以分解成两部分，一部分垂直于叶片表面，记做F压力；一部分平行于叶片表面，记做F摩擦。这样分解成二部分表达的好处是，每一部分都可以根据自身明确的物理意义进行表达。F压力表示在无粘情况下叶片对气流的压力，它垂直于叶片表面，实现叶片排与流体间功的交换。平行于叶片表面的力F摩擦是由叶片表面的粘性剪切层产生的，这一部分体积力产生损失，同时也对流体做功。实践证明在一般情况下叶片压力要远大于粘性力，即粘性力的估计精度相对而言对计算结果的精度的影响较小。同时粘性力和表面压力都还可以分解成轴向、周向和径向三个分量。图11为二维叶片力分解示意图，三维情况同理。图11二维叶片力分解示意图4.3叶片对气流切向摩擦力的模化气流流过叶片表面由于粘性作用受到叶片表面摩擦产生损失，这部分损失我们可以用简单的壁面摩擦阻力系数来模化。在空气动力学中壁面摩擦阻力系数的定义式是：（4.1）图12壁面阻力系数定义即面积为S的平板有速度为V∞、密度ρ的气流流过，平板对气流的阻力为Xf,就可以推算出平板的壁面阻力系数Cf。如果Cf已知，就可以推算出气流流过平板所受到的摩擦阻力。只要不发生显著的脱体现象，曲面情形的摩擦阻力和平板情景差不多[8]。因此平面计算结果可以近似应用到叶片摩阻计算中。我们可以把叶片近似看成一个两面都有气流流过的平板，如果已知叶片表面的壁面摩擦阻力系数，则气流流过叶片所受到的阻力就可以计算出来。叶片表面的摩擦阻力系数可以用平板的壁面摩擦阻力系数近似计算[8]，层流流动状态下根据附面层布拉休斯解可以推导出层流时平板壁面阻力系数的计算公式：适用于Rel＜5×105～106,其中：,紊流时：适用于5×105＜Rel＜5×107,适用于5×107＜Rel＜5×109l为叶片弦长,ν为气体运动粘性系数。叶片表面的摩擦阻力系数可以参考上面公式来选取，但因为上述公式是在零压力梯度下推导出来的，而在真实的压气机中叶片是工作在很大的逆压力梯度下，而且随着压气机工作状态的不同逆压力梯度的程度也不同，所以叶片表面的摩擦阻力系数的选取还要根据不同的设计水平和经验来确定。叶片的摩擦力应包括叶盆和叶背两个面，见图13，因此一个叶片的摩擦阻力应该等于：（4.2）图13叶片对气流粘性摩擦力的模化示意图然后要将叶片的摩擦力转化为体积力加到每一个网格控制体上，实际计算时叶片槽道内沿轴向和径向都会有几个网格，单个网格控制体上的体积力应该等于：（4.3）其中：β叶型角——叶片的叶型角、SABT网格单元体子午面的面积、Z叶片数、IMM1周向总的网格数。然后根据网格当地的叶片构造角和叶片倾角将粘性摩擦力沿径向、周向和轴向三个方向分解——。4.4叶片对气流法向压力的模化在不考虑粘性作用的情况下，叶片对气流的作用力应该垂直于叶片表面，因此在叶片几何角度已知的情况下，如果给定叶片力中的一个分量，其他两个分量就可以根据叶片几何角求出，见图14。图14叶片法向压力的三个分量（4.4）如何确定叶片压力周向分量？需要利用角动量守恒方程首先求出叶片对气流总的周向力Fθ，这个总的周向力是粘性摩擦力和无粘压力的周向分量的总合，4.3节中已经列出了粘性摩擦力的周向分量计算方法，用总的周向力减去周向粘性摩擦力，就可以得出叶片法向压力周向分量。总的周向分量Fθ可以根据周向角动量方程推导出来。图15用角动量方程求叶片力周向分量Fθ示意图从角动量方程：m为控制体中流体的质量，r为半径。即对于定常流流过控制体的质量不随时间变化，即=0则为控制体的平均密度，可以用控制体八个节点上密度值平均，为控制体体积。其中：所以：上面的公式适用于压气机流路半径变化不大的情况，对于流路变化大的情况应用下式：（4.5）其中为子午速度，为图15所示的网格长度。对于非定常情况计算公式为：（4.6）（4.5）和（4.6）式中其它参数都可以由欧拉方程求解器直接计算，只有Vθ出口需要给定，Vθ出口可以利用公式（4.7）计算：（4.7）公式（4.7）中的仍可由欧拉求解器直接计算，β叶型角作为程序输入数据也是已知的，只有落后角需要进行模化，影响落后角的因素有很多，比如叶型、攻角、马赫数等，准确地确定落后角可以大大地提高计算精度。落后角的计算见4.6节。叶片对气流总的周向力Fθ确定后，利用以下公式计算叶片无粘压力周向分量:（4.8）然后再根据（4.4）式求出无粘压力另外两个分量、。这样叶片粘性摩擦力和无粘压力都模化成体积力了。4.5堵塞判断随着出口反压的逐渐降低叶片排的某些截面会发生堵塞，发生堵塞截面的损失将大幅增长，如果这时仍用4.3节中的叶片摩擦阻力系数方法来计算损失就会偏小，因此在堵塞发生时程序就会加大摩擦阻力系数来计算堵塞状态下的损失，堵塞时叶片表面的摩擦阻力系数根据经验来选取。体积力模型因为计算时没有加入真实叶片，无法考虑叶片本身对压气机通道造成的堵塞，计算的堵塞流量偏大，因此必须通过别的方法来判断叶片排各截面的气流是否发生堵塞。准确判断叶片各截面是否发生堵塞，对损失和最大流量的计算十分重要的。本文采用轴向质量通量最大值的方法来判断堵塞[6]。所谓轴向质量通量最大值是指叶片某截面可能通过的最大流量，轴向质量通量最大值确定后，将它与叶片截面进口实际流量比较，如果叶片截面进口实际流量大于或等于轴向质量通量最大值，就说明该截面已经处于堵塞状态。叶片某个截面可以通过的轴向质量通量最大值可以根据质量连续原则推导出来。根据质量连续，叶栅进口处流量等于叶栅喉道处流量，即：图16堵塞判断模型示意图（4.9）其中：、、分别为叶片截面进口处的密度、轴向速度和面积的轴向投影；Prel、Trel叶片喉道处的相对总压和总温（静子为绝对总压和总温）A*为有效喉道面积喉道处的流量函数叶片堵塞实际上是叶片喉道处流速达到音速，即喉道处=1.0，因此叶片某截面可以通过的最大流量应该等于：忽略喉道上游的任何损失，喉道处的相对总温和总压等于叶片进口的值，即,，可得到（4.10）其中F相当于流函数，,是马赫数为1.0时的值，=1.281。式（1-1）中喉道面积A*可以根据叶片的真实叶型计算，也可以根据经验公式（4.11）近似计算。图17叶型喉道面积计算示意图（4.11）轴向质量通量最大值确定后，将它与叶片截面进口实际流量比较，如果叶片截面进口实际流量大于或等于轴向质量通量最大值，就说明该截面已经处于堵塞状态，气流损失将大幅增长。4.6落后角计算（4.7）式中Vθ出口的计算需要用到叶片落后角，确定落后角一般要依靠长期经验的总结，本文采用落后角经验计算公式为：(4.12）式中δ为落后角，为最大挠度位置，为出口叶型角，b/t为稠度，θ叶型弯角,为修正系数，，为速度修正系数，为轮毂比修正系数，为出口气流角修正系数，是强制修正系数。

第五章进气畸变条件下风扇定常与非定常流场数值模拟为验证计算模型，对某大涵道比发动机的风扇转子进行多种进口畸变场的数值模拟，首先计算均匀进气状态不同转速条件下风扇转子流场及特性，并与试验结果相对比。然后分别进行了稳态总压、总温畸变和动态温度脉冲与总压周期振荡条件下风扇畸变流场数值模拟。对模型模拟畸变流场、预估风扇畸变特性及非定常计算能力做出了全面考核。5.1算例风扇转子的主要气动参数所有算例中均采用某航空大涵道比发动机风扇/增压级部件的风扇转子，该转子有完整的设计参数和均匀进气与稳态总压畸变条件下风扇特性试验结果，为检验程序提供了可靠依据。表1为风扇转子主要气动参数。表1风扇主要气动参数（摘自文献[9]）叶片数32片展弦比2.6轮毂比0.35压比1.64效率0.91图18中指定了计算区域中的四个轴向位置，将在这些位置上提取气流参数用来进行试验结果与计算结果的对比。图18计算区域5.2计算网格子午面网格周向网格图19计算网格，叶片前后缘用黑实线表示计算网格为图19所示的轴对称网格，轴向121个网格，周向60个网格，径向30个网格，叶片排内沿轴向取9个网格，叶片排沿整个一周总的网格节点数为16200个（9*30*60）。计算网格沿周向和径向都是等距的，沿轴向也尽量保证网格的均匀变化，因此网格的均匀性和正交性都能得到很好的保证，这种规则网格有利于数值计算的收敛，对应于大的计算时间步长，大幅度加快了数值计算的速度。5.3均匀进气条件下数值模拟在计算畸变流场之前，先进行了均匀进气条件下风扇特性计算，用以验证程序计算准确性。共算了六个转速下的等转速特性，图20为计算结果与试验结果对比图。从整体来看计算风扇特性与试验结果相符，计算的流量、压比和效率值与试验结果十分符合，每条等转速线相应的最大质量流量点与试验值基本相符。=1.0和=0.958转速近喘点压比特性线趋势与试验对比可以看出，靠近不稳定边界时压比增加的趋势并未减缓，而试验结果则出现了压比增幅明显减缓的趋势，这是因为靠近不稳定边界时风扇的某些截面可能已经出现分离或旋转失速，风扇的增压能力下降，目前的计算模型还不能反映这一现象，因此如果想在近喘点得到更好计算结果，还要进一步改进近喘点模型。计算的效率-流量特性与试验结果也相当符合，除了效率的绝对值与试验十分接近外，特性线的形状也很相似，最高效率点的位置基本与试验重合。在计算模型中没有加入堵塞模型时，越靠近堵塞状态计算的效率越高，堵塞模型加入后计算的效率特性线才出现拐头，更加接近真实情况，所以认为本文采用的堵塞模型能够较为准确地判断堵塞，堵塞后的选取的摩擦阻力系数也比较合理。(a)风扇压比-流量特性(b)风扇流量-效率特性图20均匀进气条件下计算的风扇特性与试验结果对比由图21可以看出，用本文计算的落后角与商用全三维粘性程序计算的落后角在最高效率点符合得非常好，接近喘点状态两者的落后角也很相似，只在根部差别略大，接近堵点状态计算的落后角略偏小，总的来看本模型中经验公式计算的落后角在大部分状态下精度还是可以满足要求的。(a)近堵点(b)最高效率点(c)近喘点图21计算的风扇转子落后角沿叶高分布与全三维粘流程序计算结果对比(a)近堵点(b)最高效率点(c)近喘点图22计算的风扇转子压比沿叶高分布与全三维粘流程序计算结果对比由图22可见在最高效率点和近喘点计算的压比沿径向分布和商用全三维粘性程序计算计算结果在趋势上是一致的，在近堵点本文计算的结果偏高，因为商用程序计算风扇特性时近堵点压比比试验要低很多，而本文计算的近堵点压比与试验结果相差很小，所以认为本文计算的近堵点压比沿叶高分布应该与真实情况更接近。图23风扇转子中计算的相对马赫数分布，风扇=1.0设计转速、流量371公斤/秒工作，模型在进口处捕获一道激波图24风扇转子叶尖区域进口到出口计算的静压和相对气流角分布。进口是第一个方块符号，出口是最后一个。结果显示激波是在前缘三个网格点被捕获在=1.0转速情况下转子上相对马赫数的分布如图23所示。马赫数沿前缘的分布显示了拟激波的存在。由于沿轴向方向转子上只有9个网格分布，从所得马赫数等值线很难知道在转子前缘是否会出现激波。但叶尖的相对气流角分布和静压分布都绘出时，激波的存在变得明显。如图24所示，在此区域内相对气流角变化不大，而相对静压突增（在前三点以内）。捕捉拟激波的重要性就在于，它是跨音叶片排中产生压升和损失的主要机制。捕捉激波的现象与Gong的计算结果非常一致。通过以上对均匀进口流场数值模拟结果的分析和与试验结果的对比，可以认为本文所采用的计算模型能够很好地刻画风扇的相关流场，并能以较高精度预估风扇转子的特性曲线。5.4稳态进气畸变条件下模拟计算5.4.1稳态压力畸变稳态压力畸变是发动机进口最常见的一种畸变形式，飞机大攻角飞行、侧风、进气道与压气机不匹配都可能引发进口稳态总压畸变，它也是压气机畸变试验中必做的一个科目。本节将计算风扇转子在不同转速和畸变指数下的畸变流场，计算风扇畸变特性并与均匀进气计算结果和试验结果相比较，详细地分析总压畸变流场的流动结构，总结其流动特点，验证模型模拟畸变流场的能力。5.4.1.1给定进口畸变场风扇稳态压力畸变试验时在距离风扇2米位置处安装畸变栅产生进口总压畸变，在风扇进口安装压力测量装置测量总压不均匀度，畸变栅的形式见文献[10]。因试验报告中未提供风扇进口的畸变图谱，模拟计算时进口畸变场低总压区根据畸变栅的形状确定在146度的范围内，因畸变栅插入较深，近似认为畸变类型为周向畸变。分别计算了三个状态的畸变流场，每个状态对应的转速和畸变指数见表2。表2三个稳态总压畸变状态的转速和畸变参数状态相对换算转速畸变指数(%)低压区总压（pa）11.034.9593198.9121.04.693748.7830.9583.80595006.83低压区总压计算根据下面公式推算。式中：——畸变指数——低压区总压平均值——进口总压平均值进口总压畸变场云图和总压沿周向分布见图25。图25模拟计算时采用的计算域进口总压畸变场形式5.4.1.2计算结果与试验结果对比计算的畸变特性、均匀特性和试验特性都画在图26中，由图26可见当进口存在总压畸变时风扇性能在明显下降，流量降低、压比减小，随着转速的降低这种差距逐渐减小，模拟计算很好地反映出了这种趋势。计算结果与试验符合得相当好，各个转速下流量和压比的降低幅度都与试验十分相近，准确地反映出不同转速下进口总压畸变对风扇性能的衰减，说明程序可以定量地反映进气畸变条件下风扇特性的变化。图26均匀进气和进口总压畸变条件下计算特性与试验对比5.4.1.3进口总压畸变流场流动结构分析压气机内畸变流场最大特点就是流动畸变和压气机的耦合使得流动被重新分布。风扇上游总压和静压的分布见图27，在畸变区总压低于非畸变区，静压沿整个周向变化不大。流量系数见图28，在畸变区流量系数明显偏低。由图29可见流量系数在风扇进口的分布相对其在风扇上游位置的分布有了改变，在风扇进口流量系数的非均匀性被消减，流动系数的分布型也发生了变化，这种分布现象的原因是，在低流动系数区域压气机的压升高于高流动系数区域，这种扰动向上传播改变了流量系数的分布。风扇出口的压比分布见图30，其形状和图28中风扇进口流量系数分布正好相反，因为在低流量系数的区域风扇增压能力强，压比高，因此在畸变区的压比高于非畸变区，这与平行压气机理论是一致的。图31显示了计算域进口和风扇出口总压对比，显然通过风扇后总压不均匀性被削弱了。正是因为总压低的区域通过压气机后产生较高的压升，所以总压畸变通过压气机后总是被衰减。图27风扇上游换算总压、静压沿周向分布图28风扇上游流量系数沿周向分布图29风扇上游和转子进口处流量系数图30转子出口压比周向分布图31风扇上游和风扇出口换算总压沿周向分布(a)风扇上游(b)风扇进口(c)风扇出口图32不同轴向位置处的总压云图(a)风扇上游(b)风扇出口图33不同轴向位置处的总温云图图32中列出了不同轴向位置处总压云图，可以看出流动具有很强的三维性，在叶尖处低压区产生了较大相位移动，而叶根处则基本没有相位变化。图33显示即使进口处的畸变仅仅是总压畸变，但由于加功量的周向分布因总压畸变而产生不均匀性，出口处的总温分布也会呈现明显的不均匀性。风扇转静子在进口总压畸变条件下流场数值模拟计算为验证程序计算多排叶片的能力，在风扇转子叶片后加上一排静子叶片进行两排叶片的总压畸变场模拟计算，计算选取的相对换算转速为=1.0，进口畸变场与5.4.1.1节中状态2相同。图34为转静子联算的子午网格图。图35为不同轴向位置处总压沿周向的分布，可以清楚地看到总压畸变沿轴向的衰减，总压畸变在经过转子后被衰减，经过静子叶片后又进一步衰减，说明静子叶片也具有衰减总压畸变的能力。图36为风扇不同轴向位置处总压云图，由图可见畸变流的三维性非常明显，风扇上游的周向畸变在经过风扇转子后，叶中的低总压区已经几乎被抹平了，但叶尖和叶根的低压区还很明显，而且有一定的相位移动。气流流过静子后叶尖区的低压区也被抹平了，低压区只存在在叶根区附近。图34风扇转静子联算子午网格图图35风扇不同轴向位置处换算总压沿周向分布(a)风扇上游(b)转子出口(c)静子出口图36风扇不同轴向位置处总压云图5.4.1.5稳定性预估工作是目前流体力学界公认的难点问题，在相关的文献[7][15]中可以看到目前判断流场稳定性的一种方法是在流场中加入一个波长为若干个动叶栅距的扰动，然后根据扰动的发展来判断失速发生，目前这种方法还在研究阶段，只能够定性地判断失速。作者试图利用这种方法来判断失速，但利用目前的模型几乎算不出失速现象。分析认为如果要利用这种方法来预测失速，目前的模型中还需进一步改进，因为在近失速点损失大幅度提高，叶片增压能力下降，其特性线在小于某一流量时压比会出现拐头，模型要能够算出这种现象才能够利用加扰动的方法来判断稳定性。由于时间关系模型没有加入模拟失速的部分，所以无法用加扰动的办法来判断失速。作者提出用更简单和常用的办法来判断失速，即利用叶栅的扩散因子是否超标来判断失速。扩散因子是表征叶栅逆压力梯度的参数，叶栅只有在逆压力梯度足够大时才能发生分离失速，所以用扩散因子来判断失速是有一定的物理依据的。程序中采用的扩散因子定义式为：其中W1、W2为转子叶片进、出口相对速度，为扭速，为稠度。首先计算了均匀进气条件下设计转速近喘点343kg/s流量时沿风扇叶片扩散因子最大值为0.63。然后计算了相同流量下畸变指数为4.6%的总压畸变进口时风扇转子整个叶高最大扩散因子沿周向的分布，计算结果见图48。相同流量下进口总压畸变流情况下在畸变区扩散因子最大超过了0.7，如果假设风扇可承受的最大扩散因子是0.63的话，那么在进口有总压畸变的情况下，风扇就不可能在相同的流量下还保持稳定，它的稳定流量比均匀进口情况下要小。显然在进口总压畸变时风扇的压比降低、流量减小，而且稳定边界点的流量也比均匀情况下要小，因此其稳定裕度较均匀进口情况将有所降低。图37进口总压畸变条件下风扇转子扩散因子沿周向分布因为风扇试验时没有试到喘点，本文无法验证利用扩散因子的方法是否能够定量地判断失速点。但作者认为在周向平均体积力模型加三维欧拉方程的计算模型中利用扩散因子来判断失速还是可行的，因为叶栅通道内的周向平均体积力模型与S2程序十分类似，它和S2程序一样不考虑叶片通道内部周向的不均匀，所不同的是S2程序是沿压气机整个环形通道的轴对称，而本文程序是沿单个叶片通道的轴对称。在S2程序中扩散因子是衡量叶栅逆压力梯度最有效的参数，其有效性经多年的研究和使用所验证。稳态总温畸变运输机使用反推力装置，战斗机吸入其他机群排出的气体，舰载机吸入蒸气弹射装置产生的热蒸气，直升机接近地面悬停时燃气的吸入等都引起温度畸变问题，此外低压压气机进口的总压畸变往往衍生高压压气机的总温畸变。为验证模型模拟稳态总温畸变能力，在风扇转子上游加入一周向稳态温度畸变进行模拟计算，畸变区扇形角121度，畸变区温升87.41k,畸变指数18.2%。温度畸变指数定义公式如下：式中：——畸变指数——高温区总温平均值——进口总温平均值风扇上游畸变场总温云图和总温沿周向分布见图38。总温云图总温沿周向分布图38稳态总温畸变模拟计算时风扇上游总温畸变场(a)风扇上游(b)风扇进口(c)风扇出口图39风扇不同轴向位置处的总温云图图39列出了风扇不同轴向位置处总温云图，可以看出温度畸变沿流程的变化，在风扇上游标准方波形的周向畸变在风扇进口处高温区沿周向向低温区扩散，风扇出口截面周向总温畸变已经演变成了周向和径向的混合畸变。图40为风扇上游和出口总温对比，可以看出通过风扇后总温不均匀性有少量的衰减，畸变区沿周向略有偏移。畸变区的高总温被衰减可以这样解释：畸变区较高的总温与非畸变区相比因物理转速相同而对应着较低的换算转速，因此加功量小，温升小，缓解了畸变区总温高的趋势。图41为风扇不同轴向位置处总压云图，图42为风扇上游和出口总压对比，由这两张图可以看到畸变区的总压升确实比非畸变区小，所以虽然进口处的畸变仅仅是总温畸变，出口处的总压分布也呈现明显的不均匀性，解释了风扇进口温度畸变诱导高压进口压力畸变的现象。图40风扇上游和风扇出口总温对比(a)风扇上游(b)风扇进口(c)风扇出口图41不同轴向位置总压云图图42风扇上游和风扇出口总压周向分布图43风扇上游和风扇进口轴向速度对比图43显示了从风扇上游位置到风扇进口轴向速度的变化，低总温区的流动被加速，而高总温区基本没有。平行压气机理论可以很好的解释这一现象：具有高总温来流的压气机部分必须在低质量流量下工作，而其它部分则须在高质量流量下工作，畸变区的增压能力低于非畸变区，而出口静压几乎是均匀的，因此流动畸变区域在进入风扇之前必须相对非畸变区减速增压。模拟计算与Gong的计算结果十分相似，而且与平行压气机理论定性的分析一致，可以认为本模型能够很好地反映稳态温度畸变与压气机的相互作用。5.5动态进气畸变条件下模拟计算本模型的一大特点就是能够模拟压气机对动态畸变的响应，为了验证模型这一功能我们将模拟风扇转子在进口总压周期振荡和进口温度快速瞬变条件下的畸变流场。5.5.1进口周期压力振荡动态总压畸变起源于进气道的外部和内部两个方面，地面旋风、阵风的吸入属进气道外部起因，进气道工作不稳定、嗡鸣、以及激波与附面层干扰等则属内部起因，它们可以产生周期性压力震荡和脉冲；战斗机打机关炮时也会产生周期压力振荡。文献[11]中对周期性振荡型压力畸变的研究认为，叶轮机经受这类动态畸变后，各级工作不再匹配。引起级的工作失配的主要机理是压力脉动的幅值衰减和相位错移，在各个频率下的幅值衰减和相位差关系就是动态系统的幅频特性和相频特性，其变化特点是频率越高幅值比越小、相位差越大。本文模拟计算了三种不同频率、相同幅值的周期压力振荡下风扇的响应，来了解该风扇对周期压力振荡的频率特性，并检验程序计算非定常流动的能力。模拟计算设定总压振荡畸变区扇形角98度，畸变区总压振荡幅值12650Pa,最大畸变指数4.621%，压力振荡频率分别为10Hz，20Hz，40Hz，计算时在风扇上游边界强加一个周期振荡的压力波，图44为上游边界处某一时刻的压力云图和三种振荡频率下上游边界畸变区总压随时间的变化。图44模拟计算时风扇上游边界总压云图及其畸变区总压随时间的变化通过对计算结果进行分析得出如下结论：风扇对这种动态的压力振荡表现出了良好的幅频特性，见图45，随着频率的增大，风扇出口压力振动幅值与进口之比即幅值比减小，这个结果与文献[11]的实验和计算结果是一致的。这说明在压力振荡幅值相等的情况下，低频的压力振荡经过压气机后衰减的幅度要小，对压气机性能的影响更大。模拟计算的相频特性不明显，风扇进出口压力脉动的相位差很小，频率的影响也不大，分析其原因认为，主要是因为计算的叶排数太少，一排叶片对压力振荡的影响力不足以使振荡波的相位发生明显的变化。图46在风扇特性图上绘制出不同压力振荡频率下风扇工作点运动轨迹，由图可见低频压力振荡时工作点偏离均匀进气工作点最多，流量的变化范围也最大，低进口总压时刻风扇的工作点越接近稳定工作边界，而不同的频率下的压比变化范围相差不是很大。和上面的结论一致，在压力振荡幅值相同的情况下，频率越低对风扇的影响越大，越偏危险。数值模拟的结果可见，模型很好地反映出风扇对周期压力振荡的幅频特性，证明了程序具有模拟动态畸变的能力。图45动态压力振荡下风扇的幅频特性图46不同压力振荡频率下风扇工作点运动轨迹5.5.2进口温度快速瞬变飞机在发射武器后发动机吸入燃气，使压气机进口温度随时间突升。一架垂直起落飞机，当喷口转到垂直向下排气时，通过地面反射，发动机进口的下半部分吸入燃气，也是一个温度瞬变的过程。它们的变化近似于一个三角形波，一般用温升率表示瞬变特征。对涡轮喷气发动机的试验结果表明[11]，压气机对温度瞬变的稳定性响应主要取决于温升率(△T/△t)。当压气机转速一定时，进口温度突升，就使换算转速降低，而压比变化不大，在特性图上压气机工作点会沿着等压比线向左移动，接近稳定工作边界。图47模拟计算时风扇上游边界总温云图及其畸变区总温随时间的变化模拟计算分别在压气机进口加入相同温升不同温升率的2组温度脉冲，目的是对比不同温升率对压气机的影响。根据目前掌握资料武器发射时温升率基本上在500-2500k/s的范围内[12]，因此计算的温升率在这个范围内选取：1000k/s、2000k/s，温升为26.7k，高温区仍在周向108度的范围内。进口总温分布云图和畸变区总温随时间的变化见图47。图48为在这种动态总温情况下风扇流量随时间的变化，可以看到随着温度的上升流量减小，但温度上升到最大值时流量并未下到最小值，而是有一个滞后，如果假设温升率为1000k/s的温度脉冲周期为T1，温升率为2000k/s的温度脉冲周期为T2，温升率为1000k/s时温度最大值与流量最小值之间约相差0.17T1，而温升率为2000k/s情况下两者之间相差约0.29T2，说明温升率越高，流量相应的滞后越严重。同时图45中可以看到，虽然两个温度脉冲的温升是相同的，但低温升率的温度脉冲引起的流量下降更严重些，这一点很好地验证了相同温升情况下低温升率更加危险的结论。图49为风扇压比随时间的变化，总体上温度脉冲对风扇压比的影响不大，风扇压比略有降低，最低压比点也相对于最高温度点有一定滞后，其规律与流量变化相同，所不同的是不同的温升率情况下压比的降低几乎相同。图50在风扇特性图上描绘了发生温度脉冲过程中风扇工作点的运动轨迹，风扇进口局部区域温度的动态升高使风扇换算转速降低，工作点沿着等压比线向左侧移动，使喘振裕度降低，极限情况是有可能跨过喘振边界发生喘振。同时温升率越低工作点移动得越远，越危险，这些现象与以往试验和理论研究的结论都是一致的。