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

凹形運(yùn)動體三自由度入水運(yùn)動特性數(shù)值研究

2018-10-15 06:29:08路麗睿魏英杰王海洋宋武超
艦船科學(xué)技術(shù) 2018年9期
關(guān)鍵詞:方向

路麗睿,魏英杰,王 聰,王海洋,宋武超

(1. 哈爾濱工業(yè)大學(xué) 航天學(xué)院,黑龍江 哈爾濱 150001;2. 中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所,江蘇 蘇州 215000)

0 引 言

入水是指運(yùn)動體由空氣介質(zhì)穿越自由液面進(jìn)入到水介質(zhì)的過程。該過程涉及到空泡流動、湍動渦等復(fù)雜現(xiàn)象,具有強(qiáng)非線性和非定常特性。同時(shí)多相流場與運(yùn)動體運(yùn)動存在耦合影響,運(yùn)動體運(yùn)動特性較之單一水下運(yùn)動較為復(fù)雜。

19世紀(jì)初,Gekle等[1-2]對回轉(zhuǎn)體垂直入水開展了數(shù)值模擬研究,獲得了空泡發(fā)展演化規(guī)律和深閉合位置與弗汝德數(shù)間的關(guān)系,建立了空泡閉合射流模型。Quan等[3]基于Navier-Stokes方程對入水空泡生成、發(fā)展、頸縮和閉合過程進(jìn)行模擬。Reinhard等[4]開展了二維對稱模型入水過程空泡形成和擴(kuò)張的數(shù)值研究,建立了基于Wagner理論和水動力占優(yōu)假設(shè)的入水空泡預(yù)測數(shù)值模型,分別研究了拋物線型和對稱楔形對入水空泡形成和擴(kuò)張的影響規(guī)律。Iranmanesh等[5]基于VOF自由液面追蹤技術(shù)開展了低弗汝德數(shù)下水平圓柱入水的數(shù)值仿真研究,研究了圓柱直徑、入水速度對入水空泡流的影響。葉取源等[6]分別采用非線性自由面理論和歐拉-拉格朗日混合邊界元方法模擬了錐體和圓盤入水空泡發(fā)生、發(fā)展、面閉合和深閉合等流動現(xiàn)象,定性地分析了運(yùn)動體質(zhì)量和速度對入水空泡演化過程的影響。Gong等[7]在SPH方法基礎(chǔ)上,引入非反射邊界條件(non-reflection boundary),開展了二維平面楔形體入水過程的數(shù)值模擬研究。張珂等[8 – 9]采用單相格子Boltzmann流動模型開展了圓柱體垂直入水二維和三維數(shù)值模擬,結(jié)果表明單相格子Boltzmann流動模型可以較為精確地捕捉到入水過程中的部分流動特征。王聰?shù)萚10]對球體高速入水開展數(shù)值仿真,研究了空氣壓力對入水空演化的影響。馬慶鵬等[11]分別對球體和錐頭圓柱體入水空泡及流場環(huán)境影響開展了模擬研究。王平等[12]開展了楔形體波浪中入水?dāng)?shù)值研究,研究了波浪對楔形體受力及入水姿態(tài)的影響。

凹形運(yùn)動體作為一種復(fù)雜結(jié)構(gòu)運(yùn)動體,具有廣泛的工程應(yīng)用背景。目前僅有路中磊等[13 – 14]以火箭助推器入水為背景開展過“凹體”入水相關(guān)研究,其研究主要關(guān)注入水過程流場特性,對運(yùn)動特性關(guān)注較少。除火箭助推器外,導(dǎo)彈尾罩亦屬于凹形運(yùn)動體。在潛射導(dǎo)彈冷發(fā)射過程中,導(dǎo)彈點(diǎn)火時(shí)將尾罩從一側(cè)彈開,其運(yùn)動具有典型的多自由度特性,研究尾罩入水運(yùn)動特性對避免尾罩砸艇具有重要意義。本文以此為背景,開展基于重疊網(wǎng)格技術(shù)的凹形運(yùn)動體三自由度入水過程運(yùn)動特性的數(shù)值研究,分析初始垂向速度、水平遷移速度對凹形運(yùn)動體三自由度入水過程運(yùn)動特性及水動力特性的影響。

1 數(shù)值計(jì)算理論及方法

1.1 基本控制方程

入水運(yùn)動是一個(gè)非定常問題,流動參數(shù)均隨時(shí)空變化,主要表現(xiàn)在質(zhì)量擴(kuò)散、動量黏滯耗散等方面。在流體流動過程中,混合相滿足質(zhì)量守恒和動量守恒方程。

連續(xù)性方程為:

動量守恒方程為:

1.2 多相流模型

低速入水空泡氣液兩相介質(zhì)互不摻混,空泡壁面處兩相流速基本一致,因此可采用均質(zhì)平衡多相流理論[15]。本文選用的VOF多相流模型是一種適用于固定歐拉網(wǎng)格的界面捕捉模型,該模型可以較好描述流動過程中多相流體的界面位置。VOF模型將兩相流體視為一種密度可變的混合流體,通過逐一求解計(jì)算流域內(nèi)每一控制單元的控制方程得到混合流體參數(shù),并根據(jù)各相流體所占的體積分?jǐn)?shù)去實(shí)現(xiàn)對流體間界面的追蹤。混合相密度和動力粘度分別如下:

1.3 湍流模型

湍流是一種具有強(qiáng)非線性的流動狀態(tài),同時(shí)也是最為普遍的流動狀態(tài)。本文采用Realizable,該模型將渦粘系數(shù)()和湍動能()與耗散率()聯(lián)系起來,分別建立湍動能()與耗散率()的輸運(yùn)方程,如式(5)和式(6)所示,最終使整個(gè)控制方程封閉。

1.4 重疊網(wǎng)格技術(shù)

重疊網(wǎng)格的基本思想是對復(fù)雜流域的分解與網(wǎng)格的重新組合[16]。重疊網(wǎng)格技術(shù)將流域分解為若干個(gè)子區(qū)域,各子區(qū)域計(jì)算網(wǎng)格在計(jì)算中相互獨(dú)立,網(wǎng)格之間存在重疊、嵌套和覆蓋關(guān)系。在重疊區(qū)域中,各子區(qū)域計(jì)算網(wǎng)格通過流場信息的插值與映射進(jìn)行數(shù)據(jù)交換,從而建立起了各子區(qū)域網(wǎng)格之間的耦合關(guān)系,為子區(qū)域流場計(jì)算提供邊界條件。重疊網(wǎng)格的主要過程包含尋點(diǎn)、挖洞、建立網(wǎng)格重疊關(guān)系3個(gè)部分。在尋點(diǎn)過程中,網(wǎng)格A落入網(wǎng)格B中的節(jié)點(diǎn)稱為洞內(nèi)點(diǎn),洞內(nèi)點(diǎn)在數(shù)值求解過程中不參加計(jì)算,與洞內(nèi)點(diǎn)緊鄰的點(diǎn)稱為邊界插值點(diǎn),邊界插值點(diǎn)通過插值與映射接收來自其他子區(qū)域網(wǎng)格的數(shù)據(jù)信息,并將所接受的數(shù)據(jù)作為該區(qū)域求解的邊界條件。如圖1所示,網(wǎng)格A通過其邊界插值點(diǎn)接收B網(wǎng)格內(nèi)求解得到的流場信息,同樣,網(wǎng)格B也從其外邊界插值點(diǎn)獲得A網(wǎng)格的流場數(shù)據(jù),并將該數(shù)據(jù)作為求解邊界條件。

1.5 數(shù)值計(jì)算方法

本文基于Fluent 17.2平臺,采用有限體積法離散流動控制方程;壓力場與速度場的耦合選用Coupled算法;壓力場的空間離散采用PRESTO!格式;各相體積率離散采用CICSAM格式;各求解變量的離散采用一階迎風(fēng)格式;對流項(xiàng)采用QUICK離散格式;相界面的幾何重構(gòu)采用Geo-Reconstruct格式。通過C語言程序,嵌入Fluent用戶自定義函數(shù)定義入水運(yùn)動體慣性矩及邊界壓力,實(shí)現(xiàn)模型三自由度運(yùn)動。

早年嫁過來時(shí),采蓮掮著鋤頭跟福生做農(nóng)活,不惟下湖,還得去崖上。崖頭,就是高處的地坡。正午坐在地頭上喝水,偶爾一搭眼,就看到天宇下的崖坡上,一人一牛在耕地,遠(yuǎn)看兩個(gè)黑點(diǎn),背后一抹駱駝云,一趟一趟,不厭其煩地劃著圈……此后,她的夢里始終留存著這樣一組畫面。農(nóng)夫,耕牛,在天宇下的土坡上,永遠(yuǎn)沒有盡頭地轉(zhuǎn)悠,一簇牛角上的紅布條像火苗似的燃燒著。而那一聲吆牛號子,那份游蕩與戚然交織的天籟啊,簡直就是入心入肺了。[1]12

2 計(jì)算模型及網(wǎng)格劃分

2.1 幾何模型

本文計(jì)算所用凹形運(yùn)動體幾何模型如圖2所示,模型幾何參數(shù)及屬性見表1。

表 1 模型幾何參數(shù)及屬性Tab. 1 The parameter of geometric model

2.2 計(jì)算域及邊界條件

為簡化計(jì)算,可假設(shè)凹形運(yùn)動體在z方向不存在速度分量,且不存在繞y軸的轉(zhuǎn)動速度分量。故本算例的數(shù)值計(jì)算選擇在三維、三自由度情況下進(jìn)行。其計(jì)算域及邊界條件如圖4所示,xy平面為對稱平面,計(jì)算域?qū)挒?0D,高為7D,厚為10D。以觸水時(shí)刻作為時(shí)間零點(diǎn),觸水點(diǎn)為坐標(biāo)原點(diǎn),重力方向?yàn)閥負(fù)方向。

2.3 網(wǎng)格劃分

圖4給出了計(jì)算域網(wǎng)格和子域網(wǎng)格的劃分方式,2套網(wǎng)格均采用正六面體結(jié)構(gòu)化網(wǎng)格,以較好控制網(wǎng)格質(zhì)量和網(wǎng)格數(shù)量。為保證空泡區(qū)域計(jì)算精度,并保持較好的過度效果,對凹形運(yùn)動體附近區(qū)域采用局部加密處理,運(yùn)動體向遠(yuǎn)場過渡區(qū)域采用漸疏處理。

2.4 數(shù)值方法驗(yàn)證

為了驗(yàn)證數(shù)值方法的有效性,基于上述計(jì)算方法及邊界條件,對Aristoff[17]所開展的球體入水實(shí)驗(yàn)進(jìn)行計(jì)算。實(shí)驗(yàn)研究所選球體直徑D=25.4 mm,球體密度為 2.3 g/cm3,入水速度 v0=2.17 m/s。本節(jié)數(shù)值計(jì)算采用1:1的計(jì)算域模型,運(yùn)動參數(shù)初始條件設(shè)置與Aristoff實(shí)驗(yàn)狀態(tài)完全一致,計(jì)算邊界條件及重疊網(wǎng)格劃分方法與2.2節(jié)和2.3節(jié)所述相同。圖5和圖6分別給出了計(jì)算所得空泡形態(tài)和運(yùn)動軌跡與試驗(yàn)的對比圖,從圖中可以看出,數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)具有較好的一致性。

3 結(jié)果分析

3.1 凹形運(yùn)動體三自由度入水運(yùn)動過程

導(dǎo)彈尾罩在拋落時(shí),為了避免尾罩垂向下落砸艇,往往會在導(dǎo)彈尾罩一側(cè)安裝拉斷式鉸鏈,使尾罩彈開后向一側(cè)飛離,如圖7所示。因此本文對凹形運(yùn)動體三自由運(yùn)動初始方向進(jìn)行如圖7所示定義,水平遷移速度沿x正方向,垂向速度沿y負(fù)方向,旋轉(zhuǎn)方向正負(fù)由右手定則給出,初始角速度方向?yàn)檎?/p>

圖8給出了凹形運(yùn)動體以135°傾角、垂向初始速度 v0=7.5 m/s、水平遷移速度 u0=2 m/s、角速度 ω0=2.0 rad/s狀態(tài)入水過程氣液兩相體積分?jǐn)?shù)變化云圖。在凹形運(yùn)動體入水過程中,下邊緣撞擊自由液面,流體質(zhì)點(diǎn)獲得來自運(yùn)動體的動能而產(chǎn)生向外排開的運(yùn)動,如圖8中t=0 ms所示,由于運(yùn)動體凸面與凹面幾何特性的差異,導(dǎo)致凹面尖端更易誘導(dǎo)空泡形成。隨著入水時(shí)間的增加,噴濺液體在凹面內(nèi)部閉合,形成隨動空泡,如圖8中t=5 ms所示。在閉合瞬間,噴濺液體在凹面上半側(cè)產(chǎn)生抨擊作用而形成力矩,運(yùn)動體旋轉(zhuǎn)方向發(fā)生改變,定義該運(yùn)動為回旋運(yùn)動。運(yùn)動體的回旋導(dǎo)致由上緣誘導(dǎo)形成的二次空泡被橫向拉伸,空泡橫向發(fā)展,如圖8中t=10 ms至t=15 ms所示。隨著空泡的橫向發(fā)展,迎流面一側(cè)水流體在重力作用下下落,與背流面一側(cè)空泡界面相遇,二次空泡完成閉合,如圖8中t=20 ms至t=25 ms所示。二次空泡受慣性力作用繼續(xù)做橫向運(yùn)動,逐漸與自由液面邊界分離而形成了獨(dú)立空泡,如圖8中t=30 ms至35 ms所示。隨動空泡亦由于凹形運(yùn)動體旋轉(zhuǎn)運(yùn)動逐漸脫落形成獨(dú)立空泡,如圖8中t=15 ms所示。獨(dú)立空泡形態(tài)脫落后受運(yùn)動體凹面旋轉(zhuǎn)運(yùn)動形成的剪切力作用而逐漸卷曲拉長,如圖 8 中 t=15 ms至 t=35 ms所示。

對凹形運(yùn)動體所受到的水動力系數(shù)做如下定義:

圖9給出了凹形運(yùn)動體入水過程運(yùn)動特性曲線,圖10給出了該過程水動力特性曲線,圖11給出了凹形結(jié)構(gòu)體位置姿態(tài)隨時(shí)間變化簡圖。從圖9中可以看出在入水之前,凹形運(yùn)動體y方向速度以較小變化率減小,角速度亦呈現(xiàn)較為平緩的正向增加趨勢;在入水后的極短時(shí)間內(nèi),y方向速度迅速減小,角速度也達(dá)到峰值,之后迅速減小為0并反向增大,稱該階段為沖擊階段;之后,y方向速度變化率逐漸減緩,速度持續(xù)減小并逐漸趨于0,角速度變化率在t=5~13 ms之間出現(xiàn)短暫的微幅增大,t=10 ms之后又逐漸減小,角速度達(dá)到負(fù)方向峰值后亦逐漸減小。對比圖10可知,在入水之前,凹形運(yùn)動體靠近自由液面,其下緣與自由液面間形成氣墊,受空氣墊阻力作用,y方向速度緩慢減小,由于空氣墊阻力與凹形運(yùn)動體所示重力形成了與角速度方向相同的力矩,導(dǎo)致角速度亦呈現(xiàn)緩慢增加趨勢。在觸水后,進(jìn)入沖擊階段,沖擊階段主要包含入水沖擊和噴濺沖擊2個(gè)部分,在入水沖擊過程中,凹形運(yùn)動體受到入水沖擊力作用,y方向速度迅速減小。y方向入水沖擊力在質(zhì)心處產(chǎn)生與角速度方向相同的力矩,x方向沖擊力在質(zhì)心處產(chǎn)生與角速度方向相反的力矩。由于y方向沖擊力遠(yuǎn)大于x方向沖擊力,且其值遠(yuǎn)大于空氣墊阻力,因此合力矩與角速度方向相同,運(yùn)動體角速度迅速增加。隨著入水時(shí)間的增加,運(yùn)動體背流面所形成的液體噴濺達(dá)到凹形運(yùn)動體內(nèi)壁面上側(cè)(圖8中t=0~5 ms之間),此時(shí)入水沖擊階段已經(jīng)完成,運(yùn)動體主要受噴濺沖擊力作用,進(jìn)入噴濺沖擊階段。噴濺沖擊力在質(zhì)心處所產(chǎn)生的力矩方向與旋轉(zhuǎn)方向相反,致使角速度迅速減小并反向增加,形成回旋運(yùn)動。噴濺沖擊階段結(jié)束后,凹形運(yùn)動體下緣形成了空泡隨動,且由于空泡內(nèi)壓力小于凸面處所受環(huán)境壓力,導(dǎo)致迎流面與背流面壓差形成力矩作用,該力矩方向與運(yùn)動體回旋方向相同,導(dǎo)致角速度變化率短暫增大。在t=13 ms附近,隨動空泡脫落,運(yùn)動體處于全沾濕狀態(tài),其旋轉(zhuǎn)受到流體阻力和轉(zhuǎn)動慣性力共同作用,最終在t=18 ms時(shí)角速度達(dá)到峰值并開始減小。從圖11還可以看出,由于運(yùn)動體撞水面積較小,使得入水沖擊較噴濺沖擊弱;在t=14 ms附近,運(yùn)動體在x方向的迎流面積達(dá)到最大,x方向阻力達(dá)到峰值。結(jié)合圖8中t=10~16 ms階段,由于運(yùn)動體上緣與空泡尾腔分離,空泡界面與運(yùn)動體邊緣之間存在張力,張力的拖曳作用在y方向的分量使得y方向阻力開始微幅增大,并在t=16 ms時(shí)達(dá)到峰值。

3.2 垂向速度對凹形運(yùn)動體運(yùn)動特性的影響

考慮不同垂向速度,圖12給出了凹形運(yùn)動體入水過程質(zhì)心平移特性變化規(guī)律,分別取初始垂向速度v0=5.0 m/s,7.5 m/s,10.0 m/s ,初始水平遷移速度 u0=1 m/s,初始角速度 ω0=2.0 rad/s。從圖 12(a)中可以看出隨著垂向速度的增加,y方向速度趨于0時(shí),其最大位移逐漸增大,但增大幅度隨速度的增加而減小。由圖12(b)可知,在沖擊階段,隨著垂向速度的增加,速度衰減程度逐漸增大。從圖12(c)中還可發(fā)現(xiàn),隨著垂向速度的增大,“Z”型路徑逐漸平坦,表明垂向速度增加對運(yùn)動體回?cái)[存在抑制作用。

為分析上述現(xiàn)象原因,圖13給出了運(yùn)動體x方向和y方向的水動力系數(shù)變化曲線。從圖13(b)可以看出,隨著入水速度的增加,y方向入水沖擊力以及噴濺沖擊力峰值增量大幅加大,導(dǎo)致沖擊引起的速度衰減幅度增加。在入水沖擊之后,y方向受到的阻力亦隨著速度增加而加大,沖擊與阻力對y方向最大位移具有負(fù)激勵作用,隨著速度的增加,負(fù)激勵增大,導(dǎo)致位移增加幅度減小。

圖14給出了垂向速度對旋轉(zhuǎn)運(yùn)動特性的影響規(guī)律,從圖中可以看出,隨著垂向速度的增大,入水初期角位移峰值逐漸前移,角速度峰值逐漸增大。回旋運(yùn)動階段角速度亦隨著垂向速度的增大而增大。由于回旋運(yùn)動主要由噴濺沖擊引起,噴濺沖擊力對回旋運(yùn)動角速度具有正激勵作用,從圖13(a)可以看出,隨著垂向速度的增加,噴濺沖擊力大幅增加,導(dǎo)致回旋運(yùn)動角速度增大。

3.3 水平遷移速度對凹形運(yùn)動體運(yùn)動特性的影響

圖15給出了水平遷移速度對凹形運(yùn)動體平移運(yùn)動的影響規(guī)律,初始水平遷移速度分別取u0=2.0 m/s,3.0 m/s,4.0 m/s,初始垂向速度取 v0=6.0 m/s,初始角速度ω0=2.0 rad/s。從圖15(a)中可以看出,不同水平遷移速度下,凹形運(yùn)動體運(yùn)動具有相同的運(yùn)動趨勢,且隨著水平遷移速度增加,y方向速度趨于0時(shí),其位移最大值增加。當(dāng)u0=3 m/s時(shí),回?cái)[運(yùn)動引起了質(zhì)心的上漂。從圖15(b)中可發(fā)現(xiàn),隨著水平遷移速度的增加,入水沖擊階段,y方向速度衰減幅值逐漸減小。觀察圖15(c)可發(fā)現(xiàn),隨著速度的增加,質(zhì)心軌跡回?cái)[幅值增大,質(zhì)心運(yùn)動軌跡的“Z”型特征隨著水平遷移速度的增加而逐漸顯著。

圖16給出了水平遷移速度對x和y方向水動力系數(shù)的影響規(guī)律,從圖中可發(fā)現(xiàn),沖擊階段y方向受到的入水及噴濺沖擊作用與入水瞬間速度衰減幅值呈現(xiàn)相同趨勢。隨著水平速度的增加,x方向的入水沖擊力增大,但噴濺沖擊力卻減小。出現(xiàn)該現(xiàn)象是由于水平速度增加,噴濺到達(dá)運(yùn)動體背流面時(shí)與運(yùn)動體的相對速度減小,導(dǎo)致噴濺沖擊作用減小。在沖擊階段之后,y方向所受到的阻力以及x方向阻力均隨水平遷移速度的增加而增大。

4 結(jié) 語

本文基于重疊網(wǎng)格技術(shù),開展了凹形運(yùn)動體入水過程運(yùn)動特性數(shù)值研究,分析了初始垂向速度、水平遷移速度對凹形運(yùn)動體入水過程運(yùn)動特性及水動力特性的影響,得到如下結(jié)論:

1)在凹形運(yùn)動體三自由度入水過程中,噴濺沖擊作用對凹形運(yùn)動體角運(yùn)動具有重要影響,運(yùn)動體在水中旋轉(zhuǎn)方向與噴濺沖擊形成的力矩方向相同。

2)凹形運(yùn)動體入水運(yùn)動存在2個(gè)運(yùn)動占優(yōu)階段,在入水前期為平移運(yùn)動占優(yōu),入水后期為旋轉(zhuǎn)運(yùn)動占優(yōu)。

3)凹形運(yùn)動體在下落過程中,質(zhì)心運(yùn)動軌跡會產(chǎn)生回?cái)[而呈現(xiàn)“Z”型特征,垂向速度的增加對質(zhì)心回?cái)[運(yùn)動具有抑制作用,水平遷移速度的增加對其具有促進(jìn)作用。

猜你喜歡
方向
2023年組稿方向
方向
青年運(yùn)動的方向(節(jié)選)
2022年組稿方向
2022年組稿方向
2021年組稿方向
如何確定位置與方向
2021年組稿方向
2021年組稿方向
大自然中的方向
主站蜘蛛池模板: 国产自在线播放| 亚洲二区视频| 秋霞一区二区三区| 青青青草国产| 丁香五月婷婷激情基地| 国产高清在线精品一区二区三区| 久久综合激情网| 97视频在线观看免费视频| 手机在线看片不卡中文字幕| 欧洲亚洲一区| 亚洲国产成人自拍| 国产尹人香蕉综合在线电影| 色视频国产| 国产精品极品美女自在线网站| 99re经典视频在线| 91精品视频播放| 欧美区在线播放| 午夜日本永久乱码免费播放片| 91口爆吞精国产对白第三集| 亚洲国产成熟视频在线多多 | 新SSS无码手机在线观看| 91午夜福利在线观看| 欧美精品在线看| 亚洲中文字幕在线观看| 精品国产aⅴ一区二区三区| 欧美成a人片在线观看| 亚洲一区国色天香| 亚洲美女一级毛片| 九九九精品成人免费视频7| 小说 亚洲 无码 精品| 99视频精品全国免费品| 99re视频在线| 成年网址网站在线观看| 亚洲熟妇AV日韩熟妇在线| 免费看a毛片| 亚洲国产精品VA在线看黑人| 亚洲首页在线观看| 亚洲国产精品无码AV| 99久久精品国产自免费| 日本三区视频| 日韩毛片免费| 午夜福利在线观看成人| 一级毛片免费高清视频| 欧美激情福利| 午夜国产精品视频| 国产免费a级片| 日韩高清欧美| 中文精品久久久久国产网址 | 久久夜色精品| 久996视频精品免费观看| 国产成人综合日韩精品无码首页 | 毛片在线看网站| 久夜色精品国产噜噜| 香蕉eeww99国产在线观看| 色婷婷在线影院| 国产一二视频| 午夜一区二区三区| 男女精品视频| 亚洲国产清纯| 国产网站免费观看| 日韩小视频在线播放| 九九热在线视频| 免费看黄片一区二区三区| 一级成人a做片免费| 日韩精品毛片| 日本人妻一区二区三区不卡影院| 午夜在线不卡| 国产精品毛片在线直播完整版| 国内熟女少妇一线天| 色国产视频| 成人午夜久久| 日韩精品少妇无码受不了| 日韩一区精品视频一区二区| 手机在线免费不卡一区二| 免费观看国产小粉嫩喷水| 国产在线自在拍91精品黑人| 人妻精品久久无码区| 久久semm亚洲国产| 亚洲熟女偷拍| 国产一区二区福利| 免费又爽又刺激高潮网址| 久久青草热|