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

雙粗糙元對高超音速邊界層感受性的影響

2024-01-01 08:13:12史明方徐立丹唐小軍王振清
哈爾濱工程大學學報 2023年10期
關鍵詞:模態區域影響

史明方,徐立丹,唐小軍,王振清

(1.內蒙古科技大學 土木工程學院,內蒙古 包頭 014000;2.哈爾濱工程大學 航天與建筑工程學院,黑龍江 哈爾濱 150001;3.中國空間技術研究院 北京衛星制造廠有限公司,北京 100094)

高超音速武器以其優越的突防性能和損傷能力引起各國的重點關注[1]。在高超音速武器的設計過程中,高超音速邊界層流動狀態會直接影響高超音速武器的性能[2],因為邊界層層流轉捩湍流會導致邊界層內摩擦阻力和壁面加熱顯著增加。鑒于此,高超音速邊界層的轉捩控制是高超音速武器設計的關鍵技術之一。感受性是轉捩的第1 個階段,定義了邊界層內擾動波的初始模態和幅值等參數。感受性階段受到很多因素的影響,簡單的可以分為來流擾動和壁面擾動[3]。其中來流擾動對感受性的影響研究較早,并已經取得了一系列的成果。

壁面粗糙元作為壁面擾動的重要組成部分,成為近些年的研究熱點。在高超音速飛行器表面不可避免地存在壁面粗糙元,而壁面粗糙元又會對高超音速流場及邊界層造成顯著的影響,因此開展壁面粗糙元對邊界層感受性作用機理的研究是十分必要的。董昊[4]總結了壁面粗糙元對高超音速邊界層轉捩的影響,指出壁面粗糙元對高超音速邊界層感受性的作用機理十分復雜,包括壁面粗糙元和邊界層內擾動波的相互作用,涉及瞬態增長和尾跡等現象。Saikishan 等[5]研究了壁面粗糙元造成的壓力梯度對邊界層轉捩的影響。Clark 等[6]利用開源軟件SU2 優化了渦流發生器的形狀,以達到抑制邊界層內Mack 模態演化的目的。壁面粗糙元可以依據相對光滑壁面的位置進行分類,在本文中,將高于光滑壁面的粗糙元稱為K型粗糙元,相應的,將低于光滑壁面的粗糙元稱為D 型粗糙元,K 型粗糙元的研究相對較早。段志偉等[7]分析了壁面粗糙元形狀對轉捩模式的影響。Ryota 等[8]結合fhf 直接數值模擬(direct numerical simulation,DNS)和非線性拋物化穩定性方程(nolinear parabolizd stability equations,NPSE)方法分析了三維圓柱狀粗糙元高度對邊界層轉捩的影響,并討論了邊界層轉捩時壁面粗糙元的臨界高度。Vignes 等[9]分析了渦流發生器對圓柱形粗糙元周圍高超音速流場氣動特性的影響,發現渦流發生器有助于降低圓柱形粗糙元造成的氣動阻力。D型粗糙元(主要集中在矩形腔結構)對高超音速邊界層感受性的影響成為研究熱點。李闖等[10]采用直接數值模擬方法討論了空腔結構對6.0 Ma 高超音速邊界層內擾動波演化的作用機理,發現空腔深度的增加會改變邊界層內第2 模態的演化機理。劉俊等[11]系統地研究了空腔結構內波系的傳播規律,討論了壓力振蕩反饋模型的異同。李強等[12]利用激波風洞分析了展向凹槽對高超音速平板邊界層轉捩的影響,發現凹槽會促進轉捩,并且隨著凹槽深度和寬度的增加而增強。上述學者所分析的壁面粗糙元均為獨立粗糙元,與之相對應的分布式粗糙元的研究在獨立粗糙元的基礎上也在逐步開展。董昊等[13]利用數值模擬和油膜干涉技術分析了壁面粗糙元高度和間距對轉捩位置的影響,發現高度的影響要大于間距的影響。Irimpan 等[14]利用分布式粗糙元來表征燒蝕后的熱防護系統,發現分布式粗糙元可以降低壁面加熱率17%~42%。Zhou 等[15]指出波浪形壁面可以有效地抑制邊界層內第2模態不穩定性,并且抑制效應隨波浪形壁面的延長、加深和向下游移動而增強。Tong等[16-17]研究了正弦狀壁面粗糙元對高超音速邊界層轉捩的影響。

目前,獨立粗糙元和分布式粗糙元的研究都已取得了一定成果,但雙粗糙元對高超音速邊界層感受性影響的研究相對較少,并且大部分研究集中在分析雙粗糙元間距的影響,缺乏雙粗糙元組合類型影響的研究。本文對高超音速鈍錐繞流流場進行直接數值模擬,討論了雙粗糙元組合類型對高超音速流場及邊界層的影響;向流場中引入有效振幅脈沖慢聲波,借助快速傅里葉變換方法分析了雙粗糙元組合類型對脈沖慢聲波作用下邊界層內擾動波傳播過程的影響,以揭示雙粗糙元組合類型對高超音速邊界層感受性的作用機理。

1 控制方程及數值方法

為了分析雙粗糙元對脈沖有限振幅自由來流擾動作用下高超音速流場及邊界層內流動及邊界層的影響,本文采用直接數值模擬方法對高超音速流場進行直接數值模擬,揭示雙粗糙元的組合類型對高超音速邊界層感受性的影響。

對高超音速流場的數值模擬,控制方程為二維N-S方程:

式中:U為解向量;F和G為對流項;Eυ與Fυ為黏性項。解向量的具體表達式為:

對流項的表達式為:

黏性項的表達式為:

式中:ρ、u、v、T和p分別為密度、x方向速度、y方向速度、溫度和壓強;導熱系數k可以通過假設普朗特數為常數獲得;E為單位流體的總能,可以表示為單位流體的內能e與動能的和,具體為:

對于量熱氣體,氣體內能e可計算為:

式中:比熱比γ=1.4;τ為切應力,一般將空氣看作牛頓流體:

式中μ為黏性系數,可以通過Sutherland公式獲得:

式中:T0和μ0分別為參考溫度和黏性系數;Ts=110.4 K。

本文考慮空氣為完全氣體,補充狀態方程為:

式中:R為通用氣體常數,R=286.94 J/(kg·K)。此外,本文將笛卡爾坐標系下控制方程進行坐標變換,轉換到曲線坐標系下,以便于對復雜外形下的流場進行數值求解。

為了分析雙壁面粗糙元組合形式對高超音速邊界層感受性的影響,本文以自由來流擾動作用下高超音速鈍錐繞流流場為研究對象,鈍錐鼻部半徑r=1 mm,半錐角θ=5°。本節以K型粗糙元(上游)和D型粗糙元(下游)的組合為例,給出了本文的計算模型及網格,計算網格分布為801×251,計算模型如圖1 所示。K 型粗糙元(上游)和D 型粗糙元(下游)的雙壁面粗糙元組合形式表示為Kf Db,下標f和b分別代表上游粗糙元和下游粗糙元,其余3 種組合形式以此類推。壁面粗糙元形狀方程(以弧線ABCDE為例)為:

圖1 計算模型及網格Fig.1 Computing model and grid

點A和點E分別為壁面粗糙元的起點和終點;點B和點D是壁面粗糙元形狀的拐點(xⅠf或xⅠb和xⅡf或xⅡb),分別稱為第1拐點(xⅠf或xⅠb)和第2拐點(xⅡf或xⅡb);點C為壁面粗糙元中點(xcf或xcb)。需要說明的是,K型粗糙元和D型粗糙元在同一位置處(壁面粗糙元中點相同)關于光滑壁面對稱。計算模型的幾何參數在表1中給出。本文使用的變量將分別利用頭部半徑r、來流參數ρ∞、u∞、T∞、p∞,u∞/r,r/u∞、k∞和μ∞對長度、密度、速度、溫度、壓強、頻率、時間、導熱系數和黏性系數進行無量綱化處理,自由來流條件在表2中給出。

表1 壁面粗糙元幾何參數Table 1 Geometric parameters of roughness element mm

表2 自由來流條件Table 2 Conditions of freestream

自由來流參數,如擾動類型、持續時間等,會顯著地改變高超音速邊界層感受機理[18-19]。較大幅值的自由來流擾動與邊界層的作用更加強烈,導致感受性階段更復雜,甚至改變轉捩路徑。擾動波的作用時間也會顯著地影響邊界層感受性機理,短時間高強度的自由來流擾動又稱為有限振幅脈沖擾動,對高超音速邊界層感受性的影響機理與長時間低強度的自由來流擾動存在很大的差異。本文在分析壁面粗糙元對高超音速邊界層感受性的影響時,向流場中引入脈沖慢聲波擾動,其表達式為:

式中:ε、k和f分別為擾動波的幅值、波數和頻率,ε=6×10-2,k=3.144 6×10-4和f=50π,因此本文自由來流擾動的周期為4。本文中擾動量是指非穩態高超音速流場的瞬時值與穩態流場值的差值,擾動量利用上標“'”表示,如溫度擾動值利用“T'”表示。在分析壁面粗糙元對高超音速邊界層感受性的作用機理時向高超音速流場引入脈沖擾動波,即向流場中引入半個周期(Δt=2)的脈沖慢聲波擾動。

高階精度有限差分方法可以穩定準確地對高超音速流場進行數值模擬。本文基于N-S方程采用S-W(steger-warming)通量矢量分裂方法將對流項分裂為正、負2 個流通量矢量。利用五階加權本質無振蕩(weighted essentially non-scillatory,WENO),格式對上述正、負流通量矢量進行離散;對于粘性項采用六階精度中心差分格式進行離散;時間項利用三步三階TVD 型Runge-Kutta 格式進行推進。本文采用的數值方法已經廣泛使用在高超音速流場及邊界層的數值模擬中[20-21],并已對數值方法和網格獨立性進行過驗證,結果表明計算過程具有良好的穩定性,并可以保證計算結果的準確性。

2 雙粗糙元對穩態高超音速流場的影響

壁面粗糙元會改變高超音速邊界層感受性機理,為揭示雙粗糙元對高超音速邊界層感受性的影響,首先需對雙粗糙元作用下穩態高超音速流場進行分析。圖2 為壁面粗糙元作用下高超音速流場流向速度分布,可以看出壁面粗糙元顯著地改變了流場結構,在流場中形成壓縮波和膨脹波。獨立K型粗糙元導致流場中形成壓縮波-膨脹波-壓縮波波系;而獨立D 型粗糙元導致流場中形成膨脹波-壓縮波-膨脹波波系。雙粗糙元會導致流場中形成更復雜的波系,上游粗糙元造成的壓縮波或膨脹波會與下游粗糙元造成的壓縮波或膨脹波相“融合”,例如KfKb型粗糙元對導致流場中形成壓縮波-膨脹波-壓縮波-膨脹波-壓縮波波系,Kf粗糙元造成的下游壓縮波和Kb粗糙元造成的上游壓縮波發生了融合;相似的,KfDb型粗糙元則在流場中形成了壓縮波-膨脹波-壓縮波-膨脹波波系。需要指出的是K型粗糙元形成的波系結構與Xue等[22]分析在前臺階對高超音速流場的影響是相似的,同時也指出前臺階對高超音速流動的影響局限在一定范圍內。

圖2 穩定高超音速流場流向速度云圖Fig.2 Contours of velocity along x-axis in syeady hypersonic flowfield

壁面粗糙元對高超音速流場結構產生了一定影響,壁面粗糙元也會對高超音速邊界層內流動狀態的影響,圖3 給出了壁面雙粗糙元對穩態高超音速邊界層內壁面壓力的影響。可以看出壁面粗糙元對粗糙元附近壁面壓力分布產生了顯著的影響。以獨立K 型粗糙元為例,壁面壓力在K 型粗糙元前緣區域顯著增加,在第1拐點處達到最大值,之后壁面壓力沿流向方向逐漸減小至第2 拐點,在粗糙元中心恢復到光滑壁面水平。在第2拐點至壁面粗糙元后緣區域,壁面壓力逐漸增加至光滑壁面狀態。壁面壓力的變化情況是由壁面曲率變化導致的,本文所采用的壁面粗糙元由3 次多項式控制,并且壁面粗糙元中點處與光滑壁面平行,因此在壁面曲率增大的區域壁面壓力增加,并在壁面粗糙元中心處恢復到光滑壁面狀態。D型粗糙元對壁面壓力的影響情況相似。此外,壁面粗糙元的影響范圍有限,僅局限在壁面粗糙元附近,并且K 型粗糙元的影響范圍要大于D 型粗糙元。K 型粗糙元的出現導致其上游和下游區域的壓強分布產生了變化,該現象在D 型粗糙元中并不明顯。對于雙粗糙元,雖然壁面粗糙元的影響范圍有限,但本文中下游粗糙元均位于上游粗糙元的影響范圍之內,因此在K 型粗糙元下游的壁面粗糙元區域的壁面壓強均發生了一定程度的減小,而在D 型粗糙元下游的壁面粗糙元區域的壁面壓強均發生了一定程度的增長。DfDb型粗糙元在除兩粗糙元之間的區域內,壁面壓強分布與獨立粗糙元相同。

圖3 高超音速邊界層內壁面壓力比較Fig.3 Comparisons of wall pressure in hypersonic bound‐ary layer

壁面摩擦系數表征了邊界層內的剪切狀態,并且壁面摩擦力會顯著地影響高超音速飛行器的性能。鑒于此,對邊界層內壁面摩擦系數進行分析是有必要的。圖4給出了邊界層內壁面摩擦系數沿流向的分布。計算結果表明壁面粗糙元顯著地改變了邊界層內的流動狀態,在K型粗糙元的中心處,壁面摩擦系數到達極大值。這是因為對于穩態流場來說,邊界層外邊界的速度是相同的,因此壁面摩擦系數與邊界層厚度呈反比。在邊界層內粘性的作用下,K型粗糙元中心處邊界層的厚度最薄,壁面摩擦系數達到極大值。同理,壁面摩擦系數在D型粗糙元的前緣和后遠處達到極大值,在D型粗糙元中心處達到了極小值。當壁面摩擦系數小于零時,說明在區域內發生了回流,產生了流動分離。在本文中,獨立K型粗糙元并未造成流動分離,而獨立D型粗糙元均導致邊界層流動分離。這說明D型粗糙元更容易導致流動的分離。在雙粗糙元中包含D型粗糙元時,在D型粗糙元區域內邊界層都發生了流動分離。

圖4 高超音速邊界層內壁面摩擦系數比較Fig.4 Comparisons of friction coefficient in hypersonic boundary layer

此外,KfKb型雙粗糙元也造成了流動分離,流動分離發生在雙粗糙元之間。因為在Kf粗糙元后緣下游至Kb粗糙元上游區域內,壓力顯著增加,逆壓梯度很大,導致流動在該區域內發生分離。流動分離現象會顯著地改變邊界層內的流動結構,進而影響到邊界層內的壁面摩擦力和熱流密度。獨立粗糙元導致壁面摩擦系數沿流向的減小值要超過增加值,也就是說壁面粗糙元可以起到減小阻力的作用。也有學者得出相似的結論,Brian等[23]指出分布式粗糙元可以在保持邊界層內摩擦阻力基本不變的情況下穩定邊界層內擾動波。因此合理地布置壁面粗糙元可以在一定程度上達到減阻,甚至降熱的目的。結合Zhou 等[15]指出壁面粗糙元可以穩定邊界層的作用,壁面粗糙元的合理布置對高超音速飛行器的設計具有一定的指導意義。因此,對雙粗糙元作用下高超音速邊界層感受性機理的研究是有意義的。

3 雙粗糙元對非穩態高超音速流場的影響

為了分析雙粗糙元對高超音速邊界層感受性的影響,向高超音速流場中引入脈沖慢聲波擾動。圖5 以KfKb型粗糙元為例給出了脈沖擾動波在高超音速流場中的傳播過程。

圖5 脈沖擾動作用下高超音速溫度擾動流場Fig.5 Contours of transient temperature fluctuation at different time

高超音速流場可以分為4 個區域:R1區域-激波外區域;R2區域-激波區域;R3區域-邊界層與激波之間的區域;R4區域-邊界層區域,如圖5(h)所示。脈沖擾動波與弓形激波的作用更加強烈,在脈沖慢聲波經過弓形激波后會在高超音速流場中形成“拖拽波”,拖拽波的形成是因為脈沖慢聲波在經過激波后的速度不同。脈沖慢聲波在與斜激波相互作用時,并未在流場中形成明顯的拖拽波,而是在流場中依次形成低溫-高溫結構。壁面粗糙元影響了脈沖慢聲波在高超音速流場中R3和R4區域的傳播。

第2 節已經指出壁面粗糙元改變了穩態高超音速流場結構,在流場中形成一組新的波系,并導致邊界層發生流動分離。可以看出雙粗糙元造成的壓縮波和膨脹波波系改變了R3區域內擾動波的傳播,在流場中明顯的看出3道壓縮波對流場內擾動波的放大作用,如圖5(k)所示。自由來流擾動在經過KfKb型粗糙元造成的壓縮波時被放大,并在離開這3道壓縮波時恢復到光滑壁面狀態。脈沖慢聲波在經過激波后進入到邊界層內,激發出邊界層內擾動波,導致邊界層內依次形成了低溫斑和高溫斑結構,如圖5(g)中矩形區域所示。邊界層外邊界處形成“繩索”狀高溫斑-低溫斑結構,如圖5(g)中橢圓形區域所示,Ma等[24-26]指出這是邊界層內慢聲波結構的典型結構。邊界層內擾動波在向下游傳播的過程中逐漸耗散,但在邊界層內形成殘余反射波,反射波在邊界層內進行多次反射,最終在邊界層內粘性作用下耗散掉,如圖5(l)所示。

為了更清晰地揭示壁面粗糙元對高超音速邊界層感受性的影響,圖6給出了邊界層內不同位置處壁面壓強擾動值隨時間的變化情況。邊界層內越下游區域越晚感受到邊界層內擾動波,擾動波在向下游傳播過程中,持續時間越長,在x=0.90 處,擾動波持續時間Δt=3.8;而在x=18.89 處,擾動波持續時間Δt=17.9。邊界層內擾動波持續時間都超過了自由來流擾動波的時間(Δt=0.5T=2)。說明邊界層內反射波在邊界層內多次反射后才逐漸耗散消失。此外,擾動波向下游傳播過程中幅值逐漸減小,擾動波波形越復雜。邊界層內擾動波在進入到高超音速流場后呈現“正弦波”分布,在向下游傳播的過程中,在粘性的作用下,擾動波的幅值逐漸減小,波形逐漸變得復雜。邊界層內擾動波在通過第1粗糙元(Kf或Df)時,邊界層內擾動波在通過Kf粗糙元前緣(x=0.90)處被放大,而Df粗糙元前緣處壁面壓強擾動并未受到影響,因為D 型粗糙元的影響范圍較小。在Kf粗糙元中點(x=1.39)處,邊界層內擾動波稍被放大。邊界層擾動波在通過第1粗糙元后緣(x=1.79)處,Kf、Kf Kb和Kf Db粗糙元降低了該處邊界層擾動波的幅值,而Df、Df Db和Df Kb粗糙元放大了該處擾動波的幅值。Kb粗糙元也放大了擾動波的幅值,但Db粗糙元并未對該處產生影響。在兩粗糙元中點(x=1.91)處,Kf和Db粗糙元抑制了邊界層內擾動波的傳播,因此,Kf Db粗糙元對邊界層內擾動波傳播的抑制作用更明顯。同理,Df和Kb粗糙元促進了邊界層內擾動波的傳播,因此,Df Kb粗糙元對邊界層內擾動波傳播的促進作用更明顯。Kf Kb和Df Db雙粗糙元對該點擾動波的傳播基本沒影響,這可以說明同類型的上游粗糙元和下游粗糙元對兩粗糙元中點處的影響是相同的,二者的影響相互抵消。下游粗糙元對邊界層內擾動波的作用機理與上游粗糙元相同,這里不再贅述。壁面粗糙元影響的局限性導致在兩粗糙元下游區域的擾動波的傳播并未受到影響。

圖6 各流向位置處壁面壓強擾動幅值-時間對比Fig.6 Time history trace of wall pressure fluctuation at various streamwise locations

為了更準確地分析高超音速邊界層內擾動波的傳播過程,本文采用快速傅里葉變換方法將邊界層內的時域信號轉換到頻域內,以分析邊界層內不同模態擾動波的演化過程。快速傅里葉變換公式為:

式中:f為自由來流擾動波的頻率;|p'|為任意瞬時壓力擾動;|p'n|和φn分別為局部擾動幅值和相角,下標整數n代表了擾動場的波模態,當n=0代表平均流場變形,n=1代表基頻擾動,n=2代表二階諧頻擾動,以此類推。第2節已經指出壁面粗糙元對高超音速流場的影響具有局限性,因此本節僅給出x=0.90~10.01范圍內邊界層內擾動波的演化情況,如圖7所示。可以看出邊界層內擾動波在向下游傳播過程中擾動波的波形變得越來越復雜,擾動波的幅值在逐漸減小,在x=10.01處稍有增加。此外,邊界層內擾動波的主導模態在逐漸增大,如在x=0.90處,主導模態f=0.033;在x=1.79處,主導模態f=0.028;在x=2.40處,主導模態f=0.025,并且主導模態在邊界層內擾動波所占比例在逐漸減小。也就是說高超音速邊界層內擾動波存在模態競爭現象,主導模態進行了模態遷移。

圖7 邊界層內擾動波頻率振幅Fig.7 Frequency-amplitude diagram of perturbations in boundary layer

壁面粗糙元對邊界層內不同模態擾動波傳播過程的影響相同,即在K 型粗糙元前緣區域和D 型粗糙元的后緣區域內,邊界層內擾動波的幅值被放大;相應的在D型粗糙元前緣區域和K型粗糙元的后緣區域內,邊界層內擾動波的傳播被抑制。值得注意的是Kf獨立粗糙元抑制了其下游區域內擾動的演化。雙粗糙元的組合類型對邊界層內不同模態擾動波的傳播具有一定的影響,Kf Kb型雙粗糙元導致壁面粗糙元下游區域內擾動波的主導模態顯著降低,而且這種抑制效應還要強于Kf獨立粗糙元。下游的Kb粗糙元增強了上游的Kf粗糙元對Kf Kb雙粗糙元下游邊界層內擾動波的抑制效應。這就證明了K型粗糙元對下游的邊界層內擾動波演化的抑制具有協同效應。

4 結論

1)雙壁面粗糙元會顯著地改變穩態高超音速流場及邊界層內的流動狀態,導致流場中形成一組新的波系,并導致邊界層內出現流動分離。如Kf Kb型粗糙元會導致流場中形成壓縮波-膨脹波-壓縮波-膨脹波-壓縮波波系,并在兩粗糙元之間形成流動分離。

2)雙壁面粗糙元導致流場中形成新的波系顯著地影響了高超音速流場中自由來流擾動的傳播過程,使擾動波在壓縮波作用下發生了一定的變形;當自由來流擾動離開壁面粗糙元造成的波系后,新的波系影響也隨之消失。

3)雙壁面粗糙元也改變了邊界層內擾動波的演化過程,在K 型粗糙元前緣區域和D 型粗糙元的后緣區域內,邊界層內擾動波的幅值被放大;在D型粗糙元前緣區域和K 型粗糙元的后緣區域內,邊界層內擾動波的傳播被抑制。壁面粗糙元對邊界層內擾動波的作用也局限在一定范圍內。

4)邊界層內擾動波在向下游傳播過程中主導模態在逐漸減小,并存在模態競爭現象。雙粗糙元的組合類型對邊界層內不同模態擾動波的傳播具有一定的影響。其中,KfKb型雙粗糙元顯著抑制了其下游區域內主導模態擾動波的演化,并且2個K型粗糙元對其下游邊界層內擾動波演化的抑制具有協同效應。

猜你喜歡
模態區域影響
是什么影響了滑動摩擦力的大小
哪些顧慮影響擔當?
當代陜西(2021年2期)2021-03-29 07:41:24
擴鏈劑聯用對PETG擴鏈反應與流變性能的影響
中國塑料(2016年3期)2016-06-15 20:30:00
關于四色猜想
分區域
國內多模態教學研究回顧與展望
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
區域
民生周刊(2012年10期)2012-10-14 09:06:46
主站蜘蛛池模板: 色综合综合网| 国产91精品调教在线播放| 色香蕉影院| 99视频免费观看| 国产网站免费看| 亚洲天堂网在线观看视频| 青草视频久久| 999福利激情视频| 超碰精品无码一区二区| 国产精品自在自线免费观看| 狠狠色噜噜狠狠狠狠色综合久 | 国产日本一区二区三区| 国产福利小视频高清在线观看| 天天做天天爱天天爽综合区| 亚洲综合第一页| 日韩大乳视频中文字幕 | 极品性荡少妇一区二区色欲| 538国产视频| 五月婷婷导航| 婷婷亚洲综合五月天在线| 国产美女一级毛片| 青草精品视频| 亚洲成a人片77777在线播放| 亚洲无码电影| 91蝌蚪视频在线观看| 日韩国产黄色网站| 青草精品视频| 亚洲欧美精品在线| 午夜国产小视频| 国产AV无码专区亚洲A∨毛片| 亚洲欧美日韩另类在线一| 国产主播喷水| 国产精品福利尤物youwu| 无码aⅴ精品一区二区三区| 日韩高清在线观看不卡一区二区 | 久久久久亚洲AV成人网站软件| 国产亚洲精品91| 青青草国产在线视频| 东京热高清无码精品| 青青青国产视频| 国产三区二区| 亚洲日韩图片专区第1页| 亚洲不卡影院| 91成人试看福利体验区| AV在线天堂进入| 国产黄网永久免费| 91久久大香线蕉| 国产成人AV综合久久| 性做久久久久久久免费看| 国产精品99久久久| 在线观看视频一区二区| 免费高清a毛片| a在线亚洲男人的天堂试看| 久久精品视频一| 精品视频一区二区观看| 无码一区18禁| 国产丰满成熟女性性满足视频| 911亚洲精品| 国产乱人视频免费观看| 国产69精品久久久久妇女| 久久五月视频| 欧美α片免费观看| 国产精品美乳| 99热这里只有精品免费| 中日无码在线观看| 亚洲精品不卡午夜精品| 久久亚洲日本不卡一区二区| 国产第八页| 亚洲精品无码人妻无码| 国产乱肥老妇精品视频| 亚洲精品大秀视频| 四虎综合网| 国产99欧美精品久久精品久久| 高清欧美性猛交XXXX黑人猛交| 久久香蕉欧美精品| 曰韩人妻一区二区三区| 国产成人做受免费视频| 91久草视频| 国产SUV精品一区二区| 伊人精品成人久久综合| 国产福利免费视频| 精品国产中文一级毛片在线看 |