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

基于CFD的艦船遭遇畸形波數值模擬的驗證研究

2018-04-25 02:27:12張本輝石愛國
船舶力學 2018年4期

蔡 烽,張本輝,吳 明,楊 波,王 驍,石愛國

(海軍大連艦艇學院 航海系,遼寧 大連 116018)

0 引 言

畸形波(freak wave)是海洋中高且陡的大波,其持續時間很短,但出現的偶然性和巨大的破壞性,對船舶航運和海洋工程結構物等極具威脅,已引發多起海上事故[1],因此畸形波越來越引起人們的關注,它的發生機理及工程應用問題已成為當前海洋物理學界和船舶水動力學界的一個研究熱點問題[2-5]。海軍艦艇需要在各種海況下航行,海浪環境深刻地影響艦艇的作戰效能,因此如何更加準確地把握遠洋海區風浪環境的非線性規律和特點,對畸形波進行監測和預報,進一步計算艦船遭遇畸形波時的耐波性能、載荷響應以及預報其是否存在傾覆等風險,對完善現有的船舶設計和操縱規范、保障艦艇航行安全以及提升海軍戰斗力具有顯著作用。

關于畸形波的數值模擬方法,從考慮波浪調制不穩定性[6-7]的非線性波動方程(例如三階、四階非線性薛定諤方程)出發,可以研究畸形波的發生機理,但是據此很難把握畸形波發生的時空條件,因此不宜用于開展艦船遭遇畸形波相關的耐波性試驗?;贚onguet-Higgins模型是實驗室模擬產生畸形波的有效手段和常用方法,簡單易用,可以實現畸形波在實驗室的定時定點生成[8]。黃國興[9]采用人工干預組成波隨機初相位的方法,使部分組成波初相位相同,可以得到包含畸形波的波列,但模擬的效率比較低,且不能控制畸形波的生成時間和生成地點。Kriebel[10]采用一個基本隨機波列和一個瞬態波列線性疊加的雙波列疊加模型模擬了畸形波;裴玉國[11]采用三波列疊加模型優化了畸形波的模擬。這兩種模擬方法都可以實現畸形波的定時定點生成,但基于瞬態波列與隨機波列的疊加,瞬態波列的能量所占的比例會影響整個模擬波列的譜的結構。對于大尺度的畸形波,采用瞬態波列方法將影響波浪序列的統計特性。劉贊強[12]采用改進的相位調制方法來模擬畸形波,既滿足波浪序列的統計特性又可保持目標譜的結構,且模擬效率較高,是一種有效的模擬畸形波的方法。因此,本文提出了計及航速航向的相位調制方法來模擬畸形波(頂浪航行時僅需要考慮航速的影響),進而利用CFD(計算流體力學)實現了艦船頂浪狀態下定時定點遭遇畸形波,并與Bennett等人的水池試驗數據對比進行驗證,吻合較好。

1 英國南安普頓大學艦船遭遇畸形波水池實驗

對于畸形波的數值模擬、演化規律及其對近岸結構物響應的研究,國內外學者做了大量的工作,但對于艦船遭遇畸形波的研究還相對較少。國外的研究艦船遭遇畸形波的團隊主要有兩個:德國工業大學的Clauss[13]和英國南安普頓大學的Bennett[14],國內則偏重于理論的研究[15],本文利用Bennett等人的水池實驗結果對數值模擬結果進行驗證,下面對水池試驗相關情況進行簡要介紹。

1.1 水池試驗條件

拖曳水池的長、寬、深分別為60 m×3.7 m×1.86 m,最大可用拖曳速度為4.5 m/s?;尾ň劢沟奈恢迷诰嚯x造波機25 m處,可生成兩維非規則海浪的周期范圍為T≥0.7 s。造波機的槳板是由用戶自定義輸入技術所控制的;獲得的波浪測量結果顯示消波區的反射低于10%。

1.2 船模相關參數

圖1 英國利爾德級護衛艦的型線圖和船模示意圖Fig.1 The body plan and the experimental set-up of leander class frigate hull

水池試驗采用的艦船為英國利安德級護衛艦,曾是20世紀60~80年代英國海軍的主力戰艦。Andrew和Lloyd在1981年時曾用此船模研究過搖蕩運動和甲板上浪。船??s尺比為1:43.62,總長為2.6 m,型線圖和船模示意如圖1所示,主要的船型參數如表1所示。

表1 利安德級護衛艦的主要型值數據Tab.1 Principal particulars of Leander class frigate

1.3 基于相位優化技術生成畸形波

Bennett等人采用了相位優化技術生成畸形波,相位優化技術由Clauss等人[16]提出,波浪的計算機控制優化過程如圖2所示。

給定目標畸形波的聚焦時間、上跨零點周期、最大波高以及峰前波陡,構造約束條件和目標函數;初始的相位是隨機生成的,利用序列二次規劃算法對相位進行修正,波浪序列自動生成、測量、評估和修正,直到滿足目標值。該相位優化技術的優勢是利用物理的水池可以將波浪非線性效果考慮在內,實際上同時完成自我驗證的過程。

圖2 計算機控制的波浪生成優化過程Fig.2 Computer controlled optimization of waves

1.4 Bennett對畸形波的定義

Bennett等人采用的畸形波定義主要由兩個條件組成:

式中:HR為最大波高,ηR為峰值,Hs為有義波高,AI被稱為畸形指數,CI被稱為波峰指數。

2 畸形波海浪環境生成的數值模擬研究

采用CFD方法模擬艦船在波浪中的運動時,假設在固定坐標系中,航速為U0,浪舷角為χ;則長峰非規則波的波高方程為[17]:

設在位置 x=xc、y=yc,時刻 t=tc時生成畸形波,調制初相 θi使部分(或者全部)組成波在 x=xc、y=yc,t=tc時 ηi(xc,yc,tc)為正,則在此疊加的波高會增大。 令組成波數M=M1+M2,則(3)式可以寫為:

在此令后M2個組成波的合成波波面η2(x,y, t)在預定位置處聚焦出現大波,需要調制后M2個組成波的初相位 θi,使 ηi(xc,yc,tc)。

(1) 當 ki(xccos χ+ycsin χ-U0cos χtc)-ωitc<0 時,令整數 N=int[(ki(xccos χ+ycsin χ-U0cos χtc)-ωitc)/2 π ],此時 N<0,(6)式可以寫為:

調制 θi,使-π/2<ki(xccos χ+ycsin χ-U0cos χtc)-ωitc-2Nπ+θi<π/2,這樣 cos( ki(xccos χ+ycsin χ-U0cos χtc)-ωitc-2Nπ+θi)>0,此時 ηi(xc,yc,tc)>0,η2(xc,yc,tc)>0,由于-2π<ki(xccos χ+ycsin χ-U0cos χtc)-ωitc-2Nπ<0,θi在下述區間隨機取值:

(2) 當 ki(xccos χ+ycsin χ-U0cos χtc)-ωitc≥0 時,令整數 N=int[ (ki(xccos χ+ycsin χ-U0cos χtc)-ωitc)/2 π ],此時 N≥0,(6)式可以寫為

調制 θi,使-π/2<ki(xccos χ+ycsin χ-U0cos χtc)-ωitc-(2N+1 ) π+θi<π/2,這樣 cos( xccos χ+ycsin χ-U0cos χtc)-ωitc-(2N+ 1) π+θi>0,此時 ηi(xc,yc,tc)>0,η2(xc,yc,tc)>0,θi的確定方法與情況(1)中所述的相同,在此不再贅述。

艦船頂浪航行時,浪舷角為χ=0;則方程(3)可簡化為:

基于此數學模型,可以對艦船定時定點遭遇畸形波的波浪環境進行數值模擬,為利用CFD方法進行相關研究奠定基礎。

對于CFD而言,采用Bennett等人的相位優化技術還存在一定的技術困難。即使初始條件相同,相位優化解并不依賴于初始隨機相位,優化之后的相位選擇并不唯一,海浪的隨機特性并沒有完全丟失,數值模擬試驗的相位無法與水池實驗的相位完全一致,導致生成的畸形波海浪環境必然存在著差異。因此,盡可能地對相位進行篩選,使各項波浪指數與水池實驗值接近,以期與水池試驗具有可比性。另外,船體周圍的波浪場實際上是由入射波(即畸形波)、艦船反射波和輻射波構成,本節主要考慮畸形波的定時定點生成,在加載船模進行耐波性相關數值計算時,Fluent軟件可以自動將艦船與流場的流固耦合運動考慮在內。

3 低速頂浪航行時數值模擬和水池試驗對比

3.1 基于CFD進行船模水池實驗的方案

利用Fluent前處理軟件Gambit生成的利安德級護衛艦船體如圖3所示。

圖3 利安德級護衛艦船體Fig.3 Leander class frigate hull

Bennett等人通過控制造波機的波浪生成以及拖車的運動,確保船模的中部剛好遭遇畸形波波峰,本文在利用CFD進行相關的數值模擬時,在預定時刻令聚焦位置剛好位于船模中部,從而實現了艦船定時定點遭遇畸形波,并在此記錄波形。

3.1.1 計算域

計算域設置為長方體,如圖4所示,計算域尺寸為:入口距船首1L,出口距船尾2L,頂部邊界距水線0.5L,底部邊界距水線1.5L,左、右邊界距船中縱剖面1L。

圖4 頂浪流場的計算域Fig.4 The computational domain of head-sea filed

3.1.2 網格劃分

為充分發揮結構及非結構網格優勢,對計算域分為船體近流場區及遠流場區進行處理,近流區域采用非結構網格進行網格生成,船體面網格采用三角形網格,從船體表面網格以一定比例外推而生成近密外疏的非結構四面體網格。遠流場區完全布設結構網格,由自由面向上下兩底邊漸疏,在保證網格質量的同時保持合適的網格數量,整個計算域網格劃分效果如圖5所示。

網格總數為1 381 643,網格質量分布如表2所示。

3.1.3 邊界條件設置

計算域及邊界條件名稱如圖6所示。

具體設置如下:

入口邊界(in):速度入口(velocity-inlet)條件,設置入口處流體的速度、入口處水的體積分數、湍動能k,耗散率ε;

圖5 網格劃分效果Fig.5 Meshing effect

表2 網格質量評價Tab.2 The evaluation of mesh quality

出口邊界(out):壓力出口(pressure-outlet)條件,設定靜壓力、底面位置和自由面高度,回流的湍動能k及耗散率ε則采用Fluent的推薦值;

圖6 計算域及邊界名稱Fig.6 Computational domain and name of boundary conditions

上下邊界(top、bot)—速度入口,給定三個方向流速(u=U0、v=w=0)及水的體積分數(0);

左右邊界(port、stab)—滑移的壁面;

船體(ship)—有剪切力無滑移的壁面。

3.1.4 數值模擬的初始條件

數值波浪水池采用的靶譜為Jonswap譜,其水池試驗參數如表3所示。

表3 水池試驗的主要參數Tab.3 The main parameters in experiment

聚焦時刻為tc=5 s,聚焦位置位于船舶的重心處即xc=0 m。

3.2 數值模擬結果及分析

3.2.1 利安德護衛艦頂浪航行流場分析

Bennett等人利用相位優化技術在水池中定時定點生成了畸形波,為檢驗數據采集系統的不確定度,特進行了三次重復實驗,其波面的測量結果如圖7所示。

圖7 測量過程中波面的不確定度Fig.7 Wave profiles showing uncertainty in experimental measurements

由圖 7可知,畸形波的聚焦時間為tc=6.2 s,AI接近2.36,CI接近1.37,說明利用相位優化技術生成可以定時定點生成畸形波,另外,測量過程的相對誤差較小。

為了分析艦船反射波、輻射波對生成的畸形波波浪場的影響,不加載艦船的波浪場中,在船模中心位置設置浪高儀進行時歷監測,并作為參考值。在圖5的三維計算域中,加載了船模之后,分別在距離船模中心橫向距離y=0.5 m、1.5 m、2.5 m處設置浪高儀,進行波浪時歷監測,并將測量結果同參考值進行比對,對圖8所包含的極值大波進行特征統計,其畸形特征參數如表4所示。

由圖 8和表4可知,加載了船舶搖蕩之后,對入射波浪場(畸形波)產生了干擾影響,與參考值對比可知,這種干擾影響相對較小。干擾主要是由船舶在行進過程中產生的反射波和繞射波造成的,隨著遠離船模,干擾作用在不斷減小。

圖8 不同橫向位置處測得的波浪時歷Fig.8 Wave profiles measured in different transverse locations

表4 不同橫向位置處測得波浪時歷的畸形參數比對Tab.4 The freak-wave parameter of wave profiles measured in different transverse locations

3.2.2 利安德護衛艦頂浪航行搖蕩運動數值模擬

當船舶流場發展比較充分的時候,取船模搖蕩運動比較穩定后的某一時刻(t=1.5 s)作為記錄起點,數值模擬船模在包含畸形波的非規則波浪中做縱搖和垂蕩運動,當艦船遭遇畸形波時的時刻(t=4.5~6 s,每隔0.25 s取一幅圖片)的瞬時態勢如圖 9所示。

圖9 U0=0.6 m/s時艦船頂浪遭遇畸形波的過程Fig.9 The process of ship encountering with freak waves at forward speed U0=0.6 m/s

由圖 9可知,艦船在包含畸形波的非規則波浪中頂浪航行時,當遭遇的波浪比較大時,會出現砰擊、上浪等強非線性現象,會對船舶產生一系列不利于航行的搖蕩響應,甚至危及船舶的航行安全。

圖9中利安德級護衛艦船模在包含畸形波的非規則波浪航行時,測得的縱搖時歷和垂蕩時歷分別如圖10和圖11所示。

由圖10和圖11可知,當艦船遭遇比較大的波浪時,垂蕩和縱搖值相對較大,特別遭遇畸形波時,更加明顯,與圖9對應可知,產生了砰擊和上浪現象。

采用下跨零點法對U0=0.6 m/s工況下所采集的縱搖時歷和垂蕩時歷進行處理,求得其最大縱搖值和最大垂蕩值(與水池實驗一致,此處采用的是谷—峰最大值)如表5所示。

圖10 U0=0.6 m/s工況下所采集到的縱搖時歷Fig.10 The pitch profile obtained at forward speed U0=0.6 m/s

圖11 U0=0.6 m/s工況下所采集到的垂蕩時歷Fig.11 The heave profile obtained at forward speed U0=0.6 m/s

表5 低速頂浪航態下數值模擬值與水池實驗值的比對Tab.5 The comparison between CFD numerical simulation values and wave-tank experimental values at forward low-speed

由表 5可知,CFD數值模擬的初始條件同水池試驗基本一致,則獲取的ξmax(m)和θmax(deg)比較接近,誤差值可以接受,證明了低速頂浪航行情況下本文數值模擬的有效性。

4 中速頂浪航行數值模擬和水池試驗對比

當艦船以不同的航速遭遇畸形波時,其運動響應又會有何不同,結合Bennett等人的水池實驗情況,本節設計了另一航速的試驗,試驗的主要參數如表6所示。

表6 試驗的主要參數Tab.6 The main parameters in experiment

在圖5的三維計算域中,距離船模中心橫向距離y=2.5 m處設置浪高儀進行波浪時歷監測,如圖12所示。

圖12 橫向位置y=2.5 m處測得的波浪時歷Fig.12 Wave profiles measured in transverse location y=2.5 m

利用下跨零點法對圖12中的波浪時歷所包含的極值大波進行特征統計,將畸形參數與水池實驗給定值進行比對,如所表7示。

表7 水池試驗和CFD數值模擬參數的比對Tab.7 The comparison between CFD numerical simulation and wave-tank experimental values

在U0=1.4 m/s的工況下,艦船遭遇包含畸形波的非規則海浪所對應的縱搖和垂蕩時歷分別如圖13和圖14所示。

圖13 U0=1.4 m/s工況下所采集到的縱搖時歷Fig.13 The pitch profile obtained at forward speed U0=1.4 m/s

圖14 U0=1.4 m/s工況下所采集到的垂蕩時歷Fig.14 The heave profile obtained at forward speed U0=1.4 m/s

采用下跨零點法對U0=1.4 m/s工況下所采集的縱搖時歷和垂蕩時歷進行畸形參數統計,求得其最大縱搖值和最大垂蕩值(與水池實驗一致,此處采用的是谷-峰最大值)如表8所示。

表8 中速頂浪航態下CFD數值模擬值與水池實驗值的比對Tab.8 The comparison between CFD numerical simulation and wave-tank experimental values at forward middle-speed

由表 8可知,CFD數值模擬的初始條件同水池試驗一致,則獲取的ξmax(m)和θmax(deg)比較接近,誤差值可以接受,證明了中速頂浪航行情況下本文數值模擬的有效性。

將表5和表8對比可知,隨著航速的增大,垂蕩值增大,縱搖值減少,CFD數值模擬得出的趨勢與水池實驗的趨勢相同,得出這樣的結論與艦船在畸形波中的搖蕩運動有關。當艦船遭遇畸形波時,船頭有可能被掩埋在波高較大的海浪中,當船從波峰中穿過時會出現明顯的上浪現象。當高速遭遇到海浪時,船頭掩埋和甲板上浪會更加嚴重;在一定程度上使縱搖運動減弱;然而,垂蕩運動并不因為上浪嚴重而受到影響,而是隨著航速的增大而增大。

5 試驗結果與船舶設計規范的對比

由于畸形波的突發特性,可以在毫無預兆的情況下出現,意味著艦船有隨時遭遇畸形波的可能,特別是在海況比較惡劣的情況下,導致海難事故時有發生,因此在船舶設計階段,就需要將最為惡劣的情況考慮在內,以防悲劇的發生。

Bennett等人對水池試驗的搖蕩運動時歷使用中心差分方法得到垂蕩加速度和縱搖加速度,并與英國勞安德船級社的規范值進行比對,發現縱搖加速度遠小于規范值,而垂蕩加速度則有可能大于船級社的規范值(特別是航速較高的情況下)。在利用CFD進行數值模擬時,U0=0.6 m/s和U0=1.4 m/s兩種工況下的垂蕩加速度如圖15所示。U=1.4m/s U=0.6m/s

圖15 兩種工況下的垂蕩加速度對比Fig.15 The comparison of the heave acceleration under two different work conditions

求取圖15中垂蕩加速度曲線的最大值,并將其與重力加速度進行歸一化,進而與英國勞安德船級社的規范值(2010)對比,如表9所示。

由表 9 可知,當 U0=0.6 m/s時,amax(g)<0.39,并沒有超越船級社的規范值,但是當U0=1.4 m/s時,amax(g)>0.39,超越了船級社的規范值,這在船舶設計的過程中需要予以考慮。

表9 試驗得到的垂蕩加速度與船級社規范值的比對Tab.9 The comparison between heave acceleration value with Lioyd’s Register Rules

6 結 論

本文基于隨機波浪的Longuet-Higgins模型,在相位調制方法的基礎上考慮了航速航向的影響,調制部分組成波的初相位,可以實現艦船定時定點遭遇畸形波。主要可以得出以下結論:

(1)利用CFD方法對頂浪情況下艦船定時定點遭遇畸形波進行了數值模擬,并將數值模擬結果與英國南安普頓水池實驗值進行比對,吻合較好,證明了本文數值模擬的有效性。

(2)當艦船遭遇相同的畸形波時,隨著航速的增大,垂蕩值增大,縱搖值減小,CFD數值模擬得出的趨勢與水池試驗的結論一致。

(3)當艦船高速遭遇畸形波時,某些動力響應值(比如垂蕩加速度)超出了船級社的規范值,有可能存在一些安全隱患,在船舶設計過程中是應該考慮在內的。

參 考 文 獻:

[1]董艷秋,紀 凱,黃衍順.波浪中船舶橫搖穩性的研究[J].船舶力學,1999,3(2):1-6.Dong Yanqiu,Ji Kai,Huang Yanshun.Study on the ship stability in rolling in waves[J].Journal of Ship Mechanics,1999,3(2):1-6.

[2]劉首華,牟 林,劉克修,等.畸形波研究的進展及存在問題[J].地球科學進展,2013,28(6):665-673.Liu Shouhua,Mu Lin,Liu Kexiu,et al.Research advance and problems in freak waves[J].Advances in Earth Science,2013,28(6):665-673.

[3]高志一,于福江,徐福祥,等.畸形波生成條件預報方法研究進展[J].海洋通報,2011,30(3):351-356.Gao Zhiyi,Yü Fujiang,Xü Fuxiang,et al.Progress in the method of freak wave condition forecasting[J].Martine Science Bulletin,2011,30(3):351-356.

[4]Cui Cheng,Zhang Ningchuan,Pei Yuguo,Liu Qianlie.Numerical study on generation and evolution of freak waves[J].Journal of Ship Mechanics,2012,12,16(12):51-56.

[5]沈玉稿.畸形波的數值模擬及其與海洋結構物相互作用研究[D].上海:上海交通大學,2013.Shen Yügao.Numerical simulation of freak wave and the interactions between freak wave and offshore structures[D].Shanghai:Shanghai Jiao Tong University,2013.

[6]Zhang Yunqiu,Zhang Ningchuan.Numerical simulation and mechanism analysis of freak waves[J].Acta Oceanologica Sinica,2007,26(5):116-124.

[7]Lo E,Mei C C.A numerical study of water-wave modulation based on a higher-order nonlinear Schroedinger equation[J].Journal of Fluid Mechanics,1985,150:395-416.

[8]劉贊強.畸形波模擬與其在核電取水結構物作用探討[D].大連:大連理工大學,2011.Liu Zanqiang.Freak wave simulation and its interaction with nuclear power station ingarage structure[D].Dalian:Dalian University of Technology,2011.

[9]黃國興.畸形波的模擬方法及基本特性研究[D].大連:大連理工大學,2002.Huang Guoxing.Research on the simulation method and basic characteristics of freak waves[D].Dalian:Dalian University of Technology,2002.

[10]Kriebel D L.Efficient simulation of extreme waves in a random sea[C]//Abstract for Rogue Waves 2000 Workshops.Brest,France,2000:1-2.

[11]裴玉國.畸形波的生成及基本特性研究[D].大連:大連理工大學,2007.Pei Yüguo.The generation of freak waves and its behaviors[D].Dalian:Dalian University of Technology,2007.

[12]劉贊強,張寧川,俞聿修.改進的相位調制法模擬畸形波:I-理論模型與驗證[J].水動力研究與進展,A輯,2010,25(3):383-390.Liu Zanqiang,Zhang Ningchuan,Yü Yüxiu.The generation of freak waves based on a modified phase modulation method:I-Theory and validation[J].Chinses Journal of Hydrodynamics,Series A,2010,25(3):383-390.

[13]Clauss G F,M Klein,Kauffeldt A.Limiting loads and motions of ships in extreme sea states[C]//13th Congress of Int1.Maritime Assoc.of Mediterranean IMAM 2009,12-15 Oct.2009.Istanbul,Turkey,2009.

[14]Bennett S S,Hudson D A,Temarel P.The influence of forward speed on ship motions in abnormal waves:Experimental measurements and numerical predictions[J].Journal of Fluids and Structures,2013,154-172.

[15]董艷秋,楊冠聲,陳學闖.關于一種危害船舶安全的波浪—畸形波的探討[J].船舶力學,2003,7(2):33-38.Dong Yanqiu,Yang Guansheng,Chen Xuechuang.Study on the great damage caused by freak wave[J].Journal of Ship Mechanics,2003,7(2):33-38.

[16]Gunther F G,Christian E S,Marco Klein.Generation of rogue waves with predefined steepness[C]//Proceedings of OMAE 2006,June 4-9,2006.Hamburg,Germany,2006.

[17]吳 明.不規則波中艦船搖蕩運動的數值模擬及預報研究[D].大連:海軍大連艦艇學院,2013 Wu ming.The research on prediction and simulation of ship motions in irregular waves[D].Dalian:Dalian Naval Academy,2013.

主站蜘蛛池模板: www.亚洲色图.com| 久久公开视频| 精品国产福利在线| 在线视频97| 国产九九精品视频| 久久青草视频| 日韩成人在线视频| 国产性生交xxxxx免费| 国产一区二区免费播放| 中国国语毛片免费观看视频| 亚洲人免费视频| 午夜精品久久久久久久无码软件 | 国产a在视频线精品视频下载| 香蕉精品在线| 69国产精品视频免费| 天天躁狠狠躁| 无码专区第一页| 亚洲免费黄色网| 日本a级免费| 91最新精品视频发布页| 亚洲无码视频图片| 少妇精品在线| 国产91全国探花系列在线播放| 日韩AV手机在线观看蜜芽| 伊人大杳蕉中文无码| 专干老肥熟女视频网站| 欧美一级在线播放| 免费日韩在线视频| 婷婷六月在线| 中文字幕无码电影| 亚洲女同一区二区| 午夜啪啪福利| 伊人91在线| 亚洲精品成人福利在线电影| 小蝌蚪亚洲精品国产| 波多野结衣视频网站| 亚洲黄色片免费看| 国产91丝袜| 国产在线91在线电影| 欧美一级在线看| 伊人福利视频| 精品国产aⅴ一区二区三区| 免费高清a毛片| 亚洲天堂视频网| 伊人久久大线影院首页| 国产va在线| 亚洲精选高清无码| 97国产成人无码精品久久久| 国产91在线|日本| 亚洲精品视频免费看| 国产亚洲欧美在线专区| 亚洲码一区二区三区| 欧美视频在线观看第一页| 亚洲成人手机在线| 国产精品伦视频观看免费| 免费三A级毛片视频| 亚洲综合香蕉| 国产亚洲日韩av在线| 免费一级α片在线观看| 国产肉感大码AV无码| 无码一区二区波多野结衣播放搜索| 激情无码字幕综合| 久久黄色毛片| 国产99视频在线| 亚洲男人天堂网址| 97国产在线视频| 国产十八禁在线观看免费| 国产aaaaa一级毛片| 国产精品林美惠子在线播放| 五月天综合婷婷| 熟女视频91| 国产精品免费久久久久影院无码| 国产小视频免费观看| 中文字幕有乳无码| 成人福利在线免费观看| a毛片在线播放| 人人艹人人爽| 国产视频欧美| 婷婷午夜影院| 国产精品久久久久久搜索| 久久一级电影| 国产在线观看99|