• <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中的模糊控制占空比擾動法
    国产av精品麻豆| 国产探花极品一区二区| 国产永久视频网站| 日韩制服骚丝袜av| 性色avwww在线观看| 免费看av在线观看网站| 久久精品国产亚洲av涩爱| 亚洲精品乱码久久久久久按摩| 观看av在线不卡| 国产男女内射视频| 我要看黄色一级片免费的| 香蕉精品网在线| 国产伦理片在线播放av一区| 街头女战士在线观看网站| 男女边吃奶边做爰视频| 五月天丁香电影| 老师上课跳d突然被开到最大视频| 性高湖久久久久久久久免费观看| 自拍偷自拍亚洲精品老妇| 国国产精品蜜臀av免费| a 毛片基地| 在线观看一区二区三区激情| 你懂的网址亚洲精品在线观看| 日韩精品有码人妻一区| 午夜免费鲁丝| 亚洲精品国产av蜜桃| 欧美激情国产日韩精品一区| 欧美精品一区二区免费开放| 国精品久久久久久国模美| 啦啦啦啦在线视频资源| 亚洲国产日韩一区二区| 欧美日韩精品成人综合77777| 偷拍熟女少妇极品色| 亚洲欧美日韩另类电影网站 | 人人妻人人澡人人爽人人夜夜| 国产午夜精品久久久久久一区二区三区| 激情五月婷婷亚洲| 水蜜桃什么品种好| 午夜免费观看性视频| 成人毛片60女人毛片免费| 国产欧美另类精品又又久久亚洲欧美| 国产精品久久久久久精品电影小说 | 亚洲三级黄色毛片| 国产爱豆传媒在线观看| 男女边吃奶边做爰视频| 干丝袜人妻中文字幕| 亚洲欧美成人精品一区二区| 精品人妻视频免费看| 99九九线精品视频在线观看视频| 久久久久久久大尺度免费视频| 亚洲,欧美,日韩| 免费人妻精品一区二区三区视频| 久久97久久精品| 欧美成人精品欧美一级黄| 精品少妇黑人巨大在线播放| 日韩中字成人| 国产成人a∨麻豆精品| 亚洲不卡免费看| 伦理电影大哥的女人| 久久精品国产自在天天线| 免费观看性生交大片5| 国产无遮挡羞羞视频在线观看| 欧美日本视频| 亚洲熟女精品中文字幕| 精品少妇久久久久久888优播| 欧美3d第一页| 亚洲欧美日韩东京热| 日本色播在线视频| 一区二区三区免费毛片| 亚洲第一av免费看| 亚洲av中文av极速乱| 伦理电影免费视频| 国产精品蜜桃在线观看| 九草在线视频观看| 我的女老师完整版在线观看| 美女cb高潮喷水在线观看| av在线老鸭窝| 中文天堂在线官网| 中国国产av一级| 建设人人有责人人尽责人人享有的 | 亚洲精品乱码久久久久久按摩| 免费av中文字幕在线| 欧美xxxx性猛交bbbb| 免费看不卡的av| 日产精品乱码卡一卡2卡三| 舔av片在线| 婷婷色av中文字幕| av福利片在线观看| 国产亚洲欧美精品永久| 欧美xxxx性猛交bbbb| 日日摸夜夜添夜夜添av毛片| 少妇人妻 视频| videossex国产| 观看免费一级毛片| 交换朋友夫妻互换小说| 国产综合精华液| 国产精品伦人一区二区| 人体艺术视频欧美日本| 97在线视频观看| 一级黄片播放器| 人人妻人人爽人人添夜夜欢视频 | 丰满人妻一区二区三区视频av| 一本久久精品| 80岁老熟妇乱子伦牲交| 亚洲成人中文字幕在线播放| 一本久久精品| 国产精品不卡视频一区二区| 最新中文字幕久久久久| 日韩亚洲欧美综合| 在线观看免费视频网站a站| 亚洲精品视频女| 亚洲国产精品一区三区| 伊人久久国产一区二区| 22中文网久久字幕| 嫩草影院新地址| 免费黄网站久久成人精品| 夫妻性生交免费视频一级片| 国产精品免费大片| 蜜桃亚洲精品一区二区三区| 日本免费在线观看一区| 久久久久人妻精品一区果冻| 久久精品国产a三级三级三级| 美女脱内裤让男人舔精品视频| 大码成人一级视频| 久久久久网色| 91狼人影院| 国产一区二区在线观看日韩| 国内精品宾馆在线| 人人妻人人爽人人添夜夜欢视频 | 久久人人爽人人爽人人片va| 日韩一区二区三区影片| 国产精品一区二区在线观看99| 丰满人妻一区二区三区视频av| 男人爽女人下面视频在线观看| 国产在线男女| 少妇人妻 视频| 2021少妇久久久久久久久久久| 精品一品国产午夜福利视频| 一区二区三区精品91| 亚洲欧美日韩卡通动漫| 日韩电影二区| 日产精品乱码卡一卡2卡三| 久久青草综合色| 久久久久久久久久久丰满| 日本免费在线观看一区| 舔av片在线| 狠狠精品人妻久久久久久综合| 色网站视频免费| 国产男女内射视频| 亚洲精品456在线播放app| 黑丝袜美女国产一区| 五月玫瑰六月丁香| 狂野欧美白嫩少妇大欣赏| 国产精品99久久久久久久久| 国产成人免费观看mmmm| 久久人人爽人人片av| 亚洲欧美日韩东京热| 七月丁香在线播放| 麻豆乱淫一区二区| 黄色一级大片看看| 日韩av免费高清视频| 国产 精品1| 伊人久久精品亚洲午夜| 一级片'在线观看视频| a级毛色黄片| 少妇高潮的动态图| 精品一区在线观看国产| 新久久久久国产一级毛片| 国内少妇人妻偷人精品xxx网站| 国产精品精品国产色婷婷| 内地一区二区视频在线| 午夜日本视频在线| 丰满少妇做爰视频| 久久av网站| 国产精品秋霞免费鲁丝片| 日韩亚洲欧美综合| 亚洲中文av在线| 亚洲精品色激情综合| 韩国av在线不卡| 欧美最新免费一区二区三区| 偷拍熟女少妇极品色| 又爽又黄a免费视频| 日本-黄色视频高清免费观看| 少妇人妻 视频| 一级毛片aaaaaa免费看小| 色哟哟·www| 亚洲精品乱码久久久久久按摩| 国产成人91sexporn| 欧美日韩视频高清一区二区三区二| 人妻少妇偷人精品九色| 国产亚洲av片在线观看秒播厂| 久久久久久久久久成人| 男女边摸边吃奶| 各种免费的搞黄视频| 卡戴珊不雅视频在线播放| 欧美丝袜亚洲另类| 国产又色又爽无遮挡免| 久久久久久久久久久免费av| 色网站视频免费| 中文欧美无线码| 亚洲精品美女久久av网站| 天天操日日干夜夜撸| 黄色怎么调成土黄色| 国产亚洲精品第一综合不卡| 国产一区二区激情短视频 | 熟女少妇亚洲综合色aaa.| 久久精品久久久久久久性| 成人亚洲欧美一区二区av| 伦理电影免费视频| 国产免费又黄又爽又色| 观看av在线不卡| 亚洲一区中文字幕在线| 最近中文字幕2019免费版| 黑人猛操日本美女一级片| 国产深夜福利视频在线观看| 日韩av在线免费看完整版不卡| 丁香六月欧美| 欧美激情 高清一区二区三区| 七月丁香在线播放| 免费在线观看影片大全网站 | 国产成人精品久久二区二区免费| 熟女av电影| 久久天堂一区二区三区四区| 又大又黄又爽视频免费| 国产精品国产三级国产专区5o| 捣出白浆h1v1| 国产亚洲午夜精品一区二区久久| av线在线观看网站| 69精品国产乱码久久久| 一级a爱视频在线免费观看| 日韩熟女老妇一区二区性免费视频| 熟女少妇亚洲综合色aaa.| 久久久久久久国产电影| 久久性视频一级片| 国产精品99久久99久久久不卡| 女性生殖器流出的白浆| 黄色a级毛片大全视频| 国产xxxxx性猛交| 亚洲人成电影免费在线| 一级黄片播放器| 黑丝袜美女国产一区| 夫妻午夜视频| 欧美精品高潮呻吟av久久| 麻豆国产av国片精品| 女性被躁到高潮视频| 国产一区二区三区综合在线观看| 亚洲av电影在线观看一区二区三区| 国产日韩欧美亚洲二区| 国产精品国产三级专区第一集| 一区二区三区乱码不卡18| a级片在线免费高清观看视频| 国产成人欧美在线观看 | 亚洲免费av在线视频| 亚洲中文字幕日韩| 免费女性裸体啪啪无遮挡网站| 看免费av毛片| 国产91精品成人一区二区三区 | 国产一区二区三区av在线| a级毛片黄视频| 成人免费观看视频高清| 国产精品一区二区免费欧美 | 精品福利永久在线观看| 一本久久精品| 成在线人永久免费视频| 色视频在线一区二区三区| 在线观看免费视频网站a站| 少妇人妻 视频| 亚洲精品在线美女| 亚洲人成电影观看| 欧美国产精品va在线观看不卡| 在线天堂中文资源库| 一级片'在线观看视频| 一边摸一边做爽爽视频免费| 午夜福利,免费看| 久久性视频一级片| 亚洲欧美日韩高清在线视频 | 18在线观看网站| 青青草视频在线视频观看| 欧美成人精品欧美一级黄| 成年美女黄网站色视频大全免费| 制服人妻中文乱码| 十八禁高潮呻吟视频| 婷婷色麻豆天堂久久| 中文字幕av电影在线播放| 91成人精品电影| 大话2 男鬼变身卡| 欧美大码av| 丝瓜视频免费看黄片| 男女免费视频国产| 激情五月婷婷亚洲| 狠狠精品人妻久久久久久综合| 大话2 男鬼变身卡| 国产精品免费大片| 妹子高潮喷水视频| 大话2 男鬼变身卡| 人人澡人人妻人| 久久久久久久国产电影| 国产精品.久久久| 免费不卡黄色视频| 欧美国产精品va在线观看不卡| 久久精品aⅴ一区二区三区四区| 亚洲精品日韩在线中文字幕| 国产成人av激情在线播放| 91字幕亚洲| 亚洲欧美精品自产自拍| 亚洲欧洲日产国产| 一边摸一边做爽爽视频免费| 欧美日本中文国产一区发布| 少妇人妻 视频| 无限看片的www在线观看| 一区二区日韩欧美中文字幕| 国产黄频视频在线观看| 啦啦啦在线免费观看视频4| videos熟女内射| 久久免费观看电影| 日日爽夜夜爽网站| 久久久精品免费免费高清| 女人精品久久久久毛片| 午夜免费成人在线视频| 少妇人妻久久综合中文| 精品人妻1区二区| 黄色片一级片一级黄色片| 国产成人欧美| 国产成人影院久久av| 曰老女人黄片| 亚洲黑人精品在线| 日韩 亚洲 欧美在线| 国产有黄有色有爽视频| 无限看片的www在线观看| 亚洲人成77777在线视频| av国产精品久久久久影院| 欧美黑人欧美精品刺激| 青春草亚洲视频在线观看| 欧美黑人精品巨大| 熟女少妇亚洲综合色aaa.| 国产精品久久久av美女十八| 爱豆传媒免费全集在线观看| 国产在线观看jvid| 免费不卡黄色视频| 无遮挡黄片免费观看| 精品国产一区二区久久| 丝瓜视频免费看黄片| 亚洲人成电影免费在线| www.自偷自拍.com| 久久这里只有精品19| 久久久久久人人人人人| 欧美日韩福利视频一区二区| www.自偷自拍.com| av片东京热男人的天堂| 人人妻,人人澡人人爽秒播 | 欧美日韩黄片免| 只有这里有精品99| 另类精品久久| 婷婷色综合大香蕉| 国产精品久久久久久精品古装| 久久精品亚洲av国产电影网| 美女视频免费永久观看网站| 别揉我奶头~嗯~啊~动态视频 | 亚洲精品美女久久av网站| 精品一区二区三区av网在线观看 | 国产精品久久久久久精品古装| 亚洲五月婷婷丁香| 黄色 视频免费看| 国产一区二区三区av在线| 在线观看www视频免费| 又大又爽又粗| 天堂俺去俺来也www色官网| 这个男人来自地球电影免费观看| 国产日韩欧美视频二区| 青草久久国产| 中文字幕亚洲精品专区| 免费在线观看日本一区| 在线天堂中文资源库| av国产精品久久久久影院| 看免费成人av毛片| 久久国产亚洲av麻豆专区| 女人久久www免费人成看片| av在线播放精品| 久久鲁丝午夜福利片| 一区二区三区四区激情视频| videosex国产| 丰满少妇做爰视频| 大话2 男鬼变身卡| 精品福利观看| 大香蕉久久成人网| 国产高清视频在线播放一区 | 少妇猛男粗大的猛烈进出视频| av有码第一页| 国语对白做爰xxxⅹ性视频网站| 成人亚洲欧美一区二区av| 女人爽到高潮嗷嗷叫在线视频| 丝袜美腿诱惑在线| 中文乱码字字幕精品一区二区三区| 天堂中文最新版在线下载| av国产久精品久网站免费入址| 最近中文字幕2019免费版| 婷婷色综合www| 性色av乱码一区二区三区2| 超碰成人久久| 久久免费观看电影| 免费看av在线观看网站| 视频区图区小说| 亚洲av在线观看美女高潮| 菩萨蛮人人尽说江南好唐韦庄| 丁香六月天网| 欧美中文综合在线视频| 免费在线观看完整版高清| 我的亚洲天堂| 国产三级黄色录像| 一本综合久久免费| cao死你这个sao货| 丰满饥渴人妻一区二区三| 成人手机av| 午夜视频精品福利| 最新在线观看一区二区三区 | 日本wwww免费看| 国产野战对白在线观看| 一级片'在线观看视频| 考比视频在线观看| 人人妻人人添人人爽欧美一区卜| 欧美日本中文国产一区发布| 亚洲欧美一区二区三区久久| 大型av网站在线播放| 久热爱精品视频在线9| 超色免费av| 丝袜在线中文字幕| 免费不卡黄色视频| 在线观看免费视频网站a站| 国产高清国产精品国产三级| 日韩,欧美,国产一区二区三区| 中文字幕高清在线视频| 国产免费福利视频在线观看| 国产精品一区二区精品视频观看| 国产精品一国产av| 午夜免费成人在线视频| 日韩av免费高清视频| av视频免费观看在线观看| 黄片播放在线免费| 性高湖久久久久久久久免费观看| 欧美性长视频在线观看| 亚洲av国产av综合av卡| 午夜视频精品福利| 国产成人欧美在线观看 | 亚洲人成电影观看| 在线 av 中文字幕| 肉色欧美久久久久久久蜜桃| 91精品国产国语对白视频| 久久亚洲国产成人精品v| 精品久久蜜臀av无| 侵犯人妻中文字幕一二三四区| 91精品三级在线观看| 免费av中文字幕在线| 在线观看一区二区三区激情| 国产精品偷伦视频观看了| 一区二区三区乱码不卡18| 国产一区二区三区综合在线观看| 伦理电影免费视频| 久久精品熟女亚洲av麻豆精品| 99re6热这里在线精品视频| 多毛熟女@视频| 天天躁狠狠躁夜夜躁狠狠躁| 性少妇av在线| 亚洲国产成人一精品久久久| 成人国产av品久久久| 少妇 在线观看| 日本av免费视频播放| 亚洲图色成人| 午夜久久久在线观看| 尾随美女入室| 大片免费播放器 马上看| 最近最新中文字幕大全免费视频 | 国产欧美日韩综合在线一区二区| 在线观看一区二区三区激情| 丰满迷人的少妇在线观看| 性高湖久久久久久久久免费观看| 欧美精品啪啪一区二区三区 | 免费不卡黄色视频| 80岁老熟妇乱子伦牲交| 高清不卡的av网站| 久久精品国产亚洲av高清一级| 一边摸一边做爽爽视频免费| 午夜老司机福利片| 亚洲精品久久成人aⅴ小说| 黑人欧美特级aaaaaa片| 国产av精品麻豆| 九色亚洲精品在线播放| 免费在线观看黄色视频的| tube8黄色片| av天堂久久9| videosex国产| 精品福利观看| 国产熟女午夜一区二区三区| 国产精品久久久av美女十八| 久久综合国产亚洲精品| 国产男女超爽视频在线观看| 美女国产高潮福利片在线看| 久9热在线精品视频| 波多野结衣av一区二区av| 一级毛片电影观看| 成人亚洲精品一区在线观看| 欧美精品亚洲一区二区| 女警被强在线播放| 精品一品国产午夜福利视频| 脱女人内裤的视频| 满18在线观看网站| 老司机靠b影院| 可以免费在线观看a视频的电影网站| 国产91精品成人一区二区三区 | 高潮久久久久久久久久久不卡| 午夜福利在线免费观看网站| 久热爱精品视频在线9| 国产极品粉嫩免费观看在线| h视频一区二区三区| 亚洲色图综合在线观看| 高清黄色对白视频在线免费看| 国产精品久久久久久精品古装| 成人国产av品久久久| 黑人巨大精品欧美一区二区蜜桃| 免费久久久久久久精品成人欧美视频| 国产欧美亚洲国产| 下体分泌物呈黄色| 日韩伦理黄色片| 巨乳人妻的诱惑在线观看| 50天的宝宝边吃奶边哭怎么回事| 亚洲九九香蕉| 精品国产超薄肉色丝袜足j| 日本vs欧美在线观看视频| 一区二区日韩欧美中文字幕| 热re99久久国产66热| 一级,二级,三级黄色视频| 国产欧美日韩一区二区三区在线| 免费久久久久久久精品成人欧美视频| 国产一区二区激情短视频 | 成人午夜精彩视频在线观看| 久久 成人 亚洲| 欧美精品人与动牲交sv欧美| 好男人电影高清在线观看| 午夜久久久在线观看| 美女高潮到喷水免费观看| 丝袜脚勾引网站| 视频区图区小说| 欧美变态另类bdsm刘玥| 国产成人精品久久二区二区免费| 亚洲精品自拍成人| 97在线人人人人妻| 久久综合国产亚洲精品| 久久精品国产亚洲av高清一级| 亚洲国产欧美日韩在线播放| 一边亲一边摸免费视频| 性高湖久久久久久久久免费观看| 一级a爱视频在线免费观看| 亚洲国产av新网站| 1024香蕉在线观看| 丰满少妇做爰视频| 国产色视频综合| 国产又爽黄色视频| 日日爽夜夜爽网站| 日韩av在线免费看完整版不卡| 亚洲国产精品国产精品| 亚洲国产精品一区三区| 肉色欧美久久久久久久蜜桃| 美女视频免费永久观看网站| 国产激情久久老熟女| 久久精品人人爽人人爽视色| 欧美成人午夜精品| 日韩人妻精品一区2区三区| 午夜视频精品福利| 欧美国产精品一级二级三级| 色网站视频免费| 韩国高清视频一区二区三区| 777米奇影视久久| 黄色一级大片看看| 久久国产精品影院| 成人国产av品久久久| tube8黄色片| 久久青草综合色| 国产精品 国内视频| 十八禁高潮呻吟视频| 国产精品一区二区免费欧美 | 男女边吃奶边做爰视频| 久久国产精品大桥未久av| 在线观看www视频免费| 成人三级做爰电影| 视频区欧美日本亚洲| 亚洲欧美激情在线| 大码成人一级视频| 看十八女毛片水多多多| 日韩中文字幕欧美一区二区 | 久久精品国产亚洲av涩爱| 美女国产高潮福利片在线看| 精品一区二区三区四区五区乱码 | 免费久久久久久久精品成人欧美视频| 久久鲁丝午夜福利片| 一级a爱视频在线免费观看| 亚洲免费av在线视频| 欧美日本中文国产一区发布| 亚洲精品一卡2卡三卡4卡5卡 | 亚洲精品国产色婷婷电影| 肉色欧美久久久久久久蜜桃| 中文字幕亚洲精品专区| 搡老乐熟女国产| 大话2 男鬼变身卡| 黑人巨大精品欧美一区二区蜜桃| 黑人欧美特级aaaaaa片| 天天操日日干夜夜撸| 女人高潮潮喷娇喘18禁视频| 我的亚洲天堂| 日韩av不卡免费在线播放| tube8黄色片| 少妇猛男粗大的猛烈进出视频| 只有这里有精品99| 久久精品熟女亚洲av麻豆精品| 亚洲国产最新在线播放| 亚洲一卡2卡3卡4卡5卡精品中文| 成年动漫av网址|