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

基于慣性半徑相似岸橋模型的地震試驗

2016-11-23 08:14:54王貢獻胡吉全
振動、測試與診斷 2016年5期
關鍵詞:有限元結構模型

李 哲, 王貢獻, 胡吉全, 王 東

(武漢理工大學物流工程學院 武漢,430063)

?

基于慣性半徑相似岸橋模型的地震試驗

李 哲, 王貢獻, 胡吉全, 王 東

(武漢理工大學物流工程學院 武漢,430063)

岸橋結構進行縮尺模型設計時,各構件(梁)的幾何尺寸可由原形的尺寸按相似比推導出來, 梁截面厚度卻不能和長度按同一比尺縮小,導致梁截面尺寸不能完全相似。筆者采用截面慣性半徑相似的方法控制抗彎剛度設計梁截面尺寸。以該方法為參考制作了一臺岸橋縮尺模型并進行地震試驗,縮尺模型試驗測量值與數值計算結果接近,結構各考察點加速度、應變出現極值的大小以及所對應的時間都比較接近。該縮尺模型能替代完全相似模型,所得試驗值能近似反映實際結構在地震中的動態響應,驗證了截面慣性半徑相似模型設計方法的可行性,為工程中大型鋼結構縮尺模型的設計及后續的地震試驗研究提供了參考。

慣性半徑; 相似理論; 岸橋結構; 縮尺模型; 振動臺試驗

引 言

縮尺模型試驗現在已經成為大型結構動力學研究過程中不可或缺的部分。以結構動力相似理論為依據設計制作與實際結構尺寸相似而按一定比尺縮放的試驗模型[1-2],對模型進行一系列動力試驗,將所得的試驗數據按一定的相似比進行換算,即可得到能近似反映實際結構的動態響應參數。

岸橋模型地震試驗中最主要的問題在于其結構中梁截面的厚度不能按統一相似比縮放。岸橋梁截面厚度平均在10~25 mm之間,如果尺寸相似比定為1∶50,則縮尺模型中梁截面厚度應為0.2~0.5 mm才能滿足梁截面完全相似。0.2~0.5 mm厚度角鋼的焊接在實際加工中是極難實現的,工程上一般采用加厚的角鋼(2~3 mm)進行焊接,這樣得到的縮尺模型和原型就不能完全相似,無法通過縮尺模型預測原型的動態特性。現有的岸橋相似模型設計過程中,為應對梁截面不能完全相似的情況紛紛采取一些措施,如文獻[3]設計制作了一個相似比為1∶15的岸橋畸變模型(其梁截面厚度沒有按統一的相似比例進行縮放),采用有限元預測系數法,得到畸變模型動響應預測系數,實現用畸變模型預測原型結構的動態響應。雖然該方法得到的結果準確可靠,但是畸變預測系數必須經過大量的有限元計算推導,還需要相應的試驗進行驗證,操作起來復雜繁瑣、周期長。文獻[4]設計制作了1∶20的岸橋縮尺模型,采用質量補償法對模型進行修正來應對模型不完全相似的問題,在此基礎之上他們進行了一系列的地震模擬試驗。這套方案根據動力相似原理,對試驗模型進行了質量修正,計算出所需添加的質量,并將質量塊通過焊接、捆綁等方式加在模型上,質量塊的安放位置主要依據試驗人員的經驗,對于安放了質量塊的地方產生局部應力加大以及剛度變化問題未進行考慮,使得部分測量結果產生較大偏差。

現有的研究表明,結構的變形主要是彎曲變形的可選用抗彎剛度作為相似比控制參數,而對于剪切變形為主的結構,選用抗剪剛度作為相似比控制參數更為妥當。大量地震災害分析表明,岸橋結構地震激勵下的破壞形式主要是梁的彎曲變形。筆者經過大量的有限元數值計算以及簡單梁架結構的振動試驗,得出岸橋結構在地震激勵下以立柱的彎曲振動為主,由此可以選用抗彎剛度作為相似比控制參數。現采用截面慣性半徑相似的方法,設計出岸橋結構彎曲梁的截面參數,使岸橋理想試驗模型結構與實際加工出的試驗模型截面慣性半徑誤差達到最小。以岸橋模型相似比尺的量綱分析為基準,結合梁截面慣性矩公式推導出縮尺模型梁截面的尺寸,制作1∶15的岸橋縮尺模型。建立岸橋結構1∶15完全相似縮尺模型的有限元數值模型、以加工出來的試驗縮尺模型為基準的有限元模型以及岸橋原型有限元模型。對試驗縮尺模型進行振動臺地震模擬試驗,配合有限元數值計算來驗證截面慣性半徑相似模型設計方法的正確性。

1 模型相似比尺分析

筆者研究對象為某港口正在使用的一款俯仰式集裝箱岸橋起重機,工作時臂架放下總長達66 m,凈高54 m、整機重量約為630 t。根據振動臺臺面大小(1.5 m×1.5 m)、振動臺最大搭載質量(2 t)以及試驗場所的空間限制, 尺寸相似比取1∶15。簡化了的岸橋結構包括:海陸側門腿立柱、橫梁及撐桿、前后大梁、拉桿等,結構示意圖見圖1。選取與岸橋原型結構相同的材料Q345鋼,其材料彈性模量取值E=206 GPa,泊松比取值ν=0.3,材料密度取ρ=7 850 kg/m3。

圖1 集裝箱岸橋結構模型Fig.1 Structural model of the quayside container crane

結構動力學性能最直接的表現是:頻率、變形、位移[5]。在岸橋結構地震模擬試驗中,縮尺模型相關參量應選取:彈性模量E、密度ρ、幾何尺寸l、時間t、位移u、速度υ、加速度a、重力加速度g、應力σ、頻率ω。在線彈性范圍內[6],采用量綱分析法得到這些參量的函數表達式

(1)

用C來表示岸橋模型與原型之間參量的相似比,如前文所述,幾何比尺為:Cl=1/15。縮尺模型的制作材料和原形相同,所以密度比尺Cρ=1、彈性模量比尺CE=1,將這三項作為基本比尺,推導時間比尺Ct、質量比尺Cm、位移比尺Cu、速度比尺Cv、加速度比尺Ca、頻率比尺Cω、應變比尺Cε、應力比尺Cσ。梁彎曲振動的微分方程為

(2)

其中:A為梁的截面面積;I為梁的截面慣性矩。

由式(2)推導得出:

(3)

其中:CA為截面面積比尺;CI為截面慣性矩比尺。

由上述推導可得岸橋結構動力模型相似關系如表1所示。

表1 縮尺模型相似關系

2 截面慣性半徑相似模型設計

岸橋結構在地震下的響應以梁的彎曲振動為主,梁的彎曲剛度可以通過截面慣性半徑來控制(I=r2A)[7]。岸橋原型梁截面和縮尺模型的梁截面示意圖分別如圖2(a,b)所示。正確設計梁截面參數,使原型梁截面慣性半徑與縮尺模型梁截面慣性半徑的相似誤差達到最小,對提高縮尺模型的精度至關重要。

圖2 岸橋原型梁截面與縮尺模型梁截面Fig.2 Cross-sectional shapes of bending beam components of prototype and scale mode

由截面慣性矩公式可得原型與縮尺模型繞y軸的梁截面慣性矩

(4)

(5)

其中:Iyp,Iym分別為原型與縮尺模型繞y軸的梁截面慣性矩;L,H,t1,t2分別為原型的梁截面長度、寬度、左右厚度及上下厚度;Lm,Hm,tm為縮尺模型的梁截面長度、寬度及厚度。

原型與縮尺模型繞z軸的梁截面慣性矩Izp與Izm也可按截面慣性矩公式給出,形式與式(4),(5)一致。結合截面慣性半徑的公式r2=I/A可以得到

(6)

(7)

其中:ryp,rym,rzp,rzm為原型與縮尺模型繞y,z軸的梁截面慣性半徑;Ap,Am為梁截面面積;Cry,Crz為繞y,z軸的梁截面慣性半徑比尺,均為1/15。

在縮尺模型截面尺寸設計中,首先確定加工鋼板厚度tm,再基于式(6)及(7)計算出梁截面的參數Lm和Hm。由于岸橋結構梁構件多、方程復雜、計算量大,故使用通用數學與工程計算軟件maple編寫計算程序,方程沒有精確解。因此,基于工程經驗在計算前設置好各參數的取值范圍,得到的截面參數Lm和Hm的近似解。

以海測立柱為例,梁截面尺寸為:L=0.900 m,H=1.400 m,t1=0.010 m,t2=0.010 m;梁截面慣性半徑比尺Cry=Crz=1/15;鋼板厚度取tm=0.003 m。將已知參數代入計算程序得到近似解Lm=0.060 24 m,Hm= 0.096 22 m。模型加工精度為毫米級,取Lm= 0.060 m,Hm= 0.096 m。

采用上述設計方法得到的岸橋主要梁截面尺寸如表2所示。(表中Cry′,Crz′為Cry,Crz的倒數,與最初設定值15有細微差別,可忽略不計)

表2 1∶15縮尺模型彎曲梁構件的截面參數

Tab.2 Cross-sectional parameters of bending beam components of 1∶15 scale model

彎曲梁名稱截面參數/mm慣性半徑比尺(倒數)LmHmtmC'ryC'rz前大梁5757314.8214.82后大梁63102314.8315.05海測立柱6096314.9515.61陸測立柱6096314.9515.61海測橫梁63112314.8314.59陸測橫梁63123314.8314.92梯形梁6262315.0915.09

3 岸橋結構地震試驗分析

3.1 模型仿真分析

使用有限元仿真計算驗證截面慣性半徑相似模型設計方法的可行性。按原型尺寸建立岸橋有限元模型M1、按1∶15相似比嚴格縮放的完全相似縮尺模型M2(完全相似縮尺模型)以及一個按截面慣性半徑相似模型設計方法設計截面尺寸的縮尺模型M3(彎曲剛度相似縮尺模型),岸橋結構模型的有限元模型示意圖見圖3。在采用有限元計算前需對有限元模型進行修正,配合縮尺模型(見圖4)模態試驗以及加載試驗、采用參數型修正法對岸橋有限元模型進行修正,從而提高有限元數值模型精度。在有限元模型修正過程中,通過反復調整各部件材料屬性,使得有限元模型與試驗縮尺模型前幾階頻率振型一致、加載后應力值大小接近[8],且應通過結構阻尼試驗為有限元模型設置合理的阻尼參數[9-10]。

圖3 集裝箱起重機有限元模型Fig.3 Finite element model of the container crane

圖4 試驗縮尺模型Fig.4 Scale model of the container crane

首先對有限元模型M1,M2,M3進行模態計算,其計算結果如表3所示,其中模型M2的頻率值經過相似比轉換后與模型M1完全一致,M3頻率值與M2非常接近,但存在著較小的誤差,誤差在可接受的范圍內。模型只比較頻率一般是不全面的,應該還要比較振型,考慮到振型圖占用篇幅較多,且第1振型占其地震反應的主要部分[11],這里采用文字描述岸橋模型前3階振型特征。1階模態主要振型特征:前大梁揚起、門框沿大梁方向彎曲;2階模態主要振型特征:門腿沿大梁方向彎曲;3階模態主要振型特征:后大梁揚起、門框沿大梁方向彎曲。岸橋模型M1,M2前3階振型特征完全一致, M2與M3振型之間存在細微區別,可忽略不計。模型M2與M3頻率、振型產生誤差的主要原因包括:截面慣性半徑不能完全相似;雖然岸橋結構地震主要以彎曲變形為主,但模型的振型中可能包含構件的剪切變形。模型M2,M3前6階固有頻率接近且前3階振型特征相似,說明模型M3可以代替完全相似縮尺模型,即采用截面慣性半徑相似方法在有限元仿真方面是可行的。

表3 岸橋模型固有頻率

Tab.3 The first six frequencies of models

階數f/Hz誤差/%M1M2M3M3與M212.543238.69738.8240.3323.029645.89443.9564.2233.269349.85347.1255.4743.526753.61550.8575.1453.639455.13952.9423.9864.072961.77364.0843.74

進行地震時程分析,取EL-Centro南北向地震加速度記錄(20 s,Δt=0.02 s)、Taft東西向地震加速度記錄(20 s,Δt=0.02 s),加速度峰值分別調整為 0.1g,0.2g,對于1∶15縮尺模型而言,地震載荷的時間軸調整為原型的1/15(20 s調整為1.34 s),而峰值加速度則被調整為1.5g(對應原型的0.1g),3g(對應原型的0.2g)。圖5中A1,A4為加速度考察點,A2,A3為位移考察點,S1~S12為應力應變考察點。

圖5 模型測量點Fig.5 Measurement nodes of the model

經對比分析時程計算結果,發現完全相似模型M2的時程計算結果與岸橋原型M1完全一致,這符合實際情況,也驗證了數值模型的正確性。故后面分析中將模型M2的值作為基準,對比分析時,只需將M3與M2或者試驗值與M2進行對比,這樣可以避免相似比換算時產生的誤差,結果也更加直觀。加速度峰值調整為 0.1g時,EL,Taft波下考察點A1的加速度時程曲線如圖6,7所示 (注:實線代表模型M2,虛線表示M3。在EL波下考察點A1的加速度方向為小車運行方向,即南北方向;在Taft波下考察點A1的加速度方向為大車運行方向,即東西方向)。

圖6 考察點A1的加速度時程曲線(EL 0.1 g)Fig.6 Acceleration time-history curve at A1(EL 0.1 g)

圖7 考察點A1的加速度時程曲線(Taft 0.1 g)Fig.7 Acceleration time-history curve at A1(Taft 0.1 g)

加速度峰值調整為 0.2g時,EL,Taft波下考察點A1的加速度時程曲線如圖8~9所示。

圖8 考察點A1的加速度時程曲線(EL 0.2 g)Fig.8 Acceleration time-history curve at A1(EL 0.2 g)

圖9 考察點A1的加速度時程曲線(Taft 0.2 g)Fig.9 Acceleration time-history curve at A1(Taft 0.2 g)

考察點A2在地震波EL,Taft加速度峰值調整為 0.2g時位移時程曲線如圖10~11所示。

圖10 考察點A2的位移時程曲線(EL 0.2 g)Fig.10 Displacement time-history curve at A2(EL 0.2 g)

圖11 考察點A2的位移時程曲線(Taft 0.2 g)Fig.11 Displacement time-history curve at A2(Taft 0.2 g)

觀察圖6~11可以發現,由模型M3計算得到的地震響應與M2的結果十分接近,說明在有限元時程計算中采用截面慣性半徑相似方法設計的有限元模型M3代替完全相似模型M2是可行的。在相同地震激勵、不同的加速度峰值調整下考察點加速度時程曲線形狀相似,0.2g時加速度峰值變化為0.1g時的1.5~2倍,與前期預計的情況相符。

3.2 振動臺模型試驗

本試驗旨在證明采用截面慣性半徑相似方法設計的試驗縮尺模型能代替完全相似試驗縮尺模型(由于完全相似試驗縮尺模型無法加工制造,所以與試驗值比較的對象為經過修正的有限元模型M2的計算值)。試驗位于武漢理工大學港口裝卸實驗室,所采用的振動臺可進行水平方向、垂直方向的任意單方向或組合雙向振動試驗,臺面長寬均為1.5 m,可提供加速度范圍為±50 m·s-2,最大位移±200 mm,最大速度均為0.8 m·s-1,最大承載力2 t,頻率范圍0.1~100 Hz,3 dB帶寬。試驗中,2個加速度傳感器、2個位移傳感器、12個應變片(S1~S12)被布置在模型門架結構上如圖5所示。按照M3的梁截面參數制造一臺試驗用縮尺試驗模型,對其進行振動臺地震模擬試驗來驗證此相似模型是否可以代替完全相似模型,加工出的試驗縮尺模型見圖4。

采用文獻[12]中介紹的模態測試方法來獲取岸橋1∶15縮尺試驗模型的固有頻率。M2,M3有限元模態計算結果和縮尺模型試驗值及相應的誤差如表4所示。

表4 試驗模型固有頻率

為節約篇幅,只給出考察點A4的加速度時程曲線 (在EL地震波下,加速度峰值調整為 0.2g)。如圖12所示,按截面慣性半徑相似方法加工出的模型能夠較好地預測原型的地震響應。圖中實線為模型M2計算值,虛線為縮尺模型試驗實測值。

觀察對比完全相似模型

M2的計算結果和縮尺

圖12 考察點A4的加速度時程曲線Fig.12 Acceleration time history curve at A4

模型試驗實測值,在不同地震載荷下、不同加速度峰值調整下岸橋結構上各考察點加速度、應變出現極值的大小以及所對應的時間都比較接近。最大加速度出現在岸橋大梁前后兩端,這與試驗所測結果一致,也符合實際地震情況下岸橋結構加速度響應;最大的應變主要集中在前大梁與拉桿連接處以及后大梁的中部位置;原型結構的應變分布情況與表5中的應變分布情況接近,說明按照M3截面參數設計加工出的縮尺模型能得到正確可靠的試驗數據。

表5 不同測點的最大動態應變

Tab.5 Maximum dynamic strain at various nodes

測點應變(E×10-6)M2試驗值誤差/%試驗值與M2S282.578.64.72S4156.8169.37.97S6206.7211.52.32S9138.3125.39.39S10271.3275.41.51S12307.9313.11.69

完全相似模型M2計算值與試驗實測值之間存在微小差異的原因有:a.梁的截面慣性半徑相似不可能完全被滿足,與設計的15有偏差;b.有限元模型與試驗模型存在一些細小的差別;c.傳感器測量誤差;d.雖然岸橋結構地震主要以梁的彎曲變形為主,但結構的振型中可能包含構件的剪切變形。采用截面半徑相似方法設計加工的試驗模型的可行性得到了進一步的驗證,采用該方法設計加工出的縮尺模型能代替完全相似模型。

4 結 論

1) 采用截面慣性半徑相似方法對岸橋梁截面參數進行設計,加工出的1∶15岸橋縮尺模型能代替完全相似模型,準確地預測出原型結構的動態特性和地震響應。

2) 岸橋縮尺試驗模型地震試驗實測值中頻率、振型、加速度、位移以及應變與數值計算結果有較高的相似度,最大誤差為9.39%,滿足工程需要,該方法能為后續的研究提供精確可靠的試驗縮尺模型。

[1] 陳喆,陳國平. 相似理論和模型試驗的結構動響應分析運用[J]. 振動、測試與診斷, 2014,34(6):995-1001.

Chen Zhe, Chen Guoping. Research of dynamics response based on similarity theory and model test [J]. Journal of Vibration, Measurement & Diagnosis, 2014,34(6):995-1001.(in Chinese)

[2] 張瑩,孫廣俊,李鴻晶.開洞填充墻對混凝土框架柱地震損傷影響分析 [J]. 振動、測試與診斷,2014, 34(5):932-938.

Zhang Ying, Sun Guangjun, Li Hongjing. Numerical simulation on seismic damage of reinforced concrete frame columns considering the influence of infill walls with opening [J]. Journal of Vibration, Measurement & Diagnosis, 2014, 34(5):932-938. (in Chinese)

[3] 李哲,胡吉全,王東.地震載荷作用下岸橋結構單參數畸變相似模型研究[J]. 振動與沖擊. 2014,33(20):174-179.

Li Zhe, Hu Jiquan, Wang Dong. Distortion model of container cranes subjected to seismic load[J]. Journal of Vibration and Shock, 2014, 33(20): 174-179.(in Chinese)

[4] Jacobs L D, Des Roches R, Roberto T L. Seismic behavior of a jumbo container crane including uplift[J]. Earthquake Spectra, 2011, 27(3):745-773.

[5] 何川,汪洋,方勇,等. 土壓平衡式盾構掘進過程的相似模型試驗[J]. 土木工程學報,2012,45(2):162-169.

He Chuan, Wang Yang, Fang Yong, et al. Similarity model test of earth-pressure-balanced shield tunneling process[J]. China Civil Engineering Journal, 2012,45(2):162-169. (in Chinese)

[6] 金玉龍,吳天行. 集裝箱碼頭岸橋結構的動力相似分析與試驗驗證[J]. 上海交通大學學報:自然科學版, 2012, 46(10): 1609-1615.

Jin Yulong, Wu Tianxing. Dynamic similarity analysis and experimental verification on a quayside container crane[J]. Journal of Shanghai Jiaotong University:Science, 2012, 46(10): 1609-1615.(in Chinese)

[7] Jin Yulong, Li Zengguang. Theoretical design and experimental verification of a 1/50 scale model of a quayside container crane[J]. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2012, 226(6):1644-1662.

[8] 蘇忠亭,徐達,李曉偉,等. 步兵戰車車體結構有限元模型修正 [J]. 振動、測試與診斷,2014, 34(6):1148-1155.

Su Zhongting, Xu Da, Li Xiaowei, et al. Finite-element model updating for infantry combat vehicle car-body[J]. Journal of Vibration, Measurement & Diagnosis, 2014, 34(6): 1148-1155.(in Chinese)

[9] 張教超,王敏慶,馬建剛. 結構阻尼的聲衰減時間測量方法 [J]. 振動、測試與診斷,2013, 33(2):330-333.

Zhang Jiaochao, Wang Minqian, Ma Jiangang. Approach for determining structural damping by measurement of sound attenuation time[J]. Journal of Vibration, Measurement & Diagnosis, 2013, 33(2): 330-333.(in Chinese)

[10]張安付,閆孝偉,盛美萍,等. 自由阻尼結構損耗參數的換算方法 [J]. 振動、測試與診斷,2015, 35(3):572-578.

Zhang Anfu, Yan Xiaowei, Sheng Meiping, et al. Conversion method for loss factor of unconstrained damping structures[J]. Journal of Vibration, Measurement & Diagnosis, 2015, 35(3): 572-578.(in Chinese)

[11]劉世忠. 系桿拱橋的橫向振動特性分析 [J]. 振動、測試與診斷,2015, 35(2):245-250.

Liu Shizhong. Analysis of transverse vibration frequency of tied arch bridge[J]. Journal of Vibration, Measurement & Diagnosis, 2015, 35(2): 245-250.(in Chinese)

[12]孫熙平,王元戰,趙炳皓. 環境激勵下高樁碼頭物理模型模態實驗[J].振動、測試與診斷,2013,33(2):263-268,340.

Sun Xiping, Wang Yuanzhan, Zhao Binghao. Modal experiment of physical model for high-piled wharf under ambient excitation[J]. Journal of Vibration, Measurement & Diagnosis, 2013, 33(2): 263-268,340.(in Chinese)

doi:10.16450/j.cnki.issn.1004-6801.2016.05.025

10.16450/j.cnki.issn.1004-6801.2016.05.024

國家自然科學基金資助項目(51275369)

2015-07-08;

2015-09-02

U448.21; TH13

李哲,男,1986年5月生,博士生。主要研究方向為機械工程。曾發表《地震載荷作用下岸橋結構單參數畸變相似模型研究》(《振動與沖擊》2014年第33卷第20期)等論文。

E-mail:172042756@qq.com

簡介:王貢獻,男,1976年12月生,博士、副教授。主要研究方向為結構動力學。

E-mail:wgx@whut.edu.cn

猜你喜歡
有限元結構模型
一半模型
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
3D打印中的模型分割與打包
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 国产精品视频999| 日本午夜影院| 久久精品视频亚洲| 国产综合精品一区二区| 91在线中文| 秋霞午夜国产精品成人片| 国产天天射| 亚洲男人天堂2018| 亚洲自拍另类| www.99精品视频在线播放| 真人高潮娇喘嗯啊在线观看| 亚洲精品卡2卡3卡4卡5卡区| 在线观看亚洲精品福利片| 亚洲综合国产一区二区三区| 亚洲婷婷丁香| 国产一级毛片高清完整视频版| 白浆视频在线观看| 免费大黄网站在线观看| 国产xxxxx免费视频| 国产成人8x视频一区二区| 婷婷色狠狠干| 国产成人乱码一区二区三区在线| 九九免费观看全部免费视频| 国产福利一区二区在线观看| 99久久亚洲精品影院| 久久婷婷五月综合97色| 国产sm重味一区二区三区| 国产免费看久久久| 国产欧美视频综合二区| 黄色三级网站免费| 人人澡人人爽欧美一区| 91久久精品国产| 国产成在线观看免费视频| 国产在线视频自拍| 四虎成人精品在永久免费| 亚洲天堂久久新| 国产精品视频公开费视频| 国产精品自拍合集| 国产在线八区| 国产精品污污在线观看网站| 欧美人人干| 亚洲精品第一页不卡| 亚洲欧洲AV一区二区三区| 黑人巨大精品欧美一区二区区| 久久99热这里只有精品免费看| 日韩午夜片| 97青草最新免费精品视频| 久久人人爽人人爽人人片aV东京热| 久久 午夜福利 张柏芝| 欧美日韩北条麻妃一区二区| 青青草91视频| 精品久久综合1区2区3区激情| 尤物特级无码毛片免费| 午夜少妇精品视频小电影| 成人毛片在线播放| 国产一级毛片高清完整视频版| 亚洲欧美日韩中文字幕在线一区| 国产精品久久久精品三级| 中文无码精品a∨在线观看| 青青操国产| 色综合天天综合| 国产区在线看| 国产免费精彩视频| 国产精品区视频中文字幕| 国产97色在线| 免费jizz在线播放| 永久免费无码成人网站| 国产精品三级av及在线观看| 国产伦精品一区二区三区视频优播| 欧美日韩动态图| 国产在线无码av完整版在线观看| 亚洲高清中文字幕在线看不卡| 亚洲综合第一页| 国产欧美视频在线| 2021最新国产精品网站| 国产后式a一视频| 亚洲色婷婷一区二区| 国产午夜人做人免费视频| 91青青视频| 中文字幕1区2区| 国产精品短篇二区| 四虎影视8848永久精品|