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

        ?

        深海頂張式立管組合參激共振的非線性振動(dòng)分析

        2015-06-05 15:30:53唐友剛朱龍歡李楊青
        關(guān)鍵詞:渦激浮體立管

        唐友剛,朱龍歡,李楊青

        (1. 天津大學(xué)水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室,天津 300072;2. 天津大學(xué)建筑工程學(xué)院,天津 300072)

        深海頂張式立管組合參激共振的非線性振動(dòng)分析

        唐友剛1,2,朱龍歡1,2,李楊青1,2

        (1. 天津大學(xué)水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室,天津 300072;2. 天津大學(xué)建筑工程學(xué)院,天津 300072)

        研究了深海頂張式立管參數(shù)激勵(lì)和渦激共同作用下的非線性振動(dòng)特性. 考慮平臺(tái)升沉運(yùn)動(dòng)激勵(lì)和渦激力建立立管振動(dòng)方程,采用多尺度方法求解立管振動(dòng)方程的近似解析解. 考慮和型組合參激共振Ω≈ω1+ω2情況研究立管的振動(dòng)特性,計(jì)算得到了立管的幅頻響應(yīng)曲線,分析了平臺(tái)升沉運(yùn)動(dòng)對(duì)深海立管非線性振動(dòng)的影響. 結(jié)果表明:當(dāng)參激頻率滿足和型組合參激共振條件時(shí),立管振動(dòng)響應(yīng)中頻率為1/2參激頻率的亞諧波成分明顯;且由于內(nèi)共振關(guān)系的存在,立管1階模態(tài)被激發(fā),其幅值遠(yuǎn)大于2階模態(tài)幅值;隨著平臺(tái)升沉運(yùn)動(dòng)幅值的增大,立管橫向振動(dòng)幅值顯著增大,這表明平臺(tái)運(yùn)動(dòng)對(duì)于立管彎曲振動(dòng)有重要影響.

        頂張式立管(TTRs);組合參激共振;多尺度法;渦激振動(dòng)

        立管是海洋油氣開發(fā)中連接海底井口和上部浮體的結(jié)構(gòu)[1].當(dāng)水深達(dá)到1,500,m以上時(shí),海洋環(huán)境尤為惡劣.立管頂部承受平臺(tái)升沉作用,內(nèi)部有高壓的油氣通過,外部有復(fù)雜的浪、流等環(huán)境荷載作用,此環(huán)境下立管會(huì)發(fā)生碰撞、波頻振動(dòng)、渦激振動(dòng)、參激振動(dòng),這些都能引起立管的疲勞損傷和斷裂失效,不僅給工程本身帶來巨大經(jīng)濟(jì)損失,而且還會(huì)對(duì)自然環(huán)境造成嚴(yán)重的次生災(zāi)害[2].

        關(guān)于深海細(xì)長結(jié)構(gòu)的參激振動(dòng)響應(yīng)問題,國內(nèi)外學(xué)者進(jìn)行了大量理論和實(shí)驗(yàn)研究.Patel等[3-4]研究了低軸力條件下張力腿的不穩(wěn)定區(qū)域,用半解析的方法討論了張力腿在復(fù)合激勵(lì)下的動(dòng)力響應(yīng)問題;Park等[5]用有限元方法研究了海洋細(xì)長結(jié)構(gòu)物在參數(shù)激勵(lì)和波流激勵(lì)聯(lián)合作用下的動(dòng)力響應(yīng),并計(jì)算得到了馬蒂厄不穩(wěn)定參數(shù)區(qū)域;Chatjigeorgiou[6]針對(duì)參數(shù)激勵(lì)頻率接近立管非線性共振范圍的情況,研究參數(shù)激勵(lì)對(duì)立管振動(dòng)的影響,還討論了阻尼對(duì)于立管參數(shù)激勵(lì)振動(dòng)穩(wěn)定性的影響;Chatjigeorgiou等[7]在考慮內(nèi)流和參數(shù)激勵(lì)的情況下,針對(duì)參數(shù)激勵(lì)頻率接近兩倍1階固有頻率的特殊情況,得到立管振動(dòng)響應(yīng)的近似解析解.國內(nèi)學(xué)者朱克強(qiáng)[8]最早研究了考慮船舶的升沉和橫蕩運(yùn)動(dòng)情況時(shí)立管的振動(dòng)響應(yīng),用有限差分法求解了振動(dòng)方程,得到了立管在參數(shù)激勵(lì)和莫里森波浪力共同作用下的時(shí)間歷程響應(yīng);楊和振等[9]忽略波流力,考慮平臺(tái)浮體簡諧升沉運(yùn)動(dòng)作為參數(shù)激勵(lì),研究深海立管參數(shù)激勵(lì)振動(dòng)響應(yīng);邵衛(wèi)東等[10]考慮立管頂端動(dòng)邊界條件及立管軸力隨水深的變化,給出了考慮浮體升沉及張緊環(huán)運(yùn)動(dòng)時(shí)頂張力管的固有特征.

        本文將立管上部浮體的升沉運(yùn)動(dòng)簡化為參數(shù)激勵(lì),考慮立管大長徑比引起的幾何非線性,綜合渦激力的影響,將立管振動(dòng)處理為參激-強(qiáng)迫激勵(lì)振動(dòng)系統(tǒng),建立深海立管振動(dòng)方程,采用多尺度方法求立管振動(dòng)的解析解,得到立管升沉運(yùn)動(dòng)激勵(lì)和渦激共同作用下非線性振動(dòng)的響應(yīng)特點(diǎn),并研究浮體升沉幅值對(duì)立管振動(dòng)的影響.

        1 立管分析模型

        立管為大長細(xì)比結(jié)構(gòu),內(nèi)有不可壓縮無旋流體經(jīng)過,管外承受的載荷包括順流向拖曳力、垂直于流速的渦激力以及上部浮體升沉運(yùn)動(dòng)誘發(fā)的軸向力.考慮底部約束為鉸支,立管坐標(biāo)系原點(diǎn)位于底端,軸豎直向上,建立參數(shù)激勵(lì)下深海立管模型,如圖1所示.

        圖1 參數(shù)激勵(lì)下深海立管模型示意Fig.1 Schematic diagram of deep-water riser under parametric excitation

        圖1 中,L為立管總長度,v為均勻流速,上部浮體運(yùn)動(dòng)為簡諧運(yùn)動(dòng),其位移為ua( t)=u0cos(Ωt),u0為位移幅值,上部浮體的運(yùn)動(dòng)會(huì)引起立管上端部的簡諧升沉運(yùn)動(dòng).

        單位長度流體和立管的受力情況如圖2和圖3所示.

        圖2 單位長度流體受力Fig.2 Forces acting on the unit length of fluid

        圖3 單位長度立管受力Fig.3 Forces acting on the unit length of riser

        考慮單位長度流體和立管法向上受力的平衡條件[11-12],可得

        式中:m、ma和ml分別為單位長度的立管質(zhì)量、附連水質(zhì)量和流體質(zhì)量;vl為立管管內(nèi)流體的速度;w為立管法向的位移;m0=(m+ma+m1) g ;S為剪切力;ε為局部應(yīng)變;χ為彎曲曲率;φ為立管切線方向與水平方向的夾角;EI、EA分別為立管彎曲剛度和軸向剛度;F為作用在立管單位長度上的水動(dòng)力.F包括由漩渦泄放產(chǎn)生的渦激升力fL和立管在渦激力方向運(yùn)動(dòng)產(chǎn)生的拖曳力fD,即

        式中:wρ為周圍海水的密度;D為立管外徑;sω為渦激頻率;LC、DC分別為升力系數(shù)和流體阻尼系數(shù).

        僅考慮立管在垂直來流平面的運(yùn)動(dòng),立管局部應(yīng)變?chǔ)欧譃殪o應(yīng)變?chǔ)?和動(dòng)應(yīng)變?chǔ)?兩部分,且任意點(diǎn)的切向位移u( z, t)為法向位移的高階小量O( w2),即

        此外,曲率、剪力、應(yīng)變及轉(zhuǎn)角之間有如下關(guān)系[7]:

        同時(shí),結(jié)構(gòu)所受的軸向張力eT分為靜頂張力e0T和動(dòng)張力e1T兩部分,其中靜頂張力e0T為定常值,且張力與應(yīng)變呈線性關(guān)系,即

        將式(2)~(11)代入式(1)中并進(jìn)行化簡,得到立管法向的運(yùn)動(dòng)微分方程

        應(yīng)用伽遼金法和振型疊加原理對(duì)式(12)進(jìn)行處理,設(shè)解的形式為

        式中:αu(z )、αw(z)滿足立管邊界條件,且僅與位置坐標(biāo)z有關(guān);gi( z)、fi( z)分別為立管軸向振動(dòng)和彎曲振動(dòng)的振型函數(shù);ua( t)、wa( t)分別為施加于立管頂端的軸向和橫向位移邊界條件,即ua( t)= u0cos(Ωt),wa( t)=0.將式(13)和(14)代入式(12),考慮結(jié)構(gòu)阻尼的作用,在等式左邊添加2μiq˙項(xiàng),其中μi為結(jié)構(gòu)阻尼系數(shù),結(jié)合振型的正交性進(jìn)行化簡,得到關(guān)于qi( t)的微分方程為

        式中iω為立管的第i階固有頻率.本文用解析方法求解式(15),僅考慮上部浮體的升沉運(yùn)動(dòng)和立管的橫向振動(dòng),因此可忽略式中含a()w t的乘積項(xiàng)和jkp q項(xiàng).由于阻尼項(xiàng)的存在,高階模態(tài)迅速收斂,因此僅保留前兩階模態(tài)進(jìn)行分析,由此整理式(15)得

        其中

        2 多尺度法求解模態(tài)振動(dòng)微分方程

        考慮整個(gè)振動(dòng)系統(tǒng)為弱非線性系統(tǒng),可將方程(16)和(17)中各非線性項(xiàng)系數(shù)看作小量,運(yùn)用多尺度法求解模態(tài)振動(dòng)方程.其中渦激升力項(xiàng)為周期變化的外激勵(lì),本文考慮渦激升力頻率遠(yuǎn)離立管固有振動(dòng)頻率,即外激勵(lì)為非共振的情形.

        式中e為無量綱小參數(shù).設(shè)方程(16)和(17)的解的形式為

        將式(27)代入式(16)和(17),并考慮式(26),令e同冪次項(xiàng)系數(shù)相等,并統(tǒng)一省略各變量中的“^”,得到式(28)~(31),即

        式中D0、D1分別是對(duì)T0、T1求偏導(dǎo)數(shù)的運(yùn)算符號(hào). 系數(shù)表達(dá)式分別為

        設(shè)式(28)和(29)解的形式為

        式中:cc代表前面表達(dá)式的共軛;參數(shù)iΛ的表達(dá)式為

        將式(35)和(36)代入式(30)和(31),消除永年項(xiàng),得到關(guān)于11()A T和22()A T的方程組.此振動(dòng)系統(tǒng)中存在幾種參激共振的情形,包括12Ωω≈,22Ωω≈,.這里Ω為參激頻率,1ω和 2ω分別為橫向彎曲振動(dòng)1階和2階固有頻率.本文考慮組合參激共振情形,對(duì)12+Ωωω≈的情況進(jìn)行詳細(xì)討論.同時(shí),該立管模型兩端為鉸支邊界條件,前兩階固有頻率之間存在2倍關(guān)系,即前兩階模態(tài)間存在內(nèi)共振關(guān)系,由此可以得到可解性條件為

        代入式(38)和(39)中,分離虛部與實(shí)部,為求振動(dòng)的穩(wěn)態(tài)響應(yīng),令10a′=,即可得到前兩階近似解的表達(dá)形式.

        3 算例及結(jié)果

        本文算例以深海中1,500,m長的生產(chǎn)立管為例,計(jì)算和分析立管的振動(dòng)特性.對(duì)上文討論的解析近似解進(jìn)行編程,得到立管參數(shù)激勵(lì)下非線性振動(dòng)的幅頻響應(yīng)曲線,運(yùn)用數(shù)值解法求解響應(yīng)的時(shí)間歷程曲線,并做了對(duì)比分析.立管參數(shù)如表1所示.

        表1 深海1,500,m立管參數(shù)Tab.1 Parameters of 1,500,m deep-water riser

        將表1中參數(shù)代入式(38)和(39)中求解,得到組合參數(shù)激勵(lì)下立管振動(dòng)幅值隨解諧參數(shù)σ的變化情況,即幅頻響應(yīng)曲線,如圖4和圖5所示.

        組合參激共振情況下,由圖4和圖5分析可得,在σ趨近零的范圍內(nèi),兩階振幅均隨頻率增大而減小,且1階振幅1a明顯大于2階振幅2a,頻響曲線中的實(shí)線和虛線部分分別表示穩(wěn)定解和不穩(wěn)定解,這是由解的形式中兩個(gè)不同分支決定的.

        圖4 立管的1階模態(tài)幅頻響應(yīng)曲線Fig.4 Amplitude-frequency response curve of the 1st modal riser

        圖5 立管的2階模態(tài)幅頻響應(yīng)曲線Fig.5 Amplitude-frequency response curve of the 2nd modal riser

        對(duì)式(38)、(39)分離虛、實(shí)部的結(jié)果進(jìn)行直接數(shù)值求解,可以得到1階模態(tài)坐標(biāo)1q的時(shí)間歷程曲線,見圖6.

        圖6 1階模態(tài)坐標(biāo)1q的時(shí)間歷程曲線Fig.6 Time history curve of the 1st modal cordinate1q

        再對(duì)圖6所示振動(dòng)歷程進(jìn)行傅里葉變換,得到圖7所示振動(dòng)響應(yīng)譜.

        由圖7頻譜分析可得,響應(yīng)成分主要包括3個(gè):其一為參激頻率Ω的1/2的亞諧振動(dòng),即按照頻率0.296,9,rad/s進(jìn)行振動(dòng)的成分,此時(shí)立管產(chǎn)生明顯的亞諧參激振動(dòng);其二為渦激振動(dòng)部分,即按照渦激升力頻率1.609,0,rad/s進(jìn)行振動(dòng)的成分;其三為按照立管1階固有頻率0.199,2,rad/s進(jìn)行振動(dòng)的成分,產(chǎn)生該成分的主要原因是由于本文立管模型為兩端鉸支的梁結(jié)構(gòu),其橫向振動(dòng)的前兩階固有頻率之間存在內(nèi)共振關(guān)系,即212ωω≈,使2階固有頻率相對(duì)1階固有頻率而言恰好滿足參激共振條件,激發(fā)了立管1階模態(tài)的振動(dòng),由此說明此時(shí)立管1階模態(tài)仍為主要激發(fā)模態(tài),且在參激共振和內(nèi)共振共同作用下,亞諧振動(dòng)諧波占主要部分,這表明了組合參數(shù)激勵(lì)對(duì)于立管振動(dòng)的重要影響..從圖4和圖5可得,0σ=時(shí),

        圖7 1階模態(tài)振動(dòng)響應(yīng)譜Fig.7 Spectrum of the 1st modal vibration

        改變頂部升沉幅值,得到1,maxa,如表2所示.由表2中數(shù)據(jù)可知,在參激頻率不變的情況下,立管1階模態(tài)幅值最大值隨浮體升沉幅值的增加而增加,當(dāng)浮體升沉幅值達(dá)到3,m時(shí),幅值響應(yīng)發(fā)生很明顯的增大.立管處于參數(shù)共振區(qū)時(shí)易發(fā)生失穩(wěn),需控制立管振動(dòng)響應(yīng).

        表2 1階模態(tài)幅值的最大值1,maxaTab.2 Maximum response value of the 1st modal1,maxa

        4 結(jié) 論

        (1) 理論分析結(jié)果與數(shù)值模擬分析結(jié)果基本吻合,說明該理論分析結(jié)果是合理的,分析方法是可靠的.

        (2) 當(dāng)立管頂部參激頻率接近立管1階固有頻率與2階固有頻率之和時(shí),1階響應(yīng)明顯大于2階響應(yīng),為主要響應(yīng)成分.

        (3) 參數(shù)激勵(lì)使響應(yīng)中產(chǎn)生明顯的亞諧參激振動(dòng)成分,并且結(jié)構(gòu)本身的內(nèi)共振關(guān)系使2階固有頻率相對(duì)1階固有頻率滿足參激共振條件,在立管振動(dòng)系統(tǒng)中產(chǎn)生重要作用.

        (4) 在參激頻率不變的情況下,立管1階模態(tài)幅值最大值隨浮體升沉幅值的增加而增加,并且當(dāng)升沉幅值較大時(shí),立管橫向振動(dòng)幅值顯著增大.因此在實(shí)際設(shè)計(jì)中,應(yīng)該避免立管頂部發(fā)生較大的升沉幅值.

        (5) 本文選取渦泄頻率遠(yuǎn)離立管橫向振動(dòng)固有頻率,而當(dāng)外激頻率與某階固有頻率接近時(shí),渦激共振與參激共振同時(shí)發(fā)生,可能導(dǎo)致立管振動(dòng)大幅加劇,使立管處于危險(xiǎn)狀態(tài),因此在實(shí)際工程中,還需要關(guān)注兩種諧振同時(shí)發(fā)生的情形,確保立管安全工作.

        [1] 王一飛. 深海立管渦激振動(dòng)疲勞損傷預(yù)報(bào)方法研究[D]. 上海:上海交通大學(xué)船舶海洋與建筑工程學(xué)院,2008. Wang Yifei. Prediction Methods for the VIV-induced Fatigue Damage in Deep Sea Risers [D]. Shanghai:School of Naval Architecture,Ocean and Civil Engineering,Shanghai Jiao Tong University,2008(in Chinese).

        [2] 陳云水,王德禹. 深水立管的若干結(jié)構(gòu)力學(xué)研究進(jìn)展[J]. 中國海洋平臺(tái),2007,22(5):1-5. Chen Yunshui,Wang Deyu. Some advance of research on structural mechanics of deep water risers[J]. China Offshore Platform,2007,22(5):1-5(in Chinese).

        [3] Patel M H,Park H I. Dynamics of tension leg platform tethers at low tension (Part I):Mathieu stability at large parameters[J]. Marine Structures,1991,4(3):257-273.

        [4] Patel M H,Park H I. Combined axial and lateral responses of tensioned buoyant platform tethers[J]. Engineering Structures,1995,17(10):687-695.

        [5] Park H I,Jung D H. A finite element method for dynamic analysis of long slender marine structures under combined parametric and forcing excitations[J]. Ocean Engineering,2002,29(11):1313-1325.

        [6] Chatjigeorgiou I K. On the parametric excitation of vertical elastic slender structures and the effect of damping in marine applications[J]. Applied Ocean Research,2004,26(1/2):23-33.

        [7] Chatjigeorgiou I K,Mavrakos S A. Nonlinear resonances of parametrically excited risers numerical and analytical investigation for12Ωω?[J]. Computers and Structures,2005,83(8/9):560-573.

        [8] 朱克強(qiáng). 駁船升沉和橫蕩運(yùn)動(dòng)下海洋立管動(dòng)力響應(yīng)[J]. 中外船舶科技,2002(3):1-4. Zhu Keqiang. Dynamic response of submarine riser connected to a vessel with heave and lateral motion[J]. Shipbuilding Science and Technology,2002(3):1-4(in Chinese).

        [9] 楊和振,李華軍. 參數(shù)激勵(lì)下深海立管動(dòng)力特性研究[J]. 振動(dòng)與沖擊,2009,28(9):65-69. Yang Hezhen,Li Huajun. Vibration analysis of deep-sea risers under parametric excitations[J]. Journal of Vibration and Shock,2009,28(9):65-69(in Chinese).

        [10] 邵衛(wèi)東,唐友剛,樊娟娟,等. 考慮浮體升沉及張緊環(huán)運(yùn)動(dòng)深海立管固有振動(dòng)特征研究[J]. 海洋工程,2012,30(2):8-13. Shao Weidong,Tang Yougang,F(xiàn)an Juanjuan,et al. Study of natural vibration characteristics of deep-water riser considering heave motion of platform and tensionring’s motion[J]. The Ocean Engineering,2012,30(2):8-13(in Chinese).

        [11] Paidoussis M P,Li G X. Pipes conveying fluid:A model dynamical problem[J]. Journal of Fluids and Structures,1993,7(2):137-204.

        [12] 李楊青. 深海頂張力立管非線性振動(dòng)特性研究[D]. 天津:天津大學(xué)建筑工程學(xué)院,2013. Li Yangqing. Study on the Nonlinear Vibration Behaviors of Top-Tensioned Risers in Deep Water[D]. Tianjin:School of Civil Engineering,Tianjin University,2013(in Chinese).

        (責(zé)任編輯:樊素英)

        Nonlinear Vibration Analysis of Combination Parametric Resonance for TTRs in Deep Water

        Tang Yougang1,2,Zhu Longhuan1,2,Li Yangqing1,2
        (1. State Key Laboratory of Hydraulic Engineering Simulation and Safety,Tianjin University,Tianjin 300072,China;2. School of Civil Engineering,Tianjin University,Tianjin 300072,China)

        The nonlinear vibration characteristics induced by parameters and vortex were investigated for toptensioned risers (TTRs) in deep water. The vibration equations were established with consideration of both the parametric excitation caused by the heave motion of platform and vortex induced force. Approximated analytical solution to the equations was calculated using the method of multiple scales. Considering the case of combination parametric resonance,Ω≈ω1+ω2, the frequency response curves of the riser were obtained and the influence of the heave motion of platform on the nonlinear vibration of the riser was analyzed. The results show that when the parametric excitation frequency satisfies the combination parametric resonance conditions, the sub-harmonic resonances (response frequency equal to 1/2 parametric frequency) is obvious. Because of the existence of internal resonance, the vibration response of the first mode is excited and the amplitude is much larger than that of the second mode. As the heave motion amplitude increases, the amplitude of the lateral vibration response of TTRs increases obviously, which indicates that the heave motion of platform has important influence on the lateral vibration of TTRs。

        top-tensioned risers(TTRs);combination parametric resonance;method of multiple scales;vortex induced vibration

        P756.2

        A

        0493-2137(2015)09-0811-06

        10.11784/tdxbz201402022

        2014-02-18;

        2014-09-02.

        國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)資助項(xiàng)目(2014CB0468);國家自然科學(xué)基金重點(diǎn)資助項(xiàng)目(51239008).

        唐友剛(1952— ),博士,教授.

        唐友剛,tangyougang_td@163.com.

        猜你喜歡
        渦激浮體立管
        浮體結(jié)構(gòu)沉浮過程周圍水流特性研究
        人民長江(2023年6期)2023-07-25 12:24:14
        不同間距比下串聯(lián)圓柱渦激振動(dòng)數(shù)值模擬研究
        物探船硬浮體陣列自擴(kuò)變量分析與應(yīng)用
        超大型浮體結(jié)構(gòu)碰撞損傷研究
        渦激振動(dòng)發(fā)電裝置及其關(guān)鍵技術(shù)
        常見高層建筑物室內(nèi)給水立管材質(zhì)解析
        有限流動(dòng)水域浮體受力及側(cè)傾研究
        盤球立管結(jié)構(gòu)抑制渦激振動(dòng)的數(shù)值分析方法研究
        電子制作(2018年14期)2018-08-21 01:38:42
        深水鋼懸鏈立管J型鋪設(shè)研究
        柔性圓管在渦激振動(dòng)下的模態(tài)響應(yīng)分析
        美女被躁到高潮嗷嗷免费观看| 日本一区二区三区熟女俱乐部 | 日韩一区二区中文天堂| 337p日本欧洲亚洲大胆色噜噜 | 日本精品熟妇一区二区三区| 国产精品久久久在线看| 国产精品久久久国产盗摄| 国产又爽又黄的激情精品视频| 亚洲日韩AV无码美腿丝袜| 国产三级av在线播放| 亚洲最大中文字幕在线| 国模丽丽啪啪一区二区| 中文字幕日韩高清| 精品一区二区三区久久久| 色小姐在线视频中文字幕| 免费人成在线观看视频高潮| 人妻少妇精品视频一区二区三区| 性一交一乱一伦一视频一二三区| 日韩精品午夜视频在线| 国产日产一区二区三区四区五区| 麻豆╳╳╳乱女另类| 午夜精品久久久久久久久久久久| 性做久久久久久久| 91中文在线九色视频| 日日噜噜夜夜狠狠va视频v| 久久久久成人片免费观看蜜芽| 色婷婷资源网| 精品国产车一区二区三区| 亚洲精品中文字幕视频色| 国产后入清纯学生妹| 国产成人av性色在线影院色戒 | 成年女人18毛片毛片免费| 蜜臀av一区二区三区| 91九色免费视频网站| 久久国产精品-国产精品| 色老头一区二区三区| 免费在线观看视频专区| 国产香蕉视频在线播放| 麻豆高清免费国产一区| 中国人妻沙发上喷白将av| av男人的天堂亚洲综合网|