亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        實用粘彈性阻尼器耗能結構非平穩(wěn)地震響應的快速求解

        2018-09-10 11:39:16李創(chuàng)第朱騰飛柏大煉葛新廣
        廣西科技大學學報 2018年4期

        李創(chuàng)第 朱騰飛 柏大煉 葛新廣

        摘 要:本文為提高實用粘彈性阻尼器耗能減震結構基于非平穩(wěn)地震激勵的響應分析效率,對實用粘彈性阻尼器單自由度耗能系統(tǒng)隨機地震響應的數(shù)值分析方法進行了系統(tǒng)研究.首先,采用設置支撐的六參數(shù)實用粘彈性阻尼器進行建模;然后,利用傳遞函數(shù)法直接得到耗能減震結構系統(tǒng)瞬態(tài)響應精確解;最后,基于虛擬激勵法獲得了快速求解調制非平穩(wěn)激勵下耗能減震結構的時域瞬態(tài)響應.通過均勻非平穩(wěn)和非均勻非平穩(wěn)算例分析表明:該方法對于此類問題的計算具有效率高、工程應用強的特點,為設置粘彈性耗能減震結構在非平穩(wěn)地震激勵下的快速求解提供了參考.

        關鍵詞:實用粘彈性阻尼器;耗能減震系統(tǒng);非平穩(wěn)地震響應;快速求解

        中圖分類號:TU311.3 DOI:10.16375/j.cnki.cn45-1395/t.2018.04.003

        0 引言

        工程上,為減小結構體系的地震動響應,常采用增加結構阻尼的方法,其中粘彈性阻尼裝置可以有效提供阻尼被廣泛采用[1-3].粘彈性阻尼器有著非常復雜的力學性能,且其性能易受外部激勵和環(huán)境因素的影響,因此,研究人員為了精確描述其本構關系,提出了多種恢復力模型,主要有:Kelvin模型[3]、Maxwell模型[3]、分數(shù)導數(shù)模型[3]、六參數(shù)模型[4]等.其中,Mazza F和Vulcano A提出了一種六參數(shù)實用粘彈性阻尼器模型.該實用粘彈性阻尼器是由兩支Maxwell單元并聯(lián)一支Kelvin單元組成,其本構方程從力學性能說是Kelvin模型[3]和Maxwell模型[3]的復合,能更準確的模擬粘彈性材料的蠕變及松弛性能,模型計算參數(shù)更易于測試數(shù)據(jù)擬合.此外,該模型在原結構上可進行擴階求解,因而已廣泛應用于實際工程當中.在實際工程中,阻尼器需要與支撐串聯(lián)安裝,常用的支撐形式有人字形支撐和對角式支撐,近年來還出現(xiàn)了剪刀式和肘桿式等新型支撐[5],支撐的水平等效剛度直接影響到阻尼器的減震效果.因此,采用設置支撐的六參數(shù)實用粘彈性阻尼模型分析了耗能結構的動態(tài)特性,具有良好的工程應用價值.

        一般粘彈性耗能減震變頻結構常用擴階精確法和非擴階近似法作為分析方法.擴階精確法針對廣義Maxwell[6-9]、GHM[10]、分數(shù)導數(shù)Kelvin[11]等易擴階粘彈性近似模型,利用擴階復模態(tài)法獲得結構響應解析解,因擴階方程組物理意義不明確、變量個數(shù)劇增、計算效率低等缺陷,難以用于實際工程的分析.非擴階近似法主要是模態(tài)應變能法[1,12]和取結構基頻的強行振型解耦法[13],但因其在多自由結構振動分析時忽略模態(tài)交叉項的影響,使其精度和使用范圍受到限制.

        虛擬激勵法[14]將平穩(wěn)振動分析轉化為簡諧振動分析和確定性時間歷程分析,在計算步驟簡化的基礎上仍保持理論上的高度精確性.虛擬激勵法與鐘萬勰[15]提出的精細積分法二者相結合的方法是運動方程的時域數(shù)值求解中常用的方法,簡化了計算過程的同時效率得到很大提升.Su等[16]針對均勻調制非平穩(wěn)激勵下非線性結構的響應分析,提出了一種快速等效線性化方法.傳統(tǒng)的運動方程在時域內(nèi)的求解存在矩陣的階數(shù)較高,不便于大型的時程響應分析.非齊次的運動方程均需計算指數(shù)矩陣以及矩陣求逆,而矩陣求逆容易造成計算精度降低、出現(xiàn)數(shù)值不穩(wěn)定以及需要考慮逆矩陣不存在的情況,從而限制了使用范圍.非平穩(wěn)隨機振動分為均勻非平穩(wěn)激勵和非均勻平穩(wěn)激勵,兩者都可以表示為同源平穩(wěn)激勵關于時間的函數(shù)[17-18],可以利用虛擬激勵法化為確定性的時間歷程分析.

        本文應用傳遞函數(shù)法,直接在原結構基礎上獲得設置支撐六參數(shù)實用粘彈性阻尼器耗能減震結構在任意激勵和非零初始條件下位移、速度和阻尼器受力時域瞬態(tài)響應的精確解,避免了結構擴階、計算指數(shù)矩陣以及矩陣求逆.結合虛擬激勵法,快速求解在非平穩(wěn)地震激勵下(含均勻調制非平穩(wěn)和非均勻調制非平穩(wěn))耗能減震結構的非平穩(wěn)響應.

        1 結構運動方程及其系統(tǒng)響應

        1.1 實用粘彈性模型阻尼器本構關系

        實用粘彈性阻尼器[PQ(t)]的分析模型如圖1所示,該模型由兩支Maxwell單元并聯(lián)一支Kelvin單元組成,該實用粘彈性阻尼器受力[PQ(t)]與相對位移[xQ]的本構關系為[3]:

        [PQ(t)=0tQ(t-τ)xQ(τ)dτ=k0xQ(t)+c0xQ+P0Q(t)] (1)

        [P0Q(t)=0thQ(t-τ)xQ(τ)dτ] (2)

        [hQ(t)=Q(t)-Q(+∞)=k1e-k1c1t+k2e-k2c2t] (3)

        式中:[Q(t)]為實用粘彈性阻尼器的松弛模量函數(shù),[k0]為實用粘彈性阻尼器的平衡剛度,[c0]為實用粘彈性阻尼器的平衡阻尼,[hQ(t)]實用粘彈性阻尼器的松弛函數(shù);[k1、k2]和[c1、c2]分別表示實用粘彈性阻尼器各支Maxwell單元的剛度和阻尼.

        1.2 結構原始運動方程

        設m、k、c分別為帶支撐六參數(shù)實用粘彈性阻尼器單自由度耗能減震結構的質量、剛度和阻尼;將支撐與實用粘彈性阻尼器的整體串聯(lián)系統(tǒng)作為等效阻尼器[PG(t)].在地震動[xg(t)]的作用下,結構相對于地面的位移為[x],結構模型如圖2所示,其運動方程為:

        [mx+cx+kx+PG(t)=-mxg(t)] (4)

        根據(jù)前期研究,等效阻尼器受力[PG(t)]與相對位移[x]的本構關系為[19]:

        [PG(t)=kGx(t)+0 thG(t-τ)x(τ)dτ] (5)

        式中:[kG]和[hG(t)]分別為等效阻尼器[PG(t)]的平衡剛度、松弛函數(shù),可根據(jù)文獻[19]直接得出.

        根據(jù)前期研究,由傳遞函數(shù)法可直接求解出耗能減震系統(tǒng)響應解析式可統(tǒng)一表達為[19]:

        [Slt=j=15ρljyjt , l=1, 2, 3] (6)

        其中:[ρlj]——相應的響應系數(shù),其值可由文獻[19]直接得出;[S1t]——結構的位移響應,[S2t]——結構的速度響應,[S3t]——等效阻尼器受力.

        [yj(t)]為標準一階系統(tǒng)對地震激勵的響應,其表達式為:

        [yj(t)-Sjyj(t)=-xg] (7)

        [yj(t)=-0teSj(t-τ)xg(τ)dτ] (8)

        傳遞函數(shù)法的最大優(yōu)點是不需要擴階,本節(jié)應用傳遞函數(shù)法,直接在原結構基礎上獲得設置支撐六參數(shù)實用粘彈性阻尼器耗能減震結構在任意激勵和非零初始條件下位移、速度和阻尼器受力時域瞬態(tài)響應的精確解,避免了因擴階方程組物理意義不明確,變量個數(shù)劇增,計算效率低等所帶來的缺陷.

        2 耗能減震系統(tǒng)結構非平穩(wěn)響應的快速分析

        為考慮地震的強度非平穩(wěn)和頻率非平穩(wěn),采用Priestley提出的演變譜模型[20]來分析均勻調制和非均勻調制非平穩(wěn)地震響應,它可以用下式表示:

        [xg(t)=-∞+∞a(ω,t)eiωtdα(ω)] (9)

        其協(xié)方差函數(shù)可表示為:

        [Cxg(t1,t2)=-∞+∞eiω(t1-t2)a(ω,t1)a?(ω,t2)Sxf(ω)dω] (10)

        式中:[Sxf(ω)]為0均值平穩(wěn)隨機過程的自譜密度;[i=-1];“*”表示取復共軛;[α(ω)]是一個正交增量過程;[a(ω,t)]是一滿足[a(ω,t)=a?(-ω,t)]的調制函數(shù).

        由式(6)知,耗能減震系統(tǒng)的一般響應[Sl(t)](含結構位移、速度及阻尼器受力)的非平穩(wěn)協(xié)方差函數(shù)為:

        [E[Sl(t)Sl(t+τ)]=j=15k=15ρljρ?lkE[yj(t)y?k(t+τ)]] (11)

        [E[yj(t)y?k(t+τ)]=0t0t+τeSj(t-η) eS?k(t+τ-ξ)Cxg(η,ξ)dηdξ] (12)

        將式(10)代入式(12)可得:

        [E[yj(t)y?k(t+τ)]=-∞+∞ Yj(ω,t)Y?k(ω,t+τ)Sxf(ω)dω] (13)

        [Yj(ω,t)=0t esj(t-η) a(ω,η)eiωηdη , (j=1, 2, …, 5)] (14)

        式(14)為標準一階系統(tǒng)對激勵[a(ω,t)eiωt]的響應.根據(jù)林家浩等[18]提出的虛擬激勵法原理,若構造如下虛擬激勵:

        [g(t)=Sxf(ω)a(ω,t)eiωt] (15)

        則產(chǎn)生的虛擬響應必定為:

        [Yj(ω,t)=Sxf(ω)Yj(ω,t)] (16)

        則由式(16)和式(13)可分別求得:

        [Yj(ω,t)Y?k(ω,t+τ)=Yj(ω,t)Y?k(ω,t+τ)Sxf(ω)] (17)

        [E[yj(t)y?k(t+τ)]=-∞+∞ Yj(ω,t) Y?k(ω,t+τ)dω] (18)

        由式(14)—式(16)知:

        [Yj(ω,t)=0teSj(t-η)g(η)dη] (19)

        對式(19)取時間步長[Δt=ti+1-ti],進行數(shù)值離散化得:

        [Yj(ω,ti+1)=TYj(ω,ti)+titi+1eSj(ti+1-η)g(η)dη] (20)

        式中:[T=eSjΔt].

        在很小的時間步長[Δt]內(nèi),通常可以認為激勵是線性變化的,則式(20)可表示為:

        [Yj(ω,ti+1)=TYj(ω,ti)+A1g(ti)+A2g(ti+1)] (21)

        式中:

        [A1=1-TS2jΔt+TSj;A2=T-1S2jΔt-1Sj] (22)

        由式(21)可看出,在零初始條件的情況下,即[Yj(ω,0)=0],可以推導出[ti=iΔt]時的[Yj(ω,ti)]表達式:

        [Yj(ω, t1)=A1g(t0)+A2g(t1)] (23)

        [Yj(ω,ti)=Ti-1A1g(t0)+Ti-2A2g(t1)+…+T0A3g(ti-1)+A2g(ti) , i≥2] (24)

        式中:

        [A3=TA2+A1] (25)

        上式中,若用[Bi,0, Bi,1, …, Bi,i]表示[g(t0) , g(t1) , …, g(ti)]的系數(shù),則式(24)可表示為:

        [Yj(ω,ti)=WiJi] (26)

        其中:

        [Wi=Bi,0Bi,1…Bi,i ] (27)

        [Ji=g(t0)g(t1)…g(ti) T] (28)

        由以上推導可知,[Yj(ω,ti)]對應系數(shù)的計算是一個遞推過程,[ti]時刻所對應的系數(shù)[Bi,0, Bi,1, …, Bi,i]只和結構自身有關,且可用[ti-1]時所對應的系數(shù)[Bi-1,0 , Bi-1,1, …, Bi-1,i-1]表示.

        [B1,0=A1B1,1=A2 , i=1 ] (29a)

        [B2,0=TB1,0 B2,1=TA2+A1 , i=2 B2,2=A1,1 ] (29b)

        [Bi,0=TBi-1,0 Bi,1=TBi-1,1 , 3≤i Bi,l=Bi-1,l-1 , 2≤l≤i ] (29c)

        根據(jù)式(29)建立的系數(shù)遞推關系式,各時刻對應的系數(shù)均可通過遞推式得出.本文方法未計算指數(shù)矩陣,不需求逆計算,與文獻[21]提出的顯式時域分析法相比,縮短了計算時間,提高了計算效率.

        將式(26)代入式(18)便可以簡潔的計算出[E[yj(t)y?k(t+τ)]],即:

        [E[yj(t)y?k(t+τ)]=-∞+∞ Yj(ω,t) Y?k(ω,t+τ)dω] (30)

        地震工程中,一般研究結構系統(tǒng)的響應方差,即[τ=0]時,上式可簡化為:

        [σ2yy=-∞+∞ Yj(ω,t) Y?j(ω,t)dω] (31)

        上式為無窮廣義積分,可根據(jù)數(shù)值積分方法來進行求解.在實際計算中,為考慮積分的計算問題,積分上下限一般取有限值,故而設積分區(qū)間為[-b,b].若采用等間距梯形積分公式來計算,式(31)可寫為:

        [σ2yy=Δωn=0pYj(ωn,t)Y?j(ωn,t)] (32)

        式中:[p]為離散頻點數(shù),[Δω=2b/p],[ωn=nΔω(n=0,1,…,p)].

        3 算例

        對如圖2所示設置帶支撐的六參數(shù)實用粘彈性阻尼器的某單層鋼結構建筑進行地震響應分析,結構的基本參數(shù)為:質量[m=42 500 kg],剛度[k=145.43×105 N/m],結構基本周期[T0=0.339 s],阻尼比[S0]分別取0.02、0.04、0.08、0.2;支撐剛度[kb=1.5 k],實用粘彈性阻尼器的基本參數(shù)為:Kelvin單元的平衡剛度和阻尼分別為[k0=0.36×105 N/m]、[c0=0.37×105(N?s)/m],Maxwell單元阻尼器兩分支單元的剛度和阻尼分別為[k1=42.08×105 N/m],[c1=0.83×105(N?s)/m];[k2=6.87×105 N/m],[c2=2.15×105(N?s)/m].平穩(wěn)地震動[xf(t)]譜密度函數(shù)取為Kanai-Tajimi譜:

        [Sxf(ω)=ω4f+4ξ2fω2fω2(ω2-ω2f)2+4ξ2fω2fω2?S0]

        其計算取值為:[ωf=10.9 rad/s],[ξf=0.96];[S0=0.015 54 m2/s3].

        調制函數(shù)分別取為階躍型和非均勻調制函數(shù)Spanos-Solomos型[22],計算參數(shù)取為:

        [a(t)=1]

        [a(ω,t)=a1(ω,t)=2 ω5π t e-0.5×0.15+ω225π2t]

        在階躍型調制非平穩(wěn)地震激勵作用下,結構的位移、速度和阻尼器受力隨時間的響應方差如圖3—圖5所示.由計算結果可以看出:在均勻平穩(wěn)激勵下,結構位移、速度和阻尼器受力在經(jīng)過一個時間點后會趨于穩(wěn)定;阻尼比越大結構的3種響應均越?。蛔枘岜仍酱?,達到平穩(wěn)值的時間越短;同一阻尼比下,結構速度最先達到平穩(wěn),阻尼器阻尼力速度其次,結構位移最后達到平穩(wěn)值.

        在Spanos-Solomos型非均勻調制非平穩(wěn)地震激勵作用下,結構的位移、速度和阻尼器受力響應方差如圖6—圖8所示.由計算結果可以看到:在非均勻激勵下,結構的位移、速度和阻尼力受力具有峰值效應,且同一結構響應的峰值發(fā)生時間不隨阻尼比的變化而變化;同一阻尼比下,速度最先達到峰值,阻尼器阻尼力其次,結構位移最后達到峰值;阻尼比越大3種響應對應的值就越小.

        為了研究支撐剛度對結構響應的影響,研究了結構基本參數(shù)不變時,而支撐剛度分別為0、0.5 k、1 k、1.5 k、2 k,阻尼比取[S0=0.1]時,在2類非平穩(wěn)激勵下的3種響應.

        階躍型調制非平穩(wěn)地震激勵作用下,結構的位移、速度和阻尼器受力隨時間的響應方差如圖9—圖11所示.由計算結果可以發(fā)現(xiàn):在均勻平穩(wěn)激勵下,結構位移、速度和阻尼器受力在經(jīng)過一個時間點后會趨于穩(wěn)定;支架的剛度越大,結構的3種響應均越??;支架剛度達到結構剛度的1.5倍以上時,結構位移和速度響應隨支架的變化基本趨于穩(wěn)定;支架剛度對于阻尼器阻尼力的影響較小.

        在Spanos-Solomos型非均勻調制非平穩(wěn)地震激勵作用下,結構的位移、速度和阻尼器受力響應方差如圖12—圖14所示.由計算結果可以看到:在非均勻激勵下,結構的位移、速度和阻尼力受力具有峰值效應,且同一結構響應的峰值發(fā)生時間不隨支架剛度的變化而變化;支架剛度達到結構剛度的1.5倍以上時,結構位移和速度響應隨支架的變化基本趨于穩(wěn)定;支架剛度對于阻尼器阻尼力的影響較小.

        4 結論

        本文結合虛擬激勵法對設置帶支撐實用粘彈性阻尼器單自由度耗能減震系統(tǒng)基于非平穩(wěn)地震響應作用下進行了研究.研究表明:

        1)在平穩(wěn)激勵作用下,阻尼器阻尼越大,結構響應越快達到平穩(wěn),說明阻尼對于減震的有效性.

        2)在平穩(wěn)激勵作用下,其他條件一致下,結構速度響應最先達到平穩(wěn),位移最慢達到平穩(wěn),說明阻尼對于結構響應的不同步性.

        3)在非平穩(wěn)激勵作用下,結構的3種響應均存在峰值效應,且結構速度響應最先達到峰值,結構位移最晚達到峰值,但峰值發(fā)生時間不隨阻尼而改變,說明在非平穩(wěn)激勵下,結構的響應與結構頻率有關,且阻尼對于結構不同影響存在相位差.

        4)支撐剛度對于結構的位移和速度的響應影響較大,對阻尼器阻尼力影響較小.

        參考文獻

        [1]CHRISTOPOULOS C,F(xiàn)ILIATRULT A.Principle of passive supplemental damping and seismic isolation[M]. Pavia:IUSS Press,2006.

        [2]中華人民共和國住房和城鄉(xiāng)建設部.建筑抗震設計規(guī)范:GB 50011-2010[S].北京:中國建筑工業(yè)出版社,2010.

        [3]周云.粘彈性阻尼器減震結構設計[M].武漢:武漢理工大學出版社,2006.

        [4]MAZZA F,VULCANO A.Control of the earthquake and wind dynamic response of steel-framed buildings by using additional braces and/or viscoelastic dampers[J]. Earthquake Engineering and Structural Dynamics,2011,40(2):155-174.

        [5]汪志昊,陳政清.高層建筑結構中粘滯阻尼器的新型安裝方式[J].世界地震工程,2010,26(4):135-140.

        [6]YAMADA K.Dynamic characteristics of SDOF structure with Maxwell element[J].Journal of Engineering Mechanics,2008,134(5):396-404.

        [7]PALMERI A,RICCIARDELLI F. State space formulation for linear viscoelastic system with memory[J].Journal of Engineering Mechanics,2003,129(7):715-724.

        [8]SINGH M P,VERMA N P. Seismic analysis and design with Maxwell dampers[J].Journal of Engineering Mechanics,2003,129(3):273-282.

        [9]李創(chuàng)第,李暾,尉宵騰,等. Maxwell阻尼耗能結構非平穩(wěn)地震響應解析分析[J].振動與沖擊,2016,35(19):172-180.

        [10]TRINDADE M,BENJEDDOU A,OHAYON R. Modeling of frequency-dependent viscoelatic materials for active passive vibration damping[J]. Journal of Vibration and Acoustics,2002,122:169-174.

        [11]SORRENTINO S,F(xiàn)ASANA A. Finite element analysis of vibration linear systems with fractional derivative viscoelastic modes[J]. Journal of Sound and Vibration,2007,299:839-853.

        [12]SOONG T T,DARGRUSH G F. Passive energy dissipation systems in structural engineering[M]. England:John Wiley and Ltd,1997.

        [13]OU J P,LONG X,LI Q S. Seismic response analysis of structures with velocity-dependent dampers[J].Journal of Constructional Steel Research,2007,63:628-638.

        [14]林家浩.隨機振動的虛擬激勵法[M].北京:科學出版社,2004.

        [15]鐘萬勰.結構動力方程的精細時程積分法[J].大連理工大學學報,1994(2):131-136.

        [16]SU C,HUANG H,MA H. Fast equivalent linearization method for nonlinear structures under nonstationary random excitations[J]. Journal of Engineering Mechanics,2016,142(8):04016049.

        [17]張書天,方同,孫木楠.非均勻調制激勵性系統(tǒng)的非平穩(wěn)響應分析[J].力學學報,1996,28(1):104-108.

        [18]DEODATIS G,SHINOZUKA M.Auto-regressive model for non-stationary stochastic processes[J].Journal of Engineering Mechanics,1988,114 (11):1995-2012.

        [19]李創(chuàng)第,劉鵬,葛新廣,等.六參數(shù)實用粘彈性阻尼器單自由度減震系統(tǒng)非平穩(wěn)響應分析[J].廣西科技大學學報,2018,29 (2):110-118,125.

        [20]PRIESTLEY M B.Power spectral analysis of non-stationary random processes[J].Journal of Sound and Vibration,1967,6(1):86-97.

        [21]蘇成,徐瑞.非平穩(wěn)激勵下結構隨機振動時域分析法[J].工程力學,2010,27(12):77-83.

        [22]SPANOS P T D,SOLOMOS G P. Markov approximation to transient vibration[J]. Journal of Engineering Mechanics,1983,109(4):1134-1150.

        Fast solution to non-stationary seismic response of energy dissipation structure with practical viscoelastic damper

        LI Chuangdi, ZHU Tengfei, BAI Dalian, GE Xinguang

        (School of Civil Engineering And Architecture, Guangxi University of Science and Technology, Liuzhou 545006, China)

        Abstract: In order to improve the response analysis efficiency of energy dissipation structure of practical viscoelastic dampers based on non-stationary seismic excitation, the numerical analysis method of random seismic response of practical viscoelastic dampers with single degree of freedom energy dissipation system is systematically studied. Firstly, a six-parameter practical viscoelastic damper with braces is used to model the behavior of viscoelastic dampers. Then, the transient response of the energy dissipation structure is solved by the transfer function method. Finally, based on the virtual excitation method, the time-domain transient responses of energy dissipation structure under unsteady excitation are obtained and two examples including uniform or non-uniform modulated random excitation are given. The examples show that the method has the characteristics of high efficiency and strong engineering application for the calculation of this kind of problem. This provides a reference for the rapid solution of the viscoelastic energy dissipation structure under non-stationary seismic excitation.

        Key words: practical viscoelastic damper; energy absorbing structure; non-stationary seismic response; damper force; fast solution

        (學科編輯:黎 婭)

        欧美成人看片一区二区三区尤物| av天堂手机在线免费| av在线一区二区三区不卡| 日韩内射美女片在线观看网站| a级大胆欧美人体大胆666| 人妻久久999精品1024| 日本加勒比一区二区在线观看| 亚洲精品1区2区在线观看| 亚洲av成人无码精品电影在线| 久久精品无码一区二区乱片子| 超短裙老师在线观看一区| 一本色道久久88加勒比一| 欧美一区二区三区红桃小说| 狠狠躁夜夜躁AV网站中文字幕| av手机天堂在线观看| 午夜福利理论片在线观看播放| 日韩av精品国产av精品| 亚洲高潮喷水中文字幕| 麻豆国产精品伦理视频| 欧美人与善在线com| 欧美丰满大屁股ass| 久久久精品中文无码字幕| 中文字幕综合一区二区三区| 奇米影视第四色首页| 亚洲日韩专区在线视频| 黄色三级国产在线观看| 亚洲综合极品美女av| 一本一道av无码中文字幕| 国产成人综合日韩精品无| 99青青草视频在线观看| 综合色区亚洲熟妇另类| 欧美日韩国产免费一区二区三区欧美日韩| 精品一区二区中文字幕| 国产精品18久久久白浆| 97久久精品午夜一区二区| 在线无码精品秘 在线观看| 久久这里都是精品99| 午夜精品久久久久久99热| 动漫在线无码一区| 亚洲六月丁香色婷婷综合久久| 人妻尝试又大又粗久久|