風(fēng)扇壓氣機進氣畸變周向平均體積力模型三維非定常流

STYLEREF北航章目錄北京航空航天大學(xué)碩士學(xué)位論文PAGEivPAGEiii摘要壓氣機進口的均勻流動條件是工程設(shè)計所期望的，但是當(dāng)把發(fā)動機安裝接到飛機上時，在飛行過程中，就不能保證在各種狀態(tài)下壓氣機或風(fēng)扇的進口流場都是均勻的了。壓氣機進口畸變帶來的直接嚴重后果是壓氣機的穩(wěn)定裕度下降和增壓比、效率降低，因而加強進氣畸變對壓氣機性能和穩(wěn)定性影響的研究就顯得非常重要。本文采用三維歐拉方程加體積力模型的方法，發(fā)展了一套三維非定常進氣畸變快速數(shù)值模擬程序，藉此可以實現(xiàn)風(fēng)扇/壓氣機在進氣畸變條件下的全環(huán)面三維非定常數(shù)值模擬。該方法從物理基本原理出發(fā)，能夠較好地反映工作于非均勻流場中的壓縮系統(tǒng)以及葉片通道內(nèi)部的非定常三維流動本質(zhì)。程序用極小的計算機資源即可以實現(xiàn)對風(fēng)扇全環(huán)面三維非定?；兞鲌龅哪M。利用該程序模擬計算了某航空發(fā)動機風(fēng)扇進口穩(wěn)態(tài)總壓畸變條件下特性和流場結(jié)構(gòu)，并與試驗結(jié)果進行了對比，結(jié)果表明所能夠以較高的精度預(yù)估氣流畸變對壓氣機/風(fēng)扇特性的影響。另外還計算了穩(wěn)態(tài)總溫畸變，動態(tài)總壓振蕩和溫度脈沖對該風(fēng)扇流動結(jié)構(gòu)的影響，計算結(jié)果與目前畸變試驗和理論研究結(jié)論一致，說明該程序可以很好地模擬穩(wěn)態(tài)和動態(tài)總溫、總壓畸變對壓氣機流動結(jié)構(gòu)和特性的影響。關(guān)鍵詞：風(fēng)扇/壓氣機；進氣畸變；周向平均體積力模型；三維非定常流

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本文工作內(nèi)容……………….8第二章壓縮系統(tǒng)的建?！?02.1管道中的流動……………….102.2葉片排中的流動…………….10第三章三維非定常歐拉方程求解器……………….123.1網(wǎng)格生成…………………….123.2基本方程………………….123.3計算格式…….…………….133.4程序的流程圖………….………………..14葉片通道周向平均體積力模型建立…….………………….184.1周向平均體積力模型原理………………..184.2葉片力的分解…………….………………194.3葉片對氣流切向摩擦力的?；?94.4葉片對氣流法向壓力的模化……………..214.5堵塞判斷……..…………..234.6落后角計算………..………………….….25第五章進氣畸變條件下風(fēng)扇定常與非定常流場數(shù)值模擬………….265.1算例風(fēng)扇轉(zhuǎn)子的主要氣動參數(shù)…………..265.2計算網(wǎng)格…………………..265.3均勻進氣條件下數(shù)值模擬………………..275.4穩(wěn)態(tài)進氣畸變條件下模擬計算……….…315.4.1穩(wěn)態(tài)壓力畸變………………..…..315.4.1.1給定進口畸變場…………………….….315.4.1.2計算結(jié)果與試驗結(jié)果對比……………..325.4.1.3進口總壓畸變流場流動結(jié)構(gòu)分析……..335.4.1.4風(fēng)扇轉(zhuǎn)靜子在進口總壓畸變條件下流場數(shù)值模擬計算……………..355.4.1.5利用擴散因子進行穩(wěn)定性預(yù)估的初步研究…………..365.4.2穩(wěn)態(tài)總溫畸變…………………...375.5動態(tài)進氣畸變條件下模擬計算………………….…….405.5.1進口周期壓力振蕩………………405.5.2進口溫度快速瞬變……………….425.6結(jié)論……………………….45研究總結(jié)與展望……………………….47參考文獻……………….49攻讀碩士期間發(fā)表的論文及所取得的研究成果………………..……51致謝…………………….52PAGE52PAGE51第一章緒論1.1選題意義壓氣機進口的均勻流動條件是工程設(shè)計所希望的，當(dāng)壓氣機裝在壓氣機試驗臺上進行試驗時，試驗臺的進氣裝置保證了壓氣機進口流場的均勻性，當(dāng)發(fā)動機在地面試車臺上試車時，氣流經(jīng)過工藝進氣道而流入發(fā)動機的壓氣機，這時壓氣機的進口流場也是相當(dāng)均勻的。但是當(dāng)把發(fā)動機安裝到飛機上時，在飛行過程中就不能保證在各種狀態(tài)下壓氣機或風(fēng)扇的進口流場都是均勻的了。飛機爬升、大偏航角、側(cè)風(fēng)和進氣道不起動工況等都會造成壓氣機進口處總壓、總溫、速度和流向不均，從而引起壓氣機性能的惡化。進氣畸變對壓氣機工作的影響有兩個方面，一是空氣動力學(xué)方面的影響，另一方面是氣體力學(xué)方面的影響。前者的影響使壓氣機的性能和穩(wěn)定性變壞，后者的作用是葉片上非定常力引起機械振動。本文主要研究壓氣機對進氣畸變的氣動響應(yīng)。在航空發(fā)動機性能評估中，一個影響其性能的重要的因子就是進口流場畸變，而進氣畸變對發(fā)動機性能影響最直接、最重要的表現(xiàn)就是對壓氣機性能和穩(wěn)定性的影響。壓氣機進口畸變帶來的直接嚴重后果是壓氣機的增壓比、效率降低和穩(wěn)定裕度下降。因而加強進氣畸變對壓氣機性能影響的研究就顯得非常重要?，F(xiàn)在高性能航空燃氣渦輪發(fā)動機的研制和發(fā)展過程中，為了實現(xiàn)高性能的發(fā)動機設(shè)計指標(biāo)，壓氣機設(shè)計者總是希望能夠?qū)兣c發(fā)動機相互作用的物理本質(zhì)進行深入的認識和理解，壓氣機設(shè)計必須同時進行性能和抗畸變能力的綜合考慮，設(shè)計者必須掌握能夠預(yù)測不同設(shè)計參數(shù)情況下壓氣機抗畸變能力的計算工具。研究進氣畸變對壓氣機性能影響主要有兩條途徑，第一條途徑是試驗研究，在壓氣機運轉(zhuǎn)過程中通過在壓氣機前方安裝畸變發(fā)生器來模擬各種狀態(tài)下壓氣機進口畸變場，然后對進口畸變場和壓氣機性能進行詳細測量，獲得各種畸變條件下壓氣機特性，但模擬整個使用條件范圍內(nèi)的全科目試驗成本非常昂貴，而且試驗研究本身迄今為止對加深物理問題認識的幫助甚少。在國際上基本上也是通過大量試驗來獲取壓氣機設(shè)計參數(shù)與抗畸變能力的工程關(guān)聯(lián)，但是有關(guān)的試驗數(shù)據(jù)無法在公開的文獻中獲得，另一方面有限的研究結(jié)果都屬于敏感情報，有關(guān)各國都極力保密。近些年來計算流體力學(xué)的發(fā)展為畸變問題研究開辟了另一個途徑——數(shù)值模擬，我們希望通過數(shù)值模擬來獲得壓氣機在進氣畸變條件下的特性，在設(shè)計階段就對壓氣機的抗畸變能力有所了解，利用數(shù)值計算來指導(dǎo)設(shè)計，提高壓氣機設(shè)計的成功率。進口畸變與多級壓氣機相互作用問題是一個復(fù)雜的三維、非定常、非線性問題，數(shù)值模擬方法既要保證抓住畸變問題的物理本質(zhì)，其計算量又要能夠被工程應(yīng)用所接受。在過去半個多世紀壓氣機發(fā)展的過程中，有許多研究者通過采用各種近似簡化假設(shè)，發(fā)展可許多預(yù)測壓氣機畸變性能的近似模型，但因其計算體系廣泛采用了一維、二維流動假設(shè)，無法反映真實流動的三維性，而壓氣機眾多設(shè)計新理念是基于三維流動過程控制和設(shè)計，因此這些計算系統(tǒng)不能反映出這些新的設(shè)計思想的效果，無法評估發(fā)動機新的設(shè)計理念對發(fā)動機抗畸變能力的影響。接下來很容易想到的是用非定常的RANS來直接模擬，它幾乎不依賴于經(jīng)驗，能夠計算葉片槽道內(nèi)、外很細微的非定常流動結(jié)構(gòu)，但用非定常RANS數(shù)值模擬來研究畸變流場面臨的一個主要困境就是非定常擾動尺度范圍過大。壓氣機在畸變條件下各種不同尺度的非定常擾動都將會顯示出來，一方面，進口氣流畸變是以壓氣機直徑為特征尺寸的大尺度非定常擾動，求解畸變流場就要求計算壓氣機整個環(huán)面所用通道中的非定常流動。另一方面，為了正確估算粘性作用，還需要求解比邊界層尺度還小的擾動，對于高雷諾數(shù)小尺度擾動的長度可以達到幾十個微米的量級，精細網(wǎng)格和全環(huán)面多通道計算使得計算網(wǎng)格數(shù)量驚人的巨大，它對于計算機資源的高要求就不止包含內(nèi)存，也包含CPU計算時間。因此利用非定常RANS求解畸變流場需要極為苛刻的計算機硬件能力，其計算所費機時是目前工程設(shè)計工作無法接受的。本課題的目的是尋找一個即能真實地反映壓氣機內(nèi)畸變流動的三維、非定常物理本質(zhì)，又適合于工程應(yīng)用的快速壓氣機畸變數(shù)值模擬方法，可以利用它來探索畸變對壓氣機性能影響的機理，研究不同設(shè)計參數(shù)對壓氣機抗畸變能力的影響，在壓氣機設(shè)計階段就可以對其畸變性能作出預(yù)估，對壓氣機的性能和抗畸變能力做出綜合考慮，降低設(shè)計風(fēng)險，提高壓氣機設(shè)計成功率。1.2課題背景在進氣畸變數(shù)值模擬研究方面，已經(jīng)比較成熟的技術(shù)是利用一維非線性模型加平行壓氣機理論來研究有畸變進氣的壓縮系統(tǒng)，美國的Arnold工程發(fā)展中心和俄羅斯中央航空研究院都利用該方法發(fā)展了非常完善的計算軟件。所謂一維模型，也稱stage-by-stage模型。它將整個壓氣機按照級沿著軸線方向劃分成為一個一個的計算單元（一維）。針對每一個計算單元，列出氣流運動的控制方程，通過對所有單元控制方程的聯(lián)立求解分析，獲得壓氣機性能。而平行壓氣機理論，是將壓氣機沿周向分解成許多獨立段，每一段構(gòu)成一個子壓氣機，不同的子壓氣機之間相互并行工作，具有相同的出口邊界條件?；冞M氣的流動條件則表現(xiàn)為：在不同的子壓氣機進口具有不同的進氣條件，所有并行的子壓氣機都按無畸變的壓氣機特性工作，他們之間通過出口邊界條件聯(lián)系起來[1]。平行壓氣機理論可以成功地解釋一些物理現(xiàn)象，但也具有很大的局限性，該模型不能考慮到不同的進口流場區(qū)域在壓氣機通道內(nèi)流動過程中的氣流摻混，而且一維模型由于對非定常流動的模擬方式過于簡單，無法反映壓縮系統(tǒng)內(nèi)部三維的流動本質(zhì)，直接影響到預(yù)測結(jié)果的可靠性。圖1平行壓氣機理論示意圖早在1978年，Greitzer和Strand在研究施加在周向渦流上的畸變[2]時就發(fā)現(xiàn)畸變流場表現(xiàn)為三維特性，為了確定進氣畸變對壓氣機的影響，需要深入認識進氣畸變下壓氣機流場細節(jié)。但是由于計算技術(shù)和試驗技術(shù)發(fā)展的限制，針對畸變進氣條件下發(fā)動機三維流場的細節(jié)認識及其理論計算研究工作一直進展不大。近年來數(shù)值方法和計算機科學(xué)的發(fā)展使計算流體力學(xué)有可能對畸變問題進行模擬研究。似乎我們可以通過時間精確的三維RANS方程來求解壓氣機內(nèi)的非定?；兞鲌觯珕栴}在于畸變問題是個擾動尺度相差幾個數(shù)量級的非定常問題，大尺度的擾動可以達到米的量級，想要捕捉到這種擾動在壓氣機中的傳播和衰減就要求能夠解壓氣機的整個環(huán)形通道，而不是象我們現(xiàn)在通常所做的非定常計算時每個葉排只求解幾個通道。另一方面所要求解的最小尺度比邊界層的尺度還要小，這樣才可以正確的估算粘性作用，因此計算網(wǎng)格尺度必須非常小，對于高雷諾數(shù)可以達到幾十個微米的量級。精細的網(wǎng)格加上整臺壓氣機全環(huán)面所有通道都要進行非定常計算，這樣的計算量對計算機硬件能力和計算時間的需求是工程設(shè)計工作所無法接受的。例如，由Unno等人（2001）做的模擬風(fēng)扇—吊架結(jié)構(gòu)用了大約500萬個網(wǎng)格在大型的并行（166個CPU）計算機上計算[3]。目前僅有美國NASA和英國帝國理工學(xué)院兩家有此計算能力并在進行相關(guān)真實流動情況下全壓氣機/發(fā)動機流場的非定常流動計算，但即使使用高密度計算機群，其海量的計算使一個簡單算例的計算也要進行幾周甚至幾個月。這樣的工作周期對工程分析和設(shè)計而言是遠水解不了近渴，而且計算結(jié)果是詳細的非定常流場數(shù)據(jù)，如何從中提取有用的相關(guān)物理信息，尚需復(fù)雜的后處理分析工作[4]。針對低維模型無法反映流動的三維特性，而直接求解RANS方程又因計算機硬件和計算時間的限制而無法應(yīng)用于工程研究，1985年Adamcyzk給出了他的通道平均多級流動模型[5]。這種方法直接求解壓氣機中的三維流動，但并不直接計算葉片對氣流的作用力，而是通過體積力模型來模擬，因為不需要直接求解葉片力，計算量比RANS要小很多。通道平均多級流動模型主要思想是忽略葉片通道周向不均勻性，物理上可以理解為周向有無窮多個葉片，這樣葉柵通道內(nèi)的某些流動信息將丟失，但因為進氣畸變這種非定常擾動的尺度遠遠大于葉柵柵距，我們關(guān)注的是葉片排對非定常三維非均勻流的響應(yīng)，葉柵通道內(nèi)部小尺度的非定常響應(yīng)往往不是我們關(guān)心的重點，因此這種忽略小尺度擾動的通道平均模型應(yīng)用于進氣畸變模擬是可以接受的。葉片對氣流的作用可以用激盤和半激盤模型或者通道平均體積力模型來實現(xiàn)。1997年Joo和Hynes基于圓盤激勵模型[6]，1998年Gong等人基于通道平均體力模型[7]來模擬計算進氣畸變，這兩種方法都需要通過均勻的軸對稱特性來推斷在非均勻擾動作用下的葉片排性能。通道平均體力模型比激盤模型的優(yōu)勢為可以計算葉片排內(nèi)部的流動，體積力源項可以沿徑向和軸向有不同的分布，在數(shù)值解法上體積力模型擺脫了激勵圓盤模型中難應(yīng)付的人工邊界。盡管如此體積力模型仍然存在一個問題，體力的源項在空間的分布基于均勻的軸對稱特性，但在在非設(shè)計點比如接近失速時，葉片負荷將會經(jīng)歷大的弦向和徑向的再分布，這就還需要一些特殊的處理[4]。1998年Gong等將通道平均體積力模型應(yīng)用于畸變這種非定常、三維問題的研究中[7]，他的作法是不把葉片放入流場中參與計算，用無粘的三維非定常歐拉方程求解管道流動，葉片的作用以非定常體積力形式加入到方程源項中。這種方法的關(guān)鍵技術(shù)是如何準確地建立體積力源項與當(dāng)?shù)貏恿W(xué)量的關(guān)聯(lián)，尤其是在非設(shè)計點。Gong分別采用兩種方法來建立體積力源項，對于低速大展弦比壓氣機仍采用通過均勻的軸對稱特性來獲得體積力分布，即通過每排葉片各截面的壓力升特性等參數(shù)換算出體積力源項，但是源項沿葉片弦向的分布則需要依賴經(jīng)驗來獲得。而對于跨音級高速小展弦比的壓氣機他選擇采用另一種方法，這種方法首先利用CFD計算或測試獲得壓氣機的三維流場，對流場進行周向流動推力平均后得到軸對稱流場，然后對軸對稱流場應(yīng)用動量方程求出體積力，經(jīng)過一系列的不同狀態(tài)下的CFD計算后，總結(jié)出體積力源項與當(dāng)?shù)貏恿W(xué)參數(shù)的函數(shù)關(guān)系，這種方法的優(yōu)點是對經(jīng)驗的依賴較前一種方法要少。Gong利用這個模型模擬計算了NASA的單級跨音壓氣機Stage35的穩(wěn)態(tài)總壓、總溫和動態(tài)總壓畸變流場，成功地捕捉到了高速壓氣機獨有的二種特性：阻塞和激波，驗證了周向平均體積力模型能夠?qū)簹鈾C畸變引起的流動再分配及畸變在壓氣機中的傳遞有很好的刻畫。圖2和圖3為Gong計算總壓和總溫畸變條件下壓氣機內(nèi)速度系數(shù)沿流程的變化，可以明顯地看到畸變流場的三維流動特性。圖2Gong利用體積力模型計算的穩(wěn)態(tài)總壓畸變條件下速度系數(shù)沿流程的變化圖3Gong利用體積力模型計算的穩(wěn)態(tài)總溫畸變條件下速度系數(shù)沿流程的變化2001年Xu等人提出了粘性體力模型[4]，并利用該模型來模擬畸變流場。Xu的粘性體積力模型與Gong的周向平均體積力模型的不同是他把葉片放入流場中參與計算，無粘葉片力直接由方程計算出來的，而粘性影響才作為當(dāng)?shù)卦错椩诖志W(wǎng)格上進行模擬的，被模化的只是葉片槽道中的粘性力作用，因此這種模型相對于粘性力和非粘性葉片力都需要?；哪Ｐ蛯?jīng)驗的依賴程度要小，而且它有能力求解在葉柵通道內(nèi)的詳細流動，可以解由于葉片掃過而產(chǎn)生的擾動，并可以很精確的給出葉片的非定常力。但是因為方程求解時要計算葉片非粘性力作用，計算耗時要比體積力模型高一個量級。Xu發(fā)現(xiàn)一般情況下粘性力相對于總的葉片力是很小的，它的非定常擾動則更小，只需要一個很簡單的阻力系數(shù)模型就可以很好地模擬時均和動態(tài)粘性作用，但對于有分離的情況，阻力系數(shù)的變化將會變的很大，這時就需要更詳細的模擬。Xu在研究模擬粘性體積力的方法時發(fā)現(xiàn)如果直接計算無粘的葉片力，而模擬粘性力將會十分有利，因為無粘非定常力通常對擾動有很強的動力學(xué)響應(yīng)，因此模擬無粘非定常力非常困難。還有，由于無粘葉片力是總的非定常葉片力的主要部分，模擬整體葉片力的任何誤差的影響將會比單獨模擬粘性影響重要得的多。換句話說，與周向平均的體積力模型相比粘性體力模型的模擬可以承受更大的誤差[4]。圖4為Xu利用粗網(wǎng)格粘性體積力模型計算的進氣畸變壓力場與精細網(wǎng)格RANS計算計算結(jié)果對比，兩者計算靜壓和總壓場十分相近。粘性體積力模型RANS圖4Xu利用粗網(wǎng)格粘性體積力模型計算的進氣畸變壓力場與精細網(wǎng)格RANS計算結(jié)果對比比較Gong的周向平均體積力模型和Xu的粘性體力模型，周向平均體積力模型由于在管道內(nèi)只去解通道平均的流動，沒有把葉片幾何形狀考慮進去（全部組合進源項中），因此可以采用更規(guī)則的網(wǎng)格，在計算速度上比粘性體力模型高效一個量級，更適合于工程應(yīng)用。然而，通過通道平均葉柵通道內(nèi)的流動有關(guān)信息將會丟失，并且不可能去研究葉片對非定常擾動的響應(yīng)。粘性體力模型是較周向平均體積力模型進行了改進的，由于其在葉柵通道內(nèi)解非粘性流動，因此它有能力在葉柵通道內(nèi)解非定常流動并可以很精確的給出葉片的非定常力。綜上所述，利用數(shù)值方法來模擬非定常進氣畸變流動的計算模型主要包括高度簡化的低維模型如平行壓氣機理論模型和激盤模型還有周向平均體積力模型、粘性體積力模型及RANS直接數(shù)值模擬等。2001年Xu等人為了解決葉輪機械內(nèi)部非定常流動數(shù)值模擬所面對的非定常尺度跨度大的問題，提出了一個對各種尺度的非定常流動都適用的通用框架[4]，框架示意圖見圖5。框架最上端是低維模型，它對計算機資源的需求最小，但對流動做的簡化最多，計算中被忽略的流動結(jié)構(gòu)需要依靠經(jīng)驗來建模，簡化模型越多，對經(jīng)驗的依賴也就越大，經(jīng)驗總結(jié)的正確性決定了它計算的精度?？蚣茏畹锥耸侵苯訑?shù)值模擬，它對計算機資源有極高的要求，所有的流動細節(jié)都通過控制守恒方程直接計算，幾乎不依賴于經(jīng)驗?？蚣苤虚g的方法則介于兩者之間。這個框架更像是廣義大渦模擬，在某一預(yù)設(shè)的臨界值以上尺度的非定常流通過流動的控制守恒方程直接計算，低于臨界值的小尺度擾動將會通過不同的模型——通常是由經(jīng)驗、數(shù)值、解析方法或者這些方法的結(jié)合來模擬。當(dāng)我們要對某一個特殊應(yīng)用選擇一個合適的模型時，就要在經(jīng)驗和計算機資源之間做一個折衷。圖5各種計算模型對計算機資源和經(jīng)驗的需求層次1.3本文工作內(nèi)容本文工作的目的是建立一個適用于工程研究應(yīng)用的畸變流場計算程序，要求能夠真實地反映壓氣機內(nèi)畸變流動的三維、非定常物理本質(zhì)，并且計算速度快，在壓氣機設(shè)計階段結(jié)可以對其畸變特性做出預(yù)估，為方案篩選提供依據(jù)。目前高性能壓氣機設(shè)計更加注重三維性，寬弦、掠型、端彎等設(shè)計先進的設(shè)計手段在設(shè)計中頻繁應(yīng)用，簡單的低維模型模型不能模擬出這些三維因素對壓氣機抗畸變能力的影響。而對于畸變問題直接求解三維RANS方程、大渦模擬(LES)和直接數(shù)值模擬(DNS)則需要耗費大量計算機資源和時間，在目前的情況下還不能應(yīng)用于工程設(shè)計階段。鑒于這種情況，本文采用了周向平均的體積模型來模擬計算畸變流場，周向平均的體積模型居于非定常數(shù)值模擬框架的中間，它的優(yōu)點是對經(jīng)驗和計算機資源的依賴都相對較少，可以采用更規(guī)則的網(wǎng)格，計算速度快，在模型建立合理的前提下，數(shù)值更準確，更適合于工程應(yīng)用。而且，因為在整個壓氣機流場內(nèi)求解的是三維非定常流體控制方程，因此能夠真實地反映壓氣機中畸變流場的三維、非定常流動本質(zhì)，為研究進氣畸變和壓氣機相互作用的物理機理提供了可能。本文的主要工作是構(gòu)造一個三維非定常歐拉方程加葉片通道周向平均體積力模型的進氣畸變條件下壓氣機流場計算程序。其中主要包括三部分內(nèi)容：第一部分是三維非定常歐拉方程求解器的調(diào)試。本文所使用的歐拉方程求解器具有三階時間精度和兩階空間精度的歐拉方程求解器，可以計算雙涵道多級壓氣機非定常流場，作者主要負責(zé)程序的閱讀、修改和調(diào)試及驗算工作。第二部分是建立周向平均體積力模型，這是本文工作的重點，體積力模型建立的好壞直接關(guān)系到程序計算正確性及精度。與Gong相同我們都采用了周向平均的體積力模型來模擬葉片排對氣流的作用。但是體積力模型的構(gòu)造方法與Gong不同，Gong的體積力模型需要預(yù)先輸入各個葉片排的基元特性然后推算出體積力，或者根據(jù)一系列CFD計算結(jié)果歸納出體積力與當(dāng)?shù)貏恿W(xué)參數(shù)的函數(shù)關(guān)系，而本文的體積力的模化方法是將葉片對氣體的作用力分解為粘性摩擦力和無粘壓力兩部分分別模擬，粘性摩擦力用簡單的壁面摩擦阻力系數(shù)來?；瑹o粘壓力則通過氣流流過葉片后的周向動量變化，利用角動量守恒方程計算。同時因為計算時沒有加入真實的葉片排，無法考慮葉片排對壓氣機通道的堵塞，在程序中加入了堵塞判斷模型，當(dāng)判斷出葉片某一截面發(fā)生堵塞時，程序?qū)⒋蠓忍岣吣Σ磷枇ο禂?shù)，使近堵點的計算更接近真實情況。本文的體積力構(gòu)造方法主要基于葉型基本設(shè)計參數(shù)，所需參數(shù)與通流設(shè)計所要求的參數(shù)略同，因而特別適合于設(shè)計階段的方案評估。第三部分為程序驗算，選取某航空發(fā)動機單排風(fēng)扇轉(zhuǎn)子作為計算對象，計算了該轉(zhuǎn)子均勻進氣條件下相對換算轉(zhuǎn)速=0.6～1.0轉(zhuǎn)速下風(fēng)扇特性，并與試驗進行了對比；還分別計算了進口總壓穩(wěn)態(tài)畸變條件下風(fēng)扇流場特性，也與試驗進行了對比；計算了進口總溫穩(wěn)態(tài)畸變條件下風(fēng)扇流場、進口動態(tài)溫度脈沖條件下風(fēng)扇流場、進口動態(tài)壓力振蕩條件下風(fēng)扇流場，充分驗證了程序計算壓氣機在各種畸變類型包括動態(tài)畸變流場的能力。

第二章壓縮系統(tǒng)建模對于進口流場畸變這個三維、非定常問題本文選擇用非定常三維歐拉方程加周向平均體積力模型來建立計算模型。該計算模型與Gong所采用的模型[7]相似：在整個壓氣機通道中利用三維無粘非定常歐拉方程來求解管道流動，不把葉片放入流場中參與計算，葉片作用以體積力的形式加入到方程源項中。2.1管道中的流動氣流流經(jīng)壓氣機通道時，受到葉片對它的法向壓力和切向摩擦力、上下流道的摩擦力和氣流本身的粘性力，在正常的流動情況下，粘性力主要發(fā)生在附面層中，主流區(qū)中粘性力相對很小，因此本文忽略主流區(qū)的粘性作用，同時壓氣機上下壁面對氣流的摩擦力也忽略不計，在這種情況下壓氣機前后的無葉片段和葉片排間隙段都可以直接利用無粘的歐拉方程求解。在壓氣機進口管道、出口管道和葉片排間隙中采用三維非定常歐拉方程可以如此表示：（2.1）其中：ρ為密度；p為壓力；e為內(nèi)能；Vx、Vr、Vθ為軸向、徑向、周向速度2.2葉片排中的流動在葉片排區(qū)域內(nèi)葉片對氣流的法向壓力和切向摩擦力也不直接計算，而是用一定的方法進行?；；蟮娜~片力以體積力的形式被加入到歐拉方程的源項中，因此在葉片排區(qū)域求解的是帶有體積力源項的歐拉方程，方程具體形式如下：（2.2）其中Fx、Fr、Fθ分別為體積力源項在軸向、徑向、周向的分量，為體積力對氣流做的功。

第三章三維非定常歐拉方程求解器3.1網(wǎng)格生成程序有自己的網(wǎng)格生成程序，在已知壓氣機流道尺寸的情況下可以自動生成3維H形網(wǎng)格，程序所采用的網(wǎng)格有三個特點：第一，因為畸變流場模擬要求計算壓氣機整個環(huán)形通道內(nèi)的流場，因此其網(wǎng)格也必須沿整個壓氣機通道生成，而不是象通常三維計算時只生成一個或幾個葉片通道區(qū)域的網(wǎng)格；第二，程序在計算葉片區(qū)作用時葉片排對氣流的作用力不是直接計算而是通過在葉片區(qū)域內(nèi)加入體積力源項，程序在生成網(wǎng)格時不需要考慮葉片的葉型，周向網(wǎng)格沿壓氣機軸線是軸對稱的，葉片排前后緣軸向位置利用兩條徑向網(wǎng)格線給定（見圖6）；第三：因為求解歐拉方程時不求解粘性力，近壁區(qū)網(wǎng)格沒有必要加密，因此采用的計算網(wǎng)格比較均勻，有利于數(shù)值計算的收斂。同時由于稀疏網(wǎng)格在顯式計算格式中對應(yīng)于大的計算時間步長，大幅度加快了數(shù)值計算的速度。圖6轉(zhuǎn)靜子附近的計算網(wǎng)格3.2基本方程程序采用的流體控制方程是帶有源項的三維非定常歐拉方程在絕對柱坐標(biāo)系下的積分形式，在通道和葉片軸向間隙區(qū)等沒有葉片作用的區(qū)域源項為零。(3.1)其中：ρ為密度；p為壓力；e為內(nèi)能；H為總焓；Vx、Vr、Vθ為軸向、徑向、周向分速度；U為葉排旋轉(zhuǎn)周向速度；nx、nr、nθ分別為x、r、θ三個坐標(biāo)方向的單位矢量；Fx、Fr、Fθ分別為體積力源項在軸向、徑向、周向的分量。3.3計算格式求解歐拉方程所用的計算格式與Denton采用的Scree格式相同。在空間采用有限體積積分，變量儲存在有限體積的八個角點上，而所求參數(shù)的變化量在體積的中心：(3.2)此處為位于有限體積中心的守恒量（密度、軸向動量、徑向動量、周向動量矩和總能量）的變化量。FLUX為流經(jīng)有限體積表面的通量（流量，三個方向的動量和能量）。在計算通量時，假設(shè)各流動參數(shù)在有限體積表面呈線性分布，因此此格式的空間離散精度為二階精度。時間采用一側(cè)后向差分，其精度視所用的項數(shù)，可以是二階或三階以上。在本程序中使用非標(biāo)準的三階格式，即：(3.3)計算時需要儲存前兩步時間的殘差，由于本程序所要求的網(wǎng)格密度不大，因此造成的額外的內(nèi)存要求不足為道。所求得的即為五個守恒變量之一在有限體中心的殘差，或時間變化量，將其均分到有限體的八個角點上（線性假設(shè)），即可得到有關(guān)守恒變量隨時間的變化。進口邊界條件：進口給定總溫、總壓和氣流流入角度，進口靜壓利用公式（3.4）三點外差計算。(3.4)其中系數(shù)α在0.05～0.2之間取值。圖7進口靜壓外差示意圖出口邊界條件：出口給定靜壓，靜壓沿徑向滿足簡單徑向平衡條件。程序的流程圖程序的主要流程包括生成網(wǎng)格、讀入進、出口氣動參數(shù)和葉片幾何參數(shù)，建立初始流場，求解控制方程和最后計算結(jié)果輸出。主程序流程圖見圖8。圖8主程序流程圖求解控制方程的SETFLUX子程序是程序的主體，主要功能包括：給定進、出口邊界條件，構(gòu)造葉片排內(nèi)代替葉片力的體積力源項，分別求解連續(xù)方程、動量方程和能量方程，并對流場進行光順。SETFLUX子程序的主要流程見圖9。

圖9SETFLUX子程序流程圖

第四章葉片通道周向平均體積力模型建立體積力源項模化方法可以千差百異，如根據(jù)葉片基元級特性推算出體積力，這種方法在已知試驗特性的情況下可能獲得很好的結(jié)果，但在壓氣機方案設(shè)計階段無法使用，而且獲得葉片基元級特性需要在每一排葉片后沿葉高安裝大量溫度、壓力測頭，尤其在多級壓氣機試驗中實現(xiàn)非常困難。另一種方法是利用一系列CFD計算歸納出體積力與當(dāng)?shù)貏恿W(xué)項的函數(shù)關(guān)系，這種方法需要進行大量的CFD計算，使用起來比較麻煩，將耗費大量的計算時間，而且總結(jié)出的函數(shù)關(guān)系在不同的壓氣機之間并不通用。本文采用物理上較為直觀的模化方法，基于葉型的基本設(shè)計參數(shù)得到體積力，所需參數(shù)與通流設(shè)計所要求的參數(shù)相同，因而這種模化方法特別適合于設(shè)計階段對各種方案抗畸變能力進行評估。4.1周向平均體積力模型原理周向平均體積力模型建立在無窮多葉片假設(shè)的基礎(chǔ)上，所謂無窮多葉片假設(shè)是將葉片通道內(nèi)流動假設(shè)為局部軸對稱，即在每個葉片通道內(nèi)葉片力沿周向均勻地作用在流體上，在無窮多葉片假設(shè)的基礎(chǔ)上我們才能將葉片力用周向平均的體積力模擬[7]。通常情況下進氣畸變的擾動都屬于大波長擾動，其波長遠大于葉片排柵距，因此研究畸變場時，忽略葉片通道內(nèi)細小流動結(jié)構(gòu)，著重強調(diào)葉片排對非定常三維非均勻流的響應(yīng)，在此前提下引入無窮多葉片或無窮大稠度的假設(shè)是合理的。圖10葉片力周向平均后轉(zhuǎn)化為氣流的體積力源項葉片通道周向平均體積力模型的基本思想是將葉片對氣流的作用力在葉片槽道內(nèi)沿周向平均，平均后的葉片力不再以葉片力的形式存在，而是以歐拉方程中體積力源項的形式存在，其大小就是沿槽道周向平均了的葉片力，見圖10。在求解葉片排區(qū)域時并不把真實的葉片放入流場中，葉片對氣流的作用通過周向均布的體積力源項來實現(xiàn)。所要構(gòu)建的體積力源項應(yīng)該像真實的葉片一樣隨著當(dāng)?shù)氐牧鲃訁?shù)的變化而變化。需要指出的是，周向平均的尺度只是限于葉片通道。因此雖然對應(yīng)于葉片通道的周向范圍，流動參數(shù)是均勻的；但是在大于葉片通道的尺度上，體積力仍允許有變化。葉片對周向流動畸變的響應(yīng)就是由于這些變化而產(chǎn)生。4.2葉片力的分解首先將葉片對氣流的作用力進行分解，葉片對氣流中的力可以分解成兩部分，一部分垂直于葉片表面，記做F壓力；一部分平行于葉片表面，記做F摩擦。這樣分解成二部分表達的好處是，每一部分都可以根據(jù)自身明確的物理意義進行表達。F壓力表示在無粘情況下葉片對氣流的壓力，它垂直于葉片表面，實現(xiàn)葉片排與流體間功的交換。平行于葉片表面的力F摩擦是由葉片表面的粘性剪切層產(chǎn)生的，這一部分體積力產(chǎn)生損失，同時也對流體做功。實踐證明在一般情況下葉片壓力要遠大于粘性力，即粘性力的估計精度相對而言對計算結(jié)果的精度的影響較小。同時粘性力和表面壓力都還可以分解成軸向、周向和徑向三個分量。圖11為二維葉片力分解示意圖，三維情況同理。圖11二維葉片力分解示意圖4.3葉片對氣流切向摩擦力的?；瘹饬髁鬟^葉片表面由于粘性作用受到葉片表面摩擦產(chǎn)生損失，這部分損失我們可以用簡單的壁面摩擦阻力系數(shù)來?；Ｔ诳諝鈩恿W(xué)中壁面摩擦阻力系數(shù)的定義式是：（4.1）圖12壁面阻力系數(shù)定義即面積為S的平板有速度為V∞、密度ρ的氣流流過，平板對氣流的阻力為Xf,就可以推算出平板的壁面阻力系數(shù)Cf。如果Cf已知，就可以推算出氣流流過平板所受到的摩擦阻力。只要不發(fā)生顯著的脫體現(xiàn)象，曲面情形的摩擦阻力和平板情景差不多[8]。因此平面計算結(jié)果可以近似應(yīng)用到葉片摩阻計算中。我們可以把葉片近似看成一個兩面都有氣流流過的平板，如果已知葉片表面的壁面摩擦阻力系數(shù)，則氣流流過葉片所受到的阻力就可以計算出來。葉片表面的摩擦阻力系數(shù)可以用平板的壁面摩擦阻力系數(shù)近似計算[8]，層流流動狀態(tài)下根據(jù)附面層布拉休斯解可以推導(dǎo)出層流時平板壁面阻力系數(shù)的計算公式：適用于Rel＜5×105～106,其中：,紊流時：適用于5×105＜Rel＜5×107,適用于5×107＜Rel＜5×109l為葉片弦長,ν為氣體運動粘性系數(shù)。葉片表面的摩擦阻力系數(shù)可以參考上面公式來選取，但因為上述公式是在零壓力梯度下推導(dǎo)出來的，而在真實的壓氣機中葉片是工作在很大的逆壓力梯度下，而且隨著壓氣機工作狀態(tài)的不同逆壓力梯度的程度也不同，所以葉片表面的摩擦阻力系數(shù)的選取還要根據(jù)不同的設(shè)計水平和經(jīng)驗來確定。葉片的摩擦力應(yīng)包括葉盆和葉背兩個面，見圖13，因此一個葉片的摩擦阻力應(yīng)該等于：（4.2）圖13葉片對氣流粘性摩擦力的?；疽鈭D然后要將葉片的摩擦力轉(zhuǎn)化為體積力加到每一個網(wǎng)格控制體上，實際計算時葉片槽道內(nèi)沿軸向和徑向都會有幾個網(wǎng)格，單個網(wǎng)格控制體上的體積力應(yīng)該等于：（4.3）其中：β葉型角——葉片的葉型角、SABT網(wǎng)格單元體子午面的面積、Z葉片數(shù)、IMM1周向總的網(wǎng)格數(shù)。然后根據(jù)網(wǎng)格當(dāng)?shù)氐娜~片構(gòu)造角和葉片傾角將粘性摩擦力沿徑向、周向和軸向三個方向分解——。4.4葉片對氣流法向壓力的?；诓豢紤]粘性作用的情況下，葉片對氣流的作用力應(yīng)該垂直于葉片表面，因此在葉片幾何角度已知的情況下，如果給定葉片力中的一個分量，其他兩個分量就可以根據(jù)葉片幾何角求出，見圖14。圖14葉片法向壓力的三個分量（4.4）如何確定葉片壓力周向分量？需要利用角動量守恒方程首先求出葉片對氣流總的周向力Fθ，這個總的周向力是粘性摩擦力和無粘壓力的周向分量的總合，4.3節(jié)中已經(jīng)列出了粘性摩擦力的周向分量計算方法，用總的周向力減去周向粘性摩擦力，就可以得出葉片法向壓力周向分量。總的周向分量Fθ可以根據(jù)周向角動量方程推導(dǎo)出來。圖15用角動量方程求葉片力周向分量Fθ示意圖從角動量方程：m為控制體中流體的質(zhì)量，r為半徑。即對于定常流流過控制體的質(zhì)量不隨時間變化，即=0則為控制體的平均密度，可以用控制體八個節(jié)點上密度值平均，為控制體體積。其中：所以：上面的公式適用于壓氣機流路半徑變化不大的情況，對于流路變化大的情況應(yīng)用下式：（4.5）其中為子午速度，為圖15所示的網(wǎng)格長度。對于非定常情況計算公式為：（4.6）（4.5）和（4.6）式中其它參數(shù)都可以由歐拉方程求解器直接計算，只有Vθ出口需要給定，Vθ出口可以利用公式（4.7）計算：（4.7）公式（4.7）中的仍可由歐拉求解器直接計算，β葉型角作為程序輸入數(shù)據(jù)也是已知的，只有落后角需要進行模化，影響落后角的因素有很多，比如葉型、攻角、馬赫數(shù)等，準確地確定落后角可以大大地提高計算精度。落后角的計算見4.6節(jié)。葉片對氣流總的周向力Fθ確定后，利用以下公式計算葉片無粘壓力周向分量:（4.8）然后再根據(jù)（4.4）式求出無粘壓力另外兩個分量、。這樣葉片粘性摩擦力和無粘壓力都?；审w積力了。4.5堵塞判斷隨著出口反壓的逐漸降低葉片排的某些截面會發(fā)生堵塞，發(fā)生堵塞截面的損失將大幅增長，如果這時仍用4.3節(jié)中的葉片摩擦阻力系數(shù)方法來計算損失就會偏小，因此在堵塞發(fā)生時程序就會加大摩擦阻力系數(shù)來計算堵塞狀態(tài)下的損失，堵塞時葉片表面的摩擦阻力系數(shù)根據(jù)經(jīng)驗來選取。體積力模型因為計算時沒有加入真實葉片，無法考慮葉片本身對壓氣機通道造成的堵塞，計算的堵塞流量偏大，因此必須通過別的方法來判斷葉片排各截面的氣流是否發(fā)生堵塞。準確判斷葉片各截面是否發(fā)生堵塞，對損失和最大流量的計算十分重要的。本文采用軸向質(zhì)量通量最大值的方法來判斷堵塞[6]。所謂軸向質(zhì)量通量最大值是指葉片某截面可能通過的最大流量，軸向質(zhì)量通量最大值確定后，將它與葉片截面進口實際流量比較，如果葉片截面進口實際流量大于或等于軸向質(zhì)量通量最大值，就說明該截面已經(jīng)處于堵塞狀態(tài)。葉片某個截面可以通過的軸向質(zhì)量通量最大值可以根據(jù)質(zhì)量連續(xù)原則推導(dǎo)出來。根據(jù)質(zhì)量連續(xù)，葉柵進口處流量等于葉柵喉道處流量，即：圖16堵塞判斷模型示意圖（4.9）其中：、、分別為葉片截面進口處的密度、軸向速度和面積的軸向投影；Prel、Trel葉片喉道處的相對總壓和總溫（靜子為絕對總壓和總溫）A*為有效喉道面積喉道處的流量函數(shù)葉片堵塞實際上是葉片喉道處流速達到音速，即喉道處=1.0，因此葉片某截面可以通過的最大流量應(yīng)該等于：忽略喉道上游的任何損失，喉道處的相對總溫和總壓等于葉片進口的值，即,，可得到（4.10）其中F相當(dāng)于流函數(shù)，,是馬赫數(shù)為1.0時的值，=1.281。式（1-1）中喉道面積A*可以根據(jù)葉片的真實葉型計算，也可以根據(jù)經(jīng)驗公式（4.11）近似計算。圖17葉型喉道面積計算示意圖（4.11）軸向質(zhì)量通量最大值確定后，將它與葉片截面進口實際流量比較，如果葉片截面進口實際流量大于或等于軸向質(zhì)量通量最大值，就說明該截面已經(jīng)處于堵塞狀態(tài)，氣流損失將大幅增長。4.6落后角計算（4.7）式中Vθ出口的計算需要用到葉片落后角，確定落后角一般要依靠長期經(jīng)驗的總結(jié)，本文采用落后角經(jīng)驗計算公式為：(4.12）式中δ為落后角，為最大撓度位置，為出口葉型角，b/t為稠度，θ葉型彎角,為修正系數(shù)，，為速度修正系數(shù)，為輪轂比修正系數(shù)，為出口氣流角修正系數(shù)，是強制修正系數(shù)。

第五章進氣畸變條件下風(fēng)扇定常與非定常流場數(shù)值模擬為驗證計算模型，對某大涵道比發(fā)動機的風(fēng)扇轉(zhuǎn)子進行多種進口畸變場的數(shù)值模擬，首先計算均勻進氣狀態(tài)不同轉(zhuǎn)速條件下風(fēng)扇轉(zhuǎn)子流場及特性，并與試驗結(jié)果相對比。然后分別進行了穩(wěn)態(tài)總壓、總溫畸變和動態(tài)溫度脈沖與總壓周期振蕩條件下風(fēng)扇畸變流場數(shù)值模擬。對模型模擬畸變流場、預(yù)估風(fēng)扇畸變特性及非定常計算能力做出了全面考核。5.1算例風(fēng)扇轉(zhuǎn)子的主要氣動參數(shù)所有算例中均采用某航空大涵道比發(fā)動機風(fēng)扇/增壓級部件的風(fēng)扇轉(zhuǎn)子，該轉(zhuǎn)子有完整的設(shè)計參數(shù)和均勻進氣與穩(wěn)態(tài)總壓畸變條件下風(fēng)扇特性試驗結(jié)果，為檢驗程序提供了可靠依據(jù)。表1為風(fēng)扇轉(zhuǎn)子主要氣動參數(shù)。表1風(fēng)扇主要氣動參數(shù)（摘自文獻[9]）葉片數(shù)32片展弦比2.6輪轂比0.35壓比1.64效率0.91圖18中指定了計算區(qū)域中的四個軸向位置，將在這些位置上提取氣流參數(shù)用來進行試驗結(jié)果與計算結(jié)果的對比。圖18計算區(qū)域5.2計算網(wǎng)格子午面網(wǎng)格周向網(wǎng)格圖19計算網(wǎng)格，葉片前后緣用黑實線表示計算網(wǎng)格為圖19所示的軸對稱網(wǎng)格，軸向121個網(wǎng)格，周向60個網(wǎng)格，徑向30個網(wǎng)格，葉片排內(nèi)沿軸向取9個網(wǎng)格，葉片排沿整個一周總的網(wǎng)格節(jié)點數(shù)為16200個（9*30*60）。計算網(wǎng)格沿周向和徑向都是等距的，沿軸向也盡量保證網(wǎng)格的均勻變化，因此網(wǎng)格的均勻性和正交性都能得到很好的保證，這種規(guī)則網(wǎng)格有利于數(shù)值計算的收斂，對應(yīng)于大的計算時間步長，大幅度加快了數(shù)值計算的速度。5.3均勻進氣條件下數(shù)值模擬在計算畸變流場之前，先進行了均勻進氣條件下風(fēng)扇特性計算，用以驗證程序計算準確性。共算了六個轉(zhuǎn)速下的等轉(zhuǎn)速特性，圖20為計算結(jié)果與試驗結(jié)果對比圖。從整體來看計算風(fēng)扇特性與試驗結(jié)果相符，計算的流量、壓比和效率值與試驗結(jié)果十分符合，每條等轉(zhuǎn)速線相應(yīng)的最大質(zhì)量流量點與試驗值基本相符。=1.0和=0.958轉(zhuǎn)速近喘點壓比特性線趨勢與試驗對比可以看出，靠近不穩(wěn)定邊界時壓比增加的趨勢并未減緩，而試驗結(jié)果則出現(xiàn)了壓比增幅明顯減緩的趨勢，這是因為靠近不穩(wěn)定邊界時風(fēng)扇的某些截面可能已經(jīng)出現(xiàn)分離或旋轉(zhuǎn)失速，風(fēng)扇的增壓能力下降，目前的計算模型還不能反映這一現(xiàn)象，因此如果想在近喘點得到更好計算結(jié)果，還要進一步改進近喘點模型。計算的效率-流量特性與試驗結(jié)果也相當(dāng)符合，除了效率的絕對值與試驗十分接近外，特性線的形狀也很相似，最高效率點的位置基本與試驗重合。在計算模型中沒有加入堵塞模型時，越靠近堵塞狀態(tài)計算的效率越高，堵塞模型加入后計算的效率特性線才出現(xiàn)拐頭，更加接近真實情況，所以認為本文采用的堵塞模型能夠較為準確地判斷堵塞，堵塞后的選取的摩擦阻力系數(shù)也比較合理。(a)風(fēng)扇壓比-流量特性(b)風(fēng)扇流量-效率特性圖20均勻進氣條件下計算的風(fēng)扇特性與試驗結(jié)果對比由圖21可以看出，用本文計算的落后角與商用全三維粘性程序計算的落后角在最高效率點符合得非常好，接近喘點狀態(tài)兩者的落后角也很相似，只在根部差別略大，接近堵點狀態(tài)計算的落后角略偏小，總的來看本模型中經(jīng)驗公式計算的落后角在大部分狀態(tài)下精度還是可以滿足要求的。(a)近堵點(b)最高效率點(c)近喘點圖21計算的風(fēng)扇轉(zhuǎn)子落后角沿葉高分布與全三維粘流程序計算結(jié)果對比(a)近堵點(b)最高效率點(c)近喘點圖22計算的風(fēng)扇轉(zhuǎn)子壓比沿葉高分布與全三維粘流程序計算結(jié)果對比由圖22可見在最高效率點和近喘點計算的壓比沿徑向分布和商用全三維粘性程序計算計算結(jié)果在趨勢上是一致的，在近堵點本文計算的結(jié)果偏高，因為商用程序計算風(fēng)扇特性時近堵點壓比比試驗要低很多，而本文計算的近堵點壓比與試驗結(jié)果相差很小，所以認為本文計算的近堵點壓比沿葉高分布應(yīng)該與真實情況更接近。圖23風(fēng)扇轉(zhuǎn)子中計算的相對馬赫數(shù)分布，風(fēng)扇=1.0設(shè)計轉(zhuǎn)速、流量371公斤/秒工作，模型在進口處捕獲一道激波圖24風(fēng)扇轉(zhuǎn)子葉尖區(qū)域進口到出口計算的靜壓和相對氣流角分布。進口是第一個方塊符號，出口是最后一個。結(jié)果顯示激波是在前緣三個網(wǎng)格點被捕獲在=1.0轉(zhuǎn)速情況下轉(zhuǎn)子上相對馬赫數(shù)的分布如圖23所示。馬赫數(shù)沿前緣的分布顯示了擬激波的存在。由于沿軸向方向轉(zhuǎn)子上只有9個網(wǎng)格分布，從所得馬赫數(shù)等值線很難知道在轉(zhuǎn)子前緣是否會出現(xiàn)激波。但葉尖的相對氣流角分布和靜壓分布都繪出時，激波的存在變得明顯。如圖24所示，在此區(qū)域內(nèi)相對氣流角變化不大，而相對靜壓突增（在前三點以內(nèi)）。捕捉擬激波的重要性就在于，它是跨音葉片排中產(chǎn)生壓升和損失的主要機制。捕捉激波的現(xiàn)象與Gong的計算結(jié)果非常一致。通過以上對均勻進口流場數(shù)值模擬結(jié)果的分析和與試驗結(jié)果的對比，可以認為本文所采用的計算模型能夠很好地刻畫風(fēng)扇的相關(guān)流場，并能以較高精度預(yù)估風(fēng)扇轉(zhuǎn)子的特性曲線。5.4穩(wěn)態(tài)進氣畸變條件下模擬計算5.4.1穩(wěn)態(tài)壓力畸變穩(wěn)態(tài)壓力畸變是發(fā)動機進口最常見的一種畸變形式，飛機大攻角飛行、側(cè)風(fēng)、進氣道與壓氣機不匹配都可能引發(fā)進口穩(wěn)態(tài)總壓畸變，它也是壓氣機畸變試驗中必做的一個科目。本節(jié)將計算風(fēng)扇轉(zhuǎn)子在不同轉(zhuǎn)速和畸變指數(shù)下的畸變流場，計算風(fēng)扇畸變特性并與均勻進氣計算結(jié)果和試驗結(jié)果相比較，詳細地分析總壓畸變流場的流動結(jié)構(gòu)，總結(jié)其流動特點，驗證模型模擬畸變流場的能力。5.4.1.1給定進口畸變場風(fēng)扇穩(wěn)態(tài)壓力畸變試驗時在距離風(fēng)扇2米位置處安裝畸變柵產(chǎn)生進口總壓畸變，在風(fēng)扇進口安裝壓力測量裝置測量總壓不均勻度，畸變柵的形式見文獻[10]。因試驗報告中未提供風(fēng)扇進口的畸變圖譜，模擬計算時進口畸變場低總壓區(qū)根據(jù)畸變柵的形狀確定在146度的范圍內(nèi)，因畸變柵插入較深，近似認為畸變類型為周向畸變。分別計算了三個狀態(tài)的畸變流場，每個狀態(tài)對應(yīng)的轉(zhuǎn)速和畸變指數(shù)見表2。表2三個穩(wěn)態(tài)總壓畸變狀態(tài)的轉(zhuǎn)速和畸變參數(shù)狀態(tài)相對換算轉(zhuǎn)速畸變指數(shù)(%)低壓區(qū)總壓（pa）11.034.9593198.9121.04.693748.7830.9583.80595006.83低壓區(qū)總壓計算根據(jù)下面公式推算。式中：——畸變指數(shù)——低壓區(qū)總壓平均值——進口總壓平均值進口總壓畸變場云圖和總壓沿周向分布見圖25。圖25模擬計算時采用的計算域進口總壓畸變場形式5.4.1.2計算結(jié)果與試驗結(jié)果對比計算的畸變特性、均勻特性和試驗特性都畫在圖26中，由圖26可見當(dāng)進口存在總壓畸變時風(fēng)扇性能在明顯下降，流量降低、壓比減小，隨著轉(zhuǎn)速的降低這種差距逐漸減小，模擬計算很好地反映出了這種趨勢。計算結(jié)果與試驗符合得相當(dāng)好，各個轉(zhuǎn)速下流量和壓比的降低幅度都與試驗十分相近，準確地反映出不同轉(zhuǎn)速下進口總壓畸變對風(fēng)扇性能的衰減，說明程序可以定量地反映進氣畸變條件下風(fēng)扇特性的變化。圖26均勻進氣和進口總壓畸變條件下計算特性與試驗對比5.4.1.3進口總壓畸變流場流動結(jié)構(gòu)分析壓氣機內(nèi)畸變流場最大特點就是流動畸變和壓氣機的耦合使得流動被重新分布。風(fēng)扇上游總壓和靜壓的分布見圖27，在畸變區(qū)總壓低于非畸變區(qū)，靜壓沿整個周向變化不大。流量系數(shù)見圖28，在畸變區(qū)流量系數(shù)明顯偏低。由圖29可見流量系數(shù)在風(fēng)扇進口的分布相對其在風(fēng)扇上游位置的分布有了改變，在風(fēng)扇進口流量系數(shù)的非均勻性被消減，流動系數(shù)的分布型也發(fā)生了變化，這種分布現(xiàn)象的原因是，在低流動系數(shù)區(qū)域壓氣機的壓升高于高流動系數(shù)區(qū)域，這種擾動向上傳播改變了流量系數(shù)的分布。風(fēng)扇出口的壓比分布見圖30，其形狀和圖28中風(fēng)扇進口流量系數(shù)分布正好相反，因為在低流量系數(shù)的區(qū)域風(fēng)扇增壓能力強，壓比高，因此在畸變區(qū)的壓比高于非畸變區(qū)，這與平行壓氣機理論是一致的。圖31顯示了計算域進口和風(fēng)扇出口總壓對比，顯然通過風(fēng)扇后總壓不均勻性被削弱了。正是因為總壓低的區(qū)域通過壓氣機后產(chǎn)生較高的壓升，所以總壓畸變通過壓氣機后總是被衰減。圖27風(fēng)扇上游換算總壓、靜壓沿周向分布圖28風(fēng)扇上游流量系數(shù)沿周向分布圖29風(fēng)扇上游和轉(zhuǎn)子進口處流量系數(shù)圖30轉(zhuǎn)子出口壓比周向分布圖31風(fēng)扇上游和風(fēng)扇出口換算總壓沿周向分布(a)風(fēng)扇上游(b)風(fēng)扇進口(c)風(fēng)扇出口圖32不同軸向位置處的總壓云圖(a)風(fēng)扇上游(b)風(fēng)扇出口圖33不同軸向位置處的總溫云圖圖32中列出了不同軸向位置處總壓云圖，可以看出流動具有很強的三維性，在葉尖處低壓區(qū)產(chǎn)生了較大相位移動，而葉根處則基本沒有相位變化。圖33顯示即使進口處的畸變僅僅是總壓畸變，但由于加功量的周向分布因總壓畸變而產(chǎn)生不均勻性，出口處的總溫分布也會呈現(xiàn)明顯的不均勻性。風(fēng)扇轉(zhuǎn)靜子在進口總壓畸變條件下流場數(shù)值模擬計算為驗證程序計算多排葉片的能力，在風(fēng)扇轉(zhuǎn)子葉片后加上一排靜子葉片進行兩排葉片的總壓畸變場模擬計算，計算選取的相對換算轉(zhuǎn)速為=1.0，進口畸變場與5.4.1.1節(jié)中狀態(tài)2相同。圖34為轉(zhuǎn)靜子聯(lián)算的子午網(wǎng)格圖。圖35為不同軸向位置處總壓沿周向的分布，可以清楚地看到總壓畸變沿軸向的衰減，總壓畸變在經(jīng)過轉(zhuǎn)子后被衰減，經(jīng)過靜子葉片后又進一步衰減，說明靜子葉片也具有衰減總壓畸變的能力。圖36為風(fēng)扇不同軸向位置處總壓云圖，由圖可見畸變流的三維性非常明顯，風(fēng)扇上游的周向畸變在經(jīng)過風(fēng)扇轉(zhuǎn)子后，葉中的低總壓區(qū)已經(jīng)幾乎被抹平了，但葉尖和葉根的低壓區(qū)還很明顯，而且有一定的相位移動。氣流流過靜子后葉尖區(qū)的低壓區(qū)也被抹平了，低壓區(qū)只存在在葉根區(qū)附近。圖34風(fēng)扇轉(zhuǎn)靜子聯(lián)算子午網(wǎng)格圖圖35風(fēng)扇不同軸向位置處換算總壓沿周向分布(a)風(fēng)扇上游(b)轉(zhuǎn)子出口(c)靜子出口圖36風(fēng)扇不同軸向位置處總壓云圖5.4.1.5穩(wěn)定性預(yù)估工作是目前流體力學(xué)界公認的難點問題，在相關(guān)的文獻[7][15]中可以看到目前判斷流場穩(wěn)定性的一種方法是在流場中加入一個波長為若干個動葉柵距的擾動，然后根據(jù)擾動的發(fā)展來判斷失速發(fā)生，目前這種方法還在研究階段，只能夠定性地判斷失速。作者試圖利用這種方法來判斷失速，但利用目前的模型幾乎算不出失速現(xiàn)象。分析認為如果要利用這種方法來預(yù)測失速，目前的模型中還需進一步改進，因為在近失速點損失大幅度提高，葉片增壓能力下降，其特性線在小于某一流量時壓比會出現(xiàn)拐頭，模型要能夠算出這種現(xiàn)象才能夠利用加擾動的方法來判斷穩(wěn)定性。由于時間關(guān)系模型沒有加入模擬失速的部分，所以無法用加擾動的辦法來判斷失速。作者提出用更簡單和常用的辦法來判斷失速，即利用葉柵的擴散因子是否超標(biāo)來判斷失速。擴散因子是表征葉柵逆壓力梯度的參數(shù)，葉柵只有在逆壓力梯度足夠大時才能發(fā)生分離失速，所以用擴散因子來判斷失速是有一定的物理依據(jù)的。程序中采用的擴散因子定義式為：其中W1、W2為轉(zhuǎn)子葉片進、出口相對速度，為扭速，為稠度。首先計算了均勻進氣條件下設(shè)計轉(zhuǎn)速近喘點343kg/s流量時沿風(fēng)扇葉片擴散因子最大值為0.63。然后計算了相同流量下畸變指數(shù)為4.6%的總壓畸變進口時風(fēng)扇轉(zhuǎn)子整個葉高最大擴散因子沿周向的分布，計算結(jié)果見圖48。相同流量下進口總壓畸變流情況下在畸變區(qū)擴散因子最大超過了0.7，如果假設(shè)風(fēng)扇可承受的最大擴散因子是0.63的話，那么在進口有總壓畸變的情況下，風(fēng)扇就不可能在相同的流量下還保持穩(wěn)定，它的穩(wěn)定流量比均勻進口情況下要小。顯然在進口總壓畸變時風(fēng)扇的壓比降低、流量減小，而且穩(wěn)定邊界點的流量也比均勻情況下要小，因此其穩(wěn)定裕度較均勻進口情況將有所降低。圖37進口總壓畸變條件下風(fēng)扇轉(zhuǎn)子擴散因子沿周向分布因為風(fēng)扇試驗時沒有試到喘點，本文無法驗證利用擴散因子的方法是否能夠定量地判斷失速點。但作者認為在周向平均體積力模型加三維歐拉方程的計算模型中利用擴散因子來判斷失速還是可行的，因為葉柵通道內(nèi)的周向平均體積力模型與S2程序十分類似，它和S2程序一樣不考慮葉片通道內(nèi)部周向的不均勻，所不同的是S2程序是沿壓氣機整個環(huán)形通道的軸對稱，而本文程序是沿單個葉片通道的軸對稱。在S2程序中擴散因子是衡量葉柵逆壓力梯度最有效的參數(shù)，其有效性經(jīng)多年的研究和使用所驗證。穩(wěn)態(tài)總溫畸變運輸機使用反推力裝置，戰(zhàn)斗機吸入其他機群排出的氣體，艦載機吸入蒸氣彈射裝置產(chǎn)生的熱蒸氣，直升機接近地面懸停時燃氣的吸入等都引起溫度畸變問題，此外低壓壓氣機進口的總壓畸變往往衍生高壓壓氣機的總溫畸變。為驗證模型模擬穩(wěn)態(tài)總溫畸變能力，在風(fēng)扇轉(zhuǎn)子上游加入一周向穩(wěn)態(tài)溫度畸變進行模擬計算，畸變區(qū)扇形角121度，畸變區(qū)溫升87.41k,畸變指數(shù)18.2%。溫度畸變指數(shù)定義公式如下：式中：——畸變指數(shù)——高溫區(qū)總溫平均值——進口總溫平均值風(fēng)扇上游畸變場總溫云圖和總溫沿周向分布見圖38。總溫云圖總溫沿周向分布圖38穩(wěn)態(tài)總溫畸變模擬計算時風(fēng)扇上游總溫畸變場(a)風(fēng)扇上游(b)風(fēng)扇進口(c)風(fēng)扇出口圖39風(fēng)扇不同軸向位置處的總溫云圖圖39列出了風(fēng)扇不同軸向位置處總溫云圖，可以看出溫度畸變沿流程的變化，在風(fēng)扇上游標(biāo)準方波形的周向畸變在風(fēng)扇進口處高溫區(qū)沿周向向低溫區(qū)擴散，風(fēng)扇出口截面周向總溫畸變已經(jīng)演變成了周向和徑向的混合畸變。圖40為風(fēng)扇上游和出口總溫對比，可以看出通過風(fēng)扇后總溫不均勻性有少量的衰減，畸變區(qū)沿周向略有偏移?；儏^(qū)的高總溫被衰減可以這樣解釋：畸變區(qū)較高的總溫與非畸變區(qū)相比因物理轉(zhuǎn)速相同而對應(yīng)著較低的換算轉(zhuǎn)速，因此加功量小，溫升小，緩解了畸變區(qū)總溫高的趨勢。圖41為風(fēng)扇不同軸向位置處總壓云圖，圖42為風(fēng)扇上游和出口總壓對比，由這兩張圖可以看到畸變區(qū)的總壓升確實比非畸變區(qū)小，所以雖然進口處的畸變僅僅是總溫畸變，出口處的總壓分布也呈現(xiàn)明顯的不均勻性，解釋了風(fēng)扇進口溫度畸變誘導(dǎo)高壓進口壓力畸變的現(xiàn)象。圖40風(fēng)扇上游和風(fēng)扇出口總溫對比(a)風(fēng)扇上游(b)風(fēng)扇進口(c)風(fēng)扇出口圖41不同軸向位置總壓云圖圖42風(fēng)扇上游和風(fēng)扇出口總壓周向分布圖43風(fēng)扇上游和風(fēng)扇進口軸向速度對比圖43顯示了從風(fēng)扇上游位置到風(fēng)扇進口軸向速度的變化，低總溫區(qū)的流動被加速，而高總溫區(qū)基本沒有。平行壓氣機理論可以很好的解釋這一現(xiàn)象：具有高總溫來流的壓氣機部分必須在低質(zhì)量流量下工作，而其它部分則須在高質(zhì)量流量下工作，畸變區(qū)的增壓能力低于非畸變區(qū)，而出口靜壓幾乎是均勻的，因此流動畸變區(qū)域在進入風(fēng)扇之前必須相對非畸變區(qū)減速增壓。模擬計算與Gong的計算結(jié)果十分相似，而且與平行壓氣機理論定性的分析一致，可以認為本模型能夠很好地反映穩(wěn)態(tài)溫度畸變與壓氣機的相互作用。5.5動態(tài)進氣畸變條件下模擬計算本模型的一大特點就是能夠模擬壓氣機對動態(tài)畸變的響應(yīng)，為了驗證模型這一功能我們將模擬風(fēng)扇轉(zhuǎn)子在進口總壓周期振蕩和進口溫度快速瞬變條件下的畸變流場。5.5.1進口周期壓力振蕩動態(tài)總壓畸變起源于進氣道的外部和內(nèi)部兩個方面，地面旋風(fēng)、陣風(fēng)的吸入屬進氣道外部起因，進氣道工作不穩(wěn)定、嗡鳴、以及激波與附面層干擾等則屬內(nèi)部起因，它們可以產(chǎn)生周期性壓力震蕩和脈沖；戰(zhàn)斗機打機關(guān)炮時也會產(chǎn)生周期壓力振蕩。文獻[11]中對周期性振蕩型壓力畸變的研究認為，葉輪機經(jīng)受這類動態(tài)畸變后，各級工作不再匹配。引起級的工作失配的主要機理是壓力脈動的幅值衰減和相位錯移，在各個頻率下的幅值衰減和相位差關(guān)系就是動態(tài)系統(tǒng)的幅頻特性和相頻特性，其變化特點是頻率越高幅值比越小、相位差越大。本文模擬計算了三種不同頻率、相同幅值的周期壓力振蕩下風(fēng)扇的響應(yīng)，來了解該風(fēng)扇對周期壓力振蕩的頻率特性，并檢驗程序計算非定常流動的能力。模擬計算設(shè)定總壓振蕩畸變區(qū)扇形角98度，畸變區(qū)總壓振蕩幅值12650Pa,最大畸變指數(shù)4.621%，壓力振蕩頻率分別為10Hz，20Hz，40Hz，計算時在風(fēng)扇上游邊界強加一個周期振蕩的壓力波，圖44為上游邊界處某一時刻的壓力云圖和三種振蕩頻率下上游邊界畸變區(qū)總壓隨時間的變化。圖44模擬計算時風(fēng)扇上游邊界總壓云圖及其畸變區(qū)總壓隨時間的變化通過對計算結(jié)果進行分析得出如下結(jié)論：風(fēng)扇對這種動態(tài)的壓力振蕩表現(xiàn)出了良好的幅頻特性，見圖45，隨著頻率的增大，風(fēng)扇出口壓力振動幅值與進口之比即幅值比減小，這個結(jié)果與文獻[11]的實驗和計算結(jié)果是一致的。這說明在壓力振蕩幅值相等的情況下，低頻的壓力振蕩經(jīng)過壓氣機后衰減的幅度要小，對壓氣機性能的影響更大。模擬計算的相頻特性不明顯，風(fēng)扇進出口壓力脈動的相位差很小，頻率的影響也不大，分析其原因認為，主要是因為計算的葉排數(shù)太少，一排葉片對壓力振蕩的影響力不足以使振蕩波的相位發(fā)生明顯的變化。圖46在風(fēng)扇特性圖上繪制出不同壓力振蕩頻率下風(fēng)扇工作點運動軌跡，由圖可見低頻壓力振蕩時工作點偏離均勻進氣工作點最多，流量的變化范圍也最大，低進口總壓時刻風(fēng)扇的工作點越接近穩(wěn)定工作邊界，而不同的頻率下的壓比變化范圍相差不是很大。和上面的結(jié)論一致，在壓力振蕩幅值相同的情況下，頻率越低對風(fēng)扇的影響越大，越偏危險。數(shù)值模擬的結(jié)果可見，模型很好地反映出風(fēng)扇對周期壓力振蕩的幅頻特性，證明了程序具有模擬動態(tài)畸變的能力。圖45動態(tài)壓力振蕩下風(fēng)扇的幅頻特性圖46不同壓力振蕩頻率下風(fēng)扇工作點運動軌跡5.5.2進口溫度快速瞬變飛機在發(fā)射武器后發(fā)動機吸入燃氣，使壓氣機進口溫度隨時間突升。一架垂直起落飛機，當(dāng)噴口轉(zhuǎn)到垂直向下排氣時，通過地面反射，發(fā)動機進口的下半部分吸入燃氣，也是一個溫度瞬變的過程。它們的變化近似于一個三角形波，一般用溫升率表示瞬變特征。對渦輪噴氣發(fā)動機的試驗結(jié)果表明[11]，壓氣機對溫度瞬變的穩(wěn)定性響應(yīng)主要取決于溫升率(△T/△t)。當(dāng)壓氣機轉(zhuǎn)速一定時，進口溫度突升，就使換算轉(zhuǎn)速降低，而壓比變化不大，在特性圖上壓氣機工作點會沿著等壓比線向左移動，接近穩(wěn)定工作邊界。圖47模擬計算時風(fēng)扇上游邊界總溫云圖及其畸變區(qū)總溫隨時間的變化模擬計算分別在壓氣機進口加入相同溫升不同溫升率的2組溫度脈沖，目的是對比不同溫升率對壓氣機的影響。根據(jù)目前掌握資料武器發(fā)射時溫升率基本上在500-2500k/s的范圍內(nèi)[12]，因此計算的溫升率在這個范圍內(nèi)選?。?000k/s、2000k/s，溫升為26.7k，高溫區(qū)仍在周向108度的范圍內(nèi)。進口總溫分布云圖和畸變區(qū)總溫隨時間的變化見圖47。圖48為在這種動態(tài)總溫情況下風(fēng)扇流量隨時間的變化，可以看到隨著溫度的上升流量減小，但溫度上升到最大值時流量并未下到最小值，而是有一個滯后，如果假設(shè)溫升率為1000k/s的溫度脈沖周期為T1，溫升率為2000k/s的溫度脈沖周期為T2，溫升率為1000k/s時溫度最大值與流量最小值之間約相差0.17T1，而溫升率為2000k/s情況下兩者之間相差約0.29T2，說明溫升率越高，流量相應(yīng)的滯后越嚴重。同時圖45中可以看到，雖然兩個溫度脈沖的溫升是相同的，但低溫升率的溫度脈沖引起的流量下降更嚴重些，這一點很好地驗證了相同溫升情況下低溫升率更加危險的結(jié)論。圖49為風(fēng)扇壓比隨時間的變化，總體上溫度脈沖對風(fēng)扇壓比的影響不大，風(fēng)扇壓比略有降低，最低壓比點也相對于最高溫度點有一定滯后，其規(guī)律與流量變化相同，所不同的是不同的溫升率情況下壓比的降低幾乎相同。圖50在風(fēng)扇特性圖上描繪了發(fā)生溫度脈沖過程中風(fēng)扇工作點的運動軌跡，風(fēng)扇進口局部區(qū)域溫度的動態(tài)升高使風(fēng)扇換算轉(zhuǎn)速降低，工作點沿著等壓比線向左側(cè)移動，使喘振裕度降低，極限情況是有可能跨過喘振邊界發(fā)生喘振。同時溫升率越低工作點移動得越遠，越危險，這些現(xiàn)象與以往試驗和理論研究的結(jié)論都是一致的。