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

斜拉橋面內(nèi)豎向固有振動模型及特性影響的有限差分分析

2023-07-31 07:41:26陳柯帆李源賀拴海王康殷怡萍宋一凡
關(guān)鍵詞:模態(tài)振動結(jié)構(gòu)

陳柯帆 ,李源 ,2?,賀拴海 ,2,王康 ,殷怡萍 ,宋一凡 ,2

(1.長安大學(xué) 公路學(xué)院,陜西 西安 710064;2.舊橋檢測與加固技術(shù)交通行業(yè)重點實驗室(長安大學(xué)),陜西 西安 710064)

斜拉橋美觀、經(jīng)濟、跨越能力強,近年來備受橋梁工程師青睞[1].與此同時,斜拉橋整體結(jié)構(gòu)復(fù)雜,柔度大、阻尼低、剛度不足導(dǎo)致其非線性行為極為突出[2].尤其當(dāng)拉索局部模態(tài)與斜拉橋整體模態(tài)頻率比處于“1∶2”或“1∶1”等固定比例區(qū)間時[3-6],容易在環(huán)境激勵下,引發(fā)拉索劇烈振動,給橋梁安全運營帶來了極大隱患.因此,除了目前僅有的有限元方法外[7-9],如何建立準(zhǔn)確的斜拉橋整體動力模型來便捷而又準(zhǔn)確地計算斜拉橋豎向整體模態(tài)參數(shù),對推動斜拉橋應(yīng)用與發(fā)展至關(guān)重要[10].

近年來,國內(nèi)外學(xué)者圍繞斜拉橋整體動力學(xué)建模方法,尤其對于斜拉橋主梁在多點彈性支撐作用下的力學(xué)行為模擬與分析做了大量研究工作.吳慶雄等[11]進行了單索-梁結(jié)構(gòu)和二索-梁結(jié)構(gòu)模型固有振動試驗,建立了多索-梁結(jié)構(gòu)動力學(xué)模型,討論了斜拉索對索梁結(jié)構(gòu)面內(nèi)固有振動特性的影響;Cao等[12]和李專干等[13]通過建立主梁的分段函數(shù),將主梁和拉索等效為若干獨立梁段,基于拉索錨固處的邊界條件得到了剛塔柔梁斜拉橋整體動力學(xué)建模的運動方程,討論了結(jié)構(gòu)對稱性對動力特性的影響;趙文忠等[14]通過分段函數(shù)求解了三索結(jié)構(gòu)的動力方程,研究了拉索一階和二階頻率比條件對共振的影響;Cong 等[15]、Kang 等[16-18]、蘇瀟陽等[19]建立了多梁彈簧動力學(xué)模型,運用傳遞矩陣法給出了不同體系下斜拉橋整體動力學(xué)模型動力學(xué)微分控制方程,提出了不同體系斜拉橋豎向剛度評估方法.在此基礎(chǔ)上,該課題組還對索拱結(jié)構(gòu)[20-21]、懸索結(jié)構(gòu)[22]的面內(nèi)外固有振動模態(tài)參數(shù)進行了系統(tǒng)分析與研究.

受限于現(xiàn)有動力學(xué)建模方法與整體模態(tài)參數(shù)計算過程的煩冗,目前關(guān)于斜拉橋整體豎向模態(tài)的理論計算方法大多高度簡化甚至忽略索力水平投影、振動時拉索間的耦合影響作用、主梁彎曲剛度等因素對結(jié)構(gòu)振動特性的影響,或是迭代解析方法煩瑣、復(fù)雜,不利于在實際斜拉橋工程中推廣與應(yīng)用.

針對此問題,本文通過離散斜拉橋多點彈性支承梁的集中質(zhì)量參數(shù)體系,建立了一種新的斜拉橋面內(nèi)豎向整體動力學(xué)模型.該模型考慮了斜拉索對主梁的豎向彈性支承作用與對水平梁截面的軸力影響,引入微梁段兩側(cè)的剪力以模擬主梁彎曲剛度、拉索間振動耦合等影響作用,通過微梁段間的彎矩平衡和有限差分法,修正了不同結(jié)構(gòu)體系斜拉橋面內(nèi)豎向整體動力學(xué)模型的運動方程,結(jié)合特征值法給出了斜拉橋面內(nèi)豎向模態(tài)頻率及振型計算方法.通過對比參考文獻中動力學(xué)模型算法案例分析結(jié)果與某斜拉橋的實測值,進一步驗證了本文關(guān)于斜拉橋面內(nèi)豎向運動參數(shù)體系建模方法的適用性和正確性.本文計算方法無須建立大量細化的有限元模型,運用MATLAB、Excel 等軟件按編碼流程即可準(zhǔn)確、快速地估算斜拉橋面內(nèi)豎向模態(tài)頻率及振型,簡化了計算過程,便于工程應(yīng)用.

1 多點彈性支撐梁的離散模型

約定下標(biāo)“B”和“C”分別表示梁和索;下標(biāo)“i”和“j”分別表示索和梁序號(i∈[1,I],j∈[1,J]).為便于區(qū)別有索區(qū)梁段與無索區(qū)梁段,對于Ci#拉索錨固處的有索區(qū)梁段表示為Bji#梁段.建立如圖1(a)所示的多索-主梁模型.由于橋面質(zhì)量遠大于拉索質(zhì)量,本文忽略了拉索振動對主梁振動的影響,將Ci#拉索豎向視為Bji#梁段的彈性支承,水平向的索力投影等效為Bji#梁段的軸向荷載HBCi,定義He為支座水平力,如圖1(b)所示.為了精確模擬和求解具有分布質(zhì)量、荷載和多點彈性支撐作用下斜拉橋主梁的動力行為,不考慮主梁的縱向運動,本文將主梁進一步簡化為J個間距d相同、彼此鉸接、帶有I個豎向彈力支承和I個不同軸力的集中質(zhì)量參數(shù)體系,主梁的質(zhì)量和荷載均被視為作用在這些集中質(zhì)量點之上,獨立梁段間彼此通過理想鉸連接,如圖1(c)所示.

圖1 斜拉橋主梁離散模型的簡化流程Fig.1 The reduction process of main beam of cable-stayed bridges

圖1(b)中,Ci#拉索對Bji#梁段豎向彈性支承系數(shù)kBCi為[7,10,23-24]:

式中:ECi、ACi分別為Ci#拉索的彈性模量及截面積;θCi表示Ci#拉索軸線與主梁大里程方向夾角;lCi表示 Ci#拉索上下端錨固點軸向距離.由于拉索振動中的索力增量遠小于拉索初始索力,對于拉索索力的水平投影本文僅考慮初始索力[25-26].因此,主梁上水平軸力HBCi和支座水平力HBe(下標(biāo)e表示邊界)滿足:

式中:SCi表示Ci#拉索初始索力.圖1(c)中,獨立梁段彼此鉸接,通過引入離散梁段的左右側(cè)剪力以模擬具有分布質(zhì)量的主梁豎向彎曲剛度、拉索間振動耦合作用在振動過程中的相互影響.梁段間的受力如圖2(a)所示,梁段處的受力如圖2(b)和(c)所示.

圖2 微梁段受力示意圖Fig.2 The force schematic of the micro-beam segment

圖2 中,F(xiàn)B(j-1,j)、FB(j,j+1)表示梁段左右兩側(cè)的剪力;NB(j-1,j)、NB(j,j+1)表示梁段左右兩側(cè)的軸力-表示相鄰梁段間的節(jié)段左右兩側(cè)受到的彎矩作用表示相鄰梁段運動夾角;aBj表示梁段運動加速度.考慮系統(tǒng)初始為平衡狀態(tài),基于D’Alembert 原理可以得到Bj#梁段在豎向的動力平衡方程:

式中:VBj(t)表示Bj#梁段與時間相關(guān)的豎向振動位移變化因子,后文中簡寫為VBj;kT表示斜拉橋主塔對主梁面內(nèi)豎向自由運動的剛度彈簧系數(shù),需依據(jù)圖紙和規(guī)范,以及不同結(jié)構(gòu)體系下塔-梁處的邊界條件進行取值.δ(j-ji)為狄拉克(Dirac)函數(shù),由式(5)~式(6)定義:

值得注意的是,靠近邊界的梁段需根據(jù)結(jié)構(gòu)體系邊界條件進行求解運算,如附表1所示.

附表1 不同主梁邊界條件下的振動方程系數(shù)表達式Additional Table 1 Expressions of vibration equation coefficients under different boundary conditions of main beams

2 數(shù)學(xué)表達與驗證

2.1 基于有限差分法的方程優(yōu)化

假設(shè)質(zhì)量體系分布較密,梁段間相對位移較小,則其振動的幾何關(guān)系滿足以下關(guān)系式:

圖2(a)中,顯然梁段間存在剪力與彎矩平衡:

由上式可得梁段間剪力、彎矩與振動位移間關(guān)系:

對于無索區(qū)梁段(即j≠ji),左右側(cè)截面受到同一方向常軸力影響,可以得到:

對于有索區(qū)梁段(即j=ji),左右側(cè)截面軸力突變量為拉力的水平投影:

整合式(7)~式(15)并代入式(4)后,可以得到Bj#梁段的振動方程:

整合B1#~BJ#梁段方程后可以得到斜拉橋面內(nèi)豎向固有振動方程:

式中:Γ為等效剪力效應(yīng)系數(shù)矩陣,表征了剪力效應(yīng)對主梁運動的影響,與斜拉橋邊界條件有關(guān);Ξ為等效軸力效應(yīng)系數(shù)矩陣,表征了軸力效應(yīng)對主梁運動的影響,與拉索數(shù)量和錨固位置有關(guān),矩陣具體形式如附錄1 所示皆為相同形式的J維列向量.為避免贅述,在此僅展示形式:

式(20)中:ψBj單位與頻率一致,是Bj#梁段的局部模態(tài)頻率方程,表征了Bj#梁段參與整體模態(tài)的模態(tài)頻率,由κDj、κCi、κSj、κT構(gòu)成.κDj表示索力水平投影,即主梁軸力對Bj#梁段局部模態(tài)頻率的影響;κCi表示拉索豎向彈性支承作用的影響;κSj表示主梁軸力的影響,系數(shù)λj與斜拉橋邊界條件有關(guān),具體形式如附表1所示;κT表征了不同結(jié)構(gòu)體系主塔的影響.

2.2 基于特征值法的結(jié)構(gòu)模態(tài)分析

參數(shù)體系的自由振動方程組——式(17)實際上是一個J維的齊次方程組,方程有解的前提是系數(shù)矩陣行列式為0.因此,不計結(jié)構(gòu)阻尼,根據(jù)式(17)構(gòu)造結(jié)構(gòu)特征矩陣形式如下:

式中:MBj表示質(zhì)量對角矩陣;KBj表示主梁的等效剛度矩陣.根據(jù)簡諧振動理論,構(gòu)造頻率表達式:

式中:eig 表示求解矩陣特征值與特征向量.采用以上計算公式和流程求解斜拉橋的面內(nèi)豎向固有振動模態(tài)參數(shù),在確定質(zhì)量矩陣并根據(jù)相應(yīng)結(jié)構(gòu)邊界條件選擇和編輯相應(yīng)軸力與剪力系數(shù)矩陣后,僅需借助MATLAB、Excel 等工具就能簡單地計算和求解固有振動頻率和振型,無須進行大量細化的有限元建模分析.本文依托MATLAB編制了運行算法程序,其流程如圖3所示.

圖3 斜拉橋面內(nèi)豎向固有振動模態(tài)參數(shù)計算流程Fig.3 The computation solution process of a cable-stayed bridge’s in-plane vertical natural vibration modal properties

2.3 參考文獻算例驗證

以2020年Cong 團隊的雙索-梁結(jié)構(gòu)理論解析研究結(jié)果為對比對象[14,28],該研究通過簡化斜拉橋為多索-梁動力學(xué)模型,考慮結(jié)構(gòu)構(gòu)件間的幾何非線性邊界條件,基于傳遞矩陣法研究了多索斜拉結(jié)構(gòu)的固有振動特性與面內(nèi)外振動響應(yīng).其研究模型如圖4所示.

圖4 文獻[14,28]雙索-梁結(jié)構(gòu)示意圖Fig.4 Schematic of a double-cable-beam structure in the references [14,28]

本文代入了文獻[14,28]參數(shù),選取主梁劃分節(jié)段參數(shù)d=1 m,即J=300,I=2,運用MATLAB 軟件根據(jù)圖3 計算流程編寫計算程序,討論和對比多索斜拉結(jié)構(gòu)的面內(nèi)固有振動頻率及以主梁為主要振型的前五階振動模態(tài),如表1所示.

表1 兩索-梁結(jié)構(gòu)面內(nèi)豎向固有振動頻率Tab.1 In-pane vertical natural vibration frequencies of the two-cable-beam structure Hz

根據(jù)表1,本文解析法得到的頻率數(shù)值平均絕對誤差僅為0.2‰,精確度超過了原文解析方法絕對誤差1.4%.分別按照本文有限元數(shù)值模擬與解析法求解該結(jié)構(gòu)以主梁為主的前五階固有振型,如圖5所示.

圖5 兩索-梁結(jié)構(gòu)的豎向振型Fig.5 Vertical modal shapes of the two cable-beam structure

圖5 顯示兩種方法得到的結(jié)構(gòu)振型一致性良好,進一步驗證了本文計算方法的有效性和準(zhǔn)確性.

3 斜拉橋固有振動特性及影響性分析

3.1 斜拉橋模態(tài)分析

以我國西北地區(qū)某混凝土斜拉橋為對象開展固有振動特性及其影響參數(shù)分析.該橋全長166.8 m(39 m+88.8 m+39 m),采用三跨雙臺、雙塔、雙索面對稱布置,墩塔處固結(jié),為半漂浮支承體系,其立面圖如圖6所示.鋼筋混凝土主梁由節(jié)段預(yù)制雙箱梁和預(yù)制行車道板組合形成,箱梁高1.2 m,橋面凈寬8.5 m.在兩箱梁間錨固板處設(shè)橫系梁一道,縱向長約0.22 m,為方便引用,匯總主梁各截面參數(shù)設(shè)置如表2 所示;全橋現(xiàn)有48 根斜拉索,從小里程邊跨至大里程邊跨方向以C1#~C24#對單索面拉索依次編號,參數(shù)如表3所示(僅示出一側(cè),另一側(cè)參數(shù)與之相近,斜拉索彈性模量經(jīng)恩斯特公式修正后取200 GPa).

表2 主梁參數(shù)設(shè)置表Tab.2 Parameters of the beam

表3 斜拉索參數(shù)設(shè)置表Tab.3 Parameters of stay cables

圖6 橋梁立面圖(單位:cm)Fig.6 Bridge elevation drawing(unit:mm)

根據(jù)橋梁結(jié)構(gòu)形式,采用商業(yè)有限元軟件對該橋進行動力特性分析,得到該結(jié)構(gòu)前三階自振頻率、振型特征,如圖7所示.

圖7 有限元法得到的斜拉橋前三階固有振動模態(tài)參數(shù)Fig.7 The first-three order natural modes of the cable-stayed bridge by the finite element model

為進一步對比和驗證本文解析公式的正確性,采用實橋?qū)崪y、有限元分析、本文解析三種方法對該橋進行自振特性分析.設(shè)置參考算例參數(shù)詳情如下所示:

1#算例(Referred Case 1#,RC1),實體結(jié)構(gòu)有限元法:根據(jù)實際橋梁結(jié)構(gòu)形式,采用商業(yè)有限元軟件對該橋進行動力特性分析,識別該結(jié)構(gòu)自振頻率、振型特征等.

2#算例(Referred Case 2#,RC2),現(xiàn)場實測法:在橋梁邊跨0.4L(L表示跨徑)截面及中跨跨中截面布設(shè)加速度傳感器,應(yīng)用脈動激勵法進行橋梁結(jié)構(gòu)的振動試驗,識別大橋前3 階整體模態(tài)的動力特性參數(shù),采用DHSAS 頻譜分析及模態(tài)分析軟件對其進行快速傅里葉變換得到相應(yīng)的功率譜圖,再對其作進一步的頻譜分析可得到橋梁結(jié)構(gòu)的自振頻率、阻尼比.現(xiàn)場動載試驗布置如圖8所示.

圖8 現(xiàn)場動載試驗Fig.8 On-site dynamic load test

3#算例(Referred Case 3#,RC3),簡化模型解析法:根據(jù)本文動力簡化模型及振動方程,塔梁連接處按照圖紙取kT=5.3 × 109N/m,大小里程結(jié)合墩采用簡支邊界條件,位于主梁同一截面的雙索考慮為動力彈簧的并聯(lián)關(guān)系,基于式(17)采取微梁段長度d1=0.1 m,運用MATLAB 軟件根據(jù)圖3 計算流程編寫計算程序,對該結(jié)構(gòu)的面內(nèi)豎向模態(tài)頻率和振型進行了計算和分析.

匯總以上工況下得到的該橋以主梁面內(nèi)豎向為主要振型的前三階自振頻率,如表4所示.

表4 斜拉橋面內(nèi)豎向固有振動頻率Tab.4 In-plane vertical natural vibration frequencies of the cable-stayed bridge Hz

表4 中,若以實橋測得的模態(tài)參數(shù)為標(biāo)準(zhǔn),有限元法得到的結(jié)果的絕對誤差平均值為3.3%,而本文解析法得到的結(jié)果的絕對誤差平均值為2.7%.此外,相較于其他兩種結(jié)果,實測頻率值整體偏小,這是由于該橋建設(shè)時間較長,結(jié)構(gòu)剛度在通行運營中有所下降.匯總RC1 和RC3 工況下得到的該斜拉橋前三階豎向固有振型,如圖9所示.

圖9 RC1、RC3工況下的斜拉橋前三階面內(nèi)豎向振型Fig.9 The first three order modes of cable-stayed bridge under the conditions of RC1 and RC3

圖9 中兩種工況下結(jié)構(gòu)面內(nèi)豎向前三階振型一致性良好,上述情況進一步說明了本文方法計算斜拉橋面內(nèi)豎向固有振動模態(tài)參數(shù)的精確性和適用性.

3.2 梁軸力對斜拉橋豎向模態(tài)的影響

為研究索力水平投影為梁提供的軸向力對斜拉橋面內(nèi)豎向振動模態(tài)參數(shù)的影響,引入μa表示軸向力對Bj#梁段動平衡方程的放大系數(shù),由式(28)定義:

按照本文解析法研究成果,基于RC3參數(shù)設(shè)置,圖10展示了軸力對斜拉橋固有振型的影響,圖11展示了軸力對斜拉橋固有振動頻率的影響.

圖10 軸力對斜拉橋豎向固有振型的影響Fig.10 Influence of the axial force on the vertical natural modes of the cable-stayed bridge

圖11 梁軸力對斜拉橋面內(nèi)豎向頻率的影響Fig.11 Influence trends of the beam’s axial force on the in-plane vertical frequencies of the cable-stayed bridge

圖12 μa與結(jié)構(gòu)前五階頻率關(guān)系Fig.12 Relation between the first-five order in-plane vertical frequencies with μa

圖12 顯示,μa<28 時,V1~V5 隨軸力增大而呈現(xiàn)線性緩慢下降.當(dāng)μa達到30 左右時,V1 發(fā)生頻率躍遷現(xiàn)象,此時盡管V2頻率數(shù)值低于V1,但V1的振型是結(jié)構(gòu)的一階振型.當(dāng)μa達到34.4時,V1~V3值產(chǎn)生共軛對稱解,其頻域信號值幅值相同而相位不同,而V2 在μa接近38.7 時急速下降為0,此時斜拉橋一階固有振動頻率為0,結(jié)構(gòu)整體失穩(wěn).表明從軸力角度考慮,該橋目前暫無整體失穩(wěn)風(fēng)險,只有當(dāng)軸力達到現(xiàn)有軸力38 倍以后會出現(xiàn)整體失穩(wěn)現(xiàn)象.此外,隨著結(jié)構(gòu)軸力的增大,結(jié)構(gòu)基頻將發(fā)生躍遷,結(jié)構(gòu)體系內(nèi)易產(chǎn)生非線性內(nèi)共振,需進行更深入的研究.

3.3 斷索對斜拉橋面內(nèi)豎向模態(tài)的影響

從式(20)與式(22)中可發(fā)現(xiàn),斜拉橋拉索為主梁提供了軸力與彈性支撐作用(kBCi),為結(jié)構(gòu)提供了有效剛度.因此,結(jié)構(gòu)的固有振動模態(tài)頻率與每一根拉索息息相關(guān).以此實體結(jié)構(gòu)為背景,模擬單拉索出現(xiàn)斷裂的極端情況對結(jié)構(gòu)固有振動頻率及振型的影響,圖13顯示了該工況下低階頻率變化規(guī)律.

圖13 不同位置的拉索斷裂對結(jié)構(gòu)V1~V5頻率的影響Fig.13 Influence on the V1~V5-order modal frequencies when cables of different positions broke down

拉索的索力水平投影對主梁產(chǎn)生的軸壓力降低了整體剛度,而其豎向彈性支承作用則增大了整體剛度,兩者作用下使得斷索對結(jié)構(gòu)基頻影響較小,如圖13 所示.此外,左右側(cè)索面斷索后影響效應(yīng)變化規(guī)律基本一致,其中中跨拉索斷裂對V1、V2和V5影響較大,而邊跨拉索對V3 和V4 影響較大,對比圖10中結(jié)構(gòu)前五階面內(nèi)豎向振型,表明斷索后的影響效應(yīng)與該拉索對應(yīng)的主梁質(zhì)點的振型參與系數(shù)相關(guān).拉索發(fā)生損傷甚至斷裂會改變結(jié)構(gòu)整體頻率,因此,當(dāng)此情況發(fā)生時需進一步考慮結(jié)構(gòu)因局部-整體模態(tài)耦合發(fā)生非線性共振的問題.

4 結(jié)論

1)本文方法考慮了主梁截面的變剛度、變軸力作用,建立的斜拉橋整體模態(tài)動力學(xué)模型更貼近工程實際結(jié)構(gòu).通過代入已有文獻中理論模型算例參數(shù)并對比其理論計算結(jié)果,本文方法計算結(jié)果平均誤差率為0.2‰;代入某斜拉橋參數(shù)并對比其面內(nèi)豎向固有振動頻率實測值,本文方法計算結(jié)果誤差率為2.7%,進一步說明本文解析方法具有較高的精確度和較好的適用性.

2)斜拉橋低階面內(nèi)豎向固有振動頻率隨著軸力增加而降低,振型基本無變化,當(dāng)軸力增加到一定值后,結(jié)構(gòu)的低階面內(nèi)豎向整體模態(tài)頻率將發(fā)生頻率躍遷、分叉現(xiàn)象.

3)斜拉橋發(fā)生斷索對整體結(jié)構(gòu)固有振動頻率影響較小,斷索對各階頻率值影響效應(yīng)與拉索錨固處主梁的對應(yīng)階次振型參與系數(shù)相關(guān).

4)本文方法為快速計算斜拉橋整體模態(tài)頻率以避免斜拉橋整體-局部模態(tài)耦合而發(fā)生非線性內(nèi)共振問題提供了有效幫助.下一步將考慮拉索、主塔與主梁間幾何非線性邊界條件,通過建立更加細化的斜拉橋整體動力學(xué)模型開展斜拉橋的相關(guān)非線性共振研究.

附錄1

為簡化表達,定義Pp和Pp,q表達式為:

Ξ為推導(dǎo)得到的子項系數(shù)矩陣,如式(a-3)所示,其與拉索數(shù)量、錨固位置有關(guān).不同主梁邊界條件下的振動方程系數(shù)表達式見附表1.

附錄2

系數(shù)矩陣的具體形式與邊界條件相關(guān),本文簡支-簡支表達式如式(b-1)所示.

系數(shù)矩陣的具體形式與拉索數(shù)量、錨固位置、邊界條件相關(guān),本文簡支-簡支表達式如式(b-2)所示.

猜你喜歡
模態(tài)振動結(jié)構(gòu)
振動的思考
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
振動與頻率
論結(jié)構(gòu)
中華詩詞(2019年7期)2019-11-25 01:43:04
中立型Emden-Fowler微分方程的振動性
論《日出》的結(jié)構(gòu)
國內(nèi)多模態(tài)教學(xué)研究回顧與展望
創(chuàng)新治理結(jié)構(gòu)促進中小企業(yè)持續(xù)成長
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
UF6振動激發(fā)態(tài)分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 亚洲伊人久久精品影院| 亚洲第一黄片大全| 国产高清无码麻豆精品| 欧美伊人色综合久久天天| 精品国产一二三区| 2021天堂在线亚洲精品专区| 免费 国产 无码久久久| 五月婷婷丁香综合| 国产免费久久精品99re不卡| 亚洲成a人片77777在线播放| 91精品久久久久久无码人妻| 视频二区亚洲精品| 国产情侣一区二区三区| 东京热av无码电影一区二区| 88av在线| 日本不卡视频在线| 色AV色 综合网站| 亚洲性日韩精品一区二区| 亚洲国产系列| 欧美一级色视频| 亚洲第一成年免费网站| 国产清纯在线一区二区WWW| 在线国产91| 成人国产精品视频频| 亚洲欧洲自拍拍偷午夜色| 美女亚洲一区| 98超碰在线观看| 欧美午夜一区| 91成人在线免费观看| 国产美女人喷水在线观看| 久久不卡精品| 亚洲国产91人成在线| 91网在线| 国产精品夜夜嗨视频免费视频| 日韩午夜伦| 无码视频国产精品一区二区| 亚洲日韩高清无码| 精品在线免费播放| 久久香蕉国产线看观| 99精品影院| 亚洲视频无码| 国产日韩丝袜一二三区| 五月婷婷精品| 中文纯内无码H| 中文字幕免费视频| 欧美日韩精品一区二区视频| 日本欧美中文字幕精品亚洲| 国产麻豆aⅴ精品无码| 日韩在线播放中文字幕| 久久精品这里只有国产中文精品| 高潮爽到爆的喷水女主播视频| 伊在人亚洲香蕉精品播放| 成人免费网站在线观看| 视频一本大道香蕉久在线播放| 亚洲天堂久久新| 99热在线只有精品| 波多野结衣一二三| 成年片色大黄全免费网站久久| 中文字幕精品一区二区三区视频 | 538精品在线观看| 中文字幕在线看视频一区二区三区| 亚洲精品成人片在线播放| 中文国产成人精品久久一| 国产91精品久久| 欧美亚洲一区二区三区导航 | 亚洲国产成人超福利久久精品| 手机看片1024久久精品你懂的| 91色爱欧美精品www| 国产99热| 国产精欧美一区二区三区| 国产午夜不卡| 91娇喘视频| 伊人网址在线| 中字无码精油按摩中出视频| 国语少妇高潮| 91亚洲视频下载| 欧美精品另类| 亚洲天堂.com| 亚洲AV无码一区二区三区牲色| 91精品国产麻豆国产自产在线| 亚洲人成网址| 播五月综合|