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

基于IDA的大跨網架結構豎向地震動力性能研究

2017-01-10 07:27:10張輝東戚國營
天津城建大學學報 2016年6期
關鍵詞:結構分析

張輝東,戚國營,閆 立

(天津城建大學 a. 天津市土木建筑結構防護與加固重點實驗室;b. 土木工程學院,天津 300384)

土木工程

基于IDA的大跨網架結構豎向地震動力性能研究

張輝東a,b,戚國營a,b,閆 立a,b

(天津城建大學 a. 天津市土木建筑結構防護與加固重點實驗室;b. 土木工程學院,天津 300384)

大跨空間網架結構具有跨度大、豎向剛度低、阻尼小、對豎向地震作用比較敏感等特點,因此,在豎向地震作用下,對該類結構體系非線性動力性能的研究尤為重要.基于有限元法,針對我國某地區一實際大跨網架結構,對其進行豎向地震激勵下的動力性能研究,采用增量動力分析(IDA)方法對該結構進行了抗震性能評價,得到了該結構的薄弱位置及倒塌極限.研究表明:該網架結構具有良好的豎向抗震性能;隨著結構塑性發展,IDA曲線的離散性明顯;結構概率分位50%,對應的倒塌極限位移為0.337,m;薄弱構件主要分布于下弦,尤其是下弦節點附加質量區域,且桿件的塑性發展集中于附加質量區域.

大跨網架結構;抗震性能;增量動力分析;倒塌

大跨空間結構以其優點發展迅速,廣泛用于車站、機場、會堂、體育館等大面積屋蓋公共建筑中,其中應用最廣泛的當屬網架結構.大跨網架結構跨度大、阻尼小、平面外剛度低,對豎向地震響應較敏感,其產生的局部或整體倒塌對人民生命財產造成的損失極大.因此,對豎向地震激勵下該類結構薄弱部位及倒塌問題的研究顯得尤為重要.

針對網架結構的豎向抗震分析,張毅剛等[1]利用子空間迭代法研究了網架結構的自震特性,利用反應譜法和時程分析法研究了豎向地震作用下的內力分布規律;丁萬尊等[2]在考慮阻尼的情況下,利用擬夾層板法分析了豎向地震激勵下網架結構的動力問題;沈祖炎等[3]認為網架結構的豎向振型為主要振型,豎向激勵對加速度響應具有放大作用;三向地震作用下,豎向地震處于主導地位;孫夢涵等[4]認為不可忽視豎向地震對網架變形的影響,是否考慮下部支承剛度對結構的豎向地震響應影響較小;黃興淮等[5]提出網架結構在豎向強震作用下的倒塌破壞機理,認為弦桿首先發生屈曲,進行內力重分布,某一關鍵桿件的失效會產生多米諾骨牌效應,最終導致結構整體倒塌;韓慶華等[6]對大跨格構式拱結構進行增量動力分析(IDA),分析結果與實驗數據一致,驗證了IDA方法在空間結構體系性能評價方面的準確性;孟凡林等[7]利用IDA方法對大跨網架結構進行了單條地震動響應研究,獲得結構的倒塌極限及薄弱區域,為網架結構的多條地震動IDA提供了一種參考.

增量動力分析(incremental dynamic analysis,簡稱IDA)方法,是一種建立在非線性動力時程分析基礎之上的動力參數分析理論.這種方法最早由Bertero[8]提出,之后由Vamvatiskos和Cornell[9]進行了系統性研究和改進,并被FEMA350[10]和FEMA351[11]標準采用,現已廣泛用于基于性能的結構抗震評價.IDA理論的基本原則是調整地震動記錄的強度幅值IM,以此作為輸入激勵進行非線性動力時程分析,獲得不同IM值下的結構響應值DM.IDA涵蓋了結構在地震激勵下的整個變化過程,對結構從彈性到屈服、從屈服到塑性發展,直至結構倒塌作出全面評估,反映結構在不同強度地震作用下的承載力、剛度及位移變化,以此評價結構的整體抗震性能.

近年來,IDA理論作為一種強有力的非線性分析方法[12]被國內廣泛接受,成為解決工程問題的有力手段,其分析數據可以精確地反映出結構隨地震動強度的改變而發生的變化[13],能較好地評價結構的抗震性能.目前,基于IDA方法的研究主要集中于框架、框剪等結構形式,對于網架等空間結構在豎向地震作用下的性能分析較少.本文以某實際網架結構為例進行增量動力分析,研究該結構在豎向地震激勵下的剛度變化、最大節點位移、薄弱部位以及倒塌問題,評估該結構的抗震性能.

1 模型方程

1.1 結構的運動方程

根據動力學相關理論[14],該網架結構的動力方程表示為

式中:[M]為質量矩陣;[C]為阻尼矩陣;[K]為結構剛度矩陣;{üg}為地震激勵加速度;{u}為節點位移向量;{˙}為節點速度向量;{ü}為節點加速度向量.

其中,結構剛度[K]為常量,當結構進入塑性后,其整體剛度與運動的時間歷程{u(t)}及材料的非線性相關,隨著塑性的發展程度不斷變化.考慮地震動的隨機性和彈塑性響應中結構剛度的變化,將式(1)表示為

式(2)中,[K(t)]表示t時刻結構的剛度矩陣,其他參數意義如式(1).

1.2 阻尼模型

阻尼屬于結構的基本動力屬性,用來描述結構振動過程中的綜合耗能能力,是影響結構動力問題的關鍵因子之一.因此,阻尼模型的合理選取對于結構的動力分析十分關鍵.非線性增量動力分析中,常采用瑞利(Rayleigh)阻尼,考慮到高階振型對阻尼的影響,曾提出多種修正的Rayleigh阻尼模型,如張輝東等[15]提出了一種適用于網架結構的節點-構件阻尼模型.考慮到本文的研究重點和分析精度要求,選擇瑞利阻尼模型[14]

式中:α、β分別表示質量比例阻尼和剛度比例阻尼,均為未知系數.根據振型正交性,α、β與振型阻尼比之間的關系為

由式(4)可以看出:任意給定兩個阻尼比,就可確定α、β系數.阻尼比取值直接影響著彈塑性分析中的地震響應,GB50011—2010《建筑抗震設計規范》[16]規定:該類型網架的阻尼比取值0.02.曹資等[17]認為空間結構的阻尼比隨著振型的增加而增大,隨結構的塑性發展而增大,最大值可達0.05;此外,彈塑性分析中考慮材料非線性相當于考慮了塑性阻尼,同時考慮材料非線性和塑性阻尼會放大結構的阻尼作用[18].因此,在滿足分析要求的前提下,一般指定結構的阻尼比相等,本文阻尼比取值為0.03.

另外,Rayleigh阻尼與結構的振型相關,選擇參考振型應考慮對結構動力響應有重要影響的頻率分量,一般情況下,取前2階自振頻率確定阻尼系數[19-20].

2 IDA方法

IDA方法分為單記錄IDA法和多記錄IDA法:前者屬于確定性分析,而地震動的隨機性和不確定性決定了該法并不能準確評估結構在未知地震作用下的破壞特性和趨勢;后者應選擇足夠多的地震動記錄,且記錄包含結構可能遭受的最強地震動部分,符合此類條件的多條IDA曲線才能全面、準確地評估結構的抗震性能.

IDA方法是指對已選用的某條地震動記錄根據需要進行調幅,得到單調增大的一組地震動記錄,以此作為輸入激勵,對結構進行非線性動力時程分析,得到結構響應參數,直到結構倒塌失穩為止.根據一系列IM值和DM值,繪制對應于此條地震動的IDA曲線.重復以上步驟,即得到多條IDA曲線.關鍵步驟如下.

(1)地震動記錄選取.時程分析中,地震動記錄的選擇對結構的性能分析十分關鍵.根據場地類別和地震分組等條件,選擇頻譜特性相接近的地震動記錄.

為滿足IDA方法對地震強度變化的要求,需調整地震動記錄.初始地震波強度峰值的調整公式為

式中:ai(t)、Aimax為第i次調整后的地震波曲線及強度峰值;a(t)、Amax為初始地震波曲線及強度峰值;φi為調幅系數.Vamvatiskos和Cornell[9]提出一種“Hunt&fill”原則確定φi;呂大剛等[21]曾提出一種“折半取中”的原則確定φi.本文采用等步調幅.

地震動持時是影響時程分析的重要因素.持時不同,輸入能量不同,結構損傷不同,結構響應有所區別.地震動持時應包含地震動最強記錄時段;文獻[16]規定:地震波的持續時間不宜小于結構自振周期的5倍和15,s.本文輸入地震波持時20,s.

(2)IM、DM參數的選擇.對于空間結構體系,最常采用峰值加速度(PGA)作為強度參數IM[6].結構的抗震性能水平和抗破壞能力可由結構某參數或自定義指標來確定,例如層間位移角、最大節點位移等.

(3)極限狀態點的確定準則.強震作用下,桿件塑性發展,結構剛度退化,結構響應發散,即結構失效倒塌;此外,可從最大位移角度評價結構極限狀態;FEMA規定[11]:結構的倒塌破壞準則,即單條IDA曲線調幅前后兩點的連線斜率小于0.2倍的彈性斜率時,則取調幅前的點作為極限點.本文極限狀態點的確定參考此類準則.

3 網架實例

3.1 網架模型

采用河北某實際網架作分析算例:結構7度(0.10g)設防,Ⅱ類場地,地震分組第2組;正放四角錐雙層網架,縱向長度136,m,橫向長度120,m,網架高度2.5,m,網格大小4,m×4,m;高強螺栓節點連接;采用下弦多點支撐,下部支撐采用箱型和H型截面鋼柱,柱距分別為8,24,28,m,柱網布置見圖1.鋼材采用Q345鋼,桿件截面及荷載分布見表1,其中附加質量以噸(t)計.3.2 荷載

圖1 柱網布置及節點附加質量分布

表1 桿件規格及荷載分布

對于網架結構,荷載均作用于節點,桿件不承受橫向荷載,將上弦均布荷載等效為節點處的集中力,表示為

式中:qi表示均布荷載;A表示節點所承擔荷載的面積.對表1均布荷載,據式(6)等效為恒載FG=AqG=16×0.40=6.4,kN,活載FQ=AqQ=16× 0.50=8.0,kN.

根據桿件截面參數和材料參數,程序自動計算桿件的自重,網架節點的重量一般取桿件自重的30%,.

3.3 數值模型

研究該網架結構的豎向地震響應,可將下部支撐簡化為三向鉸支[4],對本文研究內容的影響在可控范圍內.有限元分析模型中,鋼材采用雙線性隨動硬化模型;采用桿單元模擬,每個桿件離散為一個分析單元,桿件兩端彎矩釋放.有限元模型見圖2.考慮非線性重力荷載的初始條件,采用非線性直接積分的方法進行彈塑性分析,積分方法采用Newmark-β法.該結構阻尼采用瑞利阻尼,阻尼比取為0.03.分析過程需同時考慮結構的幾何非線性和以塑性鉸為表征的材料非線性.桿件中部設置延性軸力鉸,鉸長度取值

圖2 有限元模型

式中:LP為鉸長度,塑性變形長度范圍;L0為桿件計算長度.塑性鉸代表桿件出現塑性變形(桿件易損),易損準則為

式中:A表示毛截面面積;σcr表示屈服極限應力值,與長細比λ有關;Pcr表示極限軸力.

4 結果與分析

綜合考慮地震動震級、場地土和結構周期等因素,從PEER強地震動數據庫[22]中選取9條地震動記錄作為輸入樣本,見表2.文獻[16]規定,7度設防的時程分析加速度取值:小震35,cm/s2,中震105,cm/s2,大震220,cm/s2.通過調幅分別取PGA為0.035g、0.105g、0.22g、0.3g、0.4g、0.5g、0.6g、0.7g、0.8g等直至結構倒塌.

表2 地震動參數

4.1 結構的整體響應分析

采用PGA作為輸入參數IM,豎向最大節點位移作為DM,通過數值分析繪制結構在不同輸入樣本下的IDA曲線,如圖3所示.

圖3 豎向最大節點位移

由圖3可知:不同地震動對應的IDA曲線差異性較大,尤其隨著地震動強度的增加,離散性愈明顯,這主要由地震動的隨機性所決定.但從IDA曲線簇來看,該結構在不同地震作用下的響應趨勢是一致的,大體經歷彈性、彈塑性、塑性三個階段.各個IDA曲線初始部分斜率較大,均存在一個線彈性階段[23],因為小震作用下,桿件沒有塑性發展,結構處于彈性階段.隨著地震動強度的增大,部分桿件由彈性狀態進入彈塑性狀態,結構剛度降低,表現為曲線斜率不同程度的減小.隨著輸入強度的繼續增大,塑性發展加劇,結構的破壞繼續發展,IDA曲線斜率急劇降低,即在很小IM增幅的情況下,損傷指標DM產生不成比例的增加,明顯存在斜率突變點(剛度退化點),即預示著結構倒塌破壞.

豎向地震作用下,結構在PGA為0~0.22g范圍內未出現剛度退化極限,說明7度大震作用下,該網架結構沒有倒塌破壞,其設計是安全的.

在上述分析數據基礎上,利用統計方法得到豎向地震作用下,結構概率分位值10%,、50%,、90%,對應的倒塌極限位移分別為0.234,0.337,0.440,m,見圖4.50%,分位代表該結構在可能遭遇地震下的平均響應水平,即結構最大位移達到0.337,m時即可視為倒塌破壞,10%,和90%,分位代表結構極限位移的離散程度.4.2 基于單條IDA的響應分析

圖4 倒塌極限位移值

以Kobe波為例(見圖5),介紹該網架的豎向地震響應:最大節點位移、桿件軸力、易損桿件分布及結構薄弱區域等.

圖5 Kobe地震動記錄

4.2.1 基于Kobe波的響應分析

大震PGA=220,cm/s2作用下,位移最大的節點(505)的加速度和位移時程曲線見圖6-7.取代表性桿件#3127、#3159、#3158為研究對象,得到其軸力時程曲線見圖8-10,節點及桿件定位見圖11標示.

圖6 505節點加速度時程曲線

圖7 505節點z向位移時程曲線

由圖5-7可以看出:最大節點位移時程曲線與加速度時程曲線基本一致,說明在大震作用下,該結構整體始終處于彈性階段;位移峰值的響應稍滯后于加速度峰值,這是結構振動過程中由阻尼作用產生的合理響應.

圖8 #3127軸力時程曲線

圖9 #3159軸力時程曲線

圖10 #3158軸力時程曲線

分析圖8-10知:在大震激勵過程中,#3127雖有鉸出現,但桿件塑性變形極小,可近似認為處于彈性階段,其軸力時程曲線與地震波曲線(見圖5)基本一致;#3159在地震激勵6.5,s左右軸力突變,桿件發生塑性變形,承載能力急劇減小,之后隨地震動強度的變化重新回到彈性工作狀態;#3158在地震作用6.7,s左右于A點出現軸力極值,桿件在A點發生塑性破壞,軸力不再隨地震激勵的改變而改變,經歷3.5,s之后于點B處變為0,說明該桿件斷裂,退出工作.

4.2.2 薄弱構件分布

程序中材料的非線性由塑性鉸體現,鉸出現的先后次序及分布表明了結構中的易損桿件及薄弱區域分布.地震激勵PGA=220,cm/s2作用下,結構下弦于2.47,s首現失效桿件(鉸),自2.47~5.00,s時段,地震強度較小,結構桿件均處于彈性狀態(無鉸出現);從5~8,s時段內,包含地震動最強部分(見圖5),且前期輸入能量不斷累積,下弦多個桿件發生塑性變形.當輸入樣本PGA=500,cm/s2,時,結構于12.48,s倒塌破壞.圖11-12為PGA=220,500,cm/s2(第12,s)作用下結構薄弱桿件的分布.由圖11-12可大致看出該結構的易損桿件、薄弱區域的分布,為該類型結構的設計和加固提供一種參考.豎向地震作用下,該網架薄弱桿件大部分位于上、下弦桿,腹桿極少存在塑性變形;與上弦桿相比,下弦失效桿件數量較多,且多集中于附加質量節點區域;節點附加質量越大,產生地震力越大,桿件發生塑性變形的可能性相應增加.此外,薄弱桿件的分布比較集中,結構中某一單元的失效會導致連鎖反應,初始破壞單元的失效,引起與之相互作用的單元發生連鎖失效[5],最終導致結構的整體失穩.因此,網架結構設計時,要特別注意附加質量節點區域桿件的設計.

圖11 PGA為220,cm/s2,上、下弦薄弱桿件分布

圖12 PGA為500,cm/s2,(第12,s)上、下弦薄弱桿件分布

5 結 論

(1)基于IDA方法的動力性能研究可以較好地評估網架結構的抗震性能,為豎向地震作用下網架結構的整體倒塌破壞研究提供一種方法.

(2)該網架結構具有良好的抗震性能,在設防烈度7度的不同地震動大震激勵下,雖有不同程度、不同數量塑性鉸的出現,但均沒有整體倒塌,該結構設計是合理、經濟和安全的.

(3)隨著塑性發展,結構的非線性得以體現,IDA曲線的離散性明顯.結構概率分位50%,對應的倒塌極限位移為0.337,m,即最大位移達到0.337,m時即可視為結構倒塌破壞.

(4)豎向地震作用下,網架結構的弦桿為薄弱桿件,且下弦桿塑性鉸數目較上弦桿多;桿件的塑性發展集中于附加質量區域.在此類型網架設計時,要注意附加質量區域桿件的設計,以防該區域桿件的塑性發展導致局部失穩及整體倒塌的危險.

[1] 張毅剛,藍惆恩. 網架結構在豎向地震作用下分析方法[J]. 建筑結構學報,1985,6(5):2-14.

[2] 丁萬尊,陳慶云. 正交正放類網架結構動力響應的擬夾層板分析法[J]. 建筑結構學報,1997,18(1):18-25.

[3] 沈祖炎,陳學潮,陳揚冀. 網架結構模型抗震實驗研究[J]. 土木工程學報,1994,27(1):30-40.

[4] 孫夢涵,范 峰,支旭東,等. 考慮下部支承四角錐網架結構動力響應分析[J]. 土木工程學報,2014:47(12):10-16.

[5] 黃興淮,徐趙東,楊明飛. 多維地震下大跨網格結構倒塌分析與抗倒塌分析[J]. 東南大學學報,2012,42(1):109-113.

[6] 韓慶華,盧 燕,徐 穎,等. 基于IDA的格構式結構倒塌能力研究性能分析[J]. 土木工程學報,2015,48(3):1-7.

[7] 孟凡林,孟祥瑞,徐德英,等. Sap2000在大跨空間結構抗震分析設計中的應用[J]. 吉林建筑工程學院學報,2011,28(4):1-4.

[8] BERTERO V V. Strength and deformation capacities of buildings under extreme environments[J]. Structural Engineering and Structure Mechanics,1977,53(1):29-79.

[9] VAMVATISKOS D,CORNELL C A. Incremental dynamic analysis[J]. Earthquake Engineering and Structure Dynamics,2002,31(3):491-514.

[10] FEMA. Recommended seismic design criteria for new steel moment-frame buildings:FEMA-350[R]. Washington,D.C.:Federal Emergency Management Agency. SAC Joint Venture,2000.

[11] FEMA. Recommended seismic evaluation and upgrade criteria for existing welded steel moment-frame buildings:FEMA-351[R]. Washington,D.C.:Federal Emergency Management Agency.SAC Joint Venture,2000.

[12] VAMVATISKOS D,CORNELL C A. Applied incremental dynamic analysis[J]. Earthquake Spectra,2004,20(2):503-522.

[13] VAMVATISKOS D,CORNELL C A. Direct estimation of seismic demand and capacity of multi-degree of freedom systems through incremental dynamic of single degree of freedom approximation incremental[J]. Journal of Structural Engineering,2005,131(4):589-599.

[14] CHOPRA A K,謝禮立,呂大剛,等. 結構動力學:理論及其在地震工程中的應用[M]. 北京:高等教育出版社,2005:276-277,339-341.

[15] 張輝東,王元豐. 空間網殼結構的節點-構件阻尼模型研究[J]. 土木工程學報,2015,48(2):55-61.

[16] 中華人民共和國住房和城鄉建設部.建筑抗震設計規范:GB50011—2010[S]. 北京:中國建筑工業出版社,2010:32-315.

[17] 曹 資,薛素鐸,王雪生,等. 空間結構抗震分析中的地震波選取和阻尼比取值[J]. 空間結構,2008,14(3):4-8.

[18] 楊志勇,黃吉峰,田家勇. 罕遇地震彈塑性靜、動力分析方法中結構阻尼問題的探討[J]. 地震工程與工程振動,2009,29(6):116-120.

[19] 何 利,葉獻國. 不同阻尼模型對框架結構地震反應影響因素[J]. 合肥工業大學學報(自然科學版),2011,34(7):1 028-1 030.

[20] 李 田. 結構時程動力分析中的阻尼取值研究[J]. 土木工程學報,1997,30(3):68-73.

[21] 呂大剛,于曉輝,王光遠. 基于單地震動記錄IDA方法的結構倒塌分析[J]. 地震工程與工程振動,2009,29(6):34-39.

[22] Pacific Earthquake Engineering Research Center. PEER ground motion database[DB/OL]. [2014-10-25]. http:// ngawest2. berkeley. edu/.

[23] 李 磊,鄭山鎖,李 謙. 基于IDA的型鋼混凝土框架的地震易損性分析[J]. 廣西大學學報,2011,36(4):536-541.

Analysis on the Vertical Seismic Performance of a Large Span Grid Structure Based on IDA Method

ZHANG Huidonga,b,QI Guoyinga,b,YAN Lia,b
(a. Tianjin Key Laboratory of Civil Buildings Protection and Reinforcement;b. School of Civil Engineering,TCU,Tianjin 300384,China)

The large span space grid structure is highly sensitive to the vertical component of the seismic action due to the large span,low vertical stiffness and its small damping property,therefore,a research of nonlinear dynamic behaviors on this structural system subjected to the vertical seismic component is particularly important.Based on the finite element method,the nonlinear dynamic behaviors of a real large span grid structure located in a region of China are investigated under earthquake actions.The seismic performance of the structure is evaluated using the Incremental Dynamic Analysis(IDA)method.The weak location and collapse ultimate point of the structure is obtained.Results show that the structure has a good seismic performance;the difference among the IDA curves is remarkable as plasticity development of members.It is observed that the collapse limit displacement of 50%,probability is 0.337 m,and that weak members occur at the lower region where the masses are assigned.

large span grid structure;seismic behavior;incremental dynamic analysis;collapse

TU393.3

A

2095-719X(2016)06-0416-07

2015-10-25;

2016-03-02

國家自然科學基金青年基金(51108301)

張輝東(1976—),男,山東青島人,天津城建大學教授,博士.

猜你喜歡
結構分析
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
隱蔽失效適航要求符合性驗證分析
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統及其自動化發展趨勢分析
論《日出》的結構
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
中西醫結合治療抑郁癥100例分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 欧美精品H在线播放| 奇米精品一区二区三区在线观看| 国产91精品调教在线播放| 国产精品免费入口视频| 亚洲中文制服丝袜欧美精品| 国产精品熟女亚洲AV麻豆| 国产精品免费p区| 91亚洲视频下载| 国产欧美日韩免费| 国产黑丝视频在线观看| 无码国产偷倩在线播放老年人| 日韩精品无码免费专网站| 拍国产真实乱人偷精品| 91精品国产一区自在线拍| 欧美日本激情| 99久久国产综合精品2020| 国产好痛疼轻点好爽的视频| 亚洲日韩精品综合在线一区二区| 亚洲一区二区约美女探花| 视频一本大道香蕉久在线播放 | 色九九视频| 国产一区二区精品福利| 欧美福利在线| 日韩在线中文| 四虎在线观看视频高清无码| 国产一级小视频| 精品少妇人妻无码久久| 亚洲资源站av无码网址| 成人综合在线观看| 国产熟睡乱子伦视频网站| 国产精品第三页在线看| 日本手机在线视频| 国产色爱av资源综合区| 四虎永久免费地址| 国产成人亚洲精品无码电影| 欧美成人一级| 久久综合一个色综合网| 伊人查蕉在线观看国产精品| 四虎永久免费网站| a天堂视频| 亚洲成aⅴ人在线观看| 午夜日韩久久影院| 91精品在线视频观看| 激情六月丁香婷婷| 亚洲视频在线青青| 日韩AV无码免费一二三区| 91在线激情在线观看| 亚洲第一成人在线| 久草性视频| a毛片基地免费大全| 777午夜精品电影免费看| 日韩人妻无码制服丝袜视频| 久草中文网| 人妻中文字幕无码久久一区| 一级毛片免费的| 久久免费成人| 免费一级全黄少妇性色生活片| 试看120秒男女啪啪免费| 精品91在线| 99精品视频九九精品| 亚洲色图欧美视频| 91精品国产自产在线观看| 久久亚洲黄色视频| 欧美a√在线| 人妻无码中文字幕一区二区三区| 青草午夜精品视频在线观看| 亚洲第一中文字幕| 亚洲成在线观看| 91精品国产自产在线老师啪l| 国产又粗又爽视频| 欧美在线导航| 动漫精品中文字幕无码| 一本大道香蕉中文日本不卡高清二区 | 亚洲午夜天堂| 国产性生交xxxxx免费| 亚洲美女一级毛片| 亚洲午夜天堂| a毛片在线播放| 无码'专区第一页| 乱人伦99久久| 欧美一区福利| 伊人福利视频|