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

受限空間欠采樣條件下航空發動機關鍵截面流場重構方法研究進展

2023-11-27 03:33:14徐笳森鄭培英婁方遠
航空發動機 2023年5期
關鍵詞:發動機

徐笳森 ,鄭培英 ,張 軻 ,婁方遠

(1.南京航空航天大學航天學院,南京 211106;2.中國航發沈陽發動機研究所,沈陽 110015;3.清華大學航空發動機研究院,北京 100084)

0 引言

航空發動機及燃氣輪機的性能評估需要測量關鍵截面的流體熱力參數,針對航空發動機內部復雜、非定常、強不均勻流場,探針的測點布局及基于離散測點的數據分析方法對發動機性能評估有較大影響[1]。當前國內外針對探針沿周向和徑向布局(包括探針數目,位置等)缺乏系統的理論指導,還沒有一套通用的周向探針優化方法;針對離散探針數據的處理方法及誤差機理尚未明確;基于離散數據對發動機關鍵截面流場實現重構的研究還處在初期階段。

王振華等[1]強調探針布局對航空發動機及燃氣輪機性能評估非常重要,然而當前國內外針對探針沿周向布局(包括探針數目,位置等)依然廣泛采用基于工程經驗的探針布局方案。在航空發動機的研究工作中,Saravanamuttoo[2]和Pianko 等[3]分別采用沿周向等距、等節距的探針布局,對發動機部件性能以及氣動性能進行了探究。此外,針對離散探針數據的處理方法,當前國內外通常基于均值計算。然而,Livesey等[4]、Pianko 等[3]以及Greitzer 等[5]的研究表明,面積平均、質量平均、總焓平均都不能完全代表所有非均勻流場的氣動參數或熱力參數;其后,Prasad[6]和Cumpsty 等[7]的研究表明,功平均以及推力平均同樣無法完全代表所有非均勻流場參數。綜上所述,基于傳統的探針布局方案及均值處理方法,無法精確評估航空發動機關鍵截面0 維流場參數。而通過李紅麗等[8]、Stummann 等[9]以 及Chilla 等[10-11]等 諸 多 學 者 的 研 究 發現,關鍵截面流場參數的模糊估計對航空發動機整機及部件性能評估會造成較大誤差;Methel 等[12]、白磊等[13]、Seshadri 等[14-15]在進行優化探針位置的研究中,發現提高關鍵截面流場參數的采樣精度可使整機及部件性能評估的準確性有較大提升。例如近期由英國劍橋大學惠特實驗室開展的試驗研究[10]顯示:當采用傳統周向均布的探針布局方案及面積均值計算方法評估多級壓氣機等熵效率時,會造成高達2.8 個點的不確定性和2.8個點的誤差。

為解決這一難題,針對航空發動機內強3 維、強不均勻流場,本文采用公開試驗數據與高精度仿真相結合的方式系統研究航空發動機關鍵截面復雜流場關鍵特征及演變規律;在此基礎上通過解析航空發動機內復雜流場的主導特征,探索在受限空間欠采樣條件下航空發動機內流復雜流場重構技術及其應用。

1 發動機關鍵截面典型流場特征

在航空發動機及燃氣輪機中,由于受進氣畸變、多級壓氣機級間相互作用、燃燒室熱斑的影響,使發動機進口、壓氣機出口、燃燒室出口、渦輪出口截面的流場具有非定常、強不均勻特性,沿著周向、徑向變化。為探明航空發動機及燃氣輪機內流動特征,首先基于公開文獻,開展系統的管流(進排氣系統[16])與環流(壓氣機[17-19]、燃燒室[20]、渦輪[21])氣動特征研究,在此基礎上,采用高保真度仿真方法對進排氣系統、多級軸流壓氣機、燃燒室、渦輪開展仿真分析。

1.1 典型進氣道出口截面流場特征

對于大涵道比民用航空發動機的進氣系統而言,在進氣道或短艙長度小、無進氣畸變的情況下,航空發動機進口流場不均勻度較低。對于軍用發動機而言,常采用S 彎進氣道,進氣道較長,在發動機進口截面會產生總壓及預旋畸變,沿周向具有同軸系固有頻率相同的分布特征。在此基礎上,通過高保真度仿真開展S 型進氣道內流場特征研究。仿真分析結果表明,在S 型進氣道靠近出口的下壁面存在逆壓梯度,原本吸入的附面層加上壁面的粘性作用使得越靠近壁面的氣流速度越低。在逆壓力梯度與下壁面外擴曲線的雙重作用下,流體的沖擊力無法抵抗逆壓力,將會在下壁面某一處產生回流,回流的流體在分離點處與氣流的沖擊力達到平衡,又會向壁面法向方向卷起而再次順流,因而形成大片分離區。靠近上壁面的氣流在上下壁面的壓差作用下沿著側壁偏轉,而靠近下壁面的回流區為上壁面發展的二次流提供了空間,二次流又會反過來抬起回流區,二次互相加強。所以在進氣道出口截面表現為帶有“對渦”形式的總壓畸變與旋流畸變,沿流向進氣道總壓、流線分布圖譜的發展歷程如圖1所示。

圖1 沿流向進氣道總壓、流線分布圖譜的發展歷程

1.2 典型多級軸流壓氣機出口截面流場特征

對于主要由多級軸流壓氣機組成的航空發動機壓縮系統而言,壓氣機內流動具有非定常、強3 維、強不均勻特征。公開研究結果[11-12]表明,決定壓氣機周向流場均勻特征的因素主要為上下游靜子尾跡的疊加及相互作用。此外,壓氣機內流道-流道間的流場具有非周期性,因此基于離散探針數據及周期性流場假設的傳統壓氣機性能評估方法具有較大誤差。近期,Chilla 等[11]針對RR 公司8 級核心機壓氣機內部流場進行了研究,50%葉高截面總壓場分布如圖2 所示。為了精確模擬多級軸流壓氣機級與級之間的相互作用,采用半周非定常RANS仿真分析方法。由于級間相互作用包括上下游靜子尾跡疊加及勢場的相互作用,多級軸流壓氣機末級出口截面的總壓場沿著周向是強不均勻的。盡管如此,壓氣機流場沿周向不均勻分布特征卻有據可循,而且通常由少數幾個波數主導。其中,壓氣機前部和中部的主導波數主要是上游靜子葉片數,而壓氣機后部的主要波數受到出口支板數的影響。以第6 級靜子S6 出口50%葉高截面總壓場分布為例(如圖3 所示),總壓力場的周向不均勻分布特征主要由第5級靜子S5,第6級靜子S6以及之間的相互作用主導。主導波數包括S5、S6,以及S6-S5。

圖2 RR公司8級軸流壓氣機50%葉高截面總壓場分布

圖3 RR公司多級壓氣機S6出口截面50%葉高總壓分布

此外,基于歷史試驗數據研究了3 級軸流壓氣機中第1 級和中間級出口總壓場的周向不均勻度[17]。例如,第1、2 級靜子S1、S2 出口截面88%葉高處不同流道的總壓沿著周向的分布以及相對應的空間頻譜分別如圖4、5 所示。顯而易見的是,S1 和S2 出口總壓場沿周向存在顯著的不均勻性。不同流道出口的總壓分布不盡相同。因此無論是采用單一流道試驗測量結果還是基于單一流道及周期性邊界條件假設的仿真結果,都會造成較大的壓氣機性能評估誤差。整體而言,由于多級軸流壓氣機級間相互作用,尤其是靜子尾跡的相互疊加及同勢場的相互作用,不同流道出口的總壓分布不盡相同且分散度較大。然而在試驗測量中,通常是基于不同流道內流動一致的假設,沿不同流道不同節距位置布置探針,因此會造成較大的試驗性能評估誤差,以本研究中的3 級軸流壓氣機為例,由于流場不均勻及離散探針數據造成的壓氣機等熵效率評估誤差可達到2.3 個百分點,出口總壓不均勻度對壓氣機等熵效率評估影響如圖6所示。

圖4 第1級靜子S1下游88%葉高截面總壓場分布

圖5 第2級靜子S2下游88%葉高截面總壓場分布

圖6 出口總壓不均勻度對壓氣機等熵效率評估影響

1.3 典型燃燒室出口截面流場特征

燃燒室出口溫度分布的不均勻性是影響航空發動機渦輪工作狀態的重要參數之一,對渦輪的效率和壽命有著重要的影響。由于受燃燒噴嘴的離散特性和摻混氣體的影響,燃燒室出口平面的溫度測量顯示出大的徑向和圓周溫度變化,流動具有強不均勻性。通常認為:

(1)燃燒室徑向溫度不均勻由襯里冷卻氣流主導。

(2)燃燒室周向不均勻由噴嘴及冷卻氣流及主流摻混主導。

通常采用溫度分布因子(Temperature Distribution Factor,TDF)量化燃燒室出口處的溫度不均勻性程度。2 個廣泛采用的參數包括徑向TDF(RTDF)和整體TDF(OTDF)。RTDF是測量周向平均溫度場不均勻性的參數,而OTDF是測量最熱斑與平均溫度的差異的模式因子。定義為

式中:T3為燃燒室進口溫度;T4為燃燒室出口溫度,上標area代表區域平均,cir代表周向平均。

Qinetiq 公司測試的軍用航空發動機全環形燃燒室出口溫度分布如圖7 所示[20]。該燃燒室出口平均溫度為2072 K,且沿著徑向、周向具有強不均勻、非周期特征。基于公開文獻[21-23]數據的典型燃燒室出口徑向溫度分布如圖8 所示。從圖中可見,燃燒室徑向溫度分布同燃燒室設計及內襯冷氣流量高度相關。不同燃燒室的出口徑向溫度分布因子形狀分散度較大[24]。

圖7 軍用發動機全環形燃燒室出口溫度分布

圖8 航空發動機燃燒室典型出口徑向溫度分布

1.4 典型渦輪出口截面流場特征

對于航空發動機及燃氣輪機中的渦輪部件而言,由于燃燒室噴嘴以及冷卻射流的離散性,渦輪入口或燃燒器出口溫度場本質上有著強烈的不均勻性。渦輪內流場不均勻主要由以下3種因素造成:

(1)燃燒室出口不均勻溫度場,即熱斑在渦輪中的遷移。

(2)渦輪葉排之間的相互作用,包括靜子尾跡的疊加及相互作用。

(3)高壓渦輪冷卻氣流同主流路的摻混。

本文基于高保真仿真方法對1 個典型雙級高壓渦輪開展了半周非定常仿真分析[25]。渦輪為GE 公司設計的E3-雙級高壓渦輪。研究內容包括:

(1)燃燒室出口溫度場及熱斑對渦輪內流場特征的影響。

(2)渦輪內葉排相互作用對渦輪內流場特征的影響。

(3)高壓渦輪冷卻氣流摻混對渦輪內流場特征的影響。

數值模擬所采用的湍流模型為k-ε模型。渦輪的流體域離散化后,在所有端壁/葉片表面上都實現了10~50范圍內的y+值。高壓渦輪計算域、網格劃分及冷卻氣流孔位置如圖9 所示。渦輪冷卻采用源項進行模擬。

圖9 高壓渦輪計算域、網格劃分及冷卻氣流孔位置

高壓渦輪進口完美周期性熱斑總溫分布及真實整機環境總溫分布如圖10 所示。通過對比具有完美進口周期溫度場(圖10(a))及具有非周期性進口溫度場(圖10(b))的2個仿真結果展開分析。高壓渦輪進口完美周期性熱斑總溫分布及真實整機環境總溫分布如圖11所示。從圖中可見,在進口為完美周期性熱斑情況下,在高壓渦輪入口觀察到的交替熵的重復模式對應于周期性入口溫度分布;相比之下,當燃燒室出口溫度具有非周期特征時,高壓渦輪入口處的熵分布周期性較低。隨著熱斑在高壓渦輪葉片行中遷移,熵模式在這2 種情況下都變得更加復雜,這主要是由于葉片行相互作用。對于這2 種情況,在高壓渦輪出口處仍然可以看到熱斑的足跡,這與公開文獻中報告的情況一致。

腰椎管狹窄癥作為脊柱外科臨床工作的重要組成部分,往往表現出行走間歇性跛行,其診斷往往需要借助于臨床影像學(CT、MRI等)檢查。Amundsen T等[1]通過臨床癥狀與影像學相關指標分析發現,腰椎管狹窄的程度與患者實際臨床癥狀并沒有顯著的相關關系。相對于腰椎管狹窄癥患者長期存有的腰部疼痛,基于臨床癥狀診斷顯得意義不足。

圖10 高壓渦輪進口完美周期性熱斑總溫分布及真實整機環境總溫分布

圖11 高壓渦輪進口完美周期性熱斑總溫分布及真實整機環境總溫分布

此外,進口完美周期性熱斑情況下50%葉高截面總溫按空間和頻域的分布如圖12 所示。從圖12(a)中可見,在高壓渦輪入口處具有完美周期性熱斑的基線情況下①~⑤的中跨時間平均總溫度分布:①位于高壓渦輪入口處,沿圓周方向的溫度變化為61.5 K,相當于沿周向平均總溫度的7.9%。熱斑在第1 排葉片中幾乎沒有衰減,只有約12.3°的輕微相移。與定子排相比,熱斑在轉子葉排遷移過程中衰減得更多。

圖12 進口完美周期性熱斑情況下50%葉高截面總溫按空間和頻域的分布

總體而言,熱斑在接近高壓渦輪出口時明顯衰減,環向溫度變化的絕對值降至39.6 K,占當地平均總溫度的8.7%。此外,當入口熱斑穿過高壓渦輪時,圓周溫度不均勻性的驅動機制也在演變。在高壓渦輪入口(截面①),中跨處的周向不均勻性主要由波數為18 的熱斑主導,從而形成極佳的正弦形狀。然而,在高壓渦輪出口處,更多影響因素包括高壓渦輪葉排間相互作用對圓周溫度不均勻性作出了貢獻。

為了研究圓周溫度不均勻性的驅動機制,進行了空間快速傅立葉變換(Fast Fourier Transform,FFT),進口非周期性熱斑情況下50%葉高截面總溫按空間和頻域的分布如圖13所示。從圖13(b)中可見,正如預期的那樣,熱斑(波數Wn=18)壓倒性地控制了高壓渦輪入口溫度圓周不均勻性。在S1 的出口(截面②),可以觀察到S1 的影響,由1 個波數為46 的小峰表示。在1 級轉子R1 的出口(截面③),也可以觀察到S2 的勢場,由Wn=48 的小峰表示。值得注意的是,頻譜圖沒有任何與轉子葉片排相關的波數。

圖13 進口非周期性熱斑情況下50%葉高截面總溫按空間和頻域的分布

這是因為流場是時間平均結果,所以通過平均消除了轉子小波數的存在,模擬了使用穩態探頭測量總溫度的場景。

從圖13 中可見在高壓渦輪入口處具有非周期性熱斑的Qinetiq 情況下截面①~⑤的時間平均總溫度分布。燃燒器與燃燒器之間不對稱的存在不會影響熱斑遷移的總體趨勢,例如熱斑的幅度在轉子上比在葉片排或定子上衰減得更多。然而,燃燒器到燃燒器的不對稱性會在空間域中產生更分散的能量分布(圖13(b))。盡管熱斑分量(Wn=12)是高壓渦輪入口溫度圓周不均勻性的主要貢獻者,但是當入口熱斑在高壓渦輪上遷移時,定子排的尾流和電勢會導致溫度圓周不均勻(圖13(b))。此外,在高壓渦輪出口處有1個很強的低階分量(Wn=2),這很可能是由所選的半環模擬策略引起的混疊分量。

2 基于多波束的流場重構方法

理論上空間周期為2π 的半徑為R的任意截面處的穩態流場可以用不同波數的無限序列來描述

式中:x(θ)為沿圓周方向的流動特性;co為信號的直流分量;Wn,i為第ith個波數;Ai和φi為信號的幅度和相位。

此外,定義ai=Aicosφi和bi=Aisinφi,式(3)可以轉換為

2.1 多波束近似方法

因此,壓氣機中的圓周流動可以用N個主導波數來近似

這是重建圓周流場的重要一步,因為式(3)中未知系數的數量從無窮大減少到了2N+1個。

為了求解包含2N+1 個未知數的方程,需要在不同圓周位置測量的最少相同數量的數據點,θ=(θ1,θ2,θ3,…θm)。該系統可以描述為

式中:A為設計矩陣,維度為m×(2N+1);F為包含2N+1個未知系數的向量;x為包含來自不同圓周位置的所有m個測量數據點的向量。

A、F和x的數學表達式為

求解式(4)中描述的N個感興趣的波數。向量x中數據點的個數必須大于未知系數的個數,或者m≥2N+1。然而,在實踐中,由于x(θ)的不確定性,重建信號顯得至關重要,它包含誤差,評估重建信號的置信度。這些都需要x(θ)中的額外數據點。因此,至少需要2N+2個測量點來表征N個感興趣的波數。然而,這會產生1 個方程多于未知數的超定系統。在本研究中,最小二乘擬合方法用于求解式(6)中的未知系數。

有了所有主導波數的幅度和相位,可以使用式(5)來重構圓周流場。使用多波束近似方法利用空間欠采樣數據重構航空發動機典型部件關鍵截面流場的技術路線如圖14所示。

圖14 基于多波束方法航空發動機內復雜流場重構技術路線

2.2 條件數

采用多波束近似方法重建的流場容易因探頭測量或探頭定位的不確定性而導致x(θ)誤差。在式(6)描述的線性系統中,重構信號中的誤差受設計矩陣A的條件數k的影響。設計矩陣的條件數定義了F中相對誤差相對于x中的相對誤差的上限,即

設計矩陣條件數的變化范圍可以為1~∞。具有大條件數的系統會導致重構信號中的過大誤差。矩陣條件數的計算公式有多種,本文采用二范數進行向量和矩陣范數計算。條件數的計算公式為

式中:A+為方陣A的Moore-Penrose偽逆。

在本研究中,設計矩陣A的條件數由探頭位置θ和感興趣的波數Wn確定。知道感興趣的波數后,設計矩陣的條件數描述了探針在捕獲感興趣波數方面的分布情況。這是選擇探頭位置的最重要參數,是下一節討論的重點。因此,以確保條件數較小的方式選擇圓周位置θ,這是利用遺傳算法來實現的。選擇200 的種群大小和200 的迭代次數,優化的執行將在幾秒鐘內完成。值得注意的是,由于遺傳算法的隨機性,每次優化運行都可能導致不同的最小值。雖然這里沒有詳細說明,但本文將討論條件數的影響。

為了評估重建信號的置信度,使用皮爾遜相關系數或皮爾遜r,并計算

式中:x(θ)為測量的真實信號;xfit,j(θ)為重建信號;皮爾遜r的范圍在0~1 之間,對于重構良好的圓周流場,預測的流動特性應與所有測量位置的實際值一致,并且皮爾遜r的值接近1,反之亦然。

本文主要對方法論進行概述,詳細方法論證詳見文獻[26-27]。

3 試驗驗證

相關探針布局及流場重構技術在評估航空發動機多級壓氣機出口截面總壓分布、燃燒室出口溫度因子、高壓渦輪出口截面總溫分布的應用中得到初步試驗驗證。

3.1 多級軸流壓氣機總壓場重構

該探針布局優化方法及流場重構技術在普渡大學葉輪機械試驗室3 級軸流壓氣機中得到驗證[17]。該3 級軸流壓氣機能夠模擬航空發動機核心機末級壓氣機的流動特征(馬赫數及雷諾數)。采用該方法優化的探針布局,僅用<20%的流場信息實現了第1、2 級靜子下游總壓場的重構,多級軸流壓氣機出口總壓場場重構結果同真實值對比如圖15 所示。重構流場計算得到的均值與真實值之間的誤差小于0.1%,第1、2 級靜子下游不同葉高處基于重構流場計算得均值總壓與真實值比較見表1。

表1 第1、2級靜子下游不同葉高處基于流場重構的總壓值與真實值對比

圖15 多級軸流壓氣機出口總壓場場重構結果同真實值對比

3.2 燃燒室出口總溫場重構與溫度因子評估

該探針布局優化方法及流場重構技術在Qinetiq公司測試的典型軍用航空發動機全環形燃燒室上進行了驗證[24]。典型全環形燃燒室出口溫度試驗結果與重構結果之間對比如圖16所示。基于公開文獻[9]數據,該燃燒室出口平均溫度為2072 K,且沿著徑向、周向具有強不均勻特征,如圖16(a)所示;采用該方法優化的探針布局,僅用沿周向10 個位置的測量數據實現了燃燒室出口溫度場的重構,如圖16(b)所示。基于重構溫度場計算得到的徑向溫度分布(RTDF)與真實值(由周向144測點計算獲得)誤差小于0.5%,如圖16(c)所示。

圖16 典型全環形燃燒室出口溫度試驗結果同重構結果之間對比

3.3 高壓渦輪出口總溫場重構

該探針布局優化方法及流場重構技術在美國GE公司設計的E3高壓渦輪上進行了驗證[20]。E3高壓渦輪為2 級帶冷卻結構。采用該方法優化的探針布局,僅用周向8 支探針的測量數據實現了高壓渦輪出口溫度場的重構,重構流場計算得到的均值同真實值之間的誤差小于0.3%,如圖17所示。

圖17 E3高壓渦輪出口溫度場重構結果與真實值對比

4 總結

(1)由于受整機環境下幾何特征對探針位置的限制以及整機環境下3 維流動等因素影響,通過有限的測點數據準確獲取整機關鍵測試截面參數非常困難。目前應用比較廣泛的航空發動機關鍵截面探針布局方案針對發動機不同截面流場特性,在進口截面多采用沿周向均勻分布探針布局方案,在級間關鍵截面采用與上游靜子階梯型布置/等節距方案,但這種布局往往不能體現發動機整機的流場特征。為了解決這一難題,針對航空發動機內強3 維、強不均勻流場,介紹了一種受限空間欠采樣條件下航空發動機內流復雜流場重構技術。

(2)依據整機環境下不同部件截面的幾何特性,結合文獻及數值仿真,開展系統性的管流(進排氣系統)與環流(壓氣機、燃燒室、渦輪)氣動特征研究,形成整機環境下各部件關鍵截面特征提取方法,總結各部件關鍵截面參數變化規律。

(3)探明了主導航空發動機內流場不均勻特征的內在機理,即航空發動機內流場不均勻特征往往由少數幾個波數主導。基于以上發現,介紹了“多波束近似”的方法,該方法在重構多級壓氣機出口截面總壓分布、燃燒室出口熱斑分布、高壓渦輪出口截面總溫分布取得初步試驗驗證。

猜你喜歡
發動機
元征X-431實測:奔馳發動機編程
2015款寶馬525Li行駛中發動機熄火
2012年奔馳S600發動機故障燈偶爾點亮
發動機空中起動包線擴展試飛組織與實施
奔馳E200車發動機故障燈常亮
奔馳E260冷車時發動機抖動
新一代MTU2000發動機系列
2013年車用發動機排放控制回顧(下)
VM Motori公司新型R750發動機系列
發動機的怠速停止技術i-stop
主站蜘蛛池模板: 国产网友愉拍精品视频| a级毛片网| 精品成人一区二区三区电影| 一级毛片在线直接观看| 亚洲成人黄色在线观看| 国产精品精品视频| 国产国产人在线成免费视频狼人色| 国产精品久久自在自线观看| 伊人激情久久综合中文字幕| 国产精品自在线拍国产电影| 国产在线精品人成导航| 国产极品粉嫩小泬免费看| 国产高潮视频在线观看| 色哟哟国产精品| 国产理论最新国产精品视频| 就去吻亚洲精品国产欧美| 特级精品毛片免费观看| 日本午夜三级| 精品福利国产| 99久久国产综合精品2020| 日韩毛片免费观看| 国产AV毛片| 91精品国产无线乱码在线 | 亚洲另类国产欧美一区二区| 亚洲精品天堂自在久久77| 久久亚洲美女精品国产精品| 亚洲视频三级| 丁香亚洲综合五月天婷婷| 欧美在线黄| 黄色网页在线播放| 亚洲天堂久久| 国产在线精品美女观看| 午夜福利亚洲精品| 久久精品亚洲中文字幕乱码| 99视频全部免费| 亚洲无限乱码| 国产无遮挡猛进猛出免费软件| 伊人久久青草青青综合| 国内视频精品| 国产91精品最新在线播放| 日本高清有码人妻| 国产精品女熟高潮视频| 女同久久精品国产99国| 国产一区成人| 97人人做人人爽香蕉精品| 99r在线精品视频在线播放| 中文字幕亚洲第一| 国产一区二区三区免费| av在线无码浏览| 日韩亚洲综合在线| 久久精品波多野结衣| 亚洲人精品亚洲人成在线| 亚洲人成网站18禁动漫无码| 亚洲最黄视频| 自拍中文字幕| 国产精品偷伦在线观看| jizz国产视频| 国产中文在线亚洲精品官网| 真实国产乱子伦高清| 欧美国产在线一区| 久久免费精品琪琪| 国产高清不卡| 97青草最新免费精品视频| 亚洲成a人片77777在线播放 | 99re这里只有国产中文精品国产精品| 日韩欧美视频第一区在线观看| 日韩小视频在线观看| 中文纯内无码H| 婷五月综合| 色欲不卡无码一区二区| 99激情网| 欧美一级色视频| 青青青国产视频手机| 久久亚洲高清国产| 这里只有精品在线播放| 精品福利一区二区免费视频| 亚洲热线99精品视频| 日韩麻豆小视频| 久久精品国产精品青草app| 国内精品久久人妻无码大片高| 亚洲三级电影在线播放| 欧美午夜网|