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

艦面飛機尾噴流對進氣道溫度場影響的仿真分析

2021-10-21 12:40:38王霄程健慧沈天榮許保成孟軒
航空學報 2021年8期
關鍵詞:發動機模型

王霄,程健慧,沈天榮,許保成,孟軒

1. 國防科技大學 空天科學學院,長沙 410073

2. 航空工業沈陽飛機設計研究所,沈陽 110035

艦載戰斗機在艦面準備起飛時,若發動機產生的高溫、高速尾噴流直接噴向后方,會對后方較大范圍內的工作人員和設備造成極大的威脅,因此需要對發動機尾噴流進行有效控制。目前,航母艦載機對發動機尾噴流的控制措施主要是安裝偏流板,將尾噴流引向空中,以起到保護其后方人員和設備的作用[1]。但偏流板的存在又會帶來另外一個問題,高溫尾噴流經偏流板反射后的部分高溫氣體會回流至艦載機機腹下方和兩側,當高溫氣體回流至艦載機進氣道附近位置時,容易經進氣道被吸入發動機內,造成發動機進口溫度畸變嚴重,甚至誘發喘振,對艦載機的安全使用造成極大的威脅。

飛機尾噴流沖擊偏流板既有超聲速尾噴流的高速流動問題,又有噴流反射后回流的低速流動問題,同時還存在高速流與低速流之間的剪切流動問題。此外,高溫尾噴流與偏流板作用后,存在復雜的湍流脈動和強渦旋流場,使得噴流流場具有復雜的激波和旋渦結構,增加了流場計算和分析研究的難度。

目前,國外對于相關問題的研究手段主要分為噴流場試驗和仿真計算2個方面。針對偏流板的試驗方面,能查到的公開資料僅有俄羅斯和美國的試驗研究,俄羅斯為研究航母艦載機噴出的高溫高壓氣流對甲板和周圍飛機的影響做了大量試驗。其曾經建造了一個圍阱作為試驗臺,研究如何將熱氣流從一個柵狀噴流偏流板沿著圍阱引向旁處,試驗中由于熱氣流非但無法被導走,反而反射回來作用在飛機上導致試驗失敗,因此俄羅斯轉而對偏流板的位置、形狀和耐高溫覆蓋層開展了進一步的研究。美國在艦載機偏流板方面能查到的公開資料是其使用的為水冷式噴流偏流板,完全升起后能承受40 824 kg的噴氣推力,美國尼米茲級航母上有4組偏流板,其中有3組是MK7MODO型,每組6塊板,偏流板尺寸為14 ft×6 ft (1 ft=0.304 8 m),安裝在第1~3號彈射器的后面,其他資料不詳。國外針對噴流流場的仿真計算大多未考慮偏流板存在時的情況,主要研究的是噴流壓力場和噪聲場,以及射流沖擊平板后的壓力分布等[2-14],針對帶進氣道的偏流板前回流溫度場沒有相關仿真結果。

國內針對發動機噴流流場的研究主要集中在射流壓力場的模擬和偏流板后方溫度場的求解上。馬彩東等[15]通過流固耦合的傳熱算法研究了不同偏流板轉角對偏流板周圍流場的影響,獲得了45°為相對理想的偏轉角度的結論。何慶林等[16]以國外某型艦載機和噴氣偏流板為研究對象,采用三維雷諾平均Navier-Stokes方程、k-ε湍流模型和離散坐標模型對偏流板噴流流場進行了三維數值模擬,結果表明,偏流板對尾流場起到了很好的偏轉作用,尾噴流對噴氣偏流板背部區域的設備和人員幾乎沒有影響。趙留平[17]針對艦載機發動機噴管高溫高壓流動特性開展了仿真分析,得到了超聲速流動下中度欠膨脹和高度欠膨脹的特征和氣體射流場隨馬赫數的增加所產生的變化特性。郭凱和王強[18]針對噴管出口截面與偏流板不同距離、偏流板不同偏轉角和不同折角情況,分析了偏流板表面熱力學參數及其周圍流場的分布,結果表明偏流板與噴管出口的水平距離越遠,偏轉角越小,偏流板上的壓力與溫度就越低。張群峰等[19]研究了偏流板不同傾角時噴流對進氣道溫升的影響,結果表明偏流板傾角越大,回流對進氣道影響越小。Gao等[20]研究了艦載機發動機噴流沖擊偏流板后的流場結構,比較了不同湍流模型對射流流場結構的求解結果,Shear Stress Transport (SST)k-ω模型能更好地對可壓縮射流問題進行求解,通過仿真計算獲得了偏流板后方艦面甲板上危險區域。

綜上,目前各國針對艦載機起飛時超聲速射流沖擊偏流板問題的主要研究方向集中在射流壓力場和噪聲場上,針對雙發尾噴流沖擊偏流板后回流的流動機理和回流溫度場特征的研究很少。因此,針對上述的問題開展相關的仿真計算分析是十分必要的。

本文首先通過數值仿真計算對某型艦載機雙發尾噴流沖擊偏流板后的流動機理和溫度場特征進行了分析,然后通過非定常計算獲得了高溫超聲速尾噴流撞擊偏流板后的動態流動特性,分析了高溫回流的強渦流動結構以及流向,給出了進氣道出口的溫升率,最后針對發動機轉速不對稱、尾噴管到偏流板的距離、來流風速等參數對進氣道出口溫度畸變強度影響展開研究,為艦面環境機艦適配性分析提供參考。

1 建模背景及計算精度校核

1.1 數值計算模型

本文以某型艦載機起飛狀態時艦面環境作為仿真計算的模型,左、右進氣道按飛機順航向區分,圖1為計算模型的俯視圖和右側發動機對稱面的剖視圖,定義尾噴管入口直徑為D,偏流板尺寸為10D×5D,右側尾噴管到偏流板中心的距離大致為3D。偏流板相對地面抬起45°,相對于飛機對稱面偏轉5°。

圖1 計算模型Fig.1 Simulation model

1.2 數值方法驗證

為了驗證偏流板干擾下的發動機噴流流場湍流模型和網格設置的精確度,考慮到仿真流場關鍵部分在于射流場和沖擊場,因此選擇了相似的射流沖擊流場[21]進行數值仿真方法的校驗。

選擇射流沖擊平板的計算模型如圖2所示。噴口直徑d=25.4 mm,沖擊距離h/d=8,平板直徑為2 m。

圖2 校驗計算模型Fig.2 Checking calculation model

來流總溫與環境總溫均為288.15 K,環境壓力為1個標準大氣壓,噴管射流入口總壓為253 315 Pa,落壓比NPR=2.5。邊界條件設置為壓力出口。湍流模型分別采用了Spalart-Allmaras(S-A)、SSTk-ε、Realizablek-ε、Re-normalization Groupk-ε(RNGk-ε)模型。

為了檢驗網格無關性,保證射流區域網格長寬比基本不變的條件下,分別劃分了數量為3.7 萬(噴管徑向網格數為10)、7.3 萬(噴管徑向網格數為16)、21 萬(噴管徑向網格數為24)、90 萬(噴管徑向網格數為44)的網格。圖3展示了總網格數為90萬的網格。

圖3 校驗計算網格Fig.3 Checking calculation grid

圖4、圖5為采用不同網格總量的網格計算時得到的流場壓力分布對比,圖中Cp為壓力系數,r為徑向,p為壓力,下標“0”表示來流參數??梢钥闯鼍W格量越大,射流的激波-膨脹波交替的流場結構計算得越精細。同時,對比平板表面的徑向壓力分布發現,網格量降低到3.7 萬時,平板徑向壓力分布才相對試驗結果出現明顯的偏差。

圖4 不同網格數量時平板表面徑向壓力分布Fig.4 Radial pressure distribution on flat plate surface with different number of grids

圖5 不同網格數量時射流中心軸線壓力分布Fig.5 Pressure distribution of jet central axis with different number of grids

圖6展示了不同網格量時計算得到的子午面馬赫數(Ma)分布對比,網格越稀疏,耗散越快,達到偏流板時速度越低。

圖6 不同網格數量時子午面馬赫數分布Fig.6 Mach distribution of meridian plane with different number of grids

圖7、圖8展示了采用總量90 萬的網格時,不同湍流模型得到的流場壓力分布對比,可以看出:S-A模型獲得的平板表面徑向壓力分布結果與試驗結果最接近。

圖7 不同湍流模型下平板表面徑向壓力分布Fig.7 Radial pressure distribution on flat plate surface with different turbulence models

圖8 不同湍流模型下射流中心軸線壓力分布Fig.8 Pressure distribution of jet central axis with different turbulence models

不同湍流模型計算得到的子午面馬赫數分布如圖9所示,S-A湍流模型射流到達偏流板時速度低,但是射流的擴散角度較大。

圖9 不同湍流模型下子午面馬赫數分布Fig. 9 Mach distribution of meridian plane with different turbulence models

綜上,由于S-A模型獲得的平板表面徑向壓力分布結果與試驗結果最接近,本文選擇S-A湍流模型完成數值仿真計算,同時設置徑向網格數為44。

1.3 湍流模型和邊界設置

經過湍流模型的校驗后,計算中采用理想氣體模型,黏性系數采用Sutherland模型,湍流模型采用S-A模型。來流總溫與環境總溫均為288.15 K,環境壓力為1個大氣壓(101 325 Pa),計算遠場大小為220 m×200 m×50 m,計算域邊界設置為壓力入口和壓力出口邊界條件,如圖10 所示。通過設置進氣道出口和噴管入口的邊界條件,實現發動機不同狀態時進氣道和噴流的一體化模擬,進氣道出口設置為質量流量邊界,噴管入口設置為總溫總壓邊界,發動機不同狀態時對應具體參數見表1。求解軟件采用Fluent。

圖10 邊界條件Fig.10 Boundary conditions

表1 發動機不同狀態時的參數

1.4 計算網格設置

本文計算采用ICEM CFD軟件進行網格劃分,全模型使用六面體網格,根據對湍流模型的校驗結果,設置噴管徑向網格數為44,網格總量在2 000萬 左右。圖11展示了計算域和模型噴管出口的面網格劃分,在噴流的核心流位置進行了加密處理。

1.5 溫度場數據處理說明

溫度場數據處理按國軍標定義。面平均溫升為

ΔT2FAV=T2FAV-Tt∞

式中:T2FAV為進氣道出口平面平均總溫;Tt∞為來流總溫。

溫度畸變強度為

δT2FAV=(T2FAV-Tt∞)/Tt∞

溫升率為

T=ΔTimax/Δτm

式中:ΔTimax為最大溫升測量值;Δτm為從溫度躍升到溫升達到極值的時間。

2 尾噴流回流溫度場的回流機理

2.1 穩態仿真計算結果分析

為了摸清飛機尾噴流回流溫度場的回流機理,結合艦載機真實工作環境,針對飛機起飛滑跑前的發動機最大狀態開展了穩態數值仿真計算分析,得到的進氣道溫度畸變強度計算結果見表2,進氣道出口溫升圖譜見圖12,仿真結果表明艦面狀態兩側進氣道均吸入了一定的高溫氣體,但左側進氣道吸入的高溫氣體明顯高于右側進氣道。

圖12 進氣道出口溫升(順航向)Fig.12 Inlet outlet temperature rise(follow course)

表2 進氣道出口穩態溫度畸變強度

穩態計算得到的發動機噴流流線見圖13,可以清晰地看出雙發尾噴流相互碰撞后到達偏流板經反射后回流的整個過程:飛機地面靜止狀態、偏流板打開時,高溫高壓的尾噴流從噴管中射出以后,沿著周向具有一定的擴散角度,但雙發尾噴流互相阻滯其沿偏流板側向溢流,流態與單發射流的沖擊流場明顯不同,兩股噴流中間的通道被高溫氣流堵死。

圖13 兩側發動機噴流流線Fig.13 Jet flow streamlines of two engines

圖14給出了總溫等溫面(500 K)的速度梯度分布情況,著色變量為速度,可以清晰看到噴流交匯處的速度與主噴流相比降低了很多,因此尾噴流撞擊偏流板后,能量高的部分沿著偏流板向上逃逸,小部分由于撞擊阻滯在偏流板前的低能流不斷堆積,受偏流板向左后方的5°偏角影響,導致部分高溫氣流沿飛機左前方回彈。

圖14 總溫等溫面速度梯度Fig.14 Velocity gradient of total temperature isothermal surface

結合圖15給出的地面流線和溫度分布圖,并對比圖16中左右兩側發動機對稱面x方向速度分量,找到了左側進氣道相對右側進氣道溫升嚴重的原因:與偏流板碰撞后的高溫回流氣體向飛機左前方運動,左側發動機對稱面的初始反射速度可達-140 m/s,隨著距偏流板距離的增加,氣流的回流速度逐漸降低,在進氣口附近(L/D=15,L為距離偏流板的實際距離),受到進氣道的抽吸,氣流改變方向。高溫氣體團左側回流速度高,右側回流速低,剪切力帶來氣團的右旋,同時受進氣道抽吸作用的影響,在機頭前卷起一個高溫回流渦,導致左、右兩側進氣道均吸入了一定的高溫氣體,由于高溫氣流沿飛機左側回彈,左側進氣道的溫升更嚴重。

圖15 噴流反射地面流線圖Fig.15 Diagram of jet ground streamlines

圖16 噴管鉛垂對稱面回流速度Fig.16 Recirculation velocity of nozzle symmetry plane

2.2 動態仿真計算結果分析

考慮到發動機尾噴流與偏流板碰撞后的回流場的產生與發展具有典型的非定常流動特性,為了觀察高溫回流流場的生成和發展,以及獲得進氣道出口平面的溫升率,在2.1節的基礎上展開了非定常數值仿真計算研究,非定常計算采用了與定常計算同樣的網格、來流條件和湍流模型,非定常時間步長為0.000 1 s,每個時間步內迭代次數為50 步。非定常計算的初始流場為發動機暖機狀態對應的定常計算收斂流場,非定常計算時發動機噴管進口的總溫和總壓條件從暖機狀態到發動機最大狀態按階躍形式給出。

動態仿真計算得到的進氣道出口溫度畸變強度隨時間變化結果如圖17所示,進氣道出口溫度畸變強度隨時間先增加后降低,左右兩側進氣道出口均會出現一個極值,左側進氣道非定常計算得到的出口溫度畸變極值與定常計算結果基本一致,非定常計算得到了左側進氣道溫度畸變強度在0.70 s左右開始上升,在1.73 s達到極值狀態,溫升率接近100 K/s。

圖17 進氣道出口溫度畸變強度隨時間變化Fig.17 Variation of total temperature distortion intensity of inlet outlet with time

為了觀察0.70~6.00 s區間進氣道吸入高溫回流導致進氣道出口的溫度畸變逐漸增加到極值又降低趨于穩定這一過程,圖18展示了總溫350 K等值面的高溫氣團的回流運動軌跡,以及對應時刻進氣道出口溫升圖譜,可以看到射流沖擊偏流板后,回流場逐漸向來流方向移動和擴大,由于存在側偏角,同時受到進氣道抽吸影響,回流在發展過程中,具有明顯的不對稱性?;亓髟娇拷M氣口,受到抽吸作用的影響越大,在越過進氣口時(t=1.70 s)時,進氣道出口溫度畸變達到了極值狀態,對應進氣道出口高溫區的位置與穩態計算也基本一致,回流一方面繼續前傳,一方面被位于其一側的進氣口吸引,形成了卷吸渦,靠近回流中心的左側進氣道受影響更大,這與圖18顯示的左側進氣道出口溫度高于右側進氣道的結果相吻合。

圖18 回流高溫氣團隨時間的運動及進氣道出口溫升Fig.18 Movement of return high temperature air mass with time and inlet outlet temperature rise

動態仿真計算得到的進氣道溫度畸變極值為35.5%,與定常計算36%的結果基本一致,進氣道出口溫度圖譜也基本一致,右側進氣道在非定常計算中捕捉到了高于定常計算10%左右的溫度畸變強度,但畸變強度極值仍低于左側進氣道的仿真結果,考慮到實際工作中,一般根據進氣道出口的極值溫度到發動機可抗極值溫度之間的范圍來評估發動機抗溫度畸變的穩定工作裕度,定常計算捕捉到了非定常仿真的溫度畸變極值,因此可采用穩態仿真計算開展后續偏流板反射回流場的分析。

3 參數影響分析

3.1 噴流不對稱影響分析

艦載機在起飛過程中,由于發動機的差異性,左右兩側發動機的噴流狀態往往不完美對稱,相對于理論值會存在一定的偏差,針對第2節對基準流場的研究發現,左右兩側噴流的阻滯作用導致高溫廢氣的堆積,考慮到艦載機的安全性,采用穩態仿真計算分析左右兩側噴流不對稱狀態時候,進氣道吸入高溫氣體的風險,為飛行員的操縱提供一定的指導。仿真計算在2.1節的基礎上開展,保持左側發動機噴流流量不變,右側發動機噴流流量降低了5%左右,計算結果見表3,與表2的計算結果相比,左側進氣道溫升相對基準狀態變化不大,右側進氣道溫升提高。

表3 非對稱噴流進氣道出口穩態溫度畸變強度

圖19給出了非對稱噴流狀態的地面流線圖,噴流回流流場初始高溫區的大小和位置與圖15基準流場對稱噴流狀態接近,不同的是,隨著高溫回流逐漸向機頭方向流動,非對稱噴流狀態受進氣道吸力影響更嚴重,高溫氣團回流到進氣口,導致右側進氣道吸入了高溫氣體。結合圖20中非對稱噴流狀態和對稱噴流狀態回流場不同剖面的總溫云圖發現,由于右發噴流流量降低,導致其對左發噴流的擠壓作用變弱,高溫回流氣體團開始向上卷起,高溫回流受進氣道抽吸作用更明顯,對進氣道影響更大,導致右側進氣道也吸入了更多的高溫氣體。

圖19 非對稱狀態噴流反射地面流線Fig.19 Ground streamlines of asymmetric jet

圖20 回流場不同剖面總溫對比(左:對稱噴流,右:非對稱噴流,順航向)Fig.20 Comparison of total temperature on different sections of backflow field (left: symmetric jet, right: asymmetric jet, follow course)

3.2 尾噴口到偏流板距離影響分析

艦載機真實工作環境中,為了滿足不同的起飛條件,一般配備相對飛機起飛站位不同距離的偏流板,因此,本文針對發動機噴口距離偏流板不同的起飛距離展開了定常數值仿真研究,兩側發動機均為最大狀態。在2.1節的基礎上,將偏流板沿飛機軸向前后分別移動了一定的距離,進氣道出口溫升的仿真計算結果見圖21,仿真結果表明,隨著尾噴口到偏流板距離的增加,左側進氣道出口溫度畸變強度先增加再降低。

圖21 進氣道出口溫度畸變強度隨尾噴口到偏流板距離的變化Fig.21 Variation of temperature distortion intensity at inlet outlet with distance from nozzle tail to jet blast deflector

圖22通過對比不同尾噴口到偏流板距離時同一總溫等溫面發現,2D站位與3D站位相比,尾噴流達到偏流板上后,由于兩股噴流間距離較大,大部分噴流可沿著通道向板上逃逸,高溫氣體回流量明顯減少。3D站位由于噴流距離更長,噴流經長距離膨脹,到達偏流板時兩股噴流基本相連,導致噴流中心以下的高溫氣流被噴流主流封閉了向上的空間,形成大量噴流回流。當尾噴管到偏流板的距離進一步增加到6D時,折返射流的強度已經大大減弱,回流速度降低,從而很難到達進氣口附近,并且其溫度在回流至進氣口的路徑中由于摻混作用已經有很大下降,使得進氣道出口的總溫升明顯下降。

圖22 尾噴口到偏流板不同距離時的回流溫度場Fig.22 Temperature field of backflow at different distances from nozzle tail to jet blast deflector

因此,偏流板的距離決定了回流場整體的強度與分布,以及進氣口位置影響了回流的狀態,從而決定了進氣道抽吸流場與回流場的耦合特性。

3.3 風速影響分析

艦載機真實工作環境中,航母一般具有25~30 kn(1 kn=0.514 m/s)的行進速度,即相對飛機坐標系來流風速在15 m/s左右,為此,針對不同行進速度下發動機尾噴流對飛機周邊流場的影響開展了定常仿真計算研究,尾噴流距離偏流板距離選擇了3.2節中進氣道受影響程度最大的3D距離,計算結果見圖23,左右兩側進氣道出口溫度畸變強度均隨著來流速度的增加先增加后減小,來流風速在20 m/s時,左右兩側進氣道的溫升均達到了極值。

圖23 進氣道出口溫度畸變強度隨來流風速變化Fig.23 Variation of temperature distortion intensity at inlet outlet with inflow wind speed

從圖24可以看到,高溫尾噴流經偏流板反射后沿著左側進氣道的左前方向前傳播,受到來流壓迫,高溫回流形成上洗渦流,隨著來流速度的增加,將這個渦不斷地向后壓縮,距離進氣道進口越來越近,導致了更多的高溫氣流被吸入左右兩側進氣道,當來流風速達到一定值時,剛好將這個高溫回流渦壓縮到進氣道進口附近,導致左右側進氣道吸入了大量高溫氣體,此時如果風速繼續增大,回流影響區繼續向后收縮,將不會對進氣道流場產生影響。

圖24 20 m/s風速時噴流反射地面流線Fig.20 Jet ground streamlines with 20 m/s wind speed

因此來流風速對進氣道吸入高溫氣體的影響與回流高溫氣體團距離進氣口位置距離相關,回流高溫氣團作為低能流本身具備一定能量,隨著風速的升高,對進氣道的影響先增加再減弱。

4 結 論

本文采用經過試驗數據驗證的數值仿真計算方法完成了飛機尾噴流沖擊偏流板后回流溫度場的穩、動態仿真計算分析。

1) 通過仿真分析發現雙發尾噴流相互干擾阻滯導致高溫回流尾氣在偏流板前堆積,受到噴流的回流引射向飛機機頭方向運動,在發動機的抽吸作用下,導致進氣道吸入回流出現高溫溫升,但左側溫升更明顯。

2) 通過非定常仿真計算,得到了進氣道出口溫度畸變隨時間先增加后減小,最后穩定的變化趨勢,進氣道出口的溫升率極值可達100 K/s。

3) 通過進一步的仿真研究,得到了噴流條件、風速、尾噴口到偏流板的距離等參數變化對進氣道溫度畸變強度的影響規律,獲得了尾噴口到偏流板的距離對回流場整體的強度與分布起決定作用,以及進氣口的位置影響了進氣道抽吸流場與回流場的耦合特性這一結論。

猜你喜歡
發動機模型
一半模型
元征X-431實測:奔馳發動機編程
2015款寶馬525Li行駛中發動機熄火
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
新一代MTU2000發動機系列
發動機的怠速停止技術i-stop
新型1.5L-Eco-Boost發動機
主站蜘蛛池模板: 国产成人麻豆精品| 一区二区三区四区在线| 日本少妇又色又爽又高潮| 国产精品污污在线观看网站| 天堂成人在线| 五月婷婷亚洲综合| 欧美日韩北条麻妃一区二区| 国模粉嫩小泬视频在线观看| 97亚洲色综久久精品| 亚洲第一区在线| 国产欧美在线| 国产成+人+综合+亚洲欧美| 亚洲国产午夜精华无码福利| 国产成人狂喷潮在线观看2345| 婷婷成人综合| 老色鬼久久亚洲AV综合| 国产一区二区三区免费| 欧美成人午夜影院| 91偷拍一区| 欧美日韩亚洲国产| 国产偷国产偷在线高清| 日韩成人午夜| 日韩精品无码免费专网站| 激情在线网| 91网在线| 亚洲中文字幕在线精品一区| 亚洲成人在线网| 日韩欧美中文| 久久国产高潮流白浆免费观看| 日韩精品一区二区深田咏美| 日韩精品无码免费一区二区三区| 国产白浆视频| 成人小视频网| 无码区日韩专区免费系列| 青青久视频| 亚洲九九视频| 免费在线不卡视频| 成人福利一区二区视频在线| 国产国拍精品视频免费看| 国产精品分类视频分类一区| 色综合久久88| 欧日韩在线不卡视频| 亚洲 日韩 激情 无码 中出| 2021国产v亚洲v天堂无码| 国产毛片片精品天天看视频| 99国产精品国产高清一区二区| 99999久久久久久亚洲| 国产不卡国语在线| 国产国语一级毛片| 色综合久久无码网| 91色国产在线| 久久国产精品娇妻素人| 911亚洲精品| 国产精品理论片| 国产一级毛片网站| 操国产美女| 亚洲国产亚洲综合在线尤物| 午夜久久影院| 欧美爱爱网| 四虎成人精品| 成人噜噜噜视频在线观看| 亚洲精品国产成人7777| 国产主播一区二区三区| 午夜福利亚洲精品| 欧美成人亚洲综合精品欧美激情| 浮力影院国产第一页| av大片在线无码免费| 午夜福利在线观看成人| 亚洲手机在线| 亚洲妓女综合网995久久| 国产精品精品视频| 先锋资源久久| 亚洲无码视频一区二区三区| 一级毛片在线免费视频| 国产精品对白刺激| 国产黄网永久免费| 凹凸国产分类在线观看| 国产精品一区二区在线播放| 国产经典三级在线| 日韩 欧美 小说 综合网 另类| 波多野结衣中文字幕一区| 亚洲一区二区日韩欧美gif|