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

高速三體船的阻力預報方法研究

2016-05-04 05:53:55周廣利艾子濤豆鵬飛黃德波
船舶力學 2016年7期
關鍵詞:船舶模型

周廣利,艾子濤,鄧 銳,豆鵬飛,黃德波

(1.哈爾濱工程大學 船舶工程學院,哈爾濱 1500012;2.英輝南方造船廠(廣州番禹)有限公司 技術中心,廣州 511431)

高速三體船的阻力預報方法研究

周廣利1,艾子濤2,鄧 銳1,豆鵬飛1,黃德波1

(1.哈爾濱工程大學 船舶工程學院,哈爾濱 1500012;2.英輝南方造船廠(廣州番禹)有限公司 技術中心,廣州 511431)

高速航行的三體船,其航行姿態會隨著傅汝德數(Fn)的增加而發生變化。常見的CFD計算預報實船阻力常常以設計浮態為基準、固定船舶航行姿態的方法進行,也是導致計算偏差的原因之一。文中選定一艘三體船模型,針對自由和約束兩種展開CFD數值模擬和阻力計算,得到其差別。為了進一步了解其興波情況和摩擦阻力的變化規律,在0.3<Fn<0.6的范圍內對處于自由和約束狀態下的模型流場相關細節進行了分析。最后,對該三體船型進行了模型試驗驗證。

三體船;航行姿態;CFD;阻力預報

0 引 言

三體船作為近年來比較受關注的船型,它具有一些單體船難以比擬的優點。寬闊的甲板面積有利于總體布置,瘦長的船型以及適當的側體位置使得高速航行時興波阻力較小,耐波性能良好[1]。優良的水動力性能使得三體船應用前景廣闊,國內外關于三體船興波阻力的理論計算、模型試驗以及數值模擬都取得了一定的進展,國內黃德波[2-5],王中[6-8],盧曉平[9]等均在多體船領域進行了相關研究。

三體船在高速航行時,船體受到的流體動壓力作用顯著增加,可以支托起部分船體重量,縱傾力矩也隨著改變,船體發生較明顯的升沉和縱傾,進而其航行姿態發生變化[10-11]。目前國內對于考慮航行姿態變化的船舶阻力預報研究不多,而這也是高速船型阻力預報需要解決的一個問題。通過實驗手段研究船舶航行姿態的方法主要有圖譜法與分部計算法。在圖譜法中會預先設定一組吃水和一組縱傾,交叉組合出各種吃水與縱傾狀態,進行不同速度下的阻力實驗,通過繪制曲線得到不同航行姿態下的水動力變化,這種方法工作量大,十分繁瑣,不宜操作;分部計算法是對這一問題的簡化,操作簡單但精度不高[12]。Joe和Fred[13]通過試驗手段對5512船體模型進行了計及航行姿態變化的阻力試驗,并取得了一定的結果。阿根廷學者Azcueta Repetto[14]對系列60船型和肥大型鈍艏船型在考慮航行姿態變化的情況下進行了阻力預報,對約束模型和自由模型分別進行了考察。

對于固定船舶航行姿態(約束)狀態的傳統數值模擬方法很難滿足高速船型的阻力預報精度,特別是Fn>0.3時,傳統的CFD預報結果與模型試驗結果會產生較大的偏差[15-19]。但隨著CFD商業軟件的發展,對模型在自由狀態下的數值模擬受到越來越多的關注,預報精度有了極大提高,將該方法運用于對高速三體船的阻力預報成為可能。

1 計及航行姿態的CFD模擬方法

1.1 VOF法

VOF模型[20]是基于兩種或多種流體(或相)是互不滲透的這一事實。對于引入模型中的每一相,引入一個稱為單元相體積分數的變量。基本原理是取坐標系o-xyz,x朝前,y朝左舷,z軸朝上,o-xy面位于靜水面,o-xz面為舯橫剖面。通過研究單元網格中流體與網格體積比的函數F,在有流體質點取值為1,無流體質點取值為0。這樣,在一個網格單元內,F的平均值就表示了單元內流體所占比例。假設流場中任意點(x,y),定義函數f(x,y,t):

1.2 船體運動

fr為定義成下式所示的延遲因子,表示定義的作用在模型上的力,為外部作用的力矩,對于六自由度模型,釋放時間ts和延遲時間tr組成的延遲函數定義如下:

延遲函數fr(t)應用在由于流體流動和重力產生的作用在模型上的力、力矩。釋放時間是指在計算開始前的一小段允許流體流動初始化的時間,其長短取決于模型網格的數量等因素。在釋放時間內,在流場啟動的瞬間,模型周圍的壓力場與速度場會發生劇烈變化,作用在模型上的力和力矩會使模型產生大幅振動,稱為瞬時效應,而應用于力和力矩的延遲時間是為了減緩這種瞬時效應。

2 CFD阻力預報

2.1 自由狀態下船模的阻力計算

自由狀態下模型的拖曳試驗無法展現流場細節,為了了解航行姿態變化對船體周圍流場細節的影響,數值模擬過程中釋放了沿z軸方向的運動與繞y軸方向的轉動,即升沉與縱傾兩個自由度,在0.3<Fn<0.6的范圍內對該三體船模型進行了數值模擬。船模的主尺度參數如表1所示,其中L為設計水線長,B為設計水線寬,D為設計吃水,Δ為排水量,S為浸濕面積。

表1 三體船模型主尺度參數Tab.1 marine particular of trimaran model

本文中基于六自由度求解RANS方程,湍流模型采用標準的k-ε模型,采用VOF方法捕捉自由液面,每個船模速度點計算時間為100 s。對于高速船的CFD數值模擬,船體及其附近的網格質量和網格密度對計算精度有很大的影響。為了避免流域內的網格扭曲程度過高,本文采用船行波加密區,向靠近船體的方向逐步加密。在離船體較近的區域,網格足夠密集,可以避免由于船體運動引起船體附近區域細小網格的扭曲;在離船體相對較遠的區域,網格則較為稀疏,是因為船體在運動過程中的幅度相對網格長度不是很大,這部分區域作為變形區域。流場采用隨體坐標,坐標原點位于船體中心處。船體周圍矩形區域為剛體網格,隨船體一起運動。最終生成網格約240萬,船體周圍網格如圖1~4所示。

圖1 計算初始時刻船體周圍橫向網格Fig.1 Transverse grids around the hull in the initial calculating moment

圖2 計算結束時刻船體周圍縱向網格Fig.2 Longitudinal grids around the hulls of the finished calculating moment

圖3 計算初始時刻船體周圍垂向網格Fig.3 Transverse grids around the center hull of the initial calculating moment

圖4 計算結束時刻船體周圍垂向網格Fig.4 Vertical grids around the hull of the finished calculating moment

數值模型可以沿z方向和繞y軸自由運動。在釋放這兩個自由度的瞬間,船體周圍的壓力場會發生急劇變化,船體表面的受力亦會發生較大的變化,隨著力矩變化,船體不斷調整航行姿態。最終,外力實現平衡,船體升沉及縱傾達到穩定。從圖1~4中可以看出在計算初始時刻與計算結束時刻船體周圍橫向與垂向網格基本沒有發生變形。

圖5 船體波形等高線圖(Fn=0.509)Fig.5 Contour lines of the ship(Fn=0.509)

從波形等高圖(圖5)中可以看出:當Fn=0.509時,靠近中體首部的興波較強,其幅度大于船體其他位置;而中體及側體尾部的波幅較小。這兩點與模型試驗時觀測到的宏觀現象十分相似。從速度云圖(圖6)上可以看到:當Fn=0.509時,中體及側體首部的邊界層很薄。主體球鼻處出現明顯的速度梯度且數值均大于來流速度,在設計吃水線以上出現首部射流,其速度均小于來流速度。側體首部處靠近中體一側與中體間興波干擾明顯,側體外側速度分布較為均勻。尾部中體與側體間的興波干擾較為明顯,但靠近中體處速度變化的位置略高于靠近側體速度變化的位置。對于側體外側,設計水線以下部分亦受到興波干擾。

圖6 流場速度分布云圖(Fn=0.509)Fig.6 Velocity profiles of the flow field(Fn=0.509)

2.2 與約束狀態下的對比分析

圖7 自由狀態下船體受力隨時間變化情況Fig.7 Change of the hull resistance with time in free condition

與約束狀態不同,自由狀態下的船體可以在z方向運動和繞y軸自由轉動,船體周圍的速度場和壓力場也會發生劇烈的變化,導致船體受力產生大幅度的波動。CFD通過若干迭代時間步內,加載的運動函數實時計算指定船體表面網格,整船受力可以通過船體表面各個網格節點的受力求和得出,不同航速下船體受力(沿y軸的力矩)隨時間變化情況如圖7所示。

由圖中可見,流場開啟瞬間,船體所受的繞y軸的力矩較大,使得船體在產生升沉的同時亦產生一定大小的縱傾力矩。船體吃水增加使得浮力再次與重力平衡,此時船體的縱傾角發生改變,船體受力達到新的平衡。迭代一定時間后,船體受力慢慢趨于穩定,船體升沉與縱傾也達到該航速下的動態平衡。圖8所示為通過VOF法捕捉自由液面,得到了船體縱傾與升沉過程隨時間的變化。

在低速范圍內,船舶的浮態不會發生較大變化,常規固定航行姿態的CFD模擬方法比較貼近船舶實際的航行狀況。但隨著航速的增加,船舶的浮態也隨之改變,其穩定航行時的姿態較初始的正浮狀態會發生較大的變化,這種姿態變化對船體興波及其濕表面積變化產生一定的影響。表2給出了船模濕表面積和摩擦阻力隨航速的變化情況,其中Δ f表示摩擦阻力的增量。

圖8 自由狀態下船模縱傾與升沉過程隨時間變化(Fn=0.566)Fig.8 Change of trim and heaving with time in free condition(Fn=0.566)

表2 濕表面積變化引起的摩擦阻力變化Tab.2 Change of frictional resistance with the wet-surface area

圖9 沿船體的浪高曲線(Fn=0.566)Fig.9 Wave height curve along the hull surface(Fn=0.566)

從表2可知,由于航行姿態改變,船舶實際的濕表面積均大于正浮狀態下的濕表面積且隨傅汝德數的增加而增加。因此,針對約束狀態下的船模CFD阻力計算結果,其摩擦阻力部分是偏小的。

通過對比約束和自由狀態下沿船體的浪高曲線可以發現,對于中體、側體內側和側體外側,兩者都具有十分相似的波形。自由狀態下船首的波高略高于約束狀態;對于船體后段,約束狀態下模型的興波波高高于自由狀態。分析發現,自由狀態下模型的船體首尾壓力差大于約束模,這也是自由狀態模型的興波阻力較大的原因,另外,自由與約束狀態下模型在水面上興起波浪的波峰位置與波高亦有差異,自由狀態模型的興波范圍大于約束狀態。

3 約束與自由狀態下的船模阻力試驗

本次試驗的目的主要是考察分析三體船在約束與自由狀態下的阻力差別,驗證數值模擬計算正確性。考慮到實船常用航速范圍,并為了在較大航速范圍內分析三體船阻力特點,確定模型試驗的速度與CFD計算點對應。該試驗在哈爾濱工程大學拖曳水池完成(圖10),阻力測量使用四自由度適航儀,測力傳感器為應變式,量程 100 N,精度0.1%。試驗中,船模處于約束狀態時,其縱、橫蕩,縱、橫搖,升沉以及搖首運動均受到限制;自由狀態時,其縱、橫蕩,橫搖及搖首運動受到限制,而縱搖與升沉運動則是自由的。試驗采用雙升沉桿拖帶方式,其中前升沉桿位于船模重心處,作為拖點,測力傳感器亦位于拖點處,后升沉桿則作為導向桿使用。在整個試驗過程中,船模的總排水量及壓載分布均保持不變。

圖11中給出了約束模型、自由模型CFD計算與實驗結果對比情況。由圖中可以看出,自由模型和約束模型CFD阻力預報結果分別與兩種狀態下的模型試驗結果吻合良好,自由模型與約束模型CFD預報結果與對應模型試驗結果誤差均在3%以內。無論是模型試驗還是CFD計算,自由模型的阻力值均大于約束模型,且偏差在5%~10%以內。分析其原因:一方面,約束模型較自由模型在自由表面的興波小,興波阻力被低估;另外,根據重力場中不可壓縮流體的伯努利方程V2/2+P/ρ+gz=c,當速度增大時,模型垂向位置z減小,自由模型吃水增加,導致其濕表面積增加,相對于浮態不發生變化的約束模型來說,按照自由模型進行摩擦阻力計算亦更為準確。

對于高速航行中的三體船,因船體運動而引起其周圍流場的變化,船體所受到的垂向力與力矩通過不斷的變化和調整最終達到穩定狀態,即三體船以某種新的姿態航行。常規的CFD理論計算并未考慮模型的升沉、縱傾變化對其阻力性能的影響,這也是采用常規CFD方法計算高速三體船阻力與模型試驗結果產生偏差的重要原因之一。在實際航行中,三體船模型的浮態也會隨著航速不同而發生變化,三體船首部的射流以及主、側體間劇烈的興波干擾也會引起船體濕表面積的較大改變,在CFD計算中船體濕表面積的變化一般被忽視,這也常常造成船體摩擦阻力被低估。由此可知,在預報三體船這種高性能船阻力時,航行姿態必須加以考慮。

圖10 約束模型的拖曳試驗Fig.10 Towing test of the constraint model

圖11 阻力預報結果對比Fig.11 Comparison between the results of resistance prediction

4 結 論

本文基于六自由度求解RANS方程進行了三體船模型在自由和約束狀態下的數值模擬,從流場細節(船體橫向截面的速度云圖、濕表面積統計和沿船長的浪高統計等)分析了自由和約束狀態下模型阻力產生差別的原因,通過對自由和約束狀態下的模型阻力試驗進行了驗證說明,現得到如下結論:

(1)對于高速船型如三體船,隨著傅汝德數的增加,其航行姿態會發生變化,這種變化對船體興波阻力以及摩擦阻力會產生一定的影響,常規固定航行姿態的CFD數值模擬方法對船體興波阻力及摩擦阻力的計算均偏小。

(2)在高航速時,自由狀態下的數值模擬與船模實際航行狀況更為接近,計算結果互相吻合。另外,該數值模擬方法可以實現模型的多自由度運動,操作比較簡單,對高速艦船的快速性預報具有重要的意義,同時也能彌補線性理論在船舶高航速強非線性條件下的不足以及獲得模型試驗中難以得到的船體周圍流場信息。

(3)由數值模擬可以精確得到三體船沿主體、側體內側和側體外側沿船體的浪高曲線,避免了通過試驗測量的麻煩以及測量誤差,可為設計工作者進行船型優化以減小主、側體間的興波干擾提供有價值的參考。

[1]酈 云,盧曉平.高速三體船阻力性能研究[J].船舶力學,2007,11(2):191-198. Li Yun,Lu Xiaoping.An investigation on the resistance of high speed trimarans[J].Journal of Ship Mechanics,2007,11 (2):191-198.

[2]黃德波,周廣利,鄧 銳.關于船舶阻力阻力性能試驗結果換算的1+K問題[C].中國造船工程學會船舶力學學術委員會、水動力學重點實驗室、船舶振動噪聲重點實驗室、江蘇省綠色船舶重點實驗室,2013.

[3]黃德波.近年我國船舶阻力方面的若干研究[C].中國造船工程學會船舶力學學術委員會,2010.

[4]黃德波,張雨新,鄧 銳,李 佳.單體與三體高速船舶粘性流場數值模擬[J].哈爾濱工程大學學報,2010(6):683-688. Huang Debo,Zhang Yuxin,et al.Numerical simulation of viscous flow around high speed monohull and trimaran ships [J].Journal of Harbin Engineering University,2010(6):683-688.

[5]韓開佳,黃德波.三體船的興波阻力計算[J].哈爾濱工程大學學報,2000(1):6-10. Han Kaijia,Huang Debo.Wavemaking resistance calculation of trimaran[J].Journal of Harbin Engineering University, 2000(1):6-10.

[6]王 中,盧曉平,王 瑋.非線性興波數值方法在高速三體船側體布局方案比較中的應用[J].船舶力學,2010,14(8): 863-871. Wang Zhong,Lu Xiaoping,Wang Wei.Application of the nonlinear wave making numerical method in the high-speed trimaran side hull position optimization[J].Journal of Ship Mechanics,2010,14(8):863-871.

[7]王 中,盧曉平,詹金林.高速三體船的水動力學和船型研究新進展[J].船舶力學,2011,15(7):813-826. Wang Zhong,Lu Xiaoping,Zhan Jinlin.New development on the investigation of high speed trimaran hydrodynamics and hull form[J].Journal of Ship Mechanics,2011,15(7):813-826.

[8]王 中,盧曉平,王 瑋.三體船興波阻力計算的自由面網格快速生成[J].哈爾濱工程大學學報,2010(4):409-413. Wang Zhong,Lu Xiaoping,Wang Wei.Fast free-surface mesh generation for the calculation trimaran wavemaking resistance[J].Journal of Harbin Engineering University,2010(4):409-413.

[9]盧曉平,王 鵬,詹金林.高速三體船片體構型數學表達和興波阻力計算[J].中國艦船研究,2013(1):13-19. Lu Xiaoping,Wang Peng,Zhan Jinlin.A calculation method of wave-making resistance and the demihull designing for high speed trimarans[J].Chinese Journal of Ship Research,2013(1):13-19.

[10]韓翔希,趙成璧,唐友宏,林 慰,曹藝齡.高速船航態模擬與阻力預報CFD方法應用[J].科學技術與工程,2013, 21:6176-6183. Han Xiangxi,Zhao Chengbi,et al.Navigation state simulation and resistance prediction of high speed vessels based on CFD simulation[J].Science Technology and Engineering,2013,21:6176-6183.

[11]姚朝幫,董文才.基于Fluent計及浮態變化的中高速船阻力預報方法研究[C].第七屆船舶力學學術委員會全體會議論文集,中國造船工程學會船舶力學學術委員會,2010.

[12]倪崇本,朱仁傳,繆國平,等.計及航行姿態變化的高速多體船阻力預報[J].水動力學研究與進展A輯,2011(1): 101-107. Ni Chongben,Zhu Renchuan,et,al.The resistance prediction for high speed multi-hull vessels with consideration of hull gesture variation during voyage[J].Chinese Journal of Hydrodynamics,2011(1):101-107.

[13]Joe L,Fred S.Resistance,Sinkage and trim,wave profile,and nominal wake tests and uncertainty assessment for DTMB model 5512[C]//Proc 25th ATTC.Iowa City,1998.

[14]Rodrigo Azcueta Repetto.Computation of turbulent free-surface flows around ships and floating bodies[D].PHD Thesis. Technischen Univ.Hamburg-Harburg,2001.

[15]陳 康,黃德波.CFD技術在三體船阻力性能研究中的應用[J].哈爾濱工程大學學報,2006(3):362-366. Chen Kang,Huang Debo.Application of CFD technology to the research of resistance performance of trimaran[J].Journal of Harbin Engineering University,2006(3):362-366.

[16]陳 康,黃德波,李云波.三體船阻力計算的改進方法研究[J].哈爾濱工程大學學報,2009(2):126-131. Chen Kang,Huang Debo,Li Yunbo.Discussion of methods for improving resistance calculation of trimaran hulls[J].Journal of Harbin Engineering University,2009(2):126-131.

[17]趙 藤,孫 鵬,于 雷.基于CFD的三體船流場數值模擬可靠性研究[J].華中科技大學學報(自然科學版),2012, 11:81-83. Zhao Teng,Sun Peng,Yu Lei.Research on reliability of fluidity numerical simulation for trimaran vessel under the CFD [J].Hua Zhong Univ.of Sci&Tech.(Natural Science Edition),2012,11:81-83.

[18]周利蘭,高 高,尹 巍.三體船興波問題的數值計算[J].船舶力學,2012,16(8):853-859. Zhou Lilan,Gao Gao,Yin Wei.Numerical calculation of the wave-making problem for trimarans problem for trimarans [J].Journal of Ship Mechanics,2012,16(8):853-859.

[19]周廣利,黃德波,鄧 銳,馬 勇,由世洲.三體船阻力性能的模型系列試驗研究[J].哈爾濱工程大學學報,2010,31 (5):576-584. Zhou Guangli,Huang Debo,et al.Investigating resistance of a range of trimaran designs[J].Journal of Harbin Engineering University,2010,31(5):576-584.

[20]王詩洋,王 超,常 欣,黃 勝.CFD技術在船舶阻力性能預報中的應用[J].武漢理工大學學報,2010,21(77):80-93. Wang Shiyang,Wang Chao,et al.Application of CFD technology to the research of resistance performance[J].Journal of Wuhan University of Technology,2010,21(77):80-93.

Resistance prediction method research of high speed trimaran

ZHOU Guang-li1,AI Zi-tao2,DENG Rui1,DOU Peng-fei1,HUANG De-bo1
(1.College of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China;2.Technical Center, Afai Southern Shipyard(Panyu Guangzhou)Ltd,Guangzhou 511431,China)

For the High speed trimarans,navigation attitude will change with the increase of Froude number. In order to design the attitude of the float state of structure as a benchmark,there is a certain error for fixed shipping normal CFD resistance calculation forecasting arc resistance.Aimed at tourism body form,using the free model state and validating the resistance test of the fixed model state analysis,get the difference.In order to further understand the laws for wave-making and frictional resistance,0.3<Fn<0.6 within the scope of the state of free model and fixed model respectively through numerical simulation,the flow field of the details is analyzed.Finally,the model test was carried out to verify the results.

trimaran;hull attitude;CFD;resistance prediction

U661.31

:Adoi:10.3969/j.issn.1007-7294.2016.07.003

1007-7294(2016)07-0805-11

2016-03-11

工信部高技術船舶科研計劃資助項目(G014612002)

周廣利(1969-),男,博士,副教授,通信作者,E-mail:zhouguangli@hrbeu.edu.cn;艾子濤(1989-),男,碩士,助理工程師。

猜你喜歡
船舶模型
一半模型
計算流體力學在船舶操縱運動仿真中的應用
基于改進譜分析法的船舶疲勞強度直接計算
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
船舶!請加速
BOG壓縮機在小型LNG船舶上的應用
船舶壓載水管理系統
中國船檢(2017年3期)2017-05-18 11:33:09
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产新AV天堂| 国产99欧美精品久久精品久久| 重口调教一区二区视频| 91精品久久久无码中文字幕vr| 国产成人免费| 国产麻豆永久视频| 亚洲精品无码AⅤ片青青在线观看| 久久综合色88| 亚洲欧美日韩天堂| 国产精品免费久久久久影院无码| 国产亚洲精品精品精品| 久久成人国产精品免费软件| a级毛片免费看| 一级不卡毛片| 亚洲热线99精品视频| 久久青草免费91线频观看不卡| 在线观看国产小视频| 国产91av在线| 一级毛片在线播放免费观看| 国产99视频免费精品是看6| 欧美亚洲国产精品第一页| 日韩成人在线一区二区| 欧美成人二区| 视频二区中文无码| 免费AV在线播放观看18禁强制| av在线5g无码天天| 日韩精品资源| 中文字幕不卡免费高清视频| 国产区在线观看视频| 成年人免费国产视频| 亚洲综合色吧| 色悠久久综合| 久草网视频在线| 18黑白丝水手服自慰喷水网站| 国产极品美女在线播放| 亚洲一区黄色| 国产欧美成人不卡视频| 久久亚洲黄色视频| 狠狠干综合| 婷婷综合亚洲| 国产国产人成免费视频77777| 蜜芽一区二区国产精品| 高潮毛片免费观看| 国产剧情一区二区| 91在线精品免费免费播放| 免费毛片在线| 亚洲区第一页| 农村乱人伦一区二区| 美女无遮挡免费网站| 日韩av手机在线| 亚洲va在线观看| 久久99国产精品成人欧美| 久久青草免费91线频观看不卡| 波多野结衣一区二区三区AV| 国产亚洲视频免费播放| 亚洲日韩日本中文在线| 成人年鲁鲁在线观看视频| 热99re99首页精品亚洲五月天| 欧美性色综合网| a级毛片视频免费观看| 国产成年无码AⅤ片在线| 精品一区二区三区水蜜桃| 3344在线观看无码| 国产人在线成免费视频| 久久精品丝袜| 伊人成人在线| 亚洲无限乱码一二三四区| 午夜毛片福利| 国产精品免费电影| 凹凸精品免费精品视频| 国产无码精品在线播放| 亚洲三级成人| 五月婷婷伊人网| 91在线播放免费不卡无毒| 真实国产乱子伦视频| 2020最新国产精品视频| 狠狠干欧美| 久久综合色88| 亚洲成a人片| 国产成人久久综合777777麻豆| 婷五月综合| 无码专区在线观看|