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

跨聲速軸流壓氣機多葉排反問題優化方法

2016-05-07 06:33:42劉昭威吳虎唐曉毅
西北工業大學學報 2016年1期
關鍵詞:優化設計

劉昭威, 吳虎, 唐曉毅

(西北工業大學 動力與能源學院, 陜西 西安 710072)

?

跨聲速軸流壓氣機多葉排反問題優化方法

劉昭威, 吳虎, 唐曉毅

(西北工業大學 動力與能源學院, 陜西 西安710072)

摘要:為提升跨聲速壓氣機的氣動性能,以葉輪機械三維黏性反問題設計理論為基礎,發展了適用于壓氣機多葉排三維反問題優化設計方法。描述了壓氣機多葉排流場計算所采用的數值求解技術,詳細推導了反問題優化設計方法和流程。為驗證方法的正確性,運用德國宇航中心單級跨聲速壓氣機R030-SUKS/31的實驗數據與計算結果進行對比。在對原型葉片表面載荷進行分析的基礎上,重新修改葉片表面載荷分布,并通過反問題設計方法計算得到新的葉片幾何構型。結果表明,通過修改葉片表面載荷分布,運用反問題設計方法得到的壓氣機,其進口相對馬赫數有所提升,壓氣機流量和壓比分別提高了3.5%和2.0%,氣動性能有明顯提升,驗證了方法的正確性和有效性。

關鍵詞:反問題設計方法;載荷分布;優化設計;跨聲速壓氣機

隨著航空發動機總體性能的提升,風扇以及增壓級的設計難度不斷增大。為了提高發動機的推重比,降低單位燃油消耗率,風扇和壓氣機的設計壓比不斷提高,壓氣機設計朝著高負荷,高葉尖馬赫數的方向發展。當前較為先進的航空發動機風扇和壓氣機前幾級均采用超聲速或跨聲速設計,其最大的特點就是葉片通道內部存在激波,氣體流動的三維效應較強,流動損失的產生原因復雜。因此,跨聲速壓氣機的氣動優化設計十分重要。

20世紀80年代初,Denton[1]首次將全三維計算流體力學(computational fluid dynamics,CFD)技術應用于葉輪機械設計領域。經過多年的發展,CFD在葉輪機械設計中的應用已日臻成熟。其主要作用是對葉輪機械內部流場進行數值模擬和性能預測,并以此為依據對葉片幾何構型進行優化設計,以期得到較好的氣動性能。這種已知葉片幾何構型求解氣動性能的過程通常被稱為正問題求解。在實際設計過程中,為了得到一個氣動性能較好的葉片幾何構型,往往需要對葉片型面坐標進行不斷調整,反復計算,這需要耗費大量的時間和計算資源。因此,設計人員提出了一種全新的設計方法,通過給定葉輪機械的氣動性能參數來求解葉片的幾何構型,即反問題設計方法。

反問題設計方法從20世紀80年代中期開始起步,近期已經成為國際上研究的熱點問題。最初的反問題求解的是二維勢流方程,隨著正問題求解方法的不斷改進,反問題設計也逐步擴展到求解Euler方程和Navier-Stokes方程。Dang[2]采用4階 Runge-Kutta時間推進、基于網格中心的有限體積法求解二維Euler定常流場,發展了反問題所適用的滲透邊界條件,并建立了一套二維反問題設計方法體系。Ghaly等[3]通過給定目標靜壓,計算虛位移的方法進行葉片幾何構型更新,并應用于壓氣機葉柵的改型設計,取得了良好效果。隨著三維CFD技術的逐步成熟,三維反問題設計方法不斷發展。Qiu等[4]基于非均勻有理B樣條曲線(NURBS)和最小二乘的思想,提出了一種統一的葉片中弧面生成方法,增強了方法的魯棒性。

國內對反問題設計方法的研究起步并不晚,西北工業大學周新海等[5]基于有限體積法求解Euler方程進行了跨聲速葉柵的反問題求解方法研究。北京理工大學楊策等[6]采用規定葉片表面無量綱目標速度分布,通過比較目標速度分布和計算獲得的速度分布來修改葉片壓力面和吸力面的坐標,最終

獲得滿足要求的葉片形狀。中科院工程熱物理所的王正明等[7-8]從1985年就開始研究二維反問題,并在2000年發展了黏性的三維反問題設計方法,該方法基于非正交曲線坐標系下完全守恒型Navier-Stokes方程和B-L湍流模型,運用MacCormark顯式時間推進,并運用該方法進行葉型設計。楊金廣[9]采用有限體積法求解三維Navier-Stokes方程,并基于滲透邊界條件以及葉片表面載荷分布來改變葉片形狀,建立了三維黏性葉輪機械反問題設計體系。

目前的葉輪機械反問題設計主要針對平面葉柵和單葉排葉片幾何構型,將反問題應用于多葉排乃至多級葉輪機械設計的研究較少。國外方面,van Rooij等[10]利用反問題設計方法研究多級壓氣機多葉片排優化改型設計,取得了良好效果。國內對于葉輪機械多葉排反問題的研究較為欠缺,從公開發表的文獻來看,尚未有較為顯著的研究成果。

本文在課題組前期的研究基礎上,采用反問題設計方法,對單級雙葉片排跨聲速壓氣機進行反問題改型設計研究。文中介紹所使用的反問題設計方法的計算模型,以德國宇航中心(Deutsches Zentrum für Luft-und Raumfahrt,DLR)單級跨聲速壓氣機R030-SUKS/31的實驗數據為依據驗證求解器的正確性,并運用葉片返回實驗驗證反問題求解的精度和有效性。最后針對該跨聲速壓氣機流場特點對葉片表面載荷分布進行優化改進,運用反問題設計方法求出新的葉片幾何構型,從而優化葉片通道中的流動結構和激波強度,并通過對比流場及轉子出口氣動性能參數,驗證方法的有效性。

1模型和計算方法

1.1控制方程及數值求解方法

反問題設計方法所求解的控制方程與正問題相同,本文中采用求解積分形式的Navier-Stokes方程,其具體形式如下:

(1)

式中,Ω為控制體體積,?Ω為控制體邊界,n為控制體邊界外法向量,W為通量變量,Fl為對流通量,Fv為黏性通量。本文采用基于網格中心點的有限體積法對控制方程進行離散求解,運用Jameson所發展的改進JST格式進行通量求解,時間推進采用5步混合Runge-Kutta法,CFL數達到3.6,湍流模型采用S-A模型,為了保證求解的穩定性,只采用了當地時間步長一項加速收斂措施。

在葉輪機械多葉排流場數值模擬中,動、靜葉片排之間的摻混面(mixing-plane)的處理較為重要,摻混面不僅需要正確傳遞上、下游葉排流場信息,還必須盡量消除特征波反射對數值求解的影響。Giles[11]基于線性Euler方程提出了一整套無反射邊界條件理論,Chima[12]將該套理論應用于葉輪機械多葉排摻混面的處理,取得了良好效果。本文所采用的摻混面數值計算方法均基于上述理論。作者在前期工作中基于上述理論,自主開發了一套壓氣機全三維黏性流場數值計算程序[13],能夠較為準確地預測跨聲速軸流壓氣機的氣動性能,本文的工作均是該套計算程序的基礎上開展的。

1.2反問題計算方法

反問題設計方法的本質是通過給定葉片周圍流場的氣動參數分布從而得到葉型幾何數據,氣動參數的分布由設計者給定。本文中所給定的氣動參數為葉片表面沿軸向的載荷分布,即葉片吸力面與壓力面之間的壓力差 ,Dang[2]對葉片載荷與葉片對氣流的加功量之間關系進行了詳細推導。反問題計算的核心是反問題邊界條件。本文所采用反問題邊界條件基于無滑移物面邊界條件,保證葉片幾何構型更新前后葉片表面的特征值守恒。通過所給定的葉片表面載荷 與迭代計算過程中的載荷之差,求解壁面的虛擬移動速度 。

(2)

式中正、負號分別代表葉片的上、下表面,ρ和c分別為鄰近壁面處的氣體密度和聲速。令葉片中弧面的虛擬法向速度Wn與葉片上下表面相等

(3)

用虛擬速度乘以相應的時間步長,則可以求出葉片中弧面的虛擬位移δf

(4)

由于時間步長與實際流場中的參數變化相關,因此選取時間步長Δt為鄰近壁面的當地時間步長,nθ為葉片中弧面法向量的切向分量,計算過程中,保證葉片的軸向和徑向坐標不改變,只改變切向坐標,則更新后的葉片中弧面坐標可表示為

(5)

計算時周期性的更新中弧面,選取中弧面上某一條徑向網格線作為積疊線,保持其切向坐標不變,更新中弧面時從該條積疊線開始,順流動方向或者逆流動方向分別求解更新后的中弧面坐標。

保持原始葉片厚度不變,將其疊加至更新后的中弧面,最終求出更新后的葉片吸、壓力面坐標。運用更新后的葉片幾何構型重新進行網格劃分,繼續進行計算,直至流場收斂。改進后的反問題求解流程如圖1所示。

圖1 反問題設計計算流程

2計算結果及分析

2.1正問題求解驗證

為了保證反問題計算的正確性,首先需對正問題求解器進行驗證。本文首先應用所發展的正問題計算程序對德國宇航中心單級跨聲速壓氣機R030-SUKS/31近峰值效率工況點進行氣動計算,并將計算結果與實驗數據[14]進行對比,表1為該壓氣機的實驗工況。

表1 R030-SUKS/31跨聲速壓氣機級實驗工況

計算給定的收斂標準為連續方程通量變量殘差收斂至 ,進、出口流量偏差小于0.5%。表2為R030-SUKS/31單級壓氣機峰值效率點附近氣動參數計算值和實驗值的對比,可見計算結果與實驗值符合較好,計算精度能夠滿足工程實用要求。

表2    R030-SUKS/31跨聲速壓氣機級氣動性能

2.2葉片返回實驗

為了驗證反問題求解的正確性,本文進行了反問題葉片返回實驗。將原始葉片的載荷分布作為目標載荷,以保證反問題解的存在性,在保證葉片厚度分布不變的基礎上通過人為改變中弧面生成新的葉片,將新葉片作為初始幾何,運用反問題計算使之穩定收斂到原始葉片幾何構型,證明收斂解的唯一性以及反問題求解方法的可靠性。反問題計算的收斂標準以及進、出口邊界條件與正問題求解相同。求解初始階段給定前文中正問題計算所得到的葉片表面載荷分布,葉片的厚度分布以及一個初始的中弧面坐標。

圖2為反問題計算最終收斂后3個不同展向位置葉片截面型線之間的對比,左側為轉子,右側為靜子。由圖中可以看出,在不同的展向位置,葉片截面型線均與原始葉型符合較好,可見反問題計算的結果與原始葉型幾何符合較好。同時,反問題計算結果很好的滿足了給定的葉片表面載荷分布。葉片返回實驗驗證了反問題計算程序的精度和可靠性。

圖2 葉片返回計算結果與原始葉型對比

2.3多葉排反問題改型優化設計

反問題設計方法的最大優勢在于其直觀性,設計人員可以根據設計意圖,給定氣動參數分布,求得滿足氣動參數分布的葉片幾何構型。本文采用反問題設計方法,對德國宇航中心單級跨聲速壓氣機R030-SUKS/31峰值效率點附近的葉片表面載荷進行優化改型設計,力求提升該跨聲速壓氣機氣動性能。

跨聲速壓氣機轉子葉片彎角較小,與亞聲速壓氣機的增壓原理不同,跨聲速壓氣機轉子主要依靠葉片通道內的激波進行增壓,且波前馬赫數越大,激波強度越高,增壓能力越強,但過高的激波強度會造成較大的激波損失。而激波在葉片通道內的位置與強度能夠較為直觀的體現在葉片表面的載荷分布上。

圖3為R030-SUKS/31單級跨聲速壓氣機轉子、靜子葉片載荷分布等值線圖。圖中等值線分布的疏密程度代表了載荷梯度的大小。從圖中的葉片表面載荷分布可以看出,轉子的葉片載荷從葉根到葉尖逐漸增大,在50%葉高以上的部分,載荷沿流動方向先增大后減小,葉片通道內激波與吸力面和壓力面相交的位置與載荷變化較為劇烈的位置相對應。而靜子葉片載荷沿徑向的分布變化不大,該壓氣機靜子葉片采用可控擴散葉型,載荷主要分布在葉片前半部分。

圖3 原始葉片表面載荷分布等值線圖

由于該單級壓氣機峰值效率已經達到90%,該工況點效率的提升空間有限,因此對其優化改型的主要目標是在保持效率基本不變,增大該工況下的流量和壓比。在對轉子葉片載荷進行調整時,首先保證葉片通道內激波與葉片的相對位置不發生明顯變化,在此基礎上適當增加來流相對馬赫數,一方面增加了激波強度和壓比,另一方面也增大了流量。為達到上述目的,轉子葉片載荷的調整策略為:在靠近葉尖附近,提高葉片前緣載荷以減小葉片吸力面的靜壓,提高進口相對馬赫數,減小葉片中部的載荷峰值大小以減小壓力面的逆壓梯度,降低附面層分離的風險,由于原型壓氣機轉子靠近葉片尾緣處的載荷較小,在改型時提高了尾緣處的載荷以增加加功量;轉子葉片中部以下的載荷調整保持原始載荷分布形式不變,對載荷沿軸向的分布曲線進行光順,使得載荷沿流動方向變化較為平緩,減小逆壓梯度導致的附面層損失。靜子葉片的調整策略為:不改變原型的載荷分布形式,減小各個截面上的載荷峰值大小,同時對載荷曲線進行光順,使載荷更為均勻的加載到整個葉片。

圖4 改型前、后3個不同半徑處葉片型線對比

基于上述多葉排反問題載荷調整策略,本文采用樣條曲線對原葉片表面軸向載荷分布進行重新構造。圖5為3個不同徑向位置處壓力載荷曲線的對比,圖中實線為改型設計所構造的載荷曲線,虛線為原型載荷曲線,圓圈為反問題計算得到的載荷曲線。從圖中可以看出,改型后的載荷分布與給定的載荷定解條件吻合較好,反問題計算結果很好的滿足了改型設計意圖。與原型壓氣機載荷曲線相比,改型后的載荷曲線在保證了原有載荷分布形式的基礎上更為光滑,變化更為平緩。圖4為3個截面葉片型線改型前、后對比,可以看出,在載荷定解條件的驅動下,葉片型線較原型發生了改變。

圖6為原型壓氣機與優化改型后的壓氣機在3個不同徑向位置對馬赫數等值線圖對比。通過對比可以看出,在80%和50%2個徑向位置處,由于改變了轉子葉片前緣附近的載荷分布,使得葉片進口區域相對馬赫數有所提升,相應提高了進口流量,同時按照激波理論,波前馬赫數提高使得激波強度增強,因此壓氣機轉子的增壓能力有所增強。在20%徑向位置處的葉片載荷調整量較小,故流場馬赫數分布變化不大。由于在構造靜子葉片的載荷時保證了原有的載荷分布形式,因此靜子葉片通道內流動的馬赫數分布未發生較大變化,由于上游轉子葉片出口流場發生了改變,影響到了下游的靜子葉片,為了保證給定的載荷分布,靜子葉片的型線發生了一定改變。

圖5 原始葉片構型、反問題輸入以及反問題計算結果葉片3個不同半徑處載荷分布曲線對比

圖6 改型前、后葉片3個不同半徑處相對馬赫數分布等值線圖對比

圖7為改型后的葉片載荷分布等值線圖,與原型對比可以看出,轉子葉片葉尖處葉片前緣載荷有所提升,葉片中部的載荷峰值有所降低,激波位置與原型沒有發生較大變化。葉片中部的載荷沿流動方向分布更加均勻,靠近尾緣處的載荷有所提升。靜子葉片的載荷分布形式未發生較大變化,葉片前部載荷最大值有所減小,但高載荷區域面積較原型有所增大。從葉片載荷分布來看,通過反問題計算,葉片表面載荷分布基本滿足了初始設計意圖。

圖7 改型后葉片表面載荷分布等值線圖

參數原型改型流量/(kg·s-1)17.3217.93總壓比1.5521.584效率0.9050.901

表3為改型前、后壓氣機出口截面氣動性能參數對比,可以看出,由于提高了進口相對馬赫數,轉子的流通能力增強,流量比原型提高了3.5%,同時進口相對馬赫數提高使得激波強度增大,提升了轉

子的增壓能力,壓比較原型提高了2.0%。隨著流量和壓比的提高,壓氣機轉子負荷增大,因此絕熱效率有所下降,但下降幅度較小。從工程實用的角度出發,整個單級壓氣機的氣動性能有明顯提升。

3結論

本文以反問題設計理論為基礎,發展了跨聲速壓氣機多葉排三維黏性反問題改型設計方法,在運用實驗結果驗證了該程序的正確性后,將其運用到單級跨聲速壓氣機多葉排改型設計上,通過對比,得出以下結論:

1) 本文所發展的正問題求解器所計算出的單級跨聲速壓氣機氣動特性與實驗值符合較好,證明了正問題求解器特別是摻混面處理方法的可靠性。

2) 從葉片返回實驗的結果可以看出,反問題計算的能夠較好的還原原始葉片幾何構型,驗證了本文所采用的反問題計算方法是正確有效的。

3) 通過對葉片表面軸向載荷分布進行合理調整,能夠有效提升壓氣機進口相對馬赫數,從而增大壓氣機流通能力和激波強度,提高壓氣機的流量和壓比。文中算例改型后的壓氣機流量提高了3.5%,壓比增加了2.0%,而效率基本不變,改型設計提高了壓氣機轉子的氣動性能。

參考文獻:

[1]Denton J D. The Calculation of Three-Dimensional Viscous Flow through Multistage Turbomachines[J]. ASME Journal of Turbomachinery, 1992, 114(1): 18-26

[2]Dang T Q. A Fully Three-Dimensional Inverse Method for Turbomachinery Blading in Transonic Flows[J]. ASME Journal of Turbomachinery, 1993, 115(2): 354-361

[3]Ahmadi M, Ghaly W. Aerodynamic Inverse Design of Turbomachinery Cascades Using a Finite Volume Method on Unstructured Meshes[J]. Inverse Problems in Science and Engineering, 1998, 6(4):281- 298

[4]Qiu X, Dang T. 3D Inverse Method for Turbomachine Blading with Splitter Blades[R]. ASME-2000-GT-0526

[5]周新海, 朱方元. 求解葉柵跨音速流動反問題的有限體積方法[J]. 工程熱物理學報, 1985, 6(4):331-335

Zhou Xinhai, Zhu Fangyuan. Finite Volume Method to Solve The Inverse Problem for Transonic Flows in Cascades[J]. Journal of Engineering Thermophysics, 1985, 6(4):331-335 (in Chinese)

[6]楊策,老大中,蔣滋康. 求解跨聲速壓氣機葉柵黏性流動反問題的數值解[J]. 推進技術, 1999,20(4):57-60

Yang Ce, Lao Dazhong, Jiang Zikang. Numerical Soluation on Viscous Inverse Problem for Transonic Compressor Cascades[J]. Journal of Propulsion Technology, 1999, 20(4): 57-60 (in Chinese)

[7]王正明,賈希誠. 正反問題數值解法相結合三維葉片的優化設計[J]. 工程熱物理學報, 2000, 21(5): 567-569

Wang Zhengming, Jia Xicheng. Optimum Design of Three-Dimensional Blades by Combination of Numerical Methods for Solving Direct and Inverse Problems[J]. Journal of Engineering Thermophysics, 2000, 21(5): 567-569 (in Chinese)

[8]Wang Zhengming, Cai Ruixian. A Three-Dimensional Inverse Method Using Navier-Stokes Equation for Turbomachinery Blading[J]. Inverse Problems in Engineering, 2000(8): 529-551

[9]楊金廣. 葉輪機械全三維黏性反方法設計技術研究[D]. 西安:西北工業大學, 2013

Yang Jinguang. Fully Three Dimensional Viscous Inverse Method for Turbomachinery Aerodynamic Design[D]. Xi′an: Northwestern Polytechnical University, 2013 (in Chinese)

[10] van Rooij M P C, Dang T Q, Larosiliere L M. Enhanced Blade Row Matching Capabilities via 3D Multistage Inverse Design and Pressure Loading Manager[R]. ASME Paper GT2008-50539, 2008

[11] Giles M B. Nonreflecting Boundary Conditions for Euler Equation Calculations[J]. AIAA Journal, 1990, 28(12): 2050-2058

[12] Chima R V. Calculation of Multistage Turbomachinery Using Steady Characteristic Boundary Conditions[R]. AIAA-1998-0968

[13] 劉昭威, 吳虎. 改進的反問題邊界條件在葉輪機械中的應用[J]. 工程熱物理學報, 2015, 36(10): 1-5

Liu Zhaowei, Wu Hu. Application of Improved Inverse Method Boundary Condition in Turbomachinery[J]. Journal of Engineering Thermophysics, 2015, 36(10): 1-5 (in Chinese)

[14] Dunker R, Rechter H, Starken H, et al. Redesign and Performance Analysis of a Transonic Axial Compressor Stator and Equivalent Plane Cascades with Subsonic Controlled Diffusion Blades[J]. Journal of Engineering for Gas Turbines and Power, 1984, 106(2): 279-287

[15] 劉昭威, 吳虎, 唐曉毅. 基于反問題設計方法的葉柵激波損失控制[J]. 推進技術, 2014,35(6):766-773

Liu Zhaowei, Wu Hu, Tang Xiaoyi. Numerical Soluation on Viscous Inverse Problem for Transonic Compressor Cascades[J]. Journal of Propulsion Technology, 2014,35(6): 766-773 (in Chinese)

Optimization of Transonic Axial Compressor Using Multi-Row Inverse Method

Liu Zhaowei, Wu Hu, Tang Xiaoyi

(Department of Aero-Engines, Northwestern Polytechnical University, Xi′an 710072, China)

Abstract:In order to improve the aerodynamic performance of transonic compressor, a compressor inverse design program is developed based on the turbomachnary three-dimensional viscous inverse design method theories. The numerical solution methods of Multi-Row transonic axial compressor flow field are introduced first and then the inverse design method and process is described in detail. The DLR R030-SUKS/31 single-stage transonic axial compressor experimental data are used to validate the calculation results. On the basis of the analysis of original blade loading distribution, the blade geometry is redesigned by applying the inverse method with the modified loading distribution. The numerical results show an increase in the compressor inlet Mach number. The mass flow rate and pressure ratio are increased respectively by 3.5% and 2.0%. Better performance of the redesigned compressor is achieved; this demonstrates the effectiveness of this method.

Keywords:inverse method, loading distribution, optimization design, transonic compressor

中圖分類號:V235.13

文獻標志碼:A

文章編號:1000-2758(2016)01-0118-07

作者簡介:劉昭威(1986—),西北工業大學博士生,主要從事葉輪機械氣動熱力學的研究。

基金項目:國家自然科學基金(51076131)資助

收稿日期:2015-04-23

猜你喜歡
優化設計
導彈舵面的復合材料設計與分析
航空兵器(2016年4期)2016-11-28 21:47:29
礦井主排水系統的優化設計與改造
科技資訊(2016年19期)2016-11-15 08:34:13
數據挖掘對教學管理的優化設計
如何實現小學數學課堂練習設計優化
文理導航(2016年30期)2016-11-12 14:56:57
淺析人機工程學在家具創作中的作用
試析機械結構優化設計的應用及趨勢
汽車行李箱蓋鉸鏈機構的分析及優化
東林煤礦保護層開采卸壓瓦斯抽采優化設計
橋式起重機主梁結構分析和優化設計
對無線傳感器網絡MAC層協議優化的研究與設計
科技視界(2016年22期)2016-10-18 15:25:08
主站蜘蛛池模板: 中文字幕欧美日韩高清| 2021天堂在线亚洲精品专区| 97se亚洲综合不卡| 国产中文在线亚洲精品官网| 亚洲福利视频网址| 国产欧美日韩资源在线观看| 久久先锋资源| 亚洲无限乱码一二三四区| 国内精品小视频福利网址| 国产综合欧美| 亚洲精品麻豆| 99re在线免费视频| 波多野结衣一区二区三视频| 午夜视频www| 香蕉蕉亚亚洲aav综合| 亚洲天堂视频在线观看免费| 高潮毛片免费观看| 在线观看国产网址你懂的| 日韩无码白| 九色视频线上播放| 国产男人天堂| 免费不卡视频| 国产精品一区二区国产主播| 国产欧美日韩一区二区视频在线| 亚洲精品自产拍在线观看APP| 国产精品乱偷免费视频| 在线免费a视频| 成年人免费国产视频| 国产成人精品一区二区秒拍1o | 国产成人精品综合| 97精品久久久大香线焦| 婷婷五月在线视频| 国产精品久久久久久影院| 国产成人精品一区二区三区| 综合色区亚洲熟妇在线| 亚洲欧美日韩精品专区| 国产在线一二三区| 成人韩免费网站| 久视频免费精品6| 自偷自拍三级全三级视频| 国产精品人成在线播放| 伊人中文网| 美女内射视频WWW网站午夜 | 欧美精品色视频| 第九色区aⅴ天堂久久香| 国产午夜福利亚洲第一| 美女一区二区在线观看| 久久久久久久蜜桃| 成人在线观看不卡| 亚州AV秘 一区二区三区| 一级毛片免费高清视频| 亚洲中文字幕精品| 伊人91在线| 成人午夜网址| 久久精品女人天堂aaa| 九月婷婷亚洲综合在线| 久久夜色精品国产嚕嚕亚洲av| 国产高清在线观看91精品| 国产综合在线观看视频| 国产一区二区免费播放| 中文字幕 91| 国产真实乱子伦精品视手机观看 | 国产在线视频自拍| 色综合婷婷| 国产成人麻豆精品| 久久一日本道色综合久久| 性69交片免费看| 亚洲区欧美区| 亚洲高清中文字幕| 国产成人资源| 在线播放国产99re| 四虎影视8848永久精品| 亚洲国产精品无码AV| 日本草草视频在线观看| 国产精品亚洲αv天堂无码| 国产精品无码翘臀在线看纯欲| 久热re国产手机在线观看| 69av免费视频| 激情午夜婷婷| 色妺妺在线视频喷水| 日韩经典精品无码一区二区| 国产精品综合色区在线观看|