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

橋梁斷面渦激振動氣動效應特性識別

2017-05-10 01:10:57葛耀君曹豐產
哈爾濱工業大學學報 2017年3期
關鍵詞:風速振動系統

許 坤, 葛耀君, 曹豐產

(1.土木工程防災國家重點實驗室(同濟大學), 上海 200092; 2.北京工業大學 道路與橋梁工程研究所, 北京 100124)

橋梁斷面渦激振動氣動效應特性識別

許 坤1,2, 葛耀君1, 曹豐產1

(1.土木工程防災國家重點實驗室(同濟大學), 上海 200092; 2.北京工業大學 道路與橋梁工程研究所, 北京 100124)

為獲取橋梁斷面渦激共振過程中的氣動效應特性,提出一種基于風洞試驗的氣動效應特性識別方法,首先提出渦振系統總體阻尼及剛度瞬時特性識別方法,在此基礎上提出從渦振系統總體阻尼及剛度中分離氣動阻尼及氣動剛度的方法.以一個開口斷面風洞試驗數據為例對該方法的準確性進行了驗證,并對該斷面渦振氣動效應特性進行研究.結果表明:識別得到的系統總體阻尼比為零的等高線與結構渦振響應曲線能很好吻合,從而驗證了該方法的準確性.在此基礎上對渦振過程中作用于結構的氣動阻尼及氣動剛度特性進行了研究,獲得了氣動阻尼及剛度隨結構振幅和約化風速的變化規律,可為后續理論模型研究提供參考.

橋梁斷面;風洞試驗;渦激振動;氣動效應;識別方法

渦激振動是一種廣泛存在于工程領域的非線性振動現象.已有健康監測資料表明,大跨度橋梁較易出現渦激振動現象[1-3],必須在設計階段對其渦振性能進行準確研究.現有渦振研究仍需依賴物理風洞試驗手段,通過彈性懸掛節段模型或全橋氣彈模型對其渦振性能進行判定[4-5].然而由于模型加工誤差或懸掛彈簧的影響,氣彈模型或節段模型不可避免地與實際結構在阻尼比等方面存在某些偏差.由于結構渦振性能受模型質量-阻尼系數影響明顯,通過風洞試驗獲得的結構渦振響應僅能反映試驗模擬質量-阻尼系數工況下的結構渦振性能,無法覆蓋其他未進行模擬的質量-阻尼系數工況.實際結構渦振可發生模態有多階,各階模態質量或模態阻尼比并不一致,若要完全掌握結構渦振性能,需對渦振可發生的各階模態進行研究,這將帶來較大的試驗工作量.

為避免上述工作量,研究人員期望通過引入某些數學模型來模擬結構渦振現象,并利用該數學模型計算不同質量-阻尼系數工況下的結構渦振響應.文獻[6]最早基于Van der Pol振子的概念,提出了經驗非線性模型來描述橋梁斷面渦激振動現象.后續又有許多學者在此基礎上提出了許多改進形式,如:經驗線性模型及廣義非線性模型等[7-10].可惜現有經驗模型都是從模擬結構自限幅振動特性出發,模型參數需通過擬合試驗獲得的結構自限幅振動響應獲得.由于結構自限幅振動響應受結構質量-阻尼系數影響明顯,因此現有渦振模型更多的是描述一種已經發生的渦振現象,而難以用來預測其他未進行試驗的工況.

導致結構出現渦激振動的主要原因是結構受到的非線性氣動效應.如果能通過試驗手段獲得結構渦振過程中所受氣動效應特性,并利用某些數學模型模擬該效應,將該數學模型與結構運動方程結合后,理論上可用來計算不同參數(質量-阻尼系數)下的結構渦振響應.從這個角度出發,或許可突破現有方法的局限,以較小的工作量完成對結構各階渦振性能的完全掌握.實現上述目的的關鍵在于結構渦振氣動效應特性識別,因此本文有針對性的提出了相應識別方法,并以實際橋梁斷面風洞試驗為背景對上述方法進行了詳細驗證.本文首先介紹渦振氣動效應識別方法;然后以一個開口斷面風洞試驗為背景對本文方法進行驗證,并對該斷面所受氣動效應特性進行研究;最后給出了氣動效應的相應研究結論.

1 氣動效應識別方法

對于橋梁斷面而言,風致氣動效應往往引起結構振動頻率和系統阻尼特性的改變,因此本文較為關注氣動效應的兩個方面:氣動剛度效應和氣動阻尼效應.

1.1 氣動剛度識別

渦激振動具有自限幅振動特性,運動過程中系統振幅、相位等會隨時間變化,因此需要對系統瞬時特性進行識別.由于氣動剛度的作用,運動過程中的系統振動頻率實時變化,因此,需對系統瞬時振動頻率進行識別.

其中P為柯西主值.

基于希爾伯特變換的系統辨識理論以及上述瞬時頻率計算公式的推導過程可參見文獻[11-13],關于這部分研究的文獻較為詳實,本文不再贅述.

上述方法雖具有理論意義,但是用于處理實際試驗數據時對數據本身要求極其嚴格,識別結果容易受到某些局部極值影響,文獻[14-15]曾對此進行過詳細的分析.為了克服上述困難,文獻[15]提出了用于瞬時特性計算的DQ方法(directlycomputequadraturemethod),并在其論文中對DQ方法、傳統HHT方法、以及其他常用方法進行了詳細比較,并驗證了DQ方法相較傳統HHT方法的優勢.因此,在本文渦振系統瞬時特性研究中,系統特性識別采用文獻[15]提出的DQ方法.

由于瞬時頻率的概念建立在窄帶振動信號基礎上,一般需首先對系統振動位移進行EMD(empiricalmodedecomposition)分解,將其分解為若干單組分固有模態函數(intrinsicmodefunction,IMF).將系統振動位移進行EMD分解后,得到的單組分IMF可表示為

其中:Ai(t)為第i個單組分成分Fi(t)的包絡項(envelope);φi(t)為瞬時相位;cosφi(t)為Fi(t)的振動承載項(carrier).

瞬時相位φi(t)關于時間t的導數,即瞬時頻率為

采用Huang的方法計算瞬時頻率時,需首先將Ai(t)項和cosφi(t)進行分離,具體分離步驟如下.

步驟1 將Fi(t)取絕對值得到|Fi(t)|;

步驟3 用3次樣條曲線插值上述局部最大值點離散序列Lmax(k),得到Fi(t)的第1次經驗包絡線e1(t);

步驟4 利用e1(t)對Fi(t)進行第一次正則化處理:

步驟5 理想情況下對Fi(t)進行正則化處理后N1(t)的值應該全部小于等于1,然而由于在某些局部位置樣條曲線可能無法完全穿過所有離散點Lmax(k),因此需要N1(t)基礎上重復步驟4,即

經過上述步驟后,即可分別得到cosφi(t)、Ai(t)分別為

實現上述計算步驟的前提條件是窄帶信號Fi(t)關于零軸對稱(零均值),并且振動信號的第1個點為某一周期的最大點(即滿足cos函數樣式),上述前提條件在實際試驗數據處理過程中可輕松實現.

獲得窄帶信號Fi(t)包絡項Ai(t)和振動承載項cosφi(t)后,瞬時相位的計算公式為

(8)

根據式(4)對瞬時相位關于時間t求導,即得到瞬時頻率ωi(t). 根據包絡項Ai(t)和瞬時頻率ωi(t)的一一對應關系,可得到隨瞬時振幅變化的瞬時頻率ωi(A).上述剛度項是渦振系統總體剛度效應,氣動剛度需在此基礎上減去模型在無風情況下得到的懸掛系統本身剛度,具體將在第2部分討論.1.2 氣動阻尼識別

渦激共振過程中,作用于結構表面的氣動力與結構運動狀態有關.結構周邊流場對結構的影響主要體現在附加氣動阻尼上.當把結構和周邊流場看作一個統一的動力系統時,系統總體阻尼特性中即包含懸掛體系機械阻尼影響又包含了結構所受附加氣動阻尼的影響.

為了獲得系統總體阻尼特性,需要用到系統“衰減-共振”或“增長-共振”瞬變段位移時程.以“衰減-共振”為例,對結構施加一個大于穩定渦振振幅的初始激勵后,結構運動會自起始狀態出發逐漸衰減至穩定振動狀態,此時的瞬變段類似于自由振動衰減過程,結構振幅可以表示為

其中Ai0為某一初始振幅,ωni為系統特征頻率,ωni與系統實際振動頻率ωi之間滿足

對式(9)兩邊取對數可得

結合式(10)可得

對式(12)關于ξi求解,可得到瞬時阻尼比為

同樣,根據瞬時阻尼比和瞬時振幅的對應關系,可得到隨振幅變化的系統總體阻尼比ξi(A).上述系統總體阻尼中包含了懸掛體系機械阻尼和結構所受附加氣動阻尼兩部分的影響,因此附加氣動阻尼效應識別需在總體阻尼比中扣除無風情況下得到的懸掛體系機械阻尼的影響.

2 開口斷面氣動效應特性

為驗證上述氣動效應識別方法用于實際橋梁斷面渦振研究的可行性,采用某一實際橋梁斷面風洞試驗數據為背景來進行研究.試驗對象為某一開口斷面斜拉橋,模型高44.9mm,寬559.3mm,長1 740mm,單位長度模型質量為10.46kg/m,豎向振動頻率為3.37Hz,阻尼比為0.4%.試驗在同濟大學TJ-2號風洞展開,模型斷面尺寸及形狀如圖1所示,試驗獲得的節段模型渦振特性如圖2所示.

圖1 開口斷面形狀(mm)

Fig.1 Configuration of the open-section bridge deck (mm)

圖2 渦振響應約化振幅

若要獲得系統氣動效應隨振幅的變化規律,需利用試驗測得的“衰減-共振”或“增長-共振”位移時程進行系統特性識別.考慮到“衰減-共振”位移時程能覆蓋較大的振幅范圍,本文開口斷面氣動效應識別采用試驗獲得的“衰減-共振”位移時程.鎖定區某一風速條件下得到的“衰減-共振”位移時程如圖3所示,圖中振幅和時間為無量綱形式.

圖3 模型“衰減-共振”位移時程

Fig.3 “Decay-to-resonance” displacement of the cross-sectional model during VIV

首先對試驗獲得的模型“衰減-共振”位移時程進行EMD分解以獲得窄帶的固有模態函數(IMFs);然后利用上一節渦振系統瞬時特性識別方法獲取系統總體阻尼及總體剛度隨約化振幅的變化規律. 對模型“衰減-共振”位移進行EMD分解后可獲得不同頻率的IMFs.對于一個試驗信號而言,引入EMD分解的目的是將其分解為不同的窄帶信號,然后對各窄帶信號的頻率特性進行識別.由于節段模型渦振位移時程本身即為窄帶信號,對其進行EMD分解后僅能獲得一個主導IMF,因此,系統特性識別時采用該主導IMF進行.此主導IMF相當于在原始位移信號的基礎上進行了時域濾波.

2.1 渦振系統剛度特性

基于上述主導IMF按照1.1節瞬時頻率識別方法進行計算,可得到不同風速條件下渦振系統瞬時振動頻率隨結構約化振幅的變化規律,如圖4所示.其中U/(fnD)=0表示彈性懸掛體系在無風情況時的振動頻率. 圖中實心矩形為根據試驗數據識別得到的系統瞬時頻率成分,實線為四次多項式擬合結果.

圖4 不同風速條件下系統振動頻率隨約化振幅的變化規律

從圖4可以看到:不同約化風速條件下系統振動頻率之間的差別相對系統振動頻率而言十分微小.同一約化風速條件下,系統振動頻率并不隨系統振幅的改變而變化,可以認為:渦振系統振動頻率不隨系統振幅變化,即某一風速條件下流場對結構的附加剛度效應可看作常數.零風速(U/(fnD)=0)條件下的結果表明,懸掛系統本身振動頻率也不隨系統振幅的改變而變化,這說明試驗過程中節段模型的懸掛彈簧始終處于彈性工作范圍.

2.2 渦振系統阻尼特性

獲得鎖定區不同風速條件下渦振系統瞬時振動頻率后,可按照第1.2節計算方法對不同風速條件下渦振系統瞬時阻尼特性進行識別.

圖5給出了不同風速條件下渦振系統總體阻尼比隨約化振幅的變化規律,圖中實心矩形為根據試驗數據識別的阻尼比,實線為四次多項式擬合結果.

圖5表明:不同約化風速條件下渦振系統總體阻尼比隨約化振幅的變化規律并不一致.在鎖定區外(U/(fnD)=7.27和U/(fnD)=12.56),系統總體阻尼比隨振幅的變化規律與零風速條件下總體阻尼比變化規律較為接近,呈現“下凹”形,即當結構振幅變小時系統總體阻尼比隨振幅的變化規律也變得較為“緩和”,在較小的約化振幅條件下系統總體阻尼比更接近于常數.

圖5 不同風速條件下系統總體阻尼比隨約化振幅的變化規律

鎖定區內(U/(fnD)=8.6~11.89),系統總體阻尼比隨約化振幅的變化規律呈現“上凸”形,即隨著約化振幅的減小,系統總體阻尼比迅速變小,鎖定區內系統總體阻尼比隨約化振幅的變化規律與鎖定區外或零風速條件下的變化規律明顯不同.

固定風速條件下系統總體阻尼比隨約化振幅的變化規律呈現明顯非線性特性,這表明渦振系統總體阻尼比同時與系統約化風速及約化振幅相關,系統總體阻尼比為約化振幅和約化風速的二元函數.

值得注意的是,零風速條件下彈性懸掛系統本身阻尼比也與系統約化振幅相關.這是由于彈性懸掛體系的阻尼主要由懸掛彈簧相鄰線圈間摩擦耗能提供,而不同振幅條件下彈簧相鄰線圈間運動周期內的摩擦耗能并不相等,因此即使彈簧處于彈性工作范圍,不同振幅條件下的系統機械阻尼比也并不一致.

由于零風速條件下彈性懸掛體系本身阻尼比與約化振幅有關,并不是一個常數,因此在識別渦振系統氣動阻尼效應時必須扣除這部分懸掛系統本身阻尼效應帶來的影響.傳統基于系統自限幅振動位移識別經驗模型參數的過程中并未考慮這部分影響,而是假設其阻尼比為常值,這也是導致傳統經驗模型無法準確模擬附加氣動阻尼隨結構振幅變化規律的因素之一.

根據不同風速條件下識別獲得的系統總體阻尼比隨約化振幅的變化規律,可得到總體阻尼比在不同約化風速和約化振幅條件下的等高線圖.系統總體阻尼比等高線圖可用來直觀地反映渦振鎖定區內外系統總體阻尼比的變化趨勢.

圖6給出了開口斷面渦振系統總體阻尼比等高線圖.圖中等高線的值采用圖5中的四次多項式擬合值.由于試驗過程中鎖定區風速間距較大,相鄰風速間距間的總體阻尼比等高線值通過插值獲得.

從圖中可以看到:在某些約化振幅和約化風速范圍內系統總體阻尼比為負值(圖中深色區域).由于系統總體阻尼比為負,處于該區域的渦振系統是不穩定的,系統振幅會逐漸增加,此時系統總體阻尼比會沿著縱軸向淺色區域發展,當發展至總體阻尼比為零的等高線處附近時,系統將趨于穩定,即穩定振動狀態.

而在深色區域外的其他區域,系統總體阻尼比為正值,當系統處于這些區域時,振動系統會逐漸發展至靜止狀態(鎖定區風速范圍以外)或發展至穩定振動狀態(鎖定區風速范圍以內).

當渦振系統處于穩定周期振動狀態時,系統總體阻尼比為零.圖6中紅色矩形為開口斷面在不同約化風速條件下測得的穩定振動幅值.從圖中可以看到:紅色矩形的分布位置與系統總體阻尼比為零的等高線能很好契合,這也進一步驗證了本文系統特性識別方法的準確性.

2.3 附加氣動效應

利用圖5中四次多項式擬合曲線,在不同約化風速條件下,用系統總體阻尼比曲線減去零風速條件下彈性懸掛系統本身阻尼比曲線,即得到渦振發生時系統氣動阻尼比隨系統約化振幅的變化規律.

得到的開口斷面渦振系統氣動阻尼比如圖7所示.可以看到:在鎖定區域內(曲面下凹區域),氣動阻尼隨約化振幅和約化風速的變化規律呈現明顯非線性特性,且氣動阻尼值為負值,這是導致結構在該區域內出現渦激振動的原因.而在鎖定區域外,氣動阻尼比隨約化振幅和約化風速的變化較為平緩,這些區域內的氣動阻尼值更接近于常數,即線性氣動阻尼效應.

圖6 系統總體阻尼比等高線圖

圖7 開口斷面渦振氣動阻尼比

Fig.7 Aeroelastic damping ratio of the open-section bridge deck during VIV

根據圖4結果,渦振系統振動頻率可看作只與約化風速有關,與振幅無關.因此各約化風速條件下的氣動剛度為常數.圖8給出了不同約化風速條件下的氣動剛度值,其中氣動剛度定義為

式中:m為模型質量;ω為渦振系統振動的圓頻率(有風條件);ω0為彈性懸掛體系振動的圓頻率(無風條件).

圖8表明:不同約化風速條件下模型所受氣動剛度效應并不相等,鎖定區后半段模型所受氣動剛度效應大于鎖定區前半段模型所受氣動剛度效應.不同約化風速條件下模型所受氣動剛度效應的不同,表明跨越鎖定區時氣動力與結構運動位移間的相位差發生了明顯變化.

圖8 開口斷面渦振氣動剛度

Fig.8 Aeroelastic stiffness of the open-section bridge deck during VIV

3 結 論

1)識別得到的系統總體阻尼比為零的等高線與試驗測得的結構渦振響應曲線能很好吻合,從而驗證了系統瞬時特性識別方法的準確性.

2)利用開口斷面試驗數據得到的渦振系統附加氣動效應結果表明:氣動阻尼隨結構振幅呈明顯非線性變化規律,且各風速下氣動阻尼隨結構振幅的變化規律并不一致,表明氣動阻尼與約化風速和結構振幅同時相關,而各風速下的氣動剛度隨結構振幅的變化規律則呈線性特性,表明氣動剛度可僅看作與約化風速相關而與結構振幅無關.

3)后續研究可嘗試利用某一數學模型對所得到的附加氣動效應隨結構振幅的變化規律進行模擬,用以計算不同結構參數下的渦振響應.

[1] LI H, LAIMA S, OU J, et al. Investigation of vortex-induced vibration of a suspension bridge with two separated steel box girders based on field measurements[J]. Engineering Structures, 2011, 33(6): 1894-1907.

[2] FUJINO Y, YOSHIDA Y. Wind-induced vibration and control of Trans-Tokyo Bay Crossing Bridge[J]. Journal of Structural Engineering, 2002, 128(8): 1012-1025.

[3] LARSEN A, ESDAHL S, ANDERSEN J E, et al. Storeblt suspension bridge-vortex shedding excitation and mitigation by guide vanes[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2000, 88(2/3): 283-296.

[4] 張文明, 葛耀君, 楊詠昕,等. 帶挑臂箱梁渦振氣動控制試驗[J]. 哈爾濱工業大學學報, 2010,42(12):1948-1952.

ZHANG Wenming, GE Yaojun, YANG Yongxin, et al.Experimental study on aerodynamic control of the vortex induced vibrations of a box girder with projecting slab[J].Journal of Harbin Institute of Technology,2010,42(12):1948-1952.

[5] 董銳, 楊詠昕, 葛耀君. 斜拉橋Π型開口斷面主梁氣動選型風洞試驗[J]. 哈爾濱工業大學學報, 2012, 44(10):109-114.

DONG Rui,YANG Yongxin,GE Yaojun. Wind tunnel test for aerodynamic selection of Π shaped deck of cable-stayed bridge[J].Journal of Harbin Institute of Technology,2012,44(10):109-114.

[6] EHSAN F, SCANLAN R H. Vortex-induced vibrations of flexible bridges[J]. Journal of Engineering Mechanics-Asce, 1990, 116(6): 1392-1411.

[7] GOSWAMI I, SCANLAN R H, JONES N P. Vortex-induced vibration of circular-cylinders: II new model[J]. Journal of Engineering Mechanics-Asce, 1993, 119(11): 2288-2302.

[8] SCANLAN R H. Bridge flutter derivatives at vortex lock-in[J]. Journal of Structural Engineering, 1998, 124(4): 450-458.

[9] LARSEN A. A generalized model for assessment of vortex-induced vibrations of flexible structures[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1995, 57(2): 281-294.

[10]D’ASDIA P, SEPE V, CARACOGLIA L, et al. A model for vortex-shedding induced oscillations of long-span bridges[C]//Proceedings of the 2nd International Structural Engineering and Construction Conference (ISEC-02). Rome: CRC Press, 2003: 2331-2336.

[11]FELDMAN M. Non-linear free vibration identification via the hilbert transform[J]. Journal of Sound and Vibration, 1997, 208(3): 475-489.

[12]FELDMAN M, BUCHER I, ROTBERG J. Experimental identification of nonlinearities under free and forced vibration using the hilbert transform[J]. Journal of Vibration & Control, 2009, 15(10): 1563-1579.

[13]HUANG N E, WU Z. A review on Hilbert-Huang transform: method and its applications to geophysical studies[J]. Reviews of Geophysics, 2008, 46(2): 2008.

[14]NUTTALL A H, BEDROSIAN E. On the quadrature approximation to the Hilbert transform of modulated signals[J]. Proceedings of the IEEE, 1966, 54(10): 1458 - 1459.

[15]HUANG N E, WU Z, LONG S R, et al. On instantaneous frequency[J]. Advances in Adaptive Data Analysis, 2009, 1(2): 177-229.

(編輯 魏希柱)

Identification of aeroelastic effects on bridge decks during vortex-induced vibration

XU Kun1,2, GE Yaojun1, CAO Fengchan1

(1.State Key Laboratory of Disaster Reduction in Civil Engineering(Tongji University), Shanghai 200092, China;2.Institute of Road and Bridge Engineering, Beijing University of Technology, Beijing 100124, China)

To identify the aeroelastic effects on bridge decks during vortex-induced vibration (VIV), a new method is proposed. The procedure of instantaneous identification of the VIV system is firstly proposed. Based on this procedure, the method for extracting the aeroelastic damping and stiffness from the total damping and stiffness of the VIV system is establised. Data of wind tunnel experiments on an open-section bridge deck undergoing vortex-induced resonance is involved for validation and for studying the features of aeroelastic effects on this deck during VIV. The zero-contour of the total damping for the VIV system compares well with the experimental VIV responses, which confirms the accuracy of this method. Based on the results, aeroelastic effects at different wind velocities and amplitudes on this open-section bridge deck are studied in detail, which can offer some valuable information for the subsequent study on mathematical modeling of this phenomenon.

bridge deck; wind tunnel test; vortex-induced vibration; aeroelastic effects; identification method

10.11918/j.issn.0367-6234.2017.03.014

2015-09-24

國家重點基礎研究發展計劃(2013CB036300)

許 坤(1988—),男,博士; 葛耀君(1958—),男,教授,博士生導師

葛耀君,yaojunge@tongji.edu.cn

U441.3

A

0367-6234(2017)03-0086-07

猜你喜歡
風速振動系統
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
中立型Emden-Fowler微分方程的振動性
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
基于GARCH的短時風速預測方法
主站蜘蛛池模板: 亚洲成人一区在线| 亚洲无码精彩视频在线观看 | 999国内精品视频免费| 久久国产V一级毛多内射| 国产精品无码制服丝袜| 九色在线视频导航91| 青青操视频免费观看| 亚洲第一成年免费网站| 91尤物国产尤物福利在线| 国产精品3p视频| 久久精品亚洲热综合一区二区| 狠狠ⅴ日韩v欧美v天堂| 亚洲人成网站在线观看播放不卡| 国产激爽爽爽大片在线观看| 青青青伊人色综合久久| 日韩资源站| 久久a级片| 国产免费黄| 国产欧美日韩视频怡春院| 亚洲国产精品无码久久一线| 精品人妻系列无码专区久久| 国产福利小视频在线播放观看| 日韩免费成人| 88av在线播放| 精品少妇人妻无码久久| 中文字幕日韩视频欧美一区| 国产幂在线无码精品| 五月天综合网亚洲综合天堂网| 色香蕉影院| 久久婷婷国产综合尤物精品| 国产成人精品亚洲日本对白优播| 国产美女视频黄a视频全免费网站| 日日拍夜夜操| 国产色爱av资源综合区| 蜜臀AV在线播放| 成人在线观看一区| 久热中文字幕在线| 久久综合九色综合97婷婷| 日韩成人在线视频| 色哟哟精品无码网站在线播放视频| 99精品热视频这里只有精品7| 四虎永久在线| 久草视频精品| www.狠狠| 亚洲精品动漫| 喷潮白浆直流在线播放| 国产中文在线亚洲精品官网| 亚洲码一区二区三区| 久久天天躁狠狠躁夜夜躁| 国产亚洲欧美日韩在线一区二区三区| 亚洲区欧美区| 日本a级免费| 色成人综合| 亚洲大尺度在线| 欧美日在线观看| 国产视频大全| 久久亚洲精少妇毛片午夜无码 | 曰韩人妻一区二区三区| 国产免费一级精品视频 | 国产偷倩视频| 国产精品林美惠子在线播放| 亚洲资源站av无码网址| 日韩AV无码一区| 一级做a爰片久久毛片毛片| 日韩欧美国产另类| 2020国产在线视精品在| 欧美狠狠干| 日韩黄色精品| 欧美性猛交一区二区三区| 国产精品无码一二三视频| 天天色天天综合| 波多野结衣无码视频在线观看| 色综合a怡红院怡红院首页| 国产在线八区| 日韩精品久久久久久久电影蜜臀| 欧美综合中文字幕久久| 免费看a级毛片| 91av成人日本不卡三区| 国产尤物在线播放| 国产精品无码AV中文| 国产免费高清无需播放器| 国产91全国探花系列在线播放|