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

    均勻流下柔性立管渦激振動(dòng)響應(yīng)及渦激力載荷特性研究

    2017-11-30 05:49:28宋磊建付世曉于大鵬張萌萌
    振動(dòng)與沖擊 2017年22期
    關(guān)鍵詞:立管幅值柔性

    宋磊建, 付世曉, 任 鐵,3, 于大鵬, 張萌萌

    (1. 上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海 200240; 2. 上海交通大學(xué) 高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心,上海 200240;3. 中國(guó)船舶及海洋工程設(shè)計(jì)研究院,上海 200011; 4. 海軍裝備研究院,北京 100161)

    均勻流下柔性立管渦激振動(dòng)響應(yīng)及渦激力載荷特性研究

    宋磊建1,2, 付世曉1,2, 任 鐵1,2,3, 于大鵬4, 張萌萌1,2

    (1. 上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海 200240; 2. 上海交通大學(xué) 高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心,上海 200240;3. 中國(guó)船舶及海洋工程設(shè)計(jì)研究院,上海 200011; 4. 海軍裝備研究院,北京 100161)

    采用模型試驗(yàn)的方法研究了均勻流下柔性立管的渦激振動(dòng)(VIV)響應(yīng)特性及渦激力載荷特性。對(duì)均勻流場(chǎng)中柔性立管的VIV響應(yīng)特性進(jìn)行了分析,而后通過(guò)歐拉-伯努利梁動(dòng)態(tài)響應(yīng)控制方程和最小二乘法求取了柔性立管順流向(IL)和橫流向(CF)的渦激力系數(shù)。研究結(jié)果表明:均勻流下柔性立管的VIV為位移和主導(dǎo)頻率不隨時(shí)間變化的穩(wěn)態(tài)響應(yīng),順流向渦激振動(dòng)的主導(dǎo)頻率為橫流向的2倍;柔性立管的激勵(lì)系數(shù)與強(qiáng)迫振動(dòng)試驗(yàn)獲得的系數(shù)不一致,無(wú)因次頻率處于激勵(lì)區(qū)間的激勵(lì)系數(shù)存在負(fù)值,激勵(lì)系數(shù)不僅和無(wú)因次頻率及無(wú)因次振幅相關(guān),還與CFamp;IL方向位移相位角相關(guān);在無(wú)因次頻率0.13~ 0.22時(shí),橫流向的附加質(zhì)量系數(shù)在1.5~ 3.0振蕩變化;而順流向的附加質(zhì)量系數(shù)在無(wú)因次頻率在0.26~ 0.42內(nèi)從-1.0迅速增大到1.2后基本保持不變。

    柔性立管; 渦激振動(dòng);渦激力; 激勵(lì)系數(shù);附加質(zhì)量系數(shù)

    在洋流作用下,立管兩側(cè)會(huì)出現(xiàn)交替的泄渦,當(dāng)泄渦頻率接近于立管某一階固有頻率時(shí),會(huì)引起立管在順流向(In-Line, IL)及橫流向(Cross-Flow, CF)發(fā)生振動(dòng),即:渦激振動(dòng)(Vortex-Induced Vibration, VIV)。VIV的出現(xiàn)會(huì)使得立管產(chǎn)生嚴(yán)重的疲勞損傷。因此準(zhǔn)確預(yù)測(cè)立管的VIV在立管的設(shè)計(jì)中占有重要地位。在立管VIV的預(yù)報(bào)中,渦激力載荷的構(gòu)建直接決定著立管VIV預(yù)報(bào)結(jié)果的準(zhǔn)確性,而當(dāng)前VIV的預(yù)報(bào)中對(duì)渦激力載荷的構(gòu)建均是基于剛性圓柱體單自由度(純CF或純IL)強(qiáng)迫振蕩試驗(yàn)建立的渦激力系數(shù)數(shù)據(jù)庫(kù),包括激勵(lì)系數(shù)數(shù)據(jù)庫(kù)和附加質(zhì)量系數(shù)數(shù)據(jù)庫(kù)進(jìn)行的。在剛性圓柱體單自由度強(qiáng)迫振蕩試驗(yàn)中假設(shè)立管CF方向的渦激振動(dòng)與IL方向的渦激振動(dòng)互不影響,圓柱體只在純CF或純IL方向上運(yùn)動(dòng),忽略二者的耦合作用。然而,真實(shí)情況下,柔性立管的CF與IL方向均會(huì)發(fā)生VIV,且兩個(gè)方向上的VIV相互耦合,這種運(yùn)動(dòng)上的耦合使得兩個(gè)方向上的渦激力也相互影響,已有的研究表明,IL方向的VIV會(huì)使得CF方向渦激力的激勵(lì)系數(shù)數(shù)據(jù)庫(kù)產(chǎn)生新的峰值[1]。Sumer等[2]指出,當(dāng)CF方向VIV振幅超過(guò)0.2D~0.3D時(shí),脫落渦會(huì)變的整齊有序,瀉渦強(qiáng)度增大,此時(shí)IL方向的渦激力會(huì)被明顯放大。正是由于柔性立管雙向VIV的這種耦合作用,使得當(dāng)前對(duì)柔性立管VIV的預(yù)報(bào)結(jié)果不準(zhǔn)確[3],在海洋立管的設(shè)計(jì)中,不得不采用高達(dá)10倍以上的安全系數(shù)來(lái)確保其結(jié)構(gòu)設(shè)計(jì)的安全性[4]。然而,隨著油氣開(kāi)采進(jìn)一步向更深的水域發(fā)展,僅靠提高設(shè)計(jì)安全系數(shù)仍然無(wú)法確保管線結(jié)構(gòu)的安全,這極大地制約了油氣開(kāi)發(fā)向更深水域的發(fā)展。因此,研究柔性立管CF與IL方向均發(fā)生VIV時(shí)的渦激力載荷特性顯得越來(lái)越重要。

    計(jì)算流體力學(xué)方法(Computer Fluid Dynamics, CFD)是研究VIV這種復(fù)雜流固耦合問(wèn)題的理想方法[5-6]。為了研究渦激振動(dòng)發(fā)生時(shí)柔性立管的水動(dòng)力特性,Yamamoto等[7]采用直接渦流法估算柔性立管的水動(dòng)力,并利用準(zhǔn)3維數(shù)值模型研究柔性立管的振動(dòng)與其水動(dòng)力之間的關(guān)系。Evangelinos等[8]利用直接數(shù)值模擬方法模擬剛性圓柱體以及柔性圓柱體的3 維流場(chǎng),并得到了流場(chǎng)作用在剛性圓柱體和柔性圓柱體上的阻力。然而上述數(shù)值模擬結(jié)果均未得到試驗(yàn)的驗(yàn)證。此外,CFD對(duì)流場(chǎng)和結(jié)構(gòu)的網(wǎng)格質(zhì)量要求十分精細(xì),其計(jì)算量特別巨大,即便是一個(gè)6 m長(zhǎng)的試驗(yàn)立管,每個(gè)流速工況計(jì)算都需耗時(shí)數(shù)月[9],因而目前采用CFD方法系統(tǒng)的研究細(xì)長(zhǎng)柔性立管的渦激力是不太現(xiàn)實(shí)的。

    鑒于CFD的局限性,Mukundan等[10-11]強(qiáng)迫振蕩試驗(yàn)獲得的激勵(lì)系數(shù)數(shù)據(jù)庫(kù)參數(shù)化,利用將經(jīng)驗(yàn)?zāi)P蚔IVA預(yù)測(cè)結(jié)果與試驗(yàn)結(jié)果之間的誤差最小化得到最優(yōu)參數(shù)的方法,獲得了CF方向上新的激勵(lì)系數(shù)數(shù)據(jù)庫(kù),新的數(shù)據(jù)庫(kù)具有更大的激勵(lì)區(qū)間,且主激勵(lì)區(qū)間與第二激勵(lì)區(qū)間混合在一起。然而,Mukundan的方法不能用于研究立管CF方向VIV高階響應(yīng)下的渦激力以及IL方向的渦激力。Huarte等[12]利用試驗(yàn)中測(cè)得的剪切流下豎直頂張力立管模型CF與IL方向上的應(yīng)變信息結(jié)合梁有限元方程得到了作用在立管CF與IL方向上的水動(dòng)力載荷,然而其沒(méi)有進(jìn)一步獲得立管CF與IL方向的渦激力系數(shù),即:激勵(lì)系數(shù)和附加質(zhì)量系數(shù)。Wu等[13-14]分別利用梁有限元狀態(tài)矢量空間方程和基于最優(yōu)控制理論的逆分析法獲得了剪切流下立管CF方向上的渦激力及其載荷系數(shù),然而二者均沒(méi)有對(duì)立管IL方向上的渦激力進(jìn)行研究。

    本文通過(guò)利用拖車(chē)拖動(dòng)兩端承受恒定預(yù)張力的柔性立管模型在拖曳水池中做勻速運(yùn)動(dòng)的方法,研究均勻流下CF與IL方向均發(fā)生VIV時(shí)的柔性立管VIV響應(yīng)特性及渦激力特性。文中利用試驗(yàn)中測(cè)得的立管模型的VIV應(yīng)變信息,采用FFT(Fast Fourier Transform)、小波分析法和模態(tài)分析法研究均勻流場(chǎng)中柔性立管CF與IL方向上的VIV響應(yīng)特性。同時(shí),利用歐拉-伯努利梁動(dòng)態(tài)響應(yīng)控制方程和最小二乘法求取VIV發(fā)生時(shí)柔性立管CF與IL方向的渦激力系數(shù),包括激勵(lì)系數(shù)和附加質(zhì)量系數(shù)。

    1 試驗(yàn)描述

    本試驗(yàn)在拖曳水池中進(jìn)行,水池尺寸為192 m×10 m×4.2 m(長(zhǎng)×寬×深),立管模型橫置于水池中,其中心距離水面1.5 m,試驗(yàn)裝置示意圖如圖1所示。立管模型利用萬(wàn)向節(jié)與端部裝置進(jìn)行連接,可在IL與CF方向上彎曲并可以沿軸向方向運(yùn)動(dòng)。端部裝置可以為立管模型提供恒定預(yù)張力。試驗(yàn)中均勻來(lái)流的模擬方法是將立管橫置于拖曳水池中,由拖車(chē)帶動(dòng)立管在拖曳水池中勻速前進(jìn),從而形成相對(duì)均勻來(lái)流。

    圖1 試驗(yàn)?zāi)P秃?jiǎn)圖Fig. 1 Schematic of the test apparatus

    試驗(yàn)中所使用的立管模型為縮尺模型,模型的外徑為30 mm,有效長(zhǎng)度為7.9 m,模型具體參數(shù)如表1所示。表1中的結(jié)構(gòu)阻尼比為立管模型在空氣中的結(jié)構(gòu)阻尼比。

    表1 立管模型基本參數(shù)

    立管模型表面布置有88個(gè)光纖光柵應(yīng)變傳感器,分別布置于模型的CF1、CF2、IL1以及IL2四個(gè)方向,用于測(cè)量模型在此四個(gè)方向上的應(yīng)變,即εCF1、εCF2、εIL1和εIL2,如圖2所示。其中CF1和CF2方向上的傳感器分別為19個(gè),IL1和IL2方向上分別為25個(gè)。傳感器在CF與IL方向上均勻分布,CF與IL方向上相鄰測(cè)點(diǎn)的間距分別為0.42 m和0.315 m。光纖光柵傳感器的采樣頻率為250 Hz。

    圖2 立管模型表面光纖光柵傳感器布置方式Fig. 2 Arrangement of strain gauges on the surface of the riser model

    2 數(shù)據(jù)處理

    2.1 渦激振動(dòng)產(chǎn)生的彎曲應(yīng)變

    試驗(yàn)中受預(yù)張力的立管模型,由光纖光柵應(yīng)變傳感器測(cè)量的立管模型CF方向各測(cè)點(diǎn)處的應(yīng)變?chǔ)臗F1和εCF2均包含兩部分:軸向張力產(chǎn)生的拉伸應(yīng)變?chǔ)臗F-T以及VIV產(chǎn)生的彎曲應(yīng)變?chǔ)臗F

    εCF1=εCF-T+εCF
    εCF2=εCF-T-εCF

    (1)

    對(duì)式(1)進(jìn)行簡(jiǎn)單的變換, 便可得到CF方向由VIV引起的彎曲應(yīng)變?chǔ)臗F

    εCF=[εCF1-εCF2]/2

    (2)

    立管模型IL方向各測(cè)點(diǎn)處由光纖光柵應(yīng)變傳感器測(cè)量的應(yīng)變?chǔ)臝L1和εIL2包含三部分:張力產(chǎn)生的拉伸應(yīng)變?chǔ)臝L-T、VIV產(chǎn)生的彎曲應(yīng)變?chǔ)臝L以及平均拖曳力產(chǎn)生的彎曲應(yīng)變?chǔ)舖b

    εIL1=εIL-T+εIL+εmb
    εIL2=εIL-T-εIL-εmb

    (3)

    由于立管模型IL方向的VIV為均值為零的周期性振動(dòng),故由VIV引起的彎曲應(yīng)變?chǔ)臝L的時(shí)間平均值為零,即

    (4)

    此外,由平均拖曳力產(chǎn)生的彎曲應(yīng)變?chǔ)舖b不隨時(shí)間變化,因此在試驗(yàn)穩(wěn)定段內(nèi)平均彎曲應(yīng)變?chǔ)舖b滿足

    (5)

    由式(3) 可以得到

    εmb+εIL=[εIL1-εIL2]/2

    (6)

    對(duì)式(6)的兩端求時(shí)間的平均值,結(jié)合式(5)和式(4)可得到立管模型IL方向各測(cè)點(diǎn)處由VIV引起的彎曲應(yīng)變?chǔ)臝L

    (7)

    2.2 柔性立管渦激力載荷

    將試驗(yàn)中所用的立管模型理想化為受張力的歐拉-伯努利梁[15],在忽略立管模型軸向扭轉(zhuǎn)變形的基礎(chǔ)上,立管CF與IL方向上的VIV可以通過(guò)歐拉-伯努利梁的結(jié)構(gòu)動(dòng)力學(xué)方程進(jìn)行描述。定義坐標(biāo)系O-XYZ:坐標(biāo)系的原點(diǎn)位于立管模型的某一端點(diǎn);Z軸沿立管的軸向方向;X軸沿來(lái)流方向,即IL方向;Y軸與X軸和Z軸相互垂直,即CF方向,如圖1所示。在此坐標(biāo)系下,立管模型CF和IL方向VIV的控制方程可表示為

    (8)

    (9)

    式中:EI為立管模型的彎曲剛度;m為立管模型單位長(zhǎng)度質(zhì)量;C為立管模型在空氣中的結(jié)構(gòu)阻尼;T為VIV發(fā)生時(shí)立管兩端軸向張力的時(shí)間平均值;fx(z,t)和fy(z,t)分別為使得立管在IL與CF 方向上發(fā)生VIV的渦激力;x(z,t)和y(z,t)分別為立管模型IL與CF方向VIV的位移;?y(z,t)/?t與?2y(z,t)/?t2分別為立管模型CF方向VIV的速度和加速度;?x(z,t)/?t與?2x(z,t)/?t2分別為立管模型IL方向VIV的速度和加速度。

    立管IL與CF方向的VIV位移x(z,t)和y(z,t)可由式(2)和式(7)求得的立管CF與IL方向的VIV彎曲應(yīng)變結(jié)合模態(tài)分析法[16]求得。立管CF與IL方向的VIV速度和加速度可由對(duì)VIV位移求時(shí)間的一階偏導(dǎo)和二階偏導(dǎo)獲得。

    2.3 柔性立管渦激力載荷系數(shù)

    本節(jié)中將給出柔性立管CF與IL方向渦激力載荷系數(shù),即:激勵(lì)系數(shù)與附加質(zhì)量系數(shù)的求解方法。以立管CF方向?yàn)槔?/p>

    (10)

    式中:y0(z)為橫截面z處的VIV位移響應(yīng)幅值;ω為立管模型CF方向VIV響應(yīng)圓頻率。

    在此情況下,立管橫截面z處的渦激力載荷fy(z,t)可表示為

    fy(z,t)=f0(z)sin(ωt+φ(z))

    (11)

    式中:f0(z)為渦激力載荷的幅值;φ(z)為橫截面z處的渦激力與此橫截面處VIV位移之間的相位角。將式(11)展開(kāi),可得

    fy(z,t)=f0(z)sin(φ(z))cos(ωt)-
    [-f0(z)cos(φ(z))sin(ωt)]

    (12)

    式中,f0(z)sin(φ(z))cos(ωt)為與速度同相位的渦激力,可將其表示為

    (13)

    式中:CLe(z)為橫截面z處的激勵(lì)力系數(shù);ρ為流體密度;U為流速;D為立管模型的水動(dòng)力直徑。

    式(12)中的-f0sin(φ)cos(ωt)項(xiàng)為與加速度同相位的渦激力,即附加質(zhì)量力,可以表示為

    (14)

    式中,CLa(z)為橫截面z處的附加質(zhì)量系數(shù)。

    由式(13)和式(14)可得橫截面z處的渦激力載荷fy(z,t)表示為

    (15)

    若已知立管橫截面z處的VIV渦激力載荷以及VIV

    響應(yīng)速度和加速度后,可采用最小二乘法求解立管橫截面z處的激勵(lì)系數(shù)CLe(z)和附加質(zhì)量系數(shù)CLa(z),具體過(guò)程為:

    將根據(jù)式(8)計(jì)算得到的VIV渦激力載荷fy(z,t)看作為t時(shí)刻立管橫截面z處渦激力載荷的測(cè)量值,將fp(z,t)看作同一時(shí)刻立管橫截面z處渦激力載荷的預(yù)測(cè)值,如式(16)表示

    (16)

    根據(jù)最小二乘法基本原理,激勵(lì)系數(shù)CLe(z)和附加質(zhì)量系數(shù)CLa(z)的選取要使得整個(gè)采樣時(shí)間內(nèi)渦激力載荷的預(yù)測(cè)值與測(cè)量值誤差的平方和為最小

    (17)

    式中:n為整個(gè)采樣時(shí)間內(nèi)的采樣點(diǎn)個(gè)數(shù);e2(z)為整個(gè)采樣時(shí)間內(nèi)立管橫截面z處渦激力載荷預(yù)測(cè)值與測(cè)量值誤差的平方和。

    將式(17)的右邊項(xiàng)進(jìn)行展開(kāi),可得

    (18)

    則式(18)可簡(jiǎn)化為

    (19)

    若要使得激勵(lì)系數(shù)CLe(z)和附加質(zhì)量系數(shù)CLa(z)的選取使得e2(z)為最小,則

    (20)

    結(jié)合式(19)和式(20),可得

    (21)

    由式(21)可求得立管橫截面z處CF方向上的VIV激勵(lì)系數(shù)CLe(z)和附加質(zhì)量系數(shù)CLa(z),如式(22)所示

    (22)

    柔性立管IL方向上各橫截面處的VIV激勵(lì)系數(shù)Cde(z)和附加質(zhì)量系數(shù)Cda(z)可采用與CF方向同樣的方法求得。

    3 結(jié)果分析與討論

    3.1 立管模型VIV響應(yīng)特性

    流速2.6 m/s下立管模型在測(cè)點(diǎn)Z=3.95 m處CF與IL方向VIV彎曲應(yīng)變的時(shí)歷曲線、頻率譜及小波分析結(jié)果,如圖3所示。從圖3可知:①立管模型CF和IL方向上的VIV均為幅值穩(wěn)定的穩(wěn)態(tài)振動(dòng);②立管模型CF與IL方向VIV的主導(dǎo)頻率分別為16.79 Hz和33.62 Hz,立管模型CF與IL方向上VIV的主導(dǎo)頻率隨時(shí)間的發(fā)展均保持不變,且IL方向上的主導(dǎo)頻率約為CF方向上主導(dǎo)頻率的2倍,F(xiàn)u等[18-19]都發(fā)現(xiàn)了VIV的這種雙倍頻率現(xiàn)象;③立管CF方向的VIV出現(xiàn)了明顯的3倍主導(dǎo)頻率的高階響應(yīng),Vandiver等[20]也發(fā)現(xiàn)了VIV的這種高階響應(yīng)。

    圖 3 流速2.6 m/s下立管模型CF與IL方向上測(cè)點(diǎn)Z=3.95 m處VIV應(yīng)變時(shí)歷、頻率譜及小波分析結(jié)果Fig. 3 Time histories, frequency spectra and result of the wavelet analysis of strain because of VIV in the CF and IL directions at Z=3.95 m for a 2.6 m/s current

    本文只關(guān)心立管模型CF與IL方向主導(dǎo)頻率下的VIV響應(yīng)及其渦激力,故采用帶通濾波對(duì)采集到的立管模型的VIV應(yīng)變信號(hào)進(jìn)行濾波處理,從而獲得主導(dǎo)頻率下的VIV響應(yīng)。流速2.6 m/s 下立管模型CF與IL方向的濾波寬度為15.0~18.0 Hz 和32.0~35.0 Hz。

    在通過(guò)濾波得到流速2.6 m/s立管模型CF與IL方向主導(dǎo)頻率下的VIV應(yīng)變后,利用模態(tài)分析法求解立管模型CF與IL方向的VIV位移響應(yīng)。圖4為立管模型CF與IL方向VIV無(wú)因次位移RMS值([y/D]RMS和[x/D]RMS)沿立管軸向的分布。從圖4可知,立管VIV響應(yīng)的主導(dǎo)模態(tài)分別為5階和7階。如圖4所示,立管模型CF方向的VIV位移響應(yīng)遠(yuǎn)大于IL方向的VIV位移響應(yīng),而圖3顯示立管CF與IL方向的VIV應(yīng)變?cè)谕涣考?jí)。這表明同一流速下,雖然立管CF方向的VIV位移大于IL方向的VIV位移,然而兩個(gè)方向的VIV應(yīng)變?cè)谕涣考?jí),考慮到IL方向的VIV響應(yīng)頻率為CF的兩倍,故而相比于CF方向的VIV,IL方向的VIV會(huì)使得立管產(chǎn)生同量級(jí)甚至更為嚴(yán)重的疲勞損傷。

    3.2 柔性立管激勵(lì)系數(shù)

    (a)

    (b)圖4 流速2.6 m/s下立管模型CF與IL方向VIV的無(wú)因次位移RMS值沿軸向分布Fig. 4 Axial distributions of the RMS values of the non-dimensionalized VIV displacement in the CF and ILdirectionsfor a 2.6 m/s current

    圖5給出了流速2.6 m/s下立管模型CF與IL方向的VIV激勵(lì)系數(shù)CLe和Cde以及無(wú)因次幅值(A/D)CF和(A/D)IL沿立管模型軸向的分布。圖5顯示在振動(dòng)節(jié)點(diǎn)處,激勵(lì)系數(shù)會(huì)發(fā)生突變,這是由于振動(dòng)節(jié)點(diǎn)處的振動(dòng)幅值太小導(dǎo)致的?;趧傂詧A柱體強(qiáng)迫振蕩試驗(yàn)的半經(jīng)驗(yàn)頻域預(yù)報(bào)理論指出無(wú)因次頻率在0.125~0.2為VIV的激勵(lì)區(qū)[21],本試驗(yàn)中流速2.6 m/s下立管CF方向VIV的無(wú)因次頻率為0.192,因此整個(gè)立管在CF方向應(yīng)該均處于激勵(lì)區(qū);但從圖5可知,立管模型某些區(qū)域內(nèi)的激勵(lì)系數(shù)CLe為負(fù)值,這表明均勻流速下,無(wú)因次頻率處于激勵(lì)區(qū)間的立管仍存在著激勵(lì)區(qū)和阻尼區(qū),在激勵(lì)區(qū)內(nèi),能量從流場(chǎng)傳到立管,而在阻尼區(qū)內(nèi),能量從立管傳到流場(chǎng)。從圖5還可知,在激勵(lì)系數(shù)CLe為負(fù)值的區(qū)域內(nèi)立管的無(wú)因次幅值并不總是最大的,這與Gopalkrishnan所發(fā)的規(guī)律是不一致的,Gopalkrishnan指出在同一無(wú)因次頻率下,只有當(dāng)幅值增大到某一值后,其激勵(lì)系數(shù)才變?yōu)樨?fù)值。此外,還可發(fā)現(xiàn),同一流速,即同一無(wú)因次頻率下,無(wú)因次幅值相同的立管各橫截面處的激勵(lì)系數(shù)并不相同,甚至差別很大,這是由于立管各橫截面的運(yùn)動(dòng)軌跡以及CF與IL方向VIV位移相位角的不同,導(dǎo)致立管尾流區(qū)具有不同的泄渦模式,繼而引起VIV水動(dòng)力的變化[22]。這表明:柔性立管的激勵(lì)系數(shù)不僅與無(wú)因次頻率和無(wú)因次振幅相關(guān),還與CFamp;IL方向的VIV位移相位角相關(guān)。立管IL方向上VIV的激勵(lì)系數(shù)Cde沿立管軸向的分布規(guī)律與CF方向的激勵(lì)系數(shù)CLe有類(lèi)似的規(guī)律。

    (a)

    (b)圖5 流速2.6 m/s下立管模型CF與IL方向激勵(lì)系數(shù)及無(wú)因次位移幅值沿模型軸向分布Fig.5 Axial distributions of the CF and IL excitation coefficients and non-dimensionalized response amplitudes for a 2.6 m/s current

    3.3 柔性立管附加質(zhì)量系數(shù)

    圖6給出了流速2.6 m/s下立管模型CF與IL方向的附加質(zhì)量系數(shù)CLa和Cda以及無(wú)因次幅值(A/D)CF和(A/D)IL沿立管模型軸向的分布。從圖6可知,立管CF與IL方向的振動(dòng)節(jié)點(diǎn)處,附加質(zhì)量會(huì)發(fā)生突變,其原因與激勵(lì)系數(shù)的一致。圖6顯示,均勻流場(chǎng)中當(dāng)立管的CF與IL方向均發(fā)生VIV時(shí),在同一流速,即同一無(wú)因次頻率下,立管CF方向上不同無(wú)因次幅值下的附加質(zhì)量系數(shù)并不相等,IL方向亦是如此。這與半經(jīng)驗(yàn)頻域預(yù)報(bào)模型VIVANA以及Aronsen[23]認(rèn)為在一定的振動(dòng)幅值范圍內(nèi),附加質(zhì)量系數(shù)與VIV振動(dòng)幅值關(guān)系不大,只與無(wú)因次頻率的結(jié)論并不一致。

    (a)

    (b)圖6 流速2.6 m/s下立管模型CF與IL方向附加質(zhì)量系數(shù)及無(wú)因次位移幅值沿模型軸向分布Fig. 6 Axial distributions of CF and IL added-mass coefficients and non-dimensionalized vibration amplitudes for a 2.6 m/s current

    為了檢驗(yàn)本文中用于計(jì)算立管模型渦激力載荷方法的正確性,利用圖5和圖6中的激勵(lì)系數(shù)和附加質(zhì)量系數(shù),根據(jù)式(15)計(jì)算立管模型在流速2.6 m/s下CF與IL方向的渦激力。然后,根據(jù)表1中的立管模型參數(shù)在有限元軟件ABAQUS中建立有限元模型,而后利用獲得的渦激力計(jì)算立管模型的應(yīng)變響應(yīng)。圖7為流速2.6 m/s下立管模型CF與IL方向測(cè)點(diǎn)Z=3.95 m處應(yīng)變的測(cè)量值和計(jì)算值時(shí)歷以及應(yīng)變測(cè)量值和計(jì)算值的RMS沿立管軸向的分布。從圖7可知,立管CF與IL方向應(yīng)變響應(yīng)的測(cè)量值和計(jì)算值非常吻合,這表明了本文中用于計(jì)算VIV發(fā)生時(shí)立管模型渦激力載荷的方法是正確的。

    圖7 流速2.6 m/s下立管模型CF與IL方向VIV應(yīng)變的測(cè)量值及計(jì)算值在測(cè)點(diǎn)Z=3.95 m處的時(shí)歷曲線以及應(yīng)變RMS值沿立管軸向的分布Fig.7 Time histories of the measured and calculated values of CF and IL VIV strains at Z=3.95 m and axial distributions of the RMS values of the measured and calculated strains for a 2.6 m/s current

    3.4 柔性立管平均附加質(zhì)量系數(shù)

    為了研究附加質(zhì)量系數(shù)與無(wú)因次頻率之間的關(guān)系,忽略VIV振幅對(duì)附加質(zhì)量系數(shù)的影響,將附加質(zhì)量系數(shù)在整個(gè)立管模型上進(jìn)行平均,得到試驗(yàn)中各流速下CF與IL方向上的平均附加質(zhì)量系數(shù)。根據(jù)立管模型的參數(shù),建立立管模型的有限元模型,并將得到的平均附加質(zhì)量系數(shù)帶入到有限元模型中,計(jì)算立管模型CF與IL方向VIV主導(dǎo)模態(tài)對(duì)應(yīng)的固有頻率,將其與試驗(yàn)測(cè)得的主導(dǎo)頻率相比較,如圖8所示。從圖8可以看出,利用平均附加質(zhì)量系數(shù)計(jì)算得到的主導(dǎo)頻率與試驗(yàn)測(cè)量值幾乎相同。這表明雖然當(dāng)CF與IL方向同時(shí)發(fā)生VIV時(shí),立管CF與IL方向上各橫截面處的附加質(zhì)量系數(shù)不相等,但是在進(jìn)行立管的VIV預(yù)報(bào)時(shí),仍可以忽略VIV振動(dòng)幅值對(duì)附加質(zhì)量系數(shù)的影響,采用平均附加質(zhì)量系數(shù)來(lái)計(jì)算立管的固有頻率。

    (a) CF

    (b) IL圖8 不同流速下立管模型CF與IL方向VIV主導(dǎo)頻率的計(jì)算值和測(cè)量值Fig. 8 Comparison of the calculated and measured values of the dominant frequencies in the CF and IL direction fordifferent currents

    將不同流速下立管模型CF的VIV主導(dǎo)頻率進(jìn)行無(wú)因次化,給出CF方向的平均附加質(zhì)量系數(shù)與無(wú)因次頻率的函數(shù)曲線,如圖9所示。在圖9中同時(shí)給出了半經(jīng)驗(yàn)頻域預(yù)報(bào)模型VIVANA 中的附加質(zhì)量系數(shù)與無(wú)因次頻率關(guān)系曲線。從圖9可知,本文中立管CF方向的附加質(zhì)量系數(shù)與無(wú)因次頻率的關(guān)系與VIVANA給出的曲線并不相符,當(dāng)無(wú)因次頻率在0.13~0.16時(shí),VIVANA中的附加質(zhì)量系數(shù)為負(fù)值,大致保持為-0.6,而本文的附加質(zhì)量系數(shù)大于2.5;在無(wú)因次頻率0.16~0.22,VIVANA給出的附加質(zhì)量系數(shù)隨著無(wú)因次頻率的增加迅速增大到約2.0后基本保持不變,而本文的附加質(zhì)量系數(shù)在此無(wú)因次頻率區(qū)間內(nèi)在1.5~3.0振蕩變化。

    圖9 立管CF方向上平均附加質(zhì)量系數(shù)隨無(wú)因次頻率的分布Fig. 9 Distribution of the added-mass coefficient with the non-dimensional frequency in the CF direction

    圖10給出了立管模型IL方向的平均附加質(zhì)量系數(shù)與無(wú)因次頻率的函數(shù)曲線?;趧傂詧A柱體純IL方向的強(qiáng)迫振蕩試驗(yàn)結(jié)果,Aronsen給出了IL方向附加質(zhì)量系數(shù)與其無(wú)因次頻率之間的關(guān)系曲線,如圖10所示。圖10顯示:同一無(wú)因次頻率下,本文中得到的IL方向的附加質(zhì)量系數(shù)基本上大于Aronsen的值,且附加質(zhì)量系數(shù)隨無(wú)因次頻率的變化關(guān)系也不太一致。在Aronsen給出的關(guān)系曲線中,當(dāng)無(wú)因次頻率從0.26增大到0.42時(shí),附加質(zhì)量系數(shù)從-0.3緩慢的增大到0.75;而在本文中,當(dāng)無(wú)因次頻率從0.26增大到0.32時(shí),附加質(zhì)量系數(shù)從-1.0迅速增大到1.2,而后,隨著無(wú)因次頻率的增大基本保持不變。

    圖10 立管IL方向上平均附加質(zhì)量系數(shù)隨無(wú)因次頻率的分布Fig. 10 The distribution of the added-mass coefficient with the non-dimensional frequency in the IL direction

    4 結(jié) 論

    本文采用水池模型試驗(yàn)的方法研究了均勻流下柔性立管的渦激振動(dòng)響應(yīng)特性及渦激力載荷特性。文中采用FFT、小波分析法和模態(tài)分析法對(duì)均勻流場(chǎng)中柔性立管的VIV響應(yīng)特性進(jìn)行了分析,并通過(guò)歐拉-伯努利梁動(dòng)態(tài)響應(yīng)控制方程和最小二乘法求取了柔性立管IL和CF方向的渦激力系數(shù)。本文的主要結(jié)論如下:

    (1)均勻流下柔性立管的VIV為響應(yīng)幅值和主導(dǎo)頻率不隨時(shí)間變化的穩(wěn)態(tài)響應(yīng),IL方向的主導(dǎo)頻率為CF的兩倍;CF方向的VIV位移大于IL方向,但兩個(gè)方向的應(yīng)變?cè)谕涣考?jí)。

    (2)柔性立管的激勵(lì)系數(shù)與圓柱體強(qiáng)迫振動(dòng)試驗(yàn)獲得的系數(shù)不一致:無(wú)因次頻率處于激勵(lì)區(qū)間的激勵(lì)系數(shù)存在負(fù)值,激勵(lì)系數(shù)不僅和無(wú)因次頻率及振幅相關(guān),還與CFamp;IL方向位移相位角相關(guān)。

    (3)柔性立管的附加質(zhì)量系數(shù)為無(wú)因次頻率和VIV響應(yīng)幅值的函數(shù);但可以忽略VIV響應(yīng)幅值對(duì)附加質(zhì)量系數(shù)的影響,采用平均附加質(zhì)量系數(shù)計(jì)算立管的固有頻率。

    (4)當(dāng)無(wú)因次頻率為0.13~0.22時(shí),CF方向的平均附加質(zhì)量系數(shù)為1.5~3.0;當(dāng)無(wú)因次頻率為0.26~0.42時(shí),IL方向的平均附加質(zhì)量系數(shù)從-1.0迅速增大到1.2,而后基本保持不變。

    [ 1 ] MARCOLLOA H, HINWOODB J B. On shear flow single mode lock-in with both cross-flow and in-line lock-in mechanisms[J]. Journal of Fluids and Structures, 2006, 22(2): 197-211.

    [ 2 ] SUMER B M, JERGEN F. Hydrodynamics around cylindrical structures[M]. Singapore: World Scientific, 2006: 334-

    413.

    [ 3 ] CHAPLIN J R, BEARMAN P W, FONTAINE E, et al. Blind predictions of laboratory measurements of vortex-induced vibrations of a tension riser[C]∥Fluid Induced Vibration Conference. [S. l.]: Journal of Fluids and Structures, 2005: 25-40.

    [ 4 ] Design of risers for floating production systems (FPSs) and tension leg platforms (TLPs): API RP 2RD[R]. Washington: American Petroleum Institute, 1998.

    [ 5 ] KAIKTSIS L, TRIANTAFYLLOU G S, ?ZBAS M. Excitation, inertia, and drag forces on a cylinder vibrating transversely to a steady flow[J]. Journal of Fluids and Structures,2007, 23(1): 1-21.

    [ 6 ] SARPKAYA T. A critical review of the intrinsic nature of vortex-induced vibrations[J]. Journal of Fluids and Structures, 2004, 19(4): 389-447.

    [ 7 ] YAMAMOTO C T, MENEGHINI J R, SALTARA F, et al. Numerical simulations of vortex-induced vibration on flexible cylinders[J]. Journal of Fluids and Structures, 2004, 19(4): 467-489.

    [ 8 ] EVANGELINOS C, LUCOR D, KARNIADAKAS G E. DNS-derived force distribution on flexible cylinders subject to vortex-induced vibration[J]. Journal of Fluids and Structures, 2004, 14(3): 429-440.

    [ 9 ] HUANG Z. CFD simulation of riser VIV[D]. College Station: Texas Aamp;M University, 2012.

    [10] MUKUNDAN H. Vortex induced vibration and force recons-truction from field and experimental data[D]. Cambridge: Massachusetts Institute of Technology, 2008.

    [11] GOPALKRISHNAN R. Vortex-induced forces on oscillating bluff cylinders[D]. Cambridge: Massachusetts Institute of Technology, 1993.

    [12] HUARTE F J H, BEARMAN P W, CHAPLIN J R. On the force distribution along the axis of a flexible circular cylinder undergoing multi-mode vortex-induced vibrations[J]. Journal of Fluids and Structures,2006, 22(6/7): 897-903.

    [13] WU J, LARSEN C M, KAASEN K E. A new approach for identification of forces on slender beams subjected to vortex induced vibrations[C]∥ Proceedings of Offshore Mechanics and Arctic Engineering Conference. Estoril: OMAE, 2008.

    [14] WU J, MAINON P, LARSEN C M, et al. VIV force identification using classical optimal control algorithm[C]∥ Proceedings of Offshore Mechanics and Arctic Engineering Conference. Hawaii: OMAE, 2009.

    [15] 陳鐵云, 陳伯真. 船舶結(jié)構(gòu)力學(xué)[M]. 北京:國(guó)防工業(yè)出版社, 1984: 29-30.

    [16] LI L, FU S, YANG J, et al. Experimental investigation on vortex-induced vibration of risers with staggered buoyancy[C]∥ ASME 2011 30th International Conference on Ocean, Offshore and Arctic Engineering. Rotterdam: OMAE, 2011.

    [17] VIKESTAD K, VANDIVER J K, LARSEN C M. Added mass and oscillation frequency for circular cylinder subjected to VIV and external disturbance[J]. Journal of Fluids and Structures, 2000, 14(7): 1071-1088.

    [18] FU S, REN T, LI R, et al. Experimental investigation on VIV of the flexible model within full scale Re number regime[C]∥ ASME 2011 30th International Conference on Ocean, Offshore and Arctic Engineering. Rotterdam: OMAE, 2011.

    [19] FANG S M, NIEDZWECKI J M, FU S, et al. VIV response of a flexible cylinder with varied coverage by buoyancy elements and helical strakes[J]. Marine Structures, 2014, 39: 70-89.

    [20] VANDIVER J K, JAISWAL V, JHINGRAN V. Insights on vortex-induced, traveling waves on long risers[J]. Journal of Fluids and Structures,2009, 25(4): 641-653.

    [21] LARSEN C M, YTTERVIK R, PASSANO E, et al. VIVANA theory manual[M].Trondheim: Norway, 2001.

    [22] YIN D, LARSEN C M. Experimental and numerical analysis of forced motion of a circular cylinder[C]∥ Proceedings of Offshore Mechanics and Arctic Engineering Conference. Rotterdam: OM, 2011.

    [23] AROSEN K H. An experimental investigation of in-line and combined in-line and cross-flow vortex induced vibrations[D]. Trondheim: NTNU, 2007.

    Structuralresponsesandvortex-inducedforceofflexiblerisersundergoingvortex-inducedvibrationinuniformflow

    SONG Leijian1, 2, FU Shixiao1,2, REN Tie1,2,3, YU Dapeng4, ZHANG Mengmeng1,2

    (1. State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China;2. Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, Shanghai Jiao Tong University, Shanghai 200240, China;3. Marine Design amp; Research Institute of China, Shanghai 200011, China;4. Naval Academy of Armament, Beijing 100161, China)

    In this study, the structural responses and vortex-induced forces of flexible risers undergoing vortex-induced vibration (VIV) in uniform flow were investigated. The VIV response characteristics were analyzed using the measured strains. Then, using the Euler Bernoulli beam structure dynamic equation and the least square method, the vortex-induced force coefficients in cross-flow (CF) and in-line (IL) direction were computed. The results indicate that the VIV of the flexible riser under uniform current is steady vibration with the amplitude and dominant frequency independent of time and the dominant frequency in the IL direction is twice as much as that in the CF direction. The excitation coefficients of the flexible riser do not always agree with those obtained by the forced oscillation tests. In the CF direction, when the non-dimensional frequency varies from 0.13 to 0.22, the added-mass coefficient oscillates between 1.5 and 3.0. In the IL direction, the added mass coefficient increases quickly from -1.0 to 1.2 and then keeps constantly when the non-dimensional frequency varies from 0.26 to 0.42.

    flexible riser; vortex-induced vibration; hydrodynamic force; excitation coefficient; added-mass coefficient

    國(guó)家自然科學(xué)基金資助項(xiàng)目(51239007; 51279101;51490674;51490675)

    2016-03-29 修改稿收到日期: 2016-07-26

    宋磊建 男,博士生,1987年生

    付世曉 男,博士,研究員,博士生導(dǎo)師,1976年生

    TH212;TH213.3

    A

    10.13465/j.cnki.jvs.2017.22.003

    猜你喜歡
    立管幅值柔性
    一種柔性拋光打磨頭設(shè)計(jì)
    灌注式半柔性路面研究進(jìn)展(1)——半柔性混合料組成設(shè)計(jì)
    石油瀝青(2021年5期)2021-12-02 03:21:18
    高校學(xué)生管理工作中柔性管理模式應(yīng)用探索
    常見(jiàn)高層建筑物室內(nèi)給水立管材質(zhì)解析
    基于S變換的交流電網(wǎng)幅值檢測(cè)系統(tǒng)計(jì)算機(jī)仿真研究
    電子制作(2017年7期)2017-06-05 09:36:13
    正序電壓幅值檢測(cè)及諧波抑制的改進(jìn)
    深水鋼懸鏈立管J型鋪設(shè)研究
    The Power of Integration
    Beijing Review(2015年43期)2015-11-25 03:12:04
    低壓電力線信道脈沖噪聲的幅值與寬度特征
    海洋立管濕模態(tài)振動(dòng)分析
    404 Not Found

    404 Not Found


    nginx
    国产1区2区3区精品| 99精品欧美一区二区三区四区| 国产三级黄色录像| 香蕉国产在线看| 可以免费在线观看a视频的电影网站| 黄色怎么调成土黄色| 亚洲欧美日韩另类电影网站| 国产不卡av网站在线观看| 亚洲av第一区精品v没综合| 久久午夜亚洲精品久久| 高潮久久久久久久久久久不卡| 九色亚洲精品在线播放| 久久精品国产亚洲av香蕉五月 | 欧美在线一区亚洲| 日韩大片免费观看网站| 一级毛片精品| 黄色毛片三级朝国网站| 老司机亚洲免费影院| 欧美成人午夜精品| 人人妻,人人澡人人爽秒播| 欧美变态另类bdsm刘玥| 五月开心婷婷网| 国产成人精品久久二区二区91| 久久香蕉激情| 亚洲精品国产一区二区精华液| 亚洲精品久久午夜乱码| 99九九在线精品视频| 亚洲精品一卡2卡三卡4卡5卡| 每晚都被弄得嗷嗷叫到高潮| 12—13女人毛片做爰片一| 亚洲全国av大片| 乱人伦中国视频| 宅男免费午夜| 久久免费观看电影| 国产精品欧美亚洲77777| 一本久久精品| 精品人妻1区二区| 国产xxxxx性猛交| 国产精品久久久人人做人人爽| 天堂中文最新版在线下载| 一进一出抽搐动态| 真人做人爱边吃奶动态| 亚洲伊人色综图| 亚洲欧美色中文字幕在线| 一区二区三区精品91| 国产免费福利视频在线观看| 人成视频在线观看免费观看| 超色免费av| 亚洲色图 男人天堂 中文字幕| 黄片大片在线免费观看| 国产麻豆69| 首页视频小说图片口味搜索| 高清在线国产一区| 2018国产大陆天天弄谢| 91麻豆av在线| 人人澡人人妻人| 五月天丁香电影| 黄色 视频免费看| 伊人久久大香线蕉亚洲五| 国产精品久久久av美女十八| 午夜福利欧美成人| 99国产精品免费福利视频| 色94色欧美一区二区| 国产精品久久久久久人妻精品电影 | 黄网站色视频无遮挡免费观看| 每晚都被弄得嗷嗷叫到高潮| 天堂8中文在线网| 欧美一级毛片孕妇| 9色porny在线观看| 欧美 亚洲 国产 日韩一| 午夜两性在线视频| 超碰97精品在线观看| 国产野战对白在线观看| 亚洲av美国av| 国内毛片毛片毛片毛片毛片| 久久天躁狠狠躁夜夜2o2o| 99热国产这里只有精品6| 欧美黄色淫秽网站| 欧美久久黑人一区二区| 欧美乱码精品一区二区三区| 变态另类成人亚洲欧美熟女 | 每晚都被弄得嗷嗷叫到高潮| 国产真人三级小视频在线观看| 久久精品国产亚洲av香蕉五月 | 精品人妻在线不人妻| 十八禁网站网址无遮挡| 精品国产乱码久久久久久男人| 一夜夜www| 亚洲一区二区三区欧美精品| 最新在线观看一区二区三区| 黄色怎么调成土黄色| 免费在线观看视频国产中文字幕亚洲| 丰满饥渴人妻一区二区三| 激情视频va一区二区三区| 国产男靠女视频免费网站| 午夜福利一区二区在线看| 日韩三级视频一区二区三区| 亚洲国产欧美一区二区综合| 中文字幕人妻熟女乱码| 欧美日韩亚洲国产一区二区在线观看 | 国产一卡二卡三卡精品| 久久国产精品大桥未久av| 亚洲精华国产精华精| 国产成人av教育| 一区福利在线观看| 美女午夜性视频免费| 免费在线观看完整版高清| 成在线人永久免费视频| 日韩人妻精品一区2区三区| 欧美人与性动交α欧美精品济南到| 一本—道久久a久久精品蜜桃钙片| 亚洲精品国产色婷婷电影| 国产亚洲精品一区二区www | 精品亚洲乱码少妇综合久久| 亚洲,欧美精品.| 91av网站免费观看| 亚洲成a人片在线一区二区| 一级毛片精品| 丰满饥渴人妻一区二区三| 精品国产一区二区三区四区第35| 大香蕉久久网| 丝袜人妻中文字幕| 人妻久久中文字幕网| 在线观看免费高清a一片| 一边摸一边做爽爽视频免费| 国产亚洲精品第一综合不卡| 18禁国产床啪视频网站| 亚洲视频免费观看视频| 国产一区二区三区综合在线观看| 中文字幕人妻丝袜一区二区| 免费高清在线观看日韩| 老司机午夜福利在线观看视频 | 亚洲性夜色夜夜综合| 久久精品国产99精品国产亚洲性色 | a级毛片在线看网站| 99国产精品免费福利视频| 大陆偷拍与自拍| 99热国产这里只有精品6| 久久人妻av系列| 久久 成人 亚洲| 精品国产乱子伦一区二区三区| 91精品国产国语对白视频| 女人爽到高潮嗷嗷叫在线视频| 在线亚洲精品国产二区图片欧美| 51午夜福利影视在线观看| 一级毛片电影观看| av天堂在线播放| 亚洲欧美一区二区三区黑人| 国产成人av激情在线播放| 久久人人97超碰香蕉20202| 肉色欧美久久久久久久蜜桃| 中文字幕最新亚洲高清| 色在线成人网| 后天国语完整版免费观看| 成人精品一区二区免费| 国产精品国产av在线观看| 麻豆国产av国片精品| 一区二区三区乱码不卡18| avwww免费| 欧美日韩视频精品一区| 欧美精品啪啪一区二区三区| 五月开心婷婷网| 自拍欧美九色日韩亚洲蝌蚪91| 中文亚洲av片在线观看爽 | 国产黄色免费在线视频| 搡老熟女国产l中国老女人| 国产不卡一卡二| 一区二区三区精品91| 中文字幕高清在线视频| 少妇被粗大的猛进出69影院| 日韩一卡2卡3卡4卡2021年| 色老头精品视频在线观看| 国产一卡二卡三卡精品| 成人国产av品久久久| av不卡在线播放| 他把我摸到了高潮在线观看 | 免费少妇av软件| 99国产精品一区二区蜜桃av | 纯流量卡能插随身wifi吗| 91av网站免费观看| 99精品在免费线老司机午夜| 精品一品国产午夜福利视频| 无人区码免费观看不卡 | 中文字幕人妻熟女乱码| 九色亚洲精品在线播放| 久久中文字幕一级| 91大片在线观看| 国产精品影院久久| tocl精华| 亚洲色图av天堂| 欧美黑人欧美精品刺激| 国产精品久久电影中文字幕 | 久久国产精品男人的天堂亚洲| 高清毛片免费观看视频网站 | 老司机午夜十八禁免费视频| 不卡一级毛片| 男女边摸边吃奶| 在线观看免费高清a一片| 在线观看免费视频网站a站| 别揉我奶头~嗯~啊~动态视频| 久久精品亚洲熟妇少妇任你| av有码第一页| 黄色丝袜av网址大全| 亚洲伊人色综图| 麻豆av在线久日| 国产成人精品久久二区二区免费| 日本精品一区二区三区蜜桃| 人人妻人人澡人人看| 大片免费播放器 马上看| 丰满少妇做爰视频| 国产男女内射视频| 啦啦啦中文免费视频观看日本| 亚洲精品国产区一区二| 欧美日本中文国产一区发布| 日韩一卡2卡3卡4卡2021年| 精品午夜福利视频在线观看一区 | 建设人人有责人人尽责人人享有的| 99国产精品免费福利视频| 亚洲国产av新网站| 国产亚洲精品一区二区www | 精品一品国产午夜福利视频| 最近最新中文字幕大全免费视频| 久久影院123| 建设人人有责人人尽责人人享有的| 国产精品二区激情视频| 如日韩欧美国产精品一区二区三区| 啦啦啦在线免费观看视频4| 飞空精品影院首页| 欧美激情极品国产一区二区三区| aaaaa片日本免费| 757午夜福利合集在线观看| 国产淫语在线视频| 免费看十八禁软件| 国产精品一区二区免费欧美| 中文亚洲av片在线观看爽 | 日韩视频在线欧美| 蜜桃在线观看..| 欧美日韩亚洲综合一区二区三区_| 国产黄色免费在线视频| 法律面前人人平等表现在哪些方面| 国产精品国产高清国产av | 宅男免费午夜| 亚洲色图 男人天堂 中文字幕| 中文字幕av电影在线播放| 国产深夜福利视频在线观看| 天天添夜夜摸| 一本色道久久久久久精品综合| 日韩免费高清中文字幕av| 少妇精品久久久久久久| cao死你这个sao货| 久久午夜综合久久蜜桃| 精品国产超薄肉色丝袜足j| 老汉色av国产亚洲站长工具| svipshipincom国产片| 精品人妻1区二区| 香蕉国产在线看| 在线十欧美十亚洲十日本专区| 美女高潮到喷水免费观看| 国产片内射在线| 欧美精品高潮呻吟av久久| 日韩有码中文字幕| 欧美日韩黄片免| 国产成人精品无人区| 国产成人免费无遮挡视频| 另类亚洲欧美激情| 欧美另类亚洲清纯唯美| 国产精品麻豆人妻色哟哟久久| 最黄视频免费看| 又紧又爽又黄一区二区| 国产欧美亚洲国产| 亚洲色图 男人天堂 中文字幕| 国产免费现黄频在线看| 日本一区二区免费在线视频| 国产成人av激情在线播放| 国产无遮挡羞羞视频在线观看| 1024香蕉在线观看| 国产福利在线免费观看视频| 国产精品亚洲av一区麻豆| 欧美国产精品va在线观看不卡| 国产精品香港三级国产av潘金莲| 亚洲av欧美aⅴ国产| 国产日韩一区二区三区精品不卡| 亚洲精品在线美女| 香蕉丝袜av| 亚洲欧洲日产国产| 天堂俺去俺来也www色官网| 欧美激情极品国产一区二区三区| 亚洲国产av影院在线观看| 高清欧美精品videossex| 一区福利在线观看| 欧美成狂野欧美在线观看| 一区二区三区国产精品乱码| 日韩成人在线观看一区二区三区| 久久免费观看电影| 黄色a级毛片大全视频| 久久精品aⅴ一区二区三区四区| aaaaa片日本免费| 欧美 日韩 精品 国产| 日韩大码丰满熟妇| 一级a爱视频在线免费观看| 免费女性裸体啪啪无遮挡网站| 99国产精品99久久久久| 日韩熟女老妇一区二区性免费视频| 青青草视频在线视频观看| 一区二区日韩欧美中文字幕| 欧美日韩亚洲国产一区二区在线观看 | 看免费av毛片| 人人妻人人爽人人添夜夜欢视频| 9色porny在线观看| 美女高潮喷水抽搐中文字幕| 国产在线精品亚洲第一网站| www.999成人在线观看| 精品午夜福利视频在线观看一区 | 别揉我奶头~嗯~啊~动态视频| 一边摸一边做爽爽视频免费| 亚洲国产欧美一区二区综合| av免费在线观看网站| 中文字幕人妻丝袜制服| 亚洲av片天天在线观看| h视频一区二区三区| 国产淫语在线视频| 久久毛片免费看一区二区三区| a级毛片在线看网站| √禁漫天堂资源中文www| 国产精品影院久久| 国产精品久久久久久人妻精品电影 | 在线观看免费高清a一片| 在线观看免费视频网站a站| 99久久精品国产亚洲精品| 国产视频一区二区在线看| 五月开心婷婷网| 午夜成年电影在线免费观看| 国产亚洲一区二区精品| 久久精品人人爽人人爽视色| 夜夜夜夜夜久久久久| 亚洲精品在线美女| 日本a在线网址| xxxhd国产人妻xxx| 国产高清激情床上av| 成人黄色视频免费在线看| 久久精品国产综合久久久| 巨乳人妻的诱惑在线观看| a在线观看视频网站| 国精品久久久久久国模美| 少妇精品久久久久久久| 久久性视频一级片| 无限看片的www在线观看| 国产一区二区三区综合在线观看| 热99re8久久精品国产| 一级a爱视频在线免费观看| 19禁男女啪啪无遮挡网站| 久久精品国产99精品国产亚洲性色 | 后天国语完整版免费观看| 国产1区2区3区精品| 国产精品麻豆人妻色哟哟久久| 视频在线观看一区二区三区| 午夜福利视频精品| 午夜福利视频在线观看免费| 免费人妻精品一区二区三区视频| 午夜日韩欧美国产| 国产精品香港三级国产av潘金莲| 人成视频在线观看免费观看| 亚洲色图综合在线观看| 日韩 欧美 亚洲 中文字幕| 一区二区三区精品91| 丰满少妇做爰视频| 亚洲第一青青草原| 国产欧美日韩综合在线一区二区| a级毛片黄视频| 黄色视频不卡| 性少妇av在线| 水蜜桃什么品种好| 国产精品一区二区在线观看99| 国产日韩一区二区三区精品不卡| www.自偷自拍.com| 免费女性裸体啪啪无遮挡网站| 国产人伦9x9x在线观看| 91成人精品电影| 国产精品成人在线| 激情在线观看视频在线高清 | 国产精品一区二区在线观看99| 国产片内射在线| 人人妻,人人澡人人爽秒播| 久久精品91无色码中文字幕| 欧美乱码精品一区二区三区| 久久精品熟女亚洲av麻豆精品| 成人黄色视频免费在线看| 成人精品一区二区免费| 一本一本久久a久久精品综合妖精| 日本撒尿小便嘘嘘汇集6| 精品福利观看| videosex国产| 黑人猛操日本美女一级片| 久久婷婷成人综合色麻豆| 日本a在线网址| 久久国产精品人妻蜜桃| 搡老岳熟女国产| 美女主播在线视频| 亚洲午夜精品一区,二区,三区| 深夜精品福利| svipshipincom国产片| 丝袜在线中文字幕| 69av精品久久久久久 | 国产亚洲欧美在线一区二区| 国产精品久久电影中文字幕 | 岛国在线观看网站| 人人妻人人澡人人爽人人夜夜| 人人澡人人妻人| 欧美av亚洲av综合av国产av| 日本av免费视频播放| 曰老女人黄片| 免费高清在线观看日韩| 人妻久久中文字幕网| 狠狠精品人妻久久久久久综合| 亚洲一区二区三区欧美精品| 久久精品亚洲精品国产色婷小说| 99re6热这里在线精品视频| 国产单亲对白刺激| 欧美日韩亚洲高清精品| 狠狠狠狠99中文字幕| 少妇猛男粗大的猛烈进出视频| 99精品欧美一区二区三区四区| 国产男女内射视频| 免费看a级黄色片| a在线观看视频网站| 黄色毛片三级朝国网站| 亚洲精品一二三| 国产单亲对白刺激| 少妇粗大呻吟视频| 国产精品 国内视频| 精品高清国产在线一区| 极品少妇高潮喷水抽搐| 美女福利国产在线| 国产日韩欧美视频二区| 狠狠狠狠99中文字幕| 丝袜美足系列| 一本一本久久a久久精品综合妖精| 亚洲avbb在线观看| 精品人妻1区二区| 亚洲熟女精品中文字幕| 夜夜骑夜夜射夜夜干| 亚洲男人天堂网一区| 亚洲成国产人片在线观看| 国产成人欧美| 久久亚洲真实| 又黄又粗又硬又大视频| 99国产精品免费福利视频| av网站在线播放免费| 欧美日韩精品网址| 一二三四社区在线视频社区8| 国产单亲对白刺激| 韩国精品一区二区三区| 一区二区三区精品91| 岛国在线观看网站| 天天添夜夜摸| 怎么达到女性高潮| 精品一品国产午夜福利视频| 国产一区二区在线观看av| 超碰97精品在线观看| 在线看a的网站| netflix在线观看网站| 精品少妇黑人巨大在线播放| 美女扒开内裤让男人捅视频| 人人妻,人人澡人人爽秒播| 欧美日韩一级在线毛片| 久久久精品区二区三区| 亚洲av国产av综合av卡| 亚洲午夜精品一区,二区,三区| 国产亚洲精品久久久久5区| 久久国产亚洲av麻豆专区| 男女床上黄色一级片免费看| 亚洲精品国产一区二区精华液| 国产亚洲精品第一综合不卡| 国产极品粉嫩免费观看在线| 国产亚洲一区二区精品| 夜夜爽天天搞| 亚洲va日本ⅴa欧美va伊人久久| 国产一区二区三区在线臀色熟女 | 国产在线一区二区三区精| 中国美女看黄片| 激情在线观看视频在线高清 | 久久天躁狠狠躁夜夜2o2o| 亚洲午夜精品一区,二区,三区| 男女之事视频高清在线观看| 咕卡用的链子| 国产免费现黄频在线看| 一级毛片女人18水好多| 亚洲成av片中文字幕在线观看| 国产精品欧美亚洲77777| 日韩欧美国产一区二区入口| 国产熟女午夜一区二区三区| 国产精品1区2区在线观看. | 人人妻人人爽人人添夜夜欢视频| 多毛熟女@视频| 天天影视国产精品| 69精品国产乱码久久久| 国产亚洲av高清不卡| 亚洲一码二码三码区别大吗| 我要看黄色一级片免费的| 国内毛片毛片毛片毛片毛片| 久久 成人 亚洲| 在线天堂中文资源库| 啦啦啦中文免费视频观看日本| 蜜桃在线观看..| 亚洲精品国产色婷婷电影| 久久久久久免费高清国产稀缺| 高清黄色对白视频在线免费看| 不卡av一区二区三区| 91av网站免费观看| 91九色精品人成在线观看| 极品教师在线免费播放| 久久国产亚洲av麻豆专区| 国产三级黄色录像| 国产成人免费无遮挡视频| 国产欧美亚洲国产| 在线观看免费高清a一片| 欧美日韩中文字幕国产精品一区二区三区 | 精品久久久久久久毛片微露脸| 在线观看免费高清a一片| 在线天堂中文资源库| 啦啦啦中文免费视频观看日本| 99在线人妻在线中文字幕 | 日日夜夜操网爽| 欧美午夜高清在线| 国产高清国产精品国产三级| 久久狼人影院| 国产激情久久老熟女| av天堂在线播放| 桃花免费在线播放| 亚洲伊人色综图| 两人在一起打扑克的视频| 十八禁高潮呻吟视频| 69精品国产乱码久久久| 日本黄色视频三级网站网址 | 在线观看免费日韩欧美大片| 国产视频一区二区在线看| 又大又爽又粗| 日本黄色视频三级网站网址 | 久久久精品国产亚洲av高清涩受| 99国产极品粉嫩在线观看| 亚洲成a人片在线一区二区| 亚洲国产毛片av蜜桃av| 亚洲第一欧美日韩一区二区三区 | 欧美在线一区亚洲| 99国产精品一区二区三区| 女人久久www免费人成看片| 性色av乱码一区二区三区2| 超碰成人久久| 视频区欧美日本亚洲| 久久人妻av系列| 精品一品国产午夜福利视频| 久久久国产一区二区| 国产日韩欧美在线精品| 老司机午夜福利在线观看视频 | 极品教师在线免费播放| 9色porny在线观看| 美女高潮喷水抽搐中文字幕| 脱女人内裤的视频| 精品久久久久久久毛片微露脸| 夜夜骑夜夜射夜夜干| 99国产精品一区二区蜜桃av | 亚洲全国av大片| 女人精品久久久久毛片| 欧美av亚洲av综合av国产av| 国产成人啪精品午夜网站| 欧美激情 高清一区二区三区| 久久青草综合色| 一个人免费在线观看的高清视频| 亚洲自偷自拍图片 自拍| 69av精品久久久久久 | 高清毛片免费观看视频网站 | 亚洲熟女精品中文字幕| 国产深夜福利视频在线观看| 18禁裸乳无遮挡动漫免费视频| 考比视频在线观看| 咕卡用的链子| 亚洲熟女精品中文字幕| 国产主播在线观看一区二区| 久久精品国产综合久久久| 两性夫妻黄色片| 国产又色又爽无遮挡免费看| 成人精品一区二区免费| 高清毛片免费观看视频网站 | 老熟妇乱子伦视频在线观看| 精品一区二区三区四区五区乱码| 91麻豆精品激情在线观看国产 | 9191精品国产免费久久| 国产精品成人在线| 国产精品香港三级国产av潘金莲| 午夜日韩欧美国产| 欧美乱妇无乱码| 久久亚洲精品不卡| 一区二区三区乱码不卡18| 久9热在线精品视频| 久久精品国产亚洲av香蕉五月 | 巨乳人妻的诱惑在线观看| 国产1区2区3区精品| 高清毛片免费观看视频网站 | 欧美精品亚洲一区二区| e午夜精品久久久久久久| 亚洲精品美女久久久久99蜜臀| 老司机在亚洲福利影院| 亚洲人成伊人成综合网2020| 亚洲成av片中文字幕在线观看| 9热在线视频观看99| 日本撒尿小便嘘嘘汇集6| 精品人妻熟女毛片av久久网站| 国产精品1区2区在线观看. | 国产aⅴ精品一区二区三区波| 亚洲欧美一区二区三区黑人| 午夜精品国产一区二区电影| 制服人妻中文乱码| 男女下面插进去视频免费观看| 精品亚洲乱码少妇综合久久| 久久青草综合色|