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

大跨網殼地震反應時Rayleigh阻尼構建方法比較

2017-12-12 02:35:32潘旦光李雪菊
哈爾濱工業大學學報 2017年12期
關鍵詞:模態結構方法

潘旦光,程 業,李雪菊

(1.北京科技大學 土木系,北京 100083;2.土木工程防災國家重點實驗室(同濟大學),上海 200092)

大跨網殼地震反應時Rayleigh阻尼構建方法比較

潘旦光1,2,程 業1,李雪菊1

(1.北京科技大學 土木系,北京 100083;2.土木工程防災國家重點實驗室(同濟大學),上海 200092)

為研究Rayleigh阻尼系數參考頻率選取方法對結構地震反應的影響,以一個長85.2 m,寬61.8 m的網殼為例,對比分析了基于地震波頻譜特性選取第二個參考頻率的方法、Idriss方法、傳統方法、優化方法所得Rayleigh阻尼系數引起頂點位移和基底剪力計算誤差.對于地震波的頻譜特性,討論了擬加速度反應譜、擬速度反應譜、位移反應譜的峰值頻率和形心點頻率的統計范圍.數值計算結果表明:基于擬速度反應譜、位移反應譜峰值頻率和形心頻率以及Fourier譜的峰值頻率所得Rayleigh阻尼將使計算結果偏??;用小于基本周期范圍內的擬加速度反應譜的峰值頻率或形心點頻率作為第二個參考頻率計算誤差較?。划數卣鸩ǖ淖吭筋l率大于基頻的結構時,Idriss方法用平滑化Fourier譜的卓越頻率計算參考頻率效果更好;優化方法可以直接得到Rayleigh阻尼系數,避免了人為選取參考頻率的任意性,且計算精度高.

地震反應;Rayleigh阻尼;卓越頻率;參考頻率;優化方法

在地震反應過程中,阻尼是影響反應結果的重要因素[1-2].為進行強震下結構的彈塑性時程反應分析,需要構建相應的阻尼矩陣.在各種阻尼模型中,Rayleigh阻尼由于計算簡便而得到廣泛應用.

采用矩陣表達,Rayleigh阻尼矩陣C可表達為質量矩陣M和剛度矩陣K的線性組合,即

C=αM+βK,

(1)

式中α和β分別為質量和剛度比例阻尼系數.α和β這兩個系數可以通過指定兩個參考頻率(ωi和ωj)及其阻尼比(ζi*和ζj*)進行計算.顯然,Rayleigh阻尼除兩階參考頻率的阻尼比外,其他各階模態的阻尼比有一定的誤差存在.因此,所構造的Rayleigh阻尼是否合理依賴于選取的兩個參考頻率是否合理.在結構地震反應分析時,Rayleigh阻尼參考頻率的選取有三類常用方法.

第一類方法是直接從結構的動力特性角度選擇參考頻率.Idriss等[3]直接利用體系基頻ω1建立阻尼矩陣,此時所構造的Rayleigh阻尼高估了除基頻以外的所有阻尼,從而導致結構動力反應偏小[4].對于簡單的建筑結構,常用兩個低階自振頻率作為參考頻率[5].但對復雜結構而言,對結構有顯著貢獻的模態數目多達幾百階[6].此時,如果用兩個低階模態建立Rayleigh阻尼矩陣將使結構高階模態的阻尼比偏大,從而導致低估結構的地震反應.Chopra[7]定性的指出,選擇的兩階參考頻率應使對結構反應有顯著貢獻模態的阻尼比取值合理.Clough等[8]建議ωi=ω1,ωj從對結構動力反應有顯著貢獻的高階振型中選取.Youssef等[9]的計算結果表明,最優參考頻率的階數隨土層深度變化而變化.

第二類方法是令ωi=ω1,第二個參考頻率ωj由地震波的頻譜特性進行確定.Hudson等[10]令ωj為地震波卓越頻率ωe和結構基頻的奇數倍.即

ωj=nω1,

(2)

式中n為大于ωe/ω1的奇數.樓夢麟等[11]建議以輸入地震波加速度反應譜的峰值頻率作為ωj進行深覆蓋土層的地震反應計算.

第三類方法采用優化算法直接得到Rayleigh阻尼系數.楊大彬等[12]以振型的峰值應變能系數作為權重函數,建立了加權最小二乘法.潘旦光等[13-14]分別以結構峰值位移誤差和基底剪力誤差最小為目標函數,提出了一種求解Rayleigh阻尼系數的優化求解方法.

為比較上述三類Rayleigh阻尼系數計算方法的特點及對結構動力反應的影響.本文將以大跨網殼結構水平方向和豎向地震反應為例,討論不同方法所得阻尼矩陣對網殼頂點位移及基底剪力的影響.除此以外,本文還從Rayleigh阻尼構建角度討論了地震波主要頻率成分及其選取范圍的問題.

1 大跨網殼有限元模型及輸入地震波

1.1 有限元模型

結構計算模型如圖1所示,上部為85.2 m長,61.8 m寬單層網殼結構,下部由8組支架支撐整個結構.結構平面關于x軸對稱,關于y軸不對稱,z為豎向坐標.結構主要構件如下:網殼內部桿件和環梁分別為300 mm×50 mm×8 mm×12 mm和1 000 mm×600 mm×30 mm×30 mm的空心方管,下部斜撐和底柱分別為300 mm×20 mm和800 mm×35 mm的空心圓管.所有桿件材料為Q235,采用梁單元進行建模.模型的模態特征如表1所示.表中ry和rz分別表示y和z方向的振型參與質量比,sy和sz分別表示y和z方向的累積振型參與質量比.部分典型模態圖如圖2所示.并設各階模態的阻尼比為0.02.

圖1 網殼有限元模型

圖2 部分模態圖

1.2 輸入地震波

為比較不同類型地震波對Rayleigh阻尼系數計算的影響,選用表2中的3條不同場地類型地震波分別作為柱根部的水平和豎向地震輸入.輸入地震波的加速度時程如圖3所示,加速度時程的幅值統一調整為0.35 m/s2.

表1 固有頻率及振型參與質量比

表2 地震波

圖3 地震波加速度時程

2 Rayleigh阻尼系數計算方法

由模態分析可知,單層網殼結構在y方向(水平)地震輸入時,前45階模態的振型參與質量即超過90%,第一階模態的反應貢獻具有絕對統治地位.對于z方向(豎向)地震輸入時,前205階模態的振型參與質量超過90%,第三階模態的反應貢獻最大,但不具有統治地位.對于水平地震反應,結構的顯著貢獻模態少;而豎向地震反應,結構的顯著貢獻模態多,且沒有具有絕對統治地位的模態.同時,在水平地震反應時,結構的第一個顯著貢獻模態的頻率小于大部分地震波的卓越頻率,而豎向地震反應時,結構的第一個顯著貢獻模態的頻率高于部分地震波的卓越頻率,因此,討論y方向和z方向的地震反應,相當于討論了兩種結構類型下Rayleigh阻尼系數的計算問題.不同方法計算阻尼系數的差別在于兩個參考頻率選取的不同,為此,下面討論單層網殼結構y方向和z方向地震反應時三類方法所得的參考頻率及對結構地震反應的影響.

2.1 第一類方法

對于y方向地震輸入,一種是根據經驗選擇ωi=ω1和ωj=ω19(記為i=1 &j=19),另一種是直接令ωi=ω1和ωj=ω2(i=1 &j=2).對于z方向地震輸入,分別考慮兩種組合ωi=ω3和ωj=ω45(i=3 &j=45)及ωi=ω1和ωj=ω2(i=1 &j=2).并在后面的討論中將這種方法稱為傳統方法.

2.2 第二類方法

在第二類方法中,y方向和z方向地震輸入的第一個參考頻率分別選為ω1和ω3.第二個參考頻率可基于反應譜和Fourier譜進行選擇.圖4為3條地震波的位移、擬速度和擬加速度反應譜.

在利用反應譜選擇Rayleigh阻尼第二個參考頻率時,常直接選取反應譜峰值頻率[15].但這種方法存在的一個問題是:當結構的基頻高于反應譜峰值頻率時,構建的Rayleigh阻尼將使所有的高階模態阻尼比大于真實的阻尼比,從而使計算結果偏小,這一點在位移反應譜中尤其明顯,因為位移反應譜的峰值常出現在長周期.譬如,自振周期在0~6 s范圍內,El Centro、Northridge和Tianjin波的位移反應譜峰值頻率分別為0.352、0.410和0.568 Hz,這些峰值頻率都小于基頻,此時,都無法構造出合理的Rayleigh阻尼.事實上,反應譜反映的是在一個確定地震波作用下不同自振周期單自由度體系的最大反應,然后,以一定組合規則得到結構的反應.結構高階模態的周期Ti都是小于基頻周期T1的,即Ti

圖4 地震波反應譜

考慮到多自由度體系是一系列單自由度體系反應的疊加,因此,部分學者將地震波的主要頻率定義為反應譜形心的頻率[15].反應譜的形心頻率fRg-i定義如下

(3)

式中:Si(ζ,T)為反應譜,i分別取位移、擬速度和擬加速度.顯然fRg-i是積分范圍TR的函數.圖5為不同地震波反應譜fRg-i隨積分區間變化的曲線.由圖可知,隨著積分區間的增加,形心點的頻率逐步降低.這是長周期部分反應譜影響的必然結果.由于只有[0,T1]區間的反應譜參與結構反應的計算,因此,計算形心頻率的反應譜范圍定義為[0,T1]更合理.表4和表5中列出擬加速度反應譜形心頻率fRg-a、擬速度反應譜形心頻率fRg-v和位移反應譜形心頻率fRg-d都是指[0,T1]基本區間不同地震波反應譜的形心頻率.

圖5 反應譜形心頻率

圖6為不同地震波Fourier譜.Fourier譜表明地震波中含有哪些頻率分量,及哪些頻率分量振幅大.因此,將分量振幅最大的頻率稱為卓越頻率,并記為Fourier譜峰值頻率fF.實際地震波的Fourier譜呈劇烈起伏的鋸齒狀,為避免Fourier譜中個別尖刺的影響,常將平滑化方法所得的峰值頻率作為地震波卓越頻率.本文將平滑化后Fourier譜的峰值頻率記為fP.圖6中的平滑化曲線是采用矩形脈沖窗的結果,窗的帶寬取為1.2 Hz.

在得到地震波的卓越頻率后,可采用式(2)計算Rayleigh阻尼系數.當ωe取為Fourier譜的峰值頻率fF時,稱為I-1方法;當ωe取為平滑化Fourier譜的峰值頻率fP時,稱為I-2方法.且當ωe<ω1時,令n=1.

圖6 地震波Fourier譜

2.3 第三類方法

對于第三類方法,以文獻[14]提出的約束優化解法作為對比.該方法首先建立代數方程組:

(4)

根據以上討論,共形成12種Rayleigh阻尼系數的計算方法,其中傳統方法在算例比較時,采用兩種頻率組合.表3列出了各種計算方法的簡化名稱和所對應的類型.下面討論各種方法所得Rayleigh阻尼對單層網殼結構地震反應的影響.

表3 Rayleigh阻尼系數計算方法

3 數值計算結果

3.1 計算誤差公式

為比較不同方法的計算誤差,進行結構的線彈性地震反應分析.同時,以前300階模態振型分解時程分析方法的計算結果為精確解,記為r*,采用Rayleigh阻尼模型所得的近似解記為r,則Rayleigh阻尼模型計算結果的相對誤差為

(5)

在m條地震波作用下,各反應量的平均誤差可采用以下兩式進行統計:

(6)

3.2 不同Rayleigh阻尼模型計算誤差

El Centro、Northridge和Tianjin三條地震波分別沿y方向(水平)和z方向(豎向)輸入下,前述12種方法的參考頻率如表4和表5所示.在y方向地震輸入時,僅統計y方向的地震反應,z方向地震輸入時,僅統計z方向的地震反應.圖1中A點的位移計算誤差和基底剪力誤差如表6~9所示.El Centro波作用下A點位移反應的精確解及其Fourier譜如圖7所示,基底剪力的精確解及其Fourier譜如圖8所示.表和圖中的uAy和uAz分別表示A點的y方向和z方向的位移.Fy和Fz分別表示y方向和z方向的基底剪力.

表4y方向地震輸入下的參考頻率

Tab.4 Reference frequencies under the y direction seismic input Hz

表5z方向地震輸入下的參考頻率

Tab.5 Reference frequencies under the z direction seismic input Hz

表6 y方向地震輸入下uAy的相對誤差

表7 z方向地震輸入下uAz的相對誤差

表8 y方向地震輸入下基底剪力Fy的相對誤差

表9 z方向地震輸入下基底剪力Fz的相對誤差

圖7 El Centro波作用下A點的位移反應時程及Fourier譜

Fig.7 Time histories and Fourier spectra of pointAdisplacement under the El Centro wave

圖8 El Centro波作用下基底剪力反應時程及Fourier譜

Fig.8 Time histories and Fourier spectra of base reactions under the El Centro wave

由計算結果可看出:

2) 對于基底剪力Fy和Fz,以fRv、fRd、fRg-v、fRg-d為第二參考頻率的計算結果都偏小.對比反應譜的峰值頻率和形心頻率可以發現fRa>fRv>fRd,以及fRg-a>fRg-v>fRg-d.這是由于反應譜在大于2 Hz區域為加速度敏感區,小于0.3 Hz為位移敏感區,中間區段為速度敏感區[7].而由圖8基底剪力的Fourier譜可知,高階模態對基底剪力的反應也有顯著貢獻.譬如f19=2.905 Hz和f45=4.443 Hz模態的反應分別對Fy和Fz有顯著影響,而基于速度反應譜的fRv、fRg-v和位移反應譜的fRd、fRg-d的第二個參考頻都遠小于f19和f45,導致這兩階模態的阻尼比偏大而低估了高階模態的反應.基于加速度反應譜峰值頻率和形心頻率的計算誤差相對較小.

3)fF作為第二個參考頻的計算結果偏小.這是由于Fourier譜和擬速度反應譜類似,因此,fF和fRv的計算結果也類似.以平滑化的峰值頻率fP作為第二個參考頻的計算結果要略優于fF的計算結果,但無法改變計算結果偏小的特點.同時,當fF或fP小于結構的基頻時,直接用卓越頻率作為第二個參考頻率的計算誤差更大.這是由于這種方法除基頻的阻尼比外,其余模態的阻尼比都大于精確解而導致計算結果偏小.因此,對于地震波卓越頻率小于結構基頻的地震波,采用地震波卓越頻率作為參考頻率不合理.

4)I-2方法的平均誤差小于I-1方法,這表明采用平滑化Fourier譜所得卓越頻率進行Idriss方法[10]第二個頻率的計算更合理.而且,由于fF易受Fourier譜中個別低頻尖刺的影響,而使計算結果離散性大.同時,當地震波的卓越頻率小于結構的基頻時,此時,Idriss方法實際上就是基于基頻建立阻尼矩陣,Rayleigh阻尼高估了高階模態的阻尼比,由此使結構的地震反應偏小,因此,當地震波的卓越頻率小于結構的基頻時,I-1和I-2方法都是不合理的.當地震波的卓越頻率大于結構的基頻時,可采用I-2方法建立Rayleigh阻尼.

5)對于優化方法,無論是y方向地震輸入還是z方向地震輸入,由于優化算法中考慮了結構動力特性、地震波頻率特性的影響,因此,對于兩個方向的地震輸入和所有地震波的計算結果誤差都小且穩定.而且,優化方法是直接得到Rayleigh阻尼系數,避免了人為選擇兩階參考頻率的任意性,適用于不同工程結構的Rayleigh阻尼構建.

6)對于傳統方法,本文所選的y方向和z方向的頻率組合i=1 &j=19和i=3 &j=45實際上和優化分析方法的頻率基本相同,因此,對于富有經驗的計算人員,直接選取兩階合理的參考頻率用于Rayleigh阻尼計算也是可行的.但是,任意選擇前兩階模態進行Rayleigh阻尼計算所得計算結果偏小,且絕對誤差最大,是不合理的計算方法.

4 結 論

在大型復雜結構的非線性地震反應分析時,常需建立Rayleigh阻尼矩陣進行直接積分法計算.由于參與結構振動的模態多且復雜,因此,如何選取合理的參考頻率是一個需要仔細斟酌的事情.本文以一個長85.2 m,寬61.8 m的大跨屋蓋為例,對比分析了12種國內外Rayleigh阻尼系數計算方法對計算結果的影響,由理論分析和數值計算可得出以下結論:

1)對于第一類方法,工程技術人員如果對結構的動力反應有充分認識,可直接指定兩階參考頻率.但是,任意選擇前兩階模態進行Rayleigh阻尼計算,易造成計算結果偏小,是不合理的計算方法.

2)對于第二類方法,建議采用[0,T1]區間的擬加速度反應譜峰值頻率或形心頻率、I-2方法作為Rayleigh阻尼系數計算的第二個參考頻率.

3)優化方法所得Rayleigh阻尼系數是綜合考慮結構的動力特性、輸入地震波頻譜特性的綜合結果,直接得到Rayleigh阻尼系數,避免了基于經驗指定兩階頻率的任意性,且計算誤差較小,適合于工程結構的計算與分析.

[1] 沈聚敏,周錫元,高小旺. 抗震工程學[M].北京: 中國建筑工業出版社,2002.

SHEN Jumin, ZHOU Xiyuan, GAO Xiaowang. Earthquake engineering [M]. Beijing: China Architecture & Building Press, 2002.

[2] 翟長海,謝禮立,張茂花. 阻尼對工程結構等延性地震抗力譜的影響分析[J]. 哈爾濱工業大學學報,2007,38(10): 1705-170.

ZHAI Changhai, XIE Lili, ZHANG Maohua. Influence analysis damping on constant-ductility seismic resistance spectra for seismic design of structures [J]. Journal of Harbin institute of technology, 2007, 38(10) : 1705-1708.

[3] IDRISS I M, LYSMER J, HWANG R, et al. Quad 4: a computer program for evaluating the seismic response of soil structures by variable damping finite element procedures[R]. Berkeley: University of California, 1973.

[4] 鄒德高,徐斌,孔憲京.瑞利阻尼系數確定方法對高土石壩地震反應的影響研究[J].巖土力學,2011,32(3):797-803.

ZOU Degao, XU Bin, KONG Xianjing. Study of influence of different methods for calculating Rayleigh damping coefficient on high earth-rock dam seismic response[J]. Rock and Soil Mechanics, 2011,32(3):797-803.

[5] 周國良,李小軍,劉必燈,等. 大質量法在多點激勵分析中的應用、誤差分析與改進[J]. 工程力學,2011,28(1): 48-54.

ZHOU Guoliang, LI Xiaojun, LIU Bideng, et al. Error analysis and improvements of large mass method used in multi-support seismic excitation analysis [J]. Engineering Mechanics, 2011, 28(1): 48-54.

[6] 潘旦光,靳國豪,高莉莉. 大跨斜拉橋Rayleigh阻尼系數的約束優化解[J]. 振動與沖擊,2014,33(16): 34-41.

PAN Danguang, JIN Guohao, GAO Lili. Constraint optimal solution of Rayleigh damping coefficients for long-span cable-stayed bridges[J]. Journal of Vibration and Shock, 2014, 33(16): 34-41.

[7] CHOPRA A K. Dynamics of structures: theory and applications to earthquake engineering [M]. New Jersey: Englewood Cliffs, Prentice-Hall, 1995.

[8] CLOUGH R W, PENZIEN J. Dynamics of structures[M]. New York: Mc-Graw Hill Inc, 1993.

[9] YOUSSEF M A, HASHASH D P. Viscous damping formulation and high frequency motion propagation in nonlinear site response analysis[J]. Soil Dynamics and Earthquake Engineering, 2001, 22: 611-624.

[10]HUDSON M, IDRISS I M, BEIKAE M. User manual for QUAD4m: A computer program to evaluate the seismic response of soil structures using finite element procedures and incorporating a compliantbase[D]: Berkeley: University of California, 1994.

[11]樓夢麟,邵新剛.應用通用程序計算深覆蓋土層地震反應的幾個問題[J]. 振動與沖擊,2015, 34(4): 63-68, 109.

LOU Menglin, SHAO Xin’gang. Several problems in seismic response calculation of soil layer with deep deposit using general software[J]. Journal of Vibration and shock, 2015, 34(4): 63-68, 109.

[12]YANG Dabin, ZHANG Yigang, WU Jinzhi. Computation of Rayleigh damping coefficients in seismic time-history analysis of spatial structures [J]. Journal of the International Association for Shell and Spatial Structures, 2010, 51(2): 125-135.

[13]潘旦光. 直接確定Rayleigh阻尼系數的一種優化方法[J]. 工程力學,2013,30(9): 16-21.

PAN Danguang. An optimization method for the direct determination of Rayleigh damping coefficients[J]. Engineering Mechanics, 2013, 30(9): 16-21.

[14]PAN D G,CHEN G D, GAO L L. A constrained optimal Rayleigh damping coefficients for structures with closely-spaced natural frequencies in seismic analysis[J]. Advances in Structural Engineering, 2017, 20(1):81-95.

[15]樓夢麟,邵新剛.深覆蓋土層Rayleigh阻尼建模問題討論[J]. 巖土工程學報,2013, 35(7):1272-1279.

LOU Menglin, SHAO Xingang. Discussion on modeling issues of Rayleigh damping matrix in soil layer with deep deposit[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(7):1272-1279.

(編輯趙麗瑩)

ComparisonstudyonRayleighdampingconstructionmethodsforlong-spanreticulatedshellseismicresponse

PAN Danguang1,2, CHENG Ye1, LI Xueju1

(1.Department of Civil Engineering, University of Science and Technology Beijing, Beijing 100083, China;2.State Key Laboratory of Disaster Reduction in Civil Engineering (Tongji University), Shanghai 200092, China)

To analyze the effects of Rayleigh damping coefficients reference frequency selection method on structural seismic response, some methods to estimate Rayleigh damping coefficients including the frequency contents-based second reference frequency, IDRISS method, traditional method and the optimization method were compared by the errors of top displacement and base shear with a 85.2 m length, 61.8 m width reticulated shell. As for the frequency contents of earthquake waves, the statistics ranges of the peak and centroid frequencies of pseudo-acceleration-, pseudo-velocity- and displacement response spectrum were discussed. The numerical results show that: the second reference frequency specified by the peak and centroid frequencies of the pseudo-velocity- and displacement response spectrum as well as peak frequency of the Fourier spectrum make the response results smaller than the exact ones. The estimation error is reasonable, which the second reference frequency is equal to the peak and centroid frequencies of the pseudo-acceleration response spectrum in the range below the fundamental frequency. It is suitable that the reference frequency is estimated by the smoothing peak frequency of the Fourier spectrum in the IDRISS method when the predominant frequency of earthquake wave is greater than the fundamental frequency of structure. Optimization method can get the Rayleigh damping coefficient directly, eliminate the arbitrariness of reference frequency and have high calculation accuracy.

seismic response; Rayleigh damping; predominant frequency; reference frequency; optimization method

10.11918/j.issn.0367-6234.201611130

TU311.3; TU352.1+1

A

0367-6234(2017)12-0045-08

2016-11-28

土木工程防災國家重點實驗室開放基金(SLDRCE15-01)

潘旦光(1974—),男,研究員,博士生導師

潘旦光,pdg@ustb.edu.cn

猜你喜歡
模態結構方法
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
國內多模態教學研究回顧與展望
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 成人午夜亚洲影视在线观看| 伊人查蕉在线观看国产精品| 久久亚洲中文字幕精品一区| 2019国产在线| 国产69囗曝护士吞精在线视频| 国产导航在线| 久久精品午夜视频| 国产成人久久777777| 2020精品极品国产色在线观看| 成人自拍视频在线观看| 免费激情网址| 思思热精品在线8| 国产精品流白浆在线观看| 欧美一级高清片久久99| 精品国产一区二区三区在线观看 | 嫩草影院在线观看精品视频| 自拍亚洲欧美精品| 中国一级特黄视频| 亚洲国产天堂久久综合| 一边摸一边做爽的视频17国产| 国产精品乱偷免费视频| 久久亚洲欧美综合| 久久人体视频| a亚洲视频| 中文字幕 日韩 欧美| 欧美www在线观看| 91人妻日韩人妻无码专区精品| 精品一区二区三区水蜜桃| 米奇精品一区二区三区| 一区二区自拍| 国产精品青青| 五月婷婷导航| 日韩午夜片| www.av男人.com| 亚洲狼网站狼狼鲁亚洲下载| 国产91精选在线观看| 色国产视频| 亚洲三级色| 日韩123欧美字幕| 九九视频免费在线观看| 国产www网站| 亚洲中文字幕无码爆乳| 无码精品国产dvd在线观看9久| 国产一级毛片yw| 毛片a级毛片免费观看免下载| 在线观看国产小视频| 69综合网| 亚洲国产91人成在线| 99性视频| 午夜电影在线观看国产1区| 欧美色亚洲| 亚洲一区二区精品无码久久久| 欧美日韩国产系列在线观看| 午夜不卡视频| 日韩专区欧美| 亚洲av色吊丝无码| 国产视频大全| 99国产精品国产| 99伊人精品| 国产区福利小视频在线观看尤物| 在线看片中文字幕| 国产人碰人摸人爱免费视频| 亚洲色欲色欲www网| 国产成人毛片| 久久久久中文字幕精品视频| 四虎影视库国产精品一区| 亚洲男人的天堂在线观看| 亚洲欧洲日韩综合| 国产不卡在线看| 亚卅精品无码久久毛片乌克兰| a亚洲天堂| 国产一区成人| 2048国产精品原创综合在线| 亚洲免费成人网| 尤物亚洲最大AV无码网站| 一本二本三本不卡无码| 亚洲国产欧美国产综合久久| 波多野结衣一二三| 亚洲香蕉久久| 成人午夜视频免费看欧美| 久久无码高潮喷水| 国产青青草视频|