亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        艦船頂浪遭遇畸形波的波浪水池?cái)?shù)值模擬研究

        2016-07-16 00:16:42張永勝?gòu)埍据x蔡烽王驍薛亞?wèn)|李
        船舶 2016年2期
        關(guān)鍵詞:數(shù)值模擬

        張永勝?gòu)埍据x蔡 烽王 驍薛亞?wèn)|李 東

        (1.海軍大連艦艇學(xué)院 訓(xùn)練部 大連116018; 2.海軍大連艦艇學(xué)院 航海系 大連116018)

        ?

        艦船頂浪遭遇畸形波的波浪水池?cái)?shù)值模擬研究

        張永勝1張本輝2蔡 烽2王 驍2薛亞?wèn)|2李 東2

        (1.海軍大連艦艇學(xué)院 訓(xùn)練部 大連116018; 2.海軍大連艦艇學(xué)院 航海系 大連116018)

        [摘 要]畸形波對(duì)船舶航運(yùn)極具威脅,已引發(fā)多起海上事故, 因此有必要對(duì)艦船遭遇畸形波進(jìn)行數(shù)值模擬。基于長(zhǎng)峰波隨機(jī)海浪模型,在相位調(diào)制生成畸形波的同時(shí)考慮了艦船航速的影響,可以對(duì)艦船在頂浪航行定時(shí)定點(diǎn)遭遇畸形波的非線性波浪環(huán)境進(jìn)行建模,利用CFD生成數(shù)值波浪水池,并對(duì)畸形波生成過(guò)程中畸形參數(shù)、峰度、陡度的變化以及速度分量隨深度的變化進(jìn)行分析,從一定程度上反映畸形波的演化特性,并驗(yàn)證數(shù)值模擬的有效性。

        [關(guān)鍵詞]畸形波;數(shù)值波浪水池;數(shù)值模擬;相位調(diào)制

        張本輝(1988-),男,博士,研究方向:艦船操縱性與耐波性。

        蔡 烽(1973-),男,博士后,副教授,研究方向:非線性海浪。

        王 驍(1980-),男,博士,講師,研究方向:艦船操縱性。

        薛亞?wèn)|(1990-),男,碩士,研究方向:非線性海浪。

        李 東(1990-),男,碩士,研究方向:艦船耐波性。

        引 言

        畸形波是海洋中高且陡的大波,其持續(xù)時(shí)間很短,但出現(xiàn)的偶然性和巨大的破壞性,對(duì)船舶航運(yùn)和海洋工程結(jié)構(gòu)物等極具威脅[1]。因此越來(lái)越引起人們關(guān)注,其發(fā)生機(jī)理及工程應(yīng)用問(wèn)題已成為當(dāng)前物理海洋學(xué)界和船舶水動(dòng)力學(xué)界的一個(gè)研究熱點(diǎn)。越來(lái)越多的記錄中提到艦船遭遇畸形波的現(xiàn)象,對(duì)這樣的情形進(jìn)行建模研究是非常必要的,進(jìn)而合理地預(yù)報(bào)艦船的極限載荷與響應(yīng),旨在評(píng)價(jià)艦船遭遇畸形波的風(fēng)險(xiǎn)情況以及完善現(xiàn)有的船舶設(shè)計(jì)和操縱規(guī)范、進(jìn)而保障艦艇航行安全,其中較為關(guān)鍵的一點(diǎn)就是生成艦船遭遇畸形波的(數(shù)值)波浪環(huán)境。

        關(guān)于畸形波的數(shù)值模擬方法,從考慮波浪調(diào)制不穩(wěn)定性[2]的非線性波動(dòng)方程(例如三階、四階非線性薛定諤方程)出發(fā),可以研究畸形波的發(fā)生機(jī)理,但是很難控制畸形波發(fā)生的時(shí)空條件(即需要等待畸形波的出現(xiàn)且不能預(yù)測(cè)其出現(xiàn)的位置),因此不適用于開(kāi)展艦船遭遇畸形波相關(guān)的耐波性試驗(yàn)。為研究艦船遭遇畸形波時(shí)的響應(yīng),則定時(shí)定點(diǎn)生成預(yù)定波高的畸形波就非常有必要。劉贊強(qiáng)等[3]采用基于線性疊加原理,調(diào)制子波的初相位,使子波在預(yù)定位置預(yù)定時(shí)間進(jìn)行聚焦來(lái)實(shí)現(xiàn)畸形波的數(shù)值模擬,既滿足波浪序列的統(tǒng)計(jì)特性又可保持目標(biāo)譜的結(jié)構(gòu),且需要的子波個(gè)數(shù)相對(duì)較少,模擬效率較高。國(guó)內(nèi)外學(xué)者對(duì)于畸形波對(duì)于近岸結(jié)構(gòu)物響應(yīng)的研究做了大量的工作[4-6],但對(duì)于艦船遭遇畸形波的研究還相對(duì)較少。本文在劉贊強(qiáng)相位調(diào)制方法的基礎(chǔ)上,考慮了航速的影響,可以對(duì)艦船頂浪航行定時(shí)定點(diǎn)遭遇畸形波進(jìn)行建模。

        目前實(shí)際海上艦船遭遇畸形波的相關(guān)數(shù)據(jù)非常少且不易獲取,更可能由于海況惡劣,提取的數(shù)據(jù)誤差較大?;贑FD(computational fluid dynamics,計(jì)算流體力學(xué))方法能有效彌補(bǔ)海上試驗(yàn)、水池試驗(yàn)及勢(shì)流理論等方面的不足,利用CFD可以實(shí)現(xiàn)艦船在規(guī)則波或者非規(guī)則波浪中搖蕩運(yùn)動(dòng)的數(shù)值模擬及預(yù)報(bào),能夠計(jì)算波浪破碎、甲板上浪等強(qiáng)非線性問(wèn)題,在船舶水動(dòng)力學(xué)研究領(lǐng)域取得了豐碩的研究成果,因此可以利用CFD對(duì)艦船頂浪遭遇畸形波進(jìn)行數(shù)值模擬。

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

        當(dāng)模擬艦船遭遇畸形波時(shí)候,需要盡可能使模擬的海況接近于真實(shí),海上興起的實(shí)際波浪通常是短峰非規(guī)則波。Tromans等人[7]認(rèn)為畸形波海況可以近似為長(zhǎng)峰波,F(xiàn)aulkner和Williams[8]將其描述為“水墻”。因此,本文基于長(zhǎng)峰不規(guī)則波海浪模型來(lái)進(jìn)行畸形波的數(shù)值建模,某一固定的波面方程可由無(wú)數(shù)隨機(jī)的余弦波線性疊加予以描述。

        采用CFD方法模擬艦船在波浪中的運(yùn)動(dòng)時(shí),在相對(duì)坐標(biāo)系中,航速為U,頂浪航行時(shí)為正,順浪航行時(shí)為負(fù),非規(guī)則波的波高方程[9]為:

        三個(gè)方向的速度方程為:

        設(shè)在位置x=xc、時(shí)刻t=tc時(shí)生成畸形波,調(diào)制θi,使部分或者全部組成波在x=xc、t=tc時(shí)為正,在此疊加的波高會(huì)增大。令組成波數(shù)M = M1+ M2,可以寫(xiě)為:

        在此令后M2個(gè)組成波的合成波波面在預(yù)定位置處聚焦出現(xiàn)大波,需要調(diào)制后M2個(gè)組成波的初相位θi,使。

        調(diào)制θi使這樣此時(shí)則由于, θi在下述區(qū)間隨機(jī)取值:

        調(diào)制θi使,這樣,此時(shí)θi的確定方法與情況(1)中所述的相同,在此不再贅述。

        2 數(shù)值仿真及結(jié)論

        2.1仿真條件

        目標(biāo)譜為Jonswap譜,有義波高Hs= 0.1 m (縮尺比為1∶43.62,則對(duì)應(yīng)于真實(shí)海況的4.362 m),譜峰升高因子γ = 3.3,譜峰周期Tp= 1.45 s(對(duì)應(yīng)全尺度的譜峰周期為9.48 s),頻率范圍為f =0.5 fp~3.5 fp;組成波數(shù)M=70,由于從高頻向低頻調(diào)制優(yōu)于低頻向高頻調(diào)制[10],故調(diào)制70個(gè)組成波序列中后60個(gè)組成波的初相位。

        2.2計(jì)算域與網(wǎng)格劃分

        根據(jù)流場(chǎng)的有關(guān)特征,建立的三維計(jì)算域范圍為12 m×5.2 m×4.5 m,其中水深為3.9 m,將其劃分為造波區(qū)(前9 m)、消波區(qū)(后3 m),暫不考慮船模的影響,采用全結(jié)構(gòu)網(wǎng)格對(duì)計(jì)算域進(jìn)行劃分,網(wǎng)格節(jié)點(diǎn)總數(shù)為1 065 400。參考文獻(xiàn)9中有關(guān)CFD計(jì)算網(wǎng)格無(wú)關(guān)性的實(shí)驗(yàn)結(jié)論,兼顧效率與精度,網(wǎng)格布設(shè)如圖1所示。

        圖1 水池自由面網(wǎng)格加密示意圖

        本文以Flunt軟件為試驗(yàn)平臺(tái),以不可壓縮粘性流體的N-S方程為控制方程,利用軟件的二次開(kāi)發(fā)功能,采用邊界造波法和阻尼消波法生成數(shù)值波浪水池。在y = 0 m的縱剖面上設(shè)置9個(gè)浪高儀進(jìn)行時(shí)歷監(jiān)測(cè),浪高儀沿x方向的坐標(biāo)位置分別為 x=2.5 m、3 m、3.5 m、4 m、4.5 m、5 m、5.5 m、6 m、6.5 m;另外對(duì)整個(gè)波面演化進(jìn)行實(shí)時(shí)錄像,以分析畸形波生成過(guò)程中的演化情況。計(jì)算時(shí)間步長(zhǎng)為0.001 s,采樣頻率為100 Hz,假設(shè)生成畸形波的預(yù)定時(shí)間和預(yù)定位置分別設(shè)為tc= 10 s,xc= 4.5 m,當(dāng)U = 1.4 m/s時(shí)(對(duì)應(yīng)于實(shí)際航速為9.25 m/s,18 kn),數(shù)值波浪水池在t=10 s時(shí)的瞬時(shí)波面如下頁(yè)圖2所示。

        從圖2中可以看出,在預(yù)定時(shí)刻預(yù)定位置出現(xiàn)了極值大波,此時(shí)波峰值為0.142,遠(yuǎn)大于有義波高0.1 m。x = 4.5 m處浪高儀監(jiān)測(cè)的波浪時(shí)歷如圖3所示。

        根據(jù)Klinting和Sand[11]對(duì)畸形波的定義,即畸形波的波高Hj應(yīng)該滿足以下條件其中 Hj-1和Hj+1是畸形波前后相鄰波浪的波高,ηj是畸形波對(duì)應(yīng)的波峰高度(相對(duì)于水平線)。在此將α1、α2、α3、α4統(tǒng)稱為畸形波特征參數(shù)。利用上跨零點(diǎn)法對(duì)波浪時(shí)歷進(jìn)行特征統(tǒng)計(jì),則圖3中極值波的特征參數(shù)分別為2.1、2.4、2.03、0.67,上述指標(biāo)全部符合畸形波的定義,證明了本文所建模型的有效性。

        圖2 t=10 s時(shí)刻的瞬時(shí)波面

        圖3 x = 4.5 m監(jiān)測(cè)的波浪時(shí)歷

        2.3畸形波海浪的非線性演化分析

        畸形波多在未知和不可預(yù)測(cè)的條件下發(fā)生,實(shí)測(cè)資料十分缺乏,且所有的資料都是海上或者海岸某一固定點(diǎn)監(jiān)測(cè)到的波面時(shí)歷,沒(méi)有艦船移動(dòng)坐標(biāo)系下記錄到的完整畸形波發(fā)展過(guò)程的空間記錄,本文通過(guò)設(shè)置的9個(gè)浪高儀,可以監(jiān)測(cè)畸形波生成、發(fā)展演變過(guò)程,進(jìn)而探討艦船移動(dòng)坐標(biāo)系下畸形波的非線性特征。無(wú)調(diào)制的非規(guī)則波浪為正態(tài)隨機(jī)過(guò)程,受非線性的影響,其波面過(guò)程分布與高斯分布有些偏差,其偏差的程度可以用偏度(Skewness)和峰度(kurtosis)來(lái)衡量,偏度是不對(duì)稱性的綜合度量,峰度表示的是波面高度頻率分布曲線峰的尖銳程度。對(duì)于高斯分布,偏度為0,峰度為3.0。各個(gè)監(jiān)測(cè)點(diǎn)的畸形參數(shù)、最大峰高、偏度、峰度以及統(tǒng)計(jì)得到的有義波高如表1所示。

        表1 各個(gè)監(jiān)測(cè)點(diǎn)的統(tǒng)計(jì)特征參數(shù)

        由表1可知,只有x = 4.5 m的波浪時(shí)歷中所含的極值波嚴(yán)格地滿足了畸形波的要求,說(shuō)明畸形波的空間時(shí)間持續(xù)時(shí)間很短,突然出現(xiàn)又很快消失。在x = 2.5~5.5 m內(nèi)α1變化相對(duì)較?。碒j變化較?。?,基本上都在2以上,即使畸形波存在一個(gè)發(fā)展和消亡的過(guò)程,但是在某個(gè)空間時(shí)間范圍內(nèi)其波高相對(duì)較大。x=2.5~6.5 m范圍內(nèi),總體上而言α2基本上逐漸增大,而α3逐漸減小,Hj-1逐漸減少,而Hj+1逐漸增大,通過(guò)錄像可以看出,畸形波峰之前的波峰迅速減少,波谷在演化過(guò)程中由深谷迅速變淺,而畸形波峰之后的波谷逐漸變深,且后一波峰也逐漸隆起。α4表示畸形波的垂直不對(duì)稱性,在畸形波發(fā)展過(guò)程中,α4逐漸增大,說(shuō)明其波峰不斷增大,即子波能量不斷聚集,之后又逐漸消散。由畸形參數(shù)的變化可知,畸形波生成演化的整個(gè)過(guò)程中能量變化非常快。

        由偏度和峰度可知,在畸形波預(yù)定的生成位置附近,其值偏離高斯分布的程度非常大,然而x = 4.5 m處的偏度和峰度并非最大值,Matlab理論計(jì)算的結(jié)果為偏度和峰度在x = 4.5 m處偏離正態(tài)分布程度最大,與CFD數(shù)值模擬的結(jié)論略有不同。后者考慮了水的粘性、波波相互作用等,更加接近于海上實(shí)際,結(jié)果的可信度更高一些。有資料認(rèn)為[12]畸形波破壞最大時(shí)刻不是其波高最大的時(shí)刻,而在發(fā)展過(guò)程中的向前沖擊力更大,對(duì)結(jié)構(gòu)物的破壞可能更為嚴(yán)重,當(dāng)波高最大時(shí)刻應(yīng)該是“強(qiáng)弩之末”,對(duì)船舶或者結(jié)構(gòu)物的沖擊力變小,這些結(jié)果還有待船舶搖蕩試驗(yàn)的進(jìn)一步驗(yàn)證。

        2.4速度場(chǎng)分析

        波面的傳播速度是波浪最基本和重要的特征參數(shù)之一,利用Fluent軟件,可直接捕捉自由面的速度。圖2中自由面對(duì)應(yīng)的速度場(chǎng)如圖4所示。

        圖4 圖2中瞬時(shí)波面對(duì)應(yīng)的速度場(chǎng)

        從圖4可以看出,在畸形波波峰處,畸形波波峰處的波速比其它位置大很多,而波谷處的波速則相對(duì)較小,注意此處的速度場(chǎng)包含著艦船的相對(duì)航速。在指定的聚焦位置xc= 4.5 m和聚焦時(shí)刻tc= 10 s,水質(zhì)點(diǎn)沿三個(gè)方向的速度分量隨水深的變化如圖5所示。

        圖5 水質(zhì)點(diǎn)三個(gè)方向的速度隨水深的變化曲線

        從圖5可見(jiàn),沿y方向的速度分量v的量級(jí)為10-5,相對(duì)于x方向的速度分量u和z方向的速度分量w則可忽略不計(jì)。隨著深度的增加,u呈負(fù)梯度由最大值2.79 m/s迅速減少至航速;而w首先呈正梯度增大至最大值0.16 m/s,然后呈負(fù)梯度迅速減少至0。對(duì)于總的速度幅值而言,z方向的速度分量的平方之后的量級(jí)僅為x方向的速度分量的百分之一,因此在畸形波附近位置,整個(gè)波面猶如一堵向前運(yùn)動(dòng)的水墻,其瞬時(shí)能量非常集中,可能對(duì)船舶的結(jié)構(gòu)強(qiáng)度產(chǎn)生巨大沖擊力,具有巨大的破壞性,會(huì)對(duì)艦船造成致命的傷害。

        3 結(jié) 論

        本文實(shí)現(xiàn)了艦船頂浪情況下遭遇畸形波波浪水池的CFD模擬,可得出以下結(jié)論:

        (1)本文所建立模型,可有效模擬艦船頂浪定時(shí)定點(diǎn)遭遇畸形波的波浪環(huán)境,且可保持波浪序列的統(tǒng)計(jì)特性。

        (2)畸形波生成演化的整個(gè)過(guò)程中能量變化非??欤陬A(yù)定的位置附近畸形波生成時(shí)刻波浪的非線性最強(qiáng),然而偏度和峰度值的最大偏移點(diǎn)并非在畸形波發(fā)生時(shí)波峰最大的時(shí)刻,在畸形波發(fā)展過(guò)程中其向前的沖擊力可能會(huì)更大。

        (3)在畸形波波峰處的波速相對(duì)較大,x方向的速度分量對(duì)總速度幅值起主導(dǎo)作用。

        另外,波浪場(chǎng)穩(wěn)定后,載入船模進(jìn)行搖蕩試驗(yàn)時(shí),船模在網(wǎng)格中的位置不變,可以保證艦船在移動(dòng)坐標(biāo)系下定時(shí)定點(diǎn)遭遇畸形波,而船模繞固定點(diǎn)進(jìn)行六自由度搖蕩運(yùn)動(dòng)。此時(shí),則需考慮船模繞射波和輻射波對(duì)畸形波的干擾作用。艦船以某種航速頂浪遭遇畸形波時(shí)測(cè)得最大縱搖與垂蕩值以及船體某些部位的垂直彎矩,并與水池?cái)?shù)據(jù)進(jìn)行對(duì)比,能進(jìn)一步對(duì)本模型進(jìn)行驗(yàn)證。

        [參考文獻(xiàn)]

        [ 1 ] 張運(yùn)秋,張寧川,裴玉國(guó). 畸形波數(shù)值模擬的一個(gè)有效模型[J].大連理工大學(xué)學(xué)報(bào),2008 (3):406-410.

        [ 2 ] Lo E, Mei C C.A numerical study of water-wave modulation based on a higher-order nonlinear Schr?edinger equation [J]. Journal of Fluid Mechanics,1985(1):395-416.

        [ 3 ] 劉贊強(qiáng),張寧川,俞聿修. 改進(jìn)的相位調(diào)制法模擬畸形波:I -理論模型與驗(yàn)證[J].水動(dòng)力研究與進(jìn)展,A輯,2010 (3):383-390.

        [ 4 ] 沈玉稿. 畸形波的數(shù)值模擬及其與海洋結(jié)構(gòu)物相互作用研究[D]. 上海:上海交通大學(xué),2013.

        [ 5 ] 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,Istanbul,2009:12-15.

        [ 6 ] 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(5):154-172.

        [ 7 ] Tromans P S,AnaturkA R,Hagemeijer P. A new model for kinematics of large ocean waves-application as a design wave[C]. ISOPE,1991:64-71.

        [ 8 ] Faulkner D,Williams R A, Design for abnormal ocean waves[J]. Structural Mechanics,1997:187-193.

        [ 9 ] 吳明.不規(guī)則波中艦船搖蕩運(yùn)動(dòng)的數(shù)值模擬及預(yù)報(bào)研究[D].大連:海軍大連艦艇學(xué)院,2013.

        [10] 劉贊強(qiáng),張寧川,俞聿修.改進(jìn)的相位調(diào)制法模擬畸形波:II-畸形波特征參數(shù)和模擬效率的影響因素探討[J].水動(dòng)力研究與進(jìn)展(A輯),2010 (6):813-821.

        [11] Klinting P,Sand S. Analysis of Prototype Freak Waves[C]. Coastal Hydrodynamics,1987:618-632.

        [12] 俞聿修,桂滿海.海浪的群性及其主要特征參數(shù)[J].海洋工程,1998 (3):9-21.

        Wave tank numerical simulation of ships in head sea encountering freak waves

        ZHANG Yong-sheng1ZHANG Ben-hui2CAI feng2WANG Xiao2XUE Ya-dong2LI Dong2
        (1. Department of Training, Dalian Naval Academy, Dalian 116018, China; 2. Department of Navigation, Dalian Naval Academy, Dalian 116018, China)

        Abstract:The freak wave is a great threat to shipping, and has caused many marine accidents. It is necessary to carry out the numerical simulation of ships encountering freak waves. Based on the long-crested (unidirectional) random wave model, freak wave is generated by the phase modulation synchronously considering the influence of ship speed. Then the nonlinear wave environment, which simulates the ships encountering freak waves at the appointed time and the position in head waves, can be numerically modeled. The numerical wave tank that is generated by the computational fluid dynamics (CFD) can analyze the variation of the freak parameters, peak and steepness, and the velocity components varied with water depth. It can reflect the evolution characteristics of the freak waves to some extent, verify and validate the numerical simulation.

        Keywords:freak wave; numerical wave tank; numerical simulation; phase modulation

        [中圖分類號(hào)]U661.71

        [文獻(xiàn)標(biāo)志碼]A

        [文章編號(hào)]1001-9855(2016)02-0025-06

        [基金項(xiàng)目]十二五預(yù)研項(xiàng)目(51314030101),大連市科技基金(2012J21DW027),海軍大連艦艇學(xué)院科研發(fā)展基金。

        [收稿日期]2015-06-28;[修回日期]2015-12-29

        [作者簡(jiǎn)介]張永勝(1973-),男,博士,高級(jí)工程師,研究方向:海浪非線性和艦船耐波性。

        猜你喜歡
        數(shù)值模擬
        基于AMI的雙色注射成型模擬分析
        錐齒輪精密冷擺輾成形在“材料成型數(shù)值模擬”課程教學(xué)中的應(yīng)用
        基于氣象信息及風(fēng)場(chǎng)信息的風(fēng)機(jī)輪轂處風(fēng)速預(yù)測(cè)
        鉆孔灌注樁樁底沉渣對(duì)樁體承載特性影響的模擬分析
        西南地區(qū)氣象資料測(cè)試、預(yù)處理和加工研究報(bào)告
        科技資訊(2016年18期)2016-11-15 08:01:18
        張家灣煤礦巷道無(wú)支護(hù)條件下位移的數(shù)值模擬
        科技視界(2016年18期)2016-11-03 23:14:27
        張家灣煤礦開(kāi)切眼錨桿支護(hù)參數(shù)確定的數(shù)值模擬
        科技視界(2016年18期)2016-11-03 22:57:21
        跨音速飛行中機(jī)翼水汽凝結(jié)的數(shù)值模擬研究
        科技視界(2016年18期)2016-11-03 20:38:17
        姚橋煤礦采空區(qū)CO2防滅火的數(shù)值模擬分析
        雙螺桿膨脹機(jī)的流場(chǎng)數(shù)值模擬研究
        科技視界(2016年22期)2016-10-18 14:53:19
        国产成人综合亚洲国产| 国产乱人伦av在线无码| 久久青草伊人精品| 中文字幕国内一区二区| 亚洲一区二区三区偷拍厕所| 色视频线观看在线网站| 国产精品麻花传媒二三区别| 最新手机国产在线小视频| 久久精品国产亚洲av日韩一| 欧美激情肉欲高潮视频| 天美传媒精品1区2区3区| 妺妺窝人体色www聚色窝| 视频一区精品中文字幕| 国产精品久免费的黄网站| av蓝导航精品导航| 欧美中出在线| 国产午夜免费一区二区三区视频| 日韩精品久久无码中文字幕 | 可以免费观看的毛片| 精品国产一区二区三区久久狼| 青青草高中生在线视频| 亚洲女初尝黑人巨高清| 久久国产亚洲精品超碰热| 小草手机视频在线观看| 极品少妇被黑人白浆直流| 国产成人综合亚洲精品| 99热在线播放精品6| 美腿丝袜在线观看视频| 亚洲国产av玩弄放荡人妇| 无码国产精品一区二区免| 国产av一区仑乱久久精品| 亚洲av综合av一区二区三区| 国产色a在线观看| 国产AV无码无遮挡毛片| 色婷婷久久精品一区二区| 精品乱码久久久久久久| 亚洲国产精品悠悠久久琪琪| 国产女人av一级一区二区三区 | 国产乱妇无码大片在线观看| 国产成人精品三级麻豆| 亚洲国产精品美女久久久|