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

浮式結構物二階波浪力求解方法比較研究

2017-10-12 01:02:33歐紹武付世曉
海洋工程 2017年4期
關鍵詞:方法

歐紹武,付世曉

(1. 上海交通大學 海洋工程國家重點實驗室,上海 200240; 2. 高新船舶與深海開發裝備協同創新中心,上海 200240)

浮式結構物二階波浪力求解方法比較研究

歐紹武1, 2,付世曉1, 2

(1. 上海交通大學 海洋工程國家重點實驗室,上海 200240; 2. 高新船舶與深海開發裝備協同創新中心,上海 200240)

對遠場法、近場法和中場法這幾種常用的二階波浪力計算方法進行了總結,并結合算例對它們之間的差異進行了研究,重點分析了網格密度、浮體形狀、水深等因素的影響。此外,還對全QTF法和Newman近似這兩種不規則波中的二階力計算方法進行了討論,比較了不同水深條件下兩種方法求得的二階波浪力譜之間的差異。研究結果表明,遠場和中場法給出的結果基本一致,且具有較好的數值穩定性;而近場法受網格密度、浮體形狀等因素影響較大;Newman近似給出的結果在頻率較低時與全QTF法接近,但在頻率較高時存在一定的誤差。

二階波浪力;遠場法;近場法;中場法;全QTF法;Newman近似

Abstract: The commonly used methods in calculating second-order wave drift forces, such as far-field integration, near-field integration and middle-field integration, are summarized. The difference between these methods are studied by several examples, with special emphasis on element density, model shape and water depth. In addition, the drift forces in irregular sea with different water depths are calculated by both full-QTF method and Newman’s approximate, and the difference between two methods is further discussed through spectrum analysis. The results show that far-field integration and middle-filed integration give similar results with better robustness, while near-field integration is largely influenced by element density and model shape. It also shows that the results given by the Newman’s approximation are close to those by the full QTF method at low frequencies, but with certain errors when the frequency is higher.

Keywords: second-order drift force; far-field integration; near-field integration; middle-field integration; full-QTF; Newman’s approximate

對于海洋結構物而言,在不規則波浪的作用下,受到的波浪力包含一階和二階成分。一階波浪力幅值較大,與波高成線性關系且頻率與波浪相同。而二階波浪力包括定常二階力、低頻二階力和高頻二階力三部分,幅值相對較小,且與波高的平方成比例。由于二階波浪力的存在,使得波浪中的海洋結構物不僅存在波頻的搖蕩運動,而且還存在長周期的漂移運動,這一運動的頻率遠低于不規則波的特征頻率。對于錨泊的海洋浮式結構物如FPSO、半潛平臺等而言,二階波浪力中的低頻成分有可能與錨泊系統的固有頻率接近而產生共振,從而引起相當大的水平運動,在錨泊系統中產生較大的附加應力,對錨泊系統的安全構成威脅。

對二階波浪力的深入研究源于20世紀60年代。Maruo[1]在由自由面、物面以及一個遠離浮體的控制面圍城的流域內,應用動量和能量守恒方程,得到了浮體在規則波中的平均二階力的表達式。該理論之后由Newman[2]及Faltinsen等[3]加以推廣。上述方法均是在遠離浮體的控制面上積分,具有計算簡便的優點,但只能求出二階力中的定常部分,被稱為遠場積分法。后來,Molin[4], Pinkster[5], Ogilvie[6]等提出了壓力積分法,直接求出物面上壓力的所有二階分量,然后在物面上積分,求得二階力在各坐標軸上的分量。這種方法又被稱為近場積分法,它可以清晰地顯示出二階力的各種組成部分,能夠用于計算定常、低頻和高頻二階力。但與遠場積分相比,該方法計算相對復雜費時,計算結果容易受物面形狀及網格劃分的影響。近年來,Chen[7]在近場積分法的基礎上提出一種新的二階力計算方法,被稱為中場積分法。其做法是首先在浮體附近建立一個控制面,然后利用斯托克斯公式和高斯公式,將近場積分法公式中對物面的積分轉化成在控制面及平均水線面上的積分。這一方法能夠用于計算定常、低頻和高頻二階力,同時,由于控制面可以選取解析的性狀,避免了浮體表面網格劃分對速度勢求解精度的影響,具有較好的數值穩定性。

對不規則波中的浮式錨泊系統而言,低頻二階力預報是十分重要的。對于低頻二階力的計算,可以采用近場積分法或者中場積分法求解完整的二階傳遞函數矩陣,被稱為全QTF法,但這些方法需要求解二階速度勢,求解過程極其耗時。工程中常使用Newman[8]提出的利用定常二階力近似低頻二階力的方法,該方法被稱為Newman近似。

在錨泊浮式海洋結構物的運動響應預報中,選取合適的方法,對二階波浪力進行準確考慮是十分重要的。在以往的研究中,有學者曾分別對上述方法中的某一種或兩種的特點進行過討論[9-12],但這些研究均是基于某一特定的船型開展的,得出的結論存在一定的局限性。本文首先對上述幾種常用的二階力計算理論進行了梳理,隨后結合一些算例就網格劃分、物體形狀、水深等因素對這幾種方法計算結果的影響進行了分析,并在此基礎上對它們的適用性進行了討論。

1 基本理論

1.1勢流理論

為了討論浮體在波浪中的運動,這里定義三個坐標系[13],分別是:大地坐標系OXYZ,原點位于靜水面上,用于描述波浪入射角β;隨體坐標系OAXAYAZA,原點固結在浮體重心處,隨船一起搖蕩;參考坐標系OBXBYBZB,原點始終位于浮體的平衡位置,不隨浮體一起搖蕩,當浮體靜止時,其與隨體坐標系重合。三個坐標系均為右手系,Z軸垂直向上,浮體靜止時,三個坐標系的X,Y軸平行,如圖1所示。

圖1 坐標系及浪向角定義Fig. 1 Definition of coordinate system and wave direction

本文討論的浮體二階波浪力計算基于勢流理論,流體滿足無黏、無旋、不可壓縮三個基本假定,當簡諧的入射波頻率為ω時,總的速度勢可以表示為:

其中,φI是入射勢;φD為繞射勢,是固定不動的結構物在入射波作用下產生的速度勢;φR表示輻射勢,是浮體在靜水中搖蕩產生的速度勢。φp=φD+φR是由于結構物存在而產生的速度勢,又稱為擾動勢。

速度勢滿足下列定解條件:

1)Laplace方程:

2φ=0 (流域內)

2)統一的線性化的自由面條件:

3)物面條件:

其中,(n1,n2,n3)=n,(n4,n5,n6)=r×n,n表示浮體表面某點的單位外法向量;r表示從浮體隨體系原點到該點的矢徑。

4)水底條件:

5)遠方輻射條件:

假設海底是水平的,自由表面無限長,根據Airy波理論可以得到入射勢的表達式為:

其中,A為波幅,β為浪向角,h為水深,k為波數。由色散關系確定:

輻射勢和繞射勢的定解問題可以利用格林函數法求得。在柱坐標系(R,θ,z)下,距離浮體很遠(R→∞)的某點的繞射勢和輻射勢可以表示為[7]:

其中,K(θ)稱為Kochin函數,定義為:

式中:S為浮體的濕表面;σ(Q)為表面的源強。

1.2遠場法

遠場法的基本思想是在由物面S、自由面SF、半徑無限大的垂直柱面S∞和水底平面SB圍成的流域內,運用動量和能量守恒方程,得到浮體所受的平均二階力。Newman[2]給出的平均二階力的表達式為:

從表達式來看,遠場法是在一個柱面上積分,對浮體表面網格劃分不敏感,因此具有較好的數值穩定性,在工程中被廣泛應用。但同時,由于遠場法需要選取無限大控制面,故只能用于計算單浮體的平均二階力。

1.3近場法

近場法的基本思想是直接求出物面上壓力所有的二階分量,并在物面上進行積分,以求得二階波浪力(力矩)在各坐標軸上的分量。具體做法是將速度勢、壓力、瞬時波面升高、運動向量等均展開成ε的級數的形式,然后在浮體的瞬時濕表面上積分,并精確至二階[13]:

注意式中已將浮體的瞬時濕表面分成平均濕表面S0,和因波浪及浮體運動引起的濕表面變化ΔS。式中帶括號的上標(0)、(1)、(2)表示對應物理量的階次。

p(0),p(1),p(2)分別為靜壓力、一階動壓力和二階動壓力,由下面的式子確定:

注意這里為了避免混亂,引入了符號x1,x2,x3來表示x,y,z三個方向。類似的,后面還將使用α1,α2,α3來表示α,β,γ三個歐拉角。式中,向量X表示浮體表面一點在參考坐標系下的位置,定義為:

浮體表面法向N的定義與X類似,為:

N=n+α(1)×n+D·n=n(0)+εN(1)+ε2N(2)

將上面式子代入式(18),保留二階項可以得到二階漂移力的表達式如下:

在不規則海況中,入射波可以記為各種頻率的規則波的疊加,即:

若只研究差頻二階力,則可將表達式化簡成以下形式:

從近場法的推導過程可知,該方法具有普適性,能應用于單體或者多體系統,并能夠給出六個自由度上包括定常、低頻和高頻成分在內的所有二階力分量。但由于近場法需要在船體表面積分,在船體表面存在尖角的地方容易產生奇點,進而對計算結果造成影響。

1.4中場法

中場法從近場法發展而來,其基本思想是在浮體附近建立一個封閉的控制面,然后對近場法的二階漂移力計算公式使用雙參數的斯托克斯公式和高斯公式,將其轉換成在浮體表面S、自由面SF及控制面C上的積分[7],即:

F=ρ?S[Φ(·n)+()(X·n)]ds-ρg∮ΓΘkdl-

M=ρ?Sx×[Φ(·n)+()(X·n)]ds-ρg∮ΓΘ(x×k)dl-

式中:Γ為浮體與自由面的交線;ΓC為控制面與自由面的交線;Θ定義為:

式(22),(23)是通用的公式,可用于計算包括和頻和差頻成分在內的六個方向上的所有二階力。若只需計算平面內三個方向(橫蕩、縱蕩和艏搖)的平均二階力,則可將上述公式進一步簡化成在控制面及其與自由面的交線上的積分,即:

式中:〈·〉表示在一個周期內取平均值。

從表達式來看,中場法將積分轉換成了在控制面及其與自由面交線上的積分,受物體幾何形狀影響較小,能夠應用于單體或者多體的情況,并能給出六個自由度上的定常、低頻及高頻二階力。

1.5Newman近似

對不規則波中的浮式錨泊系統而言,低頻二階力預報是十分重要的。對于低頻二階力的計算,常用的方法主要有兩種:第一種是利用近場法式(19)或者中場法式(22)、式(23)直接計算完整的QTF矩陣,然后根據不規則波時歷計算低頻二階力,被稱為全QTF法;另一種方法是根據規則波下計算得到的平均二階力,使用近似的方法得到QTF矩陣,其中使用最廣泛的是Newman近似。

一般而言,ωi,ωj平面內的所有頻率均對二階力有貢獻,但二階力的幅值通常都比較小,因此只有頻率接近錨泊系統固有頻率的二階力會對浮體的運動響應有明顯影響。頻率差Δω等于固有頻率ωn在ωi,ωj平面內表示兩條線ωi=ωj±ωn,而大多數情況下,浮式錨泊系統在水平面內運動的固有頻率都遠小于波浪頻率,此時這兩條線接近QTF矩陣的對角線。Newman近似[8]認為,這種情況下QTF矩陣中的非對角元素Tωi,ωj可以用對角元Tωi,ωi來表示:

上式是原始的Newman近似公式,該式還有多種修正的表達式,其中使用最多的是[14]:

與全QTF法相比,Newman近似能夠極大地提高計算效率,在浮式錨泊系統固有頻率很低的情況下有較高的精度,但其完全忽略了二階速度勢的影響,在某些情況下的正確性有待討論。

2 數值計算

2.1數值模型

本文的數值計算主要基于Hydrostar軟件[15],共涉及到兩種主尺度相近的模型,分別是方Weigly船和半橢球體,其主要參數如表1所示。

表1 模型的主要參數Tab. 1 Models’ main characteristic parameters

2.2浮體形狀及網格劃分的影響

為了研究網格劃分對不同二階力計算方法的影響,選取網格數為400、800、1 200、2 000和4 000五種情況,進行比較分析。數值計算中的波浪頻率選取為0.05~2.5 rad/s,間隔0.05 rad/s,浪向角選取45°,浮體濕表面及中場積分法所用的控制面如圖2所示。

圖2 濕表面及控制面模型(網格數N=1 200)Fig. 2 Panel model of control surface and wet surface (element number N=1 200)

圖3給出了不同網格密度下Weigly船的一階縱蕩、橫蕩和艏搖運動RAO。其中,無符號的曲線表示縱蕩;帶“o”的曲線表示橫蕩;帶“*”的曲線表示中場法艏搖。從計算結果來看,各網格密度下的結果吻合得很好,表明網格密度對一階運動RAO沒有明顯的影響。

圖3 不同網格密度下Weigly船和橢球體的運動RAOFig. 3 Motion RAOs of weigly ship and ellipsoid with different element densities

圖4給出了不同網格密度下,分別用遠場法、近場法和中場法計算平均二階力的結果。其中,帶“x”的曲線表示遠場法;帶“o”的曲線表示中場法;無符號的曲線表示近場法。

圖4 不同網格密度下Weigly船和橢球體的平均二階力Fig. 4 Mean drift force of Weigly ship and ellipsoid with different element densities

總體而言,三種方法給出結果隨頻率的變化趨勢相同,中場法和遠場法給出的結果基本一致,且具有較好的收斂性。對于算例中的兩個模型,當網格數大于1 200時,中場法和遠場法的結果基本不再隨網格數變化而變化;但對于近場法,在波浪頻率較大時,不同的網格密度下計算結果有較大的差異,隨著網格密度的提高,近場法與其它兩種方法的結果有趨同的傾向。此外,在波浪頻率較低時,三種方法的結果基本一致,但隨著波浪頻率的增大,近場法的結果與遠場法和中場法的結果出現了一定的差異。值得注意的是,在波浪頻率為0.75~1 rad/s時,近場法給出的Weigly船的縱蕩二階力與中場法和遠場法的結果有較大的差異,其原因可能是Weigly船的表面存在尖角,進行壓力積分時這些地方容易產生奇點,進而對計算結果造成影響。相比之下,表面相對光滑的橢球體則不存在這類問題,波浪頻率較低時(1.2 rad/s以下)三種方法的結果基本一致。

2.3水深的影響

一般而言,不同水深下浮體所受的波浪力及其運動響應存在一定的差異,常用的三種二階波浪力計算方法的原理不同,其在不同水深條件下給出的預報結果可能存在一定的差異。為了比較這種差異,這里以橢球體為例,網格數選取N=2 000,水深與吃水之比選取無限大(inf),10,5,2,1.5五種情況進行分析,如圖5所示。其中,帶“*”的曲線表示遠場法;帶“o”的曲線表示中場法;無符號的曲線表示近場法。從圖5中的結果來看,在各種水深條件下,遠場法給出的結果均與中場法一致。而對于近場法,由于算法的不同,在算例中的幾種水深條件下,給出的平均二階力結果與遠場法和中場法略有差異,但隨水深變化的規律與基本一致,即隨水深的變淺而有一定程度的增大。這一現象表明,在淺水條件下預報平均二階力時,三種方法均適用。

圖5 不同水深下橢球體的平均二階力Fig. 5 Mean drift force of ellipsoid at different water depths

2.4Newman近似與全QTF法的差異

由于Newman近似簡單地根據二階傳遞函數的定常項(對角元項)得到QTF矩陣的其余部分,這種做法難免會帶來一定的誤差。

圖6 橢球體的橫蕩QTF矩陣Fig. 6 Sway QTF matrix of ellipsoid

為了研究Newman近似與全QTF法的差異,這里以橢球體為例,網格數選取N=2 000,水深/吃水選取無限大,利用中場法求得平均二階力,然后利用Newman近似求出近似的QTF矩陣,如圖6(a)所示;圖6(b)則給出了直接用中場法求出的完整的QTF矩陣。由圖可知,在對角線附近(圖中兩根虛線之間的位置),兩種方法給出的結果基本一致,表明兩種方法對二階力中頻率很低的部分預報結果是相近的。但在遠離對角線的位置,兩者存在明顯的差異,Newman近似得出的QTF矩陣較為平滑,具有明顯的規律性;而全QTF法的結果在非對角元上波動較大,存在明顯的峰谷交替現象。這一結果表明,兩種方法對頻率稍高的二階力成分預報的結果可能存在差異。

為了進一步說明兩種方法的預報結果的差異,根據式(21)計算了不同水深條件下,不規則波中橢球體所受的二階力時歷,并進一步利用FFT給出了相應的二階力譜,如圖7所示。其中不規則波選取ITTC雙參數譜,有效波高Hs=1 m,特征周期T1=4.6 s,其譜形如圖8所示。

圖7中H/d表示水深與吃水之比,其中虛線表示Newman近似的結果,帶“o”的實線表示全QTF法的結果。顯然,在頻率較低的情況下(小于0.1 rad/s),Newman近似與全QTF法給出的結果基本一致;但在稍高的頻段內(0.1~0.8 rad/s),兩種方法存在明顯的差異,大多數情況下,Newman近似的結果明顯大于全QTF法的結果。此外,在水深較小的情況下(H/d=1.5),全QTF法的結果中,0.2~0.6 rad/s范圍內的艏搖二階力明顯大于Newman近似的值,原因是淺水條件下二階速度勢對二階力有較大的貢獻,而Newman近似完全忽略了這一項的影響,從而導致較大的誤差。

圖7 橢球體的橫蕩二階力譜Fig. 7 Sway drift spectrum of ellipsoid

圖8 ITTC雙參數譜(Hs=1 m, T1=4.6 s)Fig. 8 Two-parameter ITTC spectrum (Hs=1 m, T1=4.6 s)

3 結 語

通過上述討論,可以得出以下結論:

1) 遠場積分法只能用于單浮體計算,能給出縱蕩、橫蕩和艏搖三個自由度上的平均二階力;近場積分法清楚地解釋了二階力的來源,具有普適性,能應用于單體或者多體系統,并能夠給出六個自由度上包括定常、低頻和高頻成分在內的所有二階力分量;中場積分法能應用于單體或者多體系統,也能用于計算二階力的定常、低頻和高頻部分。

2) 三種方法給出的結果在波浪頻率較低時(1.2 rad/s以下)差異不大。但隨著波浪頻率的提高,遠場法與中場法表現出更好的數值穩定性,網格密度對結果影響不大;相比之下,近場法的結果受網格密度影響比較明顯。此外,使用近場法計算具有尖角的結構時可能會出現誤差。

3) 三種方法給出的平均二階力隨水深變化的規律一致,即隨水深的變淺而有一定程度的增大,表明這三種方法對深水或淺水的情況均適用。

4) Newman近似與全QTF法給出的結果在對角線附近基本一致,但在遠離對角線的位置,兩者存在明顯的差異。從不規則波下的響應譜來看,在頻率很低(0.1 rad/s以下)的情況下,Newman近似與全QTF法給出的結果基本一致;但隨著頻率的增大,Newman近似出現了明顯的誤差。此外,在水深較淺的情況下,二階速度勢對二階力有較大貢獻,導致浮體某些自由度的低頻二階力將會發生明顯的變化,但Newman近似忽略了這一影響。

[1] MAURO H. The drift of a body floating on waves[J]. Journal of Ship Research, 1960, 4: 1-5.

[2] NEWMAN J N. The drift force and moment on ships in waves[J]. Journal of Ship Research, 1967, 11(1): 51-60.

[3] FALTINSEN O M, MICHELSEN F C. Motions of large structures in waves at zero froude number[J]. International Symposium on the Dynamics of Marine Vehicles and Structures in Waves, 1975, 90:3-18.

[4] MOLIN B. Computations of wave drift forces[C]//Proceedings of the offshore Technology Conference. 1979: OTC-3627-MS.

[5] PINKSTER J A. Low frequency second order wave exciting forces on floating structures[D]. TU Delft, Delft University of Technology, 1980.

[6] OGILVIE T F. Second-order hydrodynamic effects on ocean platforms[M]//International Workshop on Ship and Platform Motions. Berkeley: University of California, 1983.

[7] CHEN X B. Middle-field formulation for the computation of wave-drift loads[J]. Journal of Engineering Mathematics, 2007, 59(1): 61-82.

[8] NEWNAN J N. Second-order, slowly-varying forces on vessels in irregular waves[J]. Marine Vehicles, 1974: 182-186.

[9] KASHIWAGI M, ENDO K, YAMAGUCHI H. Wave drift forces and moments on two ships arranged side by side in waves[J]. Ocean Engineering, 2005, 32(s 5-6):529-555.

[10] INOUE Y, ISLAM M R. Numerical investigation of slowly varying drift forces of multiple floating bodies in short crested irregular waves[C]//Proceedings of the Tenth International Offshore and Polar Engineering Conference. International Society of Offshore and Polar Engineers, 2000.

[11] MALENICA &, OROZCO J M, CHEN X B. Some aspects of multibody interactions in seakeeping[C]//Proceedings of the Fifteenth International offshore and Polar Engineering Conterence. International Society of Offshore and Polar Engineers, 2005: ISOPE-1-05-250.

[12] NACIRI M, BUCHNER B, BUNNIK T, et al. Low frequency motions of LNG carriers moored in shallow water[C]//Proceedings of the 23rd International Conference on Offshore Mechanics and Arctic Engineering. 2004:995-1006.

[13] 劉應中,繆國平. 船舶在波浪上的運動理論[M]. 上海: 上海交通大學出版社, 1987. (LIU Y Z, MIAO G P. Motion theory of ship on waves[M]. Shanghai: Shanghai Jiao Tong University Press, 1987. (in Chinese))

[14] DNV-RP-C205, Environmental conditions and environmental loads[S]. Veritas D N, 2010.

[15] Veritas B. Hydrostar for experts user manual[S]. 2010.

A comparative study on the numerical methods of second-order drift forces on floating structures

OU Shaowu1, 2, FU Shixiao1, 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 200240, China)

U661.1

A

10.16483/j.issn.1005-9865.2017.04.013

1005-9865(2017)04-0100-10

2016-08-15

歐紹武(1991-),男,碩士研究生,主要從事多體水動力研究。E-mail:shaowuloveh@163.com

付世曉。E-mail:shixiao.fu@sjtu.edu.cn

猜你喜歡
方法
中醫特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數學教學改革的方法
河北畫報(2021年2期)2021-05-25 02:07:46
化學反應多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學習方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 国产精品冒白浆免费视频| 日本a∨在线观看| 欧美午夜网站| 欧美亚洲欧美区| 午夜欧美在线| 激情在线网| 国产91色在线| 国产97区一区二区三区无码| 国产欧美日韩va| 尤物成AV人片在线观看| 99视频在线精品免费观看6| 香蕉视频在线观看www| 亚洲毛片网站| 四虎成人在线视频| 亚洲动漫h| 免费人成在线观看成人片| 国产成人精品亚洲日本对白优播| 国产亚洲精品无码专| 国产丝袜无码一区二区视频| 搞黄网站免费观看| 免费大黄网站在线观看| 亚洲区第一页| 91久草视频| 国产哺乳奶水91在线播放| 欧美日韩精品一区二区视频| 98超碰在线观看| 久热中文字幕在线| 蜜臀AVWWW国产天堂| 2021国产精品自产拍在线| 精品综合久久久久久97超人| 午夜高清国产拍精品| 亚洲丝袜中文字幕| 少妇精品网站| 网久久综合| 亚洲av片在线免费观看| 欧美一级高清视频在线播放| 毛片网站观看| 国产高清在线精品一区二区三区| 亚洲无码视频喷水| 久久激情影院| 色婷婷狠狠干| 国产91色在线| 亚洲精品福利网站| 美女扒开下面流白浆在线试听| 亚洲成A人V欧美综合| 色妞永久免费视频| 亚洲人人视频| 国产拍揄自揄精品视频网站| 欧美中文字幕无线码视频| 国产91全国探花系列在线播放| 无码精品一区二区久久久| 国产SUV精品一区二区6| 丝袜美女被出水视频一区| 成年av福利永久免费观看| 亚洲一区二区成人| 国产精品综合久久久| 国产日韩欧美中文| 久久精品国产一区二区小说| 老熟妇喷水一区二区三区| 在线免费a视频| 欧美亚洲国产日韩电影在线| 91青草视频| 91视频精品| 亚洲国产精品日韩av专区| 国产啪在线91| 91一级片| jizz在线免费播放| 成色7777精品在线| 波多野结衣一级毛片| 久久婷婷综合色一区二区| 九九久久精品免费观看| 亚洲午夜福利在线| 91在线激情在线观看| 色婷婷啪啪| jizz国产视频| 国产成人综合日韩精品无码不卡| 国产精品香蕉在线观看不卡| 亚洲国产高清精品线久久| 最新痴汉在线无码AV| 国产青青草视频| 黄网站欧美内射| 香蕉网久久|