亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        三維覆冰斜拉索風(fēng)致振動馳振分析

        2016-05-20 02:38:31譚冬梅王凱麗瞿偉廉高遠志
        振動與沖擊 2016年7期

        譚冬梅, 王凱麗, 瞿偉廉, 韓 玲, 高遠志

        (1.武漢理工大學(xué) 土木工程與建筑學(xué)院,武漢 430072; 2.湖北省建筑科學(xué)研究設(shè)計院,武漢 430071)

        ?

        三維覆冰斜拉索風(fēng)致振動馳振分析

        譚冬梅1, 王凱麗1, 瞿偉廉1, 韓玲2, 高遠志1

        (1.武漢理工大學(xué) 土木工程與建筑學(xué)院,武漢430072; 2.湖北省建筑科學(xué)研究設(shè)計院,武漢430071)

        摘要:斜拉索偏心覆冰后,氣動外形不再穩(wěn)定,在風(fēng)的作用下,可能誘發(fā)馳振;大幅的馳振將威脅到拉索的結(jié)構(gòu)安全,因此,有必要對覆冰斜拉索的馳振穩(wěn)定性進行深入的研究。應(yīng)用FLUENT中的SST k-ω模型對三維覆冰斜拉索的繞流場進行數(shù)值模擬,得到全攻角下的阻力系數(shù)、升力系數(shù)及馳振力系數(shù),依此判定覆冰斜拉索是否發(fā)生馳振,并得到某大跨斜拉橋部分斜拉索的馳振臨界風(fēng)速。結(jié)果表明:經(jīng)過三維模擬計算的覆冰斜拉索的部分馳振力系數(shù)小于零,并且馳振臨界風(fēng)速較小,斜拉索易發(fā)生覆冰馳振;與風(fēng)洞試驗和二維模擬數(shù)據(jù)比較,三維覆冰斜拉索模擬獲得的氣動力系數(shù)比二維更接近試驗值。

        關(guān)鍵詞:覆冰;斜拉索;氣動力參數(shù);馳振;臨界風(fēng)速

        大跨斜拉橋由于過冷卻水附著在拉索表面,致使拉索表面被包裹上冰層。拉索偏心覆冰后,氣動外形不再穩(wěn)定,在風(fēng)的作用下,可能誘發(fā)馳振。馳振是由于平均升力系數(shù)突降產(chǎn)生的氣動負阻尼引起的一種橫風(fēng)向風(fēng)致振動,屬發(fā)散性振動,對結(jié)構(gòu)破壞很大。近年來,很多專家對覆冰導(dǎo)線的馳振進行了大量的研究,主要有Den Hartog[1]機制和Nigol[2]機制,前者認為氣動負阻尼是馳振的關(guān)鍵原因,后者認為是冰翼氣動力產(chǎn)生的扭轉(zhuǎn)自激失穩(wěn)。李萬平等[3-4]曾對覆冰導(dǎo)線空氣動力特性進行了測試;嚴波等[5-6]對四分裂導(dǎo)線的舞動特性和尾流馳振進行了數(shù)值模擬研究;Braun等[7]利用數(shù)值模擬分析了雙分裂、三分裂和六分裂導(dǎo)線的繞流場。

        覆冰斜拉索由于其動力特性與覆冰導(dǎo)線有差異,并且斜拉索與水平方向成較大的夾角(一般>30°),其風(fēng)攻角可為0°~360°,較易滿足馳振的起振條件[8],因此,覆冰導(dǎo)線馳振不同于覆冰斜拉索馳振。Demartino等[9-10]利用風(fēng)洞試驗研究了橋梁拉索覆冰的不同類型,分析了覆冰對拉索氣動參數(shù)的影響和覆冰拉索的馳振穩(wěn)定性;李壽英[11]對覆冰拉索進行了風(fēng)洞試驗和二維模擬研究。但是,當ReD>250時,二維模擬計算的氣動參數(shù)會出現(xiàn)錯誤值,三維模擬計算對正確預(yù)測流動特性是有必要的[12],為驗證三維覆冰拉索模擬的精確性,本文運用FLUENT軟件,對三維新月形覆冰斜拉索的繞流場進行數(shù)值模擬,獲得Den Hartog系數(shù)和斜拉索馳振臨界風(fēng)速,并與二維和風(fēng)洞試驗數(shù)據(jù)進行比較,證明計算流體力學(xué)(Computational Fluid Dynamics,CFD)可以成為替代風(fēng)洞試驗研究馳振的一種高效手段,也為進一步的覆冰拉索風(fēng)致振動研究提供數(shù)據(jù)。

        1橫向馳振原理及臨界風(fēng)速計算

        1.1氣動力參數(shù)

        氣動力參數(shù)采用如下定義:

        CL=2FL/(ρU2LB)

        (1)

        CD=2FD/(ρU2LB)

        (2)

        式中:FL為斜拉索模型的升力,來流速度方向逆時針轉(zhuǎn)動90°為升力的正方向;FD為斜拉索模型的阻力,沿來流方向為正;ρ為空氣密度,取為1.225 kg/m3;U為前方均勻來流風(fēng)速;B為斜拉索模型截面特征長度;L為斜拉索模型長度;CL、CD分別為升力系數(shù)和阻力系數(shù)。

        1.2風(fēng)致斜拉索橫向馳振原理

        偏心覆冰拉索在風(fēng)的激勵下,在其上產(chǎn)生升力,誘發(fā)馳振。一個振動系統(tǒng)是否穩(wěn)定,取決于其阻尼項的正負,當系統(tǒng)的阻尼項為正時,系統(tǒng)穩(wěn)定;反之,系統(tǒng)失穩(wěn)[13]。依據(jù)Den Hartog馳振理論,當氣動力為負阻尼時結(jié)構(gòu)發(fā)生失穩(wěn),即馳振發(fā)生的描述為:

        (3)

        1.3臨界風(fēng)速計算方法

        當氣動阻尼與結(jié)構(gòu)阻尼之和為零時為馳振發(fā)生的臨界狀態(tài),則可得馳振臨界風(fēng)速的理論公式為[14]:

        (4)

        式中:m為振子質(zhì)量;ω為結(jié)構(gòu)振動圓頻率;ξ為阻尼比。

        由于拉索的一階模態(tài)頻率最小,相應(yīng)馳振臨界風(fēng)速最小,因此,應(yīng)計算拉索一階模態(tài)相對應(yīng)的馳振臨界風(fēng)速。則第一階模態(tài)馳振臨界風(fēng)速可由下式得到:

        (5)

        式中:m1為單位長度斜拉索質(zhì)量;ω1為斜拉索一階模態(tài)圓頻率;ω1=2πf1;f1取斜拉索一階模態(tài)頻率;ξ1為斜拉索一階模態(tài)的阻尼比。

        2模型建立與數(shù)值模擬方法

        2.1模型的建立

        覆冰斜拉索選取新月形典型冰型,斜拉索直徑為120 mm,覆冰厚度為50 mm,拉索模型長度為900 mm,其模型外形尺寸見圖1。計算流域采用矩形區(qū)域,計算區(qū)域的大小為3.4 m×5.1 m,拉索的中心放置于原點處,距離上游流體入口為1.7 m,距離下游流體出口為3.4 m,

        圖1 三維覆冰拉索模型示意圖Fig.1 Sketch of 3D stay cables with iced accretion model

        距離左右流域壁面均為1.7 m,詳見圖2,圖中X軸正方向為速度來流方向,來流方向逆時針轉(zhuǎn)動90°為Y軸正方向,Z軸為拉索的長度方向。

        圖2 計算區(qū)域尺寸及坐標Fig.2 Computational domain size and coordinate

        三維拉索繞流的網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格,拉索周圍邊界層網(wǎng)格采用外O形網(wǎng)格,對邊界層網(wǎng)格進行加密,邊界層徑向網(wǎng)格增長系數(shù)為1.05,網(wǎng)格劃分見圖3。

        圖3 覆冰拉索網(wǎng)格劃分Fig.3 Mesh of iced cable

        2.2數(shù)值模擬方法

        邊界條件的定義:上游流體入口定義為速度入口(Velocity inlet),下游流體出口定義為壓力出口(Pressure outlet),上下壁面定義為對稱邊界(Symmetry),其他定義為固壁邊界(Wall)。進口速度為12 m/s,湍流強度(Turbulent Intensity)設(shè)置為3.8%,湍流黏性率(Turbulent Viscosity Ratio)設(shè)置為10。

        求解器(Solver)選用基于壓力法的求解器(Pressure Based),三維空間(3D),采用非穩(wěn)態(tài)計算方法,時間步長取0.001s。湍流模型采用SSTk-ω(Shear Stress Transportk-ω)模型,松弛因子采用默認,動量、湍動能、比耗散采用二階迎風(fēng)格式,SIMPLEC算法。

        模擬計算時取5°攻角遞增,經(jīng)過計算后,為增加計算精確度,在升力系數(shù)曲線的明顯尖峰處增加計算工況,取2°攻角遞增,風(fēng)攻角示意圖見圖4,為了使模擬更精確,針對每一工況進行建模、劃分網(wǎng)格。

        圖4 風(fēng)攻角示意圖Fig.4 Diagram of wind attack

        3三維新月形覆冰斜拉索數(shù)值模擬結(jié)果

        3.1典型風(fēng)攻角下流場分布

        為了分析拉索展向不同位置處的表面風(fēng)壓,取模型展向六個監(jiān)控截面,監(jiān)控截面分別為Z=100 mm、Z=200 mm、Z=300 mm、Z=400 mm、Z=500 mm、Z=600 mm。

        圖5給出了同一時刻90°風(fēng)攻角下監(jiān)控截面的壓力云圖,從圖5可知,不同截面脫落旋渦中心的顏色有深有淺,說明其壓力值不一致,并且脫落散發(fā)形式也有著細微的差別,隨著拉索高度變化,相應(yīng)截面的壓力云圖是存在相位差的,說明了空間拉索繞流場的三維流動特性,這在二維模擬中是無法反應(yīng)到的,因此,三維模擬比二維模擬更接近于實際流動特性。

        圖5 90°風(fēng)攻角監(jiān)控截面的壓力云圖(單位:Pa)Fig.5 Pressure contour of monitoring section of 90°wind attack(unit: Pa)

        圖6所示為0°、90°風(fēng)攻角時的速度分布,由圖6可知,有明顯的尾流渦交替現(xiàn)象。在拉索迎風(fēng)面和背風(fēng)面速度較小,迎風(fēng)的兩側(cè)速度較大,而壓力在迎風(fēng)面出現(xiàn)最大值,背風(fēng)面出現(xiàn)最小值。

        圖6 典型風(fēng)攻角的速度分布(單位:m/s)Fig.6 Velocity contour of typical wind attack(unit: m/s)

        3.2全攻角下三維覆冰拉索的氣動力系數(shù)

        圖7分別是三維FLUENT模擬獲得的0°風(fēng)攻角和90°風(fēng)攻角時阻力系數(shù)和升力系數(shù)時程曲線。在0°風(fēng)攻角時,隨著穩(wěn)定渦街和旋渦脫落的形成,阻力系數(shù)達到穩(wěn)定,約為0.48,而升力系數(shù)在0附近呈穩(wěn)定的周期變化;在90°風(fēng)攻角時,隨著周期性的旋渦脫落,阻力系數(shù)在1.6附近呈周期變化,升力系數(shù)在-0.1附近呈周期性變化。由升力系數(shù)時程曲線進行傅里葉變換,即可得到旋渦脫落頻率和Strouhal數(shù),為渦激共振研究提供參考。

        圖7 典型風(fēng)攻角的阻力系數(shù)和升力系數(shù)時程曲線(風(fēng)速:12 m/s)Fig.7 Drag and lift coefficient time-history curve of typical wind attack(wind speed:12 m/s)

        通過模擬得到12 m/s風(fēng)速下各風(fēng)攻角的阻力系數(shù)和升力系數(shù)時程曲線,進而獲得時程穩(wěn)定后的阻力系數(shù)和升力系數(shù)的平均值,其隨風(fēng)攻角的變化規(guī)律(見圖8)。平均阻力系數(shù)曲線呈兩端低中間略高的形狀,70°~100°風(fēng)攻角的平均阻力系數(shù)波動相對比較大。平均升力系數(shù)曲線呈現(xiàn)“M”狀,曲線兩側(cè)各有一尖峰,升力系數(shù)在0°、90°和180°風(fēng)攻角處約為0,在20°風(fēng)攻角時達到最大值1.0,在120°風(fēng)攻角附近達到負峰值。

        將圖8中三維模擬的阻力系數(shù)和升力系數(shù)與二維模擬、風(fēng)洞試驗值相比較,三維模擬獲得的平均阻力系數(shù)曲線處于二維模擬和風(fēng)洞試驗數(shù)據(jù)曲線的中間,略高于風(fēng)洞試驗數(shù)據(jù)曲線,對于平均升力系數(shù)而言,三維模擬的數(shù)據(jù)曲線與二維模擬、風(fēng)洞試驗相應(yīng)數(shù)據(jù)曲線變化規(guī)律一致,均呈“M”狀,但在風(fēng)攻角60°~140°范圍內(nèi),三維模擬的升力系數(shù)與風(fēng)洞試驗值接近,而二維模擬數(shù)據(jù)與風(fēng)洞試驗差別較大,三維比二維模擬更接近于試驗值,三維模擬曲線比風(fēng)洞試驗曲線平緩。

        圖8 三維模擬覆冰拉索氣動力系數(shù)(風(fēng)速:12 m/s)Fig.8 Aerodynamic coefficient of 3D simulation for iced cable(wind speed:12m/s)

        3.3全攻角下三維覆冰拉索的馳振力系數(shù)

        為了判斷覆冰拉索發(fā)生馳振的可能性,運用式(3)計算可得到平均升力系數(shù)隨風(fēng)攻角變化曲線的斜率值與相應(yīng)平均阻力系數(shù)的矢量和,即馳振力系數(shù),但對于一個特定風(fēng)攻角而言,?CL/?α可能取CL-α曲線的左側(cè)曲線斜率或者右側(cè)曲線斜率[15]。本文將兩種情況下覆冰拉索隨風(fēng)攻角變化的馳振力系數(shù)列于圖9中,當某風(fēng)攻角處左側(cè)曲線斜率和右側(cè)曲線斜率分別對應(yīng)的馳振力系數(shù)同時<0時,才認為拉索處于不穩(wěn)定狀態(tài)。從圖9可知,在風(fēng)攻角25°、30°、35°、42°、153°、175°處,兩種情況的馳振力系數(shù)均出現(xiàn)負值,說明拉索處于不穩(wěn)定狀態(tài),易發(fā)生覆冰馳振。

        圖9 三維模擬覆冰拉索馳振力系數(shù)(風(fēng)速:12m/s)Fig.9 Galloping coefficient of 3D simulation for iced cable(wind speed:12m/s)

        4覆冰斜拉索馳振臨界風(fēng)速

        由于冬季華中地區(qū)空氣濕度大,會出現(xiàn)雨雪邊降邊凍,導(dǎo)線和拉索極易受覆冰災(zāi)害的影響。本文則以某大跨斜拉橋為工程背景計算了拉索的馳振臨界風(fēng)速。表1給出了大跨斜拉橋部分斜拉索的結(jié)構(gòu)參數(shù)[16]。

        表1 某大跨斜拉橋斜拉索參數(shù)

        由阻尼比1%及表中各參數(shù),通過式(5)獲得25°、30°、35°、42°、153°以及175°風(fēng)攻角下各個覆冰斜拉索兩種情況下的馳振臨界風(fēng)速,其數(shù)值見表2中(a)~(f)的數(shù)據(jù)。從表2可知在同一風(fēng)攻角下,左側(cè)曲線斜率和右側(cè)曲線斜率分別對應(yīng)的斜拉索馳振臨界風(fēng)速差別較大,在30°、153°、175°風(fēng)攻角處,左側(cè)曲線斜率和右側(cè)曲線斜率分別對應(yīng)的馳振臨界風(fēng)速均較小,最小臨界風(fēng)速僅有7.9 m/s,易發(fā)生覆冰馳振。覆冰斜拉索馳振臨界風(fēng)速隨著質(zhì)量和一階頻率的增大而增大,一般斜拉索長度越長,越容易發(fā)生覆冰馳振。

        表2 覆冰斜拉索馳振臨界風(fēng)速

        5結(jié)論

        本文利用FLUENT軟件對三維新月形覆冰斜拉索進行數(shù)值模擬,得到了覆冰拉索的氣動力參數(shù)、馳振

        力系數(shù)以及馳振臨界風(fēng)速,通過數(shù)值模擬結(jié)果分析得到如下結(jié)論:

        (1) 文中監(jiān)控截面的流場分布說明了空間拉索繞流場的三維流動特性,三維模擬比二維模擬更接近于實際流動。與風(fēng)洞試驗和二維模擬數(shù)據(jù)相比較,三維模擬得到的氣動力系數(shù)隨風(fēng)攻角的變化規(guī)律比二維模擬更接近于風(fēng)洞試驗值。

        (2) 三維模擬的新月形覆冰拉索在25°、30°、35°、42°、153°以及175°風(fēng)攻角時左側(cè)曲線斜率和右側(cè)曲線斜率對應(yīng)的馳振力系數(shù)均出現(xiàn)負值,拉索具有發(fā)生覆冰馳振的可能性。

        (3) 以某大跨斜拉橋為工程背景,研究出其斜拉索在新月形覆冰條件下,30°、153°和175°風(fēng)攻角的左側(cè)曲線斜率和右側(cè)曲線斜率分別對應(yīng)的馳振臨界風(fēng)速均較小,最小臨界風(fēng)速僅有7.9 m/s,說明其斜拉索易發(fā)生覆冰馳振,一般拉索越長,馳振臨界風(fēng)速越小。

        參 考 文 獻

        [ 1 ] Den Hartog J P. Transmission line vibration due to sleet[J]. Transactions of the American Institute of Electrical Engineers,1932,51:1074-1077.

        [ 2 ] Nigol O,Buchan P G. Conductor galloping-Part Ⅱ:torsional mechanism[J]. IEEE Transactions on Power Apparatus and Systems,1981,100(2):708-720.

        [ 3 ] 李萬平. 覆冰導(dǎo)線群的動態(tài)氣動力特性[J]. 空氣動力學(xué)學(xué)報,2000,18(4):414-420.

        LI Wan-ping. Dynamic aerodynamic characteristicsof the galloping of bundled iced power transmission lines[J]. ACTA Aerodynamic Sinica,2000, 18(4):414-420.

        [ 4 ] 李萬平,黃河,何锃. 特大覆冰導(dǎo)線氣動特性測試[J]. 華中科技大學(xué)學(xué)報,2001,29(8):84-86.

        LI Wan-ping,HUANG He,HE Zeng. Aerodynamic characteri-stics of heavilyiced conductors[J]. Journal of Huazhong University of Science and Technology, 2001,29(8):84-86.

        [ 5 ] 嚴波,蔡萌琦,呂欣,等. 四分裂導(dǎo)線尾流馳振數(shù)值模擬研究[J]. 振動與沖擊,2015,34(1):182-189.

        YAN Bo,CAI Meng-qi,Lü Xin,et al.Numerical simulation on wake galloping of quad bundle conductor[J].Journal of Vibration and Shock,2015,34(1):182-189.

        [ 6 ] 蔡萌琦,嚴波,呂欣,等. 覆冰四分裂導(dǎo)線空氣動力系數(shù)數(shù)值模擬[J]. 振動與沖擊,2013,32(5):132-137.

        CAI Meng-qi, YAN Bo, Lü Xin,et al. Numerical investigation on aerodynamic coefficients of iced quad bundle conductor[J]. Journal of Vibration and Shock,2013,32(5):132-137.

        [ 7 ] Braun A L,Awruch A M. Aerodynamic and aeroelastic analysis of bundled cables by numerical simulation[J]. Journal of Sound and Vibration,2005,51-73.

        [ 8 ] 李壽英,黃韜,葉繼紅. 覆冰斜拉索馳振穩(wěn)定性的理論研究[J]. 振動與沖擊,2013,32(1):122-127.

        LI Shou-ying,HUANG Tao,YE Ji-hong. Theoretical analysis of galloping stability for stay cables with iced accretions[J].Journal of Vibration and Shock,2013,32(1):122-127.

        [ 9 ] Demartino C,Ricciardelli F. Aerodynamic stability of ice-accreted bridge cables[J]. Journal of Fluids and Structures,2015,52:81-100.

        [10] Demartino C,Koss H H,Georgakis C T, et al. Effects of ice ac-cretion on the aerodynamics of bridge cables[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2015,138:98-119.

        [11] 李壽英,黃韜,葉繼紅. 覆冰斜拉索氣動力的試驗與數(shù)值研究[J]. 湖南大學(xué)學(xué)報,2012,39(8):1-6.

        LI Shou-ying,HUANG Tao,YE Ji-hong. Experimental and numericalinvestigations of the aerodynamic forces on stay cables with iced accretion[J]. Journal of Hunan University,2012,39(8):1-6.

        [12] Kravchenko A G, Moin P. Numerical studies of flow over a circular cylinder at ReD=3 900[J]. Physics of Fluids, 20 00, 12(2):403-417.

        [13] 騰二甫,段忠東,張秀華. 新月形覆冰導(dǎo)線氣動力特性的數(shù)值模擬[J]. 低溫建筑技術(shù),2008(1):86-88.

        TENG Er-fu,DUAN Zhong-dong,ZHANG Xiu-hua. Numerical simulations of aerodynamic characteristics of iced conductor with crescent shape[J]. Low Temperature Architecture Technology,2008(1):86-88.

        [14] 陳政清. 橋梁風(fēng)工程[M]. 北京:人民交通出版社,2005.

        [15] 溫曉光. 斜拉索干索馳振的機理研究[D]. 長沙:湖南大學(xué),2013.

        [16] 張武劍,王波,楊金保,等. 九江長江公路大橋斜拉索振動特性研究[J]. 世界橋梁,2012,40(6):47-51.

        ZHANG Wu-jian, WANG Bo,YANG Jin-bao,et al. Study of vibration property of stay cables of Jiujiang Changjiang river highway bridge[J]. World Bridges, 2012,40(6):47-51.

        Galloping analysis of wind-induced vibration for 3D stay cables with iced accretion

        TANDong-mei1,WANGKai-li1,QUWei-lian1,HANLing2,GAOYuan-zhi1(1.School of Civil Engineering, Wuhan University of Technology, Wuhan 430072, China;2. Hubei Provincial Academy of Building Research and Design, Wuhan 430071, China)

        Abstract:The stay cables with iced accretion have the off-centered cross section, and their aerodynamic shape is no longer stable. These may lead to their galloping vibration under the action of wind. A galloping vibration with a large amplitude threatens the security of stay cables, so it is necessary to analyze deeply the galloping stability of stay cables with iced accretion. Here, the flow field of 3D iced cables was simulated by applying SST k-ω model of FLUENT. The drag, lift and galloping force coefficients under the whole attack angle were computed to estimate the possibility of their gallop. The critical galloping wind speeds of part of stay cables of a certain cable-stayed bridge were obtained. The results showed that part of galloping force coefficients computed with 3D simulation of the iced cables are less than zero, and the critical galloping wind speeds are smaller, the stay cables are easy to gallop; compared with wind tunnel tests and 2D simulation data, the aerodynamic force parameters computed with 3D simulation of iced cables are more close to the test results than those computed with 2D simulation.

        Key words:iced accretion; stay cables; aerodynamic force parameters; galloping; critical wind speed

        中圖分類號:U448.27

        文獻標志碼:A

        DOI:10.13465/j.cnki.jvs.2016.07.024

        收稿日期:2015-01-27修改稿收到日期:2015-04-02

        基金項目:國家自然科學(xué)基金資助項目(51408452);國家重點實驗室開放基金資助項目(2013B114);中央高校基本科研業(yè)務(wù)費專項資金資助(2013-IV-035)

        第一作者 譚冬梅 女,博士,副教授,1976年生

        亚洲精品一区二区网站| 91免费播放日韩一区二天天综合福利电影| 亚洲一区二区久久青草| 成人性生交大片免费看激情玛丽莎| 亚洲av高清一区二区三| 天码人妻一区二区三区| 亚洲av日韩aⅴ永久无码| 中文字幕丰满人妻有码专区| 亚洲tv精品一区二区三区| 无码人妻久久一区二区三区app | 我也色自拍俺也色自拍| 色爱情人网站| 成在人线av无码免观看麻豆| 免费国产h视频在线观看86| 日本不卡不二三区在线看| 国产乱人对白| 中文字幕无码精品亚洲资源网久久| 亚洲 国产 韩国 欧美 在线| 国产三级精品av在线| 妺妺窝人体色www聚色窝| 久久亚洲Av无码专区| 一本久久伊人热热精品中文| 亚洲av福利天堂一区二区三| 激情 人妻 制服 丝袜| 亚洲精品国产不卡在线观看| 亚洲av日韩专区在线观看| 白丝兔女郎m开腿sm调教室| 亚洲黄色一级毛片| 一本久道在线视频播放| 野花香社区在线视频观看播放| 毛片在线播放a| 亚洲综合网一区二区三区| 精品国内日本一区二区| 国产精品久久久久影院| 国产午夜无码精品免费看动漫| 国产色视频在线观看了| 午夜理论片yy6080私人影院 | 国产精品一区二区无线| 在线视频青青草猎艳自拍69| 全亚洲最大的私人影剧院在线看| 国产精品理论片|