999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

高超聲速航空發動機強預冷技術研究

2015-04-28 02:55:06鄒正平劉火星唐海龍萬敏王洪偉陳小龍陳懋章
航空學報 2015年8期
關鍵詞:發動機

鄒正平,劉火星,唐海龍,萬敏,王洪偉,陳小龍,陳懋章,*

1.北京航空航天大學 能源與動力工程學院 航空發動機氣動熱力國家級重點實驗室,北京 100191 2.先進航空發動機協同創新中心,北京 100191 3.北京航空航天大學 機械工程及自動化學院,北京 100191

伴隨著科技的進步,人類對探索和利用空間越來越重視,因此對高超聲速飛行器的需求凸顯出來。高超聲速飛行器具備即時、高精度打擊能力,比現有各種隱形技術有更好的突防能力,且既具威懾性又具實用性,因而被稱為“常規的戰略武器”,將在未來的空天戰中成為主角。另一方面,現有的各種航天器多以多級火箭為動力,分級推動航天器進入最終的空間軌道。這種方式的缺點是不能重復使用,酬載成本很高。降低有效載荷成本、重復使用一直是航天器動力追求的目標。鑒于高超聲速飛行器所具備的巨大潛在戰略、戰術及經濟價值,美、俄、歐、日等航空強國均已將其作為研究熱點,并給予高度重視和持續的資金支持,以期優先占領此未來技術制高點。美國早已圍繞高超聲速飛行器及其動力裝置,開展了大量研究,如 Hy Fly(Hypersonic Flight Demonstrator)計劃、Hy Tech(Hypersonic Technology)計劃、Hyper-X(eXperimental Hypersonic)計劃、NGLT(Next Generation Launch Technology)計劃和高超聲速轟炸機計劃等,并陸續取得關鍵性的進展[1]。如X-43A已成功進行了7倍聲速和10倍聲速飛行試驗,X37B進行了入軌飛行等。2013年,美國公布了其正在研發的SR-72高超聲速飛行器,飛行速度可達Ma=6,兼具偵察與攻擊雙重職能,計劃于2025年首飛,2030年裝備部隊。俄羅斯也在全力研制自己的高超聲速武器并給予優先發展,以期盡快趕上美國。

動力是高超聲速飛行器最為關鍵的決定性技術。美國等航空強國針對高超聲速飛行器動力進行了大量研究工作,但由于技術復雜,技術成熟度普遍較低,可用于水平起降重復使用高超聲速飛行器的動力技術更是如此,短時間還難以見到工程可用的方案。因此,動力已成為高超聲速飛行器的瓶頸所在。

若要降低航天器發射成本或實現高超聲速飛行器水平起降,首先要考慮航空渦輪發動機,在低馬赫數(0~3)飛行時,它是理想的動力,技術成熟、可靠性很高,其比沖是各類發動機中最高的,但在高馬赫數(>3)時航空渦輪發動機性能急劇惡化且可靠性降低。而沖壓發動機在高馬赫數下性能優越。因而,兼具航空渦輪發動機和沖壓發動機優點的渦輪基沖壓組合發動機(Turbine-Based Combined Cycle,TBCC)成為了高超聲速飛行器動力的主要發展方向之一。但時至今日,TBCC仍難以解決模式轉換工程技術困難和“推力縫隙”、“背死重”以及高超聲速飛行時熱防護困難等問題,嚴重制約了其發展和工程應用。

在高超聲速動力的探索過程中,預冷成為了主攻方向之一。所謂預冷是指將進入發動機的高溫空氣進行冷卻,使溫度降低到航空渦輪發動機能正常工作的溫度。各航空強國都提出過自己的預冷方案,并進行了大量研究,如美國提出的射流預冷卻技術方案,日本和俄羅斯提出的預冷器技術方案,英國提出強預冷器與閉式循環組合的技術方案。在這些方案中最值得關注的是英國的方案。

英國噴氣發動機公司(Reaction Enging Limited,REL)在20世紀80年代提出了強預冷技術方案,30多年的沉寂之后于2012年11月宣布,該方案地面試驗取得成功,實現了強預冷,可將以Ma=5飛行的發動機進口氣流滯止溫度降低到發動機能正常工作的進氣溫度。以這種強預冷方案為基礎,REL公司先后研發了兩種發動機,即彎刀(Scimitar)和軍刀(Sabre,或譯為佩刀)。彎刀發動機采用常規的渦扇發動機技術,加上強預冷系統,用于Ma=5的高超聲速飛機。軍刀發動機從原理方案看,則是彎刀發動機與火箭發動機的組合,用于“云霄塔”航天飛行器Skylon,采用吸氣式發動機將飛行器加速到Ma=5.5,再用火箭進一步推動飛行器,實現單級入軌。可將1 kg有效載荷送入地球低軌道的成本從目前的15 000英鎊/千克降至650英鎊/千克,與常規入軌方案相比,具有水平起降、重復使用、單級入軌、低酬載成本、兩次運行間隔時間短等一系列技術優勢。

國際上對于此強預冷航空發動機技術給予了極高評價,認為是“噴氣推進發明以來的第二次革命”。英國大學與科學部長 Willetts發表評價:“這項技術會徹底改變未來的空中和太空旅行”,英國業內人士預計2025年可最終在工程上實現這項技術。并認為,強預冷航空發動機技術可使整個航空發動機面貌發生根本性變化,具有非常巨大的潛在的技術優勢,可以把它看成是航空發動機領域的顛覆性技術。

本文將綜述國內外高超聲速強預冷航空發動機相關技術的研究與進展情況。

1 國外航空發動機強預冷技術研究

為解決高馬赫數飛行時,發動機進口的高滯止溫度給發動機帶來的種種不利影響,對進入發動機的空氣降溫的預冷技術日益受到重視,并得到了快速的發展。該技術即在常規發動機前增添預冷裝置,預先冷卻壓氣機進口空氣。預冷可能帶來諸多好處:在高超聲速飛行時,降低進氣溫度可擴展高馬赫數飛行范圍,改善各部件工作條件,緩解發動機機體高溫防護問題;降低進氣溫度可提高進氣密度以增大進氣質量流量,進而增大推力;高速滯止熱的適當利用可提高循環熱效率。實現預冷的途徑主要有兩種:一是在壓氣機進口噴入冷卻介質,其中的典型代表如美國MSE技術應用公司提出的射流預冷卻TBCC發動機方案(Mass Injection Pre-compressor Cooling Turbine-Based Combined Cycle,MIPCC-TBCC);二是利用預冷器,如日本的吸氣式渦輪沖壓膨脹循環發動機(Expander Air Turbo Ramjet Engine,ATREX)、俄羅斯的深冷空氣渦輪發動機(Deeply Cooled Air Turborocket,ATRDC)和英國REL公司等方案。

在射流預冷卻MIPCC-TBCC方案中,通過加裝在傳統渦輪發動機壓氣機前部的液體噴射裝置在進氣道內噴射冷卻介質,達到有效降低壓氣機進口空氣溫度的目的,如圖1所示。美國在F100發動機基礎上發展MIPCC并進行了試驗驗證,結果表明海平面推力增加了約1倍,在24.7 km的高度上最大飛行馬赫數可達3.5,可以有效拓寬原有飛行包線[2]。但是冷卻介質的噴入會使得高空高速飛行時,壓氣機進口來流氧含量下降,通常需要在壓氣機后注入氧化劑。同時該方案還存在比沖降低、發動機長度增加、來流總壓損失、進氣總壓/總溫畸變以及所需冷卻介質量大等問題。實際上為解決高溫大氣條件下起飛推力不足的問題,一個已得到實際應用的辦法就是在進氣道前方噴水。

圖1 MIPCC-TBCC發動機方案Fig.1 MIPCC-TBCC engine scheme

使用預冷器預冷的技術效率相對高,但設計難度增大,其中最為關鍵的是輕質、高效的緊湊型預冷器的設計。日本提出的ATREX發動機方案即采用該技術,通過安裝在進氣道后的空氣/氫預冷器將進口空氣冷卻[3-4](見圖2),并在預冷器前噴灑少量的酒精以有效阻止預冷器結冰[5]。地面縮尺模型試驗結果表明:當風扇進口溫度低于160 K,且預冷器總壓恢復系數為0.9時,ATREX循環發動機的推力和比沖比無預冷時分別增加了2倍和1.5倍。基于ATREX發動機方案,日本宇航研究開發機構(Japan Aerospace eX-ploration Agency,JAXA)進一步研制了高速渦輪發動機(S-發動機),并于2014年2月成功進行了在Ma=4條件下的地面試驗[6]。但是該方案發動機推重比較低,且無法保證整個飛行馬赫數范圍內都擁有較高的比沖(其比沖在Ma=1,5時分別為3 500 s和2 500 s)[7],而且用氫直接冷卻,一旦預冷器發生氫脆問題,氫燃料泄漏到主流道,容易發生爆炸等安全問題;此外,它還存在冷卻所需的氫超過燃燒所需的氫,造成氫的浪費、換熱器單位質量換熱量小以及結霜等問題。雖然存在諸多問題,S-發動機的實質性進展證明了預冷卻技術在提升渦輪發動機工作馬赫數方面具有巨大的潛力。

圖2 ATREX發動機方案Fig.2 ATREX engine scheme

俄羅斯在深冷空氣渦輪發動機(ATRDC)研究中同樣對預冷器冷卻技術進行了探索。ATRDC也是采用燃料氫對進口空氣進行深度預冷并驅動渦輪,如圖3所示。該發動機擁有比ATREX更高的增壓比,壓氣機進口溫度為98~112 K,壓比為40時,其在Ma=0~6范圍內的平均比沖可達2 500 s,推重比高達18~22,但是預冷器很重,約占整個發動機質量(不含進氣道)的40%,它同樣存在預冷所需氫的量大于燃燒所需氫的量的問題,導致約1/2的氫在渦輪中膨脹后直接被排掉,造成了嚴重的資源浪費[8]。

圖3 ATRDC發動機方案Fig.3 ATRDC engine scheme

針對傳統TBCC和上述預冷動力方案的不足,在歐洲航天局和英國政府等的支持下,英國和比利時等聯合提出了一種適用于Ma=5的高超聲速航空發動機即“彎刀”(Scimitar)方案[9],這種發動機主要采用強預冷和閉式循環等新技術,與TBCC相比,具有明顯優勢:無論是低馬赫數還是高馬赫數飛行,耗油率皆可下降18%~23%;所有部件在全部飛行馬赫數范圍內一直工作,不存在TBCC渦輪-沖壓模態之間的轉換和再啟動問題,以及渦輪或沖壓不工作時的“死重”問題;沒有TBCC模態切換時推力下降問題;由于強預冷,使高馬赫數時進口溫度大幅下降,大大減輕了材料的高溫防護問題等,且與已有預冷動力方案相比不會大幅度增加發動機質量,技術優勢非常明顯。2012年11月,研制方噴氣發動機公司(REL)宣布地面試驗達到了強預冷的目標(見圖4),即所發展的預冷技術具備在0.01 s瞬間將100 kg/s空氣的溫度由1 000℃降到-150℃的能力。這是一項意義非凡的重大技術突破,將整個改變航空發動機的面貌,有望成為未來最適用的高超聲速動力技術。

圖4 緊湊快速強換熱器試驗裝置Fig.4 Device diagram for compact precooler test

鑒于英國強預冷方案的巨大技術優勢和已取得的重大突破,該方案已作為高超聲速飛機推進技術的研究重點之一列入歐盟遠期先進推進概念和技術計劃(Long-Term Advanced Propulsion Concepts and Technologies,LAPCAT),獲得了歐盟第六框架和第七框架計劃的經費支持[10]。該項目預計將在2019年進行整機測試,英國政府于2013年7月為該項目2014—2016年的研究投資6 000萬英鎊,為實現2025年發動機投入使用的目標,預計總投資80億歐元。該項目目前已吸引歐盟多個國家的關注和直接參與。同樣美國軍方對“軍刀”和“彎刀”發動機技術高度重視。2014年1月,美國空軍實驗室(Air Force Research Laboratory,AFRL)與英國噴氣發動機公司達成了合作研究和發展協議,要求與英方共享關鍵性的緊湊快速強換熱技術,以支持美國的高超聲速飛機研發。同時,AFRL也開展了針對“軍刀”發動機的建模和性能評估工作,并于2015年初給出初步研究結果,認為該新型發動機方案原理先進可行,并且在工程實現上不存在重大技術障礙,他們將與REL繼續保持緊密合作進一步深入研究相關關鍵技術。日本在2014年S-發動機試驗后,也改變原有研究思路,計劃在下一階段的預冷發動機研究過程中,參照英國“彎刀”方案思路,采用閉式換熱介質循環和惰性冷卻介質等相結合的改進措施來進一步發展高超聲速航空發動機技術。

2 緊湊快速強換熱機理及換熱器設計技術研究

強預冷航空發動機優異性能的核心支撐是緊湊快速強換熱技術,該創新性技術的核心是超臨界換熱介質和微小尺度換熱單元的結合。處于超臨界狀態的流體,其黏性系數和擴散系數接近氣態而傳熱系數和密度等接近液態。因此,用超臨界狀態流體作為換熱介質可得到流動損失低、換熱量大的換熱效果。再結合幾何尺度為毫米或亞毫米級的基本換熱單元,可實現換熱器單位質量換熱能力遠超現有換熱技術的強換熱器。但是將超臨界與微小尺度二者相結合導致換熱器內部流動換熱機制更為復雜,必須對相關機理進行研究,弄清楚超高熱流密度情況下微小尺度流動與換熱耦合機制及參數之間的關聯、超大的微通道長徑比對超臨界流動換熱的影響機制、超強溫度梯度引起的流動結構變化及其與換熱的耦合機制等問題。

首要工作是機理等研究手段的發展,包括發展流熱耦合數值模擬程序,使之具備詳細模擬及分析緊湊快速強換熱器內部復雜流動和換熱的能力;發展超臨界介質微小尺度換熱流動實驗方法、流動測試方法并建立專用試驗平臺,解決緊湊快速強換熱器流動換熱機理試驗及性能試驗能力問題;探討該緊湊換熱技術在航空發動機設計中的應用途徑,解決強預冷高超聲速航空發動機總體性能設計問題等。限于篇幅,有關試驗和發動機總體性能等方面的內容,本文不作介紹。

2.1 高精度流熱耦合數值模擬方法及校驗

該數值模擬程序應具備的主要特征是精度高、通用性好、功能豐富、對湍流模型和轉挨模型的調試改進具有可擴展性[11]。

1)采用基于預處理的時間推進法,可求解由多孔材料、純流體以及純固體組成的耦合區域復雜流動和換熱問題。程序分別采用基于原始變量以及守恒變量的預處理方法來求解純流體域流動問題,而求解多孔區域流動與換熱問題只采用基于原始變量的預處理方法。

對于純流體域定常流場的預處理方法,其控制方程組可以寫為

式中:Ω為空間有界閉合區域;S為封閉曲面;W=[ρ ρvxρvθr ρvrρE]T為守恒變量,ρ 為密度,r為半徑,vx、vθ、vr分別為軸向、周向、徑向的速度,E為單位質量流體的總能量;Fc、Fv和Q分別為對流項、擴散項和源項。

經過預處理后的Navier-Stokes方程,用原始變量Wp表示為

式中:Wp=[p vxvθvrT]T,T為溫度,p為壓力;ˉΓ為預處理矩陣[12]。

多孔介質區域宏觀控制方程與純流體Navier-Stokes方程具有相似的方程形式,采用相似的預處理方法來解決時間推進法在求解多孔域宏觀控制方程所遇到的方程剛性問題。根據不同的實際問題,多孔域能量方程采用局部非熱平衡(Local Thermal Equilibrium,LTNE)模型或者局部熱平衡(Local Thermal Nonequilibrium,LTE)模型。

2)添加基于局部變量的γ-Reθt轉挨模型,提高了有轉挨條件下的流熱耦合數值模擬精度。

Menter&Lantry提出的γ-Reθt轉挨模型是基于k-ω剪切應力輸運(SST)湍流模型的四方程模型,這里只給出間歇因子γ輸運方程以及Reθt輸運方程。

守恒形式的間歇因子輸運方程為

式中:Tu為當地湍流度;U為當地速度;θ為動量厚度;s為流線弧長;F(λθ)代表了當地壓力梯度以及湍流度的影響;γ為間歇因子;ui為速度分量;μT為渦黏系數,μL和μ均為動力黏性系數;Pγ和Eγ分別為過渡源項和再層流化源項;λθ為流向壓力梯度。

3)添加保證通量守恒的非匹配流固交界面邊界條件,增強了程序處理復雜流熱耦合問題的能力。分區流熱耦合算法要求流固交界面上溫度和熱流通量連續,程序中采用的是流體域向固體域傳遞邊界熱流量,固體域向流體域傳遞邊界溫度的方法。其中交界面上的溫度以及熱流密度采用了有限單元形函數插值方法進行插值[13],在交界面上將固體域的溫度傳遞給流體域,同理將流體域的熱流密度傳遞給固體域,直到最終收斂。

4)添加了非定常模塊以及并行功能,實現了多孔/流體/固體非定常流熱耦合問題的快速求解。

通過一系列具有代表性的算例,分別對程序的對流項離散格式、各種湍流模型、轉挨模型、預處理方法、熱傳導求解器和流熱耦合方法等進行了詳細的校驗和分析,結果表明各個模塊功能良好,具有較高的模擬精度和可靠性。

第一個校驗算例是對固體肋強化換熱進行校驗,校驗算例來自文獻[14],平行通道上下壁面以交錯的方式各設置2個多孔/固體肋片,不可壓流體以恒定的速度以及溫度進入平行通道,以通道高度和主流進口來流條件定義的Re為100,通道內流動為層流流動。

圖5(a)和圖5(b)分別給出了平行通道采用固體肋時其上下壁面的局部努塞爾數Nu分布,并與H.Y.Li采用Simpler方法得到的結果進行了對比。圖中橫坐標為x方向相對位置,縱坐標為Nu,其中x為軸向當地位置,H為通道軸向長度。從中可以看出,由于固體肋的阻塞作用,使得流體以較高的速度流過肋片與通道之間的區域,故在下壁面設置肋片可以增強上壁面相應位置處的換熱,Nu存在著峰值。本文結果與文獻結果符合得比較好,說明本程序能夠準確地處理較為復雜的不可壓流熱耦合問題。

圖5 采用固體肋時通道上下壁面局部Nu分布(Re=100)Fig.5 Local Nu profiles at upper and lower walls of channel with solid baffles(Re=100)

圖6 可滲透壁面平行通道內層流流動Fig.6 Laminar flow in a plate channel with permeable walls

圖7 不同注入率時軸向速度分布Fig.7 Axial velocity distributions at different coolant injection rates

第二個校驗算例是可滲透壁面層流平行通道流動。圖6給出了可滲透壁面平行通道示意圖,平行通道上下壁面均為可滲透壁面,流體以恒定的縱向速度Vw從下壁面流入通道,又以Vw的速度流出上壁面,通道內流動為層流流動,且流動是不可壓的,當通道內的流動處于充分發展狀態時,通道內流體速度分布為橫向速度U(y)與縱向速度Vw的疊加[15]:

圖7給出了不同流體注入率時通道橫截面上橫向速度分布,并與解析解進行了對比。圖中縱坐標為相對位置,橫坐標為橫向相對速度大小,其中Um為通道截面平均速度。從中可以看出,隨著流體注入率F的增大,通道內的流體橫向速度分布越來越不對稱,橫向速度最大值變大,最大值位置也向上壁面移動。在不同的注入率下數值計算結果(離散點)與解析解(曲線)均符合得非常好,說明本程序計算的可靠性。

此外,還通過Rayleigh-Bénard流動、方腔周期振蕩頂蓋驅動流、多孔介質栓塞流、多孔介質中的Bénard流動等多個典型算例對該程序進行了校驗和分析,驗證結果表明,該程序對微小尺度的多孔/流體/固體多區域流熱耦合問題具有較高的模擬精度,可以滿足緊湊快速強換熱器數值模擬要求[16]。

在發展流熱耦合數值模擬程序的同時,還針對緊湊快速強換熱器的特點,對商用軟件的預測精度進行了校驗,分析湍流模型、邊界條件、浮升力以及流動加速等因素對緊湊快速強換熱器內流動和換熱的影響。圖8給出了用商用軟件計算緊湊快速強換熱器內流動得到的結果。圖中橫坐標為沿流向的幾何相對位置,縱坐標為溫度,其中d為換熱單元內徑。黑色方點為實驗數據[17],各曲線為使用不同湍流模型得到的數值模擬結果。數值模擬結果與實驗結果的對比表明,合理地選取湍流模型能夠得到更為準確的計算結果。

圖8 緊湊快速強換熱器內流動和換熱計算結果Fig.8 Computing results of flow and heat transfer in compact precooler

圖9 多孔介質強化換熱計算模型Fig.9 Computational models for heat transfer enhancement in porous media

2.2 緊湊快速強換熱機理

利用2.1節發展和校驗的數值模擬手段,本文對緊湊快速強換熱器內部流動換熱機理進行了研究,分別從換熱器內、外換熱兩個方面入手,目的在于以最小的流動損失代價得到最大的強化換熱效果。

在內換熱方面,首先對多孔介質微小尺度強化換熱機理及Re、孔隙率、達西數、多孔介質層厚度等參數的影響規律進行了研究。采用多孔材料是實現強化換熱的有效手段,已有研究表明,由于多孔材料可以擁有微米甚至納米級的微觀結構,流體在流過多孔介質時產生強烈混合,從而能夠大大改善換熱效果。限于篇幅,本文只給出了Re、達西數等參數的影響規律,分析表明高滲透率、高熱傳導系數以及高比表面積的多孔材料能更有效地增強換熱性能。雖然多孔材料能夠有效強化換熱,但是由于目前多孔材料在加工過程中的不確定性較大,難以實現關鍵參數的精確控制,距離實際應用還有較大差距。本文進一步對微小尺度圓管作為換熱基本單元情況下,超臨界介質流動換熱、微小尺度流動換熱機理以及二者的耦合作用機制進行了探索,初步得到了微小尺度圓管內部超臨界介質流動換熱的關鍵特征參數影響規律。研究結果表明以超臨界和微小尺度為特征的內換熱可以得到很高的對流換熱系數,完全可以滿足緊湊快速強換熱的要求。

在外換熱方面,主要考慮如何增大單位介質換熱面積和提高外表面對流換熱系數。緊湊快速強換熱器工作在航空發動機環境中,外換熱介質為空氣或燃氣,提高外換熱能力最為直接的方法是減小基本換熱單元的徑向尺寸,實現微小尺度換熱。而提高對流換熱系數則要從流動結構入手,通過控制流動結構的發展,實現以最小的流動損失為代價提高外表面對流換熱系數的目的。

2.2.1 多孔介質微小尺度強化換熱機理

為實現強化內換熱的目的,利用上述發展和校驗的多孔/流體/固體多區域三維流熱耦合數值模擬程序,對交錯肋片平行通道[18]、受限層流沖擊射流[19]以及層流平板發散冷卻[20]這3個多孔介質強化換熱典型應用中的流動換熱機理進行了研究,其模型如圖9所示。圖9(a)中Lt為通道總長,Lc為通道中帶肋部分總長,U為氣流速度,L為肋間距,Li為第一個肋片與進口間長度,t為肋厚,h為肋高,H為通道高度;圖9(b)中L為多孔介質層長度,h為多孔介質層厚度,B為槽寬,H為計算域厚度,vin為氣流進口速度,Tin為氣流進口溫度,Tw為底面溫度;圖9(c)中H為多孔介質層厚度,Vc為冷卻流體速度,Tc為冷卻流體溫度,Uin為主流速度,Tin為主流溫度,δv為邊界層位移厚度。

圖10 不同Re、Da下壁面的局部Nu分布Fig.10 Local Nu profiles at lower wall with different Reynolds numbers and Darcy numbers

圖10 (a)給出了設置有交錯多孔/固體肋片的平行通道內不同雷諾數時下壁面的局部Nu分布,橫坐標為相對位置,縱坐標為局部Nu(Nu=h L/λ,h為對流換熱系數,L為特征長度,λ為導熱系數)。從圖中可以看出,下壁面處的局部Nu分布呈現出周期性分布,由于多孔肋對于流體的阻礙作用,大部分流體會從多孔肋與下壁面的通道通過,會增強該處的對流換熱效果;另一方面,水平流動流體會逐漸改變方向,從而對下壁面有沖擊作用,也會增強該處的對流換熱作用;之后在兩個連續肋片之間,Nu出現極小值。另外,隨著Re的增大,Nu分布曲線整體上移,換熱能力增強。圖10(b)給出了不同Da數(Da=K/L2,K為多孔介質滲透率,L為特征長度)時通道上下壁面的局部Nu分布,橫坐標為相對位置,縱坐標為局部Nu。顯然,Da對通道壁面局部Nu分布總體趨勢影響不大,但對Nu的峰值大小有很大的影響。隨著Da的減小,采用多孔肋的換熱特性與采用固體肋片的相類似。與固體肋片相比,采用多孔肋會使得流體傾斜沖擊壁面的能力減弱,故Nu極大值會減小,同時流體能夠穿透多孔肋會改善肋片之后的對流換熱,Nu極小值峰值也會有所削弱,這說明采用多孔肋會在整體上削弱換熱,但能使得通道壁面的換熱更加均勻。隨著Da的減小,多孔肋對流體的阻礙作用越來越小,Nu峰值進一步減小。

圖11(a)給出了受沖擊壁面上鋪設有多孔介質層的受限層流沖擊射流在不同Re時通道下壁面的Nu分布,圖中橫坐標為相對位置,縱坐標為Nu,ε為孔隙率,B為槽寬,Bi為畢渥數,λs為固體導熱系數,λf為流體導熱系數,Cf為阻力系數。從中可以看出隨著Re的增大,滯止區的換熱性能逐漸增大,之后均沿流向急劇減小。另外,由于多孔介質的存在使得光通道中下壁面區域存在的二次回流區消失,故相應的二次峰值現象也消失,相應Nu分布更加平滑。圖11(b)給出了相應條件下不同Da時下壁面Nu的流向分布,橫坐標為相對位置,縱坐標為Nu。從中可以看出,隨著Da的增大,穿透多孔層沖擊到下壁面的冷卻流體增多,滯止區的Nu隨之增大;當Da小到通道內重新出現二次回流區時,下壁面Nu分布相應的會出現二次峰值現象,但不明顯。

圖12(a)給出了典型平板發散冷卻中不同冷氣注入率時的多孔域/流體域交界面無量綱溫度分布,圖中橫坐標為相對位置,縱坐標為無量綱溫度(θ=(Ti-Tc)/(Tin-Tc),其中Ti為交界面溫度,Tc為冷卻介質進口溫度,Tin為主流進口溫度)。從圖中可以看出,隨著冷氣注入率F的增加,多孔域/流體域交界面的溫度顯著降低,且隨著冷卻流體的不斷積累,越往下游邊界層厚度越大,相應的溫度下降幅度越高。在本文的計算條件下,冷氣注入率不到0.2%時便能使得相對溫度下降到0.5,說明發散冷卻的冷卻效果非常好。圖12(b)給出了其他條件不變,不同主流Re時多孔域/流體域交界面的無量綱溫度分布,圖中橫坐標為相對位置,縱坐標為無量綱溫度。從中可以看出,Re的改變對于層流發散冷卻的冷卻效果還是有比較明顯的影響;在本文的計算條件下,隨著Re的增加,多孔域/流體域交界面上的溫度不斷降低,其冷卻效果不斷減弱。對于層流發散邊界層來說,Re主要改變的是邊界層特性,且隨著Re的增加,邊界層內黏性力的作用減弱,邊界層厚度不斷降低,相應的交界面上溫度和對流換熱系數升高,冷卻效果變差。另外,在本文計算條件下,不同的數量級時Re對冷卻效果的影響也不相同,在Re較小時,層流發散冷卻的冷卻效果對Re的變化更加敏感。

圖11 不同Re、Da時下壁面Nu分布Fig.11 Nu profiles at lower wall with different Reynolds numbers and Darcy numbers

圖12 不同主流Re、冷氣注入率多孔域/流體域交界面無量綱溫度分布Fig.12 Dimensionless temperature profiles at porous/fluid interface at different Reynolds numbers and coolant injection rates

2.2.2 微小尺度圓管內部超臨界介質流動換熱規律

提高壁面對流換熱系數是提高換熱能力的主要手段之一,研究氣動熱力幾何參數對換熱規律的影響,可為緊湊快速強換熱器一維設計參數的選取提供依據。限于篇幅,本文只給出圓管直徑D、進口Re、進口質量流量密度G等參數對壁面換熱系數影響規律的分析。

圖13給出了在保證進口Re相同和進口質量流量密度G相同的情況下,圓管直徑D對圓管壁面對流換熱系數的影響,其中橫坐標為距入口的無量綱管徑D/D0(其中D0為研究中選取的特征直徑),縱坐標為壁面對流換熱系數h,在這兩種情況下均保證對不同管徑D的流體單位質量加熱量相同。從圖中可以看出,兩種情況下圓管壁面的對流換熱系數均隨D的減小而增大,并且存在一臨界直徑。

圖13 無量綱管徑對圓管壁面對流換熱系數的影響Fig.13 Heat transfer coefficient profiles at wall of tube for different dimensionless tube diameters

圖14 圓管壁面對流換熱系數隨進口Re的變化關系Fig.14 Heat transfer coefficient profiles at wall of tube with different Reynolds numbers of inlet

圖14 給出了管徑D相同時,圓管壁面對流換熱系數隨進口Re的變化關系,其中橫坐標為進口Re,縱坐標為壁面對流換熱系數h,保證進口Re不同時流體單位質量加熱量相同。從圖中可以看出,圓管壁面換熱系數隨進口Re的增大而增大。這是由于隨著進口Re的增大,速度邊界層厚度減小,而加熱量相同意味著在相同流向位置處的物性參數相同,即普朗特數Pr相同,此時的溫度邊界層厚度也隨著隨度邊界層的減小而減小,因此壁面對流換熱系數增加。上述研究結果表明,減小管徑D和提高進口Re均可以有效提高圓管壁面換熱系數h。

2.2.3 微小尺度圓管內超臨界介質流動換熱熵產分析

良好的換熱能力無疑是緊湊快速強換熱器重要的性能指標之一,但若從氣動角度考慮,由于管內被加熱的超臨界氦及管外預冷空氣均將最終進入渦輪內部進行做功,因此有效控制換熱過程中造成的管內外工質做功能力的下降,對于航空發動機用換熱器的設計而言是極為重要的,即換熱器的設計追求換熱效率與流動損失等指標的綜合性能最優。而流動和換熱等因素引起的熵產之間關系的分析有助于理解換熱過程中熵產生的原因,可為設計提供參考。本文將熵產Sgen分為兩大部分,一部分為黏性耗散引起的耗散熵產SgenF,另一部分是由于導熱過程存在溫差所引起的加熱熵產SgenH,其表達式為

式中:u為速度;ε為湍流耗散系數;λ為導熱系數;α為熱擴散系數;αT為湍流熱擴散系數;δij為克羅內克符號。式(9)等號右邊兩項分別代表統計平均流場內的黏性耗散熵產與由湍流脈動導致的黏性耗散熵產,式(10)等號右邊兩項分別代表統計平均流場內的換熱熵產與由湍流脈動導致的換熱黏性耗散熵產。

本文結合對不同管徑換熱管內超臨界氦的湍流流動數值模擬,對流場內總熵產及其兩個組成部分的情況進行分析。圖15為相同來流雷諾數下熵產隨管徑的變化趨勢,圖16為相同質量流量密度下熵產隨管徑的變化趨勢。

圖15 相同來流雷諾數下熵產隨無量綱管徑的變化Fig.15 Entropy generation profiles for different dimensionless tube diameters with the same freestream Reynolds number

圖16 相同質量流量密度下熵產隨無量綱管徑的變化Fig.16 Entropy generation profiles for different dimensionless tube diameters with the same mass flux

首先從圖15(b)中可以看出,隨管徑的增加,耗散熵產呈現大幅度減小,造成這一現象的原因在于,在相同的來流雷諾數下,管徑越小,來流速度越大,這兩點共同造成流場徑向速度梯度的增加,從而增強了平均流場內的剪切耗散作用,最終使流動耗散隨管徑的減小而增強。而加熱熵產隨管徑的減小則是由于小管徑內溫度梯度較大所致。雖然在同樣的熱流密度下小管徑內由于換熱系數增強使得平均壁面-流體溫差減小,但由于該減小程度不如管徑的減小程度,而造成溫度梯度隨管徑的增加。導致相同來流雷諾數下熵產隨管徑的減小而增大。對于圖16所示情況,由于來流速度相同,同樣會出現在管徑減小下的速度梯度增加,造成耗散熵產的加劇。而導熱溫差導致的加熱熵產也同樣是管徑減小導致的溫度梯度增大所致。從熵產分析結果可以初步看出,減小管徑D雖然可以增強換熱性能,但相應地也會造成流場各處總熵產的升高,導致流動系統能量品質的降低更為嚴重。所以,設計時管徑的選擇應綜合考慮流動損失與換熱效果。

圖17 單排微小尺度圓柱管束總壓恢復系數和對流換熱系數隨進口Re變化Fig.17 Total pressure recovery and heat transfer coefficient profiles with different Reynolds numbers of inlet at a single row of micro-size tubes

圖18 多排微小尺度圓柱管束總壓恢復系數和對流換熱系數隨進口Re變化Fig.18 Total pressure recovery and heat transfer coefficient profiles with different Reynolds numbers of inlet at micro-size tube matrix

2.2.4 微小尺度圓柱擾流換熱機理

采用緊湊排列圓柱管束來實現強化外換熱的同時會增大流動損失。因此,需要通過對緊湊排列圓柱管束中流動換熱參數影響規律的研究,得到綜合考慮流動損失和換熱性能后的最佳參數范圍,以便為緊湊快速強換熱器設計提供參考。圖17給出了單排微小尺度圓柱管束中的流動換熱規律,其中橫坐標為進口Re,圖17(a)和圖17(b)的縱坐標分別為總壓恢復系數ζ和對流換熱系數h。從圖中可以看出,隨著進口Re的增加,單排和多排管束的總壓恢復系數以及壁面平均對流換熱系數總體變化趨勢相同,即總壓恢復系數逐漸降低,對流換熱系數逐漸升高。在單排管束時,總壓恢復系數的變化趨勢基本呈線性規律,而在多排管束情況(見圖18)下,隨著Re的增加,總壓恢復系數的下降趨勢逐漸加劇。在單排管束和多排管束情況下,壁面平均對流換熱系數隨Re增加而增大,但增大趨勢逐漸變緩。這就意味著在選取進口Re時,應根據設計要求的流動和換熱性能選擇適中的范圍。

圖19和圖20分別給出了單排圓管在4個不同時刻的展向渦量與耗散熵產云圖及時均云圖,圖中Vorticity Z表示展向渦量,其中Z向為沿管長方向。由圖可知,圓柱后渦對的強剪切引起了很強的耗散熵產,并遷移到下游較遠的區域。相對而言,圓管外流動中換熱引起的熵產主要集中在圓管壁面附近(見圖21),其影響范圍相對集中,因此從減小熵產的角度來看,應主要從減小耗散引起的熵產入手。在本文研究的情況下,管束中耗散引起的熵產主要包括兩個方面,即脈動耗散引起的熵產和平均流耗散引起的熵產。

圖19 單排圓管時均和瞬時展向渦量分布云圖Fig.19 Time-averaged and transient vorticity Z distributions of a single row of tubes

圖20 單排圓管時均和瞬時耗散熵產分布云圖Fig.20 Time-averaged and transient entropy generation due to energy dissipation for a single row of tubes

圖21 單排圓管時均和瞬時換熱熵產分布云圖Fig.21 Time-averaged and transient entropy generation due to heat transfer for a single row of tubes

圖22和圖23分別給出了多排圓管情況下展向渦量和耗散熵產的非定常時均云圖。與單排管束的時均解對比可以看出,在多排管束中,由于前一排管束尾跡對下游管束的作用,使得下游管束前緣處的耗散熵產明顯增大。同時,由于管束排列比較緊湊,管束尾跡受下游圓管固壁的限制,被擠壓在較窄的范圍內,并在此范圍內產生了較大的熵產。

圖22 多排圓管情況下時均展向渦量分布Fig.22 Time-averaged vorticity Z distributions for tube matrix

圖23 多排圓管情況下時均耗散熵產分布Fig.23 Time-averaged entropy generation due to energy dissipation distributions for tube matrix

2.3 緊湊快速強換熱器設計方法發展

由于傳統的換熱器設計方法已不再適用于緊湊快速強換熱器設計,本文利用前述低維模型,發展了緊湊快速強換熱器的確定性設計方法,并利用該方法完成了用于高超聲速強預冷航空發動機不同位置的多種不同結構形式的緊湊快速強換熱器確定性設計方案。其中,毛細管式換熱器能將高溫來流在流經換熱器的瞬間冷卻至發動機可以正常工作的溫度,同時換熱器單位重量的換熱量超過100 k W/kg。圖24為毛細管式換熱器示意圖。

為保證設計方案滿足設計要求,采用數值模擬手段對設計結果進行了檢驗。利用毛細管換熱器的簡化模型進行數值模擬,模型中空氣流經10排毛細管。為減小計算網格量,每排以一根毛細管進行模擬,保證兩種工質進口流速與設計值相同,網格數總計780萬,網格示意圖如圖25所示。

圖24 毛細管式緊湊換熱器Fig.24 Micro-size tube compact precooler

圖25 毛細管換熱器數值模擬網格示意圖Fig.25 Meshes used in simulations of micro-size tube compact precooler

冷卻后空氣出口溫度的數值模擬結果與設計值對比如圖26所示,圖中橫坐標為沿毛細管高度方向的相對位置,縱坐標為空氣出口溫度。結果證明換熱器的設計方案初步滿足設計要求。

圖26 冷卻后空氣出口溫度數值模擬結果與設計值對比Fig.26 Comparison of cooled air temperature at outlet between computing result and design value

3 緊湊快速強換熱器先進加工及檢測技術

圖27 不同溫度下GH4169晶粒尺寸效應Fig.27 Grain size effect of GH4169 at different temperature

3.1 薄壁毛細管制造及檢測技術

薄壁毛細管的制造是緊湊快速換熱器的關鍵。其壁厚和均勻度嚴重影響到后續的換熱能力、換熱效率和結構可靠性。此外,薄壁毛細管后續的定位裝配過程要求兩端端口具有很高的同心度,以保證釬料均勻分布在待焊接頭處。這些都對薄壁毛細管的制造提出了非常苛刻的要求,為了成功制造薄壁毛細管,必須深入掌握其流動變形機理,包括微尺度變形中的晶粒尺寸效應、幾何尺寸效應關系,以及溫度對其尺寸效應的影響等問題。

圖27分別給出了不同溫度下應力隨晶粒尺寸的變化趨勢,以分析溫度對于GH4169晶粒尺寸效應的影響,圖中橫坐標為晶粒尺寸的倒數,縱坐標為應力。圖例中,25、200、400分別代表變形溫度,0.002等數字代表真實應變大小。

從圖中可以看出,在不同溫度、不同應變下的晶粒尺寸效應曲線中,都存在分段現象,且分段點隨著應變的增加而變大;常溫下存在兩個分段點,第一個分段點隨著應變增加變化趨勢明顯,第二個分段點和應變沒有太大關系,這表明當晶粒小到一定程度后,晶粒對應力的影響可以忽略,應力主要受應變的影響,即晶粒尺寸效應現象不明顯;當在200℃ 的時候,分段點呈微微下凹的趨勢,在同一應變下的線性趨勢最好,應力水平主要受不同應變的影響較大;在400℃下,分段點隨應變增加而增大現象明顯。

將不同溫度下屈服時的晶粒尺寸效應擬合線斜率放在一起,如圖28所示,圖中橫坐標為晶粒尺寸,縱坐標為屈服應力。可以進一步看出明顯的分段現象,另外,存在一臨界晶粒尺寸,晶粒大于該尺寸時,斜率隨溫度下降明顯,且溫度越高下降速度越慢;晶粒小于該尺寸時,斜率隨溫度上升而下降,但下降速度基本不變。這就要求在選擇成形工藝路線時,需要考慮成形溫度、臨界晶粒尺寸等因素對于成形力的影響。

圖29給出了試制的毛細管樣品掃描電子顯微鏡的照片,通過SEM照片,可以看到其內壁粗糙度滿足要求,而對于使用中影響最大的成形管壁內部孔洞缺陷,在隨機選擇的多個樣品的SEM照片中均未發現成形孔洞缺陷,表明薄壁毛細管在拉拔過程中不存在成形孔洞缺陷問題。

表1則給出了3種不同熱處理狀態(軟態、半硬態和硬態)的薄壁毛細管的單向拉伸力學性能檢測結果,結果表明不同狀態的抗拉強度不同,半硬態與硬態的薄壁毛細管抗拉強度相似,而軟態薄壁毛細管抗拉強度相對較低。圖30給出了半硬態成形下的單拉曲線圖,圖中橫坐標為位移,縱坐標為拉力。

圖28 不同溫度下屈服時的晶粒尺寸效應擬合線斜率Fig.28 Different slopes fitted at various temperatures for grain size effect of yield stress

圖29 毛細管壁厚電鏡照片Fig.29 SEM photograph for wall thickness of micro-size tube

表1 不同成型方法毛細管力學性能Table 1 Mechanical properties of micro-size tube in different forming processes

圖30 半硬態單拉試驗結果Fig.30 Results of uniaxial-tension test in middle-hard state

3.2 薄壁毛細管的組合焊接

薄壁毛細管的組合焊接,直接關系到最終產品的質量。在一般釬料的制備方法上,進行適合工況的釬料研發,掌握釬料成分對于焊接強度的作用規律;找出釬焊接頭的組織、性能隨溫度、釬縫間隙等焊接參數的演變規律;探明不同焊接熱循環過程中接頭內部殘余內應力的分布規律;組合釬焊還要求掌握密集毛細管的定位裝配方法,探究多個焊接接頭間熱力學影響規律;最后還需制定大量薄壁毛細管釬焊接頭焊接質量表征方法,合理評價焊接質量。

本文研究中對現有的BNi-2釬料進行改進,并通過釬焊時間的合理控制以保障釬焊質量,工藝試驗表明:通常情況下,釬縫越大,越容易在釬縫中央部位形成連續金屬間化合物的脆性相,導致釬縫的脆性增加;這是由于釬料向母材內部擴散不均勻導致的,本研究中通過試驗得到合理的釬焊間隙。為了改善釬縫組織和提高接頭韌性,焊后進一步在合適溫度下保溫一定時間。最終通過多次試驗確定合理的工藝參數。圖31給出了兩種不同形式焊接試驗件的照片,圖32則給出了焊接處的SEM照片,表明焊接接頭處的釬料溶解均勻,沒有凹凸不平等表面缺陷。

圖31 毛細管焊接試驗件Fig.31 A test sample of welded micro-size tube

圖32 焊接部位的釬料填充情況Fig.32 Filling result of solder at the welded position

4 結 論

國內外的研究表明,強預冷高超聲速航空發動機技術是一項具有非常巨大的潛在技術優勢和前瞻性的共性技術,特別是REL公司的地面試驗已實現了強預冷的技術指標,值得引起關注。目前國內也已開展了強預冷航空發動機相關研究,并在微尺度流熱耦合換熱機理、流熱耦合數值模擬、緊湊強換熱器設計制造等方面均取得一定進展,為下一步的研究打下了良好的理論和技術基礎。

[1] Tang M,Hamilton B A,Chase R L.The quest for hypersonic flight with air-breathing propulsion,AIAA-2008-2546[R].Reston:AIAA,2008.

[2] Balepin V,Engers R,Terry S.MIPCC technology development,AIAA-2003-6929[R].Reston:AIAA,2003.

[3] Miyagi H,Miyagawa T,Monji T,et al.Combined cycle engine research in Japanese HYPR project,AIAA-1995-2751[R].Reston:AIAA,1995.

[4] Tanatsugu N,Sato T,Balepin V.Development study on ATREX engine[J].Acta Astronautica,1997,41(2-8):851-862.

[5] Harada K,Tanatsugu N,Sato T.Development study of a precooler for the air-turboramjet expander-cycle engine[J].Journal of Propulsion and Power,2001,17(6):1233-1238.

[6] Zhang D B.Analysis of Japanese precooling high-speed turbine engines[J].Aeronautical Information,2014,1599(4):1-4(in Chinese).張東寶.日本預冷卻高速渦輪發動機研究取得重要進展[J].航空情報,2014,1599(4):1-4.

[7] Sato T,Kobayashi H,Tanatsugu N.Development study of the precooler of ATREX engine,AIAA-2003-6985[R].Reston:AIAA,2003.

[8] Vladimir B.High speed propulsion cycles,RTO-AVTVKI LS CSP-07-5052[R].Rhode Saint Genese:VKI,2007.

[9] Steelant J.Sustained hypersonic flight in Europe:technology drivers for LAPCAT II,AIAA-2009-7240[R].Reston:AIAA,2009.

[10] Steelant J.Sustained hypersonic flight in Europe:first technology achievements within LAPCAT II,AIAA-2011-2243[R].Reston:AIAA,2011.

[11] Zhang H J,Zou Z P,Li Y,et al.Preconditioned densitybased algorithm for conjugate porous/fluid/solid domains[J].Numerical Heat Transfer Part A,2011,60(2):129-153.

[12] Weiss J M,Smith W A.Preconditioning applied to variable and constant density flows[J].AIAA Journal,1995,33(11):2050-2057.

[13] Zhang B,Han J L.Multi-field coupled computing platform and thermal transfer of hypersonic thermal protection structures[J].Acta Aeronautica et Astronautica Sinica,2011,32(3):400-409(in Chinese).張兵,韓景龍.多場耦合計算平臺與高超聲速熱防護結構傳熱問題研究[J].航空學報,2011,32(3):400-409.

[14] Li H Y,Leong K C,Jin L W,et al.Analysis of fluid flow and heat transfer in a channel with staggered porous blocks[J].International Journal of Thermal Sciences,2010,49(6):950-962.

[15] Nikitin N V,Pavellev A A.Turbulent flow in a channel with permeable walls,direct numerical simulation and results of three-parameter model[J].Fluid Dynamics,1998,33(6):826-832.

[16] Zhang H J.Investigation of numerical conjugate heat transfer method and coupling mechanism for hybrid porous/fluid/solid domains[D].Beijing:Beihang University,2013(in Chinese).張紅軍.多孔/流體/固體多區域流熱耦合數值模擬方法以及耦合機制研究[D].北京:北京航空航天大學,2013.

[17] Li Y.A 3-D conjugate heat transfer solver and methodology research[D].Beijing:Beihang University,2011(in Chinese).李宇.三維流/熱耦合數值模擬程序的發展及方法研究[D].北京:北京航空航天大學,2011.

[18] Zhang H J,Zou Z P,Song S H,et al.Investigations of heat transfer enhancement in a channel with staggered porous ribs by the preconditioned density-based algorithm[J].Numerical Heat Transfer Part A,2015,67(2):1370-1385.

[19] Zhang H J,Zou Z P.Investigation of a confined laminar impinging jet on a plate with a porous layer using the preconditioned density-based algorithm[J].Numerical Heat Transfer Part A,2012,61(4):241-267.

[20] Zhang H J,Zou Z P.Numerical investigation of laminarplate transpiration cooling by preconditioned density-based algorithm[J].Numerical Heat Transfer Part A,2012,62(10):761-779.

猜你喜歡
發動機
元征X-431實測:奔馳發動機編程
2015款寶馬525Li行駛中發動機熄火
2012年奔馳S600發動機故障燈偶爾點亮
發動機空中起動包線擴展試飛組織與實施
奔馳E200車發動機故障燈常亮
奔馳E260冷車時發動機抖動
新一代MTU2000發動機系列
2013年車用發動機排放控制回顧(下)
VM Motori公司新型R750發動機系列
發動機的怠速停止技術i-stop
主站蜘蛛池模板: 啪啪永久免费av| 国产区在线观看视频| 2021天堂在线亚洲精品专区| 亚洲色图欧美| 日本福利视频网站| 91最新精品视频发布页| 十八禁美女裸体网站| 亚洲美女AV免费一区| 国产亚洲高清在线精品99| 国产成人毛片| 日本亚洲欧美在线| 都市激情亚洲综合久久| 动漫精品中文字幕无码| 成人字幕网视频在线观看| 亚洲水蜜桃久久综合网站| 91亚瑟视频| 成人福利在线免费观看| 麻豆精品国产自产在线| 国产成人AV大片大片在线播放 | 伊人福利视频| 亚洲天堂.com| 欧美国产日本高清不卡| 天天干伊人| 久久久久无码国产精品不卡| 亚洲欧美精品一中文字幕| 国产特级毛片aaaaaa| 精品99在线观看| 日本尹人综合香蕉在线观看 | 亚洲第一成年网| 国产一区二区三区夜色| 波多野结衣的av一区二区三区| 欧美国产日韩在线| 国产幂在线无码精品| 国产乱人免费视频| 最新加勒比隔壁人妻| 精品少妇人妻一区二区| 在线观看国产黄色| 国产福利影院在线观看| 日韩第八页| 久久女人网| 亚洲三级影院| 人妻出轨无码中文一区二区| 亚洲一区第一页| 亚洲天堂网在线观看视频| 韩国v欧美v亚洲v日本v| 精品日韩亚洲欧美高清a| 妇女自拍偷自拍亚洲精品| 99久久国产精品无码| 亚洲视频在线青青| 中文成人无码国产亚洲| 国产人前露出系列视频| 国产后式a一视频| 久草中文网| 综合网久久| 欧美成人精品一级在线观看| 57pao国产成视频免费播放| 亚洲一级无毛片无码在线免费视频| 日韩精品一区二区深田咏美| 久久九九热视频| 国产乱子伦视频在线播放| 国产精品页| 精品少妇人妻av无码久久| 毛片基地视频| 精品91自产拍在线| 国产高清毛片| av大片在线无码免费| 91精品国产91久久久久久三级| 国产成人狂喷潮在线观看2345| 免费国产好深啊好涨好硬视频| 亚洲区一区| 国产免费羞羞视频| 亚洲成人黄色在线观看| 亚洲高清无码久久久| 视频一区视频二区中文精品| 无码中文字幕精品推荐| 中文字幕亚洲第一| 国产97视频在线观看| 国产亚洲视频免费播放| 国产精品久久久久久久久| 91久久精品日日躁夜夜躁欧美| 国产成人亚洲无码淙合青草| 欧美国产中文|