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

        ?

        彈性邊界下圓弧拱的自由振動分析

        2016-12-15 10:41:25趙章泳邱艷宇王明洋宋春明
        振動與沖擊 2016年21期
        關(guān)鍵詞:無量圓弧振型

        趙章泳, 邱艷宇,2, 王明洋,2, 宋春明,2, 曹 侃

        (1.解放軍理工大學(xué) 爆炸沖擊防災(zāi)減災(zāi)國家重點實驗室,南京 210007; 2.南京理工大學(xué) 機(jī)械工程學(xué)院,南京 210094)

        ?

        彈性邊界下圓弧拱的自由振動分析

        趙章泳1, 邱艷宇1,2, 王明洋1,2, 宋春明1,2, 曹 侃1

        (1.解放軍理工大學(xué) 爆炸沖擊防災(zāi)減災(zāi)國家重點實驗室,南京 210007; 2.南京理工大學(xué) 機(jī)械工程學(xué)院,南京 210094)

        工程中圓弧拱的邊界不能總被簡化為理想的簡支或固支形式。為了研究彈性支承對圓弧拱自由振動特性的影響規(guī)律,將力、位移等變量無量綱化。根據(jù)平衡方程和坐標(biāo)轉(zhuǎn)換推導(dǎo)得出極坐標(biāo)下圓弧拱在水平、豎直和轉(zhuǎn)動方向支撐條件為彈性時的邊界條件方程。并采用考慮彎曲和軸向變形而忽略剪切變形及轉(zhuǎn)動慣量的自由振動的運動控制方程。運動方程在邊界條件情況下,其解僅為關(guān)于矢跨比f,細(xì)長比s和無量綱剛度陣[K]的函數(shù)。采用Runge-Kutta法和行列式搜索法求解運動方程的特征值即無量綱頻率Ωi以及特征向量即振型。通過計算發(fā)現(xiàn),與理想支撐相比,彈性支承情況下細(xì)長比s對拱自振頻率的影響要明顯下降。理想支撐情況下圓弧拱的自振頻率越高,則彈性支承對其自振頻率的影響越小。與水平和豎向彈性支承相比,轉(zhuǎn)動方向彈性支承僅對圓弧拱基頻有較大影響。

        圓弧拱;自由振動;彈性邊界;振型;自振頻率;無量綱化

        同梁、板、柱等一樣,拱結(jié)構(gòu)是工程中最重要的基本構(gòu)件之一,在不同領(lǐng)域都有著廣泛的應(yīng)用,CHIDAMPARAM 等[1]整理了一百多年來關(guān)于拱的研究。由于工程中拱結(jié)構(gòu)的橫截面通常都是左右對稱的,因此可以將拱的平面內(nèi)撓曲振動和平面外扭轉(zhuǎn)振動解耦從而分別處理[2]。HENRYCH[3]通過嚴(yán)密的數(shù)學(xué)推導(dǎo)得出了極坐標(biāo)下拱結(jié)構(gòu)的考慮彎曲、軸向、剪切變形及轉(zhuǎn)動慣量等因素時的數(shù)學(xué)模型,并且得到了忽略軸向變形和切向慣性力的情況下拱振型的閉合解。然而對于較淺的拱,其軸向變形對拱自由振動頻率具有顯著影響[4]。KARNOVSKY[5]著重對在一定程度上滿足工程要求的基礎(chǔ)上對拱結(jié)構(gòu)進(jìn)行多自由度離散。在這些研究中,絕大多數(shù)采用的都是理想的簡支和固支的邊界條件。但是在工程中結(jié)構(gòu)的邊界不能總簡化為理想邊界,因此需要研究彈性支撐邊界對拱結(jié)構(gòu)自振特性的影響。LEE等[6]曾經(jīng)對拱結(jié)構(gòu)在徑向和轉(zhuǎn)動方向彈性的邊界情況進(jìn)行了研究。然而實際工程中的彈性卻總是存在于水平、豎直方向和轉(zhuǎn)動方向的??垫玫萚7]曾利用結(jié)構(gòu)動力學(xué)中柔度系數(shù)的概念對水平彈性支撐下園拱的動力特性進(jìn)行了計算。本文將更全面的針對圓弧拱在這三種彈性邊界下的自由振動特性進(jìn)行分析。

        1 計算模型

        α=4arctan(2f)

        (1)

        式中:f和s分別為拱的矢跨比和拱的細(xì)長比。拱的邊界在三個方向上分別具有大小為Kx、Ky、KΘ的彈性支撐。極坐標(biāo)系下拱軸線上θ位置點的徑向和切向位移分別為u、v,轉(zhuǎn)角為ζ。直角坐標(biāo)系下拱軸線上θ位置點的橫向和豎向位移分別為X、Y,轉(zhuǎn)角為Θ;對應(yīng)的力為Fx、Fy、M。

        圖1 圓弧拱幾何模型Fig.1 Circular arch model

        圖2所示為極坐標(biāo)系下拱微元體受外荷載作用的內(nèi)力圖,其中N、Q、M分別為軸力、剪力、彎矩,Pr、Pθ分別為微元體所受的徑向分布荷載和切向分布荷載。所有位移和力的正方向均同圖1、圖2中所示一致。

        圖2 拱微元體受力模型Fig.2 Loads on an arch element

        WUNG[8]曾對軸向變形、剪切變形及轉(zhuǎn)動慣量在理想支座下對拱的自振頻率的影響進(jìn)行了詳細(xì)的研究和計算,表明在細(xì)長比s小于0.1的拱中,軸向變形對拱的頻率影響遠(yuǎn)大于其他兩項。因此本文亦采用忽略剪切變形和轉(zhuǎn)動慣量的影響,考慮軸向變形的幾何方程。HENRYCH[3]推導(dǎo)得到了這種情況下拱的幾何、物理和平衡方程分別如式(2)~(4)所示:

        (2)

        (3)

        (4)

        需要注意的是,由于忽略了剪切變形,剪力Q是由平衡方程導(dǎo)出而并非物理方程。

        為了避免數(shù)值計算中可能出現(xiàn)的病態(tài)現(xiàn)象以及計算結(jié)果的推廣,將所有物理量無量綱化。定義無量綱位移為位移與l之比,無量綱彎矩為彎矩與EI/l之比,無量綱力為力與EI/l2之比,轉(zhuǎn)動彈簧無量綱剛度為剛度與π4EI/l之比,線性彈簧無量綱剛度為剛度與π4EI/l3之比[6]。為了方便推導(dǎo)和計算,采用向量形式表示變量:

        (5)

        向量中元素均已無量綱化,在后文中如無特殊說明,所有變量均為無量綱變量。

        (6)

        圓拱內(nèi)力與位移的關(guān)系為:

        (7)

        式中:

        如果考慮圓拱的軸向變形而忽略剪切變形和轉(zhuǎn)動慣量,則圓拱自由振動的控制微分方程為[3]:

        (8)

        (9)

        在邊界處,θ=±α/2,根據(jù)平衡條件有:

        (10)

        對于彈簧,根據(jù)胡克定律:

        (11)

        又直角坐標(biāo)和極坐標(biāo)的轉(zhuǎn)換關(guān)系為:

        (12)

        由于[T]陣為對稱正交陣,聯(lián)立式(10)~(12)可得,在邊界處有:

        (13)

        式中:

        聯(lián)立式(7)、(13)可得:

        (14)

        式(14)即為彈性邊界下圓拱自由振動的邊界條件。

        2 振型正交性

        在振型正交性的推導(dǎo)中,由于對如同時間,質(zhì)量等變量進(jìn)行無量綱化并不影響結(jié)果,因而在振型正交性推導(dǎo)中并不采用無量綱參數(shù)。

        記:

        [dm(θ,t)]=am[dm(θ)]sin(ωmt+φm)

        [dn(θ,t)]=an[dn(θ)]sin(ωnt+φn)

        (15)

        分別為圓拱按第m、n階振型振動時的位移函數(shù)向量。此時第m階振型的慣性力向量為

        [Im(θ,t)]=

        (16)

        其中質(zhì)量陣

        (17)

        則在任意時刻t,第m階振型對應(yīng)的慣性力按照第n階振型的虛位移上所做的虛功為

        sin(ωmt+φm)sin(ωnt+φn)

        (18)

        同理在任意時刻t,第n階振型對應(yīng)的慣性力按照第m階振型的虛位移上所做的虛功為:

        sin(ωmt+φm)sin(ωnt+φn)

        (19)

        根據(jù)Maxwell-Betti定律,第一組荷載在第二組荷載所引起的變位上所做的功等于第二組荷載在第一組荷載所引起的變位上所做的功,則

        Wm,n(t)=Wn,m(t)

        (20)

        由于

        [dm(θ)]T[M]T[dn(θ)]=

        [dn(θ)]T[M]T[dm(θ)]

        (21)

        聯(lián)立式(18)、(19)、(20)可得

        其中

        H(t)=amansin(ωmt+φm)sin(ωnt+φn)

        (22)

        由于對任意時刻t,式(22)成立,且對一般結(jié)構(gòu),圓頻率為重根的情況非常少見[9]。故:

        (m≠n)

        (23)

        此式即為圓拱結(jié)構(gòu)對質(zhì)量的正交性。

        3 數(shù)值方法的收斂性及與有限元對比

        為了解決微分方程(8)在邊界條件(14)下的解,可以使用AL-KHAFAJI[10]以及LEE等[11]介紹的數(shù)值方法。通過行列式尋根法和試位法來計算拱的無量綱自振頻率Ωi,再根據(jù)Runge-Kutta計算拱的振型,其中i表示模態(tài)階數(shù)。

        圖3 數(shù)值方法收斂性分析Fig.3 Convergence analysis of the numerical method

        圖3為f=0.3,s=0.01前四階頻率隨半拱劃分段數(shù)N的收斂情況。從圖中可發(fā)現(xiàn):邊界彈簧剛度越小,收斂越快;振型階數(shù)越低,收斂越快;當(dāng)N>40時,前四階的自振頻率Ωi均能收斂。通過對不同的f、s及[K]的計算,N=40的劃分均能收斂。所以在計算中,統(tǒng)一采用了N=40的劃分方法。

        在對拱的有限元計算中,常用線單元對拱進(jìn)行離散。采用多少線單元可以得到理想的近似效果是值得注意的。通過對離散后圓拱的剛度陣和質(zhì)量陣進(jìn)行無量綱化,同樣可以得到其無量綱頻率Ωi[4]。將全拱劃分為100段,計算得到了幾種典型情況下圓拱的無量綱頻率Ωi?,F(xiàn)將兩種方法計算結(jié)果列在表1中對比。通過對比發(fā)現(xiàn),當(dāng)劃分段數(shù)達(dá)到全拱100段后,對于不同矢跨比和細(xì)長比,前四階振型均能滿足2%以內(nèi)的誤差。

        表1 本文Ωi的計算結(jié)果與有限元方法對比

        4 討 論

        4.1 彈性支撐下矢跨比對頻率的影響

        對于板拱,當(dāng)f在0.08附近時,存在著振型的轉(zhuǎn)換點,即當(dāng)f減小時,拱的一階振型會從反對稱變換成為對稱。這種變化意味著曲桿作為拱結(jié)構(gòu)的剛度的喪失,并轉(zhuǎn)化為曲梁[2]。通過本文計算,發(fā)現(xiàn)當(dāng)拱具有彈性支撐時,這種現(xiàn)象將不僅僅出現(xiàn)在淺拱中。如圖4所示為s=0.01,kx=∞、ky=100、kΘ=0時圓拱前四階自振頻率Ωi隨矢跨比f的變化曲線,圖中ia、is分別為反對振動和對稱振動的階數(shù)??梢园l(fā)現(xiàn)第一階對稱振動和第一階反對稱結(jié)構(gòu)自振頻率曲線在f=0.4時相交,之后第一振型由對稱變?yōu)榉磳ΨQ,而在理想邊界條件下,這一現(xiàn)象發(fā)生在f=0.08附近。且第二階對稱振動和第二階反對稱振動的頻率曲線在f=0.07處發(fā)生了相交。而在LEE[6]的算例中,高階的頻率曲線只是靠近而并沒有相交。說明對于高階振型,彈性邊界也會使得對稱與反振型發(fā)生互換。

        圖4 Ωi 隨矢跨比f變化曲線Fig.4 The ΩiVersus rise to length ratio f curves

        由于這種振型的互換現(xiàn)象,為了研究彈性邊界對拱自振頻率的影響規(guī)律,在計算研究中把反對稱振動和對稱振動區(qū)分開來是很有必要的。

        4.2 彈性支撐下細(xì)長比對拱頻率的影響

        圓拱的細(xì)長比s對于拱的頻率變化的影響是十分顯著的[12]??傮w來說,矢跨比越小,模態(tài)階數(shù)越高,影響越為明顯[13]。這是由于s越大,拱軸線的軸向應(yīng)變與拱截面上遠(yuǎn)離拱軸線的纖維因彎曲而發(fā)生的應(yīng)變相比,越來越不能忽略[14]。拱的軸向應(yīng)變在拱振動過程中的影響主要表現(xiàn)在使圓拱的無量綱自振頻率下降[8]。如圖5所示,在豎向彈性支撐條件下,由于圓拱邊界剛度的下降導(dǎo)致拱振動過程中軸力減小,拱軸向變形也相應(yīng)減小,因而較兩端理想支撐情況下,s對Ωi的影響也大大減小。與豎向彈性支承相比,當(dāng)水平彈性支承剛度減小時,拱中軸力將減小而彎矩將增大[15],因而s對Ωi的影響將比豎向彈性支撐更弱,在計算中也證明了這一點。

        圖5 Ωi 隨細(xì)長比s變化曲線Fig.5 The ΩiVersus slenderness s curves

        4.3 彈性支承對拱結(jié)構(gòu)振型與頻率影響的關(guān)系

        為了更深入的研究彈性邊界如何使得對稱和反對稱振動的頻率發(fā)生改變,將f=0.2,s=0.001的彈性支撐圓拱與理想邊界圓拱的前四階振型繪制于圖6中。從中可以發(fā)現(xiàn),對于一階振型,由于其模態(tài)形狀并未發(fā)生改變,因此其頻率也并無太大改變。而對其余三階振型,彈性支撐使得其振型中的“半波”數(shù)目減少,即其彎曲形式發(fā)生明顯變化,因此其頻率將發(fā)生顯著下降。這是因為根據(jù)Rayleigh原理,自振頻率的平方等于勢能與動能之比[16]。如果彈性支撐剛度的增加并未使得其模態(tài)振型發(fā)生顯著變化,則其將僅通過略微減小動能而引起自振頻率的增加。當(dāng)剛度繼續(xù)增加時,振型將會發(fā)生顯著變化(主要表現(xiàn)為半波數(shù)目的增加),此時由于勢能的增加將會使得自振頻率明顯升高。

        圖6 豎向彈性支撐對振型影響Fig.6 The influence of the general boundary onto mode shapes

        4.4 水平彈性支承對拱結(jié)構(gòu)的頻率影響

        圖7 水平方向彈性支承kx對頻率比Ωi/Ω′i的影響Fig.7 Effects of kxon Ωi/Ω′i of horizontal elastic boundary conditions

        4.5 豎向彈性支承對拱結(jié)構(gòu)的頻率影響

        圖8為f=0.2,s=0.01的雙鉸圓拱,豎直方向無量綱剛度ky對頻率比Ωi/Ω′i的影響。其形狀與規(guī)律同圖7大致相同。但是對頻率比Ωi/Ω′i有影響的ky值的范圍將明顯減小。

        圖8 豎直方向彈性支承ky對頻率比Ωi/Ω′i的影響Fig.8 Effects of kyon Ωi/Ω′i of vertical elastic boundary conditon

        為了進(jìn)一步確定這種影響現(xiàn)象,通過計算矢跨比f從0.2~0.5,細(xì)長比s從0.1~0.001的圓拱的頻率恢復(fù)曲線,其規(guī)律同圖7、圖8所示相同。在此范圍內(nèi),彈性支座對于圓拱頻率的影響隨矢跨比f和細(xì)長比s的減小而增加,隨頻率階數(shù)的升高而增加。這兩種影響都可以歸結(jié)為:彈性支座對圓拱頻率的影響的增加隨圓拱在理想支座情況下頻率的升高而降低。

        4.6 轉(zhuǎn)動方向彈性支撐對拱結(jié)構(gòu)的頻率影響

        圖 9所示為轉(zhuǎn)動方向彈性支撐對圓拱結(jié)構(gòu)的自振頻率的影響。kΘ由0到無窮,圓拱即由雙鉸拱過度到無鉸拱。從中可以看出,與水平彈性支承與豎向彈性支承相比,轉(zhuǎn)動方向彈性支撐的剛度只在很小的范圍內(nèi)影響拱結(jié)構(gòu)的頻率,且模態(tài)階數(shù)越高,轉(zhuǎn)動支撐的剛度對頻率影響越弱。

        圖9 轉(zhuǎn)動方向彈性支承kΘ對頻率比Ωi/Ω′i的影響Fig.8 Effects of kΘon Ωi/Ω′i of rotational elastic boundary conditon

        5 結(jié) 論

        本文通過將圓拱參數(shù)無量綱化,建立圓拱邊界為豎直、水平和轉(zhuǎn)動彈性支撐下的無量綱微分方程組和邊界條件。并通過數(shù)值方法求解了圓拱的頻率和振型隨彈性支撐的變化所發(fā)生的變化。

        (1) 在彈性支撐條件下,拱的基頻由對稱轉(zhuǎn)換為非對稱的臨界矢跨比將會增大。

        (2) 與理想邊界情況相比,在彈性邊界下,細(xì)長比s對Ωi的影響將大大減小。

        (3) 對于不同矢跨比f、細(xì)長比s的拱的某一階模態(tài),水平和豎向彈性支承均只在某一范圍內(nèi)顯著影響拱的Ωi。并且彈性支承對Ωi的影響范圍和影響程度隨矢跨比f和細(xì)長比s的減小、模態(tài)階數(shù)的增大而減小。彈性支承的變化對Ωi的顯著影響主要是使圓弧拱的振型發(fā)生了變化。

        在忽略轉(zhuǎn)動慣量的情況下,轉(zhuǎn)動彈性支承僅在很小范圍內(nèi)對基頻有顯著影響,對高階Ωi的影響較小。

        [1] CHIDAMPARAM P, LEISSA A W. Vibrations of planar curved beams, rings, and arches[J]. Applied Mechanics Reviews, 1993, 46(9): 467-483.

        [2] 項海帆. 拱結(jié)構(gòu)的穩(wěn)定與振動/高等結(jié)構(gòu)力學(xué)叢書[M]. 北京:人民交通出版社, 1991.

        [3] HENRYCH J. The dynamics of arches and frames[M]. Elsevier Scientific Pub. Co., 1981.

        [4] 錢七虎,孫乃光, 陳震元,等. 拱型結(jié)構(gòu)的自振頻率計算及軸向變形對自振頻率的影響[D]. 錢七虎院士論文選集, 1999.

        [5] KARNOVSKY I A. Theory of arched structures[M]. New York, Ny: Springer New York, 2012.

        [6] LEE B K, LEE J Y, CHOI K M, et al. Free vibrations of arches with general boundary condition[J]. KSCE Journal of Civil Engineering, 2002, 6(4): 469-474.

        [7] 康婷, 白應(yīng)生, 孫惠香. 水平彈性支撐圓拱的動力特性研究[J]. 力學(xué)與實踐, 2013, 35(2): 50-55. KANG TING,BAI YING SHENG ,SUNHUIXIANG.Dynamic characteristics of the horizontal elastic support circular arch[J].Mechanics in Engineering,2013,35(2):50-55.

        [8] WUNG S. Vibration of hinged circular arches[D]. Rice University, 1967.

        [9] 張相,黃東才. 結(jié)構(gòu)振動力學(xué)[M]. 上海:同濟(jì)大學(xué)出版社, 1994.

        [10] AL-KHAFAJI A W, TOOLEY J R. Numerical methods in engineering practice[M]. Holt, Rinehart And Winston, 1986.

        [11] LEE B K, WILSON J F. Free vibrations of arches with variable curvature[J]. Journal of Sound and Vibration, 1990, 136(1): 75-89.

        [12] LEE B K, LEE T E, AHN D S. Free vibrations of arches with inclusion of axial extension, shear deformation and rotatory inertia in cartesian coordinates[J]. KSCE Journal of Civil Engineering, 2004, 8(1): 43-48.

        [13] CHIDAMPARAM P, LEISSA A W. Influence of centerline extensibility on the in-plane free vibrations of loaded circular arches[J]. Journal of Sound and Vibration, 1995, 183(5): 779-795.

        [14] 金尼克. 拱的穩(wěn)定性[M]. 上海:建筑工程出版社, 1958.

        [15] FANNING P J, BOOTHBY T E, ROBERTS B J. Longitudinal and transverse effects in masonry arch assessment[J]. Construction and Building Materials, 2001, 15(1): 51-60.

        [16] 克拉夫, 彭津. 結(jié)構(gòu)動力學(xué)[M]. 北京:高等教育出版社, 2006.

        Free vibration analysis of arches under elastic support boundary conditions

        ZHAO Zhangyong1, QIU Yanyu1,2, WANG Mingyang1,2, SONG Chunming1,2, CAO Kan1

        (1. State Key Laboratory of Disaster Prevention & Mitigation of Explosion & Impact, PLA University of Science and Technology, Nanjing 210007, China;2. School of Mechanical Engineering, Nanjing University of Science and Technology, Nanjing 210094, China)

        The boundary conditions of circular arches in practical engineering can not be always simplified as ideally hinged ends or fixed ones. In order to study the influences of elastic supports on the free vibration characteristics of circular arches, the variables, such as, forces and displacements were nondimensionalized. According to the balance equations and coordinate conversion, the boundary condition equations of arches with vertical, horizontal and rotational elastic supports were derived in polar coordinates. The flexural and axial deformations were considered while shear deformations and rotation inertias were neglected to establish the governing equations of motion of arches. The solutions to the governing equations of motion under the boundary conditions were the functions with respect to rise to span ratiof, slenderness s, and dimensionless stiffness matrix [K]. Runge-Kutta method and the determinant search method were used to solve the eigenvaluesΩiof the equations of motion and eigenvectors, i.e., modal shapes. Through calculating different configurations, it was shown that under the conditions of elastic supports, the effects of slenderness s on the arches’ natural frequencies are much smaller than those under the conditions of ideal supports; the effects of horizontal and vertical elastic supports on the natural frequencies decrease with increase in natural frequencies of arches under the conditions of ideal supports; the rotational elastic supports obviously influence the fundamental frequency of arches.

        circular arch; free vibration; elastic support; modal shape; natural frequency; non-dimensionalization

        2015-04-21 修改稿收到日期:2015-10-08

        趙章泳 男,碩士生,1992年生

        宋春明 男,博士生,副教授,1979年生

        E-mail:ming1979@126.com

        O327

        A

        10.13465/j.cnki.jvs.2016.21.018

        猜你喜歡
        無量圓弧振型
        烏雷:無量之物
        關(guān)于模態(tài)綜合法的注記
        縱向激勵下大跨鋼桁拱橋高階振型效應(yīng)分析
        淺析圓弧段高大模板支撐體系設(shè)計與應(yīng)用
        劉少白
        藝術(shù)品(2020年8期)2020-10-29 02:50:02
        外圓弧面銑削刀具
        塔腿加過渡段輸電塔動力特性分析
        論書絕句·評謝無量(1884—1964)
        炳靈寺第70 窟無量壽經(jīng)變辨識
        西藏研究(2017年3期)2017-09-05 09:45:07
        結(jié)構(gòu)振型幾何辨識及應(yīng)用研究
        山西建筑(2015年14期)2015-06-05 09:37:07
        亚洲一区二区视频免费看| 无码人妻丰满熟妇啪啪7774| 91视频免费国产成人| 素人系列免费在线观看| 亚洲精品第四页中文字幕| 超碰色偷偷男人的天堂| 无码精品a∨在线观看十八禁| 久热爱精品视频在线观看久爱 | 91日本精品国产免| 日韩精品人妻中文字幕有码| 国产亚洲熟妇在线视频| 奇米影视第四色首页| 黄色资源在线观看| 亚洲中文字幕在线第二页| 亚洲日本精品国产一区二区三区| 韩日午夜在线资源一区二区| 亚洲精品视频久久| 久久综合这里只有精品| 国产黑丝美腿在线观看| 国产成年无码v片在线| 日韩啪啪精品一区二区亚洲av| av国产免费在线播放| 久久人妻av一区二区软件| 日本亚洲国产一区二区三区| 国产精品麻豆A在线播放| 国产性虐视频在线观看| 又黄又硬又湿又刺激视频免费| 国产片AV在线永久免费观看| 国产精品国产三级国a| 久久精品中文字幕无码绿巨人| 精品国产18久久久久久| 无码制服丝袜中文字幕| 伊人青青草综合在线视频免费播放 | 亚洲美女主播内射在线| 国产做爰又粗又大又爽动漫| 91av精品视频| 久久精品成人一区二区三区蜜臀| 日韩 无码 偷拍 中文字幕| 亚洲国产美女精品久久久久| 一区二区三区午夜视频在线观看| 一区二区三区国产精品乱码|