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

大跨斜拉橋Rayleigh阻尼系數約束優化解

2014-09-08 03:30:40潘旦光靳國豪高莉莉
振動與沖擊 2014年16期
關鍵詞:模態優化結構

潘旦光,靳國豪,高莉莉

(北京科技大學 土木系,北京 100083)

我國規范[1]規定,進行斜拉橋、懸索橋、單跨跨徑150 m以上等特殊橋梁地震反應分析時可采用時程分析法、多振型反應譜法及功率譜法。用直接積分法進行時程分析時必涉及阻尼矩陣的建立。由于形成阻尼機理復雜,影響因素多[2],因此無法直接利用構件尺寸、材料性質直接構建阻尼矩陣,而采用構造的方法。由于工程結構阻尼比較易測量,且能反應結構宏觀阻尼特性,因此常用阻尼比構建阻尼矩陣。阻尼可分為復阻尼及粘滯阻尼。復阻尼耗能與迫振頻率無關,這與結構試驗結果相符。阻尼為復數,易于用頻域進行分析。對非線性分析只能采用等效線性化方法,無法進行真非線性分析。與復阻尼相比,用時域內粘滯阻尼更方便[3-4]。

目前已有諸多粘滯阻尼矩陣的構造方法[5-7],其中最具代表性的為Rayleigh阻尼、Caughey阻尼及疊加振型阻尼。董軍等[8-9]將Rayleigh阻尼與疊加振型阻尼結合,構造出任意階模態阻尼比等于精確值的阻尼矩陣。但常用的為質量矩陣與剛度矩陣線性組合的Rayleigh阻尼矩陣,可表達為

C=αM+βK

(1)

式中:α,β分別為質量、剛度比例阻尼系數;M,K分別為質量、剛度矩陣。

計算α,β一般先進行結構模態分析,后指定兩階參考模態(ωi,ωj)的阻尼比等于已知值(ζi*,ζj*)。α,β計算公式為

(2)

由式(1)得各階模態阻尼比為

ζn=α/2ωn+βωn/2

(3)

進行橋梁設計時,雖事先無法獲知模態阻尼比的精確值,但可據已建橋梁阻尼比經驗值設各階模態精確值為常量。如鋼筋混凝土橋梁阻尼比為0.05。因此式(3)所得阻尼比僅在兩階參考頻率處阻尼比等于精確值,而其它頻率阻尼比存在一定誤差。由于阻尼比會影響結構動力反應,因此合理指定兩階參考頻率成為時程分析關鍵。結構基頻為影響結構動力反應最主要模態。規范[1]中建議第一個參考頻率選擇基頻。而對第二個參考頻率選擇無明確方法。樓夢麟等[10-11]通過對大跨拱橋計算認為參考頻率的選擇會顯著影響結構反應,應選振型參與系數較大模態作為參考頻率;通過對深覆土層模型分析認為計算Rayleigh阻尼系數需考慮輸入地震波頻譜特性影響。為避免人為選擇兩階參考頻率所致任意性,劉紅石[12]提出的最小二乘法計算α,β可使阻尼比誤差在平方意義上最小,但并未考慮各階模態對動力反應貢獻的差異,無法保證對結構反應有顯著貢獻模態的阻尼比取值合理。楊大彬等[13-14]針對網殼結構地震反應分析,建立基于多參考振型的加權最小二乘法用于計算Rayleigh阻尼系數。潘旦光[15-16]基于振型疊加反應譜理論,詳細討論式(2)中兩最優參考頻率及荷載空間分布、頻譜特性及結構動力特性關系,并建立Rayleigh阻尼系數的無約束優化求解方法。但對橋梁結構而言,與結構所受激勵輸入方向對應的第1階整體振動模態為控制結構動力反應的關鍵模態之一,通常以該階模態為計算Rayleigh阻尼系數的第1個參考頻率。因此,Rayleigh阻尼系數計算實際為合理選取第2個參考頻率問題。大跨橋梁參與振動模態數較多,頻率密集,不同于普通建筑結構。為此,本文以大跨橋梁為分析對象,以結構峰值位移誤差最小為目標,建立求解Rayleigh阻尼系數的優化方法。并將基頻阻尼比等于精確值作為約束條件,形成等價于指定最優第2參考頻率的約束優化求解方法。并以840 m斜拉橋為例,說明Rayleigh阻尼系數約束優化解的變化規律及計算方法精度。

1 Rayleigh阻尼系數約束優化解

地震輸入作用下多自由度體系強迫振動方程[17]為

(4)

通過模態分析獲得結構前N階頻率ωn及模態φn(n=1,2,…,N)。由模態疊加反應譜法理論知,基于Rayleigh阻尼與精確阻尼比,分別計算所得第n階模態對第k自由度最大位移反應ukn為

(5)

Sd(ζn,ωn)=

基于平方和開平方原理,第k自由度位移反應誤差為

(7)

為簡化計算,將反應譜函數用1階Taylor級數展開,得

(8)

將式(8)代入式(7),整理得

(9)

(10)

采用Lagrange乘子法求式(10)約束條件下式(9)最小值,令

(11)

將式(11)分別對X,λ求導,并令相應導數為零得代數方程組為

(12)

2 參數研究

為說明公式的計算精度及影響因素,以圖1懸臂梁為例進行討論。梁長6 m,彈性模量20 GPa,線密度600 kg/m,截面轉動慣量1.6×10-3m4。將梁分為10個單元,用集中質量模型,體系共10階模態,其中前4階模態分析結果見表1。

圖1 懸臂梁計算模型

表1 懸臂梁模態分析結果

設體系基礎運動豎向加速度為

(13)

激振頻率考慮三種情況:ω=0.7ω1, 6.0ω1,16.0ω1,即分別考慮激振頻率小于基頻、接近ω2、接近ω3。在任意激振頻率下第n階模態反應譜及導數為

Sd(ζ,ωn)=

(14)

(15)

2.1 Rayleigh阻尼模型計算精度

阻尼為結構固有特性,在線彈性范圍內阻尼比與外部荷載無關。而Rayleigh阻尼為近似阻尼,因此為使構建的阻尼矩陣動力反應計算誤差小,應使對結構動力反應有顯著貢獻模態的阻尼比誤差越小越好。不同荷載頻譜特性不同,Rayleigh阻尼合理參考頻率與荷載有關。以含所有模態在精確阻尼比下的模態疊加法計算結果為精確解,記r*。采用Rayleigh阻尼模型計算所得近似解記為r,Rayleigh阻尼模型計算結果相對誤差為

(16)

設懸臂梁所有頻率精確阻尼比為1%,表2、表3為簡諧荷載作用下用Rayleigh阻尼矩陣進行計算所得頂點A的豎向位移及支座B處彎矩相對誤差。計算中本文方法的φkn=φn(A)用前4階模態參與優化計算。為與傳統計算方法進行比較,用式(2)計算Rayleigh阻尼系數時指定兩階參考頻率考慮三種組合,分別為i=1&j=2,i=1&j=3 ,i=2 &j=3。由表2、表3看出:① 由本文方法所得Rayleigh阻尼矩陣的反應誤差均小于或等于傳統計算方法,計算精度良好;② 在不同迫振頻率下結構動力反應的顯著貢獻模態不同,當ω/ω1= 0.7,6.0時,結構動力反應由前兩階模態控制。以i=1 &j=2構建阻尼矩陣計算誤差最小,而ω/ω1=16.0時,以i=1 &j=3構建阻尼矩陣的計算誤差最小。由此可知,計算Rayleigh阻尼系數最優參考頻率應據荷載頻譜特性作相應調整。

表2 A點位移相對誤差(%)

表3 B點彎矩相對誤差(%)

2.2 優化參考頻率影響因素

為說明本文方法所得Rayleigh阻尼系數特點,不同迫振頻率下所得阻尼系數及相應阻尼比見表4。由表4看出,ω/ω1=0.7,6.0時,第2個優化參考頻率j=2;而當ω/ω1=16.0時,第2個優化參考頻率j=3。結合表2、表3知,本文方法所得Rayleigh阻尼系數為綜合考慮結構動力特性及輸入荷載頻譜的結果,使對結構有顯著貢獻模態的阻尼取值更合理,可直接獲得Rayleigh阻尼系數,無需人為選擇參考頻率,便于應用。

表4 優化Rayleigh阻尼系數及相應阻尼比

3 工程概況及輸入地震波

對簡單結構而言,最優兩參考頻率較易判斷;但對大跨橋梁而言,結構顯著貢獻模態數量較多,人為方法選擇合理參考模態較困難。為驗證本文方法計算精度,以840 m斜拉橋為例,分析地震作用下優化Rayleigh阻尼系數的合理性及計算精度。大橋為60 m+120 m+480 m+120 m+60 m雙塔雙索半漂浮體系鋼混混合斜拉橋。塔高155.1 m,橋面板以上高103 m。主梁為寬28.5 m、高3 m的單室封閉鋼箱梁。塔柱與混凝土主梁用C50混凝土,鋼箱梁用Q345qc。塔座、邊墩及輔助墩用C40混凝土,承臺用C30混凝土。拉索用112根Φj15.24 mm鋼絞線。用單脊梁式方法建立有限元模型。塔墩及主梁采用三維梁單元,索用桿單元,建立有限元模型見圖2。

圖2 斜拉橋有限元模型

有限元模型整體坐標系x方向為順橋向,y方向為豎向,z方向為橫橋向,模型共1 149個節點、1 196個單元。設橋梁各階模態阻尼比精確值為0.02,以x方向地震輸入為例,分析塔頂(A點)、主梁跨中(B點)x方向水平位移及橋墩(C點)處剪力及彎矩在不同阻尼模型下結構動力反應特點與計算精度。其它方向地震反應分析方法類似。

式(9)計算中權重系數wkn涉及第k自由度模態位移φkn,即優化的參考自由度問題。理論上參考自由度可選擇任意自由度;但對結構而言,若結構最大位移誤差小,則其余自由度動力反應的誤差通常也較小。而斜拉橋橋塔最大位移位于橋塔頂端,橋身最大位移位于橋梁跨中,故計算中參考自由度分別取A、B點,即x方向地震輸入時分別以φkn=φn(Ax),φkn=φn(Bx)為權重系數的參考自由度,并比較其對計算結果影響。

基于圖2有限元模型,選4條不同類型地震波作為輸入,將地震波峰值加速度統一調整為0.1 g,分析不同Rayleigh阻尼系數計算方法引起結構動力反應誤差。地震波加速度時程與7個阻尼比(ζ=0.005, 0.01, 0.02, 0.03, 0.05, 0.10, 0.20)的反應譜見圖3。實際地震輸入的反應譜為不規則曲線,無法建立位移反應譜顯式表達式。統計分析結果[18]表明,位移反應譜與阻尼比對數間具有線性關系。即對某確定地震波,自振頻率ωn體系的反應譜隨阻尼比ζ變化規律可表示為

Sd(ζ,ωn)=an(ωn)+bn(ωn)ln100ζ

(17)

式中:an(ωn),bn(ωn)為擬合參數,可用最小二乘法計算。

(18)

4 斜拉橋計算結果

4.1 自振特性

對斜拉橋而言,拉索應力對結構動力特性及動力反應有顯著影響。為考慮重力對拉索應力影響及拉索幾何剛度,斜拉橋模態分析分兩步:① 計算重力及初始應力下結構應力;② 進行有預應力的模態分析。斜拉橋主要振動模態見表5。x方向累積振型參與質量見圖4。由各模態振型參與質量、模態特征看出,第1階模態為順橋向振動第1個顯著貢獻模態。式(12)計算時順橋向約束模態取第1階。

圖3 地震波加速度時程及位移反應譜

表5 斜拉橋模態分析結果

圖4 斜拉橋累積振型參與質量

4.2 Rayleigh阻尼系數變化規律

由式(9)知,Rayleigh阻尼系數優化計算結果與參與優化計算模態數N有關。順橋向優化Rayleigh阻尼系數隨模態數N的變化見圖5。結合圖4計算結果知,在累積振型參與質量小于70%時(前128階模態),隨模態數的增加使累積振型參與質量發生顯著變化的模態處Rayleigh阻尼系數亦會發生顯著變化。累積振型參與質量超過70%后,部分高階振型參與質量較大(第538、539階),但Rayleigh阻尼系數未在此處發生相應變化。由于振型參與質量系數是用于度量各模態對基底剪力貢獻大小的指標,部分高階局部振型對結構基底剪力有顯著影響,但塔頂、主梁位移反應主要由低階模態控制,而Rayleigh阻尼系數優化計算方程本質為基于位移反應的優化方程,因此當累積振型參與質量超過70%后,Rayleigh阻尼系數優化結果趨于穩定。與高層建筑等結構Rayleigh阻尼系數變化規律有較大不同。

4.3 優化Rayleigh阻尼系數及優化參考頻率

不同地震波激勵下Rayleigh阻尼系數優化解穩定計算結果見表6。將所得Rayleigh阻尼系數代入式(3),即得各階模態阻尼比,除作為約束條件的模態外,第2個等于精確阻尼比對應的固有頻率階數稱為優化參考頻率j′。由表6看出,① 本文所得優化參考頻率不一定位于結構振型參與質量較大固有頻率上。因除基頻外,橋梁結構無具有顯著統治地位模態,最優參考頻率為綜合考慮多階模態動力反應結果;② 輸入地震波不同最優參考頻率也不同,說明本文結構最優參考頻率考慮輸入地震波頻譜特性影響;③ 以主梁跨中位移為優化目標的優化參考頻率值小于以塔頂為優化目標的計算結果。此因為主梁位移由少數低階模態控制(前28階模態),而塔頂位移頻率成分相對豐富(前52階模態)。因此,優化方程自動識別出的優化目標位移顯著貢獻模態有所不同。

圖5 順橋向地震輸入Rayleigh阻尼系數優化解

表6 優化Rayleigh阻尼系數及優化參考頻率

4.4 Rayleigh阻尼模型反應峰值計算誤差

由圖4知,當模態個數達539階時,體系x方向累積振型參與質量超過90%;620階模態后累積振型參與質量超過99%。為此本文以620階模態振型分解法計算所得結構動力反應量作為精確解。對傳統Rayleigh阻尼計算方法,選i=1 &j=56,i=1 &j=108及i=1 &j=539三種組合方式進行計算。本文方法考慮φkn=φn(Ax),φkn=φn(Bx)兩種情況。塔頂A點、主梁跨中B點在不同地震波輸入下水平位移峰值計算誤差見表7。El Centro波地震輸入作用下塔頂A點、主梁跨中B點水平位移時程及反應譜見圖6、圖7。

由計算結果看出,① 由誤差的平均值看,傳統Rayleigh阻尼計算方法與本文方法位移反應誤差均小于5%,均滿足工程需要。由位移FFT譜知,A、B點水平位移均由第1階模態反應控制。第1階模態阻尼比等于精確解時,結構位移反應誤差即可控制在合理范圍內。② 第2個參考頻率越大(如i=1 &j=539),兩參考頻率之間模態阻尼比越小于精確解,所得計算結果越大于精確解。本文所用優化參考頻率使對結構動力反應有影響部分模態阻尼小于精確解,部分大于精確解,由此達到位移反應誤差為最小目的。用本文方法計算所得結果誤差均小于傳統方法,表明計算精度良好。③ 由精確解的Fourier譜知,A點水平位移主要頻率成分小于2 Hz,B點水平位移主要頻率成分小于1.2 Hz。約束優化算法中考慮各模態對總位移貢獻影響,即表6中φkn=φn(Ax)的第2個參考頻率大于φkn=φn(Bx)的原因。

表7 A、B點水平位移反應峰值相對誤差(%)

圖6 A點水平位移及Fourier譜

不同地震波作用下橋墩C點剪力Fx及彎矩Mz計算誤差(Fx,Mz為整體坐標系下計算結果) 見表8。

El Centro波順橋向地震輸入作用下C點剪力Fx與彎矩Mz時程及反應譜見圖8、圖9。由計算結果看出,① 用傳統方法計算的三種組合計算誤差平均值均大于5%。且隨第2個參考頻率增大誤差增大。此因對該橋墩剪力、彎矩貢獻最大的兩階模態為第28階(頻率1.123 Hz)、第56階(頻率3.058 Hz)。本文所選三種傳統Rayleigh阻尼計算組合均使第28階模態阻尼比小于精確解。第2個參考頻率越大該階模態阻尼比越小,導致時程分析結果較精確解大。② 本文方法雖以位移反應誤差為目標函數獲得Rayleigh阻尼系數,但仍可兼顧內力反應的頻譜特性,使剪力、彎矩計算誤差均小于傳統方法,計算精度良好。③ 比較φkn=φn(Ax)與φkn=φn(Bx)兩組統計結果,以φkn=φn(Ax)為權重函數的計算結果大于精確解,以φkn=φn(Bx)時的結果小于精確解。此因φkn=φn(Bx)時所得第2個參考頻率位于第28階模態附近,而Rayleigh阻尼將使高階振動由于人為原因造成的阻尼增大而被消除。導致對結構剪力、彎矩有顯著影響的第56階模態反應小于精確解,從而所得內力計算結果小于精確解。而φkn=φn(Ax)優化所得第2階參考頻率位于第28、56階模態之間,所得Rayleigh阻尼將高估第28階模態振動同時又低估第56階模態振動,故誤差相對較小且計算結果大于精確解。當結構動力反應小于精確解時,所得計算結果較不安全。因此以塔頂位移為參考自由度進行順橋向地震激勵下Rayleigh阻尼系數的計算更合理。

表8 C點剪力、彎矩峰值相對誤差(%)

圖8 C點剪力及Fourier譜

5 結 論

本文針對橋梁結構基頻模態為結構動力反應顯著貢獻模態特點,提出Rayleigh阻尼系數的約束優化分析方法。利用懸臂梁簡諧振動動力反應,討論荷載迫振頻率對構建合理Rayleigh阻尼影響,結論如下:

(1) Rayleigh阻尼系數約束優化解可綜合考慮荷載頻譜特性及結構動力特性影響,直接獲得Rayleigh阻尼系數,使對結構有顯著貢獻模態的阻尼比合理。

(2) 大跨斜拉橋塔頂及主梁跨中位移主要由少數低階模態控制。用傳統Rayleigh阻尼計算方法分析時,所得結果計算誤差可控制在5%以內;橋墩的剪力、彎矩因其頻率成分豐富,部分高階振型模態對結構的動力反應同樣貢獻顯著;采用人為方法選擇兩階參考模態時,選取不當將引起較大計算誤差。

(3) 隨計算模態個數的增加Rayleigh阻尼系數約束優化解所得α,β趨于穩定。大跨橋梁的計算所用模態個數對應的累積振型參與質量超過70%時,可得α,β的穩定解。

(4) 本文方法所得Rayleigh阻尼系數計算誤差在統計意義上小于傳統方法。可避免因人為選擇參考頻率帶來的任意性,適合工程結構的計算與分析??紤]順橋向振動時,以橋身位移為優化目標所得Rayleigh阻尼系數易使結構動力反應小于精確解,安全性較低;以塔頂位移為參考自由度時所得結構動力反應誤差較小且大于精確解,因此塔頂位移作為參考自由度更合理。

[1] JTG/T B02-01-2008,公路橋梁抗震設計細則[S].

[2] Jeary A P. Damping in structures [J]. Journal of Wind Engineering and Industrial Aerodynamics,1997,72(1):345-355.

[3] 樓夢麟,潘旦光. 滯后阻尼在土層時域分析中的應用[J]. 同濟大學學報, 2004, 32(3): 281-285.

LOU Meng-lin, PAN Dan-guang. Hysteretic damping application in time domain analysis of soil layer[J]. Journal of Tongji University, 2004, 32(3): 281-285.

[4] 何鐘怡,廖振鵬,王小華. 關于復阻尼理論的幾點注記[J]. 地震工程與工程振動,2002, 22(1):1-6.

HE Zhong-yi, LIAO Zhen-peng, WANG Xiao-hua. Some notes on theory of complex damping[J]. Earthquake Engineering and Engineering Vibration, 2002, 22(1): 1-6.

[5] Luco J E. Anote on classical damping matrices[J]. Earthquake Engineering and Structural Dynamics, 2008, 37(4): 615-626.

[6] Adhikari S. Dampingmodelling using generalized proportional damping [J]. Journal of Sound and Vibration, 2006, 293(1/2): 156-170.

[7] 黃宗明,白紹良,賴明. 結構地震反應時程分析中的阻尼問題評述[J]. 地震工程與工程振動,1996, 16(2): 95-105.

HUANG Zong-ming, BAI Shao-liang, LAI Ming. Review on the damping in earthquake response time-history analysis of structures[J]. Earthquake Engineering and Engineering Vibration, 1996, 16(2): 95-105.

[8] 董軍,鄧洪洲,王肇民. 結構動力分析阻尼模型研究[J]. 世界地震工程,2000, 16(4): 63-69.

DONG Jun, DENG Hong-zhou, WANG Zhao-min. Studies on the damping models for structural dynamic time history analysis[J]. World Information on Earthquake Engineering, 2000, 16(4): 63-69.

[9] 克拉夫R,彭津J,著.王光遠,譯.結構動力學(第二版)[M]. 北京:高等教育出版社,2006.

[10] 樓夢麟,張靜. 大跨度拱橋地震反應分析中阻尼模型的討論[J]. 振動與沖擊, 2009, 28(5): 22-26.

LOU Meng-lin, ZHANG Jing. Discussion on damping models for seismic response analysis of long-span bridge [J]. Journal of Vibration and Shock, 2009, 28(5): 22-26.

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

LOU Meng-lin, SHAO Xin-gang. Discussion on modeling issues of rayleigh damping matrix in soil layer with deep Ddeposit [J]. Chinese Journal of Geotechnical Engineering, 2013, 35(7): 1272-1279.

[12] 劉紅石. Rayleigh阻尼比例系數的確定[J]. 噪聲與振動控制, 1999(6): 21-22.

LIU Hong-shi. Determination of rayleigh damping scale coefficient[J]. Noise and Vibration Control,1999(6):21-22.

[13] Yang D B, Zhang G Y G, Wu J Z. 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.

[14] 楊大彬, 張毅剛, 吳金志. 基于多參考振型的Rayleigh阻尼系數計算方法在單層柱面網殼中的應用[J]. 空間結構, 2011,17(3): 8-15.

YANG Da-bin, ZHANG Yi-gang, WU Jin-zhi. Application of the multiple reference modes based computation method of Rayleigh damping coefficients in single-layer cylindrical latticed shell [J]. Spatial Structures, 2011, 17(3): 8-15.

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

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

[16] 潘旦光. 地震反應分析中Rayleigh阻尼系數的優化解[J]. 工程力學,2013, 30(11):15-20.

PAN Dan-guang. An optimization solution for Rayleigh damping coefficients in seismic response analysis[J]. Engineering Mechanics, 2013, 30(11):15-20.

[17] Chopra A K. Dynamics ofstructures: theory and applications to earthquake engineering[M]. New Jersey: Englewood Cliffs, Prentice-Hall, 1995.

[18] Newmark N M, Hall W J. Earthquake spectra and design,earthquake engineering research institute[R]. California: Berkeley, 1982:29-37.

猜你喜歡
模態優化結構
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
國內多模態教學研究回顧與展望
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
基于HHT和Prony算法的電力系統低頻振蕩模態識別
主站蜘蛛池模板: 国产成人综合日韩精品无码首页| 婷婷亚洲最大| 一级毛片免费的| 亚洲日韩久久综合中文字幕| 在线综合亚洲欧美网站| 国产地址二永久伊甸园| 亚洲无限乱码| 亚洲国产成人综合精品2020| 国产人妖视频一区在线观看| 高潮毛片无遮挡高清视频播放| AV不卡无码免费一区二区三区| 免费xxxxx在线观看网站| 国产91av在线| 国产精品性| 国产在线拍偷自揄拍精品| 国产一区二区精品高清在线观看| 91在线国内在线播放老师| 一级毛片在线免费视频| 成人亚洲视频| 成人一级黄色毛片| 亚洲一级毛片在线观| 亚洲天堂精品在线| 欧美日韩国产成人高清视频| 成人免费黄色小视频| jizz在线观看| 国产91在线免费视频| 国内自拍久第一页| 国产精品护士| 香蕉伊思人视频| 制服丝袜一区| 国产午夜不卡| 午夜无码一区二区三区在线app| 青青草原国产精品啪啪视频| 国内精品免费| 91精品伊人久久大香线蕉| 在线国产91| 亚洲国产清纯| 国产午夜看片| 欧美日韩精品一区二区视频| 囯产av无码片毛片一级| 国产正在播放| 干中文字幕| 精品亚洲麻豆1区2区3区| 久久免费看片| 无码日韩视频| 国产精品久久自在自2021| 91年精品国产福利线观看久久| 免费一级毛片不卡在线播放| 国产永久在线视频| 九九九精品成人免费视频7| 黄色网站在线观看无码| 成人午夜视频免费看欧美| 久久黄色免费电影| 五月激情婷婷综合| 老司机午夜精品网站在线观看| 国产美女无遮挡免费视频网站| 青青草91视频| 1级黄色毛片| 91福利国产成人精品导航| 五月天综合网亚洲综合天堂网| 国产欧美专区在线观看| 日韩123欧美字幕| 国产十八禁在线观看免费| 精品国产成人高清在线| 欧美天堂在线| 九九九国产| 视频一区视频二区中文精品| 国产精品免费p区| 天天操天天噜| 国产99视频免费精品是看6| 夜夜操国产| 亚洲中文字幕国产av| 99视频在线观看免费| 美女被操91视频| 99偷拍视频精品一区二区| 精品无码一区二区在线观看| 欧美另类精品一区二区三区| 欧美色图第一页| 国产女人喷水视频| 亚洲乱码在线播放| 手机在线看片不卡中文字幕| P尤物久久99国产综合精品|