• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    深海采礦系統(tǒng)中懸臂式立管渦激振動(dòng)分析1)

    2022-07-10 13:13:50金國(guó)慶
    力學(xué)學(xué)報(bào) 2022年6期
    關(guān)鍵詞:渦激立管懸臂

    金國(guó)慶 鄒 麗 ,?,2) 宗 智 ,? 孫 哲 王 浩

    * (大連理工大學(xué)船舶工程學(xué)院,工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室,遼寧大連 116024)

    ? (高技術(shù)船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心,上海 200240)

    ** (中國(guó)水利水電科學(xué)研究院,北京 100038)

    引言

    海底蘊(yùn)藏著豐富的礦產(chǎn)資源,目前已知的包括多金屬結(jié)核物、多金屬硫化物、富鈷結(jié)殼和深海稀土資源等,具有重要戰(zhàn)略意義與商業(yè)開(kāi)采價(jià)值[1-2].這類(lèi)礦產(chǎn)資源大多富含于數(shù)千米水深的深遠(yuǎn)海底部,賦存狀態(tài)的特殊性和周?chē)h(huán)境的復(fù)雜性導(dǎo)致開(kāi)采難度較大[3-4].目前主流的深海采礦方式為垂直管道混輸系統(tǒng),即采用水力提升的方式將礦物通過(guò)長(zhǎng)距離的管道輸送至水面采礦船.該系統(tǒng)主要由海底采礦系統(tǒng)、垂直輸運(yùn)系統(tǒng)和水面支持系統(tǒng)三部分組成[3].如圖1 所示,海底采礦車(chē)收集礦物并通過(guò)軟管輸送到中間倉(cāng)(緩沖倉(cāng)),之后利用水力提升泵將礦物通過(guò)垂直提升管道輸送至水面母船[5-9].

    圖1 深海采礦系統(tǒng)示意圖Fig.1 Sketch of the deep-sea mining system

    礦產(chǎn)混輸系統(tǒng)中的關(guān)鍵裝備之一為具有大長(zhǎng)細(xì)比特征的垂直提升管道.不同于傳統(tǒng)的海洋立管,垂直提升管道底部與中間倉(cāng)相連接,其底部是自由運(yùn)動(dòng)的狀態(tài),因此該輸運(yùn)管道通常被稱(chēng)為懸臂式或自由懸掛式立管.類(lèi)似的立管形式在油氣資源開(kāi)采系統(tǒng)中也會(huì)遇到,當(dāng)遭遇極端海況或者回收輸油管道時(shí),輸油管道底部脫離井口,呈現(xiàn)自由懸掛式狀態(tài)[10-12].在深海采礦垂直提升管道的研究中,底部中間倉(cāng)重量可簡(jiǎn)化為一個(gè)垂直向下的作用力,立管底部可視為弱約束狀態(tài).然而,目前對(duì)此類(lèi)立管的研究主要集中在結(jié)構(gòu)動(dòng)力學(xué)行為分析[11-15],缺乏對(duì)其水動(dòng)力學(xué)性能的研究.

    與傳統(tǒng)立管類(lèi)似的是,懸臂式立管在作業(yè)過(guò)程中也面臨著渦激振動(dòng)(VIV)等水動(dòng)力學(xué)問(wèn)題.渦激振動(dòng)指的是流體經(jīng)過(guò)鈍體結(jié)構(gòu)物,從物體表面分離后在尾部交替地產(chǎn)生旋渦,進(jìn)而產(chǎn)生周期性的外力和振動(dòng).在一個(gè)質(zhì)量-彈簧-阻尼系統(tǒng)中,當(dāng)振動(dòng)頻率和瀉渦頻率一致時(shí),出現(xiàn)鎖定共振現(xiàn)象,此時(shí)結(jié)構(gòu)物會(huì)產(chǎn)生大幅度的振動(dòng),嚴(yán)重危害作業(yè)安全.此外,長(zhǎng)期的立管振動(dòng)會(huì)導(dǎo)致結(jié)構(gòu)出現(xiàn)疲勞損傷.因此,立管的渦激振動(dòng)研究是立管水動(dòng)力性能分析的重要組成部分.對(duì)于立管的渦激振動(dòng)特性的測(cè)量,傳統(tǒng)的做法是采用縮尺模型進(jìn)行水池試驗(yàn).Chaplin 等[16]針對(duì)一根長(zhǎng)13.12 m 的柔性立管,進(jìn)行了階梯流作用下的VIV 試驗(yàn)研究,觀察到了立管在一定流速下被激發(fā)出高階模態(tài)振動(dòng).Trim 等[17]開(kāi)展了水平布置的38 m 長(zhǎng)的立管分別在均勻流和剪切流作用下的試驗(yàn),分析了不同拖曳速度下的振動(dòng)模態(tài)和頻率等特征.Song 等[18-19]和高云等[20]研究了7.9 m 長(zhǎng)立管在不同流速的均勻流作用下的非均勻分布的水動(dòng)力載荷、立管截面運(yùn)動(dòng)軌跡和疲勞損傷等問(wèn)題.此外,王俊高等[21]研究了正弦振蕩來(lái)流下小尺度模型立管的VIV 特性.徐萬(wàn)海等[22]則結(jié)合拖曳水池試驗(yàn)和數(shù)值模型分析了立管的流體力系數(shù)特性.在懸臂式立管模型試驗(yàn)方面,文獻(xiàn)[10]研究了水面母船運(yùn)動(dòng)對(duì)立管動(dòng)力行為的影響.Fujarra 等[23]研究了懸臂式管道在一階振動(dòng)模態(tài)下的鎖定問(wèn)題.Wang 等[24]分析了船舶周期性運(yùn)動(dòng)對(duì)自由懸掛式立管運(yùn)動(dòng)響應(yīng)的影響.Mao 等[25]測(cè)量了懸掛式排空采油管道在不同流速下的振動(dòng)模態(tài)和振幅.然而,模型試驗(yàn)成本高且周期較長(zhǎng),此外,縮尺模型需要改變立管的材料屬性,產(chǎn)生難以避免的誤差.

    近年來(lái),數(shù)值手段成為研究立管渦激振動(dòng)的主要方式.常見(jiàn)的立管振動(dòng)預(yù)報(bào)數(shù)值方法主要分為三類(lèi):經(jīng)驗(yàn)?zāi)P?、參?shù)化模型和CFD 模型[26].其中,經(jīng)驗(yàn)?zāi)P统2捎梦擦髡褡幽P瞳@得作用在結(jié)構(gòu)上的流體力,進(jìn)一步通過(guò)梁模型求解立管變形[27-29].參數(shù)化模型則是基于一系列的剛性圓柱渦激振動(dòng)試驗(yàn)結(jié)果來(lái)確定立管運(yùn)動(dòng)幅值[30],該方法已形成相關(guān)的計(jì)算軟件,例如Shear 7 和VIVANA 等.可以看出,前兩者分別屬于經(jīng)驗(yàn)和半經(jīng)驗(yàn)方式,因此其計(jì)算成本相對(duì)較低.計(jì)算流體力學(xué)方法則是通過(guò)求解Navier-Stokes (N-S)方程,獲得立管尾流場(chǎng)的非定常流動(dòng)特征和時(shí)域內(nèi)變化的水動(dòng)力載荷.常見(jiàn)的方法是基于切片法的思想,采用CFD 模型求解立管若干節(jié)點(diǎn)的受力,通過(guò)有限元方法求解柔性立管變形和振動(dòng)[31-34],這類(lèi)準(zhǔn)三維的數(shù)值模擬手段可兼顧計(jì)算效率和精度.此外,部分研究采用高精度的大渦模擬LES 方法[35]和直接數(shù)值模擬DNS 方法研究低雷諾數(shù)下柔性立管的渦激振動(dòng)[36-37],但該方法目前所需計(jì)算資源較大,難以求解高雷諾數(shù)下大長(zhǎng)細(xì)比立管的水動(dòng)力問(wèn)題.

    綜上所述,基于CFD 模型和有限元的耦合求解算法,可有效求解具有大長(zhǎng)細(xì)比特征的海洋立管的渦激振動(dòng)問(wèn)題.其中,離散渦方法(DVM)[38]作為一種無(wú)網(wǎng)格的純拉格朗日的流體計(jì)算方法,可以高效率地求解高雷諾數(shù)下的圓柱繞流問(wèn)題,通過(guò)與有限元方法(FEM)結(jié)合,近年來(lái)被發(fā)展用于研究立管的VIV 問(wèn)題[31-32,39-40].此外,目前海洋立管研究主要集中于傳統(tǒng)的兩端鉸支的立管形式,缺乏對(duì)深海采礦系統(tǒng)中采用的懸臂式立管水動(dòng)力性能的研究.本文基于準(zhǔn)三維的DVM-FEM 耦合算法,建立懸臂式立管渦激振動(dòng)的計(jì)算模型,即立管上端固定,底部無(wú)約束且承受一定的向下的拉力.主要目的是系統(tǒng)性地研究不同流速下懸臂式立管的振動(dòng)響應(yīng),分析其振動(dòng)模態(tài)、頻率和振幅等主要特征.本文首先分別介紹了二維離散渦算法和有限元模型.其次,通過(guò)數(shù)值預(yù)報(bào)與經(jīng)典的9.63 m 長(zhǎng)立管在均勻流作用下的試驗(yàn)結(jié)果[41]的對(duì)比,驗(yàn)證了本文采用的耦合算法可以較好地預(yù)報(bào)柔性立管的振動(dòng)響應(yīng).最后,研究不同折合速度下懸臂式立管的VIV 特征,分析其與兩端鉸支立管的異同.本工作可為深海采礦混輸提升系統(tǒng)的設(shè)計(jì)和振動(dòng)行為分析提供參考.

    1 DVM-FEM 耦合算法

    本文所采用的耦合算法基于切片法的思想,將立管沿垂向劃分成若干切片,在每個(gè)切片上通過(guò)二維離散渦算法求解彈性圓柱渦激振動(dòng)問(wèn)題,并獲得該切片所對(duì)應(yīng)的節(jié)點(diǎn)的水動(dòng)力載荷.此外,基于有限元方法建立柔性立管振動(dòng)求解的控制方程,通過(guò)迭代求解獲得立管振動(dòng)的時(shí)域解.綜上,該耦合算法主要包括流體求解器和結(jié)構(gòu)求解器兩部分,所對(duì)應(yīng)的算法分別為離散渦和有限元方法.

    1.1 流體求解器:離散渦方法

    離散渦算法是一種純拉格朗日粒子算法,不需要進(jìn)行復(fù)雜的網(wǎng)格劃分,渦元在物面上生成,運(yùn)動(dòng)后的渦元主要集中在邊界層和尾流區(qū)域,因此具有較低的計(jì)算成本[38,42].根據(jù)質(zhì)量守恒和動(dòng)量守恒定律,離散渦算法的控制方程為連續(xù)性方程和N-S 方程,即

    式中,u為速度矢量,D/Dt為物質(zhì)導(dǎo)數(shù),f為質(zhì)量力,ρ 為流體密度,ν 為運(yùn)動(dòng)黏度.對(duì)式(2)兩側(cè)取旋度,可得到二維渦量輸運(yùn)方程

    式中,ω 是垂直于二維水平面方向上的渦量分量

    式中,u和v分別是速度在x和y方向上的分量.離散渦方法引入流函數(shù) ψ,速度分量可表示為

    將方程(5)代入式(4)中可得Poisson 方程

    Chorin[38]采用算子分裂法的思想,將渦量輸運(yùn)方程分解為兩部分求解,分別為對(duì)流和擴(kuò)散方程

    對(duì)流項(xiàng)

    通過(guò)Biot-Savart 定律可求解對(duì)流項(xiàng)方程,渦元對(duì)流速度由渦元之間的誘導(dǎo)速度和來(lái)流速度兩部分組成,基于Spalart 等[42]提出的渦核模型,數(shù)值計(jì)算過(guò)程中的對(duì)流速度可表示為

    式中,u∞=(u∞,0) 為來(lái)流速度,r=(x,y) 為渦元位置坐標(biāo),Γ 為渦元強(qiáng)度,σ 為渦元半徑,k為二維平面垂向單位矢量.對(duì)于擴(kuò)散項(xiàng),基于Chorin[38]提出的隨機(jī)擴(kuò)散法求解擴(kuò)散方程,即

    式中,?xi和 ?yi分別為第i個(gè)渦元在兩個(gè)方向的隨機(jī)步長(zhǎng),?t為迭代時(shí)間步,Pi和Q分別為介于 (0,1) 和(0,2π)的隨機(jī)數(shù).采用二階Runge-Kutta 法進(jìn)行時(shí)間積分,渦元i在t+?t時(shí)刻的新位置為

    求解物面壓力分布時(shí),由于物體是運(yùn)動(dòng)的,因此需在壓力求解中考慮加速度項(xiàng),即

    式中,p為物面壓力,s表示物面切向,ab為物體運(yùn)動(dòng)加速度.

    離散渦算法具體的實(shí)現(xiàn)過(guò)程可參考文獻(xiàn)[43-44].需要注意的是,渦元在運(yùn)動(dòng)過(guò)程中會(huì)進(jìn)入物體內(nèi)部,本文采用瞬態(tài)渦量守恒方法,將進(jìn)入圓柱內(nèi)部的渦元基于圓周定理映射到物體外部,保證單個(gè)時(shí)間步內(nèi)流場(chǎng)渦量是守恒的,具體實(shí)現(xiàn)步驟可參考文獻(xiàn)[44-45].

    1.2 結(jié)構(gòu)求解器:有限元方法

    海洋柔性立管具有大長(zhǎng)細(xì)比的特征,可簡(jiǎn)化為一根垂直的Euler-Bernoulli 梁模型.通過(guò)Galerkin 方法,可獲得立管結(jié)構(gòu)動(dòng)力響應(yīng)計(jì)算的有限元控制方程[31-32],即

    式中,M,C,K分別為質(zhì)量、阻尼和剛度陣,F為外力矢量,x為立管各節(jié)點(diǎn)組成的位移矢量.其中外力矢量由離散渦算法提供,阻尼采用Rayleigh 阻尼模型,即

    式中,α 和 β 為瑞利阻尼系數(shù),計(jì)算方式為

    式中,ξ 為阻尼比,ω1和 ω2分別為立管的前兩階固有角頻率,可通過(guò)模態(tài)分析獲得,本次計(jì)算模態(tài)分析采用子空間迭代法[46].

    需要注意的是,傳統(tǒng)的兩側(cè)鉸支立管高度為 2 處的張力T(z) 由頂部張力、包含管內(nèi)流體重量的立管結(jié)構(gòu)自重和浮力項(xiàng)三部分組成,即

    式中,Ttop為頂部張力,ztop為立管頂端高度,ws和wb分別為單位長(zhǎng)度立管的自重和所受到的浮力.對(duì)于底部有附加重量的懸臂式立管,其所受張力由參考文獻(xiàn)[25-47]給出

    式中,Wb為立管底部承受的向下的拉力,即深海采礦系統(tǒng)中中間倉(cāng)的重量.

    如圖2 所示,懸臂式立管可簡(jiǎn)化為一個(gè)柔性懸臂式模型,其底部連接的中間倉(cāng)可簡(jiǎn)化為一個(gè)垂直向下的力Wb,需要注意的是該力學(xué)模型忽略了中間倉(cāng)慣性力的影響.未來(lái)工作中可以將中間倉(cāng)簡(jiǎn)化為一個(gè)集中質(zhì)量節(jié)點(diǎn)進(jìn)行結(jié)構(gòu)動(dòng)力響應(yīng)分析,并且可以研究這兩種力學(xué)模型的異同.

    圖2 垂直提升管道簡(jiǎn)化模型Fig.2 Simplified model of the vertical lifting pipeline

    綜上所示,基于切片法的思想,立管被劃分成一系列的梁?jiǎn)卧?基于離散渦方法求解立管各單元節(jié)點(diǎn)所受的外力,采用有限元方程求解柔性立管結(jié)構(gòu)變形,本文采用Newmark 方法[46]迭代求解每個(gè)時(shí)間步內(nèi)的立管振動(dòng)位移、速度和加速度,采用弱耦合的形式實(shí)現(xiàn)流體和結(jié)構(gòu)求解器之間的數(shù)據(jù)交換.在接下來(lái)的內(nèi)容中,本文首先通過(guò)與經(jīng)典的立管模型對(duì)比來(lái)驗(yàn)證該耦合算法的有效性,進(jìn)一步分析懸臂式立管的VIV 特性.

    2 數(shù)值模型驗(yàn)證

    由于懸臂式立管的模型試驗(yàn)基準(zhǔn)算例較少,在這一節(jié)中仍采用頂端和底部均為固定邊界條件的模型試驗(yàn)進(jìn)行數(shù)值驗(yàn)證.基于DVM-FEM 耦合算法的計(jì)算程序由本文作者自主開(kāi)發(fā),采用Fortran 語(yǔ)言編寫(xiě).如表1 所示,物理模型試驗(yàn)[41]中的立管總長(zhǎng)度為9.63 m,立管外部直徑為20 mm,厚度為0.45 mm,頂部張力為817 N,分別計(jì)算均勻流速0.42 m/s 和0.84 m/s 兩個(gè)流速下的立管振動(dòng)響應(yīng).

    表1 柔性立管參數(shù)[32,41]Table 1 Parameters of a flexible riser[32,41]

    本文的耦合算法基于切片法,切片數(shù)量對(duì)計(jì)算成本有重要影響,選取合適的立管單元數(shù)N可在保證計(jì)算精度的基礎(chǔ)上提高求解效率.本次研究選擇了4 套立管單元數(shù)進(jìn)行敏感性分析,數(shù)量分別為30,60,90 和120.在0.42 m/s 來(lái)流速度作用下,不同單元數(shù)計(jì)算得到的立管橫向振動(dòng)位移均方根值如圖3 所示.可以看出,立管單元數(shù)量變化對(duì)其振幅和模態(tài)的影響均較小,選取N=60 時(shí)即可滿(mǎn)足立管單元數(shù)量收斂性要求,因此本文所有模擬采用的立管單元數(shù)均為60.

    圖3 不同單元數(shù)立管的橫向位移均方根值Fig.3 RMS amplitudes of the transverse displacements for the riser models with various element numbers

    兩個(gè)來(lái)流速度下的立管在不同時(shí)刻的橫向振動(dòng)位移包絡(luò)線和對(duì)應(yīng)的均方根值如圖4 所示.可以看出,在0.42 m/s 流速下,立管橫向呈現(xiàn)出2 階振動(dòng)模態(tài).高流速下激發(fā)出了更高階的振動(dòng)模態(tài),在0.84 m/s 流速作用下,立管橫向振動(dòng)呈現(xiàn)出4 階模態(tài).進(jìn)一步對(duì)比橫向振動(dòng)位移的均方根值,如圖4(b)和4(d)所示,圖中同時(shí)給出了文獻(xiàn)[32,48-49]的數(shù)值模擬結(jié)果.結(jié)果表明:數(shù)值預(yù)報(bào)的振動(dòng)模態(tài)階數(shù)與試驗(yàn)和其他數(shù)值模擬結(jié)果均是一致的,但振動(dòng)幅值相比試驗(yàn)結(jié)果是偏小的,目前算法在準(zhǔn)確預(yù)報(bào)柔性立管振動(dòng)幅值方面仍有待改進(jìn).本文數(shù)值模型預(yù)報(bào)立管振幅偏小的原因可能是多方面的,首先,DVMFEM 耦合算法基于切片法的思想,將立管簡(jiǎn)化為梁模型進(jìn)行結(jié)構(gòu)動(dòng)力響應(yīng)求解,忽略了立管本身三維效應(yīng)的影響以及尾流中三維渦結(jié)構(gòu)的影響,Huang等[48]和Wang 和Xiao[49]則是采用完全三維網(wǎng)格的方法求解流場(chǎng)和立管變形;其次,Wu 等[50]的研究表明立管流向和橫向運(yùn)動(dòng)之間存在相互作用,本文采用的耦合算法則是單獨(dú)求解每個(gè)方向的立管結(jié)構(gòu)動(dòng)力響應(yīng);此外,離散渦方法本身及其與有限元方法的耦合方式仍有待進(jìn)一步的研究與改進(jìn).

    圖4 不同流速下的柔性立管橫向振動(dòng)響應(yīng)((a),(c)振動(dòng)位移包絡(luò)線;(b),(d)位移均方根值)Fig.4 Vibration response envelopes of a flexible riser under different current speeds ((a),(c) vibration amplitude envelopes;(b),(d) RMS displacements)

    3 懸臂式立管計(jì)算結(jié)果與分析

    本節(jié)主要關(guān)注不同流速下的懸臂式立管的渦激振動(dòng)響應(yīng)特征,同樣選用上述9.63 m 長(zhǎng)的立管,其底部承受的向下的拉力設(shè)定為Wb=817 N,這與原來(lái)兩端鉸支的立管試驗(yàn)中所選擇的頂部張力在數(shù)值上是相同的,其他物理參數(shù)保持不變.因此,不同流速作用下的柔性立管的固有頻率是一致的.基于子空間迭代法可獲得立管濕模態(tài)下的前8 階固有頻率及其對(duì)應(yīng)的歸一化的振型,如圖5 所示.立管的固有頻率隨模態(tài)階數(shù)的增加而增大,本文以1 階固有頻率f1=0.78 Hz 為參考值,取折合速度Ur從4 到46,間隔為2,計(jì)算得到的實(shí)際流速區(qū)間為u∞=0.0626~0.6884,對(duì)應(yīng)的雷諾數(shù)區(qū)間為Re=1252~ 13768.折合速度計(jì)算公式為

    圖5 懸臂式立管前8 階固有頻率和振型Fig.5 First eight-order natural frequencies and modal shapes for a cantilever riser

    雷諾數(shù)計(jì)算公式為

    基于結(jié)構(gòu)物繞流的升力時(shí)歷曲線,采用快速傅里葉變換(FFT),可得到泄渦頻率fs,進(jìn)一步可計(jì)算斯特勞哈爾數(shù)

    3.1 振動(dòng)響應(yīng)分析

    首先,本文計(jì)算了不同折合速度下的柔性懸臂式立管的橫向振動(dòng),振幅均方根值隨Ur和z/L的變化情況如圖6 所示.可以看出,隨著折合速度的增加,立管振動(dòng)模態(tài)逐漸增加,在較高折合速度下高階模態(tài)被激發(fā).在Ur=4~ 46 范圍內(nèi),橫向振動(dòng)逐漸由1 階模態(tài)向4 階模態(tài)過(guò)渡.立管振幅則呈現(xiàn)出波浪式的變化特征,即在一定的折合速度范圍內(nèi),立管的振動(dòng)模態(tài)是相同的,這與Fan 等[35]對(duì)兩端鉸支柔性立管VIV 的研究結(jié)果類(lèi)似.當(dāng)立管呈現(xiàn)出相同模態(tài)時(shí),隨著Ur的增加,立管振幅逐漸增加,同時(shí)在少數(shù)折合速度下振幅沒(méi)有呈現(xiàn)出線性增長(zhǎng)的特征.當(dāng)振動(dòng)從低階模態(tài)向高階模態(tài)轉(zhuǎn)變時(shí),模態(tài)階數(shù)的躍遷過(guò)程中會(huì)伴隨著振幅的突然降低.此外,可以很明顯地觀察到無(wú)約束的立管底部表現(xiàn)出較大幅度的位移.

    圖6 不同折合速度下的橫向振幅均方根值Fig.6 RMS vibration amplitudes in the transverse direction under a wide range of reduced velocities

    針對(duì)不同振動(dòng)模態(tài)下立管的變形特征,本文分別將1 階到4 階模態(tài)下立管橫向振動(dòng)的位移均方根值進(jìn)行了比較分析,如圖7(a)到圖7(d)所示.并且列舉了四種振動(dòng)模態(tài)下某一折合速度對(duì)應(yīng)的立管瞬態(tài)振動(dòng)包絡(luò)線圖,如圖7(e)到圖7(h)所示.同樣可以看出,大部分情況下,沿立管展向各位置處的橫向位移隨折合速度增加而逐漸增大,該變化特征在3 階模態(tài)情況下最為明顯,如圖7(c)所示.此外,在同一振動(dòng)模態(tài)下,振幅最小處所對(duì)應(yīng)的立管展向位置基本上是相同的,與來(lái)流速度關(guān)系較小.

    圖7 不同折合速度下的柔性立管橫向振動(dòng)響應(yīng) ((a)~(d)位移均方根值;(e)~(h)振動(dòng)位移包絡(luò)線)Fig.7 Vibration response envelopes of a flexible riser under different reduced velocities ((a)~(d) RMS amplitudes;(e)~(h) vibration amplitude envelopes)

    在深海采礦系統(tǒng)中,懸臂式立管底部與中間倉(cāng)相連接,礦物經(jīng)輸料軟管到達(dá)中間倉(cāng),再通過(guò)垂直提升硬管輸送到水面母船.因此,立管底部的振動(dòng)對(duì)礦物輸運(yùn)的穩(wěn)定性有較大的影響.本文將中間倉(cāng)簡(jiǎn)化為有一個(gè)作用在立管底部的垂直向下的作用力進(jìn)行研究,如圖2 所示.由上述分析可知,懸臂式立管底部位置有較大幅度的位移,不同折合速度下的立管底部位移均方根值如圖8 所示.可以看出,不同模態(tài)下振幅的增長(zhǎng)速率是不同的.在同一振動(dòng)模態(tài)下,立管底部振幅與折合速度基本上呈現(xiàn)出線性關(guān)系.

    圖8 立管底部橫向位移均方根值Fig.8 RMS amplitude of the transverse displacement at the bottom for the riser

    不同折合速度下立管振動(dòng)的主導(dǎo)頻率fv如圖9 所示.圖中St=0.2 這條線表示當(dāng)圓柱在亞臨界雷諾數(shù)范圍內(nèi)振動(dòng)時(shí),其對(duì)應(yīng)的斯特勞哈爾數(shù)位于0.2 附近[51].整體上,振動(dòng)頻率隨著折合速度的增加而逐漸增大.當(dāng)振動(dòng)模態(tài)發(fā)生轉(zhuǎn)換時(shí),振動(dòng)頻率會(huì)發(fā)生突然的跳躍.在同一振動(dòng)模態(tài)內(nèi),振動(dòng)主導(dǎo)頻率在低折合速度下沿著St=0.2 增加,且頻率值低于所對(duì)應(yīng)振動(dòng)模態(tài)下的固有頻率.隨著折合速度進(jìn)一步增加,當(dāng)振動(dòng)頻率超過(guò)固有頻率后,主導(dǎo)頻率逐漸偏離S t=0.2 所對(duì)應(yīng)的頻率值,該現(xiàn)象與剛性圓柱渦激振動(dòng)是類(lèi)似的[52].

    圖9 橫向振動(dòng)的主導(dǎo)頻率比Fig.9 Dominant frequency ratio of the transverse vibration

    3.2 三階模態(tài)下橫向振動(dòng)響應(yīng)分析

    由上文分析可知,懸臂式立管發(fā)生渦激振動(dòng)時(shí),其振動(dòng)模態(tài)在一定折合速度范圍內(nèi)是相同的,因此本節(jié)選取典型的3 階振動(dòng)模態(tài)進(jìn)行討論.本文中立管以3 階模態(tài)振動(dòng)時(shí)所對(duì)應(yīng)的折合速度Ur=24~36,間隔為2,共7 個(gè)工況,所對(duì)應(yīng)的橫向振幅均方根值如圖7(c)所示.基于立管垂向不同位置振動(dòng)時(shí)歷曲線,采用快速傅里葉變換(FFT)方法對(duì)時(shí)歷數(shù)據(jù)進(jìn)行頻譜分析,可得到振動(dòng)頻率的空間分布特征,如圖10 所示.其中顏色條為歸一化的振動(dòng)能量譜密度(PSD),白虛線表示3 階振型所對(duì)應(yīng)的固有頻率,f3=4.07 Hz.可以看出,7 種折合速度對(duì)應(yīng)的振動(dòng)頻率均由3 階模態(tài)主導(dǎo),且不同折合速度下立管底部均呈現(xiàn)出較高能量的振動(dòng).此外,隨著折合速度增加,振動(dòng)頻率逐漸遠(yuǎn)離3 階固有頻率.

    同樣采用FFT 方法,可得到不同折合速度下的立管展向各位置處的泄渦頻率空間分布,如圖11 所示,其中白虛線均表示St=0.2 所對(duì)應(yīng)的頻率值.與圖10 相比,不同折合速度下泄渦頻率與振動(dòng)頻率基本上是一致的.在低折合速度下,立管各位置處的泄渦頻率均集中在St=0.2 附近,隨著速度增加,立管展向大部分位置所對(duì)應(yīng)的泄渦頻率逐漸低于S t=0.2所對(duì)應(yīng)的頻率值,盡管仍有較少位置的泄渦頻率集中在St=0.2 附近.不同折合速度下的泄渦頻率在立管展長(zhǎng)方向上并未呈現(xiàn)出明顯的分布特征,這與振動(dòng)頻率空間分布圖中呈現(xiàn)出的3 階分布特征是不同的.

    圖10 立管橫向振動(dòng)頻率空間分布 ((a)~(g) Ur=24~36)Fig.10 Spatial distribution of transverse vibration frequency for the riser ((a)~(g) Ur=24~36)

    圖11 立管泄渦頻率空間分布 ((a)~(g) Ur=24~36)Fig.11 Spatial distribution of vortex shedding frequency for the riser ((a)~(g) Ur=24~36)

    圖11 立管泄渦頻率空間分布 ((a)~(g) Ur=24~36) (續(xù))Fig.11 Spatial distribution of vortex shedding frequency for the riser ((a)~(g) Ur=24~36) (continued)

    3 階振動(dòng)模態(tài)下不同折合速度所對(duì)應(yīng)的立管橫向振幅時(shí)空演化特征如圖12 所示,所有時(shí)空演化圖選擇的時(shí)間段均為6~8 s,其中紅色和藍(lán)色表示相反的振動(dòng)方向.可以看出,隨著折合速度的增加,立管展向不同位置處的振幅增加,且在該時(shí)間段內(nèi)駐波特征逐漸增強(qiáng).由上文分析可知,同一振動(dòng)模態(tài)下,振動(dòng)頻率隨折合速度的增加而增大,因此相同時(shí)間段內(nèi)紅藍(lán)單元交替的周期逐漸減小.

    圖12 立管不同位置的振幅時(shí)空演化 ((a)~(g) Ur=24~36)Fig.12 Temporal-spatial evolution of transverse vibration amplitude along the riser span ((a)~(g) Ur=24~36)

    本文工作基于切片法的思想,在切片上采用無(wú)網(wǎng)格的二維離散渦算法計(jì)算流體力,四種折合速度下(Ur=24,28,32,36)的立管若干截面處的瞬態(tài)渦量場(chǎng)如圖13 所示.可以看出,隨著折合速度的增加,立管底部的流向位移逐漸增大.立管不同高度截面渦量場(chǎng)中紅色和藍(lán)色粒子分別表示渦元具有正向和負(fù)向的渦量.立管不同截面處的尾流場(chǎng)中的渦結(jié)構(gòu)具有不同的模式[53],其中大部分截面呈現(xiàn)出2 S 的尾渦模式.此外,近尾流場(chǎng)渦結(jié)構(gòu)比較清晰,而遠(yuǎn)尾流場(chǎng)則演化出更為復(fù)雜的渦結(jié)構(gòu)特征.

    圖13 瞬態(tài)渦量場(chǎng) ((a)~(d) Ur=24~30)Fig.13 Instantaneous vorticity field ((a)~(d) Ur=24~36)

    3.3 兩端鉸支與懸臂式立管的渦激振動(dòng)特征比較

    為了比較傳統(tǒng)的兩端鉸支立管(對(duì)應(yīng)海洋工程中的采油管道)與懸臂式立管(對(duì)應(yīng)深海采礦的提升管道)渦激振動(dòng)特征的異同,本文計(jì)算了一系列均勻流速作用下的兩端鉸支的9.63 m 長(zhǎng)立管的渦激振動(dòng)響應(yīng).立管頂部張力與試驗(yàn)保持一致,取為T(mén)top=817 N.立管濕模態(tài)下的前8 階固有頻率及其對(duì)應(yīng)的歸一化的振型,如圖14 所示.可以看出,兩端鉸支立管的各階固有頻率均大于對(duì)應(yīng)的底部拉力為817 N 的懸臂式立管的固有頻率.兩端鉸支立管的1 階固有頻率f1=1.61 Hz,本次計(jì)算取折合速度Ur=4~ 44,間隔為2.計(jì)算得到的不同來(lái)流速度作用下立管最大的橫向位移均方根值如圖15 所示,圖中不同顏色塊表示不同的振動(dòng)模態(tài).可以發(fā)現(xiàn),橫向振幅的變化規(guī)律與懸臂式立管是一致的,即在同一振動(dòng)模態(tài)下,立管振幅與折合速度基本上呈現(xiàn)出線性關(guān)系.此外,當(dāng)模態(tài)轉(zhuǎn)變發(fā)生時(shí),振幅會(huì)突然降低,并隨著折合速度增加,振幅再一次逐漸增大,這與懸臂式立管的VIV 特征也是相同的.

    圖14 兩端鉸支立管前8 階固有頻率和振型Fig.14 First eight-order natural frequencies and modal shapes for a riser hinged at both ends

    圖15 沿立管展向最大的橫向位移均方根值Fig.15 Maximum RMS amplitude of the transverse displacement along the riser span

    橫向振動(dòng)的主導(dǎo)頻率與1 階固有頻率的比值如圖16 所示.同樣地,振動(dòng)頻率的變化趨勢(shì)與懸臂式立管也是一致的,即振動(dòng)頻率隨著折合速度的增加而逐漸增大,并且當(dāng)振動(dòng)模態(tài)改變時(shí),振動(dòng)頻率會(huì)出現(xiàn)較大幅度的增加.不同的是,在相同的折合速度區(qū)間內(nèi),兩端鉸支立管存在6 階振動(dòng)模態(tài),而懸臂式立管只存在4 階模態(tài).這是由于兩端鉸支立管的1 階固有頻率較大,相同折合速度下其對(duì)應(yīng)的真實(shí)來(lái)流速度更大,這導(dǎo)致兩端鉸支立管被激發(fā)出高階模態(tài)響應(yīng).本文選取兩端鉸支立管的3 階振動(dòng)模態(tài)進(jìn)行分析,該模態(tài)下3 個(gè)不同折合速度作用下的立管不同位置的振幅時(shí)空演化如圖17 所示,選擇的時(shí)間段均為5~7 s.可以發(fā)現(xiàn),兩端鉸支的立管展向不同位置處的振幅均隨折合速度增加而增大.在低階模態(tài)響應(yīng)下,隨折合速度增加,該時(shí)間段內(nèi)駐波特征逐漸增強(qiáng).這兩個(gè)特征與懸臂式立管是一致的.

    圖16 立管橫向振動(dòng)的主導(dǎo)頻率比Fig.16 Dominant frequency ratio of the transverse vibration for the riser

    圖17 立管不同位置的振幅時(shí)空演化 ((a)~(c) Ur=16~20)Fig.17 Temporal-spatial evolution of transverse vibration amplitude along the riser span ((a)~(c) Ur=16~20)

    4 結(jié)論

    礦產(chǎn)混輸垂直提升管道是深海采礦集成系統(tǒng)中十分重要的裝備,其懸臂式的結(jié)構(gòu)面臨復(fù)雜海洋環(huán)境時(shí)的振動(dòng)特征不同于傳統(tǒng)海洋立管,目前文獻(xiàn)仍缺乏對(duì)該類(lèi)型立管的討論.本文的主要目的便是分析具有一定底部重量的懸臂式立管在不同來(lái)流速度下的渦激振動(dòng)響應(yīng).為此,本文首先建立了基于切片法思想的DVM-FEM 耦合模型,通過(guò)數(shù)值計(jì)算兩端鉸支立管的振動(dòng)響應(yīng)并與模型試驗(yàn)結(jié)果對(duì)比,驗(yàn)證了該耦合模型的有效性.最后,系統(tǒng)性地分析了不同速度的均勻來(lái)流作用下懸臂式立管的橫向振動(dòng)響應(yīng)特征,主要結(jié)論如下.

    (1)本文采用的DVM-FEM 耦合算法可以準(zhǔn)確預(yù)報(bào)柔性立管在不同速度下的振動(dòng)模態(tài),但振幅相比試驗(yàn)結(jié)果要小,該問(wèn)題同樣出現(xiàn)在前人的研究工作中[39].由于本文目的是研究深海采礦系統(tǒng)中存在的懸臂式立管的渦激振動(dòng)特征,因此目前工作主要分析不同來(lái)流速度對(duì)這類(lèi)特殊立管的影響規(guī)律,如何實(shí)現(xiàn)對(duì)立管實(shí)際振幅的準(zhǔn)確預(yù)報(bào)仍有待進(jìn)一步的研究.

    (2)對(duì)于懸臂式立管,隨著折合速度的增加,立管的振動(dòng)模態(tài)階數(shù)逐漸增加.同時(shí),振動(dòng)模態(tài)在一定的折合速度范圍內(nèi)是相同的.在同一模態(tài)下,立管不同位置的振幅隨著速度增加而增大,尤其是無(wú)約束的立管底部會(huì)發(fā)生較大的位移.此外,當(dāng)模態(tài)發(fā)生轉(zhuǎn)換時(shí),振幅會(huì)發(fā)生突降,并隨著速度的增加振幅再一次逐漸增大.振幅隨Ur和z/L的改變整體上呈現(xiàn)出波浪式的變化特征.

    (3)深海采礦垂直提升立管底部與中間倉(cāng)及輸送軟管相連接,因此底部位移對(duì)采礦系統(tǒng)整體穩(wěn)定性有重要影響.本文結(jié)果表明,在同一振動(dòng)模態(tài)下,懸臂式立管底部位移呈現(xiàn)出隨折合速度線性變化的特征.當(dāng)振動(dòng)模態(tài)發(fā)生轉(zhuǎn)變時(shí),振動(dòng)主導(dǎo)頻率存在突然跳躍的特征.

    (4)本文重點(diǎn)分析了3 階振動(dòng)模態(tài)下的懸臂式立管的振動(dòng)響應(yīng).可以看出,振動(dòng)主導(dǎo)頻率在低折合速度下低于3 階固有頻率,隨著速度增加,振動(dòng)頻率逐漸增大并遠(yuǎn)離3 階固有頻率.此外,無(wú)約束的立管底部的振動(dòng)能量較大.對(duì)于本次研究中3 階模態(tài)振動(dòng)下的懸臂式立管,其駐波特征隨著折合速度增加而逐漸加強(qiáng),需要注意的是,本文關(guān)注的立管模型的長(zhǎng)徑比較小且為低階模態(tài)響應(yīng),對(duì)于實(shí)尺度深海采礦立管以及高階振動(dòng)模態(tài)響應(yīng)下的立管,其振動(dòng)特征仍需進(jìn)一步研究.

    (5)本文對(duì)比了懸臂式立管與兩端鉸支立管的異同.可以發(fā)現(xiàn),懸臂式立管的各階固有頻率均小于對(duì)應(yīng)的兩端鉸支立管的固有頻率.在相同折合速度區(qū)間內(nèi),兩端鉸支立管被激發(fā)出更高階的振動(dòng)模態(tài)響應(yīng).此外,兩種類(lèi)型立管的VIV 振幅和頻率隨折合速度的變化規(guī)律是相同的.

    (6)垂直管道水力提升系統(tǒng)是深海采礦的主流方式,但目前仍缺乏對(duì)這類(lèi)特殊的底部弱約束立管的數(shù)值與試驗(yàn)研究工作.未來(lái)研究工作可以考慮以下幾方面:母船運(yùn)動(dòng)對(duì)懸臂式立管渦激振動(dòng)和水動(dòng)力性能的影響;中間倉(cāng)重量改變的影響;考慮中間倉(cāng)慣性力作用下的深海采礦立管模型的結(jié)構(gòu)動(dòng)力響應(yīng);海洋剪切流動(dòng)等變化的影響;考慮內(nèi)部流動(dòng)的立管振動(dòng)特征;采礦立管渦激振動(dòng)抑制措施等.

    猜你喜歡
    渦激立管懸臂
    不同間距比下串聯(lián)圓柱渦激振動(dòng)數(shù)值模擬研究
    渦激振動(dòng)發(fā)電裝置及其關(guān)鍵技術(shù)
    常見(jiàn)高層建筑物室內(nèi)給水立管材質(zhì)解析
    懸臂式硫化罐的開(kāi)發(fā)設(shè)計(jì)
    盤(pán)球立管結(jié)構(gòu)抑制渦激振動(dòng)的數(shù)值分析方法研究
    電子制作(2018年14期)2018-08-21 01:38:42
    當(dāng)液壓遇上懸臂云臺(tái) 捷信GHFG1液壓懸臂云臺(tái)試用
    探討掛籃懸臂灌注連梁的施工
    LF爐懸臂爐蓋防傾覆結(jié)構(gòu)
    深水鋼懸鏈立管J型鋪設(shè)研究
    柔性圓管在渦激振動(dòng)下的模態(tài)響應(yīng)分析
    蜜臀久久99精品久久宅男| 国产精品久久久久成人av| 国产视频内射| 在线播放无遮挡| 纵有疾风起免费观看全集完整版| 亚洲av.av天堂| 久久青草综合色| 少妇 在线观看| 青春草视频在线免费观看| 一级毛片黄色毛片免费观看视频| 女人精品久久久久毛片| 久久国内精品自在自线图片| 一个人免费看片子| 中文字幕人妻丝袜制服| 五月玫瑰六月丁香| 99久久精品国产国产毛片| 在线亚洲精品国产二区图片欧美 | 国产精品成人在线| 97在线视频观看| 嫩草影院入口| videos熟女内射| 成人免费观看视频高清| 亚洲怡红院男人天堂| 国产成人av激情在线播放 | 国产av码专区亚洲av| 久久精品久久久久久噜噜老黄| 黄色怎么调成土黄色| 国产乱来视频区| 两个人免费观看高清视频| 成人毛片60女人毛片免费| 一本色道久久久久久精品综合| 成人无遮挡网站| 国产深夜福利视频在线观看| 国产av精品麻豆| 搡女人真爽免费视频火全软件| 欧美xxⅹ黑人| 精品一区二区免费观看| 99久久精品国产国产毛片| 老女人水多毛片| 黑丝袜美女国产一区| 人人妻人人澡人人看| 免费观看性生交大片5| 一区二区av电影网| 一二三四中文在线观看免费高清| 一二三四中文在线观看免费高清| 丝袜美足系列| 大香蕉97超碰在线| 最近的中文字幕免费完整| 99热这里只有精品一区| 国产淫语在线视频| 日韩熟女老妇一区二区性免费视频| av网站免费在线观看视频| 少妇人妻 视频| 一级a做视频免费观看| a级毛片黄视频| 国产一区有黄有色的免费视频| 一级毛片黄色毛片免费观看视频| 亚洲怡红院男人天堂| 国产高清有码在线观看视频| 飞空精品影院首页| 国产高清不卡午夜福利| 亚洲国产精品专区欧美| 老司机影院成人| 青春草视频在线免费观看| 亚洲av欧美aⅴ国产| 国产男女内射视频| 丰满饥渴人妻一区二区三| 尾随美女入室| 91aial.com中文字幕在线观看| 黄片播放在线免费| 亚洲精品乱码久久久v下载方式| 91精品一卡2卡3卡4卡| av专区在线播放| 亚洲av在线观看美女高潮| 麻豆乱淫一区二区| 亚洲精品一区蜜桃| 80岁老熟妇乱子伦牲交| 国产男女超爽视频在线观看| 亚洲一区二区三区欧美精品| 精品人妻偷拍中文字幕| 校园人妻丝袜中文字幕| 国产乱人偷精品视频| 色网站视频免费| 桃花免费在线播放| a 毛片基地| 99久久精品一区二区三区| 久久久亚洲精品成人影院| 亚洲精品国产色婷婷电影| 久久精品国产亚洲网站| 国产一区二区三区综合在线观看 | 狂野欧美激情性xxxx在线观看| 麻豆乱淫一区二区| 99re6热这里在线精品视频| 午夜av观看不卡| 久久久久久久精品精品| av黄色大香蕉| 久久久久久久久久久久大奶| 人人妻人人添人人爽欧美一区卜| 亚洲国产成人一精品久久久| videosex国产| 一本色道久久久久久精品综合| 国产亚洲精品久久久com| 精品一区二区三卡| 中文天堂在线官网| 亚洲精品国产色婷婷电影| 热re99久久精品国产66热6| 热99久久久久精品小说推荐| 黄片无遮挡物在线观看| 能在线免费看毛片的网站| 高清在线视频一区二区三区| 久久人妻熟女aⅴ| 又粗又硬又长又爽又黄的视频| av女优亚洲男人天堂| 秋霞伦理黄片| 久久精品国产自在天天线| 亚洲婷婷狠狠爱综合网| 精品一品国产午夜福利视频| 91国产中文字幕| 视频中文字幕在线观看| 精品亚洲乱码少妇综合久久| 纵有疾风起免费观看全集完整版| xxxhd国产人妻xxx| av在线老鸭窝| 男女边摸边吃奶| 99九九线精品视频在线观看视频| 欧美 日韩 精品 国产| 黄色一级大片看看| 一区二区av电影网| 狂野欧美激情性bbbbbb| 免费看光身美女| 极品人妻少妇av视频| 亚洲四区av| 80岁老熟妇乱子伦牲交| 国产黄频视频在线观看| 久久国内精品自在自线图片| 久久久久久久久久成人| 少妇人妻久久综合中文| 亚洲国产精品成人久久小说| 久久鲁丝午夜福利片| 国产一区二区在线观看av| av卡一久久| 亚洲欧美色中文字幕在线| 亚洲第一区二区三区不卡| 亚洲熟女精品中文字幕| 春色校园在线视频观看| 91成人精品电影| 少妇高潮的动态图| 国产国拍精品亚洲av在线观看| 亚洲精品成人av观看孕妇| 日本黄色日本黄色录像| 一本色道久久久久久精品综合| 日韩 亚洲 欧美在线| 午夜影院在线不卡| 日本av免费视频播放| xxxhd国产人妻xxx| 极品少妇高潮喷水抽搐| av网站免费在线观看视频| 一级二级三级毛片免费看| av线在线观看网站| 婷婷色麻豆天堂久久| 国产高清有码在线观看视频| 国产免费又黄又爽又色| 黄片无遮挡物在线观看| 午夜老司机福利剧场| 秋霞伦理黄片| 亚洲伊人久久精品综合| 国产精品久久久久久精品电影小说| 熟女av电影| 日韩中字成人| a 毛片基地| 久久鲁丝午夜福利片| 黑人高潮一二区| av女优亚洲男人天堂| 99国产精品免费福利视频| 精品一区二区免费观看| 亚洲精品成人av观看孕妇| 日日摸夜夜添夜夜添av毛片| 久久99蜜桃精品久久| 国模一区二区三区四区视频| 制服人妻中文乱码| 女性被躁到高潮视频| 在线观看人妻少妇| 国产成人精品在线电影| 一区在线观看完整版| 日韩成人av中文字幕在线观看| 国产69精品久久久久777片| 欧美人与善性xxx| 日本色播在线视频| 久久精品久久精品一区二区三区| 国产一级毛片在线| 91精品三级在线观看| 欧美人与性动交α欧美精品济南到 | 欧美 亚洲 国产 日韩一| 久久久国产欧美日韩av| 国产av码专区亚洲av| 各种免费的搞黄视频| 亚洲精品国产av成人精品| av在线app专区| 国产av精品麻豆| 大码成人一级视频| 人妻少妇偷人精品九色| 久久久欧美国产精品| 免费播放大片免费观看视频在线观看| av女优亚洲男人天堂| 一个人看视频在线观看www免费| 爱豆传媒免费全集在线观看| 午夜精品国产一区二区电影| 亚洲欧美一区二区三区黑人 | 亚洲国产欧美日韩在线播放| 一级毛片aaaaaa免费看小| 黄片无遮挡物在线观看| 亚洲四区av| 国产成人freesex在线| 亚洲av福利一区| 亚洲精品中文字幕在线视频| 91久久精品电影网| 99九九在线精品视频| 蜜桃在线观看..| 一级毛片 在线播放| 青春草国产在线视频| 人体艺术视频欧美日本| 日日爽夜夜爽网站| 国产精品不卡视频一区二区| 久久久国产一区二区| 久久久久久久国产电影| 亚洲精品久久久久久婷婷小说| 国模一区二区三区四区视频| 亚洲性久久影院| 美女国产视频在线观看| 国产69精品久久久久777片| 狠狠婷婷综合久久久久久88av| 久久综合国产亚洲精品| 日韩一区二区视频免费看| 国产一区二区在线观看日韩| 国产日韩欧美视频二区| 熟女电影av网| 黄色毛片三级朝国网站| 欧美+日韩+精品| 国产欧美另类精品又又久久亚洲欧美| 99久久人妻综合| 18在线观看网站| 一本大道久久a久久精品| 不卡视频在线观看欧美| 国产日韩一区二区三区精品不卡 | 日本色播在线视频| a级毛色黄片| 制服丝袜香蕉在线| 最后的刺客免费高清国语| 亚洲国产精品999| 亚洲精品色激情综合| 欧美成人精品欧美一级黄| 精品人妻偷拍中文字幕| 五月天丁香电影| 成人免费观看视频高清| 久久精品久久精品一区二区三区| 尾随美女入室| 日韩强制内射视频| 久久久精品免费免费高清| 一级毛片 在线播放| 嘟嘟电影网在线观看| 纵有疾风起免费观看全集完整版| 国产精品久久久久成人av| 岛国毛片在线播放| 亚洲精品久久午夜乱码| 亚洲欧美成人精品一区二区| 国产精品三级大全| 日韩视频在线欧美| 日韩中文字幕视频在线看片| 精品亚洲乱码少妇综合久久| 少妇人妻 视频| 丰满饥渴人妻一区二区三| 国产免费一级a男人的天堂| 久久韩国三级中文字幕| 看免费成人av毛片| 亚洲精品久久午夜乱码| 久久久久久久久久成人| 国产精品女同一区二区软件| 另类亚洲欧美激情| 免费久久久久久久精品成人欧美视频 | 少妇熟女欧美另类| 久久久亚洲精品成人影院| 夜夜骑夜夜射夜夜干| 欧美日本中文国产一区发布| 黑人欧美特级aaaaaa片| 另类亚洲欧美激情| 国产男人的电影天堂91| 午夜福利视频在线观看免费| 免费播放大片免费观看视频在线观看| 亚洲国产成人一精品久久久| 69精品国产乱码久久久| 国产无遮挡羞羞视频在线观看| 中文天堂在线官网| a级毛片黄视频| 插逼视频在线观看| av电影中文网址| 黄色一级大片看看| 日日摸夜夜添夜夜爱| 91精品一卡2卡3卡4卡| 高清av免费在线| 午夜福利视频在线观看免费| 女人精品久久久久毛片| 中文字幕人妻熟人妻熟丝袜美| 欧美bdsm另类| 日本黄色片子视频| 亚洲精品,欧美精品| 国产成人aa在线观看| 国产有黄有色有爽视频| 一本久久精品| 日本免费在线观看一区| 99九九在线精品视频| 少妇高潮的动态图| 91精品国产国语对白视频| 99re6热这里在线精品视频| 国产高清国产精品国产三级| 赤兔流量卡办理| 啦啦啦在线观看免费高清www| 高清毛片免费看| 亚洲精品日韩av片在线观看| 成人二区视频| 制服诱惑二区| 亚洲国产欧美日韩在线播放| 69精品国产乱码久久久| 999精品在线视频| 午夜日本视频在线| 亚洲精品成人av观看孕妇| 久久久久久久久久久久大奶| 五月天丁香电影| 99视频精品全部免费 在线| 久久97久久精品| 色5月婷婷丁香| 亚洲欧洲日产国产| 日韩一区二区视频免费看| 一级片'在线观看视频| 欧美精品高潮呻吟av久久| 99热网站在线观看| 丁香六月天网| 99久久中文字幕三级久久日本| 亚洲国产欧美日韩在线播放| 成人黄色视频免费在线看| 日本黄色片子视频| 男人操女人黄网站| 爱豆传媒免费全集在线观看| tube8黄色片| 亚洲国产精品专区欧美| 老女人水多毛片| 一区二区av电影网| 丝瓜视频免费看黄片| 最黄视频免费看| 另类精品久久| 3wmmmm亚洲av在线观看| 久热久热在线精品观看| 制服诱惑二区| 三上悠亚av全集在线观看| 一区二区三区精品91| 少妇人妻精品综合一区二区| 久久99精品国语久久久| 久久久精品区二区三区| 成人手机av| 麻豆乱淫一区二区| 日本与韩国留学比较| 久久人人爽人人片av| 满18在线观看网站| 2021少妇久久久久久久久久久| 日本免费在线观看一区| 国产片特级美女逼逼视频| 永久网站在线| 黄色视频在线播放观看不卡| 精品人妻偷拍中文字幕| 人妻系列 视频| 这个男人来自地球电影免费观看 | 成人国产麻豆网| 国产一级毛片在线| 欧美bdsm另类| 高清在线视频一区二区三区| 最近最新中文字幕免费大全7| 高清在线视频一区二区三区| 精品久久蜜臀av无| 日韩电影二区| 国产精品欧美亚洲77777| h视频一区二区三区| 九色成人免费人妻av| 亚洲三级黄色毛片| av一本久久久久| av专区在线播放| 欧美少妇被猛烈插入视频| 久久久久久久亚洲中文字幕| 国产精品无大码| 热99国产精品久久久久久7| av专区在线播放| 婷婷色麻豆天堂久久| 亚洲色图综合在线观看| 国产亚洲精品第一综合不卡 | av免费观看日本| 久久人人爽人人片av| 永久免费av网站大全| 校园人妻丝袜中文字幕| 不卡视频在线观看欧美| 卡戴珊不雅视频在线播放| 黄色配什么色好看| av卡一久久| 九九在线视频观看精品| 一本一本综合久久| 22中文网久久字幕| 麻豆精品久久久久久蜜桃| 久久99蜜桃精品久久| 欧美3d第一页| 国产精品秋霞免费鲁丝片| 亚洲人成网站在线播| 人人妻人人澡人人看| 极品少妇高潮喷水抽搐| 人人妻人人澡人人看| 亚洲精品第二区| 亚洲av国产av综合av卡| 欧美3d第一页| 男女边摸边吃奶| 欧美 亚洲 国产 日韩一| 日本黄大片高清| 久久久久精品性色| av在线观看视频网站免费| 热re99久久精品国产66热6| 日韩中字成人| 日日啪夜夜爽| 精品久久蜜臀av无| 国产又色又爽无遮挡免| 久久97久久精品| 91精品一卡2卡3卡4卡| 色哟哟·www| 亚洲国产色片| 一级二级三级毛片免费看| 蜜桃在线观看..| 国产男人的电影天堂91| 乱码一卡2卡4卡精品| 两个人的视频大全免费| 插阴视频在线观看视频| videosex国产| 啦啦啦中文免费视频观看日本| 少妇的逼水好多| 日韩精品有码人妻一区| 热re99久久精品国产66热6| 欧美亚洲 丝袜 人妻 在线| 一本一本综合久久| 十八禁网站网址无遮挡| 国产一区亚洲一区在线观看| 久久精品国产亚洲av天美| 天天躁夜夜躁狠狠久久av| 亚洲色图综合在线观看| 精品人妻在线不人妻| 最黄视频免费看| 久久精品久久久久久久性| 成人毛片60女人毛片免费| 精品人妻熟女av久视频| 中文欧美无线码| 亚洲精品国产av蜜桃| 一本色道久久久久久精品综合| av专区在线播放| 国产日韩一区二区三区精品不卡 | 少妇人妻久久综合中文| 国产熟女欧美一区二区| 婷婷色综合www| 国产永久视频网站| 超碰97精品在线观看| 又大又黄又爽视频免费| 中文字幕久久专区| 天堂8中文在线网| 精品酒店卫生间| 婷婷色综合www| 日韩av在线免费看完整版不卡| 天堂俺去俺来也www色官网| 精品国产乱码久久久久久小说| 免费黄频网站在线观看国产| 国产精品欧美亚洲77777| 亚洲在久久综合| 成人无遮挡网站| 亚洲av欧美aⅴ国产| 性色av一级| 少妇人妻 视频| 国产色婷婷99| 丰满迷人的少妇在线观看| 免费看光身美女| 日韩免费高清中文字幕av| 久久久久久久久久久久大奶| 国产精品蜜桃在线观看| 亚洲一级一片aⅴ在线观看| 99久久精品一区二区三区| 日本av手机在线免费观看| 免费高清在线观看日韩| 国产老妇伦熟女老妇高清| 看非洲黑人一级黄片| 久久精品久久久久久久性| 亚洲国产av新网站| 亚洲丝袜综合中文字幕| 一级二级三级毛片免费看| 亚洲欧美成人精品一区二区| 搡老乐熟女国产| 国产成人精品福利久久| 精品久久久精品久久久| 精品国产露脸久久av麻豆| 日韩欧美精品免费久久| 亚洲精品自拍成人| 亚洲精品国产av蜜桃| 又粗又硬又长又爽又黄的视频| 人人妻人人爽人人添夜夜欢视频| videossex国产| 91午夜精品亚洲一区二区三区| 晚上一个人看的免费电影| 日韩中文字幕视频在线看片| 亚洲无线观看免费| 性色avwww在线观看| 亚洲精品乱码久久久久久按摩| 女的被弄到高潮叫床怎么办| xxx大片免费视频| 国产精品人妻久久久久久| 男男h啪啪无遮挡| 制服丝袜香蕉在线| 亚洲精品456在线播放app| 18禁裸乳无遮挡动漫免费视频| 国产精品国产av在线观看| 毛片一级片免费看久久久久| 日韩av不卡免费在线播放| 亚洲国产欧美日韩在线播放| 国产成人精品婷婷| 91精品三级在线观看| 日韩免费高清中文字幕av| 蜜臀久久99精品久久宅男| 亚洲国产精品国产精品| 一级毛片电影观看| 免费高清在线观看日韩| 免费黄色在线免费观看| 免费播放大片免费观看视频在线观看| 99久久精品一区二区三区| 日本爱情动作片www.在线观看| 日本-黄色视频高清免费观看| a级毛色黄片| 国产成人91sexporn| 热99久久久久精品小说推荐| 午夜91福利影院| 日韩一本色道免费dvd| 99热这里只有精品一区| 一级毛片电影观看| 亚洲经典国产精华液单| 亚洲五月色婷婷综合| 国产精品秋霞免费鲁丝片| 亚洲三级黄色毛片| 亚洲色图综合在线观看| av网站免费在线观看视频| 成年av动漫网址| 少妇的逼好多水| 女性生殖器流出的白浆| 国产成人aa在线观看| 亚洲国产av新网站| 国产成人午夜福利电影在线观看| 最后的刺客免费高清国语| 国产黄频视频在线观看| av在线老鸭窝| av国产精品久久久久影院| 大话2 男鬼变身卡| 久久久久久久国产电影| √禁漫天堂资源中文www| 国产成人精品在线电影| 三级国产精品片| 女人精品久久久久毛片| 人妻一区二区av| 免费观看的影片在线观看| 哪个播放器可以免费观看大片| 国产一区有黄有色的免费视频| 日本爱情动作片www.在线观看| 国产女主播在线喷水免费视频网站| 欧美日韩亚洲高清精品| 高清av免费在线| 女性生殖器流出的白浆| 一边摸一边做爽爽视频免费| 各种免费的搞黄视频| 亚洲在久久综合| a 毛片基地| 亚洲综合色惰| 高清视频免费观看一区二区| a级毛片黄视频| 草草在线视频免费看| 赤兔流量卡办理| 午夜老司机福利剧场| 欧美日韩在线观看h| 亚洲第一av免费看| 免费观看无遮挡的男女| 婷婷色麻豆天堂久久| 国产亚洲av片在线观看秒播厂| 天美传媒精品一区二区| 中文天堂在线官网| 亚洲图色成人| √禁漫天堂资源中文www| 国产精品一区二区在线观看99| 国产亚洲一区二区精品| 高清av免费在线| 国产高清不卡午夜福利| 日韩 亚洲 欧美在线| 国产熟女午夜一区二区三区 | 日本黄色日本黄色录像| 母亲3免费完整高清在线观看 | 久久鲁丝午夜福利片| 国产精品99久久99久久久不卡 | 日韩成人伦理影院| 99久国产av精品国产电影| 欧美bdsm另类| 久久亚洲国产成人精品v| 久久久久网色| 日韩制服骚丝袜av| 日日摸夜夜添夜夜添av毛片| 精品一区二区三区视频在线| www.色视频.com| 伊人久久精品亚洲午夜| 水蜜桃什么品种好| 亚洲av男天堂| 国产一区二区三区综合在线观看 | 伦精品一区二区三区| 美女国产高潮福利片在线看| 日韩不卡一区二区三区视频在线|