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

實船自航試驗數值模擬及尺度效應分析

2016-10-11 07:59:01李亮王超孫帥孫盛夏
哈爾濱工程大學學報 2016年7期

李亮,王超,孫帥,孫盛夏

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

實船自航試驗數值模擬及尺度效應分析

李亮,王超,孫帥,孫盛夏

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

以KCS船和KP505槳為計算對象,采用RANS方法和VOF模型,開展了考慮自由液面的實船自航試驗數值模擬。首先進行KCS船模自航點工況下的流場計算和KP505槳的敞水性能計算,計算結果與試驗值吻合較好,驗證了計算方法的可行性。本文分析了船體阻力、波形和速度場的尺度效應,獲得了實船的推進因子。根據數值自航試驗曲線確定了自航點,進而插值計算得到實船推進因子,并分析發現自航點轉速和伴流分數尺度效應明顯。結果表明:實船伴流分數要小于船模,自航點轉速要大于船模。

實船試驗;數值模擬;自航點;自由液面;尺度效應;伴流

船模自航試驗是預估實船性能和判斷船-機-槳匹配好壞的重要手段,對于新設計的船舶來說,還可以對若干方案進行比較,從而選型擇優,船模自航試驗的具體操作規程和方法已由ITTC(International Towing Conference)[1]給出。近年來,隨著計算機技術的高速發展,運用CFD方法已經能成功模擬船模自航試驗[2-4],其精度足以滿足工程計算的需要,相比水池試驗能節省大量的時間和成本。但由于自航試驗船模與實船之間只滿足傅汝德數和進速系數相等,雷諾數并不相等,導致模型數據換算到實船的過程中因為尺度效應的存在而產生較大的誤差,盡管有相關的經驗統計方法(ΔCT、Δω法和1978ITTC標準方法等)用來修正,但對不同類型的船舶其可靠性也是堪憂的。目前,船槳一體模型尺度的流場和水動力數值計算方法已經較為成熟[5-7],在此基礎上已經開始有學者嘗試實尺度船體的CFD計算研究,并取得了一定成果。Alejandro M.Castro等[8]基于動態重疊網格方法,開展了實尺度的KCS數值模擬自航試驗研究,對比了船模自航和實船自航推進因子和流場等各方面的不同;熊鷹等[9]提出了不考慮船體興波的自航船模推進因子計算方法,對實船推進因子做了數值預報研究工作。由于實船自航數值模擬網格數量太大、邊界層網格厚度難以保證、收斂時間長等問題,使得氣-液兩相粘性流計算較船模自航更為復雜和困難,自由液面也不易捕捉,國內在這方面的研究工作還很少,但是開展實尺度船舶自航試驗的數值模擬工作對尺度效應的修正方法研究是極具意義的,也是未來CFD方法發展的必然趨勢。

本文基于RANS方法,采用VOF方法捕捉自由液面,開展了KCS標模實船自航試驗數值模擬。首先進行模型尺度的自航狀態計算,通過與水池試驗對比驗證了該方法的可靠性。由于直接模擬實船自航計算過于復雜,數值模擬分為兩階段進行,先保證螺旋槳靜止船體以服務航速航行,待流場收斂穩定后再以強制自航法的方式預報實船自航點,并比較了本次計算結果與試驗換算值[10]和文獻[8-9]中計算值之間的差異,還分析了模型尺度和實船之間由于尺度效應導致的阻力、流場、波形及推進因子變化。

1 數值方法

1.1控制方程和湍流模型

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

1.2VOF模型

VOF(volume of fluid)方法的基本原理是通過研究網格單元中流體和網格體積比函數來確定自由面,追蹤流體的變化,而非追蹤自由液面上質點的運動。只要知道函數在流場中每個網格上的值,就可以實現對運動界面的追蹤。

將整個計算區域定義為Ω,主相流體區域記為Ω1,副相流體區域記為Ω2。VOF定義這樣一個函數:

此外,在由兩種互不相溶流體構成的流場中,流體的速度場記為V=(u,v),函數ω滿足:

在每個網格Iij上定義ω(x,t)在網格上的積分為Cij,可以得到VOF函數:

VOF函數也滿足式(2):

顯然,當C=0時,網格中全為副相流體;當C= 1的時,網格充滿主相流體;當0<C<1時,則是含有流體界面的網格,成為界面網格。

2 計算模型的建立

2.1計算對象和工況

本文研究對象為KCS集裝箱船,與之搭配的槳為KP505槳,二者模型如圖1所示。實船具體參數如表1,本文所有計算工況如表2所示。

圖1 KCS船槳模型Fig.1 Model of KCS ship and propeller

表1 KCS船和KP505槳主要參數Table 1 Main parameters of KCS ship and KP505 propeller

表2 工況介紹Table 2 Working conditions introduction

2.2網格劃分

本次計算采用混合網格方法,全部網格劃分工作在ICEM軟件中完成。由于KCS船艉收縮曲率大,劃分結構化網格困難,所以將船艉小部分流域劃分非結構化網格,以節省網格劃分工作量,而其它部分包括螺旋槳在內都采用六面體結構化網格,結構化網格和非結構化網格之間通過interface連接,流場通過interface插值進行信息傳遞。

數值計算時模型尺度和實尺度網格拓撲結構一致,船體表面網格如圖2(a)所示,區別在于首層網格厚度和三個方向上的網格節點數目。一般來說,模型尺度Y+建議控制在60左右是合適的,由式(5)可計算出第一層網格厚度為0.8~1 mm:

式中:L為特征長度。實尺度Y+控制范圍可參考的文獻不多,本文根據大量計算統計和文獻[13]建議,發現Y+值在300左右是合理的,對應第一層網格厚度為1~2 mm。結構化網格區域第一層層網格厚度可以直接用參數定義,而船艉非結構化網格則需要通過添加多層棱柱網格完成,如圖2(b)所示。螺旋槳網格如圖2(c)所示。另外值得注意的是,計算劃分網格時,船體首尾形狀復雜,流場變化劇烈,應給與適當加密;為了較好地捕捉自由液面,自由液面附近也要增加網格節點數目以提高網格分辨率。

圖2 KCS自航數值計算網格Fig.2 Computational grid for numerical self-propulsion KCS ship

2.3邊界條件設置

流場計算區域為:入口距船艏1.0LPP,出口距船艉2.0LPP,側面和底面均距船體表面1.0LPP。入口劃分為空氣速度入口和水流速度入口,二者速度大小一致,模型條件下Vm=2.196 m/s,實尺度條件下VS=12.346 m/s;出口利用Fluent用戶自定義函數(UDF)設置為壓力出口,出口垂直方向壓強按下式變化:

式中:p0為一個大氣壓,ρ為水的密度,g為重力加速度,z為垂直方向坐標值,T為吃水。流域上邊界定義為對稱面,船體表面、側面和底面均定義為無滑移壁面,螺旋槳旋轉運動采用MRF模型,壓力速度耦合迭代采用SIMPLEC方法。

3 計算結果分析

3.1船模自航點工況數值計算

為了驗證本文計算方法的準確性,首先開展船模尺度的推進因子數值預報。參考水池自航試驗流程,要獲得船模推進因子,需要分別進行KCS船模裸船體阻力計算、KP505槳模敞水性能計算和船槳一體自航點工況數值計算,自航點時螺旋槳轉速為570 r/min。為了節省計算工作量,僅在KP505槳最高效率點附近取了四個不同進速系數J=0.6,0.7,0.8,0.9分別進行敞水性能計算。

圖3 自由液面波形圖對比Fig.3 Comparison of free surface wave pattern

圖3所示為船模尺度裸船體計算自由液面波形與試驗值對比圖,從圖中可以看出二者各處波形和波高基本一致,自由液面總體捕捉效果較好,僅在尾部其波形細節存在一定差異,這可能是因為遠離船體流域網格尺度逐漸增大,波能耗散較快。

圖4所示為螺旋槳敞水性能曲線,可以看出在最大效率點附近計算值與試驗值非常接近,最大誤差不超過3%,計算精度滿足要求。

圖4 螺旋槳敞水性能Fig.4 The open water performance of propeller

表3為船模自航點計算結果與試驗值對比,表中Ct表示裸船體阻力系數,CtSP表示自航狀態下船體阻力系數,KT表示船后螺旋槳推力系數,KQ船后螺旋槳轉矩系數,J表示據 KT用等推力法查敞水曲線所得進速系數,wm表示伴流分數,tm表示推力減額。通過比較可知,除了推力系數和推力減額外,各項計算結果均與試驗值吻合良好,誤差在3%以內,推力誤差較大有可能是因為螺旋槳旋轉域與非結構網格流域相接的這對interface網格節點數目存在一定差異,從而增大了流場數據插值傳遞誤差。

表3 船模自航點計算結果與試驗值對比Table 3 Comparison of calculation resut and experimental value at model-scale self-propulsion point

3.2實船自航試驗數值計算

3.2.1實船船體阻力數值計算

在確定實船自航點之前首先要進行的是服務航速下帶自由液面的裸船體阻力計算,其計算結果如表4所示。表中試驗值摩擦阻力修正系數SFC可由式(7)計算得到,試驗總阻力系數可由式(8)計算得到

式中:形狀因子取1+k=1.1,摩擦阻力系數根據1957ITTC公式計算得到CFM=2.832×10-3,CFS= 1.378×10-3,粗糙度補貼系數由文獻[14]提供ΔCF=0.27×10-3。

表4 實尺度KCS裸船體阻力系數計算結果Table 4 Results of resistance coefficients for full-scale KCS bare hull

從表4中可以看出,這與文獻[8]計算結果誤差趨勢一致,實船總阻力系數計算值要比ITTC換算值要略大,誤差為3.83%,摩擦阻力修正系數SFC數值計算值要比計算值(式(6))要小,誤差為2.63%。文獻[13]中的計算結果顯示形狀因子1+k存在較大的尺度效應,隨著雷諾數的增大,尺度的增加,形狀因子將增大。所以如果按照船模尺度的形狀因子進行實船阻力換算,而不考慮其尺度效應,必然導致摩擦阻力修正系數過小,從而換算得到的實船總阻力系數偏大。但是從表中數據來看,試驗換算值與本文計算值和文獻[8]計算值均吻合良好,其主要原因是船模試驗確定的形狀因子是偏大的,文獻[13]中KCS模型尺度形狀因子數值計算值僅為1.05。船模形狀因子確定的主要方法有普魯哈斯卡(Prohaska)方法和第15屆 ITTC推薦方法,在 Fr為0.1 ~0.2范圍內測量得到,但是這都是以休斯假設為基礎的,即認為形狀因子1+k是與船體形狀有關的常數,實際上諸多試驗結果都表明形狀因子1+k在低速時近似為常數,而在較高航速(Fr>0.16)后,隨著Fr增大而減小,本次計算中Fr=0.26,故Fr增大導致的形狀因子誤差很大程度上抵消了尺度效應所帶來的誤差,從而使得試驗換算值和實船數值計算值非常靠近,這也表明現有的三因次換算方法形狀因子的確定方法是可取的,大大降低了尺度效應。

圖5所示為實船與船模縱剖面y/Lpp=0.006處無量綱軸向速度等值線圖,從圖中可以明顯看到實船的邊界層要比船模的更薄,船模的低速區要向下游拖延得更遠,在相同X軸位置處,實船船艉下游的速度均大于船模。同時可以看到一個有意思的現象,就是實船時船艉還形成了一個船模沒有的速度閉合區,說明此縱剖面處有一定的回流出現。另外觀察艉封板附近流域,發現實船船艉還有一股相對船模速度更大,波峰更高的急流涌現。

圖6為實船與船模槳盤面處無量綱軸向速度等值線圖,圖中也體現出了實船邊界層要比船模更薄的特點,實船的軸向速度等值線相比船模要向里收縮,這使得在槳盤面處實船具有更大的軸向速度,也就是說其軸向標稱伴流分數要小。觀察螺旋槳所在位置區域,可以發現實船與船模之間的伴流差異是很大的,實船船艉后的螺旋槳的進速要大于船模,如果根據船模測量伴流分數來設計螺旋槳,就很容易導致安裝在實船上的螺旋槳會產生推力不足的情況,所以根據伴流的尺度效應情況必須留出一定推力余量。

圖5 實船與船模縱剖面(y/Lpp=0.006)無量綱軸向速度對比Fig.5 Comparison of full-scale and model-scale nondimensional axial velocity at y/Lpp=0.006

圖6 實船與船模槳盤面(x/Lpp=0.982 5)無量綱軸向速度對比Fig.6 Comparison of full-scale and model-scale nondimensiona axial velocity at the propeller plane x/Lpp=0.982 5

圖7 實船與船模自由液面波形圖對比Fig.7 Comparison of full-scale and model-scale free surface wave pattern

圖7所示為實船與船模波高等值線圖,船艏和船舯部分的興波二者幾乎看不出差別,這和興波阻力與雷諾數無關的假設是相符合的,但是在尾部及其下流區域波形存在一定的差異,主要是雷諾數增大導致的邊界層差異增大了實船船艉流速,進而使得實船船艉急流波峰增高,其尾部波形略微有整體向后移動,這與圖5中觀察到的現象一致。

3.2.2實船自航點確定步驟

實船自航試驗的數值模擬方法和模型基本一致,唯一不同的是確定自航點的過程中不用再考慮因為雷諾數Re的不同而要進行摩擦阻力的修正,當船體阻力和螺旋槳推力達到平衡時即可確定自航點下的螺旋槳轉速、推力和轉矩,這有效避免了阻力換算和伴流尺度效應等因素帶來的誤差。確定實船自航點的具體操作步驟如下:

1)根據船模自航試驗結果,預估實船一個實船自航點。為了方便,可以將不做修正的船模自航點根據縮尺比直接換算過來作為實船的預估自航點,本文計算預估的自航點螺旋槳轉速為:

2)在預估的自航點螺旋槳轉速Ns0前后適當范圍內再各取兩個轉速,本文取Ns1=89 r/min,Ns2= 95 r/min,Ns3=107 r/min,Ns4=113 r/min,然后在服務航速VS=12.346 m/s下,分別對五個不同螺旋槳轉速進行數值模擬;

3)繪制船體阻力和螺旋槳推力隨螺旋槳轉速的變化曲線,兩曲線的交點即為該航速下的自航點。

3.2.3實船自航點數值計算結果

依據3.2.2節中實船自航點確定步驟,數值模擬得到五個不同螺旋槳轉速下螺旋槳推力、轉矩和船體阻力,其結果如表5所示,實船自航試驗曲線繪制于圖8中。由Rt=T,在實船自航曲線上通過插值可得自航點N=107.3 r/min,Rt=T=1.995× 106N,Q=2.589×106,103Ct=2.783,KT=0.161,10KQ=0.264。

表5 不同轉速實船自航數值模擬結果Table 5 Results of numerical simulations of full-scale self-propulsion at different rotating speed

表6 實船自航推進性能預報結果Table 6 Prediction results of the propulsion performance of full-scale ship

為節省計算量,本文沒有進行實尺度螺旋槳敞水性能計算,主要是因為推力受尺度作用的影響很小,幾乎可以忽略,轉矩系數變化一般也在1%左右,文獻[8]中的計算值也驗證了這一點,實槳和槳模的敞水曲線基本重合。基于以上分析,實船自航點螺旋槳的進速系數可直接利用等推力法在模型槳的敞水特性曲線上插值得到,J0=0.755,10KQ=0. 275。

圖8 實船自航試驗曲線Fig.8 Results of full-scale self-propulsion test

數值模擬最后得到的實船自航推進性能結果具體如表6所示,表中數據顯示實船計算值和船模試驗換算值之間差異較大的量為自航點轉速和伴流分數,其誤差分別為+5.82%和+8.96%。3.2.1節中,通過分析圖4和圖5可知,由于實船雷諾數與船模雷諾數相差巨大,其邊界層厚度也要更薄,實船槳盤面的進速是要大于船模的,故實船自航試驗模擬得到實船的伴流分數要小于船模,自航點轉速要大于船模是預料之中的結果。這啟示在做模型試驗時僅考慮螺旋槳載荷的尺度效應是不夠的,邊界層厚度和伴流分數的尺度效應也是必須引起注意的。

4 結論

本文選用KCS船和KP505槳為計算對象,考慮自由液面的情況下開展了實船自航試驗數值模擬和尺度效應研究分析,根據數值模擬結果主要得出了以下結論:

1)在模型自航點工況下,自由液面波形、船體阻力和螺旋槳推力轉矩值和試驗吻合良好,表明計算精度滿足要求,本文計算方法準確可靠。

2)對比實船與船模裸船數值計算結果,發現由于雷諾數的巨大差異,實船邊界層要比船模薄,槳盤面位置進速要更大,如果依據船模伴流分數來進行螺旋槳設計容易導致實船推力不足。同時實船艉流高速區向后延伸,并有一股相對船模來說速度更大,波峰更高的急流涌現,這使得實船船艉自由液面波形也略微有所整體后移。

3)運用強制自航法數值預報得到實船自航點,插值計算得到推進因子,與試驗換算值對比,誤差均在可接受范圍之內,尺度效應主要體現在伴流分數和自航點轉速,結果顯示實船伴流分數要小于船模,自航點轉速要大于船模。

4)實船自航試驗數值模擬是數值計算實船螺旋槳誘導激振力的一個必要工作,下面將以此為基礎進一步開展實船激振力的相關研究。

[1]ITTC.Testing and extrapolation methods,Propulsion,performance propulsion test[R].Recommended Procedures and Guidelines,Report 7.5-02-03-01.1.Rio de Janeiro:ITTC Secretary,2002.

[2]CARRICA P M,CASTRO A M,STERN F.Self-propulsion computations using a speed controller and a discretized propeller with dynamic overset grids[J].Journal of marine science and technology,2010,15(4):316-330.

[3]戈亮,顧民,吳乘勝,等.水面船自航因子CFD預報研究[J].船舶力學,2012,16(7):767-773. GE Liang,GU Min,WU Chengsheng,et al.CFD prediction of self-propulsion parameters for a surface ship[J].Journal of ship mechanics,2012,16(7):767-773.

[4]NOBUAKI S,YASUTAKA K,SHOTARO U,et al.Estimation of resistance and self-propulsion characteristics for low L/B twin-skeg container ship by a high-fidelity RANS solver [J].Journal of ship research,2013,57(1):24-41.

[5]ABDEL-MAKSOUD M,MENTER F,WUTTKE H.Numerical computation of the viscous flow around the series 60 CB= 0.6 ship with rotating propeller[C]//Proceedings of the 3rd Osaka colloquium on advanced CFD applications to ship flow and hull form design.Osaka,1998.

[6]劉祥珺,孫存樓.數值水池船模自航試驗方法研究[J].艦船科學技術,2011,33(2):28-31. LIU Xiangjun,SUN Cunlou.Research on ship self-propulsion model test in numerical tank[J].Ship science and technology,2011,33(2):28-31.

[7]CARRICA P M,FU Huiping,STERN F.Computations of self-propulsion free to sink and trim and of motions in head waves of the KRISO container ship(KCS)model[J].Applied ocean research,2011,33(4):309-320.

[8]CASTRO A M,CARRICA P M,STERN F.Full scale selfpropulsion computations using discretized propeller for the KRISO container ship KCS[J].Computers&fluids,2011,51(1):35-47.

[9]熊鷹,劉志華.自航船模和實船推進因子的數值預報方法研究[J].船舶力學,2013,17(1/2):14-18. XIONG Ying,LIU Zhihua.Numerical prediction of propulsion factors of propelled ship model and full scale ship[J]. Journal of ship mechanics,2013,17(1/2):14-18.

[10]LARSSON L,STERN F,VISONNEAU M.CFD in ship hydrodynamics—results of the Gothenburg 2010 workshop [C]//MARINE 2011,IV International Conference on Computational Methods in Marine Engineering.Netherlands:Springer,2013:237-259.

[11]王福軍.計算流體動力學分析——CFD軟件原理與應用[M].北京:清華大學出版社,2004:7-11. WANG Fujun.Principle and application of CFD software [M].Beijing:Tsinghua University Press,2004:7-11.

[12]MENTER F R.Two-equation eddy-viscosity turbulence models for engineering applications[J].AIAA journal,1994,32(8):1598-1605.

[13]王展智.采用混合式CRP推進器的船舶水動力特性及尺度效應研究[D].武漢:海軍工程大學,2014. WANG Zhanzhi.Study on the hydrodynamic characteristic and scale effect of a ship equipped with the hybrid CRP pod propulsion system[D].Wuhan:Naval University of Engineering,2014.

[14]National Maritime Research Institute.Tokyo 2015 a workshop on CFD in ship hydrodynamics[EB/OL].(2016-04-15).http://www.t2015.nmri.go.jp/program.html.

本文引用格式:

李亮,王超,孫帥,等.實船自航試驗數值模擬及尺度效應分析[J].哈爾濱工程大學學報,2016,37(7):901-907.

LI Liang,WANG Chao,SUN Shuai,et al.Numerical simulation and scale effect of self-propulsion test of a full-scale ship[J].Journal of Harbin Engineering University,2016,37(7):901-907.

Numerical simulation and scale effect of self-propulsion test of a full-scale ship

LI Liang,WANG Chao,SUN Shuai,SUN Shengxia
(College of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China)

In this study,we adopt the Reynolds-averaged Navier-Stokes(RANS)method and volume of fluid (VOF)model to carry out a numerical simulation of a self-propulsion test of a full-scale ship,considering a free surface on a KCS ship and KP505 propeller.First,we performed a numerical computation of a model-scale KCS under the conditions of a self-propulsion point and open-water propeller performance.The calculated results were in good agreement with the experimental data,thereby verifying the feasibility of the calculation method.Next,we analyzed the scale effect on the ship resistance,wave pattern,and velocity field,and obtained the propulsion factors of a full-scale ship.Based on the curves of the numerical self-propulsion test,we confirmed the self-propulsion point of the full-scale ship.Using the interpolation method,we identified the propulsion factors and found the scale effect of the rotational speed at the self-propulsion point and wake fraction to be obvious.The results show that the wake fraction of a full-scale ship is smaller than that of a model ship and its rotational speed at the self-propulsion point is greater.

full-scale experiments;numerical simulation;self-propulsion point;free surface;scale effect;side flow

10.11990/jheu.201507015

U661.3

A

1006-7043(2016)07-901-07

2015-07-05.網絡出版日期:2016-05-13.

國家自然科學基金項目(51309061);中央高校基本科研業務費專項資金項目(HEUCFD1515).

李亮(1990-),男,碩士研究生;王超(1981-),男,副教授.

王超,E-mail:wangchao0104@hrbeu.edu.cn.

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

主站蜘蛛池模板: 在线亚洲精品福利网址导航| 91久久大香线蕉| 香蕉视频在线精品| av手机版在线播放| 亚洲欧洲自拍拍偷午夜色无码| 欧美成人精品一级在线观看| 亚洲国产看片基地久久1024 | 国精品91人妻无码一区二区三区| 欧美中文字幕在线二区| 午夜激情福利视频| 青草精品视频| 久久亚洲国产视频| 久久综合伊人 六十路| 91精品国产福利| 人妻一区二区三区无码精品一区| 国产原创第一页在线观看| 美女高潮全身流白浆福利区| 亚洲乱伦视频| YW尤物AV无码国产在线观看| 国产精品无码久久久久久| 国产一级小视频| 国产美女自慰在线观看| 日韩av手机在线| 精品国产中文一级毛片在线看| 婷婷色婷婷| 婷婷午夜天| 国产在线一区视频| 精品自窥自偷在线看| 色偷偷综合网| 99这里只有精品在线| 国产精品99久久久| 日韩黄色精品| 国产高潮视频在线观看| 欧美日韩成人在线观看| 免费a级毛片视频| 国产精品七七在线播放| 亚洲国产欧美国产综合久久| 欧美成人精品一级在线观看| 国产在线观看第二页| 97成人在线观看| 日韩午夜伦| 国内精品自在欧美一区| 91九色最新地址| h网址在线观看| 久草网视频在线| 伊人久热这里只有精品视频99| 91破解版在线亚洲| 99久久人妻精品免费二区| 国产美女在线观看| 色婷婷综合激情视频免费看| 婷婷开心中文字幕| 2020最新国产精品视频| 无码福利日韩神码福利片| 日本久久久久久免费网络| 国产黄网站在线观看| 亚洲国产精品日韩欧美一区| 婷婷激情亚洲| 免费无遮挡AV| 在线观看亚洲天堂| 国产精品久久久久久久久久98| 精品久久高清| 国产色伊人| 日本黄色a视频| 新SSS无码手机在线观看| 91精品国产麻豆国产自产在线| 国产女人18水真多毛片18精品| 国产小视频a在线观看| 一区二区三区精品视频在线观看| 久久久久人妻一区精品| 国产精品私拍在线爆乳| 中文国产成人久久精品小说| 爱色欧美亚洲综合图区| 91探花在线观看国产最新| 91成人精品视频| 国产另类视频| 亚洲an第二区国产精品| 国产精品亚洲一区二区三区z | 久久国产精品娇妻素人| 青草91视频免费观看| www.狠狠| 亚洲狠狠婷婷综合久久久久| 亚洲A∨无码精品午夜在线观看|