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

船舶運動狀態(tài)下舵附推力鰭水動力性能分析

2017-01-10 02:11:08王文全王詩洋張?zhí)焯?/span>
船舶 2016年6期
關(guān)鍵詞:船舶模型系統(tǒng)

王文全王詩洋張?zhí)焯?/p>

(1.中國船舶及海洋工程設(shè)計研究院 上海200011;2.大連職業(yè)技術(shù)學(xué)院 大連116035)

船舶運動狀態(tài)下舵附推力鰭水動力性能分析

王文全1王詩洋1張?zhí)焯?

(1.中國船舶及海洋工程設(shè)計研究院 上海200011;2.大連職業(yè)技術(shù)學(xué)院 大連116035)

采用基于RANSE的數(shù)值波浪水池技術(shù)和數(shù)值拖曳水池技術(shù),針對75 000 t船模型開展帶有縱搖和橫搖運動的自航CFD模擬研究。首先進行靜水中的船模阻力和自航CFD模擬,通過與模型試驗結(jié)果比較,初步驗證了文中的CFD模擬方法;隨后開展不同運動狀態(tài)下是否帶有舵附推力鰭的船模自航CFD模擬,并對該種狀態(tài)下的船模阻力、運動響應(yīng)以及船后螺旋槳的推力和扭矩等進行分析,獲得相關(guān)變化特性。該項研究對于運動中船舶快速性研究和預(yù)報等,具有較好的指導(dǎo)意義和參考價值。

舵附推力鰭;縱搖;橫搖;自航

引 言

進入21世紀(jì),能源問題已引起世界各國的高度重視,成為世界范圍內(nèi)的焦點問題之一,而節(jié)能減排儼然被視作經(jīng)濟發(fā)展中的一項緊迫任務(wù)。因此,研究船舶推進節(jié)能技術(shù)對于降低船舶能源消耗、節(jié)約船舶運輸企業(yè)的燃油消耗費用,以及降低船舶對海洋環(huán)境的污染都具有重大的意義。

以往許多學(xué)者和研究人員對舵附推力鰭的節(jié)能效果分析均關(guān)注船舶在靜水中的節(jié)能效果,而在進行數(shù)值計算時并未考慮船體在波浪中航行時的六個自由度的運動。靜水工況的數(shù)值計算結(jié)果并不能反映實際船舶在海面航行時舵附推力鰭的節(jié)能效果[1-2]。

目前已發(fā)表的有關(guān)波浪中自航CFD模擬方面的文獻主要來自于美國IOWA大學(xué)的Stern團隊。Carrica等人開展了靜水中帶自由升沉、縱傾KCS船模自航以及頂浪中船模運動響應(yīng)的CFD模擬[3],計算分析了模型自航點下的船模升沉量、縱傾角和螺旋槳轉(zhuǎn)速,以及規(guī)則波頂浪中船模的阻力和運動響應(yīng)。

本文主要工作內(nèi)容包括船-槳-舵系統(tǒng)在靜水、迎浪和橫浪狀態(tài)下的水動力數(shù)值預(yù)報、計算,以及船-槳-舵-舵附推力鰭系統(tǒng)在靜水及不同運動狀態(tài)下的節(jié)能效果數(shù)值預(yù)報,這些問題均考慮流體的粘性及自由液面產(chǎn)生的影響[4]。文中,我們建立船體、槳、舵-舵附推力鰭系統(tǒng)模型,導(dǎo)入STAR CCM+軟件中進行網(wǎng)格劃分和數(shù)值計算,對船-槳-舵和船-槳-舵-舵附推力鰭系統(tǒng)進行計及自由液面的阻力、螺旋槳推力和螺旋槳轉(zhuǎn)矩的數(shù)值計算;然后設(shè)置船體強迫縱搖與橫搖模擬不同波浪參數(shù),對船-槳-舵系統(tǒng)和船-槳-舵-舵附推力鰭系統(tǒng)進行不同運動狀態(tài)下的阻力、螺旋槳推力和螺旋槳轉(zhuǎn)矩的數(shù)值計算,對比其在靜水中和不同運動狀態(tài)下的推進節(jié)能效果,并加以分析。

1 數(shù)學(xué)模型

1.1 湍流模型

湍流模型可分為空間平均模型和時間平均模型。空間平均模型即大渦模擬,時間平均模型分為湍流粘性模型和應(yīng)力模型。其中應(yīng)力模型即雷諾應(yīng)力模型(RSM),湍流粘性模型又分為線性湍流模型和非線性湍流模型[5]。考慮到船尾部槳后旋轉(zhuǎn)尾流,本文采用雷諾應(yīng)力方程湍流模型。

雷諾應(yīng)力輸運方程在直角坐標(biāo)系下的表達式為:

該方程從左到右依次為:對流項、湍流擴散項、分子粘性擴散項、剪應(yīng)力產(chǎn)生項、浮力產(chǎn)生項、壓力應(yīng)變項、粘性耗散項、系統(tǒng)旋轉(zhuǎn)產(chǎn)生項。

雷諾應(yīng)力方程模型可直接使用輸運方程求解出雷諾應(yīng)力,從而避免其他模型的粘性假設(shè)。可見,雷諾應(yīng)力方程模型模擬強旋流較其他模型有明顯優(yōu)勢[6]。

1.2 VOF方法

VOF方法是一種以計算流體占據(jù)網(wǎng)格單元體積份額的途徑來捕捉自由液面的方法,由Hirt和Nichols[7]于20世紀(jì)80年代初提出。控制體中含有水和空氣,空氣的體積分?jǐn)?shù)是αa,水的體積分?jǐn)?shù)是αw,那么水的體積分?jǐn)?shù)輸運方程為:

控制體中只充滿水和空氣這兩種物質(zhì),所以水和空氣的體積分?jǐn)?shù)之和為1,即

采用VOF方法計算自由液面所需計算機的存儲量小,有利于三維自由液面問題的模擬分析。其不足之處是不易獲得準(zhǔn)確的自由液面高度位置。

1.3 重疊網(wǎng)格方法

重疊網(wǎng)格方法(Overset grid)的最終目標(biāo)是建立各重疊網(wǎng)格間的耦合關(guān)系,為各區(qū)域流場計算提供邊界信息的傳遞服務(wù)。

如圖1所示,網(wǎng)格1 落入網(wǎng)格2物面內(nèi)的點被標(biāo)記為洞內(nèi)點而不參與流場的計算。網(wǎng)格1中與洞點相鄰的點為洞邊界點(圖1中的方點),這些點通過插值接受網(wǎng)格2傳遞的流場信息。網(wǎng)格2的外邊界點(圖1中的圓點)也通過插值相應(yīng)地接受網(wǎng)格1傳遞的流場信息, 任意變量f由相應(yīng)的六面體貢獻單元采用Trilinear 方法插值獲得。網(wǎng)格1的洞邊界點和網(wǎng)格2的外邊界插值點之間的區(qū)域為重疊區(qū)域。

圖1 重疊網(wǎng)格示意圖

重疊網(wǎng)格將計算流域分成幾個子域,各個子域的網(wǎng)格獨立生成,并通過各域網(wǎng)格之間的嵌套、重疊等關(guān)系進行信息的傳遞,在重疊網(wǎng)格邊界區(qū)域流場的信息通過插值進行耦合。重疊網(wǎng)格的生成主要包括兩個方面:一方面,對背景域網(wǎng)格進行“挖洞”操作,即屏蔽洞面內(nèi)部的區(qū)域,將洞面內(nèi)網(wǎng)格進行標(biāo)記并在CFD計算中予以舍棄;另一方面,洞邊界之間的信息傳遞通過插值進行數(shù)值計算,即“尋點”問題。

2 計算模型的建立及網(wǎng)格劃分

2.1 幾何模型的建立

本文計算所選用的是75 000 t散貨船,計算船模基本參數(shù)如表1所示。計算中使用的螺旋槳為自行設(shè)計的不等螺距螺旋槳(PBC1),螺旋槳模型的幾何要素見表2。文中,舵附推力鰭的葉剖面類型為Clark-Y翼型,展長0.065 m、安裝角度4°、弦長0.01 m,具體參數(shù)見表3。

表1 船模基本參數(shù)

表2 原型槳(PBC1)的主要幾何參數(shù)

表3 Clark-Y翼型型值表%

計算所用的船-槳-舵系統(tǒng)模型見圖2。

圖2 船-槳-舵系統(tǒng)模型

2.2 計算流域及網(wǎng)格的劃分

本文采用重疊網(wǎng)格技術(shù),將船體包含在一個小域中,然后再建立一個適當(dāng)?shù)谋尘坝颍瑢⑿∮虻倪吔缗c背景域的交界面生成Interface。該Interface邊界條件類型設(shè)置成Overset Mesh Interface,使包含船體的小域按照設(shè)定的規(guī)則運動,同時在Interface的兩個邊界面上實現(xiàn)非定常流場的數(shù)值傳遞。計算船模自航時,還要在槳軸后給內(nèi)外螺旋槳分別創(chuàng)建一個圓柱形小域。小域與槳轂同軸,其直徑為螺旋槳直徑的1.2倍,通過滑移網(wǎng)格方法模擬螺旋槳旋轉(zhuǎn)運動。流體域、船體螺旋槳和開爾文波形區(qū)域具體如圖3所示。

圖3 網(wǎng)格劃分

3 靜水中的船-槳-舵系統(tǒng)水動力分析

3.1 靜水中的裸船體阻力預(yù)報

分別計算船模在V=1.057 m/s、V=1.098 m/s、V=1.139 m/s,V=1.179 m/s、V=1.220 m/s、V=1.261 m/s、V=1.301 m/s以及V=1.342 m/s這八個速度下的阻力值,并同試驗值進行對比分析。

表4為計算結(jié)果與試驗值的比較。可見,船體阻力計算值與試驗結(jié)果吻合較好, 最大相對偏差為2.6%,這說明本文所用方法可靠性佳。

表4 裸船體阻力值同試驗值對比

3.2 靜水中的船-槳-舵系統(tǒng)水動力分析

3.2.1 船槳舵系統(tǒng)水動力性能驗證

以JBC(Japan Bulk Carrier)船模為計算對象,其中JBC船模被選為2015東京船舶水動力CFD Workshop(Tokyo 2015-A Workshop on CFD in Ship Hydrodynamics)的計算標(biāo)模之一,其幾何外形如圖4所示。船模長7 m,螺旋槳模型直徑0.203 m(實槳直徑8.12 m)。AU型槳葉剖面的五葉槳幾何外形如圖5所示,主尺度參數(shù)見參考文獻[8]。

圖4 JBC船模外形

圖5 螺旋槳外形

船模靜水中自航CFD模擬中的網(wǎng)格劃分及邊界條件設(shè)置參照3.1節(jié)中船模靜水阻力設(shè)置,得到CFD計算的靜水中自航船模阻力、螺旋槳推力與轉(zhuǎn)矩(見表5)。通過與試驗值對比分析可以看出,CFD計算的船模阻力與模型試驗結(jié)果非常接近,兩者之間的差別在1%左右。螺旋槳推力、扭矩與模型試驗結(jié)果也相當(dāng)接近,差別在3%左右。

表5 船模靜水中的自航CFD計算結(jié)果

通過對以上靜水中船模阻力和自航CFD模擬結(jié)果的分析,初步證明本文的CFD計算方法能夠進行船模靜水狀態(tài)下的自航模擬,且模擬結(jié)果具有相當(dāng)高的精度。以上研究工作為后期開展相似船型船模靜水及運動自航CFD模擬奠定了技術(shù)基礎(chǔ)。

3.2.2 船槳舵系統(tǒng)水動力分析

通過生成網(wǎng)格并進行邊界條件的設(shè)置,采用船速1.220 m/s、螺旋槳轉(zhuǎn)速10.141 r/min(下面設(shè)置一致)。計算結(jié)果得到船-槳-舵系統(tǒng)水動力結(jié)果,見表6。

表6 船-槳-舵系統(tǒng)水動力數(shù)值計算結(jié)果(靜水)

靜水狀態(tài)中,將船-槳-舵-舵附推力鰭系統(tǒng)的計算結(jié)果與船-槳-舵系統(tǒng)的計算結(jié)果比較后發(fā)現(xiàn),前者船體阻力增加0.75 N、螺旋槳的推力增加了0.99 N,螺旋槳的轉(zhuǎn)矩從數(shù)據(jù)上看沒有變化。計算結(jié)果表明安裝舵附推力鰭對螺旋槳的轉(zhuǎn)矩沒有明顯的影響,船體阻力以及螺旋槳推力均有所增加。

圖6為船體表面的自由液面的計算結(jié)果圖形,紅色部分和藍色部分的交界處為自由液面,圖7為開爾文波形圖。舵附推力鰭對于船體自由液面的分布沒有明顯效果。從自由液面位置圖可以明顯看出在船體首柱位置自由液面明顯升高,船首處波面升高最大是由于船首處為船體周圍流線的駐點所在的位置,流體壓力最大。為平衡這一壓力升高,流體液面升高,以平衡此處的壓力變化。液面升高后又逐漸下降然后再升高,從船體中部往后,液面升高漸趨不明顯。

圖6 船體表面自由液面

從圖8的十字型參考線明顯可見船-槳-舵-舵附推力鰭系統(tǒng)的槳后盤面處的速度分布更接近左右對稱,而船-槳-舵系統(tǒng)的槳后盤面處的速度分布則逆時針旋轉(zhuǎn)一定角度,表明船-槳-舵系統(tǒng)的槳后尾流的旋轉(zhuǎn)比船-槳-舵-舵附推力鰭系統(tǒng)更為強烈,從而浪費螺旋槳輸出能量[9-14]。

圖7 開爾文波形圖

圖8 螺旋槳后盤面處速度分布云圖

4 運動中的船-槳-舵系統(tǒng)水動力分析

船舶運動通過重疊網(wǎng)格技術(shù)來實現(xiàn),通過在靜水中強制船舶搖蕩來模擬其在不同波浪條件下的運動狀態(tài),計算船舶在不同運動狀態(tài)下,船-槳-舵-舵附推力鰭系統(tǒng)中螺旋槳的推力、船體的阻力、螺旋槳的轉(zhuǎn)矩。本文基于Overset網(wǎng)格技術(shù)對船槳推進系統(tǒng)旋轉(zhuǎn)與縱搖(橫搖)的耦合運動進行數(shù)值模擬,計算過程中考慮了重力的影響,對該船槳推進系統(tǒng)在縱搖運動(縱搖最大角度1°,周期10 s)、橫搖運動(橫搖最大角度5°,周期5 s)兩種工況下的水動力性能進行計算和分析。船體縱搖運動的公式定義為:

圖9為船槳推進系統(tǒng)在振蕩運動狀態(tài)下的運動示意圖,圖中標(biāo)出了下文分析中會出現(xiàn)的四個時刻所對應(yīng)的位置,即t/T=0.00(與t/T=1.00相同)、t/T=0.25、t/T=0.5、t/T=0.75。

圖9 船體縱搖(橫搖)運動示意圖

4.1 縱搖下的船-槳-舵系統(tǒng)水動力分析

通過生成網(wǎng)格并進行邊界條件的設(shè)置,采用船速1.220 m/s,螺旋槳轉(zhuǎn)速10.141 r/min(與前面設(shè)置一致)。推力和轉(zhuǎn)矩收斂曲線見圖10和圖11。

圖10 帶有舵附推力鰭的船槳舵系統(tǒng)在縱搖運動下推力值收斂曲線

圖11 帶有舵附推力鰭的船槳舵系統(tǒng)在縱搖運動下轉(zhuǎn)矩值收斂曲線

由于船體進行規(guī)律性的運動,阻力與推力都是周期性變化,選取一個周期,取平均值得到船體阻力與螺旋槳推力和轉(zhuǎn)矩。數(shù)值計算結(jié)果見表7。

表7 船-槳-舵系統(tǒng)水動力數(shù)值計算結(jié)果(縱搖)

縱搖狀態(tài)下,將船-槳-舵-舵附推力鰭系統(tǒng)的計算結(jié)果與船-槳-舵系統(tǒng)的計算結(jié)果比較后發(fā)現(xiàn),前者船體阻力增加1.11 N、螺旋槳的推力增加了0.22 N,螺旋槳的轉(zhuǎn)矩從數(shù)據(jù)上看沒有變化。計算結(jié)果表明,安裝舵附推力鰭對螺旋槳的轉(zhuǎn)矩沒有明顯的影響,船體阻力以及螺旋槳推力均有所增加。

4.2 橫搖下的船-槳-舵系統(tǒng)水動力分析

橫搖狀態(tài)下,數(shù)值計算結(jié)果見表8,推力和轉(zhuǎn)矩收斂曲線見圖12和圖13。

表8 船-槳-舵系統(tǒng)水動力數(shù)值計算結(jié)果(橫搖)

圖12 帶有舵附推力鰭的船槳舵系統(tǒng)在橫搖運動下推力值收斂曲線

圖13 帶有舵附推力鰭的船槳舵系統(tǒng)在橫搖運動下扭矩值收斂曲線

將船-槳-舵-舵附推力鰭系統(tǒng)的計算結(jié)果同船-槳-舵系統(tǒng)的計算結(jié)果比較后發(fā)現(xiàn):前者船體阻力增加1.88 N,螺旋槳的推力增加0.09 N,螺旋槳的轉(zhuǎn)矩從數(shù)據(jù)上看并沒有變化。計算結(jié)果表明,安裝舵附推力鰭對螺旋槳的轉(zhuǎn)矩沒有明顯影響,船體阻力以及螺旋槳推力均有所增加。

5 結(jié) 論

本文進行了船-槳-舵系統(tǒng)以及船-槳-舵-舵附推力鰭系統(tǒng)的水動力計算,旨在通過計算安裝舵附推力鰭的船舶在不同運動狀態(tài)下螺旋槳的推力、轉(zhuǎn)矩,并與靜水工況下的船-槳-舵系統(tǒng)螺旋槳推力、轉(zhuǎn)矩進行對比分析,獲得船舶不同運動狀態(tài)下舵附推力鰭的水動力性能,并得出如下結(jié)論:

(1)靜水狀態(tài)與運動狀態(tài)下舵附推力鰭有助于提高螺旋槳效率。

(2)船體縱搖、橫搖運動狀態(tài)對于螺旋槳的效率影響較大,且船體阻力也會變大。

(3)本文采用強迫船體產(chǎn)生橫搖、縱搖運動的方法來模擬船舶在波浪中的運動,與船舶在波浪中的真實運動有一定的區(qū)別。

在后續(xù)工作中,將嘗試運用動態(tài)流體固體相互作用(DFBI,Dynamic Fluid Body Interaction)來模擬船舶在波浪中的六自由度搖蕩運動,設(shè)置不同的波浪參數(shù),使計算結(jié)果更接近于船舶在波浪中航行的真實情況。

[ 1 ]王超.船舶特種推進器與附加節(jié)能裝置[M].哈爾濱:哈爾濱工程大學(xué)出版社,2013:40-60.

[ 2 ]胡海洲,蘇玉民,沈海龍,等.船體影響下舵附推力鰭的節(jié)能效果研究[J].艦船科學(xué)技術(shù),2016(5):67-72.

[ 3 ]Carrica P M, Fu Huiping, Stern F.Computations of selfpropulsion free to sink and trim and of motions in head waves of the KRISO Container Ship(KCS)model[J].Applied Ocean Research, 2011,33: 309-320.

[ 4 ]趙大剛.船-槳-舵-舵球式推力鰭系統(tǒng)水動力性能數(shù)值計算[D].哈爾濱:哈爾濱工程大學(xué),2012.

[ 5 ]Giulio Dubbioso, Roberto Muscari, Andrea Di Mascio.Analysis of the performance of a marine propeller operating in oblique flow[J].Computer & Fluid,2013,75:86-102.

[ 6 ]Lin Z, Wang X, Yuan X, et al.Investigation and improvement of the staggered labyrinth seal[J].Chinese Journal of Mechanical Engineering, 2015(2):402-408.

[ 7 ]倪崇本.基于CFD的船舶阻力性能綜合研究[D].上海交通大學(xué),2011.

[ 8 ]National Maritime Research Institute(NMRI).Tokyo 2015-A Workshop on CFD in Ship Hydrodymics[EB/OL].Tokyo,[December2-4,2015].http://www.t2015.nmri.go.jp/index.html.

[ 9 ]黃勝.船舶推進節(jié)能技術(shù)與特種推進器[M].哈爾濱:哈爾濱工程大學(xué)出版社,2007.

[10]黃勝,郭春雨.船舶推進節(jié)能技術(shù)研究與進展[J].艦船科學(xué)技術(shù),2007(1):27-32.

[11]郭春雨,黃勝,常欣.非定常流中槳后舵附推力鰭的數(shù)值模擬[J].哈爾濱工程大學(xué)學(xué)報,2007(3):259-262.

[12]郭春雨.螺旋槳與舵附推力鰭相互干擾水動力性能數(shù)值計算[D].哈爾濱工程大學(xué),2007.

[13]郭春雨,黃勝.面元法計算槳后舵附推力鰭水動力性能[J].哈爾濱工程大學(xué)學(xué)報,2006(4):501-504.

[14]郭春雨,黃勝.非線性渦格法預(yù)報槳后舵附推力鰭水動力性能[J].華中科技大學(xué)學(xué)報(自然科學(xué)版),2006(6):87-89.

Analysis of hydrodynamic performance of rudder with additional thrust f ns for moving ships

WANG Wen-quan1WANG Shi-yang1ZHANG Tian-tian2
(1.Marine Design & Research Institute of China, Shanghai 200011, China; 2.Dalian Vocational & Technical College, Dalian 116035, China)

The CFD simulations of the self-propulsion test with pitch and rolling for a 75 000 t model ship were carried out by the numerical wave tank technology and numerical towing tank technology based on RANSE.The resistance test and self-propulsion test simulations for the ship in calm water were firstly performed and compared with the model test results to preliminarily validate the CFD simulation method.Then, the self-propulsion tests for the model ships with or without rudders with additional thrust fins under the different motion modes were simulated to analyze the resistance and motion responses of the model ship, and the thrust and torque of the aft propeller, resulting in the relevant variation characteristics.The research can provide guidance and reference for the study and prediction of powering performance of ships in waves.

rudder with additional thrust fins; pitch; rolling; self-propulsion

U661.1

A

1001-9855(2016)06-0019-09

2016-04-18;

2016-07-19

王文全(1987-),男,博士,工程師,研究方向:船舶總體設(shè)計。王詩洋(1981-),男,碩士,高級工程師,研究方向:船舶總體設(shè)計。張?zhí)焯恚?986-),女,碩士,實驗師,研究方向:船舶性能試驗。

10.19423/j.cnki.31-1561/u.2016.06.019

猜你喜歡
船舶模型系統(tǒng)
一半模型
計算流體力學(xué)在船舶操縱運動仿真中的應(yīng)用
Smartflower POP 一體式光伏系統(tǒng)
《船舶》2022 年度征訂啟事
船舶(2021年4期)2021-09-07 17:32:22
WJ-700無人機系統(tǒng)
ZC系列無人機遙感系統(tǒng)
北京測繪(2020年12期)2020-12-29 01:33:58
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
船舶!請加速
連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini
主站蜘蛛池模板: 欧美a在线视频| 国产在线一二三区| 国产精品嫩草影院视频| 免费观看三级毛片| 青青草原国产免费av观看| 蜜芽一区二区国产精品| 亚洲精品你懂的| 麻豆国产精品| 97亚洲色综久久精品| 999精品色在线观看| 亚洲国产成熟视频在线多多 | h视频在线播放| 欧美在线综合视频| 久久香蕉国产线| 中文字幕伦视频| 午夜国产理论| 久久夜色精品国产嚕嚕亚洲av| 亚洲视频免费在线看| 午夜激情福利视频| 综合色区亚洲熟妇在线| 免费av一区二区三区在线| 欧美成人免费| 久久这里只有精品2| 色精品视频| www.av男人.com| 高清不卡一区二区三区香蕉| 情侣午夜国产在线一区无码| 亚洲制服丝袜第一页| 国产免费福利网站| 亚洲AV永久无码精品古装片| 午夜国产精品视频| 五月天综合网亚洲综合天堂网| www亚洲精品| 国产SUV精品一区二区| 国产高清精品在线91| 99久久精品免费观看国产| 久久美女精品国产精品亚洲| 成人午夜视频网站| 深夜福利视频一区二区| 美女被狂躁www在线观看| 一级毛片不卡片免费观看| 四虎精品国产永久在线观看| 精品国产香蕉在线播出| 四虎成人精品在永久免费| 国产精品久久久久久搜索 | 97视频在线观看免费视频| 亚洲男人在线天堂| 欧洲成人在线观看| 亚洲日韩高清在线亚洲专区| 久久久无码人妻精品无码| 国产福利小视频在线播放观看| 欧美另类第一页| 国产一区二区三区夜色| a网站在线观看| 日韩av无码精品专区| 2021国产在线视频| 99精品一区二区免费视频| 成年av福利永久免费观看| 亚洲男人的天堂久久香蕉网| 一级一级特黄女人精品毛片| 日韩精品欧美国产在线| 97视频精品全国免费观看| 91尤物国产尤物福利在线| 婷婷午夜影院| 欧洲一区二区三区无码| 久久a毛片| 亚洲精品天堂在线观看| 欧美日韩国产一级| 亚洲综合第一页| 呦女精品网站| 国产拍在线| 波多野结衣一区二区三区四区视频| 91精品aⅴ无码中文字字幕蜜桃| 亚洲浓毛av| 91成人在线观看视频| 日本福利视频网站| 婷婷亚洲天堂| 国产流白浆视频| 亚洲国产亚洲综合在线尤物| 暴力调教一区二区三区| 欧美亚洲国产一区| 精品撒尿视频一区二区三区|