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

基于均布式多孔表面吹氣的索結構尾流控制

2022-05-30 10:48:04高東來余海洋陳文禮
湖南大學學報·自然科學版 2022年5期

高東來 余海洋 陳文禮

摘 要:為研究基于均布式多孔表面吹氣對索結構尾流的控制效果,進行了一系列的風洞 試驗.通過粒子圖像測速(PIV)系統對雷諾數Re為1.0×104的無控和控制工況的索結構尾流 場進行了測量,分析了尾流的瞬時和時均流動特性,并利用本征正交分解(POD)和動態模態分解(DMD)對降階后的模態特性進行了分析和對比.在吹氣控制中,控制參數為無量綱的等效吹氣系數CQ.研究結果表明:隨著 CQ的增加,POD 各模態的能量分布趨于一致,流場中的擬序結構尺度趨于均一化,對旋渦脫落起控制作用的第1、2階模態受到抑制;尾流中旋渦脫落頻率被改變,分離的剪切層間的相互作用得到削弱;多個DMD模態特征值被改變,模態幅值的頻域分布發生偏移;尾流中回流區的尺度變大,湍動能和雷諾應力得到顯著削弱.

關鍵詞:索結構;流動控制;多孔材料;粒子圖像測速;模態分析

中圖分類號:U448.25 文獻標志碼:A

Control of Wakes of Cable Structures Based on?? Blowing through Uniformly Distributed Porous Surfaces

GAO Donglai1,2,YU Haiyang1,2,CHEN Wenli1,2?

(1.Key Laboratory of Smart Prevention and Mitigation of Civil Engineering Disasters of Ministry of Industry

and Information Technology(Harbin Institute of Technology),Harbin150090,China;

2.Key Laboratory of Structures Dynamic Behavior and Control of Ministry of Education

(Harbin Institute of Technology),Harbin150090,China)

Abstract:In order to investigate the control effects on a cable structure based on blowing through uniformly dis-tributed porous surfaces,we conducted a series of wind tunnel tests.The measurements of wake flow fields of the cable structure under the uncontrolled and controlled cases with Reynolds number Re1.0×104 were achieved by a PIV system.Instantaneous and time-averaged flow characteristics was analyzed,and proper orthogonal decomposi-tion(POD)and dynamic mode decomposition(DMD)were utilized to analyze step-down modal characteristics.The non-dimensional control parameter in the blowing control was the equivalent blowing coefficient,CQ.The experimen-tal results show that with the increase of CQ,POD modal energy distributions tend to be uniform,dimensions of co-herent structures in the global flow fields are transformed to be homogeneous and the dominance of the first two PODmodes is restrained.The wake vortex shedding frequency is changed and interactions between separated shear-layers are restricted.DMD modal eigenvalues are changed and spectral distributions of mode amplitudes shifted.Further-more,the recirculation region dimension in the wake is enlarged and the turbulence kinetic energy and Reynolds shear stress are significantly weakened.

Key words:cable structures;flow control;porous materials;particle image velocimetry(PIV);modal analysis

近年來,我國大跨度橋梁建設發展迅速,主跨千 米級的橋梁不斷建成.大跨度橋梁一般采用纜索承 重體系,具體的結構形式為斜拉橋與懸索橋.與中小跨度橋梁相比,大跨度斜拉橋和懸索橋的柔度大,動力效應更明顯.斜拉橋和懸索橋的橋址一般處于山 谷、江河和沿海地區,風環境較為復雜.大跨度柔性 橋梁為風敏感結構,在風作用下產生的效應十分明 顯,需對結構和構件的風效應進行評估[1].在定常氣流的作用下,結構或構件可能會由于氣動彈性失穩 而發生自激的、發散性的風致振動,如顫振和馳振;在脈動風的作用下,則可能發生限幅的強迫振動,如 抖振.對于大跨度橋梁,一種更為常見的風致振動形 式為渦激振動,該振動形式是由于鈍體繞流時的旋 渦脫落頻率和結構自身振動的某階頻率相接近時發生的鎖定現象而產生,在較低風速區間內發生,且振 幅有限,介于發散的自激振動與限幅的強迫振動之 間[2].為保證大跨度橋梁的正常使用和承載能力,需對風致振動進行抑制,既可通過抑制結構或構件自身的振動,也可通過對繞流場的調控加以實現[3].

索結構是大跨度橋梁重要的組成部分,如斜拉 橋的拉索、懸索橋的主纜和吊桿,以及下承式、中承 式拱橋的柔性吊桿等.以斜拉橋的拉索為例,其可能發生的振動形式為渦激振動、尾流馳振、風雨激振和參數共振等[2,4-5].Chen等[6]通過節段模型風洞試驗對西堠門大橋吊桿的橋塔尾流致振現象進行了分析,建立了氣動力模型和相應的吊桿運動方程.索的馳振和風雨激振的振幅過大,一般為幾倍的拉索直徑 甚至更大,可能會影響行車安全;索發生渦激振動時的振幅一般不大,但發生頻率最高,會使其耐久性受到影響,甚至會發生疲勞破壞.因此,控制并抑制大跨度橋梁索結構的風致振動具有十分重要的意義.

索結構振動的控制可分為結構措施、機械阻尼措施和氣動措施[2].結構措施一般是通過設置輔助索減小索的自由長度,增加索面剛度,以削弱振幅,其缺點是安裝困難,且對橋梁的外觀造型有影響.機 械阻尼措施一般是在拉索與橋面間設置阻尼器,耗 散振動的能量.常見的阻尼器有油阻尼器、剪切型黏滯阻尼器和磁流變阻尼器[7]等.氣動措施是指改變結構或構件的橫截面形式,改變外部氣流的流動特 性,以達到控制風致振動的目的.根據外部能量的消耗與否,流動控制方法可以分為主動控制與被動控 制[8].被動控制措施一般通過改變鈍體的外形,以改 善其流體動力學特性,從而對流動分離和尾流進行控制并抑制氣動力,尤其是降低阻力和脈動升力.被動控制方法的主要優點是節省能量并易于安裝[8].常見的被動控制方法有表面設置凸起[9-10]、安裝導流 板[11]、開槽[12]等.然而,與主動控制相比,被動控制方法往往難以達到非常明顯的控制效果,且無法對控制過程進行調節.

在主動流動控制方面,定常吸/吹氣是一種常見的控制方法.Chen等[13]通過風洞試驗研究了在雷諾數Re為3×104時定常吸氣的流動控制效果,需要指 出的是,該雷諾數位于斜拉索通常發生渦激振動時的風速所對應的雷諾數區間.研究結果表明,尾流區的旋渦脫落過程發生了變化,作用于圓柱的不穩定氣動力得到了削弱,從而抑制了圓柱的渦激振動.在控制過程中,吸氣孔在圓柱表面的分布方位角是十分重要的控制參數.Gao等[14]通過試驗研究了圓柱在迎風側吸氣與背風側吹氣同時作用下的控制效果.結果表明,作用于圓柱上的阻力和氣動力脈動幅 值均有降低,旋渦脫落頻率也發生了變化,后駐點的吹氣在尾流中產生了一對旋渦,對原始的旋渦脫落 過程有調節作用.

本文通過風洞試驗,研究了具有均布式多孔表面的圓柱在定常吹氣控制下的尾流特性,對尾流的流動特性如流線、湍動能、雷諾應力和渦量等進行了分析,并結合尾流區的頻譜特性,得到了不同吹氣控 制程度下的尾流變化過程.

1風洞試驗概況

試驗在哈爾濱工業大學大氣邊界層風洞與浪槽 聯合實驗室1號閉口回流式風洞(SMC-WT1)中進行.試驗段的寬度和高度均為505mm,長度為1000mm.試驗段的壁面為透明玻璃,可保證良好的流動可視化觀測條件.試驗段風速在0~24 m/s范圍內連 續可調,且上游位置所布置的蜂窩和柵格可保證試 驗段風速穩定,湍流度較小(約為0.30%).試驗的來 流風速U∞設為3.0m/s,對應的雷諾數Re為1.0×104.試驗所用的多孔索結構模型如圖1(a)所示,模型兩端為具有光滑表面的空心有機玻璃管,中間段為多孔段,材料為樹脂,采用3D 打印技術制作而成.模型的外 部直徑 D為50mm,內 部空心 段直徑 d為25mm,展向長度L為504 mm,其中多孔段的長度L0為64 mm.模型多孔段的內部構造如圖1(b)所示,沿模型環向均勻布置有30個氣孔,展向孔數為16個.模型表面單個氣孔直徑為3.2mm,氣孔總面積 Sh約為3858 mm2.在本試驗中同時研究了無控圓柱(光滑圓 柱)作為對比.圖1(c)所示為粒子圖像測速(particle image velocimetry,PIV)系統的示意圖,該系統主要 由高能雙脈沖激光器(Beamtech Vlite 430)、高分辨率雙曝 光 CMOS(complementary oxide semiconduc-tor)相機(pco dimax HS4)、煙霧發生器(ROSCO Al-pha900)和數字式延遲發生器(Berkeley Nucleonics Model 577)組成.試驗所用的Nd:YAG 激光器所釋 放的激光能量為0.435 J,波長為532 nm,雙路激光的時間間隔 Δt為0.1ms.CMOS相機的采樣頻率為200Hz,其與激光器之間的信號同步由數字式延遲發生器實現.示蹤粒子為平均直徑1~5 μm的油滴,通過煙霧發生器均勻散布于試驗段內.通過幀間互相關運算,可由相機采集到的圖像計算得到對應時刻的速度場,并 進一步 計算得到尾 流 面 內 渦量(ωz= ?v/?x-?u/?y)、湍動能(k=( + )/2U)和雷 諾 切應力(= v-′/U)等流動物理量.

模型的主動吹氣流動控制通過外部風機實現,風機與模型之間通過預留的氣孔與PVC軟管進行連 接.風 機的吹氣流率Q 通 過流量控 制器(Omega FMA-2613A)進行控制并保持恒定,流率范圍為24~ 216 L/min,步長為48 L/min.本文引入了無量綱參數——等效吹氣系數CQ對主動吹氣流動控制進行量化,其定義為:

式中:Uh為表面氣孔處平均吹氣流速,m/s.在風洞試 驗中,各流率所對應的主動控制工況的吹氣控制參數如表1所示.

為分析流場中主要的擬序結構,本文采用了本 征正交分解(proper orthogonal decomposition,POD)方法,將 PIV所觀測到的原始流動物理量由前15階POD模態進行重構,并分析了POD 前 4階模態的特 性.為進一步評估均布式多孔表面吹氣對索結構的尾流控制效果,采用動態模態分解(dynamic mode de-composition,DMD)方法對一些典型工況尾流場主要模態的動力學特性進行對比分析.

2結果與討論

2.1POD模態特性

模態能量是 POD模態排序的依據,低階的POD模態能量占比最高,第j階POD模態的能量占比Pj 定 義為該階模態的特征值λj與從流場中提取出的所有n階模態的特征值之和的比值,即:

在本文中,共提取了前 500階的POD模態進行計算,即式(2)中n= 500.前100階POD模態的累積 能量分布曲線如圖2所示,對于無控圓柱,低階次模態的相應能量占比較大,且隨階次的增高而降低,即低階次的POD模態代表整個流場中大尺度的擬序結構[8].無控圓柱前 4階POD模態的累積能量占比約為流場總能量的75%,與Feng等[15]的研究結果相似.從圖2中可看出,隨著等效吹氣系數CQ的增大,各工況的前100階POD模態累積能量占比呈下降趨勢,意味著索結構在均布式多孔表面吹氣控制下,流場中大尺度的擬序結構得到了有效的抑制.同時,低階次的模態能量相對降低,而高階次的模態能量相對升高,表明流場中的擬序結構尺度趨于均一化.均布 式多孔表面吹氣控制使得前 2階POD模態能量降幅 最為明顯,由于第1、2階POD模態對“2S”旋渦脫落模式中交替脫落的旋渦起控制作用[16],由此可知吹氣控制可有效控制索結構尾流中交替脫落的旋渦(即Karman 渦街),具體控制效果將在下文進行分析與討論.

流場經POD重構后的前 4階模態渦量分布如圖3所示.各工況的第1、2階POD模態渦量沿圓柱中心 線 Y/D=0呈對稱分布,第3、4階模態則呈反對稱分布.無控圓柱的前 2階模態渦量在近壁面區分布明 顯,沿順流向衰減.第1、2階模態的渦量分布的相位差約為π/2,該相位差和渦量分布值與尾流中Kar-man渦街的形成密切相關.第3、4階模態渦量分布幅 值較小,但沿順流向的分布趨勢明顯,表征尾流中脫 落的旋渦沿順流向的能量輸運過程[17].從圖3(b)和(c)中可看出,當施加均布式吹氣控制后,均布多孔表面的圓柱尾流中POD前 2階模態渦量分布得到削 弱,且向下游偏移,遠離圓柱壁面;而第3、4階模態 渦量分布得到增強.

試驗工況的前 4階POD模態系數分布如圖4所示.從圖4(a)可看出,無控圓柱的第1、2階POD模態系數幅值較大,且具有明顯的周期性,結合圖3(a)中模態渦量分布特征,進一步表明了前 2階模態對尾 流中的旋渦脫落起著主要的控制作用[17].第3、4階POD模態系數幅值較小,時間歷程也不具有明顯的周期性.當等效吹氣系數CQ 逐漸增大時,第1、2階模態系數幅值減小,且周期性不再明顯;第3、4階模態系數幅值增大,并達到與第1、2階模態系數幅值相似的水平,如圖4(b)和(c)所示.

POD模態系數的均值和均方差(mean-square-error,MSE)值如表2所示,其中MSE 值可描述模態系數的波動程度.從表2中可看出,無控圓柱的第3階POD模態均值的偏移較大,而第1、2階模態的波動程度最明顯.當CQ=0.0299時,第1~ 4階POD模態 均值進一步發生偏移,而 MSE 值均得到了削弱,尤其是第1、2階POD模態的削弱程度最大,即模態系數波動程度被明顯抑制.當CQ=0.096 7時,各模態的均值偏移幅度減小,而第3、4階模態的波動程度變大.以上結果與圖4中的現象一致,即隨著等效吹氣系數的增大,第1、2階POD模態得到了抑制,而第3、4階模態被增強.

2.2? 尾流瞬時特性

索結構尾流中一個旋渦脫落周期的瞬時面內渦量ωz分布如圖5所示,任意給定各工況的初始時刻 t0,相鄰時刻的時間間隔為T/4,其中T為各工況所對應的尾跡渦脫落周期,分別為0.08 s、0.11s和0.12 s,對應的斯托羅哈數(St)為0.208、0.152和0.139,其中無控圓柱的St 值與Fey等[16]在相同雷諾數下的研究結果相吻合.St為描述旋渦脫落周期或頻率特征的常用的無量綱參數,其定義為:

式中:fv為旋渦脫落頻率;T為旋渦脫落周期(fv=1/ T);D為鈍體的特征長度即圓柱的外部直徑;U∞為來 流風速.從圖5(a)中可以看出,無控圓柱兩側分離的剪切層之間相互作用明顯,尾跡渦脫落模式為典型的“2S”模式,即一個周期內有兩個反對稱的旋渦脫 落[18-19],這種旋渦脫落模式在鈍體尾流中所形成的一系列旋渦即所謂的Karman 渦街.隨著等效吹氣系數CQ的增大,剪切層間的相互作用減弱,渦量分布幅 值減小,隨時間的演變周期被改變,尾流中的旋渦脫 落現象已不再明顯甚至消失.

2.3 DMD模態特性

圖6給出了無控圓柱的DMD模態特征.圖6(a)所示為DMD模態特征值λj在復平面上的分布,實部為Re{λj },虛部為Im{λj },特征值λj 也被稱為Ritz 值.

從圖6(a)中可以看出,DMD模態的特征值為共軛復數對,主要分布在單位圓|λj |=1上,與Rowley等[20]和張揚等[21]的研究結果相似.DMD模態幅值分布隨無量綱頻率(St)的變化如圖6(b)所示,根據幅值大小對模態進行降序排序后可以看出,第1階DMD模態所對應的St 值即為無控圓柱實際流場中的旋渦脫落 頻率所對應的St數,與2.2 節中的結果相吻合.圖6中的結果表明,無控圓柱的DMD模態主要分布于低 頻段,且主要模態的幅值分布較為集中.

圖7所示為索結構在等效吹氣系數CQ=0.0299時的DMD模態特征.特征值λj 仍主要分布在復平面 內的單位圓|λj |=1上,如圖7(a)所示.與無控圓柱相比,該工況的模態幅值分布向低頻段進一步偏移,主 要 DMD模態的幅值得到增強,第2階模態所對應的St 值與2.2 節中的結果相同.

等效吹氣系數CQ=0.096 7時的DMD模態特征 如圖8所示.與上述工況相比,有較多模態的特征值λj分布于復平面內的單位圓|λj |=1內部,表明流場已不再以中性穩定為主,如圖8(a)所示.從圖8(b)所示的模態幅值頻域分布結果可以看出,DMD 各模態間的幅值分布趨于均勻,主要模態間的幅值差異已不再明顯,表現出明顯的寬頻分布特征.

2.4? 時間平均特性

索結構尾流的時均流線和湍動能分布如圖9所示.從圖9(a)中可以看出,無控圓柱的尾流回流區范圍約為X/D≤1.7,-0.5≤Y/D≤0.5.尾流中的湍動能較大,峰值集中分布于X/D ≈0.8處.隨著等效吹氣系數CQ的增大,索結構尾流中的回流區尺寸變大,湍動能分布得到明顯削弱,如圖9(b)和(c)所示.當CQ=0.0299時,回流區范圍約為X/D≤3.2,-0.8≤Y/D≤0.8;當CQ增至0.096 7時,回流區范圍擴大到X/D≤3.5,-1.0≤Y/D≤1.0.此現象也與圖5所示的索結構在吹氣控制時的尾流區剪切層距離增大、相互作用受到削弱相一致.

圖10所示為索結構尾流的時均雷諾切應力分布情況,在所有的工況中,雷諾切應力均關于Y/D=0呈反對稱分布.從圖10(a)中可以看出,無控圓柱尾 流中的雷諾應力分布幅值較大,且距后駐點較近,表明湍流脈動對時均流動的影響較大.當進行吹氣控 制時,雷諾應力分布幅值受到削弱,且向下游偏移,如圖10(b)和(c)所示.

3結論

本文進行了一系列風洞試驗,通過PIV系統測量了無控和均布式多孔表面吹氣的索結構尾流速度場,在此基礎上分析了各試驗工況的瞬時和時均流動物理量變化情況,并通過POD和DMD等降階模型,對主要的模態特性進行了分析和對比,得出如下主要結論:

1)與無控工況相比,控制工況的尾流場中各 POD模態階次的能量分布趨于一致,流場中擬序結構的尺度趨于均一化;模態渦量分布得到削弱,第1、2階模態對反對稱交替脫落的旋渦起到的控制作用得到抑制,而第3、4階模態表征的順流向能量輸運效應得到增強.

2)索結構尾流的旋渦脫落頻率被改變,隨著等效吹氣系數CQ的增大,剪切層間的距離增大,其相互 作用得到抑制.當CQ 值足夠大時,初始脫落的旋渦會消失.

3)當等效吹氣系數CQ增大時,繞流場的DMD模態特征發生了明顯的變化,具體表現為:更多模態的特征值λj分布在復平面內的單位圓|λj |=1中;具有高幅值的主要 DMD模態在頻域中向低頻偏移,隨著CQ的增大,進一步地表現出寬頻特征.

4)隨著 CQ的增加,索結構尾流中的回流區沿順 流向和橫流向的尺度變大;湍動能和雷諾應力得到顯著削弱;雷諾應力向下游偏移,湍流脈動對時均流動的影響得到抑制.

參考文獻

[1]陳政清,劉光棟.橋梁風工程研究的若干新進展[J].工程力學,2006,23(S2):93-111.

CHEN Z Q,LIU G D.Some recent progresses in bridge wind engi-neering research[J].Engineering Mechanics,2006,23(S2):93-111.(In Chinese)

[2]陳政清.橋梁風工程[M].北京:人民交通出版社,2005.

CHEN Z Q.Wind engineering of bridge[M].Beijing:China Com-munications Press,2005.(In Chinese)

[3]高東來,陳文禮,楊文瀚,等.基于渦動力學的大跨度橋梁風效應流動控制[J].中國科學:技術科學,2021,51(5):517-529.GAO D L,CHEN W L,YANG W H,et al.Vortex dynamics-based flow control of wind effects on long-span bridges[J].Scien-tia Sinica(Technologica),2021,51(5):517-529.(In Chinese)

[4]DE S? CAETANO E.Cable vibrations in cable-stayed bridges

[M].Zurich,Switzerland:International Association for Bridge and Structural Engineering(IABSE),2007.

[5]李壽英,向琳琳,鄧羊晨.懸索橋吊索尾流致振非定常理論分析[J].湖南大學學報(自然科學版),2020,47(11):1-8.

LI S Y,XIANG L L,DENG Y C.Theoretical analysis on wake-induced vibration of suspension bridges hangers based on un-steady theory [J].Journal of Hunan University(Natural Sci-ences),2020,47(11):1-8.(In Chinese)

[6]CHEN W L,HUANG Y W,MENG H.Wake-induced vibration

of a suspender cable in the rear of a bridge tower[J].Journal of Fluids and Structures,2020,99:103166.

[7]李惠,劉敏,歐進萍,等.斜拉索磁流變智能阻尼控制系統分析與設計[J].中國公路學報,2005,18(4):37-41.

LI H,LIU M,OU J P,et al.Design and analysis of magnetorheo-logical dampers with intelligent control systems for stay cables[J].China Journal of Highway and Transport,2005,18(4):37-41.(In Chinese)

[8]CHOI H,JEON W P,KIM J.Control of flow over a bluff body[J].

Annual Review of Fluid Mechanics,2008,40:113-139.

[9]BEARMAN P,BRANKOVI? M.Experimental studies of passive

control of vortex-induced vibration[J].European Journal of Me-chanics-B/Fluids,2004,23(1):9-15.

[10]ZHOU B,WANG X K,GUO W,et al.Control of flow past a

dimpled circular cylinder[J].Experimental Thermal and Fluid Science,2015,69:19-26.

[11]GAO D L,HUANG Y W,CHEN W L,et al.Control of circular

cylinder flow via bilateral splitter plates[J].Physics of Fluids,2019,31(5):057105.

[12]GAO D L,CHEN W L,LI H,et al.Flow around a circular cylin-der with slit[J].Experimental Thermal and Fluid Science,2017,82:287-301.

[13]CHEN W L,LI H,HU H.An experimental study on a suctionflow control method to reduce the unsteadiness of the wind loads acting on a circular cylinder[J].Experiments in Fluids,2014,55(4):1-20.

[14]GAO D L,CHEN G B,CHEN W L,et al.Active control of circu-lar cylinder flow with windward suction and leeward blowing[J].Experiments in Fluids,2019,60(2):1-17.

[15]FENG L H,WANG J J,PAN C.Proper orthogonal decomposition

analysis of vortex dynamics of a circular cylinder under synthetic jet control[J].Physics of Fluids,2011,23(1):014106.

[16]FEY U,K?NIG M,ECKELMANN H.A new Strouhal-Reynolds-number relationship for the circular cylinder in the range 47 < Re < 2×105 [J].Physics of Fluids,1998,10(7):1547-1549.

[17]江建華,鮑鋒.基于POD方法開縫圓柱繞流流場的研究[J].氣體物理,2017,2(2):28-36.

JIANG J H,BAO F.POD analysis of a slit circular cylinder near wake[J].Physics of Gases,2017,2(2):28-36.(In Chinese)

[18]WILLIAMSON C H K,ROSHKO A.Vortex formation in the wakeof an oscillating cylinder[J].Journal of Fluids and Structures,1988,2(4):355-381.

[19]MA L Q,FENG L H.Vortex formation and evolution for flow over

a circular cylinder excited by symmetric synthetic jets[J].Ex-perimental Thermal and Fluid Science,2019,104:89-104.

[20]ROWLEY C W,MEZI? I,BAGHERI S,et al.Spectral analysisof nonlinear flows[J].Journal of Fluid Mechanics,2009,641:115-127.

[21]張揚,張來平,鄧小剛,等.飛行器大攻角復雜流動的POD和DMD對比分析[J].氣體物理,2018,3(5):30-40.

ZHANG Y,ZHANG L P,DENG X G,et al.POD and DMD analy-sis of complex separation flows over a aircraft model at high angle of attack[J].Physics of Gases,2018,3(5):30-40.(In Chinese)

主站蜘蛛池模板: 欧美亚洲欧美区| 99热这里只有免费国产精品| 2022国产无码在线| swag国产精品| 老司机精品久久| 日本高清在线看免费观看| 国产99精品视频| 动漫精品啪啪一区二区三区| 色综合天天综合中文网| 国产极品美女在线| 99视频在线免费看| a欧美在线| 熟女视频91| 亚洲欧美h| 久久亚洲天堂| 中国一级特黄视频| 久久婷婷综合色一区二区| 爆操波多野结衣| 亚洲第一成年人网站| 日本91在线| 欧美不卡二区| 色哟哟国产精品| 亚洲美女一区二区三区| 五月婷婷综合色| 亚洲天堂精品在线| www.youjizz.com久久| 国产亚洲视频在线观看| 日韩av高清无码一区二区三区| 国产精品美女自慰喷水| 久久久成年黄色视频| 欧美午夜性视频| 凹凸精品免费精品视频| 欧美精品黑人粗大| 日本不卡视频在线| 成人无码一区二区三区视频在线观看| 免费黄色国产视频| 72种姿势欧美久久久大黄蕉| 青青操国产视频| 伊人久久婷婷| 欧洲av毛片| aa级毛片毛片免费观看久| 日韩欧美国产另类| 伊人久久大香线蕉影院| 国产精品私拍在线爆乳| 久草视频一区| 青草视频网站在线观看| 久久久久人妻一区精品| 黄色网页在线观看| 香蕉精品在线| 伊人久久综在合线亚洲2019| 粉嫩国产白浆在线观看| 久久精品丝袜| 91精品福利自产拍在线观看| 青青久久91| 97se亚洲综合在线韩国专区福利| 成年免费在线观看| 毛片在线播放网址| 欧美日本在线观看| 天堂av综合网| 欧美性猛交xxxx乱大交极品| 亚洲精品无码不卡在线播放| 狠狠色综合久久狠狠色综合| 亚洲精品无码AⅤ片青青在线观看| 久久久国产精品免费视频| 国产成熟女人性满足视频| 色老二精品视频在线观看| 人妻丰满熟妇αv无码| 国产女人18水真多毛片18精品| 亚洲天堂成人在线观看| 最新日韩AV网址在线观看| 亚洲美女一区| 欧美高清视频一区二区三区| 午夜爽爽视频| 久久伊人久久亚洲综合| 亚洲中文字幕国产av| 99久久国产自偷自偷免费一区| 黄色免费在线网址| 欧美三级视频网站| 亚洲天堂视频网站| 亚洲系列无码专区偷窥无码| 特级精品毛片免费观看| 最新亚洲av女人的天堂|