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

螺旋線對斜拉索渦激振動特性影響的試驗研究

2023-12-04 04:52:28劉慶寬邵林媛孫一飛常幸韓鵬王仰雪
湖南大學學報(自然科學版) 2023年11期
關鍵詞:振動模型

劉慶寬 ,邵林媛 ,孫一飛 ,常幸 ,韓鵬 ,王仰雪

[1.省部共建交通工程結構力學行為與系統安全國家重點實驗室(石家莊鐵道大學),河北 石家莊 050043;2.河北省風工程和風能利用工程技術創新中心,河北 石家莊 050043;3.石家莊鐵道大學 土木工程學院,河北 石家莊 050043]

斜拉索是斜拉橋的主要受力構件[1],因其具有長度大、剛度小和阻尼低的特點,容易發生風雨振、渦激振動和馳振等風致振動[2-4].頻繁振動輕則增加結構疲勞應力和減少使用壽命,重則直接導致斜拉索損壞,所以斜拉索風致振動問題一直備受關注[5].

風雨振[6]是斜拉索在風-雨聯合作用下發生的極具危害性的大幅發散振動,與其表面形成的上水線有密切聯系[7];Gao 等[8-9]在風洞試驗中再現了斜拉索風雨振現象,從尾流動力學角度揭示振動機理,并應用準定常理論模型描述風雨振幅值.

對于風雨振的抑振措施,在斜拉索表面纏繞螺旋線的方法應用最為廣泛[10],關于螺旋線參數對斜拉索風雨振抑振效果的研究很多.Gu 等[11]、杜曉慶等[12]和Chang 等[13]針對特定截面尺寸斜拉索,研究表面纏繞螺旋線后對斜拉索風雨振的抑制情況.李明水等[14]和李文勃等[15]通過風洞試驗驗證螺旋線對斜拉索風雨振有良好的減振或抑制效果.劉慶寬等[16]研究了表面纏繞不同直徑和間距螺旋線的斜拉索風雨振振幅,總結出能有效抑制風雨振的螺旋線參數組合.

在對斜拉索風雨振研究過程中,王修勇等[17]、儲彤[18]和祝志文等[19]發現部分橋梁斜拉索存在渦激振動現象,這是常遇風速下發生的帶有自激性質的限幅振動[20-21],雖然振幅小,但發生頻次高,會增大錨固區應力,極易引起結構發生疲勞破壞.因此,有效抑制斜拉索的渦激振動具有極其重要的工程意義.Christiansen 等[22]研究了帶有螺旋線的斜拉索在均勻流場中的振動響應和相關截面的風荷載特性.Hung等[23]通過表面凹坑、肋條和纏繞螺旋線斜拉索的振動響應對比試驗,發現螺旋線的抑振效果最好,且在所有試驗條件下都很穩定,螺旋線的存在能顯著降低斜拉索漩渦脫落的相關長度.Yagi 等[24]通過試驗證明在斜拉索表面纏繞螺旋線后,分離點沿跨度方向發生改變,會增強三維流動特性,進而抑制卡門渦街漩渦脫落,此外,螺旋結構可以優化斜拉索的氣動阻力.Liu 等[25]、劉志文等[26]采用直徑140 mm 的斜拉索進行風洞試驗,在低阻尼比條件下,表面纏繞直徑分別為2.0 mm 和3.5 mm 的螺旋線時,斜拉索仍存在明顯的渦激振動,當螺旋線直徑達到10 mm,纏繞螺距為12 倍直徑的2 根或3 根螺旋線都可以減小斜拉索渦激振動振幅.Ishihara 等[27]利用LES 方法研究了螺旋線對斜拉索渦激振動的抑制機理,分別提出了有無螺旋線的斜拉索自由振動計算模型,并通過試驗進行驗證,發現螺旋線與斜拉索直徑比為0.1時可抑制近80%的渦激振動振幅.

Gao 等[28]研究發現斜拉索風雨振的某些特性與渦激振動存在相似之處.纏繞螺旋線對控制斜拉索風雨振有顯著作用,但對斜拉索渦激振動的影響缺乏深入研究,抑制渦激振動的螺旋線參數尚可進一步細化,同時,纏繞螺旋線后斜拉索氣動力的研究目前只針對特定雷諾數,考慮到斜拉索氣動力的雷諾數效應,纏繞螺旋線后斜拉索的氣動力隨雷諾數的變化規律需要進一步澄清.因此,本文針對無螺旋線斜拉索和纏繞螺旋線斜拉索,進行節段模型測振、測壓和測力風洞試驗,詳細研究了纏繞不同根數、直徑和間距的螺旋線對斜拉索渦激振動特性、風壓分布和氣動力特性等影響規律,為相關研究和工程設計提供依據和參考.

1 風洞試驗

在實際工程中,斜拉索直徑一般為78~190 mm[16],考慮到風洞尺寸、風速大小以及阻塞度等因素,選取直徑D=150 mm的斜拉索模型為研究對象.

試驗在石家莊鐵道大學風工程研究中心STDU-1風洞進行,如圖1 所示.測力和測壓試驗在高速段完成,試驗區尺寸為5.0 m×2.2 m×2.0 m,空風洞最大風速≥80 m∕s,背景湍流度≤0.2%;測振試驗在低速段完成,試驗區尺寸為24.0 m×4.4 m×3.0 m,空風洞最大風速≥30 m∕s,背景湍流度≤0.4%[29].

圖1 STDU-1風洞Fig.1 STDU-1 wind tunnel

節段模型外部采用有機玻璃管,能夠模擬斜拉索出廠后架設到斜拉橋上的表面狀態,內部貫穿通長鋼管,有機玻璃和鋼管之間通過環向加勁肋連接,模型具有足夠的強度和剛度.為消除模型的端部效應,綜合風洞尺寸和洞壁邊界層的厚度等其他限制因素,設置節段模型受力部分長度L=1 700 mm,此外,在測振、測壓和測力試驗中,模型兩端安裝45°邊緣倒角的圓形端板.根據文獻[30]的研究,端板尺寸定為直徑750 mm(5D)和厚度10 mm.同時在測力試驗中,模型兩側還安裝補償模型,補償模型直徑與節段模型直徑相同.端板固定在補償模型上,與斜拉索之間留有1~2 mm 縫隙,保證試驗時其他部分的力不會傳到斜拉索節段模型,端板具有足夠的剛度,避免試驗中可能發生的振動.

在測振試驗中,模型兩端分別通過4 根豎向彈簧懸掛于鋼架上,模型、端板、連接件和彈簧組成振動系統,由自由振動試驗得到豎向自振頻率和系統阻尼比,振動系統自振特性的具體參數見表1.ρ為空氣密度,根據試驗時洞內空氣的溫度、濕度和氣壓計算得到.在測力和測壓試驗中,將模型通過內置鋼管固定到風洞外部的剛性框架上,模型水平安裝,與來流方向垂直.

表1 振動系統自振特性的具體參數Tab.1 Specific parameters of self-vibration characteristics of vibration system

試驗風速及湍流度由Cobra probe探頭測試,風速測量范圍為2.0~100.0 m∕s,風速測量精度為±0.5 m∕s,采樣頻率為2 000 Hz,采樣時長為30 s.模型的振動幅值由激光位移計測得,量程為(120±60)mm,測量精度為6 μm,采樣頻率為1 000 Hz,采樣時長為30 s.表面壓力由微型ESP壓力掃描閥測量,量程為±2 500 Pa,采樣頻率為330 Hz,采樣時長為60 s.氣動力由DELTA 系列六分量高頻測力天平測量,最大量程為165 N,測量精度為1∕32 N,采樣頻率為1 500 Hz,采樣時長為30 s.

斜拉索表面布置5 圈環向測壓孔如圖2(a)所示,分別在模型的L∕4、3L∕8、L∕2、5L∕8 和3L∕4 處,每圈均勻設置36 個.根據風洞試驗段截面和模型尺寸,計算出阻塞度為5.80%,采用改進的Maskell 法[31]對所有模型的結果進行阻塞度修正.螺旋線參數包括根數、直徑和纏繞間距.其中,螺旋線間距是指同一根螺旋線相鄰兩圈之間的水平距離s,如圖2(b)所示.試驗在均勻流場中進行,工況設置如表2 所示,為方便描述試驗結果,用螺旋線直徑-間距-根數的方式對工況作簡化表示,例如2-2D-3表示在斜拉索表面纏繞直徑d=2 mm,間距s=2D的3根螺旋線.

表2 工況設置Tab.2 Condition setting

圖2 試驗模型Fig.2 Test model

為方便比較振動結果,通過式(1)將試驗風速無量綱化,得到折減風速Ur.

式中:U為試驗來流風速;f為彈性懸掛斜拉索節段模型的豎向振動頻率;D為斜拉索直徑.

氣動力系數的計算如式(2)所示.

式中:FD和FL分別為作用在節段模型上的平均阻力和平均升力;L為斜拉索的軸向長度;ρ為空氣密度.

試驗所測得的原始風壓數據為模型表面各測壓點的風壓值,通過式(3)可得無量綱風壓系數.

式中:x=1~36是測壓孔編號;Px為測得的模型表面壓力信號的時間序列;Pt為參考點總壓平均值;Ps為參考點靜壓平均值.

2 渦激振動響應特性研究

2.1 螺旋線根數對斜拉索渦激振動影響分析

試驗過程中記錄模型運動時每一時刻離開平衡位置的距離A,將A的均方根(R)與斜拉索直徑D的比值作為斜拉索無量綱位移,圖3 為斜拉索無量綱位移隨折減風速的變化.當折減風速Ur=4.7時,無螺旋線斜拉索開始發生渦激振動;當Ur=5.6 時,渦激振動達到最大幅值,最大無量綱位移為0.12;當Ur=6.9時,渦激振動結束,渦激振動的折減風速區間為Ur=4.7~6.9.

圖3 斜拉索無量綱位移隨折減風速的變化Fig.3 Change of dimensionless displacement of a cable with reduced wind speed

Larsen[32]指出斜拉索的最大無量綱渦激振動振幅隨斯克頓數Sc的變化規律符合GVPO(Generalized Van der Pol Oscillator)模型,如圖4 所示,最大無量綱渦激振動振幅與Sc關系的表達式為:

圖4 最大無量綱VIV振幅隨Sc數的變化Fig.4 The variation of maximum dimensionless VIV amplitude with Sc

式中:η對應本文中最大無量綱振幅;p為振動的無量綱時刻,在[0,2π]取值;m為單位長度質量;ξa為氣動阻尼比;?為氣動非線性參數;?為氣動指數.本次試驗結果和GVPO 模型一致性較好,表明當前試驗結果是可靠的,可以作為后續對比分析的基礎.

由圖3 可知,當螺旋線直徑d=15 mm、間距s=8D時,表面纏繞單根螺旋線后,斜拉索在折減風速Ur=3.6~6.7 時出現明顯的渦激振動現象,最大無量綱位移為0.10,同無螺旋線斜拉索相比,風速鎖定區間變寬,且最大無量綱位移僅降低約16.7%.纏繞同樣直徑和間距的2 根螺旋線后,斜拉索渦激振動被大幅抑制,最大無量綱位移為0.06,相比無螺旋線斜拉索降低約47.9%.當螺旋線增加到3 根時,斜拉索的渦激振動幾乎被完全抑制,無量綱位移基本穩定在0附近.當螺旋線直徑d=16 mm,間距s=8D時,斜拉索表面纏繞2 根或3 根螺旋線時,渦激振動被完全抑制.

在斜拉索表面纏繞不同根數螺旋線后,斜拉索渦激振動振幅和渦激振動區間會發生不同程度的改變.總的來看,單根螺旋線對斜拉索渦激振動影響較小,抑制振動效果遠不如同樣參數下的2根或3根螺旋線.

2.2 螺旋線直徑對斜拉索渦激振動影響分析

鑒于實際工程中2 根螺旋線的應用較多,因此,選取表面纏繞2根螺旋線的斜拉索進行深入分析.

當螺旋線間距s=8D、根數n=2 時,斜拉索VIV 響應如圖5 所示.當纏繞直徑d≤7 mm 螺旋線,最大無量綱位移在0.10附近,較無螺旋線斜拉索略微降低,意味著傳統的控制風雨振的小直徑螺旋線不能有效抑制斜拉索的渦激振動.當螺旋線直徑d=10 mm時,最大無量綱位移依然沒有較大變化,但渦激振動區間相較于d=7 mm有變窄趨勢.隨著螺旋線直徑增加到d=15 mm,斜拉索的渦激振動特性發生明顯改變,對應折減風速Ur=5.5~6.9,最大渦激振動振幅降至0.06,渦激振動現象得到有效抑制.當螺旋線直徑d=16~20 mm 時,無量綱位移近似為0,斜拉索不發生豎向振動,表明纏繞該直徑螺旋線能完全抑制渦激振動現象.

圖5 纏繞不同直徑螺旋線的斜拉索VIV響應Fig.5 VIV response of stay cables with different diameters of helical lines

纏繞間距s=8D、根數n=2 的螺旋線后,斜拉索最大VIV 振幅如圖6 所示.無螺旋線斜拉索的振動明顯,最大振幅為25.8 mm.對于表面纏繞直徑d≤10 mm 螺旋線的斜拉索,最大振幅穩定在20 mm 附近,相較于無螺旋線斜拉索只有微小降低.在螺旋線直徑d=15 mm 時,最大振幅僅為13.7 mm,此時渦激振動被大幅抑制.當直徑d≥16 mm 時,斜拉索的最大渦激振動振幅繼續降低,基本穩定在0 附近,即渦激振動被完全抑制.

圖6 纏繞不同直徑螺旋線的斜拉索最大VIV振幅Fig.6 The maximum VIV amplitude of stay cables with different diameters of helical lines

2.3 螺旋線間距對斜拉索渦激振動影響分析

考慮到纏繞直徑大于、小于和等于15 mm 螺旋線對斜拉索渦激振動有完全不同的抑振效果,基于3種不同螺旋線直徑,對纏繞不同間距螺旋線的斜拉索渦激振動特性展開更深入研究,得到如圖7 和圖8所示的試驗結果.

圖7 纏繞不同間距螺旋線的斜拉索VIV響應Fig.7 VIV response of stay cables with different spacings of helical lines

圖8 纏繞不同間距螺旋線的斜拉索最大VIV振幅Fig.8 The maximum VIV amplitude of stay cables with different spacings of helical lines

纏繞直徑d=10 mm 的雙螺旋線后,斜拉索在不同纏繞間距下的無量綱位移如圖7(a)所示.在圖中所列的纏繞間距下,螺旋線的抑振效果均不明顯,且圖8中對應工況的最大振幅與無螺旋線斜拉索相近.

對于圖7(b)中直徑d=15 mm 的雙螺旋線,當纏繞間距為2D、5D和11D時,斜拉索最大無量綱位移在0.10 附近,對應的折減風速依次為4.9、5.1 和5.6,而當纏繞間距為8D時,斜拉索最大無量綱位移為0.06,對應的折減風速為5.9.整體來看,除8D間距外,渦激振動區間隨螺旋線間距的增大逐漸后移,結合圖8 可以看出,螺旋線抑振效果與纏繞間距并非線性相關,當纏繞直徑d=15 mm、根數n=2時,螺旋線只有在間距s=8D時有較好的抑振效果,其他間距下抑振作用均不明顯.

纏繞直徑d=16 mm 與15 mm 螺旋線時渦激振動響應規律一致,尤其在間距s=11D時的抑振效果要弱于8D,如圖7(c)所示.在螺旋線直徑d=18 mm 和20 mm 時,斜拉索的無量綱位移和最大振幅都在0附近,即在不同間距下均能完全抑制渦激振動.總的來說,纏繞螺旋線間距對渦激振動的影響并不單一,還與螺旋線直徑有關.當纏繞的螺旋線直徑較小時,間距的影響并不明顯,且對渦激振動的抑振效果影響較小;當直徑較大時,間距的大小直接影響渦激振動的抑制效果.

3 螺旋線對斜拉索風壓分布的影響

3.1 平均風壓系數

為得到螺旋線對模型渦激振動的控制機理,選取模型中間截面處的環向測點進行表面風壓分析,結果如圖9所示.

圖9 L∕2截面處的平均風壓系數Fig.9 Mean wind pressure coefficient at L∕2 section

在圖9(a)中,無螺旋線斜拉索的平均風壓系數隨雷諾數增大經歷了對稱—不對稱—對稱的變化過程,根據文獻[33-35],這是邊界層從層流狀態向湍流狀態轉捩過程中,斜拉索模型兩側相繼出現分離泡造成的.此外,180°測點處的基壓系數絕對值隨雷諾數增大逐漸減小.圖9(b)中d=2 mm 螺旋線本身直徑較小,纏繞在斜拉索表面得到的平均風壓系數跟無螺旋線斜拉索類似,具有良好的對稱性,而纏繞圖9(c)中d=15 mm 的大直徑螺旋線會遮擋住部分測壓孔,對附近的數據產生不同程度的干擾,使得纏繞該直徑螺旋線的斜拉索平均風壓系數出現不對稱的結果和某些不連續的尖點,但總體而言,纏繞d=15 mm、s=8D螺旋線的斜拉索最小風壓系數大幅度提高,抑制了升力的出現,與圖9(b)的趨勢一致.

與圖9(a)的無螺旋線斜拉索相比,纏繞圖9(b)和圖9(c)螺旋線后,斜拉索背風面負壓區變窄,基壓系數絕對值整體變小,且隨雷諾數變化的幅度變小,根據文獻[36]的研究結果,基壓系數可以定性地反映平均阻力系數的變化規律,推測纏繞螺旋線后平均阻力系數減小且雷諾數效應減弱.

3.2 風壓系數功率譜

為研究螺旋線對斜拉索表面風壓頻率特性的影響,對風壓時程進行頻域分析.雷諾數Re=1.0×105時,L∕2 截面處90°測點的風壓系數功率譜如圖10 所示,測點位于氣流分離區,其功率譜可以反映無螺旋線斜拉索的漩渦脫落特性.從圖10(a)可以看出,無螺旋線斜拉索頻率峰值對應的斯特勞哈爾數St在0.2附近,有明顯的能量分布.而纏繞直徑d=2 mm、s=8D的2根螺旋線后,斜拉索風壓系數功率譜如圖10(b)所示.與無螺旋線斜拉索相比,雖然漩渦脫落強度 稍有減弱,但功率譜仍存在明顯峰值,意味著直徑d=2 mm螺旋線不能有效抑制振動.相反,表面纏繞d=15 mm螺旋線時,斜拉索風壓功率譜不存在頻率峰值,且能量分布基本在0附近,如圖10(c)所示,說明該螺旋線組合可以對無螺旋線斜拉索規則圓截面的漩渦脫落造成破壞,有很好的抑制振動效果.

圖10 L∕2截面處90°測點的風壓系數功率譜Fig.10 Power spectrum of wind pressure coefficient at 90°measuring point at L∕2 section

4 螺旋線對斜拉索氣動力的影響

除了風致振動之外,斜拉索的氣動力也是需要考慮的問題.對于大跨徑斜拉橋,斜拉索上的風荷載會占到全橋風荷載的很大比例,采用較小氣動力的斜拉索具有重要的實際工程意義,而纏繞螺旋線后氣動力變化如何,需進一步分析.

根據振動試驗結果,選取部分可抑制渦激振動的螺旋線參數,得到斜拉索平均氣動阻力系數隨雷諾數的變化情況,如圖11 所示.值得說明的是,相關研究[37]表明,對于同樣的風速,當風向同斜拉索軸向垂直時,斜拉索上的氣動力最大,因此,這里只研究風向同斜拉索軸向垂直的情況.

圖11 平均氣動力系數Fig.11 Mean aerodynamic coefficient

對于阻力系數,在給定雷諾數范圍內,無螺旋線斜拉索的平均阻力系數在亞臨界區基本穩定在1.1附近,之后在臨界區迅速減小至0.5,在雷諾數Re=3.5×105~4.0×105時平均阻力系數小幅度增大.而纏繞螺旋線后,平均阻力系數隨雷諾數變化的幅值變小,波動區間為0.9~1.1.在亞臨界雷諾數區,圖中所有工況的平均阻力系數都小于無螺旋線斜拉索,最多可減小10.9%,同3.1 節的表面壓力結果是一致的.圖11(a)中纏繞3 根螺旋線的斜拉索平均阻力系數要大于雙螺旋線,增量約5.8%,意味著在抑振效果差別不大的螺旋線參數設置中,要優先選擇雙螺旋線.

對于升力系數,在雷諾數Re=1.0×105~2.4×105時,無螺旋線斜拉索的平均升力系數保持在0 附近,當Re=2.4×105~3.4×105時,平均升力系數迅速增大到0.7 又降低至-0.4,如圖11(b)所示.而纏繞螺旋線后,斜拉索平均升力系數的最大值均有明顯減小,波動區間為0.1~0.4,同無螺旋線斜拉索相比,纏繞螺旋線后平均氣動力系數的雷諾數效應整體減弱,但不同參數螺旋線對雷諾數效應的削弱效果并不一樣.

斜拉索的脈動氣動力系數如圖12 所示.由圖12可知,在亞臨界雷諾數區,無螺旋線斜拉索的脈動升力系數遠大于脈動阻力系數,而纏繞螺旋線后,橫方向和順風向的脈動氣動力系數均有明顯減小.值得注意的是,在圖12(b)中,纏繞d=2 mm螺旋線斜拉索的脈動升力系數要大于圖中所示的其他大直徑工況,這與2.2節中纏繞小直徑螺旋線的抑振效果不如大直徑的規律一致.

圖12 脈動氣動力系數Fig.12 Fluctuating aerodynamic coefficient

在雷諾數Re=1.0×105時,3 種典型工況的升力系數功率譜如圖13所示.由圖13可知,纏繞螺旋線后,斜拉索功率譜的頻率峰值要小于無螺旋線的值,規律與圖10 的風壓系數功率譜一致,再次證明螺旋線可以減小斜拉索旋渦脫落的強度,并且在纏繞間距s=8D的2 根螺旋線時,直徑d=15 mm 螺旋線的抑振效果要優于d=2 mm.

圖13 升力系數功率譜Fig.13 Lift coefficient power spectrum

5 結論

通過斜拉索節段模型測振、測壓和測力風洞試驗,研究了不同參數螺旋線對斜拉索渦激振動特性和氣動力特性的影響,得到的主要結論如下:

1)在螺旋線根數方面,單根螺旋線抑振效果較差,2 根螺旋線效果較好,3 根螺旋線基本能完全抑制振動.對于氣動力系數,在其他參數相同時,纏繞2 根螺旋線的斜拉索平均阻力系數要小于3 根螺旋線.所以考慮到工程需要,在抑振效果差別不大的情況下,要優先選擇2根螺旋線.

2)對于螺旋線直徑,試驗結果表明,通常用于控制風雨振的小直徑d≤10 mm 螺旋線的渦激振動抑制效果并不明顯,風壓系數功率譜仍有明顯的頻率峰值.將螺旋線直徑增加到d≥15 mm時,可有效降低斜拉索渦激振動響應幅值,同時可以減弱雷諾數效應,尤其亞臨界雷諾數區的平均阻力系數和脈動升力系數顯著減小,證明螺旋線直徑較大時可以實現渦激振動響應的控制.

3)螺旋線的纏繞間距對斜拉索渦激振動特性的影響并不單一,與螺旋線直徑有關.當纏繞的螺旋線直徑較小時,渦激振動現象明顯,這時單純改變纏繞間距并不能有效抑制渦激振動;當螺旋線直徑較大時,間距的大小直接影響渦激振動的抑制效果.

猜你喜歡
振動模型
一半模型
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
重要模型『一線三等角』
This “Singing Highway”plays music
重尾非線性自回歸模型自加權M-估計的漸近分布
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 亚洲视频一区在线| 国产精品网址你懂的| 色综合激情网| 福利一区在线| 特级aaaaaaaaa毛片免费视频| 中文字幕va| 国产福利一区视频| 精品视频一区在线观看| 72种姿势欧美久久久大黄蕉| 国产免费a级片| 老色鬼久久亚洲AV综合| 国产精品片在线观看手机版| 免费看黄片一区二区三区| 露脸国产精品自产在线播| 日韩国产精品无码一区二区三区| 久久久精品无码一区二区三区| 18禁高潮出水呻吟娇喘蜜芽| 久久久久青草大香线综合精品| 亚洲无线国产观看| 毛片基地美国正在播放亚洲| 国产激情第一页| 9cao视频精品| 国产一级裸网站| 国产一国产一有一级毛片视频| 欧美日韩国产在线播放| 特级精品毛片免费观看| 日韩欧美中文字幕在线精品| 久久久久亚洲精品无码网站| 亚洲乱码在线播放| 精品午夜国产福利观看| 久久黄色视频影| 亚洲色无码专线精品观看| 在线观看热码亚洲av每日更新| 国产真实乱人视频| 自拍欧美亚洲| 中文字幕在线永久在线视频2020| 日韩在线2020专区| 午夜视频www| 真人高潮娇喘嗯啊在线观看 | 理论片一区| 亚洲欧美色中文字幕| 99无码中文字幕视频| 欧美日韩精品一区二区在线线| 亚洲精品少妇熟女| 亚洲午夜片| 久久久久国产一级毛片高清板| 五月婷婷综合网| 99re热精品视频中文字幕不卡| 先锋资源久久| 久青草免费在线视频| 国产在线一区视频| 91香蕉国产亚洲一二三区| 欧美狠狠干| 亚洲一区精品视频在线 | 制服丝袜国产精品| 亚洲国产一成久久精品国产成人综合| 久久久久人妻一区精品| 国产99视频精品免费观看9e| 国内精品视频区在线2021| 在线色国产| 国产在线观看精品| 国产国产人免费视频成18| 久久久受www免费人成| 国产v欧美v日韩v综合精品| 精品久久777| 伊人成人在线| 欧美综合成人| 内射人妻无套中出无码| 久久国产高潮流白浆免费观看| 综合色亚洲| 欧洲日本亚洲中文字幕| 久久国产高潮流白浆免费观看| 国产91精选在线观看| 国产精品任我爽爆在线播放6080| 国产精品亚洲五月天高清| 999国产精品永久免费视频精品久久| 97视频在线精品国自产拍| 中文字幕va| 人妻丰满熟妇av五码区| 成人小视频在线观看免费| 久久香蕉欧美精品| 精品国产乱码久久久久久一区二区|