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

張力腿平臺內孤立波作用特性數值模擬

2015-10-27 10:37:48林忠義尤云祥
海洋工程 2015年5期
關鍵詞:理論

王 旭,林忠義,尤云祥

(1. 上海交通大學 海洋工程國家重點實驗室,上海 200240;2. 嘉興南洋職業技術學院,浙江 嘉興 314003)

張力腿平臺內孤立波作用特性數值模擬

王 旭1,林忠義2,尤云祥1

(1. 上海交通大學 海洋工程國家重點實驗室,上海 200240;2. 嘉興南洋職業技術學院,浙江 嘉興 314003)

依據三類內孤立波理論KdV、eKdV和MCC的適用性條件,采用Navier- Stokes方程為流場控制方程,以內孤立波誘導上下層深度平均水平速度作為入口邊界條件,建立了兩層流體中內孤立波對張力腿平臺強非線性作用的數值模擬方法。結果表明,數值模擬所得內孤立波波形及其振幅與相應理論和實驗結果一致,并且在內孤立波作用下張力腿平臺水平力、垂向力及力矩數值模擬結果與實驗結果吻合。研究同時表明,張力腿平臺內孤立波載荷由波浪壓差力、粘性壓差力和摩擦力構成,其中摩擦力很小,可以忽略;水平力的主要成分為波浪壓差力和粘性壓差力,粘性壓差力與波浪壓差力相比較小卻不可忽略,流體粘性的影響較小;垂向力中粘性壓差力很小,流體粘性影響可以忽略。

兩層流體;內孤立波;張力腿平臺;載荷特性

內孤立波是最大振幅發生在密度分層海水內部的一種常見海洋動力學現象,不僅是海洋能量級串中的一個重要環節,也是對海洋工程結構物安全性產生重要影響的海洋環境因素。我國南海內孤立波活動頻繁。1990年,在流花油田就曾發生過因內孤立波導致纜繩拉斷、船體碰撞以致拉斷和擠破漂浮軟管的事故[1]。同年,在南海陸豐油田也發生過因內孤立波導致半潛鉆井船與錨定油輪在連接輸油管道時發生困難等問題[2],因此,內孤立波已成為深海資源開發中需要考慮的海洋環境因素之一。

內孤立波在其傳播過程中可以保持波形和傳播速度不變,這是由于非線性和色散效應在一定尺度上的平衡所致,一般地可以用KdV(Korteweg- de Vries)、eKdV (extended KdV)和MCC(Miyata- Choi- Camassa)等理論模型來描述。在KdV理論中要求內孤立波是弱非線性、弱色散且兩者平衡,而在eKdV理論中只要求內孤立波是弱非線性和弱色散的[3]。為克服需要弱非線性限制條件的缺陷,Choi和Camassa[4]建立了強非線性和弱色散的內孤立波理論,稱為MCC理論[3]。但在這三類理論中弱非線性和弱色散這兩個條件僅僅為定性描述,為此黃文昊等以系列實驗為依據給出了這兩個條件的定量表征方法[5]。

作為一種半剛性半順應性平臺,張力腿平臺受到的浮力遠大于自身重力,巨大的張力腿預張力限制了平臺橫搖、縱搖和垂蕩運動,因而該平臺在惡劣海況下仍具有良好的穩定性,是我國南海深海資源開發中首選的海洋工程裝備之一[6],合理地確定各種海洋環境條件下張力腿平臺的載荷特性則是深海張力腿平臺設計和應用中的一項關鍵問題。程友良等[7- 10]將Morison公式與KdV理論結合,研究了內孤立波作用下小尺度桿件的載荷特性問題;尤云祥等[11- 12]將Morison公式與eKdV理論結合,研究了內孤立波作用下張力腿和半潛平臺的載荷與動力響應問題;宋志軍等[13]將Morison公式與KdV理論結合,研究了內孤立波作用下Spar平臺的載荷與動力響應問題。需要指出的是,在這些文獻中,關于Morison公式中慣性力和拖曳力系數都是參照表面波的方法選取的,但這種選取方法缺乏理論和實驗依據。為此,黃文昊等[14]以系列實驗為依據針對張力腿平臺給出了這兩個系數的選取方法。

目前,對內孤立波作用下深海平臺載荷問題,許多機理尚不清楚,包括各種內孤立波載荷成分的形成機理,流體粘性對內孤立波載荷的影響機理,以及利用Morison公式計算深海張力腿平臺內孤立波載荷的合理性等。計算流體力學(Computational Fluid Dynamics,簡稱CFD)方法為進一步深入認識這些問題提供了一條有效的途徑。采用CFD方法可以直接獲得內孤立波與深海張力腿平臺強非線性相互作用過程中的流場變化特性,因此可以直接獲得內孤立波作用下深海張力腿平臺載荷等水動力特性。

關輝等[15]基于KdV理論,高原雪等[16]基于MCC理論,采用CFD方法研究了內孤立波的生成傳播問題。陳杰等[17]、劉碧濤等[18]基于eKdV理論,采用CFD方法分別研究了內孤立波作用下水下潛體和海洋立管的載荷特性等問題。需要指出的是,在這些文獻中,由于沒有考慮KdV、eKdV和MCC理論的適用性條件,致使其數值模擬結果均不同程度地出現了內孤立波振幅及其波形不可控等問題。因此,選擇合適的內孤立波理論作為CFD數值模擬的依據,是研究其生成傳播及其對海上結構物強非線性作用特性問題時需要解決的關鍵問題。

有鑒于此,本文采用 Navier- Stokes方程,依據文獻[5]確定的三類內孤立波理論KdV、eKdV和MCC的適用性條件,建立振幅及其波形可控的內孤立波CFD數值模擬方法。在此基礎上,對內孤立波與張力腿平臺的相互作用特性進行數值模擬,進而分析內孤立波作用下張力腿平臺各種載荷成分的形成機理及其影響程度,張力腿平臺對內孤立波的波形及其流場的影響,以及使用Morison公式等簡化方法計算張力腿平臺內孤立波載荷的合理性等問題。

1 數值方法

建立直角坐標系oxyz,其中oxy平面位于流體靜止時兩層流體的界面上,oz軸與張力腿平臺垂向中心軸重合且垂直向上為正。內孤立波為平面前進波,界面位移為ζ,沿ox軸正方向傳播。設各層均為不可壓流體,上層流體深度與密度分別為h1和ρ1,下層流體深度與密度分別為h2和ρ2。流場計算的控制區域如圖1所示,包括內孤立波生成傳播區和消波區兩個區域,內孤立波生成傳播區的長度為18m(圖中非陰影部分區域),消波區的長度為12 m(圖中陰影部分區域)。采用速度入口方法生成內孤立波,當造波區中形成穩定的內孤立波后,對所生成內孤立波的傳播特性進行監測分析,并對張力腿平臺的內孤立波載荷進行計算。

圖1 內孤立波數值水槽示意Fig. 1 Sketch of the numerical flume for internal solitary waves

采用Navier- Stokes方程作為控制方程進行數值模擬,流場控制方程為

其中:(u1,u2,u3)為速度矢量,p為壓力,t為時間,(f1,f2,f3)=(0,0,-g)為重力矢量,g為重力加速度,ν為運動粘性系數,ρ為流體密度,當ζ

張力腿平臺壁面取為無滑移邊界條件,而計算域頂部及底部滿足如下壁面條件:

設內孤立波振幅為a,相速度為c,則其誘導上下層流體中的層深度平均水平速度分別為[4]:

張力腿平臺主要由甲板結構、立柱和浮箱等柱型結構組成,其內孤立波水平及垂向載荷均由表面摩擦力和壓差力兩個部分組成,如下式所示:

其中:(nx,ny,nz)為張力腿平臺濕表面法向矢量,方向指向物體內部,S為張力腿平臺濕表面積,第一項為作用于張力腿平臺的摩擦力,第二項為作用于張力腿平臺的壓差力。

如何選擇合適的內孤立波理論來計算入口速度是數值模擬中一個需要解決的關鍵問題。根據文獻[5]定量確定的三類內孤立波的適用性條件,入口速度的具體計算方法如下:

采用有限體積法離散動量和連續性方程,對流項采用QUICK(quadratic upstream interpolation for convective kinetics)離散格式,壓力插值格式采用體積力加權(body force weighted)方法,壓力速度耦合迭代采用PISO(Pressure Implicit with Splitting of Operators)算法,兩層流體界面的構造方法選用幾何重構法。初始時間步長為Δt=0.005 s,計算過程中根據每個時間步長的收斂情況逐漸增加時間步長以縮短計算時間。

2 結果與分析

圖2 張力腿平臺模型(水下部分)示意Fig. 2 Sketch of the submerged part for the TLP model

文獻[14]利用大型密度分層水槽,對內孤立波作用下張力腿平臺模型的載荷特性進行了系列實驗,本文結合該文獻中的相關實驗結果進行數值模擬與分析。為此,數值水槽主尺度、上下層流體密度及其深度比均與該文一致。其中,數值水槽長度為30 m,水深為1 m,上下層流體密度分別為ρ1=998 kg/m3和ρ2=1 025 kg/m3,上下層流體深度比分別選擇h1∶h2=1∶9、2∶8和3∶7三種分層工況。

模擬計算中使用的張力腿平臺模型以ISSC- TLP平臺為原型,平臺模型水下部分如圖2所示,平臺由4根立柱和4個沉箱組成;立柱間距0.383 m,半徑0.037 5 m,高度0.338 m;沉箱高度0.047 m,寬度0.033 m;平臺重心高度0.169 m,吃水深度0.156 m,排水量4.7 kg。張力腿平臺重心距速度入口端9 m。計算區域采用六面體結構化網格進行離散,每根立柱圓周上布置80個計算單元,沉箱長度方向布置90個單元,總的單元數量為1 631 968個。

2.1內孤立波數值模擬結果

首先研究數值模擬過程中流體粘性對內孤立波生成與傳播特性的影響。為此,在數值模擬中設計如下兩種情況:一種為N- S模擬,依據Navier- Stokes方程求解,考慮流體粘性的情況;另一種為Euler模擬,依據Euler方程求解,不考慮流體粘性的情況。

圖3給出了當h1:h2=3:7和ad/h=0.106時(ad為內孤立波設計振幅),采用上述兩種方式對內孤立波生成與傳播特性的數值模擬結果。由圖可知,內孤立波在向右傳播過程中,有粘、無粘兩種情況下,數值模擬所得內孤立波在其傳播過程中均保持波形穩定、振幅衰減很小(兩者振幅相對誤差均在3%以內),沒有明顯的尾波現象。因此對內孤立波的CFD數值模擬中,粘性對內孤立波生成與傳播過程的影響較小,采用基于N- S和Euler方程的兩種方法均是可行的。下文中如無特別聲明,所有數值模擬均是在有粘情況下依據Navier- Stokes方程進行計算的。

圖3 當h1∶h2=3∶7和ad/h=0.121時兩種情況下內孤立波數值模擬結果Fig. 3 The numerical results for internal solitary waves in two cases when h1∶h2=3∶7 and ad/h=0.121

圖4給出了兩種不同工況下內孤立波波形的數值模擬結果,并與相應理論和實驗結果進行了比較,圖中實驗結果取自文獻[14]。其中,工況Case A條件下h1:h2=20:80,ad/h=0.053,此時內孤立波為弱非線性和弱色散的(此時ε=0.053,μ=0.061 1),選擇KdV理論計算入口速度;工況Case B條件下h1:h2=30:70,ad/h=0.106,此時內孤立波為中等非線性和弱色散的(此時ε=0.106,μ=0.025 1),選擇eKdV理論計算入口速度。

目前實驗室條件下常用的內孤立波造波方法主要有抽板式、捶擊式和推板式等,本文依據密度分層水槽,使用推板式造波方法進行內孤立波實驗室造波,實驗通過伺服控制系統驅動兩塊推板同時作反向水平運動,使水槽內的兩層流體產生反向流動,從而在密度分層流體界面上產生上凸或下凹的內孤立波[5]。

圖4 內孤立波波形數值模擬結果與理論和實驗結果比較Fig. 4 Comparisons of the numerical results for internal solitary wave waveforms with theoretical and experimental ones

結果表明,在兩種工況下,數值模擬所得內孤立波的波形,不僅與內孤立波理論解波形一致,而且與實驗所得波形吻合,這表明依據內孤立波理論的適用性條件,采用本文所述數值模擬方法所得內孤立波的波形是準確可控的。

圖5給出了當h1∶h2=20∶80、25∶75和30∶70時內孤立波數值模擬振幅am與其設計振幅ad之間相關關系。依據三類內孤立波理論的適用性條件,圖中除h1∶h2=20∶80,ad/h=0.053依據KdV方程,其余工況均依據eKdV方程進行數值造波。橫向坐標軸和縱向坐標軸分別為無因次設計振幅和無因次模擬振幅,圈號“О”表示數值模擬振幅,虛線表示設計振幅(其斜率為1),圈號與虛線之間的垂向距離表示兩者之間的絕對誤差。結果表明,在各數值模擬工況下,數值模擬所得內孤立波振幅均與其相應設計振幅符合較好,兩者之間的相對誤差不超過4%。這表明依據三類內孤立波理論的適用性條件,采用本文所述數值模擬方法所得內孤立波的振幅同樣是準確可控的。

圖5 內孤立波振幅數值模擬結果Fig. 5 The numerical results of the wave amplitudes for internal solitary waves

2.2內孤立波載荷特性

圖6 內孤立波無因次水平力、垂向力及力矩幅值數值與實驗結果Fig. 6 Results of numerical and experimental amplitudes for dimensionless horizontal and vertical forces, as well as torques due to internal solitary waves

由圖可知,數值模擬與實驗結果吻合,兩者之間的相對誤差一般不超過12%。個別工況誤差較大的主要原因在于:實際實驗過程中具有躍層的密度分層流體在多次實驗后,密度躍層厚度會逐漸增大,與模擬中兩層流體原有假設差異增大,從而使得內孤立波數值模擬與實驗振幅之間的相對誤差擴大。

圖7給出了在Case B工況下,內孤立波無因次水平力、垂向力及力矩時歷的數值模擬及實驗結果對比。由圖可知,張力腿平臺內孤立波載荷時歷數值結果與實驗時歷結果吻合,表明采用本文所述張力腿平臺內孤立波載荷的計算方法是合理可行的。

圖7 Case B工況下內孤立波無因次水平力、垂向力及力矩時歷特性Fig. 7 The time- variant characteristics for dimensionless horizontal and vertical forces, as well as torques due to internal solitary waves for Case B

圖8 Case B工況下內孤立波無因次壓差力及摩擦力時歷特性Fig. 8 The time- variant characteristics for the pressure difference and the fractional forces due to internal solitary waves for Case B

圖9 Case B工況下內孤立波無因次波浪壓差力及粘性壓差力時歷特性Fig. 9 The time- variant characteristics for the waves and viscous pressure forces due to internal solitary waves for Case B

3 結 語

依據KdV、eKdV和MCC理論解的適用性條件,采用N- S方程為控制方程,以內孤立波在上下層流體中誘導的深度平均水平速度作為入口條件,建立了內孤立波與張力腿平臺強非線性作用的數值模擬方法。該方法對張力腿平臺內孤立波水平力、垂向力和力矩幅值及其時歷變化特性的數值模擬結果與相應實驗結果一致,可以用于內孤立波與海洋浮式結構強非線性作用的數值模擬。

研究表明,張力腿平臺內孤立波水平和垂向力均由波浪壓差力、粘性壓差力和摩擦力組成。其中,水平摩擦力、垂向摩擦力均很小,可以忽略;水平力的主要成分為波浪壓差力和粘性壓差力,粘性壓差力與波浪壓差力相比較小卻不可忽略,流體粘性的影響較小;垂向力中粘性壓差力很小,流體粘性影響可以忽略。

[1] 陳景輝.南海流花11- 1深海油田開發工程[J].中國海洋平臺,1996,11(1):43- 45.(CHEN Jinghui. Liuhua 11- 1 deep sea oil- field development project in South China Sea[J]. China Offshore Platform, 1996,11(1):43- 45. (in Chinese) )

[2] BOLE JB, EBBESMEYER CC, ROMEA RD. Soliton currents in the South China Sea:measurements and theoretical modeling[C]//The 16thOffshore Technology Conference in Houston. Texas, USA, 1994, 367- 376.

[3] HELFRICH K R, MELVILLE W K. Long nonlinear internal waves[J]. Ann. Rev. Fluid Mech., 2006, 38:395- 425.

[4] CHOI W, CAMASSA R. Fully nonlinear internal waves in a two- fluid system [J]. J. Fluid Mech., 1999, 396: 1- 36.

[5] 黃文昊,尤云祥,王旭,等.有限深兩層流體中內孤立波造波實驗及其理論模型[J].物理學報, 2013,62(8):084705- 084714.(HUANG Wenhao, YOU Yunxiang, WANG Xu, et al. Wave- making experiments and theoretical models for internal solitary waves in a two- layer fluid of finite depth[J]. Acta Phys. Sin., 2013,62(8): 084705- 084719. (in Chinese))

[6] 張智,董艷秋,唐友剛,等.1990年后世界TLP平臺的發展狀況[J].中國海洋平臺,2004,19(2):5- 11.(ZHANG Zhi, DONG Yanqiu, TANG Yougang,et al.The development of TLP after 1990[J].China Offshore Platform,2004,19(2):5- 11. (in Chinese))

[7] CHENG YL, LI JC, LIU YF. The induced flow field by internal solitary waves and its action on cylindrical piples in the stratified ocean[J]. Recent Advances in Fluid Mechanics, Qinghua- Apringer,2004,296- 299.

[8] CAI SQ, LONG XM, GAN ZJ. A method to estimate the forces exerted by internal solitons on cylinder piles[J].Ocean Engineering, 2003, 30(5):673- 689.

[9] CAI SQ, LONG XM, WANG SG. Force and torques exerted by internal solitions in shear flows on cylinder piles[J]. Applied Ocean Research,2008, 30(1):72- 77.

[10] XIE JS,JIANG YJ, et al. Strongly nonlinear internal solution load on a small vertical circular cylinder in two- layer fluids[J].Applied Mathematical Modelling, 2010,34(8):2089- 2101.

[11] 尤云祥,李巍,時忠民.海洋內孤立波中張力腿平臺的水動力特性[J]. 上海交通大學學報, 2010,44(1):12- 17.(YOU Yunxiang, LI Wei, SHI Zhongmin,et al.Hydrodynamic Characteristics of tension leg platforms in ocean internal solitary waves[J].Journal of Shanghai Jiao Tong University,2010,44(1):56- 61. (in Chinese))

[12] 尤云祥,李巍,胡天群.內孤立波中半潛平臺動力響應特性[J]. 海洋工程,2012,30(2):1- 7. (YOU Yunxiang, LI Wei, HU Tianqun, et al. Dynamic responses of a semi- submersible platform in internal solitary waves[J]. The Ocean Engineering, 2012,30(2):1- 7. (in Chinese) )

[13] 宋志軍,勾瑩,滕斌.內孤立波作用下Spar平臺的運動響應[J]. 海洋工程,2010,32(2):12- 19. (SONG Zhijun, GOU Ying, TENG Bin, et al. The motion responses of a Spar platform under internal solitary wave [J].The Ocean Engineering, 2010,32(2):12- 19. (in Chinese) )

[14] 黃文昊,尤云祥,王竟宇,等.張力腿平臺內孤立波載荷及其理論模型[J]. 上海交通大學學報,2013,47(10):1494- 1502. (HUANG Wenhao, YOU Yunxiang, WANG Jingyu, et al. Internal solitary wave loadsand and its theoretical model for a tension leg platform [J].Journal of Shanghai Jiao Tong University, 2013,47(10):1494- 1502.(in Chinese))

[15] 關暉,蘇曉冰,田俊杰.分層流體中內孤立波數值造波方法及其應用[J]. 力學季刊, 2011, 32(2):218- 223.(GUAN Hui, SU Xiaobing, TIAN Junjie. Numerical method of solitary wave generation in a stratified flow and its applications [J]. Chinese Quarterly of Mechanics, 2011, 32(2):218- 223. (in Chinese))

[16] 高原雪,尤云祥,王旭,等.基于MCC理論的內孤立波數值模擬[J]. 海洋工程, 2012, 30(4):29- 35.(GAO Yuanxue, YOU Yunxiang, WANG Xu, et al. Numerical simulation for the internal solitary wave based on MCC theory[J]. The Ocean Engineering, 2012,30(4):29- 35. (in Chinese))

[17] 陳杰,尤云祥,劉曉東,等.內孤立波與有航速潛體相互作用數值模擬[J].水動力學研究與進展,2010,25(3):344- 350.(CHEN Jie, YOU Yunxiang, LIU Xiaodong, et al. Numerical simulation of interaction of internal solitary waves with a moving submarine[J]. Chinese Journal of Hydrodynamics, 2010,25(3):344- 350. (in Chinese))

[18] 劉碧濤,李巍,尤云祥,等.內孤立波與深海立管相互作用數值模擬[J].海洋工程,2011,29(4):1- 7.(LIU Bitao, LI Wei, YOU Yunxiang, et al. Numerical simulation of interaction of internal solitary waves with deep- sea risers[J]. The Ocean Engineering, 2011,29(4):1- 7. (in Chinese))

[19] HIRT CW, NICHOLS BD. Volume of fluid(VOF) method for the dynamics of free boundaries[J]. J.Computational Physics,1981,39(1):201- 225.

[20] 韓鵬.基于VOF方法的不規則波阻尼消波研究[D].大連:大連理工大學, 2008. (HAN Peng. The study of damping absorber for irregular waves based on VOF method [D].Dalian: Dalian University of Technology, 2008. (in Chinese))

Numerical simulation for interaction characteristics of internal solitary
waves with tension leg platform

WANG Xu1, LIN Zhongyi2, YOU Yunxiang1

(1. State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China; 2. School of Jiaxing Nanyang Profession and Technology, Jiaxing 314003, China)

According to the applicability conditions for three types of internal solitary waves theories including KdV, eKdV and MCC, a numerical method based on the Navier- Stokes equation in a two- layer fluid is presented to simulate the strongly nonlinear interaction of internal solitary waves with a tension leg platform (TLP), where the velocity- inlet boundary is applied by use of the depth- averaged velocities in the upper and lower- layer fluids induced by the internal solitary waves. Results show that the waveforms and amplitudes of the internal solitary waves based on the present numerical method are in good agreements with the experimental and theoretical results, and that the numerical results for the horizontal and vertical forces, as well as torques on the TLP due to the internal solitary waves, have good agreement with experimental results. It is shown that the horizontal and vertical forces on the TLP due to the internal solitary waves can be divided into three components which are the wave and viscous pressure forces, as well as the fractional force, which is a small amount and hence can be neglected. For the horizontal force, its main components are wave pressure and viscous pressure forces. Compared with the wave pressure force, the viscous pressure force is small but can not be neglected. For the vertical force, the component of the viscous pressure force is a small amount and hence can be neglected.

two- layer fluid; internal solitory wave; tension leg platform; load characteristics

P751

A

10.16483/j.issn.1005- 9865.2015.05.003

1005- 9865(2015)05- 016- 08

2014- 05- 20

國家自然科學基金資助項目(11372184);高等學校博士點基金資助項目(20110073130003)

王 旭(1985- ),男,博士研究生,主要從事船舶與海洋工程水動力學研究。

尤云祥。E- mail:youyx@sjtu.edu.cn

猜你喜歡
理論
堅持理論創新
當代陜西(2022年5期)2022-04-19 12:10:18
神秘的混沌理論
理論創新 引領百年
相關于撓理論的Baer模
多項式理論在矩陣求逆中的應用
基于Popov超穩定理論的PMSM轉速辨識
大電機技術(2017年3期)2017-06-05 09:36:02
十八大以來黨關于反腐倡廉的理論創新
“3T”理論與“3S”理論的比較研究
理論宣講如何答疑解惑
學習月刊(2015年21期)2015-07-11 01:51:44
婦女解放——從理論到實踐
主站蜘蛛池模板: 国产69精品久久久久孕妇大杂乱| 国产农村1级毛片| 人妻丰满熟妇αv无码| 久久综合九色综合97婷婷| 国产视频自拍一区| 婷婷五月在线| 狂欢视频在线观看不卡| 丰满人妻久久中文字幕| 伊人久久大香线蕉影院| 草逼视频国产| 日本不卡在线视频| 亚洲男人天堂久久| 欧美亚洲综合免费精品高清在线观看 | 婷婷五月在线视频| 亚洲欧美不卡中文字幕| 国产高清在线精品一区二区三区 | 青青网在线国产| 日本www在线视频| 欧美精品一区二区三区中文字幕| 亚洲人成亚洲精品| 国产小视频免费| 国产在线高清一级毛片| 欧洲亚洲欧美国产日本高清| 综合天天色| 免费无码AV片在线观看国产| 亚洲国产成人综合精品2020| 99re热精品视频国产免费| 亚洲爱婷婷色69堂| 久久永久免费人妻精品| 国产成人精品无码一区二| 国产精品片在线观看手机版| 国产精品林美惠子在线播放| 强奷白丝美女在线观看| 亚洲毛片在线看| 国产三区二区| 麻豆精品在线视频| 天天操天天噜| 亚洲成aⅴ人片在线影院八| 色网站在线免费观看| 91精品啪在线观看国产60岁| 久久精品国产亚洲AV忘忧草18| 国产乱子伦视频三区| 人妻中文字幕无码久久一区| 国内精品视频| 亚洲免费人成影院| 天天躁夜夜躁狠狠躁图片| 69av在线| 亚洲日韩精品无码专区| 亚洲视频一区在线| 97成人在线视频| 91亚洲精品第一| 国产屁屁影院| 精品撒尿视频一区二区三区| 国产精品99r8在线观看| 国内精品免费| 国产一区二区三区在线观看免费| 中文字幕欧美日韩| 全午夜免费一级毛片| 国产第一页屁屁影院| 国产成人久久综合777777麻豆| 亚洲资源站av无码网址| 日韩欧美国产另类| 999福利激情视频| 久久国产精品嫖妓| 中文字幕在线一区二区在线| 先锋资源久久| 国内精自线i品一区202| 夜色爽爽影院18禁妓女影院| 蜜臀AVWWW国产天堂| 中国国产一级毛片| 亚洲综合精品香蕉久久网| 永久免费av网站可以直接看的| 91极品美女高潮叫床在线观看| 国产精品三级av及在线观看| 99在线观看国产| 国产69囗曝护士吞精在线视频| 欧美www在线观看| 欧美成人h精品网站| AV老司机AV天堂| 国产av剧情无码精品色午夜| 538国产视频| 国产人在线成免费视频|