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

雙分裂導線尾流誘發振蕩數值模擬研究

2017-03-09 10:04:15何小寶楊曉輝
振動與沖擊 2017年4期
關鍵詞:模態振動

何小寶, 嚴 波,2, 伍 川, 張 博, 楊曉輝

(1.重慶大學 航空航天學院,重慶 400044;2.重慶大學 輸配電裝備及系統安全與新技術國家重點實驗室,重慶 400044; 3.河南省電力公司電力科學研究院,鄭州 450052;4.國家電網公司輸電線路舞動仿真技術重點實驗室,鄭州 450052)

雙分裂導線尾流誘發振蕩數值模擬研究

何小寶1, 嚴 波1,2, 伍 川1, 張 博3,4, 楊曉輝3,4

(1.重慶大學 航空航天學院,重慶 400044;2.重慶大學 輸配電裝備及系統安全與新技術國家重點實驗室,重慶 400044; 3.河南省電力公司電力科學研究院,鄭州 450052;4.國家電網公司輸電線路舞動仿真技術重點實驗室,鄭州 450052)

利用Fluent流體動力學軟件模擬雙分裂導線繞流場,計算兩根子導線的空氣動力系數,得到背風側子導線的升力系數和阻力系數隨其與迎風側子導線相對位置的變化規律。進而利用ABAQUS軟件模擬典型雙分裂線路段的尾流誘發振蕩過程,獲得兩根子導線的運動軌跡、振動幅值和頻率特征。結果表明,兩子導線的位移相位相反,運動軌跡近似為水平橢圓,與現場觀測一致。模擬分析了間隔棒等間距和非等間距排布方式下雙分裂導線尾流誘發振蕩的特征。

雙分裂導線;空氣動力系數;尾流誘發振蕩;數值模擬

當風場中分裂導線的背風側子導線位于迎風側子導線的尾流馳振區時,可能會發生氣動失穩,誘發尾流馳振。分裂導線尾流馳振包括水平舞動、垂直舞動、扭轉振動和次檔距振動等[1],其中次檔距振動為主要形式。根據對二分裂和四分裂導線觀察表明:次檔距振動頻率約為1~3 Hz,振幅約為0.1~0.5 m,子導線的運動軌跡一般呈水平扁長橢圓。導線尾流馳振嚴重時會造成子導線間相互碰撞,產生“鞭擊”現象,甚至引起金具的破壞,對線路的安全運行造成嚴重威脅。

尾流區中子導線的氣動特性是研究分裂導線尾流馳振的基礎。PRICE[2]利用風洞試驗測量了光滑和膠股雙分裂導線在兩種湍流度、不同雷諾數下背風側子導線的升力系數和阻力系數與其在尾流區中位置的關系,但測量的位置限制在較窄的范圍內。DIANA等[3]在風洞中通過按次檔距振動時導線可能的運動形式來移動導線圓柱模型,測量了光滑和膠股雙圓柱的氣動參數,考慮了不同風速下導線移動頻率和幅度對氣動參數的影響。陳元坤[4]利用計算流體動力學方法,計算了不同雷諾數下雙分裂導線繞流場,研究了子導線間距、湍流強度、雷諾數等對背風側子導線氣動系數的影響,但沒有給出子導線氣動參數在整個尾流區中的變化規律。

自20世紀70年代,分裂導線尾流馳振問題在國際上就引起關注,但公開報道的研究成果較少。WARDLAW等[5]利用風洞試驗測量了光滑和絞股子導線尾流區穩定邊界,研究了分裂導線間距和湍流對穩定邊界的影響,并建立理論簡化模型分析了次檔距振動特征。TSUI等[6]利用有限元方法模擬了雙分裂導線尾流馳振問題,指出使用三維模型的必要性。HARDY等[7]通過大量現場觀測,研究了分裂導線尾流誘發振蕩特征,并分析了風向、導線張力、子導線排列方式、間隔棒排布方式和數目等對尾流馳振的影響。YANG等[8]給出了計算分裂導線次檔距振動極限幅值和最大次檔距間距的方法。李黎等[9]采用穩定性理論利用二維振子模型研究了雙分裂導線次檔距振蕩的機理。最近,本文作者課題組[10]對四分裂導線的尾流弛振進行了數值模擬研究,但忽略了背風側子導線氣動系數隨其在尾流區中位置變化的影響,模擬得到的導線運動軌跡與現場觀測結果存在較明顯的差異,存在不足。

本文首先利用Fluent流體動力學軟件模擬雙分裂導線的繞流場,獲得雙分裂導線氣動系數隨兩子導線相對位置的變化規律。利用ABAQUS有限元軟件建立典型雙分裂線路的三維有限元模型,模擬研究雙分裂導線的尾流誘發振蕩過程,分析次檔距振動的特征,以及間隔棒等間距和非等間距兩種安裝方式對次檔距振動的影響。本文方法和獲得的結果,對雙分裂導線次檔距振蕩抑制方法和技術的研究具有重要的參考價值。

1 導線空氣動力特性數值模擬

1.1 繞流場數值模擬

目前,關于導線次檔距振動和舞動等問題的研究,導線氣動參數的數值模擬和風動試驗大都采用準靜態方法,已有的研究表明,采用準靜態氣動參數能夠反映導線的主要運動特征。因此,本文仍采用準靜態分析方法。

已有的研究表明,迎風側子導線產生的尾流區邊界與尾流中心線的夾角約為±20°。假設子導線直徑為D,迎風側子導線對背風側子導線氣動特性的影響區域在流動方向約為50D[11],而在垂直于流動方向的影響范圍與其水平間距有關,一般認為在±5D范圍內影響最大[12]。

本文流場計算模型中,兩根子導線水平間距范圍取4D~50D,垂直間距范圍取-5D~5D。分析的子導線直徑為27.6 mm,尾流影響區域及計算模型參見圖1。計算來流為穩定流,風速為12 m/s。

湍流模型選擇Spalart-Allmaras模型[13],采用有限體積法,增量步為0.001 s。其中左邊界為速度入口,右邊界為壓力出口,上下邊界為對稱邊界。模擬雙分裂導線不同間距的繞流場時,迎風側子導線位置固定,通過改變間距Y和Z來改變兩根子導線之間的距離。

(a) 尾流影響 (b) 計算模型圖1 雙分裂導線尾流影響及計算模型Fig.1 Wake-induced andsimulation model of twin bundle conductor

在計算得到繞流場后,將導線模型表面的壓強進行積分,即可得到作用在導線模型上的升力和阻力。作無量綱化后得到其阻力系數和升力系數

(1)

式中:FD、FL分別為導線模型所受的阻力和升力;ρ為空氣密度;U為風速;L為導線模型的有效長度;D為導線的直徑。

迎風側子導線的氣動特性基本不受背風側子導線影響,其類似于單導線。因此,下面重點討論背風側子導線的氣動特性。如圖2所示為雙分裂導線水平間距Y為10D,垂直間距Z分別為0D、1.5D和5D時的速度場??梢?,兩子導線在風作用下均產生了渦脫落現象。當Z=0D時,背風側子導線受迎風側子導線尾流干擾嚴重,但隨著垂直距離Z的增大,這種影響逐漸減小,直至消失。

圖2 雙分裂導線子導線不同間距時的速度場Fig.2 Velocity field of twin bundled conductor with different distance

1.2 子導線氣動系數

如圖3所示為雙分裂導線不同水平間距時,子導線氣動系數隨垂直間距的變化規律??梢?,不同水平間距下,背風側子導線的阻力系數和升力系數隨垂直間距的變化趨勢相似。隨著水平間距的增大,最小阻力系數增大,而升力系數的最大值隨水平間距的增大而減小。

阻力系數曲線關于尾流中心線對稱,均為正值,且越靠近尾流中心越小。當垂直間距Z接近4D時其值趨于1.2,即接近單圓柱的阻力系數。而升力系數曲線關于尾流中心線反對稱,在尾流中心線上為0,隨著垂直間距的增大,升力系數先增大后減少,在尾流邊界處接近于0。在尾流中心線以上升力系數為負值,在中心線以下為正,即升力始終指向尾流中心。這些規律與風洞試驗結果一致,表明數值模型和結果是合理正確的。

圖3 不同水平間距時背風側子導線氣動系數隨垂直間距的變化Fig.3 Aerodynamic coefficients of leeward conductor varying with vertical interval in different horizontal interval

圖4為背風側子導線的阻力系數和升力系數在尾流區中的變化規律。從圖中可見,水平間距在30D以內,阻力系數的最小值變化較大;在此范圍之外,其最小阻力系數隨水平位置變化較??;阻力系數在垂直間距接近±4D時趨于最大值1.2。另一方面,子導線升力系數水平間距的影響范圍約為20D,垂直間距影響范圍約在±3D以內??傊?,背風側子導線的氣動特性在迎風側子導線尾流區內的變化明顯,在研究次檔距振動時不可忽略。

圖4 背風側子導線氣動系數在尾流區中的變化規律Fig.4 Aerodynamic coefficients of leeward conductor varying in the wake

2 線路模型及其固有頻率和模態

2.1 典型線路段模型

以檔距200 m的雙分裂線路段為研究對象,導線參數如表1所示。線路的初始水平張力為29.988 kN。線路上安裝3個間隔棒,每個間隔棒的質量為4.8 kg,采用非等間距和等間距兩種布置方案。非等間距布置方案按現行輸電線路設計規程[14]確定,第1個間隔棒距離左端43.33 m,其余依次增加70 m和53.33 m。等間距布置時,各間隔棒之間相距50 m。線路簡化模型如圖5所示,導線兩端固定,風垂直作用于線路。

表1 導線參數

圖5 雙分裂線路簡化模型Fig.5 Simplified model of twin double conductor line

2.2 整檔模態和次檔距局部模態

尾流誘發振蕩發生時,除了發生次檔距振動外,還可能發生整檔振動。整檔振動特征與整檔的模態和固有頻率有關,而次檔距振動與次檔距局部模態有關。利用ABAQUS軟件分析得到兩種間隔棒布置方式下該線路段在面內和面外的低階模態和固有頻率,如表2所示??梢姡g隔棒的布置方式對線路整檔模態和固有頻率的影響非常小。

表2 線路段整檔低階模態和固有頻率

Tab.2 Natural frequencies and modes of conductor line

方向模態固有頻率/Hz間隔棒非等間距間隔棒等間距面內0.4760.7041.0620.4760.7031.064面外0.3530.7060.7050.3531.0571.059

此外,間隔棒等間距布置時,線路在頻率1.418 Hz開始出現密集的次檔距局部模態;而間隔棒非等間距布置時,在1.011 Hz便出現次檔距局部模態。圖6所示為間隔棒兩種排布方式的典型次檔距局部模態和頻率。

圖6 線路段典型次檔距局部模態和頻率Fig.6 Local sub-span modes and frequencies of twin bundle conductor line

3 尾流誘發振蕩數值模擬

3.1 數值模擬方法

輸電導線屬于柔性結構,其運動過程中的變形屬于典型的幾何非線性問題。本文采用ABAQUS有限元軟件模擬尾流誘發振蕩過程。在ABAQUS軟件中,可將空間桿單元的材料性質設置為不可壓縮得到索單元來模擬導線,而間隔棒采用空間梁單元模擬。在此,導線單元尺寸取0.5 m,經檢驗,可以滿足單元收斂性。

作用在迎風側子導線上的氣動載荷僅有阻力,其阻力系數為1.2。作用在背風側子導線上的氣動載荷與其在尾流區中的位置有關,即與導線的運動狀態有關。根據式(1),作用在子導線上的氣動力由式(2)確定

(2)

式中:氣動系數CL(y,z)和CD(y,z)與導線運動過程中當前時刻的位置有關,通過其當前位置,可計算出其與迎風側子導線之間的相對距離,再利用圖4中的曲面進行插值確定其值。氣動載荷的施加利用ABAQUS用戶自定義單元子程序UEL實現[15]。

3.2 間隔棒等間距布置時次檔距振動特征

圖7為間隔棒等間距布置時,次檔距2的兩根子導線中點的位移時程曲線及運動軌跡??梢钥闯觯瑢Ь€的水平振動幅值較大,而垂直振動幅值較小,導線的運動軌跡接近橢圓,且橢圓主軸與尾流中心線的夾角較小。

圖8所示為間隔棒等間距布置時各次檔距中點子導線一個周期的運動軌跡。從圖中可見,子導線1與2的位移相位差接近180°,與現場觀察到的現象相一致[16]。此外,次檔距1中子導線1的運動沿順時針方向,子導線2的沿逆時針方向;而次檔距3中兩子導線的運動均沿逆時針方向,但這兩種運動方式中兩子導線的位移相位均相反。

圖9為間隔棒等間距布置時線路在某時刻的運動狀態。可見,此時次檔距1中點兩子導線呈張開狀態,即所謂的“呼”;次檔距2呈閉合狀態,即所謂的“吸”。同理,次檔距3的狀態為“呼”,次檔距4的狀態也為“呼”。這一現象與現場觀測結果也是一致的。

圖7 間隔棒等間距布置時次檔距2導線中點的位移時程和運動軌跡Fig.7 Time histories of displacements and motion traces at middle of sun-span 2 with equally arranged spacers

圖8 間隔棒等間距布置時各次檔距中點子導線的運動軌跡Fig.8 Motion traces of sub-conductors at the middle of sub-span with equally arranged spacers

圖9 間隔棒等間距布置時某時刻導線的運動狀態(位移放大10倍數)Fig.9 Motions status of conductor lines with equally arranged spacers at typical times (displacement magnification times: ten)

表3為間隔棒等間距布置時各次檔距中點子導線的振動幅值的均方根值(Root Mean Square,RMS)??梢?,垂直振幅的最小值為0.034 m,最大值為0.082 m;水平振幅的最小值為0.168 m,最大值為0.392 m。雖然間隔棒是等間距排布,但四個次檔距的振動幅值并不相等,且兩端次檔距的振動幅值比中間的要小。

分裂導線發生尾流誘發振蕩時,子導線間可能發生“碰線”。圖10為各次檔距中點兩子導線間距的時程曲線。由圖可知,各次檔距中點兩子導線的間距在0.4 m處上下波動。表4給出了各次檔距中點兩子導線的最大和最小間距,最小距離為0.146 m,最大距離為0.606 m,發生在次檔距2的中點。可知子導線并未發生碰線。

表3 間隔棒等間距排布時各次檔距中點兩子導線的振動幅值

Tab.3 Oscillation amplitudes of conductors at the middle of sub-spans with equally arranged spacers

圖10 間隔棒等間距排布時各次檔距中點兩子導線的間距隨時間的變化Fig.10 Space between two sub-conductors at the middle of sub-spans varying with time with equally arranged spacers

Tab.4 The maximum and minimum distance of two sub-conductors at the middle of sub-spans varying with time with equally arranged spacers

圖11所示為次檔距1兩子導線中點位移的頻譜。由圖11(a)可見, 水平位移的頻譜在0.353 Hz、0.700

Hz和1.390 Hz處出現峰值。由表2中結果可知,第一和第二個峰值分別接近于整檔面外單半波和雙半波振動模態的固有頻率,即整檔運動包含了單半波和雙半波模態響應。參見圖6可知,第三個頻率峰值與次檔距局部固有頻率1.418 Hz相接近,即線路發生了次檔距振動。

此外,由圖11(b)可知,垂直位移的頻譜在接近于面內單半波和雙半波模態對應的固有頻率附近出現峰值,而對應于次檔距固有模態的頻率峰值很小,說明發生次檔距振動主要為面外振動。

圖11 間隔棒等間距時次檔距1中點的位移頻譜Fig.11 Displacement spectra of sub-conductors at the middle of sub-span1 with equally arranged spacers

3.3 間隔棒非等間距布置時次檔距振動特征

本節討論間隔棒按現行輸電線路設計規程確定的排布方式下,其尾流誘發振蕩的特征。圖12為各次檔距中點子導線一個周期的運動軌跡??梢?,子導線1和子導線2的位移相位差接近180°,與間隔棒等間距布置時相一致。表5為各次檔距中點子導線振動幅值,子導線水平振幅仍遠大于垂直振幅。從間隔棒排布方式來看,次檔距2的長度為70 m,是四個次檔距中最大的,其振動幅值也最大;而次檔距4的長度為33.33 m,長度最小,其振動幅值也最小。因此,次檔距的長度對子導線振動幅值有較大影響。比較表3和表5中的結果可以看出,間隔棒兩種排布方式的垂直振幅相差不大,但等間距排布時的水平振幅遠大于非等間距排布的振幅。因此,采用非等間距排布方式可以大大減小尾流誘發振蕩的振動幅值。

次檔距水平振幅/m子導線1子導線2垂直振幅/m子導線1子導線210.1130.1120.0220.06020.2350.2880.0330.07030.1680.1830.0240.06340.0880.0740.0240.045

表6為各次檔距中點兩子導線的最大和最小距離。其中,最小距離為0.213 m,最大為0.502 m,子導線未發生碰線。間隔棒非等間距排布時次檔距中點距離變化范圍小于等間距排布。

表6 間隔棒非等間距排布時各次檔距中點子導線間的距離

圖13 間隔棒非等間距排布時次檔距4中點子導線2的位移頻譜Fig.13 Displacement spectra of sub-conductor 2 at the middle of sub-span 4 with equally arranged spacers

圖13為次檔距4子導線2中點位移的頻譜。水平位移的頻譜在0.356 Hz、0.499 Hz、1.044 Hz處出現峰值,第一個和第二個峰值分別對應水平單半波和垂直單半波振動模態頻率,即,此時線路整檔振動的主模態為單半波。第三個峰值對應于次檔距模態頻率1.011 Hz。

4 結 論

本文利用Fluent流體動力學軟件模擬得到雙分裂導線的氣動系數,并利用ABAQUS軟件模擬雙分裂導線的尾流誘發振蕩,得到如下結論:

(1)背風側子導線的氣動特性受其與迎風側子導線間距的影響明顯,研究分裂導線尾流誘發振蕩時不可忽略。

(2)數值模擬結果表明,尾流誘發振蕩過程中,導線的運動軌跡近似為水平橢圓,兩子導線的位移相位差反相,接近180°,且水平振動幅值比垂直振動幅值大,與現場觀察結果一致。

(3)間隔棒非等間距排布時導線的次檔距振動幅值明顯小于等間距排布情況。

(4)本文結果對雙分裂導線次檔距振蕩的抑制方法和技術研究具有重要的參考價值,其方法可以進一步用于研究不同線路的間隔棒布置方案的優化。

[ 1 ] SNEGOVSKIY D. Wake-induced oscillations in cable structures: finite element approach [D]. Moscow: University of Liege , 2010.

[ 2 ] PRICE S J. Wake induced flutter of power transmission conductors [J]. Journal of Sound and Vibration, 1975, 38(1):125-147.

[ 3 ] DIANA G, BELLOLI M. Wind tunnel tests on two cylinders to measure subspan oscillation aerodynamic forces[J]. IEEE Transactions on Power Delivery,2014, 29(3):1273-1283.

[ 4 ] 陳元坤. 分裂導線的微風振動與次檔距振動研究[D]. 武漢:華中科技大學,2011.

[ 5 ] WARDLAW R L, COOPER K R, KO R G, et al. Wind tunnel and analytical investigation into the aeroelastic behavior of bundled conductors [J]. IEEE Transactions on Power Apparatus and Systems, 1975, 94(2):642-651.

[ 6 ] TSUI Y T. On wake-induced vibration of a conductor in the wake of another via a 3D finite element method[J]. Journal of Sound and Vibration, 1986, 107(1):39-58.

[ 7 ] HARDY C, VAN DYKE P. Field observations on wind-induced conductor motions [J]. Journal of Fluids and Structures, 1995, 9(1):43-60.

[ 8 ] YANG J W, YU Z M, LUO S Q, et al. Optimization design on subspan oscillation of multi-conductor[J]. Advanced Materials Research, 2014, 940:347-350.

[ 9 ] 李黎,王騰飛,陳元坤. 分裂導線次檔距振蕩研究[J]. 振動與沖擊,2014,33(17):62-67. LI Li, WANG Tengfei, CHEN Yuankun. Research of sub-span oscillation of split conductor [J]. Journal of Vibration and Shock, 2014,33(17):62-67.

[10] 嚴波,蔡萌琦,呂欣,等. 四分裂導線尾流馳振數值模擬研究[J]. 振動與沖擊,2015,34(1):182-189. YAN Bo, CAI Mengqi, Lü Xin, et al. Numerical simulation on wake galloping of quad bundle conductor [J]. Journal of Vibration and Shock, 2015,34(1):182-189.

[11] BELLOLI M, MELZI S, NEGRINI S, et al. Numerical analysis of the dynamic response of a 5-conductor expanded bundle subjected to turbulent wind [J]. IEEE Transactions on Power Delivery, 2010, 25(4):3105-3112.

[12] RAWLINS C B. Fundamental concepts in the analysis of wake-induced oscillation of bundled conductors [J]. IEEE Transactions on Power Apparatus and Systems, 1976,95(4):1377-1393.

[13] CAI M Q, YAN B, Lü X, et al. Numerical simulation of aerodynamic coefficients of iced quad bundle conductors [J]. IEEE Transactions on Power Delivery, 2015, 30(4):1669-1676.

[14] 張殿生.電力工程高壓送電線路設計手冊[M]. 北京:中國電力出版社,2003.

[15] HU J, YAN B, ZHOU S, et al. Numerical investigation on galloping of iced quad bundle conductors [J]. IEEE Transactions on Power Delivery,2012,27(2):784-792.

[16] SOHN H K, LEE H K, CHU J H, et al. Field observation and analysis of wind-induced vibrations in four-bundled conductor transmission lines [J]. KIEE International Transactions on PE, 2003, 3(2):109-114.

A numerical simulation on wake-induced oscillation of twin bundle conductor lines

HEXiaobao1,YANBo1,2,WUChuan1,ZHANGBo3,4,YANGXiaohui3,4

(1. College of Aerospace Engineering, Chongqing University, Chongqing 400044, China; 2. State Key Laboratory of Transmission & Distribution Equipment and Power System Safety and New Technology, Chongqing University, Chongqing 400044, China; 3. Henan Electric Power Research Institute, Zhengzhou 450052,China; 4. Key Laboratory of Galloping Simulation Technology of Transmission Line, State Grid Corporation , Zhengzhou 450052, China)

Lift coefficient and drag coefficient of leeward conductors varying with relative position of upwind conductors were simulated by fluid dynamics software Fluent. The wake-induced oscillation of typical twin bundle conductor lines were numerically simulated by ABAQUS software. The moving traces, vibration amplitude and frequency characteristics of sub-conductors were obtained. It is shown that the displacement phases of the two sub-conductors are opposite, and the moving traces are approximated as a horizontal ellipse. The simulation results are consistent with the field observation results. Finally, the characteristics of wake-induced oscillation of the twin bundle conductor line with equally and unequally installed spacers are compared.

twin bundle conductor; aerodynamic coefficient; wake-induced oscillation; numerical simulation

國家自然科學基金(51277186);國家電網公司科技項目(521702140013)

2015-11-09 修改稿收到日期:2016-01-20

何小寶 女,碩士生,1990年9月生

嚴波 男,博士,教授,1965年5月生 E-mail: boyan@cqu.edu.cn

TM753;O39

A

10.13465/j.cnki.jvs.2017.04.010

猜你喜歡
模態振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
This “Singing Highway”plays music
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
車輛CAE分析中自由模態和約束模態的應用與對比
國內多模態教學研究回顧與展望
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 国产福利影院在线观看| 国产欧美日韩另类| 欧美色伊人| 99久久婷婷国产综合精| 天天婬欲婬香婬色婬视频播放| 91精品国产麻豆国产自产在线| 天堂在线www网亚洲| 天天摸天天操免费播放小视频| 国产高清国内精品福利| 国产99在线观看| 18禁不卡免费网站| 最新国产成人剧情在线播放| 国产精品男人的天堂| 在线亚洲小视频| 波多野结衣一区二区三区88| 一区二区三区四区在线| 免费人成又黄又爽的视频网站| 国产成人亚洲毛片| 亚洲无码精彩视频在线观看| 国产精品伦视频观看免费| 91在线高清视频| 久久国产高潮流白浆免费观看| 一区二区三区在线不卡免费| 五月激激激综合网色播免费| 久久久久国产精品免费免费不卡| 国产又爽又黄无遮挡免费观看| 国产精品真实对白精彩久久| 91色在线观看| 色综合久久88| 毛片在线区| 欧美不卡视频一区发布| 国产精品lululu在线观看| 久久天天躁夜夜躁狠狠| 日本成人不卡视频| 亚洲无码精品在线播放| 国产91麻豆免费观看| 国产精品久久自在自线观看| 国产jizzjizz视频| 国产黄色免费看| 亚洲成人精品久久| 最新国产在线| 亚洲色大成网站www国产| 日韩区欧美区| 色成人亚洲| 欧美视频在线观看第一页| 国产午夜人做人免费视频| 国产日韩av在线播放| 91av成人日本不卡三区| 天天操精品| 日韩一区二区三免费高清| 国产一二视频| 凹凸国产熟女精品视频| 九九这里只有精品视频| 999国产精品永久免费视频精品久久| 尤物国产在线| 97国产在线视频| 亚洲人成色在线观看| 国产极品美女在线观看| 一级高清毛片免费a级高清毛片| 国产理论一区| 国产精品蜜臀| 999国内精品视频免费| 在线永久免费观看的毛片| 美女一区二区在线观看| 国产在线观看一区精品| 一级香蕉人体视频| 精品视频第一页| 无码专区在线观看| 无码人中文字幕| 日韩精品成人网页视频在线 | 一本一道波多野结衣av黑人在线| 国产在线91在线电影| 手机永久AV在线播放| 国产成人免费手机在线观看视频| 亚洲欧美成人| 国产精品美女免费视频大全| 久久综合九九亚洲一区 | 97在线碰| 四虎永久在线| 黄色污网站在线观看| 怡红院美国分院一区二区| 玖玖精品视频在线观看|