• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    特征值理論在穩(wěn)定性預測中的應用研究進展綜述

    2022-11-05 04:17:50孫曉峰董旭張光宇王卓孫大坤
    航空學報 2022年10期
    關(guān)鍵詞:壓氣機特征值擾動

    孫曉峰,董旭,張光宇,王卓,孫大坤,*

    1. 北京航空航天大學 能源與動力工程學院,北京 100191 2. 北京航空航天大學 航空發(fā)動機研究院,北京 100191

    長期以來,中國在航空發(fā)動機領域的發(fā)展始終徘徊在舉步維艱的探索階段,在發(fā)動機設計、制造、試驗以及材料工藝方面同歐美國家存在著一定的差距,這導致了目前以殲20、運20和C919為代表的軍用民用先進飛機遲遲無法裝備“中國心”。在“兩機”專項的重大部署和國家科技發(fā)展前瞻性布局的支持引導下,與航空發(fā)動機相關(guān)的諸多工程和基礎研究領域的問題已經(jīng)得到了解決,而進一步深入發(fā)掘工程問題背后的基礎理論與方法,積極儲備面向未來的先進技術(shù),對于中國現(xiàn)有型號的研制和未來型號預研都具有重要意義。

    現(xiàn)代先進航空發(fā)動機的發(fā)展方向不僅是追求更大的推力和更高的效率,大涵道比民用航空發(fā)動機也在同時追求更低的排放和噪聲,而這些追求的背后也蘊藏了諸多工程和基礎研究領域的關(guān)鍵科學問題,不僅包括如何利用氣動熱力學基本原理實現(xiàn)高性能設計、如何通過模擬燃燒與化學反應動力學過程提高燃燒效率和減少污染物排放、以及如何利用氣動聲學相關(guān)理論降低風扇和噴流噪聲等,更是引發(fā)了各種亟待解決的復雜流固熱聲耦合問題,尤其是高推重比與高效率氣動設計帶來了葉輪機械負荷的增加與流動復雜性的增強,這勢必會造成壓氣機穩(wěn)定性裕度不足、風扇噪聲增強、葉片流致/渦致振動和燃燒室熱聲振蕩等一系列危及發(fā)動機穩(wěn)定工作和飛機飛行安全的工程問題。

    作為航空發(fā)動機最主要的3大部件,風扇/壓氣機、燃燒室和渦輪,其各自均對發(fā)動機整機性能和可靠性起著決定性的影響。對于壓縮部件,高壓比、高性能的氣動設計的工程應用主要受到穩(wěn)定性問題的制約。通常,風扇/壓氣機的工作圖譜存在著左側(cè)失穩(wěn)邊界,而發(fā)動機穩(wěn)定工作在共同工作線上,為了避免由于過渡態(tài)以及進氣畸變等造成工作點偏移越過失穩(wěn)邊界遭遇失速和喘振,工作線與失穩(wěn)邊界需要保持一定的距離,也就是所謂的穩(wěn)定性裕度。半個多世紀的研究與實踐均著眼于如何有效地將失穩(wěn)邊界左移,并在航空發(fā)動機壓縮系統(tǒng)所需的工作范圍內(nèi)避免失穩(wěn)現(xiàn)象(失速和喘振)的發(fā)生,研究不僅涵蓋了其非定常演化的特點,同時也發(fā)展了一系列的主/被動穩(wěn)定性控制方法。早期的關(guān)于葉輪機械內(nèi)部流動失穩(wěn)問題的各類解析模型也大都關(guān)注于失速喘振的起始階段,從而建立各種各樣的特征值問題。然而,為了得到流場擾動的解析解,這些模型中需要做大量的簡化和假設。模型的預測精度取決于風扇/壓氣機損失、落后角和壓比等特性曲線的準確性,而這些特性曲線常常又不具有足夠的精度,特別是在新的壓氣機型號設計階段。大量對經(jīng)驗關(guān)系的需求限制了這些模型的應用價值,特別是在實際工程中的應用。因此,發(fā)展能夠考慮復雜流場及幾何條件,不嚴重依賴經(jīng)驗關(guān)系,適用于葉輪機械設計階段,并且對流動失穩(wěn)邊界可以較為快速準確地給出清晰判據(jù)的穩(wěn)定性模型勢在必行,非常富有理論和工程指導意義。

    不過即便在設計階段保留了足夠的氣動穩(wěn)定裕度并采取了穩(wěn)定性控制措施,真實工況下運轉(zhuǎn)的壓氣機和風扇在部分中低轉(zhuǎn)速時會在近失速區(qū)發(fā)生葉片顫振,英國帝國理工大學Vahdati教授團隊與Rolls-Royce開展了長期的合作,通過復雜的CFD計算探究了這種氣彈失穩(wěn)現(xiàn)象的特性,提出了失速邊界在部分中低轉(zhuǎn)速氣彈失穩(wěn)影響下存在“Flutter Bite”現(xiàn)象[1],二者共同決定了風扇的穩(wěn)定邊界。氣動彈性穩(wěn)定性問題涉及到復雜的氣體和固體相互作用,早期學者的研究大多基于解耦假設,因此將氣動彈性穩(wěn)定性問題重點放在求解葉片表面氣動力上,進而利用特征值法及能量法對葉片運動穩(wěn)定性進行判定。由于受到計算能力的限制,先后發(fā)展了線性氣動力模型、全勢流模型、跨聲速小擾動模型以及升力面模型等用于葉片表面氣動力預測,從而利用特征值方法或者能量法對葉片氣動彈性穩(wěn)定性進行判斷。直到20世紀80年代開始利用CFD方法進行流固耦合求解,先后求解Euler方程、RANS方程,但大多數(shù)均是弱耦合方法,在此階段為了降低計算難度,Hall等[2-3]發(fā)展了處理葉輪機氣動彈性穩(wěn)定性降階模型 (Reduced Order Models),使 CFD 在氣動彈性穩(wěn)定性領域開始占據(jù)極其重要的地位。

    隨著商用航空發(fā)動機追求更低污染物排放,軍用航空發(fā)動機朝向更高推重比目標發(fā)展,民用航空發(fā)動機環(huán)形燃燒室和軍用航空發(fā)動機加力燃燒室均面臨著愈發(fā)復雜凸顯的燃燒不穩(wěn)定性問題。燃燒不穩(wěn)定性問題往往發(fā)生在高溫高壓環(huán)境下,涉及燃燒化學反應、湍流流動及復雜環(huán)境聲學等多時空尺度問題,通過實驗和高精度數(shù)值方法能夠清晰地捕捉局部火焰反應、火焰和流動干涉、火焰對聲波響應等局部關(guān)鍵機理信息,但是高成本的實驗、高難度的測試、高成本的計算資源決定了不能完全依賴這些方法用于重復性計算的參數(shù)優(yōu)化設計來支撐燃燒室的研制。因此,通過模型方法對燃燒不穩(wěn)定性發(fā)生頻率、模態(tài)進行準確預測,并發(fā)展可靠的控制手段成為了重要研究方向。

    諸多物理和工程問題在數(shù)學上都歸結(jié)為求解矩陣的特征值和特征向量的問題,這主要是由于所研究的對象往往都是以動力學微分方程組作為控制方程來描述的,并且由于非線性模型通??梢员唤茷榫€性模型,因此,求解代數(shù)微分方程的解析解和數(shù)值解就成為了自然科學乃至社會科學問題模型研究的必要手段,相應的研究對象也逐步由代數(shù)方程轉(zhuǎn)移到矩陣本身,依靠矩陣相似變換來求矩陣的特征值與特征向量,并以此來反映線性系統(tǒng)的某些物理特征。特征值理論在求解多種穩(wěn)定性問題上得到了應用,人們在求解諸如Poiseuille流動和Couette流動的穩(wěn)定性和轉(zhuǎn)捩問題時采用的方法是線性全局穩(wěn)定性分析方法[4],假定我們可以精確已知或者估計得到Navier-Stokes方程的定常流動解,如果向這個流動解中引入小擾動,那么線化小擾動方程以及相關(guān)的邊界條件就是齊次的且與exp(σt)相關(guān),當擾動方程的解存在且Re(σ)>0時,流場就是線化不穩(wěn)定的。以此為基礎,通過小擾動假設和特征值理論,可以研究多種葉輪機領域常見的穩(wěn)定性問題,因此,本文將簡要綜述特征值理論在航空發(fā)動機流動、氣彈和燃燒穩(wěn)定性預測方面的應用研究現(xiàn)狀,并給出筆者所在課題組在相關(guān)理論建模的分析與效果驗證。

    1 基于特征值理論的壓氣機流動穩(wěn)定性預測

    1.1 線化小擾動穩(wěn)定性模型

    自然界中存在多種流動不穩(wěn)定性問題,比如與重力和加熱相關(guān)的穩(wěn)定性問題、分層流穩(wěn)定性問題、與離心慣性力有關(guān)的穩(wěn)定性問題、與表面張力有關(guān)的穩(wěn)定性問題、湍流轉(zhuǎn)捩問題。各類穩(wěn)定性問題由不同的物理機制所支配。研究流動穩(wěn)定性問題的方法主要分為線性穩(wěn)定性理論和非線性穩(wěn)定性理論兩種。線性穩(wěn)定性理論是從無限小擾動的假設出發(fā),流場中的各種變量被分解為平均量和擾動量兩部分,并且忽略了小擾動的高階項,即不考慮擾動之間的相互作用,從而使方程線性化。以線性化方程為基礎建立的理論即為線性穩(wěn)定性理論。線性穩(wěn)定性理論主要關(guān)注的是流體失穩(wěn)的起始階段,而不能用來描述失穩(wěn)后轉(zhuǎn)捩的全過程,因為這個過程主要是非線性效應起決定性作用的階段。但線性理論可以說明哪種速度剖面在擾動作用下會產(chǎn)生失穩(wěn),哪些頻率的振動增長最快等,這些對我們預測失穩(wěn)的發(fā)生有著重要意義。而非線性穩(wěn)定性理論是研究有限振幅擾動波的情況,這時擾動速度與主流平均速度比較起來就不再是很小的量,因此在線化過程中省略掉的擾動量的二階項不能再被忽略,小擾動方法所計算的結(jié)果逐漸偏離實際擾動發(fā)展的過程,非線性效應開始顯現(xiàn)。按照線性穩(wěn)定性理論,一旦失穩(wěn)開始后,擾動應該按照指數(shù)關(guān)系迅速增長,而實際情況決不可能如此。

    實際上,對于湍流轉(zhuǎn)捩問題而言,流動失穩(wěn)往往是轉(zhuǎn)捩過程的開始,而并非其全部過程,穩(wěn)定性理論也不能描述其全過程。然而,目前沒有任何一種完整的理論分析能解析地描述轉(zhuǎn)捩的全過程。雖然線性穩(wěn)定性理論不能用于湍流轉(zhuǎn)捩過程的非線性階段,但卻可以說明哪種速度剖面最不穩(wěn)定,哪些頻率波動增長最快。這和葉輪機械內(nèi)部的流動失穩(wěn)問題十分類似,又略有不同。

    葉輪機械內(nèi)部流動的不穩(wěn)定現(xiàn)象屬于流動不穩(wěn)定性問題中的一種。同時,作為旋轉(zhuǎn)機械,其變工況特征、復雜的幾何形狀和固壁邊界使這一問題更加棘手。迄今為止,葉輪機流動穩(wěn)定性理論預測研究方法大致可以分為兩類:數(shù)值計算方法[7-15]和半解析半數(shù)值的理論模型方法[16-27]。以Emmons模型[16]和Stenning模型[17]為代表的線化小擾動理論模型、以Greitzer模型[21-22]和Moore-Greitzer (M-G模型[23-24])為代表的系統(tǒng)穩(wěn)定性模型均是典型的半解析半數(shù)值理論模型,其優(yōu)點在于方法判據(jù)清晰、求解快速,適合處理多級壓縮系統(tǒng)失速起始判定問題,線化小擾動穩(wěn)定性模型基于的就是二維不可壓小擾動方法,這里不妨展開簡要介紹:1955年Emmons[16]采用二維不可壓小擾動方法,給出了葉片排出口氣流堵塞面積隨攻角的變化,建立了壓氣機失速預測模型,并給出了旋轉(zhuǎn)失速物理現(xiàn)象的經(jīng)典描述。Nenni和Ludwig[19]建立了能夠考慮單排和雙排葉柵影響的二維不可壓旋轉(zhuǎn)失速穩(wěn)定性模型,同樣給出了穩(wěn)定性判定的解析表達式,并指出總壓損失系數(shù)及其導數(shù)是影響壓氣機穩(wěn)定性的最重要參數(shù)。

    而系統(tǒng)穩(wěn)定性模型基于彈簧、質(zhì)量、阻尼系統(tǒng)類比,Greitzer建立了最初的一維系統(tǒng)穩(wěn)定性模型,后來Moore和Greitzer將該模型(M-G模型[23-24])推廣到二維不可壓流,能夠進行喘振、失速的穩(wěn)定性分析。M-G模型采用了大量簡化處理,將整個壓縮系統(tǒng)作為激盤,無法包含壓氣機復雜幾何造型等因素對穩(wěn)定性的影響效果,難以適用于具有明顯三維特征的旋轉(zhuǎn)失速問題,并且其假設壓氣機都在靜壓升系數(shù)斜率為零的條件下失穩(wěn),用B參數(shù)作為判斷失穩(wěn)模式的關(guān)鍵參數(shù),這不是本文重點討論的內(nèi)容,因此不再贅述。

    采用小擾動方法分析葉輪機流動穩(wěn)定性問題在很大程度上完全符合研究對象的物理特征,圖2給出某單級風扇失速瞬態(tài)壁面靜壓信號空間傅里葉變換幅值曲線[28],可以看出:擾動波在初始階段的幅值很小,隨后開始逐漸增長直至某一時刻突然非線性放大進入失穩(wěn)狀態(tài)。因此,若需要研究失速演化的過程,或者對于某些大尺度擾動(如進口流場畸變引起的流動不穩(wěn)定性)對流場的影響,可以采用非線性穩(wěn)定性理論或者數(shù)值模擬的方式。然而,如果我們只關(guān)注失速起始階段的系統(tǒng)穩(wěn)定特性,完全可以采用線性小擾動理論來預測失速的產(chǎn)生。

    壓氣機的流動失穩(wěn)過程尤其是在失穩(wěn)起始階段,通常被認為是擾動演化的開始,此時的擾動幅值較小滿足線性假設,通過線性穩(wěn)定性分析的方式對小擾動的演化趨勢進行分析判別可以在一定程度上判斷流動的穩(wěn)定趨勢。在穩(wěn)定性模型發(fā)展的早期階段,經(jīng)典的線化模型在揭示旋轉(zhuǎn)失速發(fā)生機理層面具有重要的意義,但是大量的諸如不可壓、無黏、無旋和均勻背景流動等簡化假設決定了這些將葉片簡化為激盤或者半激盤的集中參數(shù)模型無法包含流動非均勻性和葉片力的影響。

    為了克服早期經(jīng)典穩(wěn)定性模型的這一缺陷,人們開始利用數(shù)值線性穩(wěn)定性分析方法開展失穩(wěn)預測的嘗試,美國麻省理工學院Gordon提出了一個三維不可壓數(shù)值穩(wěn)定性模型[25],在基本流軸對稱的假設基礎上,用體積力模型替代葉片對流體的作用,包含了子午流面內(nèi)流動非均勻性。不過,該模型采用的葉片力模型較為粗糙,并且不能包含壓縮性的影響,其預測精度和使用范圍受到了限制。然而,針對現(xiàn)代高負荷高速壓氣機而言,發(fā)展三維可壓縮穩(wěn)定性模型勢在必行。1996年,孫曉峰[26]從三維可壓縮線化非定常Euler方程出發(fā),通過與研究流動穩(wěn)定性所用的三維O-S方程的相應工作類比,采用無窮維系統(tǒng)的截斷技術(shù)以及有限空間波傳播界面間的模態(tài)匹配方法,發(fā)展了可以考慮任意階徑向擾動的三維可壓縮旋轉(zhuǎn)失速穩(wěn)定性模型。在此工作基礎上,本課題組[27-28]結(jié)合等價分布源法和傳遞單元法,將激波匹配關(guān)系成功引入半激盤模型中,發(fā)展了跨聲速風扇/壓氣機穩(wěn)定性預測與擴穩(wěn)理論模型。該模型不關(guān)心擾動量的具體數(shù)值,因為擾動量的大小是隨時間變化的。而失穩(wěn)是因為擾動量隨時間放大所導致,因此該模型只根據(jù)擾動量隨時間放大還是衰減來判斷系統(tǒng)穩(wěn)定與否。為了更精確地反應流場內(nèi)的物理畫面,使之包含更多的流動細節(jié)(尤其是小尺度的擾動),葉片造型和黏性的影響就要在計算模型中得到考慮,因此求解非定常雷諾平均Navier-Stokes方程的方法開始被引入到壓氣機穩(wěn)定性模型中來。

    一個動力學系統(tǒng)的穩(wěn)定性取決于其自身對于系統(tǒng)內(nèi)部或外部的任何小擾動的響應,而并非針對某一特定的擾動。對于大部分流動穩(wěn)定性問題而言,擾動總是經(jīng)歷由小到大的發(fā)展過程。如果我們更關(guān)注流動失穩(wěn)的開始階段,那么采用線化處理是可行且合理的。這也就是所謂的理論上任何流動失穩(wěn)開始問題都可以描述為一個基于Navier-Stokes方程的特征值問題。

    該預測方法的有效性在Stage37和NASA兩級風扇上得到了驗證[28],Stage37在設計轉(zhuǎn)速下葉尖馬赫數(shù)為1.4,遠遠超出了亞聲速模型的假定范圍,因此,亞模型不適合該轉(zhuǎn)速。即使90%的設計轉(zhuǎn)速下,葉尖部分的大部分流動也是跨聲速,只有在70%的設計轉(zhuǎn)速下,葉尖的進口相對馬赫數(shù)才小于1,因此可以試用文獻[27]所發(fā)展的模型來預測其失速起始點。

    模型預測的結(jié)果見圖3,圖中有兩簇曲線,一簇對應延遲時間為0的情況,一簇對應取正常延遲時間的情況;延遲時間在高亞聲速壓氣機中的作用相比低速壓氣機要明顯的多,從圖3中可以看出,加上延遲時間后,中性穩(wěn)定點要略向高流量點偏移,高階模態(tài)下的偏移要大一些;而延遲時間對相對周向傳播速度的影響似乎要更明顯一些,由此可見,該亞聲速穩(wěn)定性模型在高亞聲速條件下也能夠做出合乎實際的預測。

    NASA兩級風扇在第1級轉(zhuǎn)子采用小展弦比設計,目的是消除原有大展弦比風扇葉片的弊端,改善通道流場,提高級壓比,減少葉片數(shù)目,并用于研究小展弦比葉片的優(yōu)越性。此外其相對充分的試驗數(shù)據(jù),也可以用來評估模型對多級壓氣機穩(wěn)定性失速先兆點的預測能力。其中,NASA兩級風扇失速實驗通過在不同的展向位置同時采集徑向11處不同點的數(shù)據(jù)來研究設計轉(zhuǎn)速從50%~100%的失速邊界。在完全設計轉(zhuǎn)速情況下,質(zhì)量流量為33.248 kg/s,總壓力比為2.399,絕熱效率為0.849(圖4)。

    圖4的結(jié)果顯示了模型對NASA兩級風扇100%設計轉(zhuǎn)速下的穩(wěn)定性預測。當衰減因子趨于0,壓縮機系統(tǒng)的穩(wěn)定性由穩(wěn)定狀態(tài)變?yōu)椴环€(wěn)定狀態(tài),從而得到失速點,通過分析可以看出理論預測和實驗結(jié)果之間的誤差約5%。

    1.2 穩(wěn)定性預測通用理論

    考慮到解析模型和數(shù)值模型各自的特點,近年來,一種新的建立穩(wěn)定性模型的方式得到了發(fā)展[29-30],一方面像解析模型那樣,能夠為判斷失速邊界或者失穩(wěn)點提供清晰的判斷依據(jù)或準則,同時像數(shù)值模型那樣,能夠包含葉片造型對流場結(jié)構(gòu)的影響,并且不像解析模型那樣嚴重依賴經(jīng)驗關(guān)系。這樣一種方式,可以提高模型預測的準確性,同時在設計階段為設計人員提供一種更好的理論預測工具,考慮葉片彎掠、流道結(jié)構(gòu)等復雜造型參數(shù)對于穩(wěn)定性的定量影響,這是以往任何模型所難以比擬的??梢栽O想,通過應用這樣一種在數(shù)值解法上可以快速實現(xiàn)、物理層面上可以包含更多因素、失穩(wěn)點判斷清晰明確的模型方法,使其作為穩(wěn)定性設計而成為壓氣機設計階段中的關(guān)鍵一環(huán),而這樣一環(huán)是當前葉輪機械設計體系中所欠缺的。為了更好地解決風扇/壓氣機流動穩(wěn)定性帶來的一系列問題,有必要將穩(wěn)定性設計納入到現(xiàn)代風扇/壓氣機設計體系中(圖5[29-30])。葉輪機內(nèi)流通用穩(wěn)定性模型就是針對這一體系所做的探索實踐。

    類比層流到滯流的轉(zhuǎn)捩過程問題的處理方法[31-32],將壓氣機流場的定常數(shù)值計算結(jié)果作為某一流動狀態(tài)的平均流場,并在此基礎上引入一個包含所有物理量的小擾動,通過計算將壓氣機流動作為動力學系統(tǒng)的特征頻率來判斷其流動穩(wěn)定性。對于任意的離心或是軸流壓氣機系統(tǒng),從最基本的Navier-Stokes方程出發(fā),應用線化小擾動理論,結(jié)合浸入式邊界方法和體積力方程,發(fā)展通用的流動穩(wěn)定性模型,該模型可以用于任意進出口條件和邊界情況下的復雜流動失穩(wěn)起始點的預測。

    該方法的具體操作流程為:從控制方程出發(fā)——CFD主流流場——體積力模型——建立特征值問題——求解特征值問題。

    根據(jù)浸入式邊界理論,葉輪機葉片用分布的力源項來代替,從而,葉輪機內(nèi)部流動由帶力源項的Navier-Stokes方程來控制,即:

    質(zhì)量守恒方程:

    (1)

    式中:ρ為密度;u為速度矢量。

    動量守恒方程:

    (2)

    式中:F為單位質(zhì)量流量的葉片力矢量;Π為二階應力張量。

    能量守恒方程:

    (3)

    式中:λ為導熱系數(shù);WF表示葉片力做功。

    如1.1節(jié)所述,發(fā)展葉輪機流動通用穩(wěn)定性理論的核心目的是進行失速起始點的預測,因此,基于小擾動假設,把原始三維非定常流動分解成定?;玖骱头嵌ǔP_動之和:

    (4)

    根據(jù)多元函數(shù)泰勒展開,葉片力可以線化為

    (5)

    將式(4)、式(5)代入式(1)~式(3),減掉基本流所滿足的零階方程并且忽略高階無窮小量,則可得到線化形式的Navier-Stokes方程。將線化Navier-Stokes方程在圓柱坐標系下展開,可得:

    (6)

    式中:Φ=[ρ′u′v′w′p′]T,而A、B、C、E、G、H、M、N、Q、R以及S是由基本流所確定的系數(shù)矩陣。

    一般地,葉輪機內(nèi)部基本流在徑向、周向以及軸向都是非均勻的,假定基本流是定常的,因此,只能對小擾動量作關(guān)于時間的正則展開,即:

    (7)

    將式(7)代入到式(6),可以得到以復頻率ω為特征值的特征值方程:

    (8)

    根據(jù)齊次線性偏微分方程理論,式(8)要有非零解,則必須滿足:

    det[L(ω)]=0

    (9)

    通過求解式(9)即可得到復頻率ω=ωr+iωi,復頻率的實部ωr表征失速先兆波的頻率,而復頻率的虛部ωi則表征小擾動的演化趨勢:ωi為正,則擾動隨時間放大,系統(tǒng)不穩(wěn)定;ωi為負,則擾動隨時間衰減,系統(tǒng)穩(wěn)定。

    考慮到葉片力建模以及求解特征方程的復雜性,根據(jù)實際工程的要求以及當前的計算能力、設計體系,需要對穩(wěn)定性理論進行不同層次的簡化,得到不同層次的簡化模型。當馬赫數(shù)較低情況下,壓縮性對于流動本身影響不大,可以將流動主控方程簡化為三維不可壓縮Navier-Stokes方程組;對于近似軸對稱流動情況,可以重點關(guān)注子午平均流面;當徑向方向流動可以忽略時,研究二維軸對稱流動便能抓住問題的本質(zhì)和關(guān)鍵;為了進一步簡化計算量,準二維流線也可以作為研究主體。這些簡化方式要根據(jù)實際問題的需要而適當選擇,所發(fā)展的葉輪機流動穩(wěn)定性理論模型可被用于軸流壓氣機和離心壓氣機的流動穩(wěn)定性預測,并分析探討氣流壓縮性、葉尖間隙及葉片造型對壓氣機穩(wěn)定性的影響。

    該方法在中國科學院工程熱物理研究所的低速軸流壓氣機實驗臺和Rotor37穩(wěn)定性預測上得到了應用驗證。圖6針對前者低速壓氣機的計算結(jié)果顯示:隨著壓氣機截流過程,流量逐漸減小,所有9個模態(tài)的衰減因子不斷增長。周向傳播速度較慢的前5個模態(tài)比后4個模態(tài)更快地趨于衰減因子臨界值。模態(tài)1的衰減因子最先由負變正,代表了壓氣機系統(tǒng)由穩(wěn)定變?yōu)椴环€(wěn)定,因而模態(tài)1為最不穩(wěn)定模態(tài)。臨界點的質(zhì)量流量為2.324 kg/s,與實驗測得的失穩(wěn)開始點質(zhì)量流量的相對誤差為0.65%。

    Rotor37采用跨聲速葉型設計,據(jù)了解,目前國際上公開發(fā)表的文獻中鮮有研究一般壓氣機結(jié)構(gòu)下跨聲速流動失速起始的模型預測。因為在跨聲速流場中,葉柵中的激波對流動有重要的影響,使得很難在壓氣機系統(tǒng)中對其?;T摲椒ú捎弥芟蚱骄亩ǔA鲌龇植?,葉柵通道內(nèi)的激波所導致的損失和壓升都能包含進入平均流場中,并且經(jīng)過平均處理后,通道內(nèi)也不存在流場的跳躍式突變。因此,本模型有可能克服激波這一突躍條件帶來的挑戰(zhàn)。圖7為Rotor37的穩(wěn)定性預測結(jié)果,并在此基礎上發(fā)展了可以考慮葉片彎掠造型對穩(wěn)定性影響的預測分析方法。圖中:NTD代表無量綱延遲時間。

    最近幾年,本課題組選擇Rotor37為基準葉型,通過沿弦向移動基元葉型來獲得一系列的掠葉片,進而采用子午面模型計算預測穩(wěn)定裕度,并結(jié)合定常流場的計算結(jié)果,分析討論氣動掠設計對壓氣機穩(wěn)定性的影響[33]。

    如圖8所示[33],前掠并沒有對轉(zhuǎn)子的穩(wěn)定性造成顯著的改變。然而,后掠轉(zhuǎn)子相比Rotor37穩(wěn)定裕度有明顯降低。進一步而言,后掠程度越大,穩(wěn)定裕度隨之單調(diào)遞減。這一結(jié)果與過去的研究定性上是吻合的[34-36]。圖中:BS代表Backward Sweep;FS代表Forward Sweep。

    以北京航空航天大學單級低速壓氣機實驗臺TA36為算例,通過在流線曲率通流設計中調(diào)整葉片環(huán)量分布設計得到不同負荷分布的轉(zhuǎn)子葉片。采用子午面模型預測其失穩(wěn)邊界,結(jié)合定常和非定常流場結(jié)果,分析討論轉(zhuǎn)子葉片軸向負荷分布對壓氣機穩(wěn)定性的影響[37]。

    對3種負荷分布轉(zhuǎn)子的壓氣機進行失穩(wěn)邊界預測的結(jié)果展示在圖9中[37]??梢钥吹?,隨著節(jié)流過程進行,轉(zhuǎn)子前加載壓氣機最先失穩(wěn),其次為均勻加載壓氣機,轉(zhuǎn)子后加載壓氣機最后失穩(wěn)。這表明,該算例中的轉(zhuǎn)子軸向負荷后移對流動穩(wěn)定性而言是有利的。

    2 基于特征值理論的壓氣機氣彈穩(wěn)定性預測

    航空領域的氣動彈性(簡稱氣彈)穩(wěn)定性問題涉及到復雜的流動和結(jié)構(gòu)相互作用,預測這一問題需要分別對流動以及結(jié)構(gòu)建立相應的模型,然后將二者結(jié)合來對氣動彈性穩(wěn)定性進行預測。Theodorsen[38]首次給出了采用特征值理論對氣動彈性穩(wěn)定性進行判斷的方法,該方法采用勢流理論結(jié)合庫塔條件描述葉片所受到的流動作用力,同時采用剛體運動方程來描述翼型的振動,聯(lián)立二者可建立葉片的動力學方程,通過對方程進行特征值分析來判斷翼型在特定工況下是否會出現(xiàn)氣動彈性穩(wěn)定性問題。在今天看來,雖然Theodorsen的工作采用了極其簡化的模型來描述流動以及結(jié)構(gòu)的動力學特性,但是其采用特征值方法對于氣動彈性穩(wěn)定性進行判斷的思想?yún)s對隨后的研究工作具有深刻的啟示意義。

    Dowell[39]系統(tǒng)總結(jié)了對于單個孤立翼型的特征值分析方法,除了對氣動力進行理論建模以外,還可以通過實驗數(shù)據(jù)對氣動力進行擬合。Whitehead[40]曾對振蕩葉柵的升力系數(shù)進行測量以此作為葉片動力學方程中的輸入來對其氣彈穩(wěn)定性進行特征值分析。Bendiksen和Friedmann[41]采用特征值方法對一個二維葉柵彎扭復合顫振進行了穩(wěn)定性分析,Kielb和Kaza[42]通過采用二維梁以及有限元模型來對葉片構(gòu)建結(jié)構(gòu)模型,結(jié)合特征值方法對帶掠風扇葉片的顫振穩(wěn)定性進行了分析,二者對于不可壓縮流動的描述均采用Whitehead[40]所給出振蕩葉柵氣動力模型。Kileb 和Ramsey[43]則采用基于拉普拉斯變換的線化氣動力模型對葉柵氣動載荷進行建模,通過特征值方法給出對超聲速流動下葉柵的氣彈穩(wěn)定性變化趨勢。

    在早期計算能力不發(fā)達的年代,進行特征值分析的前提是需要將氣動作用力解析表達成頻域形式,而對于壓氣機內(nèi)的相關(guān)問題,葉片氣動載荷的建模多取自于外流孤立翼型中相關(guān)的研究結(jié)果,這也導致理論結(jié)果與真實現(xiàn)象存在較大的出入[44]。壓氣機葉片由合金制造,質(zhì)量比較大,此時顫振通常表現(xiàn)為單一的模態(tài),而外流機翼中多為復合模態(tài);另一方面,壓氣機內(nèi)部的流動受到葉片和流道復雜幾何、固壁以及上下游結(jié)構(gòu)的影響,其復雜程度遠高于外流問題,并且表現(xiàn)出明顯的流固耦合特性,對葉片氣動作用力進行準確的理論建模難度較大,因此適用于外流氣彈分析的理論解耦模型通常難以準確描述壓氣機內(nèi)部的真實流動情況。雖然部分基于解析方法的預測模型被用于設計階段的快速迭代,但是通常認為這一類方法預測的結(jié)果與實際問題有著較大差別。

    Carta[45]提出了一種能量法的判據(jù)來預測壓氣機葉片是否會出現(xiàn)顫振問題,這一方法所采用的氣動力及葉片結(jié)構(gòu)模型與Theodorsen的工作是類似的,但是區(qū)別在于Carta并沒有對控制方程的特征值進行求解,而是只提取了結(jié)構(gòu)方程部分的特征頻率,在此基礎上假設葉片在特征頻率處做簡諧運動,這一做法與Whitehead[46-47]的工作中對流體域與固體域做解耦處理是一致的,通過計算一個周期內(nèi)氣動作用力對葉片所做的功來判斷葉片與流動之間的能量傳遞情況,即氣動阻尼功,如果氣動力對葉片輸入負功,則認為流動對葉片的振動起到阻尼作用,此時葉片處于穩(wěn)定狀態(tài),否則即可能發(fā)生顫振。這一基于能量法的顫振判斷準則在壓氣機的氣動彈性問題預測中一直沿用至今。

    2.1 基于能量法的快速氣彈預測算法

    能量法通過指定葉片振動規(guī)律來對氣彈問題進行解耦,其核心問題在于如何對流動的非定常響應進行預測。在壓氣機設計階段,對于預測方法的主要要求便是能夠?qū)崿F(xiàn)高效迭代以完成快速設計,此時流場對于葉片非定常輸入的響應多采用簡化的解析方法來建模。Sun等[48]采用三維升力面理論來求解由葉片振動誘導產(chǎn)生的非定常載荷響應:在葉片表面建立無穿透邊界條件,要求來流擾動速度場與葉片振動速度在葉片表面法向和速度為零,通過將速度脈動轉(zhuǎn)化為壓力脈動的積分方程,如圖10所示,結(jié)合廣義聲類比積分公式,可將轉(zhuǎn)子散射的壓力脈動場寫為

    (10)

    假設壓力脈動存在時空周期性,則葉片上下表面壓差可展開成傅里葉級數(shù)的形式:

    (11)

    進一步結(jié)合格林函數(shù),散射的壓力脈動場可寫為積分方程的形式:

    (12)

    近年來,部分研究結(jié)果表明轉(zhuǎn)子進口段的長度以及聲處理可能會對轉(zhuǎn)子的顫振穩(wěn)定性產(chǎn)生較大的影響[48-49],為了考慮進口段長度以及聲襯邊界條件的影響,Sun等[48]結(jié)合傳遞單元方法,通過等價分布源方法來描述聲學邊界對于葉片顫振的非定常響應,同時采用邊界元方法來刻畫管道邊界的反射效應,成功建立了有限/無限長管道內(nèi)的流固聲耦合的預測模型。

    圖11給出了不同聲學邊界條件下轉(zhuǎn)子葉片氣動功Waero的變化結(jié)果,葉片振動呈現(xiàn)為一階彎扭復合模態(tài),當振動節(jié)徑ND=2時,氣動功的剛壁值小于0,即轉(zhuǎn)子在剛壁管道中是穩(wěn)定的,而當阻抗為Zs=(0.05,-3.9)的聲襯被加入管道后,氣動功進一步下降為一個絕對值比原剛壁值大兩倍的負數(shù),這也意味著氣流對轉(zhuǎn)子的阻尼效應增強。而在同一條件下,若聲襯阻抗為Zs=(0.05,-3.2),則聲襯反而會對振蕩轉(zhuǎn)子的氣動彈性響應產(chǎn)生負面影響,觸發(fā)葉片產(chǎn)生顫振。這一對比說明,壁面聲處理對于轉(zhuǎn)子顫振穩(wěn)定性產(chǎn)生重要的影響。

    與此同時,轉(zhuǎn)子進口管道的長度也可能對振蕩轉(zhuǎn)子的氣動阻尼產(chǎn)生重要的影響,圖12給出了不同管道長度下,進口距葉片前緣面距離長度對于振蕩轉(zhuǎn)子氣動功的影響結(jié)果。圖中:FLD和LID分別表示有限長和無限長;Lduct為管道長度;Ca為葉片軸向弦長。從圖中可以看到,基于有限長和無限長管道模型所得到的非定常氣動功存在明顯差異。當管道總長度為5倍葉片軸向弦長時,如圖12(a)所示,非定常氣動功會隨著進口段長度增加而由負值逐漸變?yōu)檎@也意味著進口段過長可能導致轉(zhuǎn)子出現(xiàn)顫振。而在圖12(b)所給出的管道長度為15倍葉片軸向弦長的情形下,當轉(zhuǎn)子位于管道中間位置時,非定常氣動功的數(shù)值最小,此時轉(zhuǎn)子葉片最穩(wěn)定,而當轉(zhuǎn)子逐漸靠近進口或者出口位置時,非定常氣動功的數(shù)值呈現(xiàn)出逐漸增大的趨勢直至變?yōu)檎?,這也意味著對于圖中所考慮的振動模態(tài),轉(zhuǎn)子過于靠近管道進口或者出口可能導致顫振的發(fā)生。

    2.2 氣流非定常載荷預測方法

    近年來,隨著計算方法以及軟硬件的發(fā)展,基于CFD方法求解更為復雜的非線性流動模型在氣流非定常載荷的預測中扮演著越來越重要的角色。在20世紀80~90年代,受制于計算能力的限制,Hall和Crawley[50]通過對非線性歐拉方程進行線化的方式來提升計算效率。線化后的歐拉方程可以在特征頻率處做周期性假設,對于氣彈問題,這一頻率通常是葉片的振動頻率,進而可以求得該特征頻率所對應的空間特征向量,特征向量的物理意義是周期性的流動在空間各處脈動的幅值。計算效率是評估壓氣機氣彈預測方法的核心指標之一,而基于特征頻率的處理手段可以避免對時間項進行迭代求解,因而可以大幅提升計算效率。當流場中不存在強非線性結(jié)構(gòu)時,比如流動分離以及激波,基于線化方程的計算方法能夠?qū)τ扇~片振動引起的非定?,F(xiàn)象進行合理的預測,而當流場中非線性作用占據(jù)主導時,頻域線化方法則不再適用。

    為了進一步考慮包含非線性的影響同時兼顧頻域方法在計算效率上的優(yōu)勢,He和Ning[51]提出了一種非線性諧波方法,這一方法的核心思想是在對Navier-Stokes方程進行線化的過程中保留擾動交叉項,以特征頻率對擾動進行周期性假設,在對擾動求解的過程中同時對包含非線性因素的時均流場進行同步修正直至結(jié)果收斂。相較于非線性時間推進,非線性諧波方法的計算效率仍然有1~2個數(shù)量級的提升?;陬愃频乃枷?,Hall等[52]提出了一種求解具有時空周期性非線性流場的諧波平衡方法,這一方法通過迭代求解特征頻率及其各階諧波所對應的空間特征分布來對流場進行時空重構(gòu),能夠?qū)崿F(xiàn)較高的非定常計算效率。非線性諧波法以及諧波平衡方法是目前研究壓氣機葉片氣彈問題的主要模擬工具。Clark[53]采用諧波平衡方法對壓氣機葉片的流動進行快速模擬,在此基礎上采用POD(Proper Orthogonal Decomposition)方法提取流場中的主要模態(tài)來構(gòu)建降階模型并嘗試將其應用于轉(zhuǎn)子非同步振動問題的預測。Wang和Huang[54]提出了一種針對諧波平衡方法的加速收斂方法并將其應用在跨聲速轉(zhuǎn)子的氣彈穩(wěn)定性預測。Huang等[55]采用諧波平衡方法結(jié)合影響系數(shù)法求解了葉片振動模態(tài)處的氣動作用力,在此基礎上結(jié)合僅包含指定模態(tài)的結(jié)構(gòu)動力學方程對葉片的顫振實現(xiàn)基于小擾動的特征值分析?;谥C波平衡的思想,Du和Ning[56]提出一種采用多項式來簡化時間項的方式,同樣能夠準確地對具有時空周期性的流動進行預測。Sandberg和Michelassi[57]與Casoni和Benini[58]針對軸流葉輪機械領域典型的數(shù)值模型以及頻域和其他各類降階模型給出了全面的總結(jié)。

    2.3 耦合計算的應用場景

    采用時空周期性假設的求解方法一般要利用特征頻率對時間項做傅里葉展開,而當流場中同時包含多個特征頻率時,如多葉片排干涉上下游葉片通過頻率不同,以及轉(zhuǎn)靜干涉的同時還存在葉片振動,這一類方法通常只能對不同頻率做線性疊加而無法考慮彼此之間的非線性耦合效應。當不同頻率線性疊加,多個頻率互相靠近產(chǎn)生“拍”現(xiàn)象時,計算收斂時間也會大幅延長[59]。這一類方法實際上是解耦的,即沒有考慮流場演化對于葉片振動頻率、振幅以及葉片間相角的影響。由于目前復合材料的大量應用,葉片質(zhì)量相較于早期已經(jīng)大幅降低,此時葉片振動頻率及模態(tài)與不考慮流固耦合效應的結(jié)果會存在差異,與此同時,葉片振動響應可能包含多個模態(tài),采用頻率解耦方法不符合真實的物理演化過程。Chahine等[60]的計算結(jié)果表明,隨著質(zhì)量比以及結(jié)構(gòu)剛性的降低,耦合效應的作用越來越顯著,由此導致耦合方法與解耦的能量法所給出的結(jié)果差異越來越大。在流固耦合效應起主導作用時,基于全非線性方程的耦合計算則顯得十分必要[61-62]。

    另一類必須采用耦合模型計算的振動問題是轉(zhuǎn)子非同步振動問題。近年來國外研究機構(gòu)在不同的實驗平臺上都曾測量到了葉片非同步振動的信號[63-65],這一現(xiàn)象一般呈現(xiàn)出高周向模態(tài)數(shù),并且伴隨著沿周向傳播的分離流動結(jié)構(gòu),同時壁面壓力信號存在多個非軸頻整數(shù)倍的異常成分。Sun[48]和Vahdati[49]等的工作都表明,進口管道的反射會對轉(zhuǎn)子風扇的非同步振動穩(wěn)定性產(chǎn)生明顯的影響。對于這一現(xiàn)象必須采用耦合處理的原因在于,從實驗結(jié)果來看,受到葉片振動與流動之間干涉效應的影響,流動的周向模態(tài)特性在振動前后會出現(xiàn)明顯的改變[64-65],即流固耦合效應在這一現(xiàn)象中同樣扮演著重要的角色。目前針對非同步振動問題已有較多的研究結(jié)果,在判斷葉片非同步振動的穩(wěn)定性時,絕大多數(shù)工作[66-69]仍然采用解耦的處理方式,即給定頻率和振幅并且使用氣動阻尼這一參數(shù)。氣動阻尼這一概念源于線性范疇,需要參照給定的葉片振幅做無量綱從而保證氣動阻尼與振幅無關(guān)而非同步振動往往伴隨著沿周向傳播的分離流動,表現(xiàn)出強非線性特征,在此情形下,振幅對于氣動阻尼計算結(jié)果的影響目前尚不明確。雖然葉片的振動頻率可能仍然是固有頻率,但是非同步振動的周向模態(tài)特性與流動密切相關(guān),同樣無法在計算前獲得,并且葉片的振動會顯著改變流動的周向模態(tài)特性,這也意味著無法準確地給出周期性邊界條件,因此大多數(shù)情況下需要開展全環(huán)計算。較于上述基于頻域假設的快速非定常求解算法,采用時間推進方式求解完全非線性的Navier-Stokes方程能夠包含更多復雜流動因素,因此也被廣泛應用于壓氣機氣彈問題預測,同時其對于計算量的要求也進一步提升。

    3 基于特征值理論的燃燒不穩(wěn)定性預測

    燃燒不穩(wěn)定性的發(fā)生是燃燒系統(tǒng)內(nèi)火焰非定常熱釋放和聲波的充分耦合導致的[70]。因此,燃燒不穩(wěn)定性的發(fā)生關(guān)鍵在于火焰非定常熱釋放和聲波的耦合,將復雜問題降階簡化,保留影響燃燒穩(wěn)定性發(fā)生的關(guān)鍵因素,即火焰對聲波的響應,以及影響聲波反射的燃燒室的阻抗壁面刻畫,在其他部分忽略次要因素,比如整個燃燒室內(nèi)復雜的湍流流動、燃燒反應等,而只考慮聲波傳播[71]。通過模型實驗和高精度數(shù)值模擬可以得到火焰對聲波的響應,在線性頻率內(nèi)得到火焰?zhèn)鬟f函數(shù)(Flame Transfer Function, FTF)[72],在非線性頻域范圍內(nèi),Dowling等發(fā)現(xiàn)燃燒不穩(wěn)定性非線性的主要因素在于火焰對不同幅值聲波響應的不同,因此可以通過建立包含聲波幅值影響的火焰函數(shù),也稱火焰描述函數(shù)(Flame Describing Function)[73-74]用于非線性燃燒不穩(wěn)性分析。至于燃燒室其他部分則通過建立聲學模型,結(jié)合火焰?zhèn)鬟f函數(shù)在不同界面上進行守恒方程匹配。

    燃燒不穩(wěn)定性分析本質(zhì)上是基于小擾動特征值理論的穩(wěn)定性分析,即分析燃燒系統(tǒng)有燃油當量比脈動、旋渦脫落及燃燒噪聲等小擾動存在的情況下,系統(tǒng)的響應特征,如果在施加起始小擾動后的一段時間內(nèi),系統(tǒng)的脈動幅值不斷放大,則表明燃燒系統(tǒng)是線性不穩(wěn)定的,反之,則是線性穩(wěn)定的。基于這樣的認識,我們可以對系統(tǒng)控制方程進行線化處理,得到關(guān)于小擾動量的控制方程,并求解其復特征頻率,其中復特征頻率的實部表示系統(tǒng)的熱聲共振頻率,虛部則代表系統(tǒng)不穩(wěn)定性的增長率,通過這種理論模型方法,可以快速進行燃燒不穩(wěn)定性的預測或研究控制手段[75-76]。

    國際上通用的聲網(wǎng)絡低階模型方法在認識燃燒不穩(wěn)定性發(fā)生機理及對燃燒室進行穩(wěn)定性預測方面發(fā)揮了重要作用,能夠為燃燒不穩(wěn)定性控制提供準確的目標,不穩(wěn)定模態(tài)及頻率。在考慮控制手段設計時,以上模型方法能夠直接在特征頻率分析中考慮具有集總參數(shù)特點的Helmholtz共振器[77-78]。但是,當需要對燃燒室中常常使用的兼具冷卻壁面和提供聲學耗散作用的含冷卻偏流穿孔板聲襯[79]等具有分布參數(shù)特點的控制壁面設計時,需要配合額外聲學計算[80]或?qū)嶒瀃81]。但是這種根據(jù)頻率、模態(tài)信息再通過聲學計算或?qū)嶒灲怦畹卦O計燃燒不穩(wěn)定性抑制手段往往會帶來誤差。這是因為,首先,穿孔板聲襯等聲學軟壁面在燃燒室內(nèi)會改變系統(tǒng)的熱聲共振頻率和模態(tài)分布,因此,在硬壁面條件下建立的預測模型不能準確預測帶有聲襯的燃燒室共振頻率和模態(tài)分布。另外,通過聲學計算或?qū)嶒炘O計聲襯過程中,聲襯面臨的聲學環(huán)境也很難完全實現(xiàn)其在燃燒室內(nèi)的環(huán)境。因此,為了研究聲襯控制熱聲不穩(wěn)定性效果,并準確地對聲襯參數(shù)進行優(yōu)化設計,建立一種能夠?qū)⒙曇r耦合考慮的燃燒不穩(wěn)定性模型就顯得十分重要而且必要了。

    3.1 耦合考慮聲襯壁面的三維燃燒不穩(wěn)定性預測模型

    將聲襯耦合考慮進燃燒不穩(wěn)定性模型中,難點在于如何描述燃燒室內(nèi)聲波和壁面的相互干涉作用。在軟壁面環(huán)境下,壁面阻抗、系統(tǒng)頻率、管道特征值形成復雜的耦合迭代關(guān)系[82]。為了解決這一難題,Sun等提出一種基于等價分布源思想[83]的傳遞單元方法[84],將聲襯穿孔板阻抗邊界視為等效連續(xù)分布的單極子聲源,這樣能夠保留管道特征值的正交特性,避免軟壁面管道徑向特征值的復雜迭代求解。傳遞單元方法,可以用以考慮有限長聲襯的聲傳播問題,也可以用以將聲襯作為聲學傳遞單元包含到燃燒不穩(wěn)定性模型。

    該方法在處理聲襯段聲學單元時,為了可以將之耦合到燃燒不穩(wěn)定性分析模型中,需要得到其傳遞矩陣,即聲襯段內(nèi)聲波關(guān)于界面入射波的顯示表達式。具體地,將聲襯軟壁面視為等價連續(xù)分布的單極子聲源,聲襯段的內(nèi)壓力擾動和速度擾動表達式可以寫為界面入射聲波和散射聲波之和:

    p′=p′I+p′d

    (13)

    其中:

    (14)

    聲襯背腔管道內(nèi)由聲襯單極子聲源產(chǎn)生的散射聲場為

    (15)

    p′c-p′=Zv′n

    (16)

    式中:Z為穿孔板聲阻抗;v′n為垂直于聲襯表面的聲質(zhì)點速度。通過阻抗方程以及聲襯前后端界面匹配條件即可求解得到等階聲質(zhì)點速度的表達式,再將之代入方程(13)和方程(14)即可以得到聲襯段內(nèi)聲場關(guān)于界面入射波的顯示表達式,即完成了聲襯段聲傳遞單元的建立詳細推導過程和表達式參考文獻[84]接著可以和燃燒室內(nèi)其他硬壁面條件下聲網(wǎng)絡單元建立起燃燒不穩(wěn)定性分析模型。關(guān)于模型建立的細節(jié)可以參考文獻[85]。應用該特征值分析模型,可以研究不同形式的聲軟壁面控制對燃燒不穩(wěn)定性的影響,并進行快速的參數(shù)優(yōu)化設計。

    從Rayleigh準則出發(fā)可知,控制燃燒不穩(wěn)定性的發(fā)生可以從兩個方面著手,一是削弱熱聲耦合強度[86];二是增加系統(tǒng)的聲學耗散[82,85]。下文將從增加系統(tǒng)聲學耗散角度研究穿孔板聲襯在環(huán)形燃燒室中控制熱聲不穩(wěn)定性來闡述本三維模型方法的應用。

    3.2 穿孔板聲襯控制環(huán)形燃燒室熱聲不穩(wěn)定性

    環(huán)形燃燒室由于其能產(chǎn)生更均勻的燃燒,具有結(jié)構(gòu)更緊湊的優(yōu)點,目前是航空發(fā)動機和地面燃氣輪機領域最為常用的燃燒室類型。為了滿足日益嚴苛的排放標準,環(huán)形燃燒室多采用貧油預混燃燒,具有較低的NOx(氮氧化物)排放,但更易出現(xiàn)燃燒不穩(wěn)定性。此外,由于燃燒溫度較高,燃燒室火焰筒壁面往往開有小孔,從高壓壓氣機后引入部分冷氣從小孔內(nèi)通過,形成小孔冷卻偏流,自然形成了具有聲學耗散的穿孔板聲襯結(jié)構(gòu)。穿孔板聲襯在滿足冷卻需求的情況下,還能提供聲學耗散在線性范圍內(nèi)抑制燃燒不穩(wěn)定性的發(fā)生,或盡可能的降低振蕩幅值。目前,由于巨大的計算資源和時間成本,還難以通過數(shù)值模擬如大渦模擬(LES)研究不同聲襯結(jié)構(gòu)的控制效果。而采用本文所介紹的包含軟壁面聲襯的燃燒不穩(wěn)定性模型則可以進行聲襯控制效果的研究以及參數(shù)的優(yōu)化。

    通過與有限元分析軟件COMSOL對比無熱釋放條件下系統(tǒng)的特征頻率及模態(tài)分布情況,驗證本模型的有效性,頻率對比如圖14所示。圖中模態(tài)名稱中,L表示軸向模態(tài),A表示周向模態(tài),字母前數(shù)字表示模態(tài)階數(shù),如2L0A代表第2階軸向模態(tài)混合0階周向模態(tài)。從圖中計算結(jié)果還可以發(fā)現(xiàn),聲襯的引入使得必須考慮徑向耦合模態(tài),才能準確預測頻率以及模態(tài)分布。也就是說,當我們考慮包含聲軟壁面的燃燒不穩(wěn)性預測時,有必要發(fā)展三維模型,才能夠準確預測系統(tǒng)的頻率及增長率。

    為了研究聲襯控制效果背后機理,給出系統(tǒng)內(nèi)聲壓分布,對于模態(tài)2L0A聲襯控制效果最好的聲襯位置大小L2/Lc=0,L3/Lc=0.47,系統(tǒng)內(nèi)聲壓分布如圖15(a)和圖15(b)所示。圖中:

    以上研究結(jié)果表明,對于包含聲襯的燃燒室系統(tǒng),有必要建立三維燃燒不穩(wěn)定性模型通過特征值分析理論,預測系統(tǒng)的穩(wěn)定性并研究聲襯的控制效果。聲襯的大小和位置對于其控制燃燒不穩(wěn)定性效果有重要影響,可以通過改變聲襯背腔結(jié)構(gòu)提高聲襯的控制效果。另外,值得一提的是,在實際系統(tǒng)中,由于火焰的非線性,還可以通過本模型方法結(jié)合非線性火焰函數(shù)模型[88],以研究不同聲襯結(jié)構(gòu)對降低熱聲不穩(wěn)定性極限環(huán)幅值的影響,以滿足工程實際需求。

    4 總結(jié)與展望

    隨著航空發(fā)動機設計理論和三維數(shù)值模擬技術(shù)的不斷發(fā)展,針對發(fā)動機內(nèi)流、氣彈和燃燒領域面臨的工程技術(shù)問題,各種理論模型和實驗數(shù)據(jù)庫也更加完善,設計人員普遍習慣于通過所積累的設計經(jīng)驗來預估和判斷風扇/壓氣機以及主燃/加力燃燒室的穩(wěn)定工作范圍,但基于經(jīng)驗的評估方法很難被納入到一體化設計工具中并對設計方向給予合理的引導,因此,快速準確的穩(wěn)定性評估工具對于航空發(fā)動機的設計至關(guān)重要。本文簡要闡述了基于小擾動方法的特征值理論在航空發(fā)動機氣動、氣彈和燃燒穩(wěn)定性預測方面的應用,對從經(jīng)典的線化小擾動模型到發(fā)展較為成熟的半解析模型方法的研究發(fā)展進行了綜述,并重點以本課題組近年來的研究結(jié)果為例,給出了該方法在流動失穩(wěn)預測、葉輪機設計和燃燒不穩(wěn)定性預測方面的應用案例。實踐表明,該方法可以作為發(fā)動機一體化設計分析工具成為發(fā)動機穩(wěn)定性設計過程中的一環(huán)。

    此外,鑒于壁面邊界條件可以影響小擾動在動力學系統(tǒng)內(nèi)的演化行為,在本文所介紹的理論模型基礎上,相繼有通過壁面阻抗調(diào)控方法控制流動失穩(wěn)和熱聲不穩(wěn)定性的技術(shù)手段被發(fā)展,為失速喘振和燃燒不穩(wěn)定性的控制提供了新的實現(xiàn)手段。面向未來更先進航空發(fā)動機研制的需求,基于特征值理論的穩(wěn)定性模型方法需要考慮更加復雜的壁面邊界條件、幾何影響因素、更多非定常擾動干涉機制,以及非線性因素等,以快速低成本的方法優(yōu)勢支撐行業(yè)發(fā)展可靠、高效的穩(wěn)定性評估工具。

    猜你喜歡
    壓氣機特征值擾動
    Bernoulli泛函上典則酉對合的擾動
    軸流壓氣機效率評定方法
    一類帶強制位勢的p-Laplace特征值問題
    重型燃氣輪機壓氣機第一級轉(zhuǎn)子葉片斷裂分析
    單圈圖關(guān)聯(lián)矩陣的特征值
    壓氣機緊湊S形過渡段內(nèi)周向彎靜子性能數(shù)值計算
    (h)性質(zhì)及其擾動
    小噪聲擾動的二維擴散的極大似然估計
    基于商奇異值分解的一類二次特征值反問題
    用于光伏MPPT中的模糊控制占空比擾動法
    国产成人午夜福利电影在线观看| 久久人妻熟女aⅴ| 欧美激情极品国产一区二区三区| 亚洲人成电影观看| 美女扒开内裤让男人捅视频| av不卡在线播放| 纯流量卡能插随身wifi吗| 少妇被粗大的猛进出69影院| 制服丝袜香蕉在线| 亚洲av男天堂| 欧美日韩国产mv在线观看视频| 热re99久久精品国产66热6| 两个人免费观看高清视频| 免费观看av网站的网址| 十分钟在线观看高清视频www| 久久韩国三级中文字幕| 狠狠精品人妻久久久久久综合| kizo精华| 最新的欧美精品一区二区| 免费久久久久久久精品成人欧美视频| 五月天丁香电影| 免费黄频网站在线观看国产| 国产av码专区亚洲av| 人人妻人人澡人人看| 丰满迷人的少妇在线观看| 久热爱精品视频在线9| 搡老岳熟女国产| 丝袜美腿诱惑在线| 亚洲av男天堂| 又黄又粗又硬又大视频| 欧美av亚洲av综合av国产av | 制服人妻中文乱码| 久久久久国产精品人妻一区二区| 人妻 亚洲 视频| 国产成人精品久久久久久| 色吧在线观看| 夜夜骑夜夜射夜夜干| 十八禁高潮呻吟视频| 成人亚洲精品一区在线观看| 国产国语露脸激情在线看| av在线播放精品| 婷婷成人精品国产| 国产探花极品一区二区| 国产 一区精品| 久久人人爽av亚洲精品天堂| 蜜桃国产av成人99| 日韩,欧美,国产一区二区三区| 国产黄色免费在线视频| av天堂久久9| 一边摸一边做爽爽视频免费| 欧美老熟妇乱子伦牲交| 国产精品偷伦视频观看了| 久久久久人妻精品一区果冻| 韩国高清视频一区二区三区| 中文字幕人妻丝袜一区二区 | 熟妇人妻不卡中文字幕| 波野结衣二区三区在线| 操美女的视频在线观看| 又黄又粗又硬又大视频| 日韩电影二区| 日韩av不卡免费在线播放| 男女边吃奶边做爰视频| 一级黄片播放器| 国产不卡av网站在线观看| 婷婷色麻豆天堂久久| 久久久精品94久久精品| 亚洲精品第二区| 亚洲国产看品久久| 男女边摸边吃奶| 激情五月婷婷亚洲| 国产日韩欧美视频二区| 亚洲久久久国产精品| 香蕉国产在线看| 人人妻人人添人人爽欧美一区卜| 中文字幕制服av| 伦理电影大哥的女人| 国产国语露脸激情在线看| 亚洲婷婷狠狠爱综合网| 亚洲欧美色中文字幕在线| 精品免费久久久久久久清纯 | 国产老妇伦熟女老妇高清| 美女中出高潮动态图| 亚洲av日韩精品久久久久久密 | 免费看不卡的av| 性高湖久久久久久久久免费观看| 日本wwww免费看| 亚洲成人手机| 日日啪夜夜爽| 国产成人精品无人区| 国产男女内射视频| 青春草国产在线视频| 欧美黑人欧美精品刺激| 操出白浆在线播放| 高清欧美精品videossex| 午夜影院在线不卡| 男女边摸边吃奶| 9色porny在线观看| 亚洲av综合色区一区| 日韩,欧美,国产一区二区三区| 久久精品久久久久久噜噜老黄| 亚洲国产欧美网| 一个人免费看片子| 日韩一本色道免费dvd| 久久午夜综合久久蜜桃| 国产一级毛片在线| 性高湖久久久久久久久免费观看| 一级毛片 在线播放| 老司机亚洲免费影院| 日韩电影二区| 青草久久国产| 91精品国产国语对白视频| 国产在线一区二区三区精| 在线天堂最新版资源| 色吧在线观看| 精品一区二区三区av网在线观看 | 青草久久国产| 777米奇影视久久| 亚洲,一卡二卡三卡| 久久久久精品性色| 嫩草影视91久久| 国产午夜精品一二区理论片| av有码第一页| 丝袜人妻中文字幕| 最近最新中文字幕免费大全7| 欧美日韩视频高清一区二区三区二| 欧美人与性动交α欧美软件| 国产精品国产三级国产专区5o| 极品少妇高潮喷水抽搐| 91精品三级在线观看| 一级,二级,三级黄色视频| 成人黄色视频免费在线看| 精品酒店卫生间| 久久久亚洲精品成人影院| 免费久久久久久久精品成人欧美视频| 人妻人人澡人人爽人人| 国产老妇伦熟女老妇高清| 亚洲婷婷狠狠爱综合网| 久久精品久久久久久久性| 国产精品香港三级国产av潘金莲 | 成年动漫av网址| 国产免费现黄频在线看| 欧美日韩视频精品一区| 免费观看a级毛片全部| www.熟女人妻精品国产| 黑人欧美特级aaaaaa片| 亚洲精品日本国产第一区| 久久人人97超碰香蕉20202| avwww免费| 国产成人精品在线电影| 大香蕉久久成人网| 亚洲精品乱久久久久久| 亚洲色图综合在线观看| 国产精品99久久99久久久不卡 | 日韩精品有码人妻一区| 亚洲 欧美一区二区三区| a 毛片基地| 人妻一区二区av| 午夜老司机福利片| 一级a爱视频在线免费观看| 亚洲在久久综合| 国产精品秋霞免费鲁丝片| 在现免费观看毛片| 亚洲精品美女久久久久99蜜臀 | 国产在视频线精品| 18禁国产床啪视频网站| 极品人妻少妇av视频| 久久久久精品性色| 欧美日本中文国产一区发布| 久久久久视频综合| 久久精品久久久久久久性| 午夜福利在线免费观看网站| 欧美另类一区| 亚洲美女搞黄在线观看| 一级片免费观看大全| 人人妻人人添人人爽欧美一区卜| 国产日韩欧美在线精品| 国产一区二区三区av在线| 9191精品国产免费久久| 国产午夜精品一二区理论片| 人人妻人人澡人人看| 一边摸一边抽搐一进一出视频| 亚洲男人天堂网一区| e午夜精品久久久久久久| 狂野欧美激情性bbbbbb| 看免费成人av毛片| 啦啦啦在线观看免费高清www| 欧美成人午夜精品| 黑丝袜美女国产一区| www.精华液| 女人被躁到高潮嗷嗷叫费观| 成年美女黄网站色视频大全免费| 校园人妻丝袜中文字幕| 人人妻,人人澡人人爽秒播 | 精品亚洲成国产av| 宅男免费午夜| 美女大奶头黄色视频| 欧美少妇被猛烈插入视频| 亚洲欧美一区二区三区国产| 亚洲欧美成人精品一区二区| 国产精品 国内视频| 自拍欧美九色日韩亚洲蝌蚪91| 国产日韩欧美视频二区| 久久人妻熟女aⅴ| 国产精品女同一区二区软件| 精品一区二区三卡| 久久久久精品久久久久真实原创| 侵犯人妻中文字幕一二三四区| 大香蕉久久网| av一本久久久久| 韩国高清视频一区二区三区| 丝袜喷水一区| 多毛熟女@视频| 一二三四在线观看免费中文在| 国产xxxxx性猛交| 亚洲成人av在线免费| 成人国产麻豆网| 亚洲欧美一区二区三区久久| 久久久国产欧美日韩av| 精品福利永久在线观看| 欧美国产精品一级二级三级| 精品国产超薄肉色丝袜足j| 精品酒店卫生间| 69精品国产乱码久久久| 黑人猛操日本美女一级片| 欧美 亚洲 国产 日韩一| 青春草视频在线免费观看| 国产黄频视频在线观看| av卡一久久| 国产精品久久久av美女十八| 这个男人来自地球电影免费观看 | 亚洲精品在线美女| 久久女婷五月综合色啪小说| 美女扒开内裤让男人捅视频| 日韩 亚洲 欧美在线| h视频一区二区三区| 国产探花极品一区二区| 日本色播在线视频| av线在线观看网站| 看十八女毛片水多多多| 观看av在线不卡| 久久国产亚洲av麻豆专区| 成年动漫av网址| 亚洲美女视频黄频| 久久久久久免费高清国产稀缺| 久久精品国产亚洲av高清一级| 久久久久久久精品精品| 成人午夜精彩视频在线观看| 亚洲免费av在线视频| 另类精品久久| 日韩人妻精品一区2区三区| 大片免费播放器 马上看| 亚洲欧美成人精品一区二区| 亚洲精品乱久久久久久| 黑人巨大精品欧美一区二区蜜桃| 青草久久国产| 免费黄色在线免费观看| 国产男女超爽视频在线观看| av又黄又爽大尺度在线免费看| 久久久久网色| 国产男女超爽视频在线观看| 亚洲,一卡二卡三卡| 国产精品国产av在线观看| av片东京热男人的天堂| 纯流量卡能插随身wifi吗| 亚洲av中文av极速乱| 中文字幕精品免费在线观看视频| 亚洲av日韩精品久久久久久密 | av一本久久久久| 一个人免费看片子| 日韩制服骚丝袜av| 欧美亚洲日本最大视频资源| 十八禁人妻一区二区| 亚洲四区av| 成人免费观看视频高清| 黑人欧美特级aaaaaa片| 久久婷婷青草| 在线观看三级黄色| 丝袜喷水一区| 日日撸夜夜添| 亚洲美女搞黄在线观看| 国产精品无大码| 丝袜脚勾引网站| 美女脱内裤让男人舔精品视频| 日本av手机在线免费观看| 日韩欧美一区视频在线观看| 18禁裸乳无遮挡动漫免费视频| 欧美日韩亚洲高清精品| 91aial.com中文字幕在线观看| 欧美国产精品va在线观看不卡| 在线观看免费视频网站a站| 大香蕉久久成人网| 在线观看免费高清a一片| 亚洲美女搞黄在线观看| 亚洲人成77777在线视频| 日韩制服骚丝袜av| 天天躁日日躁夜夜躁夜夜| 国产成人系列免费观看| 永久免费av网站大全| 自拍欧美九色日韩亚洲蝌蚪91| 秋霞在线观看毛片| 久久毛片免费看一区二区三区| 两性夫妻黄色片| 丁香六月欧美| 日本91视频免费播放| 欧美日韩综合久久久久久| 日韩大片免费观看网站| 国产精品国产三级国产专区5o| 大陆偷拍与自拍| 久久99热这里只频精品6学生| 男人添女人高潮全过程视频| 两性夫妻黄色片| 嫩草影院入口| 国产xxxxx性猛交| 国产极品粉嫩免费观看在线| 老汉色∧v一级毛片| 汤姆久久久久久久影院中文字幕| 在线免费观看不下载黄p国产| 免费女性裸体啪啪无遮挡网站| 男男h啪啪无遮挡| 在线精品无人区一区二区三| 人人妻人人添人人爽欧美一区卜| 欧美人与善性xxx| 中文字幕最新亚洲高清| 亚洲精品一二三| 国产精品欧美亚洲77777| 国产精品久久久av美女十八| 午夜老司机福利片| 老司机影院成人| 精品少妇久久久久久888优播| 久久 成人 亚洲| 国产亚洲精品第一综合不卡| 一级毛片黄色毛片免费观看视频| 一本色道久久久久久精品综合| 亚洲自偷自拍图片 自拍| 美女中出高潮动态图| 欧美日韩成人在线一区二区| 丁香六月欧美| 欧美精品高潮呻吟av久久| 日韩av在线免费看完整版不卡| 亚洲欧洲精品一区二区精品久久久 | 国产精品一区二区在线不卡| 人人妻,人人澡人人爽秒播 | av福利片在线| 久久久久久人妻| 观看av在线不卡| 国产一级毛片在线| 精品国产乱码久久久久久小说| 国产精品久久久久久久久免| 看非洲黑人一级黄片| 欧美成人精品欧美一级黄| 国产 一区精品| 久久人人爽人人片av| 久久天躁狠狠躁夜夜2o2o | 成人午夜精彩视频在线观看| 亚洲精品国产色婷婷电影| 成人午夜精彩视频在线观看| 在线观看国产h片| 亚洲欧美成人精品一区二区| 欧美日韩视频高清一区二区三区二| 亚洲国产精品999| 亚洲av电影在线观看一区二区三区| 久久天堂一区二区三区四区| 青青草视频在线视频观看| 欧美黄色片欧美黄色片| 在线观看国产h片| 久久亚洲国产成人精品v| 青春草亚洲视频在线观看| 国产午夜精品一二区理论片| 欧美久久黑人一区二区| 午夜福利,免费看| 国产日韩一区二区三区精品不卡| 亚洲精品美女久久久久99蜜臀 | 狠狠精品人妻久久久久久综合| av福利片在线| 亚洲伊人色综图| 99热网站在线观看| videos熟女内射| 亚洲成国产人片在线观看| 国产精品麻豆人妻色哟哟久久| 成年女人毛片免费观看观看9 | e午夜精品久久久久久久| 老司机在亚洲福利影院| 久久久久久久久久久免费av| 一区二区三区乱码不卡18| 精品国产乱码久久久久久男人| 亚洲综合精品二区| 在线天堂最新版资源| 免费在线观看视频国产中文字幕亚洲 | 精品人妻一区二区三区麻豆| 国产一区亚洲一区在线观看| 久热这里只有精品99| 亚洲国产中文字幕在线视频| av网站免费在线观看视频| 精品国产国语对白av| 国产高清不卡午夜福利| 天天躁日日躁夜夜躁夜夜| 欧美激情高清一区二区三区 | 成人手机av| 国产一级毛片在线| 操出白浆在线播放| 久久精品国产a三级三级三级| 亚洲av男天堂| 亚洲成人免费av在线播放| 亚洲免费av在线视频| 在线观看一区二区三区激情| 男女边吃奶边做爰视频| 国产又色又爽无遮挡免| 亚洲av电影在线观看一区二区三区| 亚洲,欧美,日韩| 久久青草综合色| 国产免费现黄频在线看| 国语对白做爰xxxⅹ性视频网站| 香蕉丝袜av| 高清不卡的av网站| 午夜福利视频在线观看免费| 如何舔出高潮| 国产伦人伦偷精品视频| av天堂久久9| 免费黄频网站在线观看国产| 欧美日韩福利视频一区二区| 亚洲精品在线美女| 日韩一本色道免费dvd| www日本在线高清视频| 精品国产乱码久久久久久男人| av卡一久久| 精品人妻熟女毛片av久久网站| 成年女人毛片免费观看观看9 | 成人影院久久| 精品视频人人做人人爽| 日韩精品有码人妻一区| 一级a爱视频在线免费观看| 成人免费观看视频高清| 久久久久精品国产欧美久久久 | 欧美日韩av久久| 亚洲精品久久久久久婷婷小说| 亚洲av综合色区一区| 少妇人妻 视频| 亚洲成人免费av在线播放| 熟女少妇亚洲综合色aaa.| 亚洲精品一二三| 黑人欧美特级aaaaaa片| 人妻 亚洲 视频| 亚洲av成人精品一二三区| 亚洲自偷自拍图片 自拍| 欧美激情 高清一区二区三区| 日本一区二区免费在线视频| 亚洲精品中文字幕在线视频| 制服诱惑二区| 国产精品国产三级国产专区5o| 亚洲一级一片aⅴ在线观看| 婷婷色av中文字幕| 操美女的视频在线观看| 精品视频人人做人人爽| 欧美久久黑人一区二区| 亚洲人成77777在线视频| 色婷婷久久久亚洲欧美| 大香蕉久久成人网| 国产在线一区二区三区精| 中文精品一卡2卡3卡4更新| 国产成人午夜福利电影在线观看| 日日摸夜夜添夜夜爱| 在线天堂最新版资源| 97人妻天天添夜夜摸| 日本av免费视频播放| www.自偷自拍.com| 国产男女内射视频| 亚洲在久久综合| 午夜91福利影院| 人成视频在线观看免费观看| 天堂俺去俺来也www色官网| 少妇被粗大猛烈的视频| 天天影视国产精品| 亚洲av电影在线观看一区二区三区| 美女高潮到喷水免费观看| 老司机深夜福利视频在线观看 | 欧美日韩亚洲综合一区二区三区_| 午夜激情av网站| 男人操女人黄网站| tube8黄色片| 欧美国产精品一级二级三级| 午夜福利乱码中文字幕| 大片电影免费在线观看免费| 麻豆精品久久久久久蜜桃| 欧美av亚洲av综合av国产av | 女性被躁到高潮视频| 国产高清不卡午夜福利| 久久 成人 亚洲| 欧美黑人欧美精品刺激| 搡老乐熟女国产| 国产又色又爽无遮挡免| av一本久久久久| av视频免费观看在线观看| 夫妻午夜视频| 男女无遮挡免费网站观看| 又粗又硬又长又爽又黄的视频| 一级爰片在线观看| 成人国语在线视频| 高清在线视频一区二区三区| 日韩不卡一区二区三区视频在线| 超碰97精品在线观看| 亚洲人成电影观看| 日韩中文字幕欧美一区二区 | 亚洲美女视频黄频| 男人添女人高潮全过程视频| 免费观看性生交大片5| 免费观看av网站的网址| 啦啦啦啦在线视频资源| 亚洲成国产人片在线观看| 熟女av电影| 欧美日韩国产mv在线观看视频| 亚洲激情五月婷婷啪啪| 18禁观看日本| 成人三级做爰电影| 在线天堂中文资源库| 老司机亚洲免费影院| 午夜av观看不卡| tube8黄色片| 亚洲情色 制服丝袜| 久久国产精品大桥未久av| 少妇人妻精品综合一区二区| 亚洲国产av影院在线观看| 日韩一区二区视频免费看| 黄色毛片三级朝国网站| 母亲3免费完整高清在线观看| 悠悠久久av| 色综合欧美亚洲国产小说| 十八禁网站网址无遮挡| 国产精品嫩草影院av在线观看| 国产精品一国产av| 国产一区有黄有色的免费视频| 91老司机精品| 久久精品国产亚洲av高清一级| 国产一级毛片在线| 久久久国产一区二区| 国产在线视频一区二区| 曰老女人黄片| 视频区图区小说| 80岁老熟妇乱子伦牲交| 精品国产国语对白av| 国产极品天堂在线| 无遮挡黄片免费观看| 亚洲精品国产色婷婷电影| 成人亚洲欧美一区二区av| 国产亚洲av高清不卡| 免费不卡黄色视频| 亚洲国产精品999| 国产精品久久久久成人av| 国产 一区精品| 男人添女人高潮全过程视频| 国产1区2区3区精品| 久久免费观看电影| 亚洲av在线观看美女高潮| av天堂久久9| 一级毛片黄色毛片免费观看视频| 亚洲av综合色区一区| 色综合欧美亚洲国产小说| 国产在线免费精品| 亚洲欧美精品综合一区二区三区| 性少妇av在线| 深夜精品福利| 亚洲一区中文字幕在线| 精品国产乱码久久久久久小说| 久久久精品免费免费高清| 精品国产国语对白av| 久久精品亚洲av国产电影网| 一级爰片在线观看| 亚洲精品av麻豆狂野| 亚洲国产欧美一区二区综合| 日日啪夜夜爽| 午夜免费鲁丝| 精品福利永久在线观看| 亚洲伊人色综图| 国产精品女同一区二区软件| 久久99一区二区三区| 嫩草影视91久久| 成人18禁高潮啪啪吃奶动态图| 成年美女黄网站色视频大全免费| 在线亚洲精品国产二区图片欧美| 男女边吃奶边做爰视频| 最近中文字幕2019免费版| 另类精品久久| 日本猛色少妇xxxxx猛交久久| 夫妻性生交免费视频一级片| 精品国产超薄肉色丝袜足j| av网站在线播放免费| 国产免费现黄频在线看| 亚洲综合色网址| 国产精品.久久久| 亚洲五月色婷婷综合| 爱豆传媒免费全集在线观看| 岛国毛片在线播放| 欧美成人午夜精品| 亚洲av福利一区| 看十八女毛片水多多多| 青春草国产在线视频| 啦啦啦在线观看免费高清www| 色94色欧美一区二区| 国产乱来视频区| 狂野欧美激情性xxxx| 女人久久www免费人成看片| 国产精品国产三级国产专区5o| 亚洲av福利一区| 午夜福利在线免费观看网站| 国产无遮挡羞羞视频在线观看| 国产精品免费大片| 国产精品av久久久久免费| 亚洲精品国产区一区二| 欧美少妇被猛烈插入视频| 免费观看av网站的网址| 精品国产超薄肉色丝袜足j| 亚洲国产中文字幕在线视频| 精品国产一区二区三区久久久樱花| 亚洲国产欧美日韩在线播放| 国产片内射在线| 在线看a的网站|