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

        ?

        一種單自由度體系解析解及其在車(chē)橋動(dòng)力分析中的應(yīng)用

        2017-04-21 01:06:14杜憲亭王少欽
        振動(dòng)與沖擊 2017年7期
        關(guān)鍵詞:橋梁振動(dòng)

        杜憲亭, 喬 宏, 夏 禾, 王少欽

        (1. 北京交通大學(xué) 土木建筑工程學(xué)院,北京 100044;2. 北京建筑大學(xué) 理學(xué)院,北京 100044)

        一種單自由度體系解析解及其在車(chē)橋動(dòng)力分析中的應(yīng)用

        杜憲亭1, 喬 宏1, 夏 禾1, 王少欽2

        (1. 北京交通大學(xué) 土木建筑工程學(xué)院,北京 100044;2. 北京建筑大學(xué) 理學(xué)院,北京 100044)

        假定相鄰時(shí)刻之間荷載線性變化,推導(dǎo)出低阻尼單自由度振動(dòng)體系的解析解,在此基礎(chǔ)上給出了相應(yīng)的車(chē)橋動(dòng)力相互作用系統(tǒng)建模及求解流程。系統(tǒng)模型分解為車(chē)輛、橋梁兩個(gè)子系統(tǒng),基于部件剛體假定和達(dá)朗貝爾原理推導(dǎo)車(chē)輛子系統(tǒng)運(yùn)動(dòng)方程,采用有限元法建立橋梁子系統(tǒng)模型;借助于振型疊加法將兩個(gè)子系統(tǒng)運(yùn)動(dòng)方程解耦,車(chē)輛子系統(tǒng)非正交阻尼部分的影響以及兩個(gè)子系統(tǒng)間的動(dòng)力相互作用均按非線性虛擬力處理;以一節(jié)4軸客車(chē)勻速通過(guò)32 m簡(jiǎn)支梁為例,分別采用提出的解析解法、Newmark-β法以及高斯精細(xì)積分法進(jìn)行動(dòng)力分析。結(jié)果表明,相對(duì)于Newmark-β法和高斯精細(xì)積分法,解析解法不僅具有高精度特點(diǎn),能顯著提高計(jì)算收斂的積分步長(zhǎng),同時(shí)又能避免計(jì)算復(fù)雜的指數(shù)矩陣,具有良好的工程適用性。

        車(chē)橋系統(tǒng);動(dòng)力相互作用;單自由度體系解析解;振型疊加法;非正交阻尼

        隨著鐵路的迅猛發(fā)展,列車(chē)通過(guò)橋梁時(shí)所產(chǎn)生的動(dòng)力相互作用問(wèn)題成為了研究熱點(diǎn)[1-3]。列車(chē)輪對(duì)與橋上鋼軌之間的接觸屬于高度非線性問(wèn)題的范疇,采用Newmark-β法、Wilson-θ法等經(jīng)典數(shù)值積分方法求解時(shí),需要很小的積分步長(zhǎng)才能保證其收斂性[4-9]。計(jì)算效率是制約車(chē)橋動(dòng)力相互作用分析發(fā)展的重要因素。

        翟婉明提出的新型積分方法計(jì)算效率較高,在車(chē)-線-橋分析中得到了廣泛應(yīng)用,但是,該方法要求的積分步長(zhǎng)極小。YANG等利用動(dòng)力縮聚法,縮減了和車(chē)體相關(guān)的自由度,提高了Newmark-β方法的計(jì)算效率,但是該方法需要構(gòu)造復(fù)雜的相互作用單元,使其應(yīng)用受到了限制。張純等利用Haar小波方法對(duì)車(chē)橋系統(tǒng)耦合振動(dòng)問(wèn)題進(jìn)行了求解,發(fā)現(xiàn)Haar小波方法計(jì)算時(shí)間短、計(jì)算準(zhǔn)確性高,但在其求解過(guò)程中需要將矩陣進(jìn)行擴(kuò)階,計(jì)算過(guò)于復(fù)雜。杜憲亭等[10]綜合精細(xì)積分法和高斯積分法特點(diǎn),提出了一種適用于車(chē)橋耦合振動(dòng)的預(yù)測(cè)-校正求解分析框架,能夠顯著提高分析效率,然而精細(xì)積分法涉及計(jì)算增維的指數(shù)矩陣,過(guò)程較為繁瑣,并且時(shí)間步內(nèi)運(yùn)動(dòng)插值的引入降低了分析的數(shù)值穩(wěn)定性。因此,對(duì)于車(chē)橋動(dòng)力相互作用分析而言,尋找更加高效的積分方法和相應(yīng)的系統(tǒng)建模方法及求解策略具有重要的理論意義和工程價(jià)值。

        本文基于時(shí)間步內(nèi)荷載線性變化的假定,推導(dǎo)出單自由度低阻尼體系振動(dòng)的解析解,在此基礎(chǔ)上給出了相應(yīng)的車(chē)橋動(dòng)力相互作用系統(tǒng)建模及數(shù)值求解流程。將車(chē)橋動(dòng)力相互作用系統(tǒng)模型分解為車(chē)輛、橋梁兩個(gè)子系統(tǒng),基于部件剛體假定和達(dá)朗貝爾原理推導(dǎo)車(chē)輛子系統(tǒng)運(yùn)動(dòng)方程,采用有限元法建立橋梁子系統(tǒng)模型。借助于振型疊加法將兩個(gè)子系統(tǒng)運(yùn)動(dòng)方程解耦,車(chē)輛子系統(tǒng)非正交阻尼部分的影響以及兩個(gè)子系統(tǒng)間的動(dòng)力相互作用均為非線性虛擬力處理。最后,采用解析解法、Newmark-β法以及高斯精細(xì)積分法對(duì)一節(jié)4軸客車(chē)勻速通過(guò)32 m簡(jiǎn)支梁所構(gòu)成的耦合系統(tǒng)進(jìn)行動(dòng)力分析,通過(guò)計(jì)算結(jié)果的對(duì)比驗(yàn)證所建議方法的可靠性。

        1 低阻尼單自由度體系解析解

        低阻尼單自由度體系運(yùn)動(dòng)方程為[11]

        (1)

        式中,q、ω、ξ、f、t分別為體系的位移、圓頻率、阻尼比、外荷載、時(shí)間。

        任意形式的荷載總是可以用有限數(shù)量的直線段去逼近。假定式(1)右端的荷載項(xiàng)在時(shí)間步[t,t+Δt]內(nèi)線性變化,則有

        (2)

        式中:Δt為時(shí)間步長(zhǎng);下角標(biāo)t+Δt為相鄰時(shí)刻;t+τ為兩相鄰時(shí)刻之間的任一瞬間。

        式(2)為典型的二階常系數(shù)微分方程,相應(yīng)的解由兩部分組成

        (3)

        式中,上角標(biāo)s、g分別為特解、通解。

        式(1)對(duì)應(yīng)的特征方程為

        r2(t)+2ξωr+ω2=0

        (4)

        所對(duì)應(yīng)的特征值為

        (5)

        0不是特征根,因此特解為

        (6)

        將式(6)代入式(1),得到

        (7)

        由待定系數(shù)法得到

        (8)

        (9)

        特征方程有一對(duì)共軛復(fù)根,則其通解為

        (10)

        因此,

        qt+τ[c1cos(ωDτ)+c2sin(ωDτ)]·e-ωξτ+A·τ+B

        (11)

        (12)

        c1=qt-B

        (13)

        (14)

        取τ=Δt,得到矩陣形式的相鄰時(shí)刻體系運(yùn)動(dòng)狀態(tài)遞推關(guān)系式

        (15)

        其中,

        a11=e-ωξΔt·[cos(ωDΔt)+ωξsin(ωDΔt)/ωD]

        (16)

        a12=e-ωξΔt·sin(ωDΔt)/ωD

        (17)

        a21=-e-ωξΔt·sin(ωDΔt)·ω2/ωD

        (18)

        (19)

        (20)

        (21)

        (22)

        (23)

        當(dāng)式(2)基本假定成立時(shí),應(yīng)用解析解式(15)逐步求解低阻尼單自由度體系運(yùn)動(dòng)不會(huì)產(chǎn)生數(shù)值穩(wěn)定性問(wèn)題,即時(shí)間步長(zhǎng)Δt的取值不受制于體系自身周期的影響,這與Newmark-β法等數(shù)值方法相比具有很大優(yōu)勢(shì)。將式(15)計(jì)算結(jié)果代入式(1),可以得到體系對(duì)應(yīng)的加速度。

        2 車(chē)橋動(dòng)力相互作用系統(tǒng)建模

        列車(chē)通過(guò)橋梁時(shí)引起結(jié)構(gòu)振動(dòng),而橋梁振動(dòng)又反作用于列車(chē),兩者構(gòu)成了一個(gè)相互作用的非線性耦合振動(dòng)系統(tǒng)見(jiàn)圖1。

        圖1 列車(chē)通過(guò)橋梁Fig.1 A vehicle passing through a bridge

        將車(chē)橋動(dòng)力相互作用系統(tǒng)分解為橋梁、車(chē)輛兩個(gè)子系統(tǒng)。依據(jù)有限元法,矩陣形式的系統(tǒng)運(yùn)動(dòng)微分方程為

        (24)

        (25)

        式中:上角標(biāo)B、V分別代表橋梁、車(chē)輛子系統(tǒng);FVV為車(chē)輛內(nèi)部構(gòu)件,如抗蛇行減振器等,產(chǎn)生的非線性作用力;而FVB、FBV為兩子系統(tǒng)間的動(dòng)力相互作用,反映的是系統(tǒng)間力及位移的協(xié)調(diào)關(guān)系,可以通過(guò)耦合系統(tǒng)運(yùn)動(dòng)的線性或非線性函數(shù)表達(dá),在后繼分析中可以作為非線性虛擬力處理。

        顯然,式(24)、式(25) 不能直接采用式(15)求解。

        鐵路橋梁通常由基礎(chǔ)、墩臺(tái)、支座、梁體和橋面系等部件組成,可以采用梁?jiǎn)卧?、桿單元、剛臂單元、板單元、實(shí)體單元、彈簧單元以及質(zhì)量單元等建立其有限元模型。其中,阻尼矩陣cB可以借助Rayleigh阻尼假定確定。

        采用標(biāo)準(zhǔn)模態(tài)變換

        (26a)

        (26b)

        采用式(26)對(duì)式(24)進(jìn)行坐標(biāo)變換及正交化處理,得到

        (27)

        在精度達(dá)到分析要求的前提下,通常引入如下假定對(duì)車(chē)輛子系統(tǒng)進(jìn)行簡(jiǎn)化:①列車(chē)在橋梁行駛時(shí)間較短,忽略速度波動(dòng)對(duì)動(dòng)力系統(tǒng)的影響;②車(chē)輛由車(chē)體、轉(zhuǎn)向架、輪對(duì)等剛體部件組成,且各組成部件之間通過(guò)一系和二系彈簧、阻尼器相聯(lián);③車(chē)輛橫、豎向振動(dòng)與車(chē)輛縱向振動(dòng)可以分開(kāi)考慮,前者對(duì)橋梁振動(dòng)其控制作用。借助于達(dá)朗貝爾原理可以推導(dǎo)出每節(jié)車(chē)的獨(dú)立運(yùn)動(dòng)方程。因此,式(25)可以進(jìn)一步寫(xiě)為

        (28)

        式中:l為通過(guò)車(chē)輛的節(jié)數(shù);mV,j通常為對(duì)角陣。

        式(28)不能直接應(yīng)用振型疊加法,這是由于車(chē)輛的阻尼矩陣關(guān)于振型并不正交。將方程左側(cè)阻尼矩陣力移動(dòng)至方程的右側(cè)作為虛擬力,同時(shí)引入等效阻尼比的概念,再借助于振型疊加法得到

        (29)

        等效阻尼比可參考橋梁子系統(tǒng)取值,也可采用式(30)確定

        (30)

        需要指出的是,采用輪軌非密貼模型時(shí)車(chē)輛處于自由運(yùn)動(dòng)狀態(tài),存在零頻率和剛體模態(tài),在模態(tài)求解時(shí)需要引入移軸技術(shù),這可以借于Matlab等軟件。在后繼分析中需要剔除剛體模態(tài)的影響。

        3 求解流程

        如前所述,在進(jìn)行數(shù)值求解之前需要完成以下數(shù)據(jù)準(zhǔn)備工作:

        1) 確定分析的列車(chē)運(yùn)行速度以及時(shí)間步長(zhǎng);

        2)采用有限元分析軟件建立橋梁分析模型,提取結(jié)構(gòu)的頻率、振型,計(jì)算或假定模態(tài)阻尼比,計(jì)算各階的運(yùn)動(dòng)狀態(tài)遞推系數(shù);

        3)利用Matlab軟件的符號(hào)推導(dǎo)功能推導(dǎo)車(chē)輛運(yùn)動(dòng)方程,提取車(chē)輛頻率、振型、計(jì)算等效阻尼比,計(jì)算車(chē)輛的各階運(yùn)動(dòng)狀態(tài)遞推系數(shù);

        4)輸入系統(tǒng)激勵(lì)——軌道不平順;

        5)確定車(chē)輛、橋梁子系統(tǒng)的初始振動(dòng)狀態(tài)。

        為了便于程序編制,將式(15)展開(kāi)得到

        (31a)

        (31b)

        從式(31)可知,t+Δt時(shí)刻運(yùn)動(dòng)狀態(tài)由兩項(xiàng)組成:一項(xiàng)是t時(shí)刻運(yùn)動(dòng)狀態(tài)及外力的貢獻(xiàn);二項(xiàng)是t+Δt時(shí)刻外力的貢獻(xiàn)。

        式(27)、式(29)構(gòu)成了等效的車(chē)橋動(dòng)力相互作用系統(tǒng)運(yùn)動(dòng)方程組,可以應(yīng)用式(31)進(jìn)行逐步求解。在給定t時(shí)刻動(dòng)力相互作用系統(tǒng)狀態(tài)的條件下,迭代法計(jì)算t+Δt時(shí)刻系統(tǒng)響應(yīng)的步驟見(jiàn)圖2。其中,判別計(jì)算收斂的條件為在迭代過(guò)程中橋梁、車(chē)輛子系統(tǒng)的力、位移向量的范數(shù)變化趨于穩(wěn)定。

        ‖F(xiàn)i-Fi-1‖≤E·‖F(xiàn)i-1‖

        (32)

        ‖ui-ui-1‖≤E·‖ui-1‖

        (33)

        式中:上角標(biāo)i、i-1為時(shí)間步內(nèi)兩次臨近的迭代;E為給定的相對(duì)容差,通??梢匀?0-6。

        圖2 逐步迭代求解流程Fig.2 Step-by-step iteration solution procedure

        依據(jù)上述步驟編制分析程序。

        4 算例分析

        選取一節(jié)客車(chē)勻速通過(guò)簡(jiǎn)支梁作為分析對(duì)象,其中車(chē)輛參數(shù)、橋梁參數(shù)、輪軌關(guān)系、分析車(chē)速均與文獻(xiàn)[10]一致。為驗(yàn)證解析解法的有效性,增加了與基于Newmark-β法、高斯精細(xì)積分法迭代求解結(jié)果的對(duì)比。其中,計(jì)算得到車(chē)輛子系統(tǒng)的各階頻率及等效阻尼比見(jiàn)表1。

        表1 車(chē)輛各階頻率及等效阻尼比

        選取簡(jiǎn)支梁跨中豎向加速度、車(chē)輛第一輪對(duì)豎向加速度作為分析項(xiàng)目。圖3、圖4顯示了積分時(shí)間步長(zhǎng)Δt=10-4s時(shí)三種方法求解的結(jié)果。

        圖3 跨中豎向加速度時(shí)程(Δt=10-4 s)Fig.3 Acceleration time history of bridge at mid-span(Δt=10-4 s)

        圖4 輪對(duì)豎向加速度時(shí)程(Δt=10-4 s)Fig.4 Vertical acceleration time history of wheel-set(Δt=10-4 s)

        從對(duì)比情況可知:對(duì)橋梁振動(dòng)情況而言,三種方法計(jì)算基本一致;就輪對(duì)振動(dòng)情況而言,解析解法與Newmark-β法、高斯精細(xì)積分法略有差異,這是由于車(chē)輛引入振型疊加法所致。

        圖5、圖6顯示了在積分步長(zhǎng)增大到5倍,即Δt=5×10-4s條件下,三種方法求解的結(jié)果。從對(duì)比情況可以看出,在積分步長(zhǎng)較大時(shí)高斯精細(xì)積分和解析解法依然能夠得到與原積分步長(zhǎng)一致的結(jié)果,Newmark-β法計(jì)算結(jié)果已經(jīng)發(fā)散。

        圖5 橋梁跨中豎向加速度時(shí)程(Δt=5×10-4 s)Fig.5 Acceleration time history of bridge at mid-span (Δt=5×10-4 s)

        圖6 輪對(duì)豎向加速度時(shí)程(Δt=5×10-4 s)Fig.6 Acceleration time history of wheel-set (Δt=5×10-4 s)

        通過(guò)對(duì)比發(fā)現(xiàn),與采用Newmark-β法相比較,在積分步長(zhǎng)較大時(shí)高斯精細(xì)積分和解析解法依然能夠得到與原積分步長(zhǎng)一致的結(jié)果。

        5 結(jié) 論

        依據(jù)上述理論分析和算例計(jì)算結(jié)果,可以得出如下結(jié)論:

        (1)相對(duì)于Newmark-β法,解析解法和高斯精細(xì)積分法具有求解精度高的特點(diǎn),求解車(chē)橋耦合問(wèn)題時(shí)可以采用較大的積分步長(zhǎng)。

        (2)相對(duì)于高斯精細(xì)積分法,解析解法避免了指數(shù)矩陣的復(fù)雜運(yùn)算。

        本文提出的分析方法可以應(yīng)用到實(shí)際工程中。

        [ 1 ] XIA H, DE ROECK G, GOICOLEA J M. Bridge vibration and controls: new research [M]. New York: Nova Science

        Publishers, 2011.

        [ 2 ] 杜憲亭. 強(qiáng)地震作用下大跨度橋梁空間動(dòng)力效應(yīng)及列車(chē)運(yùn)行安全研究[D]. 北京:北京交通大學(xué), 2011.

        [ 3 ] 翟婉明, 夏禾. 列車(chē)-軌道-橋梁動(dòng)力相互作用理論與工程應(yīng)用[M].北京:科學(xué)出版社, 2011.

        [ 4 ] 喬宏, 夏禾, 杜憲亭. 基于Duhamel積分的車(chē)橋耦合動(dòng)力分析方法[J]. 西南交通大學(xué)學(xué)報(bào), 2014, 49(5): 766-771. QIAO Hong, XIA He, DU Xianting. Analytical method for calculating dynamic response of coupled train-bridge system based on Duhamel integral [J]. Journal of Southwest Jiaotong University, 2014, 49(5): 766-771.

        [ 5 ] 李小珍, 馬文彬, 強(qiáng)士中.車(chē)橋系統(tǒng)耦合振動(dòng)分析的數(shù)值解法[J]. 振動(dòng)與沖擊, 2002, 21(3): 21-25. LI Xiaozhen, MA Wenbin, QIANG Shizhong. Coupling vibration analysis of vehicle-bridge system by iterative solution method [J]. Journal of Vibration and Shock, 2002, 21(3): 21-25.

        [ 6 ] ZHAI W M. Two simple fast integration methods for large-scale dynamic problems in engineering [J]. International Journal for Numerical Methods in Engineering, 1996, 39(4): 4199-4214.

        [ 7 ] 翟婉明. 非線性結(jié)構(gòu)動(dòng)力分析的Newmark預(yù)測(cè)-校正積分模式[J]. 計(jì)算結(jié)構(gòu)力學(xué)及其應(yīng)用, 1990, 7(2): 51-58. ZHAI Wanming. The predictor-corrector scheme based on the Newmarks integration algorithm for nonlinear structural dynamic analysis [J]. Chinese Journal of Computational Mechanics, 1990, 7(2): 51-58.

        [ 8 ] YANG Y B, LIN B H. Vehicle-bridge interaction analysis by dynamic condensation method [J]. Structure Engineering ASCE, 1995, 121(11): 1636-1643.

        [ 9 ] 張純, 胡振東, 仲政. 車(chē)橋耦合振動(dòng)分析的Haar小波方法[J]. 振動(dòng)與沖擊, 2007, 26(4): 77-80. ZHANG Chun, HU Zhendong, ZHONG Zheng. Vibration analysis for vehicle-bridge intercation by Haar wavelet method [J]. Journal of Vibration and Shock, 2007, 26(4): 77-80.

        [10] 杜憲亭,夏禾, 李慧樂(lè),等. 基于改進(jìn)高斯精細(xì)積分法的車(chē)橋耦合振動(dòng)分析框架[J]. 工程力學(xué), 2013, 30(9): 171-176. DU Xianting, XIA He, LI Huile, et al. Dynamic analysis framework of train-bridge system based on improved Gauss precise integration method [J]. Engineering Mechanics, 2013, 30(9): 171-176.

        [11] CLOUGH R W, PENZIEN J. Dynamics of structures [M]. 3 ed. New York: Computers & Structures Inc., 1995.

        Analytical solution to a sdof system and its application in dynamic analysisof a train-bridge system

        DU Xianting1, QIAO Hong1, XIA He1, WANG Shaoqin2

        (1. School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China;2. School of Science, Beijing University of Civil Engineering and Architecture, Beijing 100044, China)

        With the assumption that a load varies linearly within each time step, an analytical solution was derived for a single-degree-of-freedom(SDOF) vibration system with low damping. Based on this solution, the modeling of a train-bridge system and its solving procedure were deduced. The train-bridge system model consisted of a train subsystem and a bridge one. The motion equation of the train subsystem was derived using the rigid component assumption and D′Alembert’s principle, and that of the bridge subsystem was derived using the finite element method. The mode superposition method was applied to uncouple the equations of motion of the two subsystems. The effects of the non-orthogonal damping of the train subsystem and the dynamic interaction between two subsystems were treated as nonlinear virtual forces. A 4-axle vehicle passing through a simply-supported beam of 32 m span at a constant speed was taken as a case study. The dynamic analysis of the coupled system was performed using the proposed analytical method, Newmark-βmethod and Gauss precise integration method, respectively. The results showed that compared with Newmark-βmethod and Gauss precise integration method, the analytical method can not only improve the time step of numerical integration but also avoid the computation of complex exponential matrices, so it has a good applicability in engineering.

        train-bridge system; dynamic interaction; analytical solution to a single-degree-of-freedom system; mode superposition method; non-orthogonal damping

        973計(jì)劃(2013CB036203);中央高校基本科研業(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金資助(2014JBM092)

        2015-10-27 修改稿收到日期:2016-02-14

        杜憲亭 男,博士,副教授, 1978年生

        U24

        A

        10.13465/j.cnki.jvs.2017.07.006

        猜你喜歡
        橋梁振動(dòng)
        振動(dòng)的思考
        噴水推進(jìn)高速艇尾部振動(dòng)響應(yīng)分析
        This “Singing Highway”plays music
        手拉手 共搭?lèi)?ài)的橋梁
        句子也需要橋梁
        振動(dòng)攪拌 震動(dòng)創(chuàng)新
        中立型Emden-Fowler微分方程的振動(dòng)性
        加固技術(shù)創(chuàng)新,為橋梁健康保駕護(hù)航
        無(wú)人機(jī)在橋梁檢測(cè)中的應(yīng)用
        高性能砼在橋梁中的應(yīng)用
        久久精品亚洲中文无东京热| 国产精品无码无卡无需播放器| 国产69精品久久久久777| 亚洲欧洲无码一区二区三区| 精品人妻丰满久久久a| 人妻秘书被社长浓厚接吻| 狠狠cao日日橹夜夜十橹| 干出白浆视频在线观看| 国产日产欧产精品精品蜜芽| 少妇做爰免费视频网站| 国产人澡人澡澡澡人碰视频| 92精品国产自产在线观看48页 | 日本一区二区三区在线播放| 久久91精品国产一区二区| 粗大的内捧猛烈进出看视频| 精品人妻无码一区二区色欲产成人| 亚洲国产精品无码久久九九大片健| 青青草视频在线观看精品在线| 亚洲av永久无码精品漫画| 国产免费午夜a无码v视频| 一级午夜视频| 在线看不卡的国产视频| 国产一区二区三区在线综合视频| 少妇性饥渴bbbbb搡bbbb| 国产91吞精一区二区三区| 国产女主播在线免费看| 国产极品粉嫩福利姬萌白酱| 免费国产黄网站在线观看可以下载| 免费观看国产精品| 精品亚洲视频免费观看网站| 午夜视频在线观看视频在线播放| 久久综合精品国产二区无码 | 加勒比东京热中文字幕| 天天综合网天天综合色| 国产精品视频一区日韩丝袜| 国产AV秘 无码一区二区三区| 丝袜美腿丝袜美腿丝袜美腿丝袜| 日韩人妻一区二区三区蜜桃视频| 99视频在线国产| 久久国产精品免费久久久| 深夜爽爽动态图无遮无挡|