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

不同流場下含內(nèi)流立管渦激振動響應(yīng)特性及Coriolis力效應(yīng)研究

2023-02-22 14:29:42李星輝袁昱超薛鴻祥唐文勇
振動與沖擊 2023年3期
關(guān)鍵詞:效應(yīng)振動區(qū)域

李星輝, 袁昱超, 薛鴻祥, 唐文勇

(1.上海交通大學(xué) 海洋工程國家重點(diǎn)實(shí)驗(yàn)室,上海 200240;2.上海交通大學(xué) 船舶海洋與建筑工程學(xué)院,上海 200240)

海洋立管是海洋油氣集輸?shù)闹匾考m敳繌埦o式(頂張式)立管是海洋開發(fā)中常用的立管類型之一,其底部與海底井口相連,頂部與平臺張力系統(tǒng)連接,如圖1所示。頂張力使立管保持垂直狀態(tài),避免其在外力作用下失穩(wěn)。對于細(xì)長結(jié)構(gòu)物而言,海洋來流會在結(jié)構(gòu)兩側(cè)產(chǎn)生周期性脫落的漩渦,產(chǎn)生渦激振動(vortex-induced vibration,VIV),進(jìn)而誘發(fā)疲勞損傷甚至結(jié)構(gòu)破壞。立管的渦激振動問題得到了國內(nèi)外學(xué)者的持續(xù)關(guān)注,已有大量相關(guān)的試驗(yàn)研究和數(shù)值模擬。數(shù)值模擬大體上可以分為頻域和時域兩類方法。頻域方法較為便捷,適合分析簡單的線性動力學(xué)問題。時域方法應(yīng)用更為廣泛,適合分析復(fù)雜的非線性動力學(xué)問題。常用的立管的時域分析模型有計(jì)算流體動力學(xué)(computational fluid dynamics,CFD)模型、尾流振子模型和流體力分解模型等。CFD方法較為精確,但計(jì)算成本高,方法較為復(fù)雜,更加適合機(jī)理性探究。尾流振子法可以得到和試驗(yàn)數(shù)據(jù)較為接近的結(jié)果,但需要調(diào)整較多的經(jīng)驗(yàn)參數(shù)[1]。相比之下,流體力分解法更為簡便,且可以準(zhǔn)確預(yù)報立管在多種外流場和時變激勵下的渦激振動響應(yīng)。本文主要采用流體力分解模型進(jìn)行立管渦激振動響應(yīng)的時域模擬。

當(dāng)海洋平臺正常工作時,立管內(nèi)部充滿以一定速度提升的油氣等介質(zhì),內(nèi)部流動的介質(zhì)會對立管結(jié)構(gòu)產(chǎn)生內(nèi)部作用力。柳博瀚等[2]的研究表明內(nèi)流通過增加軸向動力和改變彈性管道局部曲率的方式影響管道的振動頻率。因此研究立管的動力響應(yīng)時,內(nèi)流的影響不容忽視。含內(nèi)流立管的渦激振動響應(yīng)是更為復(fù)雜的非線性動力學(xué)問題。Guo等[3]開展了含內(nèi)流立管渦激振動的試驗(yàn)研究,試驗(yàn)結(jié)果表明內(nèi)流會降低立管渦激振動的響應(yīng)頻率,增大響應(yīng)振幅,隨后學(xué)者的數(shù)值研究大多也都驗(yàn)證了這一規(guī)律。根據(jù)側(cè)重點(diǎn)的不同,立管內(nèi)流效應(yīng)的相關(guān)研究大致可分為3類:內(nèi)流本身的屬性、立管的布置形式和外部激勵等。內(nèi)流本身的屬性包括內(nèi)流流速、密度、流體性質(zhì)等。對于內(nèi)外流同時作用下的頂張式立管,當(dāng)內(nèi)流流速處于亞臨界區(qū)時,立管仍然呈現(xiàn)出周期性的渦激振動;當(dāng)內(nèi)流流速處于超臨界區(qū)時,立管的響應(yīng)呈現(xiàn)混沌、分叉等新的特性[4]。Meng等[5-6]研究了立管渦激振動響應(yīng)頻率和振幅隨內(nèi)流速度和密度的變化規(guī)律。劉曉強(qiáng)等[7]發(fā)現(xiàn)黏性內(nèi)流對立管動力特性和響應(yīng)峰值都有所影響,但在實(shí)際工程中的影響有限。馬天麒等[8]將內(nèi)流的狀態(tài)從單相拓展到多相,研究了氣液兩相內(nèi)流的效應(yīng)。立管的布置形式有頂張式、懸臂式、懸鏈線和緩波形等。劉震等[9-10]研究了內(nèi)流對懸鏈線立管和緩波形立管動力特性的影響。Liu等[11-12]分別研究了懸臂式立管和懸鏈線立管在內(nèi)流影響下的渦激振動響應(yīng)特征。立管的外部激勵主要來源于頂部浮式平臺的運(yùn)動。吳天昊等[13-14]分別通過試驗(yàn)和數(shù)值方法研究了立管在內(nèi)流與平臺運(yùn)動聯(lián)合作用下的動力響應(yīng)特征。此外,內(nèi)流的存在可能會減少立管的疲勞壽命[15-16],因此在工程實(shí)際中內(nèi)流的影響不容忽視。

盡管含內(nèi)流立管的渦激振動問題已有一定程度的研究,但目前大多數(shù)相關(guān)研究的外部流場為均勻流。而真實(shí)的立管所處海洋環(huán)境中的外部流場更加接近剪切流的形式,剪切流下立管的渦激振動呈現(xiàn)更為復(fù)雜的多頻響應(yīng)特征。從機(jī)理層面分析,內(nèi)流的效應(yīng)體現(xiàn)為慣性力、離心力和科氏力三方面。慣性力項(xiàng)和離心力項(xiàng)中,位移分別僅對時間或空間求偏導(dǎo),而科氏力項(xiàng)中位移同時對兩者求偏導(dǎo),導(dǎo)致科氏力的效應(yīng)更為復(fù)雜。為進(jìn)一步分析內(nèi)流對立管渦激振動的影響,本文分別研究了含內(nèi)流頂張式立管在均勻流和剪切流作用下的渦激振動響應(yīng)特性,并重點(diǎn)分析了內(nèi)流影響機(jī)理中相對復(fù)雜的科氏力作用。

1 含內(nèi)流立管渦激振動時域模型

本文采用Cartesian坐標(biāo)系統(tǒng)來描述立管模型,其中x軸方向?yàn)閬砹?順流)方向,y軸方向垂直于來流(橫流)方向,z軸方向豎直向上。典型的含內(nèi)流立管渦激振動模型示意圖,如圖1所示。

圖1 浮式平臺-立管-海床系統(tǒng)示意圖

1.1 渦激振動流體力載荷

流體力分解模型將渦激振動流體力分解為與速度同相位的FV和與加速度同相位的FM

(1)

式中:CV為渦激振動激勵力系數(shù);ACF和fCF分別為橫流方向渦激振動的振幅和頻率;ρf為外流密度;D為立管外徑;V為外流流速;Ca為附加質(zhì)量系數(shù),本文的附加質(zhì)量系數(shù)設(shè)為1.0。

CV的值可以從Gopalkrishnan等[17]通過試驗(yàn)繪制的渦激振動激勵力系數(shù)云圖中得到,如圖2所示。一組渦激振動頻率和振幅對應(yīng)云圖中的一個數(shù)值,正的數(shù)值為激勵力系數(shù)CV,負(fù)的數(shù)值轉(zhuǎn)化為水動力阻尼系數(shù)cf。當(dāng)響應(yīng)頻率或振幅超過云圖范圍時,采用Venugopal[18]的阻尼模型來計(jì)算水動力阻尼系數(shù)cf。

圖2 渦激振動激勵力系數(shù)云圖

1.2 含內(nèi)流立管渦激振動微分方程

不考慮重力、阻尼力和張力等外力時,含內(nèi)流管道的運(yùn)動方程[19]為

(2)

式中:m為立管單位長度的質(zhì)量;mi為單位長度內(nèi)流的質(zhì)量;Ui為內(nèi)部流體的流速;E和I分別為立管的楊氏模量和截面慣性矩。

含內(nèi)流立管的完整運(yùn)動方程為

(3)

式中:T為立管的有效張力;c為立管的阻尼。總阻尼c可分為水動力阻尼cf和結(jié)構(gòu)阻尼cs=4πmfCFξ,其中ξ為結(jié)構(gòu)阻尼比。

本文采用HHT(Hilber-Hughes-Taylor)方法求解式(3)來得到立管結(jié)構(gòu)的渦激振動響應(yīng)。HHT方法是一種改進(jìn)的Newmark-β法,廣泛應(yīng)用于動力學(xué)問題的模擬[20]。式(4)定義了每一時間分析步內(nèi)結(jié)構(gòu)位置和速度的更新過程。

(4)

1.3 模型驗(yàn)證

本文采用的立管渦激振動時域模型已被驗(yàn)證過可準(zhǔn)確預(yù)報立管在均勻流、剪切流和振蕩流等多種外流場下的渦激振動響應(yīng),詳見文獻(xiàn)[21]。為驗(yàn)證本文對內(nèi)流效應(yīng)的準(zhǔn)確模擬,與Duan等研究的數(shù)值模擬結(jié)果進(jìn)行對比。Duan等的研究中的立管模型主要參數(shù),如表1所示。外流取為流速1.6 m/s的均勻流。

表1 立管模型參數(shù)[6]

從不同內(nèi)流流速下立管的固有頻率和渦激振動均方根位移變化兩方面進(jìn)行對比,對比結(jié)果如圖3和圖4所示。

圖3 不同內(nèi)流流速下立管固有頻率對比

(a) 本文計(jì)算結(jié)果

隨著內(nèi)流流速的升高,立管的各階固有頻率都會有所降低。且當(dāng)內(nèi)流速度較低時,固有頻率降低幅度不大;當(dāng)內(nèi)流速度較高時,固有頻率顯著降低。隨著內(nèi)流流速的升高,立管的均方根位移增大。且內(nèi)流流速較低時,均方根位移增幅較小;當(dāng)內(nèi)流流速較高時,均方根位移增幅顯著。

本文預(yù)報結(jié)果和Duan等的研究有較好的一致性,證明了本文模型對立管內(nèi)流效應(yīng)模擬的準(zhǔn)確性。

2 均勻流下含內(nèi)流立管渦激振動響應(yīng)特征

1.3節(jié)中立管模型長細(xì)比較小(約為255),而實(shí)際海洋立管的長細(xì)比通常在1 000以上,故有必要分析大長細(xì)比立管在均勻流下的渦激振動響應(yīng)特征及內(nèi)流效應(yīng)。本節(jié)選取表2中的38 m長立管進(jìn)行渦激振動時域模擬(長細(xì)比約1 407),考慮了多個外流流速,并重點(diǎn)分析科氏力的效應(yīng)。38 m立管模型取自文獻(xiàn)[22]。

表2 大長細(xì)比立管參數(shù)

2.1 渦激振動響應(yīng)特性(均勻流)

立管在不同內(nèi)流流速下的固有頻率變化,如圖5所示。與圖3相比,大長細(xì)比立管的固有頻率本身較小,且隨著內(nèi)流流速升高,固有頻率降低的程度相對更為緩和。固定外流為0.25 m/s均勻流,不同內(nèi)流流速下立管無因次均方根位移曲線和立管位移時歷快速傅里葉變換(fast Fourier transform,F(xiàn)FT)結(jié)果,分別如圖6和圖7所示。隨著內(nèi)流流速的升高,立管的均方根位移增大,響應(yīng)頻率降低,這與長細(xì)比較小的立管規(guī)律相同。下文重點(diǎn)分析不同外流及內(nèi)流速度下科氏力的作用。

圖5 不同內(nèi)流流速下的立管固有頻率

圖6 不同內(nèi)流流速下的立管均方根位移(0.25 m/s均勻流)

圖7 不同內(nèi)流流速下立管時歷曲線FFT結(jié)果(0.25 m/s均勻流)

2.2 科氏力效應(yīng)(均勻流)

取0.15 m/s,0.25 m/s和0.4 m/s 3種均勻外流,固定內(nèi)流流速為5 m/s和20 m/s,分別計(jì)算在立管數(shù)值模型中保留科氏力和去掉科氏力的渦激振動響應(yīng)。其中0.4 m/s外流、20 m/s內(nèi)流下立管均方根位移的結(jié)果對比,如圖8所示。科氏力并不會引起均方根位移的顯著改變,而是微調(diào)振型形狀。為了更加清晰地表示科氏力的影響,將考慮與不考慮科氏力的無因次均方根位移曲線作差,如圖9所示。

圖8 有無科氏力情況下的立管均方根位移(20 m/s內(nèi)流)

相同外流流速下,立管的內(nèi)流流速越高,科氏力越大,均方根位移的改變量也相應(yīng)增大。沿著立管的軸向,科氏力會在一些區(qū)域使得均方根位移有所升高,而在其余區(qū)域又使得均方根位移有所降低,且呈現(xiàn)出一定的規(guī)律性。在均方根位移振型的每一個“山峰”內(nèi),靠近立管底部區(qū)域的均方根位移降低,靠近立管頂部區(qū)域的均方根位移升高。

值得注意的是,當(dāng)模態(tài)階數(shù)較高(振型“山峰”數(shù)量較多)時,最靠近頂部的振型“山峰”內(nèi)均方根位移升高的區(qū)域可能會很小甚至不存在,且從圖9(b)可見內(nèi)流流速越高,這一特性越明顯。

(a) 0.15 m/s均勻外流

3 剪切流下含內(nèi)流立管渦激振動響應(yīng)特征

為分析剪切流下的內(nèi)流效應(yīng),本章選取表2中的90 m長立管(模型取自文獻(xiàn)[23])進(jìn)行渦激振動時域模擬。外流取頂流速為0.6 m/s的剪切來流。

3.1 渦激振動響應(yīng)特性(剪切流)

不同內(nèi)流流速下立管的均方根位移和響應(yīng)頻率的變化圖,分別如圖10、圖11所示。剪切流下立管的渦激振動不再是單頻的,而是呈現(xiàn)出多頻響應(yīng)的特征,且外流確定時,不同內(nèi)流流速下立管的主導(dǎo)響應(yīng)頻率都位于一定的鎖定區(qū)間內(nèi),但各階響應(yīng)頻率的大小和占比都會有所區(qū)別。當(dāng)內(nèi)流流速為0和30 m/s時,立管的均方根分別呈現(xiàn)出第12和第13階的振型。這表明內(nèi)流流速較高時,立管的主導(dǎo)模態(tài)振型會有所升高,這是固有頻率降低導(dǎo)致的。

圖10 不同內(nèi)流流速下的立管均方根位移(0.6 m/s剪切流)

圖11 不同內(nèi)流流速下立管中點(diǎn)時歷曲線FFT結(jié)果(0.6 m/s剪切流)

3.2 科氏力效應(yīng)(剪切流)

為分析剪切流下的科氏力效應(yīng),同樣分別計(jì)算保留科氏力和去掉科氏力的立管渦激振動響應(yīng)結(jié)果。立管中點(diǎn)位移頻響譜,如圖12所示。科氏力并不會改變響應(yīng)頻率的大小,但會改變頻譜圖中的幅值分布情況。

圖12 有無科氏力情況下立管中點(diǎn)時歷曲線FFT結(jié)果

考慮與不考慮科氏力情況下立管的均方根位移曲線,如圖13所示。科氏力同樣對響應(yīng)振型產(chǎn)生了一定的調(diào)制作用,大體呈現(xiàn)出提高立管頂部區(qū)域均方根位移、降低立管中部區(qū)域均方根位移的規(guī)律。內(nèi)流流速50 m/s下科氏力的效應(yīng)更為明顯,這是由于科氏力大小與內(nèi)流速度相關(guān),內(nèi)流速度越大,科氏力效應(yīng)越明顯。

(a)

4 科氏力效應(yīng)機(jī)理

第2章通過數(shù)值模擬,直觀展示了內(nèi)流對不同流場下立管渦激振動的影響。本章將從作用機(jī)理上分析內(nèi)流效應(yīng),并重點(diǎn)關(guān)注較為復(fù)雜的科氏力作用。

4.1 科氏力作用規(guī)律

從科氏力的基本表達(dá)式F=2mv×ω出發(fā)分析其效應(yīng)。其中v和ω分別為內(nèi)流的流速和立管的角速度矢量。圖14中的實(shí)線和虛線分別為假定的立管某處位移和轉(zhuǎn)角變化的時歷曲線。將其分為A,B,C,D4個典型區(qū)域。

圖14 科氏力作用機(jī)理示意圖

區(qū)域A內(nèi)立管位移增大,轉(zhuǎn)角減小,可以判斷出科氏力方向與立管位移同向,表明科氏力做正功。同理可知區(qū)域C內(nèi)科氏力也做正功,區(qū)域B和D內(nèi)科氏力做負(fù)功。不難得到:當(dāng)立管單元位移變化趨勢與轉(zhuǎn)角變化趨勢相反時,科氏力方向與位移方向相同,科氏力做正功;反之科氏力做負(fù)功。

4.1.1 均勻流工況

結(jié)合圖9,由于均勻流下立管的響應(yīng)是單頻的,振型較為規(guī)則。振型中每個“山峰”的左半?yún)^(qū)域與圖14中的B,D區(qū)域類似,科氏力做負(fù)功,故均方根位移降低;而振型中每個“山峰”的右半?yún)^(qū)域與圖14中的A,C區(qū)域類似,科氏力做正功,故均方根位移升高。科氏力在靠近立管頂部的區(qū)域做正功,均方根位移應(yīng)該會有所升高,但從圖9(c)可見,立管頂部區(qū)域的均方根位移反而降低,這是由于存在邊界效應(yīng)。圖15為某一時刻立管的振動形態(tài),實(shí)線代表立管的振動速度方向,虛線代表科氏力方向。在頂部邊界,科氏力做正功。在臨近的左側(cè)兩個區(qū)域,盡管科氏力分別做負(fù)功和正功,但兩個區(qū)域內(nèi)科氏力的方向是相同的,使立管產(chǎn)生向下運(yùn)動的趨勢,如圖15中的虛線所示,從而減小了頂部區(qū)域的振動幅值,最終體現(xiàn)為均方根位移的降低。由圖9(b)可知,當(dāng)內(nèi)流流速為5 m/s時,立管頂部均方根位移有所升高,而當(dāng)內(nèi)流流速為20 m/s時,立管頂部均方根位移只有輕微升高,表明內(nèi)流速度越高,這種邊界效應(yīng)越明顯。

圖15 科氏力的邊界效應(yīng)

4.1.2 剪切流工況

剪切流下立管的渦激振動響應(yīng)頻率較多,振動形態(tài)類似多個模態(tài)振型的疊加。此時科氏力的效應(yīng)無法像均勻流工況一樣直接從振型角度予以解釋。故直接采用式(5)直接計(jì)算一段時間內(nèi)科氏力在某一位置做功的大小。

(5)

一段時間內(nèi)科氏力做功及立管的轉(zhuǎn)角速度均方根值(root mean square,RMS)沿管長的分布情況,如圖16所示。可見剪切流下,立管頂部的轉(zhuǎn)角速度較大。此外剪切流下立管頂部區(qū)域的均方根位移較大,相應(yīng)的振動速度較高,從而科氏力在頂部區(qū)域做功最大。且整體來看頂部區(qū)域內(nèi)科氏力所做的正功大于負(fù)功,因此頂部區(qū)域的均方根位移會有所升高。Pa?doussis的研究和本文的計(jì)算結(jié)果均表明,對于兩端簡支的管,科氏力對整個管做功為零。因此在立管的其余區(qū)域,科氏力整體做負(fù)功,均方根位移會有所降低。即剪切流下科氏力呈現(xiàn)出增大頂部區(qū)域位移、降低中部區(qū)域位移的效應(yīng),如圖13所示。

(a)

4.2 內(nèi)流反向的影響

4.1節(jié)假定的內(nèi)流流動方向是從底部到頂部,如圖1所示。本節(jié)通過改變內(nèi)流方向進(jìn)一步探究科氏力效應(yīng)機(jī)理。由式(2)可知,內(nèi)流方向的變化只會影響科氏力。將內(nèi)流反向即從頂部到底部流動,分別計(jì)算考慮與不考慮科氏力情況下立管的渦激振動響應(yīng)。

0.15 m/s均勻外流下38 m立管保留和去除科氏力的均方根位移差曲線,如圖17所示。與圖9對比可知科氏力效應(yīng)發(fā)生了以立管中點(diǎn)橫截面為中心的對稱變換。剪切外流下90 m立管保留和去除科氏力的均方根位移對比曲線,如圖18所示。與圖13對比可知科氏力的作用在內(nèi)流反向時大體呈現(xiàn)出降低立管頂部均方根位移、提高立管中部均方根位移的規(guī)律。內(nèi)流反向?qū)?dǎo)致科氏力同時反向,從而原本科氏力做正功的區(qū)域變?yōu)榱丝剖狭ψ鲐?fù)功,反之亦然。因此科氏力對振型的調(diào)制規(guī)律也發(fā)生了改變。

圖17 有無科氏力情況下的立管均方根位移差值(內(nèi)流反向)

圖18 有無科氏力情況下的立管均方根位移(內(nèi)流反向)

4.3 內(nèi)流速度的影響

(a) 0.25 m/s均勻流

5 結(jié) 論

基于含內(nèi)流立管的渦激振動響應(yīng)時域預(yù)報模型,本文分別研究了均勻外流和剪切外流下內(nèi)流對頂張式立管渦激振動響應(yīng)特性的影響。從作用機(jī)理上重點(diǎn)分析了內(nèi)流引起的復(fù)雜科氏力效應(yīng)。得到的主要結(jié)論如下:

(1) 內(nèi)流對立管渦激振動的影響體現(xiàn)在慣性力、科氏力和離心力三方面。慣性力會使立管單位長度質(zhì)量增加;離心力會降低立管的有效張力;而科氏力會改變立管的振動形態(tài),不會改變響應(yīng)頻率。

(2) 均勻外流下,隨著內(nèi)流流速的升高,立管響應(yīng)頻率降低,均方根位移增大;剪切外流下,立管呈現(xiàn)多頻響應(yīng)特征,主導(dǎo)頻率位于一定的鎖定區(qū)間,且隨著內(nèi)流流速的升高,主導(dǎo)模態(tài)階數(shù)有上升趨勢。

(3) 均勻流下,科氏力在每個振型“山峰”內(nèi)左側(cè)和右側(cè)區(qū)域分別做負(fù)功和正功,并引起均方根位移減小和增大;剪切流下,科氏力使得立管頂部區(qū)域均方根位移增大,中部區(qū)域均方根位移減小。

(4) 內(nèi)流流向改變不影響慣性力和離心力,僅影響科氏力。內(nèi)流反向時,均勻流下,科氏力在每個振型“山峰”內(nèi)的左側(cè)和右側(cè)區(qū)域分別做正功和負(fù)功,并引起均方根位移的增大和減小;剪切流下,科氏力使得立管頂部區(qū)域均方根位移減小,中部區(qū)域均方根位移增大。

(5) 內(nèi)流流速越大,均勻外流和剪切外流下科氏力的效應(yīng)都會更加明顯。為準(zhǔn)確預(yù)報立管的渦激振動,保證結(jié)構(gòu)安全性,建議實(shí)際工程中考慮內(nèi)流效應(yīng)時應(yīng)充分考慮科氏力的效應(yīng)。

綜上,本文分別從數(shù)值模擬和機(jī)理分析角度討論了內(nèi)流效應(yīng),研究結(jié)論有助于理解內(nèi)流對海洋立管渦激振動的影響規(guī)律,并可為立管實(shí)際工程設(shè)計(jì)提供一定參考。

同時,本文還存在一定的局限和不足。工程實(shí)際中的內(nèi)流并不是單相流,而是多相流的形式,建議今后的工作將內(nèi)流效應(yīng)擴(kuò)展到多相內(nèi)流效應(yīng)。另外本文的研究對象為頂張式立管,科氏力對其他形式立管(懸鏈線立管、緩波形立管等)的效應(yīng)有待進(jìn)一步研究。

猜你喜歡
效應(yīng)振動區(qū)域
振動的思考
鈾對大型溞的急性毒性效應(yīng)
懶馬效應(yīng)
振動與頻率
中立型Emden-Fowler微分方程的振動性
應(yīng)變效應(yīng)及其應(yīng)用
關(guān)于四色猜想
分區(qū)域
基于嚴(yán)重區(qū)域的多PCC點(diǎn)暫降頻次估計(jì)
電測與儀表(2015年5期)2015-04-09 11:30:52
UF6振動激發(fā)態(tài)分子的振動-振動馳豫
主站蜘蛛池模板: 亚洲天堂久久新| 久久久久久久久亚洲精品| 中文字幕佐山爱一区二区免费| 国产精品网址在线观看你懂的| 亚洲国产精品成人久久综合影院| 亚洲精品色AV无码看| 免费国产不卡午夜福在线观看| 玖玖精品在线| 色悠久久久| 99热免费在线| 欧美日韩中文国产va另类| 毛片网站在线看| 国产精品13页| 高清亚洲欧美在线看| 五月丁香在线视频| 原味小视频在线www国产| 中文字幕亚洲乱码熟女1区2区| 人妻出轨无码中文一区二区| 亚洲综合极品香蕉久久网| 国产一级视频在线观看网站| 免费观看亚洲人成网站| 日韩欧美网址| 国产女人爽到高潮的免费视频 | 99精品免费欧美成人小视频| 亚洲人成网址| 国产精品无码一二三视频| 国产一区二区三区精品欧美日韩| 亚洲精品国产精品乱码不卞| 国产精品制服| 国产高清国内精品福利| 日韩一区二区在线电影| 国产欧美日韩综合一区在线播放| 亚洲欧美不卡| 免费看a级毛片| 天堂在线亚洲| 人妻中文字幕无码久久一区| 沈阳少妇高潮在线| 欧美国产中文| 精品福利国产| www.youjizz.com久久| 国产嫖妓91东北老熟女久久一| 无码人妻热线精品视频| 99久久精品久久久久久婷婷| 亚洲第一视频网站| 国产精品女主播| 亚洲男人天堂2018| 91探花国产综合在线精品| 中文天堂在线视频| 日本高清成本人视频一区| 中文字幕va| 日韩在线2020专区| 国产精品无码久久久久AV| 精品无码一区二区三区电影| 99视频在线免费| 亚洲综合久久一本伊一区| 九色视频在线免费观看| 中文字幕欧美日韩| 亚洲一区免费看| 国产精品午夜电影| 在线看国产精品| 亚洲成网777777国产精品| 成人国产精品2021| 99国产精品一区二区| 91亚瑟视频| 国产精品对白刺激| 国产成人AV大片大片在线播放 | 国产农村1级毛片| 2020国产精品视频| 欧美成人二区| 四虎成人免费毛片| 妇女自拍偷自拍亚洲精品| 在线人成精品免费视频| 欧洲熟妇精品视频| 欧美日韩成人在线观看| 99精品在线看| 尤物午夜福利视频| 色哟哟国产精品| 亚洲一区国色天香| 99久久国产自偷自偷免费一区| 欧美在线一二区| 久久精品丝袜| 毛片在线看网站|