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

FPSO系泊系統(tǒng)數(shù)值分析方法對比研究

2023-10-25 11:41:56姜季江付世曉許玉旺
船舶力學(xué) 2023年10期
關(guān)鍵詞:系統(tǒng)

趙 冰,姜季江,付世曉,許玉旺

(1.上海交通大學(xué)a.海洋工程國家重點實驗室;b.船舶海洋與建筑工程學(xué)院,上海 200240;2.中遠(yuǎn)海運(yùn)重工有限公司,上海 200135)

0 引 言

隨著海洋油氣資源開發(fā)走向深遠(yuǎn)海,惡劣的作業(yè)海域環(huán)境對浮式結(jié)構(gòu)物關(guān)鍵系統(tǒng)的研制提出了更高的要求。通用型FPSO 是一種適用于多種海域環(huán)境的海洋油氣開發(fā)平臺,搭載多點、單點內(nèi)轉(zhuǎn)塔及外轉(zhuǎn)塔等多種系泊系統(tǒng),工作水深可達(dá)2000 m 左右。作為保證上部浮體正常定位作業(yè)的關(guān)鍵裝備,系泊系統(tǒng)的設(shè)計分析對整個平臺的安全性能尤為重要[1]。

通用型FPSO 系泊系統(tǒng)數(shù)值分析包括船體與系泊系統(tǒng)兩部分。其中,船體作為剛體,受到風(fēng)浪流等環(huán)境載荷作用,產(chǎn)生6自由度剛體運(yùn)動。通常波浪載荷基于勢流理論采用面元法進(jìn)行頻域計算,然后轉(zhuǎn)至?xí)r域;而風(fēng)流載荷需要開展模型試驗獲得,在初步設(shè)計階段可采用已發(fā)表的近似船型風(fēng)流載荷系數(shù)進(jìn)行估算。系泊系統(tǒng)的分析主要包括解析方法和數(shù)值離散方法。解析懸鏈線模型忽略動態(tài)流力作用,僅考慮重力及浮力作用,且假定剛性海底無摩擦作用,從而得到船體偏移與系泊線頂端張力之間的非線性函數(shù)關(guān)系,但是這種方法只能用于系泊系統(tǒng)的靜力計算分析[2]。另一種是離散數(shù)值模型,將系泊線視為彈性連續(xù)體,通常采用集中質(zhì)量法、有限單元法或有限差分法等進(jìn)行結(jié)構(gòu)離散。Hall和Goupee[3]采用集中質(zhì)量法對系泊線進(jìn)行了建模分析,并忽略了彎曲及扭轉(zhuǎn)剛度;Tang等[4]同樣采用集中質(zhì)量法對網(wǎng)箱的網(wǎng)衣與系泊線的動力響應(yīng)進(jìn)行了計算分析;Guo等[5]采用有限元法對浮式風(fēng)機(jī)的懸鏈?zhǔn)较挡淳€回復(fù)剛度的遲滯特性進(jìn)行了數(shù)值分析與研究;Kim 等[6]分別采用線性彈簧模型與非線性有限元模型模擬系泊系統(tǒng)對浮體運(yùn)動響應(yīng)進(jìn)行分析,其中線性彈簧模型由于無法計及非線性因素導(dǎo)致結(jié)果不夠準(zhǔn)確。離散數(shù)值模型能夠計及節(jié)點運(yùn)動產(chǎn)生的修正莫里森流載荷、海洋生物附著、海土剛度及摩擦作用等非線性因素,這種方法可以用于系泊系統(tǒng)的非線性時域動力計算。

真實海洋環(huán)境中,船體與系泊系統(tǒng)的動力響應(yīng)存在強(qiáng)烈的非線性耦合作用。在時域計算中,通常采用不同的數(shù)值分析方法來考慮二者不同程度的耦合機(jī)制,數(shù)值分析方法主要包括準(zhǔn)靜態(tài)分析方法、解耦動力分析方法以及耦合動力分析方法[7]。準(zhǔn)靜態(tài)分析方法以懸鏈線理論為基礎(chǔ),將上一時刻船體偏移位置處的系泊系統(tǒng)張力作為外力輸入,計算當(dāng)前時刻船體的偏移位置,從而得到船體6自由度運(yùn)動響應(yīng)以及對應(yīng)的系泊水平張力[2];解耦動力分析方法把系泊系統(tǒng)作為彈性體進(jìn)行結(jié)構(gòu)離散,將準(zhǔn)靜態(tài)方法獲得的船體運(yùn)動響應(yīng)作為系泊系統(tǒng)的頂端運(yùn)動邊界條件,計及波流作用下的莫里森載荷等非線性因素,開展系泊系統(tǒng)的時域動力計算[8-9];耦合動力分析方法將船體與系泊系統(tǒng)作為統(tǒng)一整體,建立全耦合控制方程組進(jìn)行求解,這種方法由于耦合程度較高,計算結(jié)果準(zhǔn)確性較高,但往往耗費大量計算時間[10-13]。此外,還有學(xué)者采用頻域分析方法對浮式結(jié)構(gòu)物的系泊系統(tǒng)進(jìn)行耦合分析[14],但是頻域方法難以考慮所有的非線性影響因素。

本文以通用型FPSO 多點系泊系統(tǒng)為研究對象,基于準(zhǔn)靜態(tài)分析、解耦動力分析以及耦合動力分析等數(shù)值方法,開展系泊系統(tǒng)總體動力響應(yīng)計算分析,通過分析系泊張力和船體運(yùn)動響應(yīng)等結(jié)果,對比三種數(shù)值方法的計算精度與效率,為通用型FPSO系泊系統(tǒng)不同設(shè)計階段的數(shù)值分析方法提供指導(dǎo)性建議。

1 數(shù)值分析方法

通用型FPSO的系泊系統(tǒng)數(shù)值分析方法主要包括準(zhǔn)靜態(tài)方法、解耦動力分析方法和耦合動力分析方法。其中,準(zhǔn)靜態(tài)方法采用懸鏈線方程模擬系泊系統(tǒng),解耦動力方法將準(zhǔn)靜態(tài)方法獲得的船體運(yùn)動響應(yīng)作為系泊系統(tǒng)有限元模型的頂端運(yùn)動邊界條件,耦合動力方法將船體運(yùn)動方程與系泊系統(tǒng)有限元動力學(xué)控制方程聯(lián)立求解。

1.1 準(zhǔn)靜態(tài)分析方法

如圖1所示,準(zhǔn)靜態(tài)方法將系泊系統(tǒng)采用解析懸鏈線模型進(jìn)行分析,獲得頂部導(dǎo)纜點不同偏移位置處的系泊張力,為船體運(yùn)動提供系泊回復(fù)作用力-[]Km{ }x。此外,系泊系統(tǒng)的動態(tài)響應(yīng)對船體運(yùn)動還存在一種阻尼作用,在準(zhǔn)靜態(tài)分析方法中將其考慮為線性阻尼-[Cm]{x?}。因此,在風(fēng)浪流等載荷作用下,系泊船體動力學(xué)控制方程為

圖1 準(zhǔn)靜態(tài)分析方法Fig.1 Quasi-static analysis method

1.1.1 準(zhǔn)靜態(tài)系泊剛度

在準(zhǔn)靜態(tài)分析方法中,系泊系統(tǒng)為船體提供的回復(fù)作用力可通過懸鏈線理論獲得。通常單一成分的懸鏈線無法滿足工程需求,系泊線由多種成分組成,對于每一種成分,均可采用懸鏈線方程模擬,在連接點處滿足張力連續(xù),形成單根多成分系泊線的懸鏈線方程[15]。

準(zhǔn)靜態(tài)系泊剛度是系泊系統(tǒng)回復(fù)力對船體位移的偏導(dǎo)數(shù)。根據(jù)懸鏈線理論計算每根系泊線在不同導(dǎo)纜點處的系泊張力,利用剛性船體參考點與導(dǎo)纜點的空間位置關(guān)系,獲得系泊系統(tǒng)在船體參考點處的合力,其中合力水平分量對船體水平位移的偏導(dǎo)數(shù)為系泊系統(tǒng)的水平面剛度。圖2給出了500 m 水深條件下一根長1200 m、直徑111 mm、單位長度濕重234.4 kg/m 的錨鏈系泊張力水平分量與上端點水平偏移的非線性關(guān)系曲線,可以看出懸鏈線理論結(jié)果與商業(yè)軟件Riflex和OrcaFlex 的離散數(shù)值模型靜力分析結(jié)果基本相同,其中Riflex 采用有限單元法,OrcaFlex 采用集中質(zhì)量法,說明懸鏈線理論可以獲得較為準(zhǔn)確的準(zhǔn)靜態(tài)系泊剛度。

1.1.2 準(zhǔn)靜態(tài)系泊阻尼

系泊系統(tǒng)對船體運(yùn)動的準(zhǔn)靜態(tài)阻尼可以表示為線性阻尼,主要影響船體的水平面運(yùn)動自由度(縱蕩、橫蕩及首搖),阻尼系數(shù)陣為

圖2 系泊張力水平分量-船體水平偏移關(guān)系曲線Fig.2 Horizontal component of mooring tension-ship offset relationship curve

其中,各個系數(shù)之間滿足關(guān)系

式中,L為船長;Km,11為船體在初始位置處的縱蕩系泊剛度;ζ為臨界阻尼比,可采用衰減運(yùn)動數(shù)值分析方法獲得(見4.2節(jié))。

1.2 解耦動力分析方法

如圖3所示,解耦動力分析方法在準(zhǔn)靜態(tài)方法的基礎(chǔ)上,將船體六自由度運(yùn)動時歷作為系泊系統(tǒng)有限元模型的頂端運(yùn)動邊界條件,同時計及莫里森流載荷、海底剛度及摩擦作用等非線性因素,開展系泊系統(tǒng)的動力計算,獲得系泊線各個節(jié)點的動力響應(yīng)。

圖3 解耦動力分析方法Fig.3 Uncoupled dynamic analysis method

系泊系統(tǒng)的動力學(xué)控制方程為

式中,[m]為系泊線質(zhì)量矩陣,[c]為阻尼矩陣,[k]為剛度矩陣,{ }u、{ }u?、{ }u? 分別為系泊線節(jié)點位移、速度及加速度列向量,{ }f為系泊線受到的莫里森載荷、海土的支持力及摩擦力等合外力列向量。系泊線頂端節(jié)點運(yùn)動邊界條件{uline-top} 為船體導(dǎo)纜點處的準(zhǔn)靜態(tài)方法獲得的運(yùn)動響應(yīng){xship-fairlead} 。

1.3 耦合動力分析方法

如圖4 所示,耦合動力分析方法將船體動力學(xué)控制方程與系泊系統(tǒng)動力學(xué)控制方程通過導(dǎo)纜點動力響應(yīng)聯(lián)立同時求解,控制方程組可以表示為

式中,{Fm} 為系泊系統(tǒng)對船體提供的作用力向量,其中力和力矩分別由各個系泊線在導(dǎo)纜點處的系泊張力合成以及對船體參考點的力矩合成;[T]為船體參考點位移與導(dǎo)纜點位移的轉(zhuǎn)換矩陣。

圖4 耦合動力分析方法Fig.4 Coupled dynamic analysis method

2 模型參數(shù)

2.1 船體模型

本文算例FPSO 是一艘大排水量船,設(shè)計排水量高達(dá)48 萬噸,船體結(jié)構(gòu)參數(shù)如表1 所示。船體設(shè)計吃水下濕表面網(wǎng)格如圖5所示,船體坐標(biāo)系原點位于船體重心在水線面的投影點處。

表1 FPSO參數(shù)Tab.1 Parameters of FPSO

圖5 船體濕表面網(wǎng)格模型Fig.5 Wet surface mesh model of the hull

圖6 多點系泊系統(tǒng)布置方式Fig.6 Layout of multi-point mooring system

2.2 系泊系統(tǒng)

通用型FPSO可搭載4×4型式的多點系泊系統(tǒng),如圖6所示,每一組系泊線與船中線的夾角為45°,每一組的線間角為5°,系泊半徑為3200 m,錨點水深為2000 m,分別編號為#1~#4、#5~#8、#9~#12 和#13~#16。每根系泊線由頂端錨鏈-中間鋼繩-底端錨鏈的形式組成,具體參數(shù)如表2所示。

表2 系泊線結(jié)構(gòu)參數(shù)Tab.2 Structural parameters of the mooring lines

2.3 環(huán)境載荷

船體波浪載荷基于勢流理論采用面元法求解波浪速度勢,再通過壓力積分獲得頻域波浪載荷,其中平均漂移力采用Newman近似法進(jìn)行計算。

風(fēng)載荷計算公式為

式中,CWX、CWY、CWN分別為Isherwood 方法[16]的縱向、橫向和首搖載荷系數(shù),ρa(bǔ)為空氣密度,Uwr為船體相對風(fēng)速的運(yùn)動速度,AT、AL分別為船體中橫剖面和中縱剖面面積,L為船長。圖7為FPSO 在不同方向θ(迎浪0°)下的風(fēng)載荷系數(shù)。

流載荷計算公式為

式中,CCX、CCY、CCN分別為OCIMF 方法[17]的縱向、橫向和首搖載荷系數(shù),ρw為海水密度,Ucr為浮船體相對流速的運(yùn)動速度,T為吃水,L為船長。圖8為FPSO在不同方向θ(迎浪0°)下的流載荷系數(shù)。

圖7 FPSO風(fēng)載荷系數(shù)Fig.7 Wind load coefficients of the FPSO

圖8 FPSO流載荷系數(shù)Fig.8 Current load coefficients of the FPSO

對于系泊線等細(xì)長結(jié)構(gòu)物,通常采用Morison 載荷模擬洋流對其作用,單位長度圓柱體運(yùn)動過程中受到的Morison載荷為

式中,D為圓柱直徑;vc、un分別為流體質(zhì)點和結(jié)構(gòu)節(jié)點的法向速度;CD為拖曳力系數(shù),文中無擋錨鏈取2.4,中間鋼絲繩取1.2;Ca為附加質(zhì)量系數(shù),文中均取1.0。

3 分析工況

表3 給出了4 組計算工況的環(huán)境參數(shù),其中波浪采用JONSWAP 譜進(jìn)行模擬,風(fēng)采用NPD 風(fēng)譜進(jìn)行模擬,流采用定常流,三者方向相同(船首迎浪),如圖9所示。

表3 計算工況Tab.3 Cases for calculation

圖9 風(fēng)浪流方向Fig.9 Incoming directions of wind,wave and current

4 結(jié)果與討論

本文采用WAMIT 進(jìn)行船體頻域水動力計算;采用SIMA 軟件進(jìn)行船體-系泊系統(tǒng)動力響應(yīng)計算,其中Simo 模塊用于計算船體運(yùn)動響應(yīng),能夠?qū)⑾挡聪到y(tǒng)采用懸鏈線方程進(jìn)行模擬,Riflex 模塊采用有限元方法模擬系泊系統(tǒng),用于計算系泊線的動力響應(yīng)。Simo 和Riflex 兩個模塊可以實現(xiàn)解耦和耦合兩種分析方法。

4.1 系泊剛度

分別采用準(zhǔn)靜態(tài)分析方法的懸鏈線理論和動力分析方法的有限元模型,對系泊系統(tǒng)在船體縱蕩方向上不同偏移位置處的靜態(tài)水平系泊張力進(jìn)行計算,結(jié)果如圖10 所示。可以看出兩種方法的計算結(jié)果基本相同,且在偏移平衡位置100 m范圍內(nèi)水平系泊張力與偏移量呈近似線性關(guān)系,超出100 m 后表現(xiàn)出微弱的非線性。因此,在線性范圍內(nèi),圖10 中曲線的斜率可近似認(rèn)為是定值,則斜率即為多點系泊系統(tǒng)在縱蕩方向上的水平系泊剛度Km,11=40 099 N/m。考慮船體縱蕩無窮附加質(zhì)量A∞=3599.1 t 時,船體在初始位置的縱蕩固有周期T11為

圖10 系泊張力水平分量-縱蕩偏移關(guān)系Fig.10 Horizontal mooring tension-surge offset relationship

4.2 系泊阻尼

首先給船體設(shè)置初始縱蕩偏移66 m,分別采用耦合動力(coupled dynamic, CD)分析方法和準(zhǔn)靜態(tài)(quasi-static, QS)方法進(jìn)行靜水衰減數(shù)值分析,結(jié)果如圖11 中黑色實線和藍(lán)色點劃線所示。可以看出兩種方法的結(jié)果差別較大,主要是因為準(zhǔn)靜態(tài)方法中未計入系泊阻尼。

為了獲得系泊系統(tǒng)對船體運(yùn)動的阻尼,根據(jù)耦合動力分析結(jié)果利用對數(shù)衰減法可以得到系泊阻尼系數(shù)ζ為0.2,進(jìn)而運(yùn)用公式(2)~(3)獲得系泊阻尼,并將其應(yīng)用于準(zhǔn)靜態(tài)分析方法中進(jìn)行相同的靜水衰減數(shù)值分析,結(jié)果如圖11 中紅色虛線所示。從圖11 中可以看出,在準(zhǔn)靜態(tài)分析方法中考慮系泊阻尼后,結(jié)果與耦合動力分析方法的差別較小,同時也說明系泊阻尼對船體運(yùn)動響應(yīng)的影響十分重要,在數(shù)值計算中應(yīng)予以考慮。

圖11 系泊船體縱蕩靜水衰減曲線對比Fig.11 Comparison of static decay curves of moored ship in surge direction

4.3 船體位移

針對表3 中的4 組分析工況,采用耦合動力分析方法(CD)與準(zhǔn)靜態(tài)方法(QS)數(shù)值計算船體在縱蕩方向上的運(yùn)動響應(yīng),并統(tǒng)計兩種方法計算結(jié)果的平均值(mean)及標(biāo)準(zhǔn)差(standard deviation),如圖12所示。從圖中可以看出,兩種方法得到的船體縱蕩位移響應(yīng)時歷整體上基本相同,因此平均位移基本相同,而準(zhǔn)靜態(tài)方法(QS)的標(biāo)準(zhǔn)差略大于耦合動力分析方法(CD),可能的原因是準(zhǔn)靜態(tài)方法(QS)中的線性系泊阻尼略小于耦合動力分析方法(CD)中的系泊系統(tǒng)動態(tài)阻尼效應(yīng)。

圖12 船體縱蕩位移響應(yīng)及統(tǒng)計結(jié)果對比(QS:準(zhǔn)靜態(tài)方法;CD:耦合動力方法)Fig.12 Comparison of surge displacement responses and statistic results(QS:Quasi-static method;CD:Coupled dynamic method)

4.4 系泊張力

由于準(zhǔn)靜態(tài)方法的船體運(yùn)動響應(yīng)與耦合動力分析方法的結(jié)果基本相同,因此可以將準(zhǔn)靜態(tài)方法的船體運(yùn)動時歷作為系泊系統(tǒng)運(yùn)動邊界條件,即采用解耦動力(uncoupled dynamic,UD)分析方法求解系泊系統(tǒng)動力響應(yīng)。圖13 給出了不同工況下#4 和#8 系泊線分別采用三種分析方法的系泊張力時歷對比結(jié)果,可以看出準(zhǔn)靜態(tài)分析方法得到的系泊張力與動力分析方法相比瞬時脈動值明顯偏小,這主要是因為在準(zhǔn)靜態(tài)方法中未考慮系泊線受到的莫里森力,而兩種動力分析方法獲得的系泊張力變化基本相同,僅在瞬時脈動值上略微存在差異。

圖13 系泊張力時歷對比Fig.13 Comparison of mooring tension time histories

為了進(jìn)一步比較三種數(shù)值計算方法的系泊張力差異,對系泊張力特征峰值進(jìn)行Weibull極值分布擬合,將概率水平為99%的特征極值進(jìn)行對比,以工況4(Hs=8 m)為例,結(jié)果如圖14所示。可以看出:對于張緊的系泊線#1~#4 和#13~#16,兩種動力分析方法獲得的系泊張力極值明顯大于準(zhǔn)靜態(tài)分析方法;對于松弛的系泊線#5~#8,三種方法結(jié)果比較接近;此外,兩種動力分析方法的結(jié)果基本相同,系泊線#1的差異最大,但僅占耦合動力分析方法結(jié)果的4.96%。因此,采用解耦動力分析方法整體上可以獲得與耦合動力分析方法同樣精度的結(jié)果。

圖14 系泊張力Weibull極值對比(工況Hs=8 m)Fig.14 Comparison of mooring tension Weibull extreme values(Hs=8 m)

4.5 計算耗時

為了比較三種數(shù)值方法的計算效率,進(jìn)一步統(tǒng)計了三種數(shù)值方法計算每個工況所耗費的時間,如表4所示。圖15所示為三種方法的耗時比率,其中準(zhǔn)靜態(tài)方法(QS)耗時最少,為耦合動力分析方法(CD)耗時的3.02%,解耦動力分析方法(UD)耗時介于準(zhǔn)靜態(tài)方法(QS)與耦合動力分析方法(CD)之間,約為耦合動力分析方法(CD)耗時的一半。因此,采用解耦動力分析方法能夠在保證計算精度的基礎(chǔ)上節(jié)省計算時間,提高計算效率。

圖15 三種數(shù)值分析方法計算耗時比率Fig.15 Time consuming ratios of the three numerical analysis methods

表4 三種方法計算耗時對比Tab.4 Comparison of consuming time in the three methods

5 結(jié) 論

本文分別采用準(zhǔn)靜態(tài)分析、解耦動力分析和耦合動力分析等數(shù)值方法開展了通用型FPSO多點系泊系統(tǒng)總體動力響應(yīng)計算分析,針對三種分析方法的船體運(yùn)動與系泊張力等結(jié)果的計算精度以及計算效率進(jìn)行了對比研究,主要得到以下結(jié)論:

(1)系泊系統(tǒng)對船體運(yùn)動除了提供回復(fù)力外,還具有一定的阻尼作用,對船體運(yùn)動響應(yīng)影響較大,因此在準(zhǔn)靜態(tài)分析方法中需要計及系泊阻尼,可采用靜水衰減耦合動力數(shù)值方法求解系泊系統(tǒng)的準(zhǔn)靜態(tài)線性阻尼系數(shù);

(2)采用計及系泊阻尼的準(zhǔn)靜態(tài)分析方法獲得的船體運(yùn)動響應(yīng)與耦合動力分析方法基本一致,但系泊張力由于未考慮莫里森力等作用而相差較大,在解耦動力分析方法中建立系泊系統(tǒng)有限元模型,考慮了莫里森力等多種非線性因素,得到的系泊張力與耦合動力分析方法基本相同;

(3)耦合動力分析方法耦合程度最高,但是計算耗時較長,準(zhǔn)靜態(tài)分析方法計算最快,但是無法計及流載荷導(dǎo)致結(jié)果存在偏差,解耦動力分析方法能夠在保證系泊張力計算準(zhǔn)確性的基礎(chǔ)上節(jié)省計算時間,從而提高計算效率;

(4)建議在系泊系統(tǒng)初步設(shè)計階段采用準(zhǔn)靜態(tài)方法確定初步設(shè)計方案,在詳細(xì)設(shè)計階段可采用解耦動力分析方法進(jìn)行安全性校核。

猜你喜歡
系統(tǒng)
Smartflower POP 一體式光伏系統(tǒng)
WJ-700無人機(jī)系統(tǒng)
ZC系列無人機(jī)遙感系統(tǒng)
北京測繪(2020年12期)2020-12-29 01:33:58
基于PowerPC+FPGA顯示系統(tǒng)
基于UG的發(fā)射箱自動化虛擬裝配系統(tǒng)開發(fā)
半沸制皂系統(tǒng)(下)
FAO系統(tǒng)特有功能分析及互聯(lián)互通探討
連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini
一德系統(tǒng) 德行天下
PLC在多段調(diào)速系統(tǒng)中的應(yīng)用
主站蜘蛛池模板: 丁香六月激情综合| 欧美日韩资源| 香蕉视频在线观看www| 欧美a级完整在线观看| 久久国产精品影院| 国产日韩丝袜一二三区| 在线亚洲天堂| 孕妇高潮太爽了在线观看免费| 久久久久久久97| 97人妻精品专区久久久久| 香蕉久久国产超碰青草| 综合社区亚洲熟妇p| 亚洲一区网站| 中国国产A一级毛片| 久久国产精品娇妻素人| 国产亚洲视频在线观看| 她的性爱视频| 九九热免费在线视频| 国产精品网址在线观看你懂的| 精品91视频| 久久久精品久久久久三级| 在线a网站| 天天综合天天综合| 欧美亚洲网| 国产极品美女在线观看| 国产自产视频一区二区三区| 真人免费一级毛片一区二区| 夜夜爽免费视频| 免费看黄片一区二区三区| 1024国产在线| 全部免费毛片免费播放| 激情综合网激情综合| av一区二区三区在线观看| 波多野结衣视频网站| 9久久伊人精品综合| 不卡的在线视频免费观看| av手机版在线播放| 91在线激情在线观看| 91精品国产一区| 欧美成人综合视频| 久久综合亚洲色一区二区三区 | 日韩高清一区 | 日本人又色又爽的视频| 精品无码日韩国产不卡av| 久草青青在线视频| 亚洲人成色在线观看| 亚洲熟女中文字幕男人总站| 色悠久久久久久久综合网伊人| 亚洲国产成人综合精品2020| 婷婷午夜天| 国产精欧美一区二区三区| 国产9191精品免费观看| 久久伊人操| 国产丝袜丝视频在线观看| 久久人体视频| 毛片久久网站小视频| 欧美天堂在线| 亚洲一本大道在线| 国产伦精品一区二区三区视频优播| 强乱中文字幕在线播放不卡| 思思99热精品在线| 激情六月丁香婷婷| 中文字幕在线日韩91| 婷婷99视频精品全部在线观看| 极品国产在线| 91人人妻人人做人人爽男同 | 美女被操黄色视频网站| 国内精品自在欧美一区| 亚洲一区二区日韩欧美gif| 国产欧美日本在线观看| 美女毛片在线| 人妖无码第一页| 国产精品人成在线播放| 99久久精品国产自免费| 香蕉视频在线观看www| 亚洲男女在线| 亚洲日本精品一区二区| 久久久久中文字幕精品视频| 日韩精品亚洲人旧成在线| 中文字幕有乳无码| 99中文字幕亚洲一区二区| 99一级毛片|