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

帶浮力塊立管的橫流-順流雙向渦激振動響應(yīng)研究

2022-07-14 13:17:22莊陳怡李宇坤袁昱超薛鴻祥唐文勇
振動與沖擊 2022年13期
關(guān)鍵詞:模態(tài)振動

莊陳怡, 李宇坤, 袁昱超, 薛鴻祥, 唐文勇

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

隨著海洋油氣開發(fā)向深水領(lǐng)域進(jìn)軍,立管需要足夠大的頂部張力以維持自身穩(wěn)定性,然而由于頂部張力器附近立管材料強(qiáng)度和接頭強(qiáng)度的限制,必須設(shè)置浮力塊對立管提供分布式浮力,以減小系統(tǒng)濕重、降低對立管張力器張力的要求。浮力塊的增設(shè)使頂張式立管的張力及重量分布發(fā)生變化,導(dǎo)致固有頻率、振型、曲率等模態(tài)特性發(fā)生改變。此外浮力塊將增大立管的水動力外徑,這將造成兩方面的影響:對于能量輸入?yún)^(qū),較高的水動力外徑會增加流體向立管結(jié)構(gòu)的能量輸入,增大激勵力;對于能量輸出區(qū),會增加流體阻尼。因此,帶浮力塊立管的動力響應(yīng)特性與裸管相比存在明顯差異。

目前,國內(nèi)外對于立管渦激振動(vortex-induced vibration, VIV)機(jī)理已有較為深入的研究[1-3],但對于布置有浮力塊裝置的立管渦激振動特性的研究較有限,主要研究方法分為模型試驗和數(shù)值模擬。Lie等[4]在挪威MARINTEK拖曳水池中,對均勻流和線性剪切流下4種不同浮力塊覆蓋率的立管振動進(jìn)行試驗,分析了立管橫流向的試驗數(shù)據(jù),評估交錯布置的浮力塊對船用立管渦激振動的影響。SHELL Exploration and Production在MARINTEK測試了均勻流下5種不同浮力塊交錯布置方案立管模型的VIV響應(yīng)情況,Rao等[5]針對“何種浮力塊覆蓋率下,浮力塊引起的激勵頻率主導(dǎo)立管的VIV響應(yīng)”這一問題,分析了橫流振動試驗數(shù)據(jù)。挪威深水計劃在MARINTEK 對6種不同浮力塊交錯布置立管進(jìn)行水動力模型試驗,Wu等[6]對橫流向試驗數(shù)據(jù)進(jìn)行分析,總結(jié)了相鄰浮力塊段與裸管段長度比、浮力塊長徑比、浮力塊覆蓋率對于立管渦激振動的影響。數(shù)值模擬方面,Wu等[7-8]以SHELL Exploration and Production在MARINTEK開展的交錯浮力塊立管試驗為數(shù)據(jù)基準(zhǔn),采用VIVANA軟件預(yù)報了浮力塊立管橫流向的渦激振動響應(yīng),基于試驗調(diào)整了激勵力系數(shù)庫,并與試驗結(jié)果進(jìn)行對比,針對另一較小尺度立管進(jìn)行了類似工作。孫友義等[9]使用SHEAR7程序?qū)в懈×K的隔水管進(jìn)行渦激振動分析,研究不同浮力塊分布形式對隔水管系統(tǒng)渦激振動疲勞特性的影響,對浮力塊分布方案進(jìn)行優(yōu)化。工程上及相關(guān)研究表明[10-12],合理的浮力塊布置可減小立管的疲勞損傷,從而延長立管的疲勞壽命,提高立管結(jié)構(gòu)安全性。

對帶有浮力塊裝置柔性立管渦激振動的研究具有工程價值,而針對此種立管渦激振動研究的有限性使其響應(yīng)特性及機(jī)理尚未清晰可知,且現(xiàn)存的渦激振動預(yù)報模型大多以橫流向渦激振動為主要研究目標(biāo)[13],忽略了橫流向與順流向振動間明顯存在的耦合效應(yīng)[14]。為了與實際工程更為緊密的結(jié)合,有必要針對布置有浮力塊的立管渦激振動建立有效的雙向響應(yīng)預(yù)報模型以深入研究其渦激振動特性。本文針對帶有浮力塊的頂張式立管,基于雙自由度受迫振動試驗數(shù)據(jù),采用有限元法和能量平衡方程建立了考慮雙向耦合效應(yīng)的渦激振動頻域預(yù)報模型,分析不同布置方式浮力塊立管的渦激振動響應(yīng),研究浮力塊布置的優(yōu)化方案。

1 立管雙向渦激振動頻域預(yù)報模型

1.1 立管振動方程

將頂張式立管簡化為受軸向張力作用的簡支梁,坐標(biāo)原點位于立管懸掛點,z軸向下為正,x軸為順流方向,y軸為橫流方向。將簡支梁在空間上進(jìn)行有限元離散,得到單元剛度矩陣K=KM+KT,其中:

(1)

(2)

梁單元集中質(zhì)量矩陣可寫為

(3)

式中:KM與KT分別為立管單元抗彎剛度矩陣和軸向張力提供的附加剛度矩陣;EI為彎曲剛度,由于浮力塊裝置的剛度極小,本文忽略浮力塊對立管彎曲剛度的影響,認(rèn)為彎曲剛度沿立管均勻分布;T為軸向張力,浮力塊的布置將改變立管各位置處的張力分布;l為單元長度;m0為立管結(jié)構(gòu)單位長度質(zhì)量,浮力塊將改變立管結(jié)構(gòu)的質(zhì)量分布;ma=CaρπD2/2為單位長度附加質(zhì)量;Ca為附加質(zhì)量系數(shù);ρ為流體密度;D為立管直徑,對于帶有浮力塊的立管段即為浮力塊直徑。

引入簡支梁頂端和底端的約束條件,按照單元節(jié)點坐標(biāo)組裝得到總剛度矩陣[KM+KT]和質(zhì)量矩陣[M]后,可得振動方程

[KM+KT]φ-ω2[M]φ=0

(4)

式中:ω為圓頻率,f=2π/ω即為固有頻率;φ為模態(tài)振型。通過求解矩陣方程即可得到簡支梁不同模態(tài)的頻率和振型。

圖1 橫流方向附加質(zhì)量系數(shù)簡化曲線Fig.1 Simplified curve of added mass coefficient in CF direction

圖2 順流方向附加質(zhì)量系數(shù)簡化曲線Fig.2 Simplified curve of added mass coefficient in IL direction

至此,可進(jìn)行立管固有頻率及振型的計算。考慮浮力塊的影響,計算中D為立管各處實際水動力直徑。橫流方向上,立管響應(yīng)頻率及振型通過迭代計算得到:

(1) 假設(shè)Ca,CF=1,求解振動方程得到立管固有頻率fCF,0;

(3) 根據(jù)橫流方向附加質(zhì)量系數(shù)簡化曲線得到附加質(zhì)量分布,更新振動方程,計算得到新的振動頻率fCF,k+1;

(4) 重復(fù)步驟(2)、(3)直至結(jié)果收斂,得到橫流向?qū)嶋H振動頻率fCF和振型φCF(z)。

由于基于順流方向附加質(zhì)量系數(shù)計算得到的振動頻率并不總是滿足為橫流響應(yīng)頻率的兩倍,常出現(xiàn)2fi,CF落于fj,IL、fj+1,IL兩者之間的情況,因此需要對順流附加質(zhì)量進(jìn)行調(diào)整。在本文模型中通過減小順流向附加質(zhì)量的方式增大固有頻率,以滿足fj,IL=2fi,CF。順流方向上立管響應(yīng)頻率及振型計算過程如下:

(2) 根據(jù)附加質(zhì)量系數(shù)分布曲線得到順流向附加質(zhì)量分布,建立振動方程,計算第i階橫流向振動模態(tài)對應(yīng)的順流向模態(tài)及頻率;

(3) 對比fj,IL和fi,CF,調(diào)整順流向附加質(zhì)量,直至滿足fj,IL=2fi,CF,并得到實際順流向振型φj,IL(z)。

1.2 渦激振動流體力模型

1.2.1 激勵力模型

采用Zheng試驗得到的橫流向和順流向渦激振動激勵力系數(shù)計算流體對立管的激勵力。與單自由度強(qiáng)迫振動試驗不同,Zheng所進(jìn)行的試驗反映出了橫流與順流振動間存在的耦合效應(yīng),并在激勵力系數(shù)上體現(xiàn)。以0°位移相位角為例,典型激勵力系數(shù)分布云圖如圖3和圖4所示。可以發(fā)現(xiàn),隨順流振幅增大,橫流向激勵中心所對應(yīng)無因次頻率減小;而順流向激勵力系數(shù)將隨橫流振幅呈現(xiàn)增大趨勢。參考文獻(xiàn)[14],將立管橫流向無因次頻率介于0.125~0.25之間的區(qū)域作為能量輸入?yún)^(qū)。根據(jù)無因次響應(yīng)頻率及幅值進(jìn)行插值可獲得相應(yīng)的橫流及順流向激勵力系數(shù),本文通過改變水動力直徑考慮浮力塊引起的激勵力系數(shù)及激勵力變化。

1.2.2 阻尼模型

當(dāng)立管結(jié)構(gòu)響應(yīng)超出能量輸入?yún)^(qū)范圍時,認(rèn)為流體作用力將轉(zhuǎn)變?yōu)樗畡恿ψ枘帷1疚氖褂肰enugopal[16]針對渦激振動提出的經(jīng)驗阻尼模型,該模型分為高流速區(qū)和低流速區(qū)。

高流速區(qū)水動力阻尼系數(shù)為

rhh=ρU2(z)Cvh/ω

(5)

低流速區(qū)水動力阻尼系數(shù)為

rhl=rsw(z)+ρDU(z)Cvl

(6)

1.3 渦激振動響應(yīng)幅值計算

參考SHEAR7方法,通過預(yù)估能量選取振動激發(fā)模態(tài)[17]。針對橫流向振動模態(tài)i及對應(yīng)順流向振動模態(tài)j,計算預(yù)估能量

(7)

其中,

(8)

(9)

選取激發(fā)振動模態(tài)時,根據(jù)最大預(yù)估能量值Πmax,計算各橫流-順流振動模態(tài)組合預(yù)估能量比Πi&j/Πmax,以0.7為能量閾值,選取橫流及順流向渦激振動激發(fā)模態(tài)。

其中,模態(tài)激勵力為

(10)

式中:Y(z)=Ay|φCF(z)|/D為立管各位置橫流向振幅;X(z)=Ax|φIL(z)|/D為立管各位置順流向振幅。針對裸管段和浮力塊段,通過改變直徑計算各段的激勵力和預(yù)估能量。

立管模態(tài)速度為

Vi(t)=Aiωisin(ωit)

(11)

模態(tài)阻尼為

(12)

針對渦激振動激發(fā)的橫流向模態(tài)i及順流向模態(tài)j建立能量平衡方程,可得模態(tài)振幅

(13)

(14)

同時,本文通過模態(tài)疊加法處理柔性立管多階模態(tài)同時被激發(fā)的現(xiàn)象。綜上,頂張式立管渦激振動響應(yīng)計算流程如圖5所示。此前已有工作驗證了此預(yù)報模型用于預(yù)報裸露立管渦激振動響應(yīng)的準(zhǔn)確性[18]。

圖5 渦激振動預(yù)報流程圖Fig.5 The flow chart of VIV prediction

1.4 疲勞損傷計算

立管渦激振動使結(jié)構(gòu)受交變彎矩作用。根據(jù)平斷面假設(shè),第i階振動模態(tài)引起彎矩為

(15)

對應(yīng)產(chǎn)生的最大彎曲應(yīng)力為

(16)

根據(jù)Miner線性累計假設(shè)及S-N曲線,可以得到第i階模態(tài)下渦激振動引起的疲勞損傷分布

(17)

式中:fi為渦激振動響應(yīng)頻率;TD為立管使用時間;A和m為S-N曲線參數(shù);SCF為應(yīng)力集中系數(shù)。

由于頂張式立管發(fā)生渦激振動響應(yīng)時往往存在多種模態(tài)成分,疲勞損傷由各模態(tài)共同產(chǎn)生。本文使用SRSS(square root of the squre stress)方法計算渦激振動疲勞損傷

(18)

2 不同浮力塊覆蓋率下立管渦激振動響應(yīng)

目前工程上認(rèn)為浮力塊交錯分布有助于降低渦激振動對立管造成的疲勞損傷[19-20]。本節(jié)選取一實尺度立管作為算例,研究25%、50%、100%浮力塊覆蓋率下立管渦激振動響應(yīng),并與裸管進(jìn)行比較。浮力塊的形狀為有一定厚度的環(huán)形塊體,包裹于立管外圍,浮力塊模型及四種覆蓋率下浮力塊分布如圖6所示,其中浮力塊間隔均布設(shè)置:25%覆蓋率情況下每隔30 m布置10 m長浮力塊,50%覆蓋率情況下每隔20 m布置20 m長浮力塊。

圖6 浮力塊模型及不同浮力塊覆蓋率下立管模型Fig.6 Buoyancy module model and riser models with different buoyancy module coverage

該算例立管工作水深1 000 m,立管外徑為0.508 m,壁厚為0.025 4 m,單位長度質(zhì)量為536.46 kg/m,彈性模量為2.1×1011N/m2,浮力塊外徑為1.193 8 m,單位長度質(zhì)量為258.73 kg/m,提供凈浮力為312.27 kg/m,立管底端張力取0.5倍裸管濕重。海流為線性剪切流,海底速度為0,海面流速為0.4~2.0 m/s。疲勞損傷計算選取DNV給出的有陰極保護(hù)的海洋結(jié)構(gòu)物S-N曲線[21]。

靜水中立管前20階固有頻率見圖7。在設(shè)置浮力塊后由于立管整體質(zhì)量增大、張力減小,立管的固有頻率降低。其中在立管全部覆蓋浮力塊的情況下,頻率約為裸管的50%。

圖7 不同浮力塊覆蓋率下立管固有頻率Fig.7 Natural frequencies of riser with different buoyancy module coverage rates

圖8為四種浮力塊覆蓋率立管模型在不同流速剪切流作用下的渦激振動響應(yīng)頻率。可以發(fā)現(xiàn),四種覆蓋率下立管模型渦激振動響應(yīng)頻率明顯落于兩條趨勢線附近,分別對應(yīng)于立管和浮力塊外徑的斯托哈爾頻率。由于浮力塊的外徑遠(yuǎn)大于立管直徑,當(dāng)渦激振動由浮力塊所控制時響應(yīng)頻率明顯低于立管控制情況。對于50%覆蓋率立管模型,立管響應(yīng)完全由浮力塊所控制,響應(yīng)頻率落于斜率較小的趨勢線附近。而對于25%覆蓋率的立管模型,響應(yīng)頻率離散分布在兩條趨勢線附近,表明在該覆蓋率下渦激振動響應(yīng)將由兩種水動力外徑共同控制,即流體通過浮力塊向結(jié)構(gòu)輸入的能量與通過立管向結(jié)構(gòu)輸入的能量相當(dāng)。

圖9為四種浮力塊覆蓋率立管模型在不同流速下的渦激振動最大疲勞損傷。全部覆蓋浮力塊時立管最大疲勞損傷略高于裸管情況,并且在順流方向上較為明顯。而25%及50%覆蓋率情況下立管最大疲勞損傷較低。25%覆蓋率立管模型的響應(yīng)頻率由裸管與浮力塊共同控制,從而使模型不易發(fā)生渦激振動,具有較好的疲勞特性;值得注意的是,25%浮力塊覆蓋率下,損傷在1.4~1.6 m/s的流速范圍內(nèi)產(chǎn)生了隨流速增大而減小的變化,這與激發(fā)模態(tài)有關(guān),從圖8可以看到,1.4 m/s流速下立管只有一個低激發(fā)頻率,而1.6 m/s流速下存在兩個低激發(fā)頻率,后者更趨向于造成較小的疲勞損傷。50%覆蓋率模型的響應(yīng)頻率由浮力塊控制,因浮力塊交錯分布擾亂了沿立管軸向漩渦泄放的一致性,從而降低了漩渦對系統(tǒng)的能量輸入,因此50%覆蓋模型的疲勞特性要明顯優(yōu)于100%覆蓋模型。可見對于頂張式立管,布置一定覆蓋率的浮力塊在減少頂部張力需求的同時可有效降低由渦激振動引起的立管疲勞損傷。

以海面流速為1.2 m/s工況為例對立管渦激振動響應(yīng)進(jìn)行具體分析,響應(yīng)位移沿立管分布如圖10所示。其中,25%覆蓋率情況下渦激振動響應(yīng)由立管及浮力塊共同控制,使得“頻率鎖定”現(xiàn)象不易發(fā)生,響應(yīng)位移僅為裸管的1/3。50%覆蓋率情況下渦激振動響應(yīng)由浮力塊控制,浮力塊間隔布置削弱了瀉渦強(qiáng)度并使結(jié)構(gòu)中輸入能量減少,且響應(yīng)模態(tài)階次降低,其響應(yīng)位移與裸管情況下較為接近。而在全部覆蓋浮力塊后,固有頻率的降低使響應(yīng)模態(tài)與裸管相比沒有發(fā)生明顯變化,而水動力直徑的增大令立管響應(yīng)位移大幅增加。

立管疲勞損傷沿管長分布如圖11所示。100%覆蓋率情況下立管的疲勞損傷最高,全浮力塊分布方案造成立管水動力外徑增大,外部流體對模型的輸入能量大幅提高,導(dǎo)致疲勞損傷增加,尤其是在高流速區(qū)可達(dá)到裸管疲勞損傷的1 000倍。50%覆蓋率下立管由于振動頻率及模態(tài)階次較低,疲勞損傷相對較小。而對于25%情況,雖然其響應(yīng)位移為四種浮力塊覆蓋率情況下最小的,但是渦激振動響應(yīng)中存在的高階模態(tài)成分導(dǎo)致疲勞損傷與50%覆蓋率下立管處在同一水平上。同一流速工況下,不同浮力塊布置方案的立管模型其響應(yīng)位移與疲勞損傷呈現(xiàn)非一致性,因此以抑制疲勞損傷為目標(biāo)進(jìn)行浮力塊布置方案優(yōu)化時,還需兼顧考慮立管位移響應(yīng)幅值。

3 浮力塊布置方案優(yōu)化

為研究不同浮力塊布置方式對頂張式立管渦激振動的影響,針對25%浮力塊覆蓋率的立管模型,對比圖12所示四種布置方案下渦激振動響應(yīng)特性。其中,方案1浮力塊均勻間隔布置,方案2布置在高流速區(qū),方案3布置在立管中段,方案4布置在低流速區(qū)。

圖12 不同浮力塊布置方案的立管模型Fig.12 Riser models with different buoyancy module arrangements

圖13~15分別為四種浮力塊布置方案立管在不同流速下的渦激振動響應(yīng)頻率、響應(yīng)位移及最大疲勞損傷。在不同流速工況下同一浮力塊布置方案將在響應(yīng)頻率上呈現(xiàn)相同的特性:其中方案1的渦激振動響應(yīng)由立管和浮力塊共同控制;方案2的響應(yīng)由浮力塊控制;方案3和方案4則由立管控制。分析表明:立管響應(yīng)頻率由處于高流速區(qū)的模型水動力直徑?jīng)Q定;處于低流速區(qū)的浮力塊對響應(yīng)頻率起不到控制作用,僅可對響應(yīng)位移產(chǎn)生影響。

方案1的響應(yīng)頻率由浮力塊與立管共同控制,不易發(fā)生渦激共振,且浮力塊的交錯分布降低了旋渦強(qiáng)度,疲勞損傷特性整體較優(yōu),同時響應(yīng)位移最小。方案2立管的渦激振動響應(yīng)由浮力塊控制,響應(yīng)模態(tài)階次及響應(yīng)頻率相對較低,疲勞損傷與方案1較為接近,具有較好的疲勞特性,但是在該方案中由于浮力塊置于高流速區(qū)中,流體將通過浮力塊向結(jié)構(gòu)中輸入大量的能量,令立管產(chǎn)生遠(yuǎn)大于其他三種方案的位移響應(yīng)。方案3和方案4的渦激振動響應(yīng)由立管控制,響應(yīng)模態(tài)階次高、頻率大,導(dǎo)致疲勞損傷程度整體較高;尤其在方案4中,立管損傷遠(yuǎn)超其他布置方案,同時由于處于激勵狀態(tài)的裸管段處于中高流速區(qū),結(jié)構(gòu)的能量輸入較大,而浮力塊位于低流速區(qū),其產(chǎn)生的流體阻尼效應(yīng)較小,導(dǎo)致立管整體出現(xiàn)較大的響應(yīng)位移。因此,從整體看,方案1為較優(yōu)的浮力塊布置方案。

4 結(jié) 論

基于本文研究,得出如下結(jié)論:

(1) 布置浮力塊將明顯改變立管渦激振動特性。立管全部覆蓋浮力塊將加劇響應(yīng)位移與疲勞損傷;而在25%及50%覆蓋率下間隔布置浮力塊會影響流體瀉渦沿管長分布的一致性,使得流體向結(jié)構(gòu)輸入的能量較少,導(dǎo)致立管渦激振動響應(yīng)幅值相對較小,并有效降低渦激振動引起的疲勞損傷。

(2) 從降低渦激振動疲勞損傷的角度出發(fā),對頂張式立管設(shè)置浮力塊時,應(yīng)盡量使其動力響應(yīng)由浮力塊和裸管段共同控制,以抑制渦激振動。同時還可在減小高流速區(qū)內(nèi)浮力塊分布長度的前提下,令渦激振動響應(yīng)由浮力塊單獨控制,這樣既能減小響應(yīng)位移,又可降低立管的振動模態(tài)階次及響應(yīng)頻率,改善結(jié)構(gòu)疲勞性能。

(3) 對于本文研究的25%浮力塊覆蓋率情況,推薦將浮力塊間隔布置,此時立管具有較好的疲勞特性且響應(yīng)位移最小。

猜你喜歡
模態(tài)振動
振動的思考
噴水推進(jìn)高速艇尾部振動響應(yīng)分析
This “Singing Highway”plays music
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
車輛CAE分析中自由模態(tài)和約束模態(tài)的應(yīng)用與對比
國內(nèi)多模態(tài)教學(xué)研究回顧與展望
高速顫振模型設(shè)計中顫振主要模態(tài)的判斷
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
UF6振動激發(fā)態(tài)分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 大陆国产精品视频| 亚洲综合香蕉| 国产亚洲精品自在久久不卡| 免费观看国产小粉嫩喷水| 久久精品免费国产大片| 久久精品国产999大香线焦| 欧美国产成人在线| a在线观看免费| 亚洲日本精品一区二区| 国产精品欧美在线观看| 欧美日韩精品一区二区视频| 在线亚洲精品自拍| 韩日无码在线不卡| 欧洲欧美人成免费全部视频| 97av视频在线观看| 婷婷色中文网| 亚洲视频免| 99re66精品视频在线观看| 亚洲免费三区| 精品无码人妻一区二区| 亚洲欧美成人| 亚洲精品免费网站| 久久精品一卡日本电影| 欧美日韩国产成人高清视频| 尤物午夜福利视频| 亚洲国产91人成在线| 日本午夜视频在线观看| 国精品91人妻无码一区二区三区| 操美女免费网站| 亚洲一级无毛片无码在线免费视频| 又黄又爽视频好爽视频| 久久精品丝袜| 视频一区视频二区日韩专区| 人与鲁专区| 国产精品流白浆在线观看| 久久先锋资源| 久久人人爽人人爽人人片aV东京热| 久久久久国产精品嫩草影院| 老司机久久99久久精品播放| 亚洲AV无码一二区三区在线播放| 97在线国产视频| 亚洲精品无码AⅤ片青青在线观看| a网站在线观看| 91黄视频在线观看| 无码AV日韩一二三区| 婷婷色中文| 亚洲国产日韩在线成人蜜芽| 四虎永久在线视频| 国产精品一区在线麻豆| 一级香蕉视频在线观看| 国产永久无码观看在线| 国产在线拍偷自揄拍精品| 日韩欧美国产另类| 中文成人在线视频| 曰AV在线无码| 欧美精品伊人久久| 亚洲黄色片免费看| 在线高清亚洲精品二区| 亚洲视频a| 国产精品视频猛进猛出| 国产成人欧美| 色偷偷综合网| 国产精品亚洲一区二区三区z| 久久中文字幕2021精品| 午夜精品福利影院| 91久久精品日日躁夜夜躁欧美| 中文字幕人妻av一区二区| 国产特级毛片aaaaaa| 日韩无码一二三区| 国产视频久久久久| 久久黄色一级片| 影音先锋丝袜制服| 91精品国产91久无码网站| 91区国产福利在线观看午夜| 亚洲AV色香蕉一区二区| 久久黄色视频影| 久久国产精品电影| 日韩在线观看网站| 精品国产成人av免费| 91啪在线| h网址在线观看| 无码AV动漫|