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

新型多立柱式Spar平臺渦激運動特性研究

2015-11-22 05:30:28蘇云龍楊建民呂海寧
海洋工程 2015年1期

蘇云龍,楊建民,呂海寧,彭 濤

(上海交通大學 海洋工程國家重點實驗室,上海 200240)

自20 世紀90年代以來,Spar 平臺憑借其優良的運動性能和穩性在深水油氣開采中得到了廣泛的應用。然而伴隨著Spar 平臺的發展,渦激運動(VIM)問題長期困擾著工程界。當平臺主體遭遇一定流速的海流作用時,由于粘性和逆壓作用,將在平臺附近產生流動分離,并在平臺后方形成周期性的旋渦脫落,當旋渦脫落頻率與平臺的固有頻率接近時將產生較大幅度的渦激運動。渦激運動會導致系泊系統及立管系統所受張力增加,使得系泊系統和立管系統疲勞壽命降低。

目前,針對常規的單立柱式Spar 平臺渦激運動的研究主要以試驗研究為主,Roddier 等[1]針對Truss Spar 平臺的硬艙分別采用三種縮尺比模型進行試驗,研究發現:對于模型有側板的情況,雷諾數范圍處于亞臨界區或是超臨界區對試驗結果影響不大。Finnigan 等[2]研究發現:剪切流和均勻流中的平臺的渦激運動特性沒有顯著區別;波與流同向時,平臺的渦激運動會被抑制;而波與流垂直或有一定夾角時,平臺渦激運動響應幅值小于平臺在單獨流與單獨波作用下的響應幅值之和,但某些情況會大于單獨流作用時的響應幅值。

除了模型試驗外,近年來也有不少學者采用計算流體動力學(CFD)的方法對單立柱式Spar 平臺的渦激運動進行研究。Thiagarajan 等[3]模擬了三維圓柱無側板和有側板情況下渦激運動響應的對比;Oakley 和Constantinides[4]模擬了Truss Spar 硬艙部分僅有側板和有側板及錨鏈、管系等附屬物的情況下渦激運動響應的對比,并與相應的模型試驗結果進行了比較;WANG 等[5]對一Cell-Truss Spar 平臺進行了數值模擬,并與相應的試驗結果進行了對比。此外,高云等[6]通過數值求解平臺和尾流陣子之間的耦合方程對無側板常規單立柱Spar 的渦激運動響應進行了研究。

針對單圓柱繞流,國內外學者已經開展了大量的工作。Achenbach[7]研究了單圓柱繞流表面壓力和摩擦力在不同Re 數時的分布情況。Schewe[8]通過風洞試驗研究了圓柱在雷諾數從亞臨界區域到超臨界區域范圍內的受力特性。單圓柱繞流的機理可用于解釋立管等的渦激振動和常規的單立柱式Spar 平臺的渦激運動等問題。

針對單方柱繞流和多柱繞流目前也已開展了大量的研究。Dutta 等[9]通過流場可視化技術,研究了方柱在0°、22.5°、30°和45°來流角下的阻力系數、斯特哈爾數及流場結構等。Sumner 等[10]通過流體染色與PIV結合的方法對兩個圓柱不同間距比、不同來流角下的流場進行了分類識別。徐有恒等[11]通過風洞試驗研究了正品字和倒品字排列三圓柱脈動壓力幅值的大小及RMS 值之間的定量關系。Lam 等[12]通過風洞試驗研究了不同間距下四圓柱方形陣列布置在0°、15°、30°和45°來流角下各柱的受力系數、斯特哈爾數等。吳七二等[13]通過CFD 方法對串列兩方柱繞流進行了研究,得到了不同間距比下柱體的受力情況。

此外,針對多柱浮式平臺的渦激運動也已經開展了不少研究。Magee 等[14]通過拖曳試驗模擬了張力腿平臺的渦激運動特性,并與以往的結果進行了對比。Goncalves 等[15]通過模型試驗對一半潛式平臺,分析了不同來流角、流速下,平臺的受力幅值、運動幅值、響應頻率等特性。國內,白治寧等[16]對四立柱深吃水半潛式平臺渦激運動進行了試驗和CFD 計算研究。

文中新型多立柱式Spar 平臺與常規的單立柱式Spar 平臺在結構形式上有著很大的不同,為了便于施工建造,平臺硬艙采用四根圓柱方形陣列布置的形式,另外采用方形中心井對立管起到保護和約束的作用,如圖1 所示。目前還沒有針對四圓柱與單方柱組合的渦激運動特性方面的研究,各柱體間復雜的相互干擾和結構自身關于不同來流角的非對稱特性,使得其繞流特性極為復雜,各柱的流動分離和渦結構將與單圓柱、單方柱或四圓柱存在很大的不同。因此針對這種特殊形式結構物的渦激運動特性開展模型試驗研究具有重要的意義。

1 新型多立柱式Spar 平臺模型及環境條件

1.1 平臺模型

新型多立柱式Spar 平臺的試驗模型如圖1 所示,縮尺比λ 為1∶ 60,據此確定模型的尺寸,保證與實體的幾何相似。模型包括硬艙、桁架及軟艙,由于上層建筑對平臺所受水動力不會產生影響,試驗中將其簡化為等效質量。平臺渦激運動的抑制裝置由一組減渦側板組成,側板沿縱向螺旋布置,每根圓柱上布置四塊減渦側板,相互間隔90°,側板高為5%d(d 為圓柱直徑),螺距比為5。減渦側板被設計為可拆卸形式,以校驗其對渦激運動的抑制效果。另外,試驗中運動與動力相似采取的是Froude 相似和Strouhal 相似準則,在有減渦側板存在時,流動分離多發生在側板邊緣處,粘性的影響不再占主導。

具體相似準則形式:

幾何相似:

運動與動力相似:

式中:L,A,▽,U,T 分別為特征線尺度、特征面積、特征體積、特征速度和時間;λ 為縮尺比;g 為重力加速度;下標m 及s 分別表示模型和實體。

關于試驗平臺模型主要尺度,如表1 所示。

圖1 多立柱式Spar 平臺模型及設計圖Fig.1 Test model and design drawing of the multi-column spar platform

表1 多立柱式Spar 平臺主要參數Tab.1 Main parameters of the multi-column spar platform

平臺的渦激運動主要是由平臺的硬艙部分引起的,而硬艙截面形狀在縱向上基本保持一致,因此平臺的數值模型可簡化為平臺硬艙截面所對應的二維流場,其具體形式如圖2 所示(以0°來流角情況為例)。圖中,D0為平臺在0°來流角下的迎流投影寬度,即硬艙寬度。

流域左側設定為速度入口,右側設定為自由出流,上下邊界設為對稱邊界。計算區域上游來流區域為8D0,上下邊界距中線各8D0,下游尾流區域為28D0,以保證結構后方渦街的充分發展。研究中要捕捉到流動分離和尾流近場部分細致結構,因此在平臺近壁面附近通過Size Function 布置較密集的網格,平臺附近及尾部一定范圍內均采用較密集的網格。另外,平臺附近區域的網格形式采用的是非結構化網格,以方便通過該區域網格的變形和重劃來實現平臺運動的模擬。計算采用SST k-ε 湍流模型,壁面處網格滿足y+≈1,共劃分181 822 個單元。計算過程中,平臺附近約2D0×2D0范圍內的網格設置為不變形,隨平臺一起運動,以保證平臺壁面附近網格的質量。

圖2 數值模型計算域大小及網格劃分示意Fig.2 Computational domain and computational mesh of numerical model

1.2 系泊系統模型

試驗采用水平等效系泊系統,共含四根系泊纜,每根間隔90°,系泊纜由鋼絲繩和彈簧組成。在保證橫蕩和縱蕩運動固有周期相同的基礎上,將原系泊系統等效為水平系泊。通過靜水衰減試驗測得平臺的橫蕩和縱蕩的固有周期與理論值106.5 s 十分接近,具體結果如表2 所示。考慮到平臺本身結構上具有一定的對稱性,試驗模擬了0°、15°、30°和45°來流角。具體的系泊系統、坐標系以及來流角規定如圖3 所示。

表2 靜水衰減試驗結果Tab.2 Calm-water decay test results

圖3 系泊系統坐標系和模型照片示意Fig.3 Layout of the mooring system,coordinate system (0deg)and different current incidences

數值計算中通過在用戶自定義函數(UDF)使用四階龍格-庫塔方法來求解平臺的運動微分方程,進而得到平臺的位移和運動速度,實現數值模型的運動模擬,以橫蕩的運動微分方程求解為例,具體如下:

應用四階龍格-庫塔求解得:

其中,m 為平臺質量;C 為結構阻尼系數;Ky為系泊系統在y 方向上的剛度;FHy為水動力在垂直流向上的分量,即升力;h 為時間步長,n 為迭代步數。由于數值模擬過程中,將壓力和粘性的貢獻以顯示方式并入N-S方程中直接進行求解,因此方程中不包含附加質量和水動力阻尼成分。

圖4 特征尺度D 的規定Fig.4 Definition of the characteristic length D

1.3 環境條件

本研究的是流致渦脫引起的渦激運動,考慮的環境條件主要是流速和來流角。而折合流速Ur是研究渦激運動時一個重要的無量綱參數,其具體表達式如式(5)所示。模型渦激運動試驗中根據折合速度來確定相應的環境工況,各試驗工況見表3 所示。

式中:U 為流速;T 為平臺靜水中橫蕩的固有周期;D 為平臺的特征尺度,具體參見圖4 所示。

表3 試驗工況表Tab.3 Table of test cases

2 渦激運動響應特性

2.1 垂直流向運動響應結果及分析

渦激運動最主要的特征是較大幅度的垂直流向的運動,即通常所說的橫蕩運動(Sway)。研究橫蕩運動時,通常關注的是其響應幅值的統計結果,橫蕩幅值對系泊系統強度設計和疲勞分析有著重要的影響。圖5給出了0°、15°、30°和45°來流角下平臺無量綱橫蕩幅值試驗結果,其中無量綱橫蕩幅值表達式:

其中,y 為垂直于流向運動,即橫蕩運動;為了便于比較不同來流角的響應結果,特征尺度統一取為D0,即0°來流角對應的特征尺度(實型值:43 m,模型值:0.717 m)。

圖5 不同來流角下橫蕩運動無量綱幅值隨Ur變化曲線Fig.5 Nondimensional amplitudes for the sway motions for different current incidences

由圖5 可得,平臺在各來流角下橫蕩運動幅值變化規律基本一致:

1)在較低速度時,各來流角的橫蕩幅值均較小;隨流速增大,橫蕩幅值迅速增大,逐漸進入共振范圍;

2)達到峰值后,由于“鎖定”效應,在一定流速范圍內橫蕩幅值始終保持平穩的狀態;

3)之后0°和15°來流角的橫蕩運動幅值明顯減小,逐漸脫離“鎖定”,而30°和45°來流角的橫蕩運動幅值僅略有減小,沒有明顯的脫離“鎖定”的趨勢。

對比各來流角結果可得:在流速較低時,各來流角橫蕩幅值相差不是很明顯。隨著流速的增大,各來流角橫蕩幅值結果差異逐漸增大。有側板時45°來流角下橫蕩幅值最大,為0.428D0;無側板時0°來流角下橫蕩幅值最大,為0.562D0。

另外,通過對比有/無側板的結果發現:減渦側板在0°來流角流速較高時對橫蕩運動有很好的抑制效果,使得幅值減為無側板時的1/4 左右。

圖6 不同來流角下橫蕩運動幅值計算結果與試驗結果對比Fig.6 Comparison of nondimensional amplitudes for the sway motions for different current incidences

圖6 給出了數值模擬得到的平臺橫蕩運動的幅值結果,及其與試驗結果的對比。由圖可得計算結果與試驗結果變化趨勢基本一致,但比試驗結果偏于保守,這主要是由于數值模擬未考慮硬艙各柱的連接部分,連接部分的存在使得硬艙各柱體在縱向上的連續性一定程度上被打破,這會對平臺的渦激運動起到抑制作用。同時,由于數值模擬是對平臺渦激運動的初步研究,采用的是二維計算,忽略了平臺硬艙底部及連接部分上下面的繞流影響,以及平臺縱向上不同高度間的繞流不完全一致產生的相互干擾等三維效應。所以計算結果偏于保守,共振區域也比試驗結果的偏大。

2.2 XY 平面內運動軌跡

根據平臺重心處的橫蕩和縱蕩運動時歷,可得到平臺在XY 平面內的運動軌跡,如圖7 所示,為便于比較,統一采用D0對x 和y 時歷結果進行無量綱化處理,其中平臺截面僅作為參照,非原始尺寸比例。

由平臺的水平面運動軌跡結果可得:平臺水平面內運動軌跡的形狀在0°來流角無側板時呈“8”字形,其余情況呈“香蕉”形。在15°、30°和45°時,運動軌跡近似與平臺硬艙截面的對角線平行,這一特點與常規的單立柱Spar 平臺存在很大的不同,而與四立柱半潛式平臺的結果[17]相類似。

圖7 不同來流角下平臺XY 平面內運動軌跡圖Fig.7 Motions in the XY plane for different current incidences

2.3 首搖運動響應結果及分析

由于文中所研究平臺結構形式上的特殊性,試驗中除了橫蕩和縱蕩外,還觀察到了較大的首搖運動,這不同于常規的單立柱Spar 平臺,而與四立柱半潛式平臺[15]類似。圖8 給出了不同來流角下平臺首搖運動幅值統計結果,其中首搖運動幅值表達式如式(7)所示,rot 為平臺首搖運動:

由圖8 可得,平臺首搖運動幅值隨Ur的變化趨勢與橫蕩幅值隨Ur的變化趨勢基本一致,即開始時幅值較小;隨著逐漸進入“鎖定”范圍,幅值迅速增大;之后0°和15°來流角有較明顯的脫離“鎖定”。

對比各來流角結果可得:有側板時各來流角的首搖幅值相差不是很明顯,約為3°;無側板時,各來流角首搖幅值存在一定差異,0°來流角下首搖幅值最大,約為6.6°。

另外,通過對比有/無側板的結果發現:減渦側板在0°來流角流速較高時對首搖運動有很好的抑制效果,這也與橫蕩運動結果相類似。

圖8 不同來流角下首搖運動幅值隨Ur變化曲線Fig.8 Amplitudes for the yaw motions for different current incidences

3 渦激運動水動力及流場特性

3.1 動力響應結果及分析

平臺在渦激運動中所受到的水動力結果對于研究平臺渦激運動特性同樣具有重要意義。試驗測量得到了各系泊纜的張力及平臺的運動時歷x(t)及y(t)。平臺系泊系統為張緊式系泊,可近似認為其為線性,相應的橫蕩和縱蕩運動方程:

其中,m 為平臺質量;C 為結構阻尼系數;Kx、Ky為系泊系統在x 和y 方向上的剛度;FHx和FHy為水動力在沿流向和垂直流向上的分量,即阻力和升力,包含了附加質量和水動力阻尼。基于運動時歷,可通過數值方法求得平臺的速度與加速度,將x 和y 方向上的系泊張力分量及加速度值代入運動方程,即可求得到平臺所受到的水動力。

根據水動力結果可求得相應的渦激運動水動力系數,其具體表達式:

其中,ρ 為流體密度;AP為浸沒部分投影面積;U 為流速。

水動力系數的統計結果對于研究渦激運動現象具有重要意義。通常,升力系數取統計均方根結果,記為C'L,阻力系數取統計平均結果,記為在圖9 中,分別給出了不同來流角,不同折合速度下的水動力系數統計結果,其中為了方便比較不同來流角的結果,橫坐標采用的是無量綱化的折合速度Ur來表征不同流速。

由圖9(a)和(b)可得:平臺C'L 的變化趨勢與橫蕩運動幅值變化趨勢基本一致。無減渦側板0°來流角時,達到最大,約為0.4;其余來流角最大均約為0.2 左右。且各來流角下均在折合速度約5 ~9 范圍內達到峰值范圍。

圖9 不同來流角下水動力系數統計結果隨Ur變化曲線Fig.9 Lift and drag coefficients for different incidence angles

圖10 不同來流角下水動力系數統計結果與試驗結果對比Fig.10 Comparison of lift and drag coefficients for different incidence angles

3.2 繞流流場特性

通過CFD 計算可以得到流場中各位置處的詳細信息,根據各柱渦結構間的干擾和平臺整體的旋渦脫落情況,可研究渦激運動中各柱間的擾動和平臺整體渦激運動的特性。這里研究的圓柱圓心距與直徑之比(l/d)約為2.31,按楊佼等[18]的計算結果,間距比為2 左右時,對于正型排布四圓柱(對應0°來流角情況),前面的圓柱渦結構會貼附到后柱上,而上下兩組圓柱各自產生不對稱泄渦。由于文中結構物還包含了中間的方柱,所以泄渦形式與四圓柱的情況有著較大的不同:各來流角下圓柱和方柱都會產生一定程度的流動分離和渦突起,各柱之間存在復雜的相互擾動,而在柱后形成共同的渦結構,類似單柱繞流的尾渦模式。圖11 給出了不同來流角下的瞬時渦量圖。以下分析中,各圓柱的編號說明參見圖4 所示。

由圖11 可得:0°來流角時,迎流面的A 柱和D 柱均有一定的渦突起產生,方柱流動分離發生在直角處,由于方柱脫渦的影響使得A 柱和D 柱渦結構不再關于來流方向對稱而是偏向上下兩側;A 柱和D 柱的渦結構與方柱的一起分別貼附到后柱B 柱和C 柱上,最終在B 柱和C 柱后產生結構整體的漩渦脫落。

15°來流角時,A 柱和D 柱后的渦結構會對方柱的產生一定的抑制,方柱上側渦結構從方柱與B 柱間經過;C 柱處于D 柱和方柱的屏蔽效應下,流動分離受到抑制。

圖11 不同來流角下的瞬時渦量圖Fig.11 Instantaneous vorticity contours for different current incidences

4 結 語

通過CFD 數值模擬和模型拖曳試驗的方法對一新型的多立柱式Spar 平臺的渦激運動特性進行了研究。所得主要結論總結如下:

1)平臺的橫蕩運動和首搖運動在各來流角下均存在較明顯的進入“鎖定”共振區間和“鎖定”共振區間,且在0°和15°來流角時還存在較明顯的脫離“鎖定”共振區間。

2)平臺的運動軌跡呈“8”字形或“香蕉”形,15°、30°和45°來流角下軌跡線近似與平臺硬艙截面的對角線平行。

3)平臺除了較大的橫蕩和縱蕩運動外,還存在較大的首搖運動,這與常規的單立柱式Spar 平臺有著很大的不同。

4)減渦側板在0°來流角且流速較高時對抑制渦激運動效果十分明顯。

6)平臺運動時的繞流在各柱均有一定的渦突起,并存在復雜的相互干擾;而整體上看,其尾渦為共同的渦結構,模式類似單柱繞流。

[1]RODDIER D,FINNIGAN T,LIAPIS S.Influence of the Reynolds number on Spar vortex induced motions (VIM):multiple scale model test comparisons [C]//Proceedings of the ASME 28th International Conference on Ocean,Offshore and Arctic Engineering.2009:OMAE2009-79991.

[2]FINNIGAN T,IRANI M,VAN DIJK R.Truss Spar VIM in waves and currents[C]//Proceedings of the 24th International Conference on Offshore Mechanics and Arctic Engineering.2005:OMAE2005-67054.

[3]THIAGARAJAN K P,CONSTANTINIDES Y,FINN L.CFD analysis of vortex-induced motions of bare and straked cylinders in currents [C]//Proceedings of the 24th International Conference on Offshore Mechanics and Arctic Engineering.2005:OMAE2005-67263.

[4]OAKLEY Jr O H,CONSTANTINIDES Y.CFD Truss Spar Hull benchmarking study[C]//Proceedings of the 26th International Conference on Offshore Mechanics and Arctic Engineering.2007:OMAE2007-29150.

[5]WANG Ying,YANG Jianmin,LU Haining.Computational fluid dynamics and experimental study of lock-in phenomenon in vortex-induced motions of a Cell-Truss Spar[J].Journal of Shanghai Jiao Tong University,Science,2009,14(6):757-762.

[6]高云,宗智,于馨.Spar 平臺渦激運動響應分析[J].中國海洋平臺,2011,26:17-22.(GAO Yun,ZONG Zhi,YU Xin.Response analysis of vortex induced spar motions[J].China Offshore Platform,2011,26:17-22.(in Chinese))

[7]ACHENBACH E.Distribution of local pressure and skin friction around a circular cylinder in cross-flow up to Re=5 ×106[J].Journal of Fluid Mechanics,1968,34(4):625-639.

[8]SCHEWE G.On the force fluctuations acting on a circular cylinder in crossflow from subcritical up to transcritical Reynolds numbers[J].Journal of Fluid Mechanics,1983,133:265-285.

[9]DUTTA S,PANIGRAHI P,MURALIDHAR K.Experimental investigation of flow past a square cylinder at an angle of incidence[J].Journal of Engineering Mechanics,2008,134(9):788-803.

[10]SUMMER D,PRICE S J,PAIDOUSSIS M P.Flow-pattern identification for two staggered circular cylinder in cross-flow[J].Journal of Fluid Mechanics,2000,411:263-303.

[11]徐有恒,張德龍,馬健.三圓柱繞流壓力的脈動幅度及其與RMS 值數值關系的實驗研究[J].水動力學研究與進展,1997,12(3):419-425.(XU Youheng,ZHANG Delong,MA Jian.An experimental study on the amplitude of fluctuating pressure of flow around a group of three cylinders[J].Journal of Hydrodynamics,1997,12(3):419-425.(in Chinese))

[12]LAM K,LI J Y,SO R M C .Force coefficients and Strouhal numbers of four cylinders in cross flow[J].Journal of Fluids and Structures,2003,18:305-324.

[13]吳七二,周岱,付功義.串列方形鈍體構筑物繞流的數值分析[J].上海交通大學學報,2012,46(1):130-135,141.(WU Qier,ZHOU Dai,FU Gongyi.Numerical analysis of flow past two bluff bodies with square-section in tandem arrangement[J].Journal of Shanghai Jiao Tong University,2012,46(1):130-135,141.(in Chinese))

[14]MAGEE A,SHEIKH R,GUAN K Y H,et al.Model tests for vim of multi-column floating platforms[C]//Proceedings of the 30th International Conference on Offshore Mechanics and Arctic Engineering.Rotterdam:[s.n.],2011:OMAE2011-49151.

[15]GONCALVES R T,ROSETTI G F,FUJARRA A L C ,et al.Experimental study on vortex-induced motions of a semi-submersible platform with four square columns,Part I:Effects of current incidence angle and hull appendages[J].Ocean Engineering,2012,54:150-169.

[16]白治寧.深吃水半潛式平臺渦激運動響應特性研究[D].上海:上海交通大學,2013.(BAI Zhining.Research on vortex induced motions of a deep draft semi-submersible platform[D].Shanghai:Shanghai Jiao Tong University,2013.(in Chinese))

[17]RIJKEN O,LEVERETTE S.Experimental study into vortex induced motion response of semi submersibles with square columns[C]//Proceedings of the 27th International Conference on Offshore Mechanics and Arctic Engineering.Estoril:[s.n.],2008:OMAE2008-57396.

[18]楊佼,武文斐,龔志軍,等.正型排布四圓柱繞流流動模式的格子Boltzmann 模擬[J].科學技術與工程,2013,13(19):5427-5430,5435.(YANG Jiao,WU Wenfei,GONG Zhijun,et al.Lattice Boltzmann simulation for flow pattern of flow around four cylinders in square arrangement[J].Science Technology and Engineering,2013,13(19):5427-5430,5435.(in Chinese))

主站蜘蛛池模板: 国产精品视频白浆免费视频| 国产精品一区不卡| 天堂中文在线资源| 国产在线八区| 国产Av无码精品色午夜| 国产麻豆另类AV| 美女潮喷出白浆在线观看视频| 免费在线国产一区二区三区精品| 人妻中文久热无码丝袜| 成年人视频一区二区| 国产精品白浆在线播放| 毛片免费在线视频| 欧美另类第一页| 91探花国产综合在线精品| 亚洲国产AV无码综合原创| 制服丝袜 91视频| 国产乱肥老妇精品视频| 国产在线观看高清不卡| 国产欧美日韩专区发布| 中文天堂在线视频| 国产又粗又猛又爽| 九九九精品视频| 99免费视频观看| 成人毛片免费在线观看| 幺女国产一级毛片| 亚洲美女高潮久久久久久久| 久久久久国色AV免费观看性色| 在线观看亚洲人成网站| 九色国产在线| 91久久国产综合精品| 自拍偷拍欧美| 久久精品中文无码资源站| 亚洲男人天堂2020| 国产一级裸网站| 国产成人无码AV在线播放动漫| 亚洲欧美精品一中文字幕| 久久人与动人物A级毛片| 国产免费黄| 青青草欧美| 亚洲色图在线观看| 国产亚洲精品yxsp| 亚洲av无码人妻| 欧洲亚洲一区| 在线色国产| 人与鲁专区| 午夜久久影院| 小说区 亚洲 自拍 另类| 国产精品黄色片| 天堂岛国av无码免费无禁网站 | 国产熟女一级毛片| 亚洲第一天堂无码专区| 久久这里只有精品免费| 日本免费新一区视频| 亚洲天堂久久| 国产毛片高清一级国语| 亚洲精品自拍区在线观看| 一级毛片基地| 国产精品综合久久久| 亚洲区一区| 国产男女XX00免费观看| 日本亚洲成高清一区二区三区| 国产AV无码专区亚洲精品网站| 亚洲AV人人澡人人双人| 无码AV高清毛片中国一级毛片| 欧美日韩一区二区在线播放| 97国产在线观看| 多人乱p欧美在线观看| 小蝌蚪亚洲精品国产| 亚洲日本韩在线观看| 99在线视频网站| 扒开粉嫩的小缝隙喷白浆视频| 91无码网站| 免费a级毛片视频| 日本免费a视频| 亚洲精品欧美重口| 99热国产在线精品99| 国产精品久久国产精麻豆99网站| 黄色成年视频| 亚洲精品在线观看91| 精品视频一区二区观看| 国产第一页亚洲| 免费人欧美成又黄又爽的视频|