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

自主水下航行器回收過程中螺旋槳推力特性分析

2017-07-10 10:28:52杜曉旭張正棟
兵工學報 2017年6期
關鍵詞:模型

杜曉旭, 張正棟

(西北工業(yè)大學 航海學院, 陜西 西安 710072)

自主水下航行器回收過程中螺旋槳推力特性分析

杜曉旭, 張正棟

(西北工業(yè)大學 航海學院, 陜西 西安 710072)

在自主水下航行器 (AUV) 水下回收過程中,螺旋槳的推力特性對AUV的穩(wěn)定性控制及安全性能等影響非常明顯。主要基于計算流體力學方法,采用多塊混合網(wǎng)格技術結合RNGk-ε湍流模型,對AUV在不同回收工況下的螺旋槳推力特性展開數(shù)值模擬。建立AUV及螺旋槳的三維幾何模型;基于多塊混合網(wǎng)格方法對計算域進行空間離散,同時采用動網(wǎng)格技術實現(xiàn)螺旋槳與回收管的相對運動及自身旋轉運動;利用有限體積法展開數(shù)值計算。基于計算數(shù)據(jù)與試驗數(shù)據(jù)的對比及數(shù)值無關性測試,表明所使用的多塊混合網(wǎng)格方法能夠對AUV回收過程螺旋槳推力特性進行準確計算,為AUV水下回收過程的控制及安全性能評估能夠提供理論參考。

流體力學; 螺旋槳; 自主水下航行器; 多塊混合網(wǎng)格; 推力特性

0 引言

自主水下航行器(AUV)作為水下自主航行的移動節(jié)點,目前在海洋各個領域中的用途愈加廣泛。目前研究者們致力于AUV的水下自主回收研究,通常使用內徑略大于AUV外徑的管狀回收裝置實現(xiàn)AUV的回收[1]。AUV在回收過程中螺旋槳與回收裝置之間的間隙隨著AUV位置的不同而不同,從而導致該間隙流場變化十分劇烈,通常稱這種復雜流場為極限流場。隨著AUV向回收裝置運動,螺旋槳由無界流場向極限流場運動,其推力特性受該流場的影響十分明顯,推力特性的改變給AUV回收過程的穩(wěn)定性控制帶來了挑戰(zhàn),因此對該運動過程的螺旋槳推力特性研究十分有意義。

近年來隨著計算機技術的推廣普及和計算方法的新發(fā)展,計算流體力學(CFD)技術取得了蓬勃的發(fā)展,通過對雷諾平均N-S(RANS)方程的數(shù)值求解來獲得黏性流場全場信息已經成為可能,因此CFD方法已經廣泛應用到螺旋槳水動力特性研究中[2-4]。蔡榮泉等[5]基于非結構化網(wǎng)格對某五葉螺旋槳定常水動力性能進行CFD計算,計算結果與試驗誤差較小。解學參等[6]基于RANS方程和RNGk-ε湍流模型采用結構- 非結構多塊混合網(wǎng)格研究了船后非均勻流場中螺旋槳的水動力性能,并將計算結果和試驗結果比較,結果表明CFD計算方法精度較高,能夠滿足工程應用要求。Guo等[7]對不同軸深的螺旋槳展開數(shù)值模擬,并與試驗結果對比,驗證了數(shù)值模擬結果的準確性。蘇石川等[8]基于兩種網(wǎng)格(分區(qū)混合網(wǎng)格及非結構網(wǎng)格)研究某散貨船螺旋槳敞水性能,對比發(fā)現(xiàn)分區(qū)混合網(wǎng)格誤差較小??傊捎肅FD方法計算螺旋槳敞水水動力特性研究已相對成熟。

但是,目前采用CFD方法對螺旋槳水動力性能的研究主要集中在敞水性能方面,而AUV水下回收過程螺旋槳處于極限流域運動,且葉梢與管壁的間隙會隨著AUV的運動而發(fā)生變化,因此螺旋槳所處流場極其復雜,其推力性能與敞水流域的推力性能有明顯的差別。在上述過程中螺旋槳推力特性的改變對AUV的穩(wěn)定性控制及安全性能都造成很大的影響,因此急需展開AUV水下回收過程螺旋槳推力特性的研究,為AUV水下回收過程的控制及安全性能評估提供定量的參考。本文基于CFD數(shù)值模擬方法,提出采用多塊混合網(wǎng)格方法離散計算空間的計算方法,結合RNGk-ε湍流模型以及動網(wǎng)格技術利用有限體積法(FVM)研究不同回收工況下螺旋槳所產生的推力,分析上述過程對螺旋槳推力特性的影響,驗證所提方法計算極限流場螺旋槳推力特性的準確性。

1 控制方程與湍流模型

1.1 控制方程

本文假設AUV在有界流場中航行時,周圍流體為連續(xù)介質,為有黏不可壓縮流體。由于螺旋槳葉梢與管壁之間的間隙時刻變化,因此槳葉周圍流體視為不可壓縮流體的三維非定常運動。因此滿足控制方程[9]為

(1)

(2)

1.2 湍流模型

本文所采用的RNGk-ε湍流模型來源于嚴格地統(tǒng)計技術,它在ε方程中加了一個條件,從而有效地提高了精度。而且RNGk-ε湍流模型充分考慮了湍流漩渦對計算的影響,從而提高了其計算復雜旋轉流場的精度。這些特點使得RNGk-ε湍流模型計算復雜流場有更高的可信度和精度。因此本文非定常湍流計算采用RNGk-ε雙方程湍流模型使RANS方程封閉[9]。

湍流脈動動能方程(k方程)為

(3)

式中:αk為湍動能的有效普朗特數(shù)的倒數(shù);μeff=μ+μt,μt為湍動黏性系數(shù);Gk為由于平均速度梯度引起的湍動能;Gb是由于浮力影響引起的湍動能;YM為可壓縮湍流脈動膨脹對縱耗散率的影響。

湍流能量耗散率(ε方程)為

(4)

其中:

(5)

式中:Pk為湍動生成項;Rε與η為RNG湍流模型特有的控件;其他常數(shù)σε=1.39,Cε1=1.42,Cε2=1.42,Cμ=0.084 5,η0=4.38,β=0.012.

2 多塊混合網(wǎng)格劃分

2.1 幾何模型

本文幾何模型如圖1所示,具體尺寸如表1所示。以AUV頭部為原點建立參考坐標系,沿AUV向前運動軸向為x軸正方向,豎直向上為y軸正方向,根據(jù)右手準則確定z軸正方向,后文中的數(shù)值模擬計算均以AUV頭部距參考坐標系原點距離為變量比較螺旋槳推力系數(shù)。

圖1 幾何模型示意圖Fig.1 Schematic diagram of geometrical model

名稱尺寸AUV外徑/m02回收管長度/m3AUV長度/m185半錐角/(°)30螺旋槳直徑/m0175

2.2 流場及網(wǎng)格劃分

本文建立圓柱體流場計算域,長為29.85 m,直徑為14 m. 流場計算域及邊界條件設置如圖2所示。入口邊界條件為速度入口,出口設為壓力出口,AUV及螺旋槳各壁面均設為無滑移壁面,采用標準壁面函數(shù)進行處理。計算域介質為水,密度ρ=1 024 kg/m3,其運動黏性系數(shù)為ν=1.003×10-6m2/s. 本文所選螺旋槳直徑D=175 mm,定義螺旋槳推力系數(shù)為Ct=T/ρn2D4,T為螺旋槳產生推力,n為螺旋槳的轉速,本文中n取900 r/min. 通過編制用戶自定義函數(shù)(UDF)實現(xiàn)AUV及螺旋槳的自定義運動。

圖2 流場計算域示意圖Fig.2 Schematic diagram of calculation domain of flow fields

本文所用FVM進行數(shù)值求解具有良好的收斂性,能夠滿足計算需求。而采用FVM求解RANS方程實際是求解離散在計算域內各網(wǎng)格節(jié)點的流場信息,包括壓力、速度以及湍流度等。因此合適的網(wǎng)格劃分是偏微分方程準確求解的基礎,也是數(shù)值模擬與分析的載體。劃分合適的高質量網(wǎng)格在CFD數(shù)值計算中起著關鍵的作用,而網(wǎng)格的形狀和精密度等對CFD計算的精度和計算效率有重要影響。目前網(wǎng)格劃分主要分為結構網(wǎng)格和非結構網(wǎng)格,結構網(wǎng)格生成質量高但適應范圍較窄,而非結構網(wǎng)格能夠很好地適應不規(guī)則的幾何模型,網(wǎng)格生成快,但是質量難以控制。針對兩種網(wǎng)格的優(yōu)缺點,分析本文的幾何模型,考慮到本文流場的復雜性,使用單一的網(wǎng)格很難達到計算要求。于是本文采用結構網(wǎng)格和非結構網(wǎng)格相結合的多塊混合網(wǎng)格方法對整個計算域進行劃分,將流場分為靜止域及運動域兩部分,靜止域為回收管以外部分,采用三維六面體結構化網(wǎng)格進行劃分,運動域為回收管以內部分,動靜域之間設立Interface交界面。將運動域分為兩塊,切出一個包含螺旋槳的圓柱部分,采用三維非結構化網(wǎng)格,而其余部分采用六面體結構化網(wǎng)格劃分。

對于包含螺旋槳的圓柱部分采用非結構化網(wǎng)格劃分,先用全三角形網(wǎng)格生成表面網(wǎng)格,再生成四面體網(wǎng)格填充計算域,在螺旋槳槳葉生成三棱柱邊界層以更精確地捕獲流場特性,同時對于螺旋槳槳葉部分進行加密,流場網(wǎng)格劃分如圖3、圖4所示。

圖3 靜止域網(wǎng)格劃分Fig.3 Meshes of stationary domain

圖4 運動域網(wǎng)格劃分Fig.4 Motion domain meshes

3 數(shù)值計算結果及分析

3.1模型驗證

在CFD計算中,網(wǎng)格的尺寸以及時間步長對計算結果的影響十分明顯,為了驗證計算的可靠性,分別對網(wǎng)格無關性以及時間步長無關性進行驗證。本文將分別計算4種不同數(shù)量的網(wǎng)格的螺旋槳推力特性曲線,通過調節(jié)螺旋槳槳葉第1層網(wǎng)格尺寸、網(wǎng)格節(jié)點數(shù)目及最大網(wǎng)格尺寸等參數(shù),分別得到252萬、378萬、504萬以及756萬網(wǎng)格,其對應的螺旋槳槳葉第1層網(wǎng)格尺度分別為0.01D、0.005D、0.001D,0.000 5D. 計算結果如圖5(a)所示,隨著網(wǎng)格數(shù)量的增多,螺旋槳推力系數(shù)曲線逐漸收斂,當網(wǎng)格數(shù)量為504萬時,其計算結果與756萬網(wǎng)格數(shù)量的計算結果幾乎沒有差別,于是認為504萬的網(wǎng)格已達到網(wǎng)格無關,同時考慮到計算成本,因此本文取504萬的網(wǎng)格作為計算網(wǎng)格。同時可以發(fā)現(xiàn),當x<1.0 m時,螺旋槳推力系數(shù)保持不變,此時AUV與回收裝置尚未接觸,可以認為此時螺旋槳推力系數(shù)為敞水推力系數(shù)。本文所用螺旋槳根據(jù)日本MAU系列螺旋槳的圖譜所設計[10],從文獻[10]可看出在同進速系數(shù)下的螺旋槳試驗推力系數(shù)KT=0.28,而此時本文推力系數(shù)數(shù)值模擬結果為Ct=0.300 2,相對誤差不超過10%,可以認為數(shù)值模擬結果是準確的。如圖5(b)所示為不同時間步長對螺旋槳推力系數(shù)的影響曲線,分別為0.05 s、0.01 s、0.005 s以及0.001 s. 從圖5(b)中發(fā)現(xiàn),隨著時間步長的減小,推力系數(shù)曲線逐漸收斂。于是認為0.005 s的時間步長滿足計算要求,后文計算中均取此時間步長。

圖5 網(wǎng)格及時間步長無關性對比結果Fig.5 Comparative results of mesh and time step size independence

3.2 AUV勻速進入回收管時螺旋槳推力特性

本節(jié)的回收管內徑為220 mm,分別研究AUV以不同速度進入回收裝置時的螺旋槳推力特性。本節(jié)分別對AUV以1.0 m/s、1.5 m/s、2.0 m/s、2.5 m/s 4種不同的速度進入回收管展開數(shù)值模擬,圖6為AUV以不同速度進入回收管時螺旋槳推力特性曲線。

圖6 AUV以不同速度進入回收管時螺旋槳推力系數(shù)曲線Fig.6 Propeller’s thrust coefficient curves when AUV enters the tube at different velocities

圖7 不同位置處螺旋槳壓力等值線圖Fig.7 Propeller’s pressure contours at different positions

如圖6所示,隨著AUV完全進管,螺旋槳推力系數(shù)逐漸增大。當x≤1.0 m時,AUV頭部尚未進入錐頭,此時螺旋槳推力系數(shù)保持不變,這是由于螺旋槳與AUV均處于無界流域中,相對運動穩(wěn)定;當1.0 m

圖7為不同位置處的螺旋槳壓力等值線圖。由圖7可以看出,在螺旋槳還未進入錐頭時,葉片的導邊處有較大的壓力梯度,可知流體動力載荷分布在這里,而葉片的其余壓力分布較為平緩,過渡均勻,上述葉片壓力分布符合敞水流域螺旋槳葉片壓力分布規(guī)律。隨著AUV頭部靠近錐頭,流場出現(xiàn)干擾,流入盤面速度增加,產生的軸向誘導速度也逐漸增加,葉片與管壁之間的干擾使得壓力分布變得更加不均勻,同時吸力面的最大壓力下降,但是壓力梯度明顯上升,因此導致螺旋槳產生的推力上升。

如圖8所示為AUV勻速進入回收管時在不同位置的切面速度云圖,陰影部分為回收裝置。當AUV還未進入回收管時,流場穩(wěn)定,AUV及螺旋槳處于無界流場,螺旋槳推進特性不受回收裝置的影響。當AUV頭部剛要進入圓柱段時,回收裝置里的水出現(xiàn)回流,同時一部分水沿著AUV與管壁間隙加速噴出,此時由于頭部進管引起的流體加速逐漸開始影響螺旋槳盤面入口速度,導致螺旋槳推力系數(shù)開始增大。隨著AUV頭部進入回收管圓柱段,AUV與管壁的間隙流場逐漸影響到螺旋槳,螺旋槳產生的誘導速度變大,吸力面與壓力面的壓力梯度增大,導致螺旋槳產生的推力系數(shù)增大,這也與圖6的螺旋槳推力系數(shù)變化趨勢是一致的。

圖8 AUV不同位置切面速度云圖Fig.8 Velocity contours at different positions of AUV

3.3 AUV勻減速進入回收管時螺旋槳推力特性

在實際情況中,為了避免AUV與回收裝置發(fā)生碰撞,AUV回收時的速度是逐漸減小的,本節(jié)對AUV以不同初始速度勻減速進入回收裝置展開數(shù)值模擬,初始速度分別為1.0 m/s、1.5 m/s、2.0 m/s、2.5 m/s. 當x=4.0 m時,其速度降為0. 圖9為AUV以不同初始速度減速進入回收管時螺旋槳推力特性變化曲線。

圖9 不同初始速度對應的螺旋槳推力系數(shù)曲線Fig.9 Propellers thrust coefficient curves corresponding to different initial velocities of AUV

由圖9可知,螺旋槳的推力系數(shù)變化規(guī)律與AUV勻速運動時幾乎是一致的。當x≤1.0 m時,隨著AUV向前運動,其運動速度線性減小,導致螺旋槳進速系數(shù)線性減小,推力系數(shù)近似呈線性增加,符合螺旋槳敞水推進特性規(guī)律;隨著AUV進入回收裝置,螺旋槳推力系數(shù)的增長速度變快。AUV速度越大,螺旋槳推力系數(shù)變化越明顯,v=2.5 m/s所對應的螺旋槳推力系數(shù),最終大于v=2.0 m/s及v=1.5 m/s所對應的推力系數(shù)。由此可見,AUV減速運動時,螺旋槳推力系數(shù)變化較之勻速運動時更加明顯,極限流場對螺旋槳推力系數(shù)的影響也更加明顯。

3.4 間隙對螺旋槳推力特性的影響

本節(jié)選用單側間隙為3%、5%、7%、9%、11%的回收管模型進行數(shù)值計算,即回收管內徑分別為212 mm、220 mm、228 mm、236 mm、244 mm,分別計算AUV以2.0 m/s的速度進入回收管時螺旋槳推力系數(shù)。圖10為不同單側間隙下AUV進入回收管時螺旋槳推力系數(shù)變化曲線。

圖10 不同間隙下AUV進入回收管螺旋槳推力系數(shù)曲線Fig.10 Propeller’s thrust coefficient curves corresponding to different clearances between AUV and tube

圖10表明,隨著AUV進入回收管,AUV與回收管之間的間隙對螺旋槳推力特性的影響是非常明顯的,尤其當螺旋槳進入錐頭之后影響更為明顯。當AUV處于回收裝置之外時,間隙對螺旋槳推力系數(shù)的影響幾乎為0,這也與實際流場是相符的。隨著間隙的增大,間隙流場之間的干擾相對減小,從而螺旋槳的葉梢流場干擾相對減弱,從而引起推力系數(shù)的減小。同時可以發(fā)現(xiàn)當間隙大于7%的時候,螺旋槳推力系數(shù)有明顯的下降趨勢。

4 結論

本文針對螺旋槳在AUV回收過程中處于極限流域這一復雜流場問題展開的數(shù)值模擬,基于螺旋槳槳葉較為復雜的幾何外形,結合結構網(wǎng)格以及非結構網(wǎng)格的優(yōu)點,提出了能夠有效計算AUV回收過程螺旋槳推力特性的多塊混合網(wǎng)格方法。同時結合動網(wǎng)格技術以及RNGk-ε湍流模型,形成了計算AUV回收時螺旋槳推力特性的計算方法。隨后對計算域網(wǎng)格無關性以及離散時間步長無關性進行驗證,提高計算精度。最后對AUV以不同速度勻速進入回收裝置、不同初始速度勻減速進入回收裝置以及不同間隙等工況下的螺旋槳推力特性進行計算。結果表明,這種基于多塊混合網(wǎng)格方法計算AUV回收過程螺旋槳推力特性是可行的,而且對于極限流域的螺旋槳水動力特性研究具有普適性。

References)

[1] 潘光, 黃明明, 宋保維. AUV回收技術現(xiàn)狀及發(fā)展趨勢[J]. 魚雷技術, 2008, 16(6):10-14. PAN Guang, HUANG Ming-ming, SONG Bao-wei. Current situation and development trend of AUV recovery technology[J].

Torpedo Technology, 2008, 16(6):10-14. (in Chinese)

[2] Funeno I. Analysis of unsteady viscous flow around a highly skewed propeller[J].Journal of Kansai Society of Naval Architects, 2002, 237: 39-45.

[3] Bhattacharyya A, Krasilnikov V, Steen S. A CFD-based scaling approach for ducted propellers[J]. Ocean Engineering, 2016, 123:116-130.

[4] 鄭小龍, 黃勝, 王超. 基于CFD的螺旋槳定常水動力性能預報精度研究[J]. 艦船科學技術, 2014, 36(12): 11-15. ZHENG Xiao-long,HUANG Sheng,WANG Chao. Study of precision of steady hydrodynamic performance prediction of propeller of based on CFD[J]. Ship Science and Technology, 2014, 36(12): 11-15. (in Chinese)

[5] 蔡榮泉,陳鳳明,馮雪梅.使用FLUENT軟件的螺旋槳敞水性能計算和分析[J].船舶力學,2006, 10(5): 41-48. CAI Rong-quan,CHEN Feng-ming,F(xiàn)ENG Xue-mei.The use of FLUENT software of propeller open water perfomance calculation and analysis[J].Journal of Ship Mechanics, 2006, 10(5): 41-48.(in Chinese)

[6] 解學參, 姜治芳, 邱遼原. 非均勻粘性流場螺旋槳非定常水動力性能研究 [J]. 船舶工程,2015, 37(6): 37-40. XIE Xue-can, JIANG Zhi-fang, QIU Liao-yuan. Research on unsteady hydrodynamic performance of propeller on non-uniform viscous flow fields[J]. Ship Engineering, 2015, 37(6): 37-40. (in Chinese)

[7] Guo C, Zhao D, Sun Y. Numerical simulation and experimental research on hydrodynamic performance of propeller with varying shaft depths[J]. China Ocean Engineering, 2014, 28(2):271-282.

[8] 蘇石川,張力,魏承印. 基于不同網(wǎng)格類型的某散貨船螺旋槳敞水性能預報與研究[J]. 江蘇科技大學學報:自然科學版, 2016, 30(1):18-22. SU Shi-chuan, ZHANG Li, WEI Cheng-yin. Prediction of the open water performance of a bulk carrier propeller based on different grid types.[J] Journal of Jiangsu University of Science and Technology: Natural Science Edition, 2016, 30(1):18-22. (in Chinese)

[9] 劉志華, 熊鷹, 葉金銘. 基于多塊混合網(wǎng)格的RANS方法預報螺旋槳敞水性能的研究[J]. 水動力學研究與進展, 2007, 22(4): 50-56. LIU Zhi-hua, XIONG Ying, YE Jin-ming. Study on the prediction of propeller open-water performance using RANS formula and multi-block hybrid meshes[J]. Journal of Hydrodynamics, 2007, 22(4): 50-56. (in Chinese)

[10] 盛振邦, 劉應中. 船舶原理(下)[M].上海:上海交通大學出版社, 2004. SHENG Zhen-bang, LIU Ying-zhong. Ship theory (2) [M]. Shanghai: Shanghai Jiao Tong University Press, 2004. (in Chinese)

Analysis of Propeller Thrust Performance in the Process ofAutonomous Underwater Vehicle Recovery

DU Xiao-xu, ZHANG Zheng-dong

(School of Marine Science and Technology, Northwestern Polytechnical University, Xi’an 710072, Shaanxi, China)

The influence of propeller thrust performance on the stability control and safety performance of autonomous underwater vehicle (AUV) during the process of AUV recovery is very obvious. Based on CFD method, the multi-block meshes method combined with RNGk-εturbulent model is used for the simulation of propeller thrust performance of AUV in different recovery cases. The three-dimensional geometric models of AUV and propeller are built. The calculation domain is dispersed by multi-blocks meshes, and the dynamic mesh technique is used to simulate the motion of AUV and propeller. The numerical simulations are performed based on the finite volume method. The comparison of numerical and experimental results,as well as numerical independence test show that the proposed numerical method can be used to calculate the propeller thrust performance during the process of AUV underwater recovery accurately, and the numerical results can provide theoretical reference to the stability control and safety performance of underwater AUV recovery.

fluid mechanics; propeller; autonomous underwater vehicle; multi-blocks mesh; thrust performance

2016-10-18

國家自然科學基金項目(1130276)

張正棟(1993—),男,碩士研究生。E-mail:18789496770@163.com

杜曉旭(1981—),男,副教授,碩士生導師。E-mail:nwpudxx@163.com

U661.31+3

A

1000-1093(2017)06-1154-07

10.3969/j.issn.1000-1093.2017.06.015

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 成人午夜免费观看| 人妻无码AⅤ中文字| 97成人在线视频| 国产激情国语对白普通话| 性欧美在线| 国产高清不卡| 99精品伊人久久久大香线蕉 | 国产综合日韩另类一区二区| 91小视频在线观看免费版高清| 欧美午夜在线视频| 日本欧美在线观看| 亚洲天堂2014| 亚洲最大看欧美片网站地址| 亚洲大学生视频在线播放| 亚卅精品无码久久毛片乌克兰| 亚洲精品自产拍在线观看APP| 国产精品乱偷免费视频| 国产福利免费视频| 韩国v欧美v亚洲v日本v| 天堂av高清一区二区三区| 日本欧美中文字幕精品亚洲| 直接黄91麻豆网站| 福利片91| 在线亚洲精品福利网址导航| 色噜噜综合网| 国产午夜福利在线小视频| 欧洲一区二区三区无码| 国产一级视频久久| 欧美日本激情| 91麻豆国产视频| 亚洲va视频| 永久免费精品视频| 露脸真实国语乱在线观看| 国产成人一二三| 色丁丁毛片在线观看| 亚洲一级毛片免费观看| 无码日韩精品91超碰| 538精品在线观看| 日日噜噜夜夜狠狠视频| 国产亚洲欧美日韩在线一区二区三区| 国产无码精品在线播放| 白浆免费视频国产精品视频| 欧美激情第一欧美在线| 一本视频精品中文字幕| 国产无码在线调教| 婷婷色一二三区波多野衣| 免费看久久精品99| 亚国产欧美在线人成| 精品1区2区3区| 午夜少妇精品视频小电影| 黄色网在线免费观看| 国产精品视频3p| 狠狠躁天天躁夜夜躁婷婷| 精品一区二区三区波多野结衣| AV不卡无码免费一区二区三区| 日韩第一页在线| 国产精品漂亮美女在线观看| 亚洲黄色视频在线观看一区| 夜夜拍夜夜爽| 国产理论精品| 好吊日免费视频| 欧美亚洲国产精品第一页| 国产日本欧美在线观看| 播五月综合| 女同国产精品一区二区| 国产精品精品视频| 久久久久亚洲AV成人人电影软件| 毛片大全免费观看| 国产真实乱子伦精品视手机观看| 欧美精品啪啪一区二区三区| 99这里只有精品6| 在线色综合| 国产丝袜无码精品| 国产欧美在线观看一区| 8090成人午夜精品| 99九九成人免费视频精品 | 91福利在线看| 亚洲精品麻豆| 亚洲无码电影| 91久久青青草原精品国产| 亚洲国产综合第一精品小说| 国产福利小视频在线播放观看|