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

淺水效應對船舶阻力及流場特性的影響分析

2017-05-10 12:34:17孫帥王超常欣支玉昌
哈爾濱工程大學學報 2017年4期
關鍵詞:船舶

孫帥, 王超, 常欣, 支玉昌

(哈爾濱工程大學 船舶工程學院,黑龍江 哈爾濱 150001)

?

淺水效應對船舶阻力及流場特性的影響分析

孫帥, 王超, 常欣, 支玉昌

(哈爾濱工程大學 船舶工程學院,黑龍江 哈爾濱 150001)

為了分析淺水條件下船舶阻力及流場的特點,本文基于混合網格技術,結合Reynolds average numerical simulation(RANS)方法,對標準船模KRISO Container Ship (KCS)開展了船模淺水效應的數值預報分析。KCS船的船艉部分使用了三種不同的非結構網格劃分方法,與水池試驗值對比得到恰當的網格劃分方法。通過改變數值水池高度探究水深變化對于船舶阻力及流場特性的影響。研究發現:當水深h小于10倍吃水時,KCS船淺水效應明顯,隨著水深的減小,船舶阻力系數、z方向抽吸力以及槳盤面伴流分數均呈增大趨勢;隨著水深的減小,艏艉壓力差逐漸增大,船艉傾現象更加明顯。

混合網格;KCS船模;淺水效應;船舶阻力;抽吸力;伴流場;RANS

當貨船航行于港口、運河、近海等位置時,貨船就處于限制航道狀態,船舶在淺水中有著與在深水中不同的水動力性能,研究淺水中船舶的阻力及流場特性的規律,無疑對船舶的性能預報、淺水修正有著重大的意義。

船舶在限制航道中航行的各項性能研究主要有試驗方法、勢流理論以及粘流(computational fluid dynamics,CFD)方法。1985 年E.Müller利用有限水深興波源勢方法,通過在船舶縱向剖面上布置源匯,提出了薄船、線性自由表面條件下船舶水壓場的計算公式[1]; I.Sahin等采用有限水深格林函數方法,系統開展了潛艇、水面艦艇和氣墊船水壓場的理論建模和數值計算工作[2-3];郭圣漢在模型試驗中采用了隨車假底試驗技術,該方法能夠在水深吃水比h/T小于4時,對Fn<0.2的工況進行較好的預報模擬,在h/T>4時可以使用到Fn=0.35的范圍[4];李一兵通過對實船與船模試驗資料的分析計算,對航行于不同內河航道中船舶的航行阻力進行了分析探討,得出對于不同航行河流應該有不同的阻力估算方法,并提出了相應的計算公式[5];繆濤等利用面元法計算了有限水深船舶在水底引起的壓力分布[6];近年來,人們開始采用CFD方法進行淺水阻力的研究,與試驗值進行對比發現,CFD方法能夠得到較為準確的阻力計算結果[7-10]。

目前淺水航行船舶的阻力及流場特性大都采用不同水深的模型水池試驗和勢流理論的方法,模型試驗結果需進行實尺度換算,這個過程需要消耗大量資源,并且需要進行試驗次數較多,而且在換算過程中很難保證換算結果的準確性;勢流理論方法忽略了粘性作用,結果準確性較差,且得不到具體的流場信息;CFD計算大多僅就淺水阻力性能進行了分析,并沒有對流場信息等作進一步討論。本文將通過CFD方法計算得到船模在深水和淺水中的阻力及流場特性,分析船舶阻力、抽吸力以及伴流場信息隨水深的變換規律,探討產生淺水效應的內在原因。

1 計算模型的建立

流體流動要受物理守恒定律的支配,基本的守恒定律包括:質量守恒定律、動量守恒定律和能量守恒定律。由于水為不可壓縮流體,熱交換很小以至于可以忽略不計,可只對質量守恒方程和動量守恒方程進行求解。詳細公式可參考文獻[11]。計算中采用k-ωSST湍流模型,該模型有效集成了k-ε和k-ω模型的優點,能夠較好地模擬存在流動分離和強逆壓梯度的復雜流動問題。

在計算過程中監測船體阻力系數,計算初始階段動量相、湍流動能和耗散率都使用一階格式,觀測阻力系數的變化,當阻力系數在比較小的范圍內變化時改為二階格式,當阻力系數穩定以后認定計算收斂。

1.1 計算對象

本文研究對象為KCS集裝箱船,如圖1所示,模型縮尺比為31.6,其船模具體參數如表1所示。

圖1 KCS船模的幾何形狀
Fig.1 The geometry of KCS

表1 KCS船模主要參數

1.2 計算域的建立及計算參數設置

流場計算區域為:入口距船艏1LPP,出口距船艉2LPP,側面和底面均距船體表面1LPP。如圖2所示。邊界條件的入口采用速度入口,速度設為2.196m/s,側面和底面采用與入口相同速度的移動壁面,船體表面使用固定壁面,頂面使用對稱面,出口采用壓力出口,設定出口壓力為10kPa。壓力速度耦合迭代采用SIMPLEC方法。

圖2 計算域及邊界條件Fig.2 The computational domain and boundary conditions

1.3 網格劃分

由于KCS船艉部型線較為復雜,艉部難于生成高質量的結構化網格,因此本次計算采用混合網格[12]方法,船艉小部分流域劃分非結構化網格,以節省網格劃分工作量,而其他部分采用六面體結構化網格,結構化網格和非結構化網格之間通過interface連接,流場通過interface插值進行信息傳遞。

流體的邊界層區域可以分為層流底層、湍流邊界層和阻流區,第一層網格的厚度對于數值模擬的結果尤為重要,因此有了無因次的第一層網格厚度y+。一般來說,模型尺度y+建議控制在60左右是合適的,對應第一層網格厚度為0.8~1 mm。結構化網格區域第一層層網格厚度可以直接用參數定義,而船艉非結構化網格則需要通過添加多層棱柱網格完成。

圖3 船體周圍加密網格Fig.3 The refined mesh around ship

對于船艉非結構網格部分,如圖4所示分別嘗試使用了三種網格劃分方法:方案1為常規方法,圖4(a),僅對船體表面網格進行了加密處理;方案2考慮到第一層網格厚度對于物體流體邊界層的模擬比較重要,圖4(b),對船體表面附近添加棱柱網格,可以加強對船艉部分邊界層細節的模擬,以利于對船艉流場計算的準確性;方案3考慮到本文計算的一個考察重點為槳盤面速度場的模擬,因此在槳盤面和槳盤面前方的位置添加密度體,圖4(c),通過使用密度體可以得到槳盤面流場更加細致伴流場信息。

圖4 三種非結構網格劃分方法Fig.4 Three unstructured grid methods

無因次化的阻力系數Ct,其公式為

(1)

式中:RT代表船模受到的阻力;U0代表船模航速,2.196 m/s;ρ為水的密度,998.2kg/m3;SW為船模在正浮狀態下的濕表面積,9.438 m2。

表2 三種網格劃分方法阻力系數對比

三種船艉網格劃分方法的阻力系數對比如表2所示。對比三種網格劃分方法可以發現,三種方法的計算準確度依次提高,這是由于棱柱網格使船體表面網格加密,更好地模擬了物面的邊界層流動,在槳盤面處添加的密度體與船體有一定的距離,密度體的添加只是使密度體位置的流場分布更加精細,對船體附近的邊界層阻力等沒有太大影響,這就導致了添加密度體對于阻力計算結果影響不大。

伴流的大小通常用伴流分數ω來表示,軸向伴流分數ωx的表達式為

(2)

式中:ωx為槳盤面的軸向伴流分數,Vx為槳盤面的軸向速度,VA為來流速度。

圖5 不同網格方案計算結果與試驗值Fig.5 The comparison of different grid methods between calculation and experiment

圖5為槳盤面不同半徑處(r為螺旋槳半徑),三種網格方案計算結果與試驗結果。將圖中三種網格方案計算結果與試驗測量值對比可以發現,方案1網格的計算結果最差,沒有模擬出槳盤面上半部分中縱剖面附近伴流分數的一個小波峰。添加邊界層后使計算結果變好,中縱剖面附近伴流分數出現變化,但不是很明顯,且在內半徑處,軸向伴流分數計算值與試驗值偏差很大。在槳盤面處添加密度體之后,出現明顯的小波峰,除了在內半徑180°左右區域內的軸向伴流分數較小之外,盤面處其他區域的計算值與實驗值吻合很好,為得到詳細準確的船艉伴流場信息,本文考慮使用方案3劃分船艉部網格。

1.4 計算工況介紹

根據國際船模試驗池會議規定,試航時不產生淺水效應的水深h必須滿足:

(3)

h>2.75×V2/g

(4)

根據表1給出的船體參數,按照式(3)計算得到的水深需要大于 5.178m,根據式(4)計算得到水深需要大于1.326m,綜合兩個公式可以暫時認定水池試驗水深大于5.178m可以消除淺水的影響,h=5.178 m的情況相當于h/T=15.15。如圖6所示將計算域中的水深分別設定為h/T=21.3、h/T=15.15、h/T=10、h/T=4、h/T=3以及h/T=2六種工況,使用Fluent軟件計算六種工況下的KCS船的阻力及流場特性,并進行對比分析。

圖6 船體吃水與水深示意圖Fig.6 The general view of draft and depth

2 計算結果分析

2.1 淺水阻力對比分析

船舶航行在淺水中會產生淺水效應,在淺水中航行時,船舶阻力受到的影響最大,因此首先對比六種工況的船體阻力系數。

阻力系數隨水深的變化曲線如圖7所示。觀察圖中的數據信息,當航行水深大于10倍吃水時,船舶阻力受水深變化很小,可以認為此時不發生淺水效應,根據式(3)國際船模水池的規定的估算公式,是足夠保證避免淺水效應的。當水深小于10倍吃水時,阻力系數隨吃水的減小單調增加,當水深減少到4倍吃水時,船舶阻力增幅變大,此時,船舶阻力隨水深變化明顯,淺水效應顯著,當水深繼續減小到2倍吃水時,阻力系數陡增,相比于3倍吃水時,增幅達12.4%,而相比于深水航行時(大于10倍吃水),增加幅度達到32.7%。可見當水深較淺時,淺水效應對船舶阻力性能有較大影響。

圖7 不同水深阻力系數及z向抽吸力Fig.7 The resistance coefficient and suction force in different depth

2.2抽吸力對比分析

船舶淺水航行時,船底與航道底部間隙變窄,船底水流速度增加,由伯努利方程可知船底部的流體流速增大、壓強減小,船底受到的z方向(垂直于水線面朝下)抽吸力增大,如圖6所示。從圖7中可以發現,z向抽吸力變化趨勢與船舶阻力隨水深的變化趨勢基本一致:當航行水深大于10倍吃水時,z向抽吸力基本不變;當航行水深小于10倍吃水時,z向抽吸力隨水深變化明顯。由此可見,船舶航行于淺水航道時,船舶阻力增加與z向抽吸力的增加息息相關。

圖8給出了不同水深時,船底中縱剖線壓力系數沿船長方向的變化,其中X/L=1.0為船艏位置。從圖中可以看出,隨著水深的減小,船艏部壓力系數逐漸增大,艉部壓力系數逐漸減小,艏艉的壓力差逐漸增大,說明隨著水深的減小,船艉尾傾現象更加明顯。

圖8 船底中縱剖線壓力沿船長方向的變化Fig.8 The change of pressure in the longitudinal lines along the length direction

2.3 伴流場對比分析

槳盤面伴流場信息是考察船舶裸船體水動力性能中一個重要指標,對于單槳船來說,螺旋槳工作于船的后方,其運轉效果很大程度上會受到船體的影響,現在的螺旋槳設計工作多數是基于標稱伴流場的螺旋槳設計,因此考察水深對槳盤面伴流場的影響就尤為重要。測量表明,與軸向伴流速度相比,周向和徑向伴流是二階小量,在螺旋槳設計中,通常不予考慮。因此本文只考慮水深變化對軸向伴流的影響。通過式(2)可以得到標稱伴流分數[13],然后以此為基礎估算船舶的實效伴流分數,實效伴流為計及了船槳之間相互干擾的船艉槳盤面處伴流,基于實效伴流分數進行螺旋槳的適伴流設計工作,適伴流設計槳的整體性能要優于圖譜設計槳,達到耗能較少,減振降噪的目的。

圖9 不同水深下的標稱伴流分數分布Fig.9 Distribution of nominal wake fraction in different depth

KCS為單槳船,槳盤面關于船中左右對稱,為方便對比,截取如圖9所示的速度云圖進行分析。對比四個水深的軸向速度可以看到在速度分布的整體趨勢基本一致,在槳盤面中心形成一個封閉的高伴流區,等值線向船體方向延伸,越靠近船艉,軸向伴流越大。從圖9(a)可以觀察到,水深為21.3T和10T時,槳盤面軸向伴流速度云圖分布基本沒有差別,說明此時沒有發生淺水效應,這與之前的分析一致。從圖9(b)可以看出,當船舶從深水航行到淺水航道,伴流等值線向外擴散,在槳盤面上方ωx=0.6的等值線有一個偏上的尖銳區域,說明水深減小以后在槳盤面處的伴流分數增大,槳盤面中心與船體之間變化明顯。對比圖9(c)航行水深為4T和3T時,隨著水深的減小,伴流等值線繼續向外擴散,曲線上方開口逐漸增大,而在上部ωx=0.6的曲線由閉合狀態改變為在槳盤面邊界靠近船體部分出現開口,整個盤面的平均伴流分數繼續增大。繼續觀察圖9(d)可以看到當航行水深為2T時,槳盤面不存在ωx=0.1曲線,表明在h=2T時槳盤面的軸線伴流分數比前面三種工況有明顯的擴大,與阻力系數曲線相對照可以發現伴流分數與阻力系數曲線是可以相互印證的,槳盤面上半部分具有明顯特征的ωx=0.6曲線與前文分析的趨勢相同,上部分的開口有所擴大。從四個水深狀況的軸向伴流分數云圖對比中可以分析看出,隨著水深的減少,軸向伴流分數有增大的趨勢,槳盤面中心與船體之間部分,伴流分數受水深影響較大,四種工況的槳盤面平均軸向伴流分數如表3所示。

表3 不同水深下槳盤面伴流分數對比

Table 3 The comparison of nominal wake fraction in different depth

工況h/T平均速度/(m·s-1)平均伴流分數121.31.59010.276215.151.58280.2793101.58360.279441.56950.285531.53560.301621.42270.352

從表3中可以看到槳盤面平均伴流分數的變化規律與前面分析過的阻力系數分布規律、槳盤面伴流場分布規律相吻合。隨著水深h的減小,伴流分數逐漸增大,軸向平均伴流分數增幅最高達27.2%。受水深減小的影響,船艉流速增加,艉傾現象明顯,船艉伴流相應增大。

2.4 船周圍流場變化分析

隨著水深的減小,船底速度逐漸增大,其增加的速度稱為回流速度。圖10表示KCS船在淺水中航行時船底流速的變化情況。隨著水深的減小,船體底部與底面邊界之間形成狹窄的水道,導致船底的流速增大;另外由于水的粘性,在船底與底面邊界形成邊界層,使過水斷面更小,導致船底的流速增大,壓力降低,從而使船體下沉,吃水增加,進而導致阻力增大。

圖10 不同水深處船中縱剖面速度分布Fig.10 The distribution of velocity in different depth

另外,從圖中可以看出,船底與底面邊界層厚度自船艏向船艉是逐漸增加的,導致相同水深時船艉部過水斷面較船艏更小,速度更大,壓力下降更甚,故船艉下沉較船艏大,產生艉傾現象。同時,隨著水深的減小,水流與船體的相對速度增大,壓力下降亦大,故艏艉壓力差將增大,艉傾現象更加明顯。且船艉與池底的間隙小,易于發生渦流,因此渦流阻力也要增大。

圖11給出了不同吃水時,船體周圍的流線分布。圖中可以看出,h=2T時,船側流線分布與深水時相比更密,由于淺水航道限制,船體排開的水主要從船體兩側通過;對于船底部流線分布,淺水航行時要比深水航行時更加不均勻,分析認為淺水航行時,船底與底面邊界形成的狹窄水道以及底面邊界的阻礙作用造成了流線分布的不均勻。

圖11 不同水深處船體周圍流線分布Fig.11 The distribution of streamline in different depth

3 結論

本文對KCS船模在不同水深下的船舶阻力及流場特性進行模擬,分析了水深變化對船舶阻力性能及流場特性的影響。分別計算了水深為h/T=21.3、h/T=15.15、h/T=10、h/T=4、h/T=3以及h/T=2六種工況,對比分析了不同工況船體的阻力系數、抽吸力、槳盤面伴流場以及邊界層速度場的特點,主要得出以下結論:

1) 當水深h小于10T時,淺水效應明顯,船舶各項性能隨水深的減小變化較為明顯,當水深h大于10T時,船舶各項性能受水深影響很小,不存在淺水效應;

2) 當水深較淺時,池底會對船體附近流場產生較大的干擾,隨著水深的減小阻力系數不斷增大,當水深超過10倍吃水時,阻力隨水深變化不大;

3) 當水深較淺時,隨著水深的減小,抽吸力不斷增加,并且增加的幅度越來越大,當水深超過10倍吃水時,z方向抽吸力基本穩定,隨著水深的減小,船艏艉壓力差增大,艉傾現象更加明顯;

4) 水深減小的過程中,槳盤面各方向伴流分數呈增大趨勢;

5) 船在淺水中航行時,隨著吃水的減小,船底流速增加,其摩擦阻力及渦流阻力增加。

本文計算過程中,沒有考慮自由葉面及船舶六自由度運動的影響,在今后計算中需進一步探討隨著水深的變化,船周圍興波及艉傾對船舶運動性能的影響。

[1]MüLLER E. Analysis of the potential flow field and of ship resistance in water of finite depth [J]. I S P, 1985(32):266-277.

[2]SAHIN I, CRANE J W, WATSON K P. Application of a singularity panel method for hydrodynamics of underwater vehicles[J]. Ocean engineering, 1997, 24(6): 501-512.

[3]SAHIN I, HYMAN M C. Simulation of three-dimensional finite-depth wave phenomenon for moving pressure distributions [J]. Ocean engineering, 2001, 28(12): 1621-1630.

[4]郭圣漢,劉希武. 船舶的淺水效應及隨車假底淺水試驗技術[J]. 橋梁建設, 1990(2): 53-56. GUO Shenghan, LIU Xiwu. Shallow water effect of ship and towing testing method in shallow water[J]. Bridge construction, 1990(2): 53-56.

[5]李一兵. 內河船舶航行阻力計算方法討論[J]. 水道港口, 2002(1): 7-11. LI Yibing. Calculation method of riverboat navigation resistance[J]. Journal of waterway and harbor, 2002(1): 7-11.

[6]繆濤, 張志宏, 顧建農, 等. 面元法求解有限水深船舶興波及水底壓力變化[J]. 計算力學學報, 2012,29(3): 464-469. MIAO Tao, ZHANG Zhihong, GU Jiangnong,et al. The calculation of shio wave and bottom pressure variation in finite depth by panel method[J]. Chiness journal of computational mechanics, 2012, 29(3): 464-469.

[7]WANG Jinbao, YU Hai, ZHANG Yuefeng,et al.CFD-based method of determining form factorkfor different ship types and different drafts[J]. Journal of marine science and application, 2016, 15(3):236-242.

[8]WU Chengsheng, ZHOU Decai, GAO Lei, et al. CFD computation of ship motions and added resistance for a high speed trimaran in regular head waves [J]. International journal of naval architecture and ocean engineering, 2011, 3:105-110.

[9]SENTHIL P M N, BINOD C. Numerical estimation of shallow water resistance of a river-sea ship using CFD [J]. International journal of computer applications, 2013, 71(5): 33-40.

[10]PHILIPP M, OULDEL M. Numerical prediction of resistance and squat for a containership in shallow water [C]∥17th Numerical Towing Tank Symposium, 2014: 352-359.

[11]王福軍. 計算流體動力學分析-CFD 軟件原理與應用[M].北京: 清華大學出版社, 2004: 124-126. WANG Fujun.The analysis of computation fluid dynamics-the based theory and application of CFD[M].Beijing: Tsing-hua University Press, 2004: 124-126.

[12]HE Miao, WANG Chao, CHANG Xin, et al. Numerical analysis of propeller wake flowfield using hybrid mesh technology [J]. Journal of marine science and application, 2012, 11(3): 295-300.

[13]王超,何苗,郭春雨,等. 一體化節能推進裝置定常水動力性能的數值模擬分析[J]. 哈爾濱工程大學學報, 2013, 34(6): 674-679. WANG Chao,HE Miao,GUO Chunyu,et al.Steady hydrodynamic performance simulation of integrative energy-saving propulsion system[J].Journal of Harbin Engineering University, 2013, 34(6): 674-679.

Analysis of ship resistance and flow field characteristics in shallow water

SUN Shuai, WANG Chao, CHANG Xin, ZHI Yuchang

(College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China)

To analyze the characteristics of ship resistance and flow field in shallow water, we adopted a hybrid grid combined with the Reynolds-averaged Navier-Stokes (RANS) method and used the KRISO container ship (KCS) model in shallow water to perform numerical prediction and analysis. We adopted three unstructured grid methods to determine the most appropriate stern mesh generation method. By changing water depth, we analyzed the effect of water depth on the hydrodynamic performance of the ship. The results show that the shallow-water effect is distinct when the water depthhis less than 10 times the draft. With a decrease in water depth, the ship resistance coefficient, the suction force in the z direction, and the propeller flow fraction show increasing trends. With the decrease of water depth, the pressure difference between the head and the tail increases gradually, and the tail inclination is more obvious.

hybrid grid; KRISO container ship; shallow-water effect; ship resistance; suction force; wake field; Reynolds-averaged Navier-Stokes

2015-12-07.

日期:2017-03-18.

國家自然科學基金項目(51379040, 51679052).

孫帥(1989-), 男, 博士研究生; 王超(1982-), 男, 副教授.

孫帥,E-mail:sunshuaidoc@163.com.

10.11990/jheu.201512026

U661.1

A

1006-7043(2017)04-0499-07

孫帥, 王超, 常欣,等.淺水效應對船舶阻力及流場特性的影響分析[J]. 哈爾濱工程大學學報, 2017, 38(4): 499-505.

SUN Shuai, WANG Chao, CHANG Xin,et al. Analysis of ship resistance and flow field characteristics in shallow water[J]. Journal of Harbin Engineering University, 2017, 38(4): 499-505.

網絡出版地址:http://kns.cnki.net/kcms/detail/23.1390.u.20170318.0716.012.html

猜你喜歡
船舶
船舶避碰路徑模糊控制系統
計算流體力學在船舶操縱運動仿真中的應用
CM節點控制在船舶上的應用
基于改進譜分析法的船舶疲勞強度直接計算
《船舶》2022 年度征訂啟事
船舶(2021年4期)2021-09-07 17:32:22
船舶!請加速
BOG壓縮機在小型LNG船舶上的應用
船舶 揚帆奮起
軍工文化(2017年12期)2017-07-17 06:08:06
船舶壓載水管理系統
中國船檢(2017年3期)2017-05-18 11:33:09
小型船舶艉軸架設計
船海工程(2015年4期)2016-01-05 15:53:30
主站蜘蛛池模板: 98超碰在线观看| 97国产成人无码精品久久久| 中国精品自拍| 国产迷奸在线看| 午夜性刺激在线观看免费| 人妻中文久热无码丝袜| 国产成人乱无码视频| 红杏AV在线无码| 久久熟女AV| 四虎国产永久在线观看| 91口爆吞精国产对白第三集| AV熟女乱| 青青操国产| 久操中文在线| 国产精品lululu在线观看| 免费不卡视频| 无码一区中文字幕| www亚洲精品| 国产精品片在线观看手机版| 亚洲天堂福利视频| 中国一级特黄视频| 91成人在线观看| 成人国产精品一级毛片天堂| 高潮爽到爆的喷水女主播视频| 国产免费人成视频网| 无码在线激情片| 伊人色综合久久天天| 精品国产aⅴ一区二区三区| 亚洲成人免费看| 一本大道无码日韩精品影视 | 思思99思思久久最新精品| 九一九色国产| 国产在线视频导航| 无码专区第一页| 亚洲区欧美区| 88国产经典欧美一区二区三区| 亚洲经典在线中文字幕| 五月婷婷丁香综合| 99伊人精品| 人人妻人人澡人人爽欧美一区| 71pao成人国产永久免费视频| 久久综合伊人 六十路| 99久久人妻精品免费二区| 国产h视频免费观看| 黄色在线不卡| 一本视频精品中文字幕| 欧美性天天| 免费看的一级毛片| 一区二区三区成人| 免费亚洲成人| 视频二区国产精品职场同事| 欧美成人一区午夜福利在线| 亚洲无限乱码| 欧美伦理一区| 午夜少妇精品视频小电影| 久久频这里精品99香蕉久网址| 91娇喘视频| 亚洲综合狠狠| 国产亚洲精品无码专| 亚洲国产精品一区二区第一页免| 99精品免费欧美成人小视频| 亚洲视频色图| 亚洲国产精品日韩av专区| 欧美成人怡春院在线激情| 久久香蕉国产线| 久久国产成人精品国产成人亚洲| 综合社区亚洲熟妇p| 五月天久久综合| 国产精品手机在线观看你懂的| 久爱午夜精品免费视频| 狠狠亚洲五月天| 午夜精品福利影院| 自拍中文字幕| 久久精品视频亚洲| 黄色免费在线网址| 欧美精品一区在线看| 亚洲视频黄| 不卡无码网| 国产精品白浆无码流出在线看| 午夜精品国产自在| 露脸真实国语乱在线观看| 波多野结衣一二三|