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

豎向混合結構阻尼矩陣近似解耦計算的誤差分析

2015-10-14 13:35:54黃維錢江梁飛飛周知
中南大學學報(自然科學版) 2015年4期
關鍵詞:模態結構

黃維,錢江,梁飛飛,周知

?

豎向混合結構阻尼矩陣近似解耦計算的誤差分析

黃維,錢江,梁飛飛,周知

(同濟大學土木工程防災國家重點實驗室,上海,200092)

將豎向混合結構等效成兩自由度結構模型,對該兩自由度模型進行不同子結構參數取值下的近似解耦計算,分析其模態阻尼矩陣對角占優程度和近似解耦計算結構響應誤差情況。研究結果表明:在不同子結構特性下,模態阻尼矩陣對角占優指數和位移誤差指數的衰退特征具有相似性。因此在豎向混合結構動力分析時,可以采用對角占優指數來判定采用近似解耦計算的可行性。

豎向混合結構;非比例阻尼矩陣;模態阻尼矩陣;對角占優;解耦誤差

近年來,高層建筑設計為滿足功能的多樣性以及造型新穎性等需求,采用了不少體型復雜、內部空間多變的建筑方案,鋼?混凝土豎向混合結構體系的應用日漸增多。鋼?混凝土豎向混合結構由下部鋼筋混凝土結構和上部鋼結構串聯組成,其沿豎向變化的抗側剛度能較好的適應結構變形需求。例如:上海環球金融中心大廈內部的核心筒采用了豎向混合的方案,其79層以下核心筒為鋼筋混凝土筒體,為減輕結構自重、增加延性,79層以上核心筒用內置鋼框架的鋼筋混凝土筒體,而在95層以上則改變為空間鋼桁架筒體形式。還有中國民生銀行大廈原是一幢地下2層,地上35層,高度135 m的鋼管混凝土框架?核心筒高層建筑。2005年結構改造中,采用鋼結構加層至45層,總高度達到175.8 m[1]。在這類結構中,上、下部分采用不同材料,各部分耗能機制不同,使得整體結構的阻尼特性難以確定,進而影響結構抗震設計。對于由2種或2種以上的不同材料串聯組成的豎向混合結構,由于不同材料在結構的不同部分提供的能量損失機制差別很大,結構特性會與兩部分相對構成有關。采用整體結構一致的單一阻尼比在具體取值上也存在困難。假設單一材料結構體系的阻尼特性是確定的、可解耦的,根據子結構概念,可以組裝形成混合結構整體運動方程,此時的整體阻尼矩陣一般不再是可解耦的,即阻尼矩陣不再具有關于實模態矩陣的正交性,經典的模態疊加法也不再適用。但由于模態疊加法用模態坐標代替原系統的物理坐標,能使原系統的運動方程解耦,且能用較少的低階模態坐標的響應來獲得較高的近似解,因此許多學者對非比例阻尼結構體系采用模態疊加法進行了研究。Foss[2]提出的狀態空間法,建立狀態方程進行復模態求解,從而在復數領域實現了非比例阻尼結構體系的解耦,但是該方法需要求解擴階矩陣的特征值和特征向量,其計算量較大,且涉及復數運算,不便于實際工程設計使用。Ma等[3?4]在非比例阻尼結構體系復模態動力特性的基礎上,在狀態空間采用坐標變換,得到非比例阻尼結構體系的實空間解耦方法。該方法使得復系數矩陣變換到實數領域,但由于假定條件的引入,使得結構的位移響應和速度響應不是獨立的。從便于工程應用的角度,人們更傾向于實數域的近似解法。采用等效阻尼比的方法得到廣泛的研究,Roesset等[5]基于結構損耗能量和總儲存能量的比值提出了結構的等效阻尼比,該表述是基于能量的阻尼性能宏觀描述,其物理意義明確。Hwang等[6]提出了模態應變能比,其將能量的形式特指為結構的彈性應變能,即結構損耗的彈性應變能與儲存的彈性應變能的比值,該表達式物理意義清楚且形式簡單,因此被廣泛的應用[7]。但采用統一的等效阻尼比改變了不同結構部分的耗能情況。Rayleigh[8]認為如果模態阻尼矩陣非對角線元素相對于對角線元素來說非常的小,可以忽略非對角線元素對方程進行解耦求解,即近似解耦方法,其實質也是一種等效阻尼比方法。Knowles[9]對非比例阻尼系統采用近似解耦法進行了誤差研究,認為在一定條件下采用近似解耦法能得到較好結果。Prandinaa等[10?11]認為在模態阻尼矩陣對角占優的情況下,采用該方法求得的結構動力響應誤差較小。然而Morzfeld等[12]認為模態阻尼矩陣對角占優不能完全確定近似解耦方法產生較小的誤差, 造成這個現象的原因是由于非對角元素對應的特征向量對結構的響應有較大的貢獻。桂國慶等[13]研究了這種近似解耦方法引起的近似誤差,導出了誤差無窮范數,并得到了精確解所在的范圍,認為近似解耦方法有一定的適用性,但該方法計算繁瑣。忽略模態阻尼矩陣非對角線元素的近似解耦方法與所導致的結構響應誤差之間的關聯關系是一個值得研究的問題。本文作者將豎向混合結構等效成兩自由度結構模型,每一自由度的動力特性由對應子結構的主頻確定。對該兩自由度模型進行不同子結構特性下采用近似解耦計算時模態阻尼矩陣對角占優程度和位移響應誤差分析,以此判斷系統的可解耦性及近似程度。

1 豎向混合結構非比例阻尼矩陣 形成

Clough等[14?15]從理論上構造了適合混合結構的阻尼矩陣,整體阻尼矩陣可由分塊阻尼矩陣合成,各子塊阻尼矩陣可按不同模型確定。由于當前應用最廣泛的阻尼模型仍是Rayleigh模型,因此,基于Rayleigh模型構建混合結構的子結構阻尼矩陣的做法具有良好的工程應用基礎。圖1所示為典型的鋼?混凝土豎向混合結構,該類結構都可以等效為兩自由度結構模型,其上部子結構(用字母表示)等效質量、等效剛度及等效阻尼矩陣分別為,對應的阻尼比為ξ;下部子結構(用字母表示)等效質量、等效剛度及等效阻尼矩陣分別為,對應的阻尼比為ξ。整體結構的質量矩陣和剛度矩陣可以由上、下部結構組裝起來,如式(1)和(2)所示;子結構的動力特性可由式(3)求得。

根據Rayleigh阻尼模型,等效兩自由度結構模型的前2階頻率分別為1和2時,則各部分子結構的Rayleigh阻尼系數可由式(4)計算得到;整體結構的阻尼矩陣可采用式(5)和式(6)計算得到。

2 近似解耦計算

一般黏性多自由度系統受迫振動的運動方程為

其中:和分別為結構的質量矩陣、阻尼矩陣和剛度矩陣。對于豎向混合結構,結構的阻尼矩陣可由式(6)得到。

對應的無阻尼系統特征方程為

可解得個特征值ω和對應的特征向量。采用質量歸一時,可得:

其中:=[1,2,…,]為特征向量矩陣;= diag[21,22,…,2]為譜矩陣。

對方程(7)進行=的模態坐標變換,并左乘矩陣T,令()=T(),可得:

其中:為模態阻尼矩陣。

當阻尼矩陣為Rayleigh阻尼矩陣,或滿足?1?1C[16]時,模態阻尼矩陣為對角矩陣,此時方程(11)為解耦的個獨立線性方程,其求解過程為標準的模態疊加法。當模態阻尼矩陣為非對角矩陣時,可以將模態阻尼矩陣分解成:

其中:d為模態阻尼矩陣中對應的對角元素組成的維方陣,d=diag[11,22,…,d];o為模態阻尼矩陣中對應的非對角元素組成的維方陣。

近似解耦方法就是忽略模態阻尼矩陣中的非對角線元素,即在方程(11)中用d代替,設運動方程的近似解為d,有:

3 矩陣對角占優指數

一般認為當模態阻尼矩陣的非對角元素相對于對角元素較小時,采用近似方法計算的動力響應誤差較小[17]。當模態阻尼矩陣滿足如下關系時,認為是對角占優矩陣[18]:

式(15)從數值上定義了矩陣的對角占優情況,但不能反映對角占優矩陣的占優程度,Meyer[19]基于對角占優矩陣的定義,推導了廣義矩陣對角占優指數計算公式:

其中:|o|為o各元素絕對值組成的矩陣;(|d?1||o|)為|d?1||o|的譜半徑。

當模態阻尼矩陣是對角矩陣時,0;而當模態阻尼矩陣為對角占優矩陣時,0<<1。越小,模態阻尼矩陣對角占優越顯著,一般認為近似解耦計算造成的誤差越小。

4 基于位移響應的誤差指數

結構位移響應在結構設計和結構分析中起著重要的作用?,F行規范仍舊以結構的位移響應作為評價指標,因此采用位移響應誤差來評價近似解耦計算的優劣是十分合適的。

定義近似解耦計算造成的位移誤差為

將方程(11)減去方程(14),有誤差響應方程為

Ma等[20?22]對有阻尼結構體系的振動問題進行了復模態分析,認為一般黏性阻尼系統的特征值和特征向量是共軛成對存在,在特定頻率諧振荷載下有阻尼結構體系能夠進行該頻率的模態阻尼振動,可以得到有阻尼結構體系的傳遞函數。

設模態坐標變換后的方程(11)受到頻率為的諧振激勵,即()=eiωt,為力的空間分布,本文取=[1,1,…,1]T;為幅值。

可以得到運動方程式(11)的穩態位移響應為

(i)為式(11)的頻響函數:

采用近似解耦法的運動方程式(14)的穩態位移響應為

d(i)為式(13)的頻響函數:

根據式(20)~(23),位移誤差響應式(19)的穩態響應可表示為

位移誤差響應式(24)體現了在諧振荷載作用下,結構采用近似解耦計算方法計算的位移響應誤差頻譜特性。由于忽略了模態阻尼矩陣中的非對角線元素,使得結構中各自由度的位移響應都與原系統響應產生偏差,因此采用2階范數能反映各自由度位移響應在結構整體響應的偏離程度。定義位移響應誤差指數:

當模態阻尼矩陣是對角矩陣時,(i)=0;當模態阻尼矩陣是非對角矩陣時,越小,結構采用近似解耦計算求得的結果誤差越小。

式(25)是根據單頻諧振荷載計算得到的位移響應誤差指數。而結構工程抗震計算中,地震波頻譜特性復雜,因此在評定采用近似解耦方法計算結構動力位移響應誤差程度時,可以選取不同頻率諧振荷載下計算得到的最大值作為該結構的位移誤差指數。

下面給出一個算例。為簡便起見,考慮2個二自由度結構體系。

模型一:

模型二:

考慮到地面運動的復雜性,本文選取9條地震波分別對兩模型進行近似解耦動力響應分析,研究其位移響應誤差情況。其中8條為天然地震波,1條為人工模擬地震波,各地震波水平分量的反應譜曲線如圖2所示。

1—Northridge;2—San Fernando;3—Turkiye;4—Taft;5—Tohoku;6—Whittie Narrows;7—EI Centro;8—Kobo;9—SHW2

經計算,模型一的位移誤差指數=0.012 3;模型二的位移誤差指數=0.029 3。圖3和圖4所示分別為兩模型不同地震作用下采用近似解耦計算得到的各自由度最大位移響應的相對誤差分布。模型一兩自由度最大位移響應的相對誤差平均值分別為0.64%和0.41%;而模型二兩自由度最大位移響應的相對誤差平均值分別為4.0%和2.5%,其相對誤差明顯比模型一的大。采用位移響應誤差指數能較準確的評定采用近似解耦計算對結構體系的影響。

圖3 模型一相對位移誤差分布

圖4 模型二相對位移誤差分布

5 豎向混合結構近似解耦計算

將豎向混合結構等效成兩自由度結構模型,每一自由度的動力特性由對應子結構的主頻確定。對該兩自由度模型進行不同特性下近似解耦計算的模態阻尼矩陣對角占優程度和位移響應誤差分析。

定義豎向混合結構的頻率比R和質量比R

式中:ωω分別為上、下部子結構的頻率,可由式(3)計算得到;MM分別為上、下部子結構的質量。

首先對豎向混合結構在不同子結構參數取值下的模態阻尼矩陣對角占優程度情況進行分析。不失一般性,取M=2 000 kg,ω=10 rad/s2,ξ=0.05,ξ=0.02。取R變化范圍為0到3,R變化范圍為0到1。

圖5所示為子結構質量比(R=0.1,R=0.5和R=8)一定時,模態阻尼矩陣對角占優指數隨頻率比的變化情況。由圖5可見:在不同質量比下,模態阻尼矩陣對角占優指數隨頻率比的變化相似。對角占優指數隨頻率比增大呈現先增大后減小的變化情況,在頻率比約為1.0時達到最大值。圖6所示為子結構頻率比(R=0.3,R=1.5和R=2.4)一定時,模態阻尼矩陣對角占優指數隨質量比的變化情況。由圖6可見:在較小質量比(小于0.2)時,質量比對角占優指數影響較明顯,之后對角占優指數隨質量比變化影響較小。

Rm:1—0.1;2—0.5;3—0.8

Rω:1—0.3;2—1.5;3—2.4

圖7所示為不同子結構特性下,模態阻尼矩陣對角占優指數隨質量比和頻率比變化的等高線分布圖。質量比對模態阻尼矩陣對角占優指數影響較小,在圖7中出現水平等高線。而模態阻尼矩陣對角占優指數隨頻率比增大呈現先增大后減小的變化,正如圖5所示那樣。當頻率比小于0.6時,對角占優指數小于0.2。頻率比在1.0~2.0范圍內,對角占優指數處于較大值,最大值出現在頻率比為1.0附近。

胡蘿卜:每100克含胡蘿卜素1.35~17.25毫克,還含有維生素B族、維生素C、脂肪及糖類和鐵、果膠、無機鹽等。

圖7 模態阻尼矩陣對角占優指數等高線

雖然模態阻尼矩陣對角占優指數在一定程度上反映了模態阻尼矩陣非對角線元素在模態阻尼矩陣中的大小程度,從感官上判斷忽略非對角線元素造成的誤差程度,但無法判定當模態阻尼矩陣對角占優指數為何值時,采用這種近似解耦方法所造成的誤差較小。

因此分別對豎向混合結構在不同子結構特性下忽略非對角線元素近似解耦計算的位移響應誤差進行分析。由于地震波頻率分布廣泛,在進行位移響應誤差分析時,將諧振頻率從0.1 rad/s2增加到100 rad/ s2,取其最大值為位移響應誤差指數。

圖8所示為子結構質量比(R=0.1,R=0.5和R=0.8)一定時,近似解耦計算位移誤差指數隨頻率比的變化情況??梢妶D8與圖5的變化規律相似:呈現先增大后減小的變化,在頻率比約為1.0時達到最大值。在頻率比一定時,位移誤差指數隨質量比的變化各異,但當質量比大于0.2時,整體上質量比對位移誤差指數的影響較小,如圖9所示。

Rm:1—0.1;2—0.5;3—0.8

Rω:1—0.3;2—1.5;3—2.4

圖10所示為不同子結構特性下,位移誤差指數隨質量比和頻率比變化的等高線分布??梢妶D10與圖7分布規律相似。質量比對位移誤差指數影響比較小,在圖10中近似呈現水平等高線。當頻率比小于0.6時,位移誤差指數小于0.02,可以認為采用近似解耦計算產生的位移誤差較小。而頻率比在0.7~2.0范圍內位移誤差指數大于0.02,最大值在頻率約為1.0處。而當頻率比大于2.0時,位移誤差指數又隨頻率比的增加緩慢減小。

圖10 位移誤差指數變化等高線

由于采用等效的兩自由度結構位移誤差指數來判定近似解耦計算的可行性,為了保證復雜結構的動力響應誤差在較小范圍,本文推薦位移誤差指數小于0.02時,可以采用近似解耦計算。對比圖7和圖10,模態阻尼矩陣對角占優指數和位移誤差指數在不同子結構特性下的衰退特性具有相似規律。而采用式(16)計算對角占優指數比采用式(25)計算位移誤差指數容易得多。經計算,第4節中的2個算例模型計算的對角占優指數分別為:模型一為0.133;模型二為0.350。因此在豎向混合結構分析中,可以采用模態阻尼矩陣對角占優指數來判定近似解耦計算的可行性。根據本文推薦的位移誤差指數取值,可以認為當模態阻尼矩陣對角占優指數小于0.2時,該結構體系可以采用近似解耦計算方法求解其動力響應。

6 結論

1)根據有阻尼結構體系的頻響函數得到的近似解耦計算的位移響應誤差指數能反映各自由度位移響應對結構整體響應的偏離程度。

2)在不同子結構特性下,模態阻尼矩陣對角占優指數和位移誤差指數的分布情況相似。模態阻尼矩陣對角占優指數和位移誤差指數隨質量比變化影響較小,而隨頻率比呈現先增大后降低的變化趨勢,在頻率比約為1.0時達到最大值。

3)在豎向混合結構分析中,可以采用模態阻尼矩陣對角占優指數來判定近似解耦計算的可行性。當模態阻尼矩陣對角占優指數小于0.2時,采用近似解耦計算方法求解豎向混合結構體系的動力響應誤差滿足一般工程分析的要求。

[1] 呂西林. 復雜高層建筑結構抗震理論與應用[M]. 北京: 科學出版社, 2007: 50?53. Lü Xilin. Seismic theory and application of complex high-rise structures[M]. Beijing: Science Press, 2007: 50?53.

[2] Foss K A. Co-ordinates which uncouple the equations of motion of damped linear dynamic systems[J]. ASME Journal of Applied Dynamics, 1958, 25(9): 361?364.

[3] Ma F, Imam A, Morzfeld M. The decoupling of damped linear systems in oscillatory free vibration[J]. Journal of Sound and Vibration 2009, 324(1/2):408?428.

[4] 陳國平. 粘性阻尼結構振動系統的實空間解耦和迭代求解[J].振動工程學報, 2000, 13(4): 559?566. CHEN Guoping. Real decoupled method and iterative solution of vibration system with viscous damping[J]. Journal of Vibration Engineering, 2000, 13(4): 559?566.

[5] Roesset J M, Whitman R V, Dobry R. Modal Analysis for Structures with Foundation Interaction[J]. Journal of the Structural Division, 1973, 99(3): 399?416.

[6] Hwang J S, Chang K C, Tsai M H. Composite damping ratio of seismically isolated regular bridges[J]. Engineering Structures, 1997, 19(1): 55?62.

[7] 呂西林, 張杰. 鋼和混凝土豎向混合結構阻尼特性研究[J]. 土木工程學報, 2012, 45(3): 10?16. Lü Xilin, ZHANG Jie. Damping behavior of vertical structures with upper steel and lower concrete components[J]. China Civil Engineering Journal, 2012, 45(3): 10?16.

[8] Rayleigh J W. The theory of sound: Vol.1[M]. 2nd ed. New York: Dover, 1945: 193?213.

[9] Knowles J K. On the approximation of damped linear dynamical systems[J]. Structural control and health monitoring, 2006, 13(1): 324?335.

[10] Prandinaa M, Mottersheada J E, Bonisolib E. An assessment of damping identification methods[J]. Journal of Sound and Vibration, 2009, 323(3/4/5): 662?676.

[11] Hwang J H, Ma F. On the approximate solution of non-classically damped linear systems[J]. ASME Journal of Applied dynamics, 1993, 60(9): 695?701.

[12] Morzfeld M, Ajavakom N, Ma F. Diagonal dominance of damping and decoupling approximation in linear vibratory systems[J]. Journal of Sound and Vibration, 2009, 320 (18): 406?420.

[13] 桂國慶, 何玉敖. 非比例阻尼結構體系近似解耦分析中的誤差研究[J]. 工程力學, 1994, 11(4): 40?50. GUI Guoqing, HE Yu’ao. Study on errors of approximate decoupling analysis for non-proportionally damped structural systems[J]. Engineering Mechanics, 1994, 11(4): 40?45.

[14] Clough R, Penzien J. Dynamics of structures[M]. Berkeley, California, USA: Computers and Structures Inc, 2003: 234?245.

[15] Chopra A K. Dynamics of structures: theory and applications to earthquake engineering[M]. 3rd ed. New Jersey, USA: Prentice Hall International, 2007: 425–429.

[16] Udwadia F E. A Note on nonproportional damping[J]. Journal of Engineering Mechanics, 2009, 135(11): 1248?1256.

[17] Morzfeld M, Ma F, Ajavakom N. On the decoupling approximation in damped linear systems[J]. Journal of Vibration and Control, 2008, 14(12): 1869?1884.

[18] Horn R A, Johnsom C R. Matrix analysis[M]. Oxford, UK: Cambridge University Press, 1985: 402?410.

[19] Meyer C D. Meyer matrix analysis and applied linear algebra[M]. Philade Phia, USA: Society for Industrial and Applied Mathematics, 2001: 660?667.

[20] Ma F, Morzfeld M, Imam A. The decoupling of damped linear systems in free or forced vibration[J]. Journal of Sound and Vibration, 2010, 329(15): 3182?3202.

[21] Morzfeld M, Ma F. The decoupling of damped linear systems in configuration and state spaces[J]. Journal of Sound and Vibration, 2011, 330(15):155?161.

[22] 曹樹謙, 張文德, 蕭龍翔. 振動結構模態分析: 理論、實驗與應用[M]. 天津:天津大學出版社, 2002: 17?37. CAO Shuqian, ZHANG Wende, XIAO Longxiang. Vibration modal analysis: Theory, experiment and application[M]. Tianjin: Tianjin University Press, 2002: 17?37.

(編輯 楊幼平)

Errors of approximate decoupling analysis for vertically mixed structures

HUANG Wei, QIAN Jiang, LIANG Feifei, ZHOU Zhi

(State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China)

Through equivalent two degree-of-freedom (2-DOF) models with different sub-structural dynamic properties, the quantifications of diagonal dominance and decoupling errors were calculated to verify the degree of approximation and feasibility of this method. The results show that with different sub-structural dynamic properties, the distribution of diagonal dominance and displacement errors of the models are similar. Therefore, the diagonally dominant index can be adopted to analyze the feasibility of decoupling approximation method in calculation of dynamic response of vertical mixed structure.

vertically mixed structure; non-proportional damping matrix; modal damping matrix; diagonal dominance; decoupling errors

10.11817/j.issn.1672-7207.2015.04.036

TU398

A

1672?7207(2015)04?1454?07

2014?04?13;

2014?06?15

國家“十二五”科技支撐計劃項目(2012BAJ13B02);國家自然科學基金資助項目(91315301-4)(Project (2012BAJ13B02) supported by the National Science and Technology Pillar Program During the 12th Five-Year Plan Period; Project (91315301-4) supported by the National Natural Science Foundation of China)

黃維,博士研究生,從事結構抗震及數值計算研究;E-mail:2008huangwei@tongji.edu.cn

猜你喜歡
模態結構
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結構
車輛CAE分析中自由模態和約束模態的應用與對比
國內多模態教學研究回顧與展望
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 亚洲人成电影在线播放| av午夜福利一片免费看| 成年片色大黄全免费网站久久| 久996视频精品免费观看| 国产精品妖精视频| 亚洲国产精品VA在线看黑人| 亚洲一区色| 9999在线视频| 亚洲第一页在线观看| 99久久无色码中文字幕| 999精品在线视频| 国产杨幂丝袜av在线播放| 亚洲日本www| 国产乱视频网站| 色欲色欲久久综合网| 国产不卡一级毛片视频| 十八禁美女裸体网站| 国产乱肥老妇精品视频| 国产成人毛片| 免费不卡视频| 人妻无码一区二区视频| 国产精品美女在线| 久久婷婷国产综合尤物精品| 亚洲国产精品久久久久秋霞影院 | 亚洲精品欧美日韩在线| 在线视频一区二区三区不卡| 久久香蕉国产线| 国产成人精品高清不卡在线 | 国产无人区一区二区三区 | 午夜视频日本| 国产成人精品一区二区不卡| 日韩成人高清无码| 成人在线不卡视频| 欧美激情第一区| 成人午夜精品一级毛片 | 园内精品自拍视频在线播放| 欧美色视频在线| 久久久久无码国产精品不卡| 国产亚洲精品yxsp| 无码免费的亚洲视频| 操国产美女| 亚洲男人在线天堂| 无码高潮喷水专区久久| 亚洲天堂伊人| 亚洲开心婷婷中文字幕| 天天躁夜夜躁狠狠躁图片| 亚洲精品不卡午夜精品| 99er精品视频| 国产农村精品一级毛片视频| 亚洲精品制服丝袜二区| 潮喷在线无码白浆| 成人国产小视频| 亚洲第一天堂无码专区| 国产主播一区二区三区| 乱人伦中文视频在线观看免费| 久久国产精品嫖妓| 最新国语自产精品视频在| 高潮毛片无遮挡高清视频播放| 一级成人a做片免费| 国产99在线| 亚洲AV成人一区二区三区AV| 熟女视频91| 国产三级a| 亚洲va在线观看| 中国国产高清免费AV片| 婷婷色中文网| 欧美色视频网站| 国产网站一区二区三区| 嫩草在线视频| 亚洲欧美激情另类| 中文字幕色在线| 久久精品中文无码资源站| 国产高清无码麻豆精品| 亚洲日本韩在线观看| 免费一极毛片| 亚洲av综合网| a在线亚洲男人的天堂试看| 日韩精品无码免费专网站| 国产精品短篇二区| 少妇露出福利视频| 免费人欧美成又黄又爽的视频| 欧美a网站|