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

受限空間內超聲速混合層生長特性

2017-11-08 01:56:41魏祥庚曹東剛吳繼平
哈爾濱工業大學學報 2017年10期
關鍵詞:生長

魏祥庚,曹東剛,秦 飛,吳繼平

(1.高超聲速沖壓發動機技術重點實驗室(國防科學技術大學),長沙 410073;2.燃燒、熱結構與內流場重點實驗室(西北工業大學),西安 710072)

受限空間內超聲速混合層生長特性

魏祥庚1,2,曹東剛2,秦 飛2,吳繼平1

(1.高超聲速沖壓發動機技術重點實驗室(國防科學技術大學),長沙 410073;2.燃燒、熱結構與內流場重點實驗室(西北工業大學),西安 710072)

為獲得受限空間內激波作用下的超聲速混合層生長規律,以支板噴射超燃沖壓發動機典型流道為研究對象,開展了2.3Ma氫氣射流與2.0Ma空氣來流所形成的超聲速混合層的生長特性研究.基于OpenFOAM計算平臺,采用大渦模擬方法,數值研究了超聲速混合層的流場結構和特征,流場結構和組分分布與實驗結果吻合較好.通過超聲速混合層組分濃度、厚度、可壓縮效應及總壓損失的分析,獲得了超聲速混合層的生長特性.研究結果表明:受限空間內超聲速混合層的生長過程具有4個典型階段,支板末端的膨脹波/激波結構會顯著減低對流馬赫數,從而降低混合層的可壓縮性,促進混合層的生長;激波與混合層的相互作用能夠增強局部湍流強度,獲得渦量增益,加快混合層的生長速率,促進混合效率,但同時會引起較大的總壓損失,降低發動機性能.發動機設計時要綜合考慮波系結構與混合層相互作用帶來的混合增強和總壓損失,實現性能優化.

超聲速混合層;激波;波系結構;受限空間;可壓縮性;大渦模擬

湍流摻混是動力系統中最為常見且極為重要的流動現象,摻混過程的優化組織可以強化燃料與氧氣的混合(非預混燃燒)或是促進燃氣與周圍氣流的熱量交換(預混燃燒),對于縮短燃燒室長度、提高燃燒效率、優化發動機能量管理有著重要的指導意義[1-2].隨著高超聲速飛行技術的發展,高速可壓流動中的湍流摻混成了該領域的研究熱點,也是高超聲速推進系統中的關鍵技術之一.德國宇航院的Papamoschou等[3]以支板噴射超燃沖壓發動機為對象,開展了激波作用下混合層結構的實驗研究,研究指出由于高速條件下氣流的可壓縮性增強,使得混合層的生長速率減弱,湍流強度與雷諾應力也相應減小,混合層的穩定性增強,并建議采用對流馬赫數Mc作為研究參數來表征混合層的可壓縮性.Furby等[4]以HyShot為研究對象,分析了混合層發展過程中的旋渦類型及演化過程.Gutmark等[5]分析了自由空間內可壓縮混合層生長速率與對流馬赫數之間的關系,指出當對流馬赫數增大時,可壓縮混合層的生長速率趨于不可壓縮混合層生長速率的20%.文獻[6-8]對比了自由空間中不同對流馬赫數條件下混合層的生長過程,研究結果均表明隨著對流馬赫數的增加,當混合層中的大尺度結構相對于當地流動的馬赫數大于1時,當地流動受到大尺度結構的阻礙而產生激波,混合層中會出現小激波結構.文獻[9-11]研究結果表明當激波打到混合層上之后能夠強化流體微團的旋轉運動,增強局部的組分輸運,進而促進混合過程.Zhang等[12]以自由空間內的可壓混合層為研究對象,討論了激波對混合層發展過程的影響,研究表明在激波/混合層的作用點附近,由于激波的壓縮作用會使得混合層厚度降低,但在作用點下游,由于渦量的增加使得混合層厚度顯著增大.

對于高超聲速飛行器而言,由于燃燒室進口氣流為超聲速,燃料與空氣需要在有限的空間和極短的時間內完成高效摻混,進而組織燃燒.與自由空間內高速可壓流動中的湍流摻混相比,高超聲速推進系統中的混合層發展不僅受其自身的可壓縮性影響,同時也受到流道幾何結構的約束,其生長特性更為復雜.此外,由于沖壓流道內往往具有復雜的波系結構,波系/混合層的相互作用也是混合層生長的重要影響因素.為了研究受限空間內超聲速混合層的生長特性,本文以支板噴射超燃沖壓發動機典型流道為研究對象,基于計算軟件(open field operation and manipulation, OpenFOAM),采用大渦模擬(large eddy simulation, LES)對2.3Ma氫氣射流與2.0Ma空氣來流所形成的超聲速混合層進行了研究,獲得了混合層的流動結構和流場特征,分析了受限空間內超聲速混合層的生長特性,有助于加深對受限空間內超聲速混合層特征的認識和對復雜波系作用下湍流摻混過程的理解.

1 物理模型和計算方法

1.1 物理模型

本文以支板噴射超燃沖壓發動機為研究對象,其流道構型如圖1所示[13],定義空氣設備噴管喉部為x=0,支板前端位于x=95 mm處,支板末端位于x=175 mm處.冷態空氣經噴管加速后以2.0Ma的速度進入燃燒室內,中心支板對來流空氣進行適當的壓縮,并為燃料噴注提供載體,氫氣通過支板內置的噴管加速后以2.3Ma的速度進入燃燒室,與空氣來流形成超聲速混合層.設備噴管的喉部高度為20.12 mm,燃燒室入口截面為35.4 mm×40.0 mm,中心支板全長80 mm,支板末端高度為2 mm,氫氣噴嘴出口高度為0.6 mm.氫氣流量和空氣流量分別為2.4、740 g/s.

圖1 流道構型示意

1.2 計算方法

大渦模擬結合了直接數值模擬的準確性和雷諾時均模式的快速性,計算精度較高,是湍流摻混燃燒數值模擬的有效方法,LES作為研究燃燒流動的有效手段已經得到了廣泛的認可,并且逐漸被用于研究各種類型發動機的燃燒流動細節.描述湍流混合過程的基本控制方程為可壓Navier-Stokes方程組,對連續方程、動量方程、能量方程、組分方程采用Favre平均并進行濾波運算后可以得到描述湍流摻混的LES控制方程組[14-15]:

式中:“-”為空間濾波;“~”為采用Favre平均后大于濾波尺度的求解參數項;上標“sgs”為小于濾波尺度的亞格子模型求解參數項;ρ為密度;uj為速度;p為壓力;δij為單位張量函數;τij為應力項;E為能量;Ym為組分m的質量分數;q為熱通量;θ為擴散通量.在LES的求解過程中,對于大尺度旋渦進行直接求解,對于比濾波尺度小的流動結構通過亞格子模型進行模擬.本文以OpenFOAM為計算平臺,使用對于可壓縮流動具有較好適用性的rhoReactingFoam求解器進行計算,采用Menon等[16]使用的oneEqEddy亞格子模型對亞格子項進行封閉:

空間離散采用二階精度的TVD格式,考慮背壓對于流動的影響,空氣來流與燃料射流均從噴管喉部開始進行計算,計算入口的邊界參數(壓力、速度、溫度、組分)全部給定,見表1.流道出口仍以超聲速流動為主,根據特征線理論,超聲速出口參數可以根據場內信息進行外插計算得到.時間離散采用二階隱式Crank-Nicholson格式,時間步長選取為1×10-8s,計算的最大CFL數不超過0.3.網格劃分采用結構化網格,對氫氣射流出口和燃燒室中心區域的網格進行了加密以便更好地捕捉混合層中的流動結構.為了更為真實地模擬受限空間內的流動狀況,采用無滑移絕熱壁面,結合大渦模擬對壁面處網格的要求,通過附面層對壁面網格質量進行了提升,壁面處第1層網格的厚度約為0.01 mm,從而使得y+≤1,并使用文獻[13]中通過實驗和計算所提取出的二維平面特征進行計算,網格總數為74萬,燃燒室計算區域和網格劃分如圖2所示.

表1 空氣來流和氫氣射流的入口邊界條件

圖2 計算網格

2 結果與分析

在支板噴射超燃沖壓發動機中,空氣來流經中心支板適當壓縮后仍以超聲速進入燃燒室,燃料經過支板內置的噴管膨脹加速后射入燃燒室,空氣來流與氫氣射流形成超聲速混合層,進而主導燃燒室內的湍流摻混過程.圖3給出了實驗紋影圖[13],計算所得的數值紋影結果以及5 000個時間步統計平均紋影結果,可以看出計算捕捉到了與實驗觀察結果相一致的主要流場特征:由于中心支板的擠壓作用,超聲速空氣來流進入流道后會在支板前緣形成斜激波,由于受到沖壓流道自身幾何構型的約束,激波在燃燒室壁面進行反射并向下游傳播;中心支板的菱形結構使得支板前部形成壓縮型面,后部形成擴張型面,因此會在轉折點(x=135 mm)處形成膨脹波;氣流進入燃燒室后,由于流道的突然擴展,會在支板后緣形成一系列的膨脹波系,氣流經過支板尾部的膨脹波系后向噴流中心偏轉,進而在噴嘴附近又形成兩道斜激波;激波在向下游的傳播過程中,不僅受流道約束會在固壁面進行反射,而且會與混合層相互作用,整個流場表現出以混合層為主導且伴隨多種波系結構的復雜特征.

圖3 實驗紋影圖、數值紋影圖與平均紋影圖對比

Fig.3 Comparisons of experimental schlieren photograph, averaged numerical schlieren image, and instantaneous numerical schlieren image

圖4給出了沿流向4個不同位置處H2的摩爾分數分布曲線.

圖4 不同位置處H2的摩爾分數分布實驗與計算對比

Fig.4 Calculated and measured H2profiles at different locations

從圖4中可以看出計算曲線與實驗數據吻合良好,驗證了計算方法的合理性與可行性.流道內混合層的發展主導著混合過程的進行,湍流摻混首先發生在氫氣射流與空氣來流相接觸的剪切層中,伴隨著剪切層的生長過程,氫氣射流與空氣來流實現質量、動量和能量的交換.在初始階段(x=178 mm)各組分分布區域較為固定且燃料區與氧化劑區域之間存在著巨大的濃度梯度,燃料經支板噴射進入燃燒室后與空氣來流之間形成速度不連續的間斷面,由于間斷面的不穩定性,剪切層產生波動進而發展成為旋渦,強化湍流.伴隨著流動的發展,氫氣射流與空氣來流相互卷吸,射流與周圍流體的摻混自邊緣逐漸向中心發展,混合層逐漸向兩側擴展,混合層厚度逐漸增大,當混合層充分發展以后(x=325 mm),中心射流與周圍流體之間大的參數梯度逐漸被抹平,主要組分的濃度曲線逐漸變得平穩.

為了進一步描述復雜波系作用下受限空間內超聲速混合層的生長特性,定義混合層厚度δ作為表征參數[17]為

δ(x)=R1(x)φ=0.5-R2(x)φ=1.5,

其含義為局部當量比φ=0.5(氫氣質量分數YH2=1.431%)的等值線到局部當量比φ=1.5(氫氣質量分數YH2=4.25%)的等值線的距離.圖5給出了混合層厚度沿流道方向的分布曲線,從圖5中可以看出,混合層的生長過程可以分為4個階段:回流生長區(175 mm

圖5 混合層厚度沿流道方向的分布曲線

圖6 x方向速度云圖

圖7給出了流場中通過壓力梯度表征的波系結構,圖8給出了流場中相應的渦量分布,當激波打到混合層上以后(如x=223 mm),由于壓縮作用會使得混合層的厚度有所減低,與此同時入射激波會迫使外側的空氣來流進入內側的燃料射流中,而且激波引起的局部逆壓梯度會強化剪切層上、下兩側流體微團的旋轉運動,誘導形成渦量增益,增強兩股氣流在垂直方向的輸運過程.入射激波引起的劇烈擾動能夠激勵剪切層K-H不穩定性,會增大局部的湍流強度,獲得渦量增益,增大混合層中旋渦結構的卷起與合并速度,促進混合層上、下兩側燃料與空氣的摻混,使得臨近作用點下游的混合層生長速率增大.需要指出的是激波/混合層的相互作用受到激波強度與混合層內外兩側氣體參數梯度的共同影響,隨著流動向下游的發展,激波強度逐漸減弱,混合層內外兩側氣體的參數梯度也逐漸減小,當激波強度較大且氣體參數梯度較小時,激波會穿透混合層,表現出折射特征,對混合層的壓縮作用也較為顯著,但整體上還是可以促進混合過程的進行,這與文獻[15]中觀測的趨勢一致.在快速生長區之后,混合層經過充分發展,進入飽和生長區,混合層的生長速率減慢,其厚度趨于穩定,標志著混合過程的基本結束.

圖7 壓強梯度云圖

圖8 渦量分布云圖

根據文獻[15]的分析,該工況中氫氣射流與空氣來流的對流Ma>1[13],這就意味著可壓縮效應比較顯著,混合層的生長速度應該比較緩慢,實現充分摻混需要更長的距離.然而,從圖5中可以看出,初始生長區中混合層的生長速率的確比較緩慢,但是經過很短的距離后混合層就開始迅速發展,進入快速生長區,混合層的生長速率顯著增大.為了分析這一現象,圖9、10分別給出了流場中時均化的馬赫數和聲速分布云圖.氫氣射流通過燃料噴嘴以2.3Ma的速度進入燃燒室,與2.0Ma的空氣來流進行摻混,從圖9可以看出,在支板末端形成了亞聲速的回流區,有利于燃料與空氣的摻混.隨著流動向下游的發展,亞聲速的回流區逐漸減小,兩股氣流通過界面處的剪切作用進行動量交換,混合區域內呈現出以跨聲速流動為主的流場特征.

圖9 時均化馬赫數云圖

圖10 時均化聲速云圖

根據聲速的定義為

可知氣體中聲速c的大小直接代表了氣體的可壓縮性的大小,氣體中的聲速越大,則氣體的可壓縮性就越小.從圖10中可以看出,氫氣射流與空氣來流通過混合層實現摻混,混合核心區的聲速相對于氫氣射流的聲速有所降低,相對于空氣來流的聲速卻有所提高,這表明在剪切層內混合氣體的可壓縮性發生了顯著地變化.為了量化可壓縮效應,圖11給出了文獻[3]提出的對流馬赫數Mc和Slessor等[18]提出衡量混合層可壓縮效應的參數Πc沿流道方向的分布.Mc和Πc的定義式分別為:

式中:U為氣體速度;c為氣體中的聲速;γ為氣體比熱比.從圖11中可以看出,在噴射截面x=175 mm處,氫氣射流(Ma=2.3)與空氣來流(Ma=2.0)的對流馬赫數約為1.3,Slessor數約為3.75,說明兩股超聲速氣流所形成的混合層在開始階段具有很強的可壓縮性,經過很短的距離后,對流馬赫數與Slessor數迅速降低,當初始生長區結束后,Mc降為0.3,降幅約為77%,Πc降為0.4,降幅約為88%,這都表明混合層的可壓縮效應顯著降低,混合層的生長速率大幅度提高.結合圖3中的流場結構可以看出,由于流道擴張在支板末端形成膨脹波系,氣流方向偏離流道方向,受到流道幾何構型和兩股氣流所形成的氣動壓縮面的約束,在膨脹波下游會形成再附激波使得氣流的整體流動方向與流道方向一致,支板末端的這種膨脹波/激波結構會使得混合層的可壓縮效應大幅度減小,從而改變混合層的生長特性.這與通常所說的高對流馬赫數混合層的生長速率較小并不矛盾,只是與自由空間相比,受限空間內混合層的對流馬赫數成為受復雜波系結構影響的空間函數,初始對流馬赫數并不能全面描述混合層的生長特性.

圖11 Mc和Πc沿流道的分布

需要強調的是,雖然受限空間內的波系結構一方面能夠降低混合層的可壓縮性,整體上增大混合層的生長速率,另一方面激勵K-H不穩定波的發展,強化局部湍流強度,獲得渦量增益,使得局部混合層的生長速率進一步增大,從而促進湍流摻混,提高混合效率,但沖壓流道中的復雜波系結構不可避免地會引起較大的流動損失,降低發動機的整體性能.以噴射截面x=175 mm處的總壓為參考,圖12給出了總壓損失沿流道方向的分布曲線.

圖12 沿流道方向的總壓損失分布曲線

從圖12中可以看出沿著流道方向總壓損失不斷增大,在流道上游,激波強度較大,引起較大的總壓損失,總壓損失的上升速率較快,隨著流動向下游的發展以及波系向下游的傳播,激波強度逐漸減小,總壓損失的上升速率趨于緩慢,在混合層生長階段,總壓損失接近9%.因此在實際的發動機設計過程中,需要對激波提高的混合效率和引起的總壓損失進行綜合考慮,實現優化設計.

3 結 論

1)受限空間內混合層的對流馬赫數是受復雜波系結構影響的空間函數,支板末端所形成的膨脹波/激波結構會使對流馬赫數迅速降低,顯著減弱可壓縮效應,從而提高混合層的生長速率,使得混合層經過較短距離后進入快速生長區.

2)激波強度與混合層內外兩側氣流的參數梯度共同影響著兩者之間的相互作用,由于激波的壓縮效應,在作用點附近混合層的厚度會略微減小,但在作用點下游,激波會強化湍流,獲得渦量增益,促進混合層的生長.

3)受限空間內的復雜波系結構在促進湍流摻混,提高混合效率方面具有積極的作用,但同時會引起較大的總壓損失,在實際發動機設計過程中需要對其進行綜合考慮,實現優化設計.

[1] KUMAR A, BUSHNELL D H, HUSSAINI M Y. Mixing augmentation techniques for hypervelocity scramjets[J]. Journal of Propulsion and Power, 1989, 5(5): 514-522.DOI:10.2514/3.23184.

[2] 張漫.RBCC燃燒室湍流噴霧燃燒及其穩定機制研究[D]. 西安: 西北工業大學, 2010.

[3] PAPAMOSCHOU D, ROSHKO A.The compressible turbulent shear layer: an experimental study[J]. Journal of Fluid Mechanics, 1988, 197: 453-477. DOI:10.1017/S0022112088003325.

[4] FUREBY C, CHAPUIS M, FEDINA E, et al. CFD analysis of the HyShotⅡscramjet combustor[J]. Proceeding of the Combustion Institute, 2011, 33(2): 2399-2405.DOI:10.1016/j.proci.2010.07.055.

[5] GUTMARK E J, SCHADOW K C, YU K H. Mixing enhancement in supersonic free shear flows[J]. Annual Review of Fluid Mechanics, 1995, 27: 375-417.DOI:10.1146/annurev.fl.27.010195.002111.

[6] LI Qibing, FU Song. Numerical simulation of high-speed planar mixing layer[J]. Computers & Fluids, 2003, 32(10):1357-1377.DOI:10.1016/S0045-7930(02)00104-4.

[7] SANDHAM N D, REYNOLDS W C. Three-dimensional simulations of large eddies in the compressible mixing layer[J]. Journal of Fluid Mechanics, 1991, 224: 133-158.DOI:10.1017/S0022112091001684.

[8] VREMAN A W, SANDHAM N D, LUO K H. Compressible mixing layer growth rate and turbulence characteristics[J]. Journal of Fluid Mechanics, 1996, 320(1): 235-258.DOI:10.1017/S0022112096007525.

[9] MARBLE F E, HENDRICKS G J, ZUKOSKI E E. Progress toward shock enhancement of supersonic combustion process[C]//AIAA/SAE/ASME/ASEE 23rd Joint Propulsion Conference. San Diego, California: AIAA, 1987: 932-950.DOI: 10.2514/6.1987-1880.

[10]LU P J, WU K C. On the shock wave enhancement of confined supersonic mixing flows[J].Physics of Fluids A, 1991, 3(12): 3046-3062.DOI:10.1063/1.857849.

[11]HUH H, DRISCOLL J F. Shock-wave-enhancement of the mixing and the stability limits of supersonic hydrogen-air jet flames[J]. PSymposium (International) on Combustion, 1996, 26(2): 2933-2939.DOI:10.1016/S0082-0784(96)80135-4.

[12]ZHANG Yunlong, WANG Bing, ZHANG Huiqiang, et al. Mixing enhancement of compressible planar mixing layer impinged by oblique shock waves[J]. Journal of Propulsion and Power, 2015, 31(1): 156-169.DOI: 10.2514/1.B35423.

[13]GERLINGER P, BRUGGEMANN D. Numerical investigation of hydrogen strut injections into supersonic airflows[J]. Journal of Propulsion and Power, 2000, 16(1): 22-28.DOI: 10.2514/2.5559.

[14]張兆順, 崔桂香, 許春曉. 湍流大渦數值模擬的理論和應用[M]. 北京: 清華大學出版社, 2008.

[15]王振國, 孫明波. 超聲速湍流流動、燃燒的建模與大渦模擬[M]. 北京: 科學出版社, 2013.

[16]MENON S, YEUNG P K, KIM W W. Effect of subgrid models on the computed interscale energy transfer in isotropic turbulence[J]. Computers & Fluids, 1996, 25(2): 165-180.DOI:10.1016/0045-7930(95)00036-4.

[17]KIM J H, YOON Y, JEUNG I S. Numerical study of mixing enhancement by shock waves in model scramjet engine[J]. AIAA Journal, 2003, 41(6): 1074-1080.DOI: 10.2514/2.2047.

[18]SLESSOR M D, ZHUANG M, DIMOTAKIS P E. Turbulent shear-layer mixing: growth-rate compressibility scaling [J]. Journal of Fluid Mechanics, 2000, 414(1): 35-45.DOI:10.1017/S0022112099006977.

Studyonthesupersonicmixinglayergrowthinconfinedspace

WEI Xianggeng1,2, CAO Donggang2, QIN Fei2, WU Jiping1

(1.Science and Technology on Scramjet Laboratory(National University of Defense Technology), Changsha 410073, China; 2.Science and Technology on Combustion, Internal Flow and Thermal-Structure Laboratory (Northwestern Polytechnical University), Xi’an, 710072, China)

The supersonic mixing layer formed by a planar thin hydrogen jet at 2.3Maand a 2.0Masurrounding airflow in a scramjet engine model is studied in order to investigate the growth characteristics with consideration of shocks. Flow field structure and features of the supersonic mixing layer are achieved by using large eddy simulation method with the OpenFOAM software. Reasonable agreements are obtained between calculation and experiment in terms of flow field structure and component distribution. The component concentration, the thickness, the compressibility effect and the total pressure loss are analyzed. Results show that four developing regions can be observed for the growth of the mixing layer. The expansion-fan/shock-wave pattern at the injector exit makes the convective mach number decrease dramatically, leading to a reduction in compressibility effects and a contribution to the development of the mixing layer. The interaction of shock/mixing layer results in local amplification of turbulence and gain of vorticity, which is beneficial to the supersonic mixing. However, the increasement in total pressure loss is unavoidable in the presence of shocks because they can bring performance losses of the scramjet. Thus a tradeoff between the enhanced mixing efficiency and the decreased total pressure recovery should be considered in the scramjet optimization design.

supersonic mixing layer; shock wave; shockwave series; confined space; compressibility; large eddy simulation

10.11918/j.issn.0367-6234.201607066

V431

A

0367-6234(2017)10-0072-06

2016-07-16

高超聲速沖壓發動機技術重點實驗室開放基金(20110302021)

魏祥庚(1979—),男,博士,副教授

魏祥庚,realysnow@nwpu.edu.cn

(編輯張 紅)

猜你喜歡
生長
野蠻生長
碗蓮生長記
小讀者(2021年2期)2021-03-29 05:03:48
生長的樹
少兒美術(2020年3期)2020-12-06 07:32:54
自由生長的家
現代裝飾(2020年11期)2020-11-27 01:47:48
美是不斷生長的
快速生長劑
共享出行不再“野蠻生長”
生長在哪里的啟示
華人時刊(2019年13期)2019-11-17 14:59:54
野蠻生長
NBA特刊(2018年21期)2018-11-24 02:48:04
生長
文苑(2018年22期)2018-11-19 02:54:14
主站蜘蛛池模板: 99爱在线| 特级欧美视频aaaaaa| 国产精品亚洲一区二区三区在线观看| 91精品国产91久久久久久三级| 国产主播一区二区三区| 日韩欧美视频第一区在线观看| 波多野结衣在线一区二区| 欧美日韩成人| 色婷婷亚洲综合五月| 亚洲精品欧美日韩在线| 97久久免费视频| 国产91视频观看| 国产欧美日韩精品综合在线| 91成人在线免费视频| 亚洲水蜜桃久久综合网站| 成人久久精品一区二区三区| 亚洲综合天堂网| a级毛片免费网站| 五月天久久综合国产一区二区| 国产精品性| 日本午夜影院| 亚洲第一视频免费在线| 国产成人久久综合一区| 国产无码高清视频不卡| jizz亚洲高清在线观看| 热久久综合这里只有精品电影| 色成人亚洲| 亚洲一级毛片在线观| 香蕉网久久| 日韩精品一区二区三区大桥未久| 国产91视频免费| 最新亚洲av女人的天堂| 美女潮喷出白浆在线观看视频| 亚洲色无码专线精品观看| a级免费视频| 亚洲黄色成人| 天天色综网| 成人国产免费| 亚洲最新地址| 又黄又湿又爽的视频| 国产精品成人AⅤ在线一二三四| 美女高潮全身流白浆福利区| 成人毛片在线播放| 幺女国产一级毛片| 国产美女91呻吟求| 91色在线观看| 欧美亚洲国产日韩电影在线| 成人在线观看不卡| 日韩欧美色综合| 99久久国产综合精品2023| 九九视频免费在线观看| 国产丝袜丝视频在线观看| 亚洲国产日韩在线观看| 国产欧美在线观看精品一区污| 67194亚洲无码| 亚洲一区二区三区麻豆| 亚洲欧洲自拍拍偷午夜色| 欧洲一区二区三区无码| 黄色三级网站免费| 欧美日韩理论| 波多野结衣中文字幕久久| 国产高清免费午夜在线视频| 2021国产精品自产拍在线| 99热在线只有精品| 日本人真淫视频一区二区三区| 香蕉久久国产超碰青草| 丝袜久久剧情精品国产| 美女免费黄网站| 米奇精品一区二区三区| 日韩福利视频导航| 日本成人精品视频| 亚洲日韩高清无码| 久久精品日日躁夜夜躁欧美| 啪啪免费视频一区二区| 欧亚日韩Av| 亚洲91在线精品| 免费人成视频在线观看网站| 在线观看免费黄色网址| 黄色网站不卡无码| 国产精品永久久久久| 波多野结衣一区二区三视频 | 欧美精品黑人粗大|