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

基于重疊網(wǎng)格的吊艙推進(jìn)器偏轉(zhuǎn)工況水動(dòng)力性能數(shù)值模擬*

2016-05-25 02:30:01徐嘉啟時(shí)立攀

徐嘉啟 熊 鷹 時(shí)立攀

(海軍工程大學(xué)艦船工程系1) 武漢 430033) (海軍航空工程學(xué)院青島校區(qū)2) 青島 266000)

?

基于重疊網(wǎng)格的吊艙推進(jìn)器偏轉(zhuǎn)工況水動(dòng)力性能數(shù)值模擬*

徐嘉啟1)熊鷹1)時(shí)立攀2)

(海軍工程大學(xué)艦船工程系1)武漢430033)(海軍航空工程學(xué)院青島校區(qū)2)青島266000)

摘要:為了更靈活高效地對(duì)吊艙推進(jìn)器偏轉(zhuǎn)工況下水動(dòng)力性能進(jìn)行數(shù)值模擬,采用重疊網(wǎng)格劃分螺旋槳計(jì)算流域和某吊艙推進(jìn)器偏轉(zhuǎn)工況下的計(jì)算流域,對(duì)4381型、4119型螺旋槳的敞水性能進(jìn)行了數(shù)值預(yù)報(bào),將重疊網(wǎng)格和滑移網(wǎng)格獲得的計(jì)算結(jié)果進(jìn)行了比較分析,驗(yàn)證了重疊網(wǎng)格方法應(yīng)用于螺旋槳水動(dòng)力性能計(jì)算的可行性.對(duì)吊艙推進(jìn)器偏轉(zhuǎn)工況下水動(dòng)力性能進(jìn)行預(yù)報(bào)和分析.結(jié)果表明,應(yīng)用重疊網(wǎng)格獲得的計(jì)算結(jié)果具有較高的計(jì)算精度,且重疊網(wǎng)格處理吊艙推進(jìn)器偏轉(zhuǎn)工況等多體相對(duì)運(yùn)動(dòng)問題更高效.

關(guān)鍵詞:重疊網(wǎng)格;螺旋槳;吊艙推進(jìn)器;偏轉(zhuǎn)工況;水動(dòng)力性能;數(shù)值模擬

0引言

近年來,重疊網(wǎng)格技術(shù)[1]發(fā)展迅速.尤其自2006年以來,其模擬多體間大幅相對(duì)運(yùn)動(dòng)的繞流流場(chǎng)的能力在空氣動(dòng)力學(xué)數(shù)值模擬中應(yīng)用愈發(fā)廣泛,包括固定翼飛機(jī)后緣襟翼的偏轉(zhuǎn)運(yùn)動(dòng)[2]、直升機(jī)旋翼與機(jī)身相互作用及旋翼變距運(yùn)動(dòng)[3]、飛機(jī)外掛物投送[4]、風(fēng)力發(fā)電機(jī)螺旋槳與風(fēng)機(jī)頭部相互運(yùn)動(dòng)[5]等,均取得與試驗(yàn)數(shù)據(jù)吻合較好.重疊網(wǎng)格技術(shù)在船舶計(jì)算流體力學(xué)仿真中的應(yīng)用也日趨廣泛,包括救生艙投放過程數(shù)值模擬、海洋結(jié)構(gòu)物與海水相互作用數(shù)值模擬[6],以及帶有螺旋槳的船舶自航工況數(shù)值模擬[7].

在上述應(yīng)用實(shí)例中的帶襟翼的固定機(jī)翼、旋翼、飛機(jī)外掛物、救生艙艙體以及船體周圍均存在運(yùn)動(dòng)邊界,解決運(yùn)動(dòng)邊界的繞流問題采用的動(dòng)網(wǎng)格技術(shù)主要有網(wǎng)格再生方法、網(wǎng)格變形方法和重疊網(wǎng)格方法.重疊網(wǎng)格方法是一種網(wǎng)格組合和嵌套策略.該方法的優(yōu)點(diǎn)在于不需要網(wǎng)格再生和變形,可處理多物體或者物體內(nèi)子部件之間大幅相對(duì)運(yùn)動(dòng)的問題.

文中結(jié)合重疊網(wǎng)格技術(shù)的上述優(yōu)點(diǎn),采用RANS方法與Realizable k-ε湍流模型首先對(duì)螺旋槳敞水性能進(jìn)行預(yù)報(bào),通過與試驗(yàn)結(jié)果的比較驗(yàn)證該方法的可行性;再用相同方法預(yù)報(bào)某吊艙推進(jìn)器偏轉(zhuǎn)工況下的水動(dòng)力性能.

1控制方程與湍流模型

1.1控制方程

湍流的控制方程引入張量符號(hào)表示,連續(xù)方程

(1)

動(dòng)量方程(雷諾平均納維斯托克斯方程(RANS方程))

(2)

(3)

依據(jù)Boussinesq提出的渦粘假定,引入湍動(dòng)粘度μt,把雷諾應(yīng)力表示成湍動(dòng)粘度的函數(shù),其張量形式為

(4)

式中:μt為湍動(dòng)粘度,Pa·s;ui為時(shí)均速度,上標(biāo)“′”代表脈動(dòng)值,上標(biāo)“ˉ”代表平均值;δij為Kronecker delta函數(shù);k為湍動(dòng)能(turbulent kinematic energy),m2/s2;ε為湍動(dòng)耗散率(turbulent dissipation rate),m2/s3.

1.2湍流模型

數(shù)值計(jì)算使用Realizable k-ε湍流模型,該模型對(duì)正應(yīng)力進(jìn)行數(shù)學(xué)約束,避免了Standard k-ε模型處理時(shí)均應(yīng)變率特別大的情形時(shí)可能出現(xiàn)正應(yīng)力為負(fù)的情況.k和ε的輸運(yùn)方程如下.

(5)

(6)

式中:σk=1.0;σε=1.2;C2=1.9;

η=(2Eij·Eij)1/2k/ε;

(7)

式(5)中,湍動(dòng)粘度μt按下式計(jì)算.

(8)

(9)

(10)

Ωij是從角速度為ωk的參考系中觀察到的時(shí)均轉(zhuǎn)動(dòng)速度張量,對(duì)無旋轉(zhuǎn)的流場(chǎng),上式中U*計(jì)算式根號(hào)中的第二項(xiàng)為零,該項(xiàng)用來表示旋轉(zhuǎn)的影響.Gk是由平均速度梯度引起的湍動(dòng)能的產(chǎn)生項(xiàng),Gk=2μtE2,E=(2Eij·Eij)1/2.

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

采用有限體積法離散控制方程和湍流模型,對(duì)流項(xiàng)和擴(kuò)散項(xiàng)均采用二階迎風(fēng)格式進(jìn)行離散,壓力速度耦合迭代采用SIMPLE方法,采用加強(qiáng)型壁面處理方式.

螺旋槳敞水性能預(yù)報(bào)采用定常計(jì)算方法;吊艙推進(jìn)器偏轉(zhuǎn)工況水動(dòng)力性能預(yù)報(bào)先采用定常方法計(jì)算,待收斂后再采用非定常方法,以加快收斂進(jìn)程,時(shí)間步長為螺旋槳旋轉(zhuǎn)3.6°所對(duì)應(yīng)的時(shí)間.

3推進(jìn)器流場(chǎng)網(wǎng)格

3.1螺旋槳流場(chǎng)網(wǎng)格

螺旋槳模型主要參數(shù)包括:直徑D=0.305 m;轂徑比d/D=0.2;4381型槳葉數(shù)z=5,盤面比AE/AP=0.725,側(cè)斜角θS=0°,無縱傾;4119型槳葉數(shù)z=3,盤面比AE/AP=0.6,側(cè)斜角θS=0°,無縱傾.計(jì)算網(wǎng)格總數(shù)均在180萬左右.

3.1.1內(nèi)外域交界面為重疊邊界的螺旋槳流場(chǎng)重疊網(wǎng)格

螺旋槳流場(chǎng)的網(wǎng)格離散從劃分流場(chǎng)開始,整個(gè)流場(chǎng)被分為2個(gè)區(qū)域,一個(gè)是螺旋槳附近的旋轉(zhuǎn)區(qū)域,另一個(gè)是遠(yuǎn)場(chǎng)的固定區(qū)域,2個(gè)區(qū)域均采用非結(jié)構(gòu)網(wǎng)格.首先,單獨(dú)生成2個(gè)區(qū)域的網(wǎng)格,2個(gè)區(qū)域外形均為圓柱狀;然后,將旋轉(zhuǎn)區(qū)域圓柱外表面設(shè)置為重疊邊界,并建立旋轉(zhuǎn)區(qū)域與固定區(qū)域網(wǎng)格的重疊網(wǎng)格;之后,選擇重疊網(wǎng)格的插值計(jì)算方法,在此應(yīng)用線性插值方法;最后,進(jìn)行初始化過程,該過程主要包括挖洞(從固定區(qū)域中挖去與旋轉(zhuǎn)區(qū)域重疊的區(qū)域)與建立重疊網(wǎng)格關(guān)系.

將一個(gè)相對(duì)于絕對(duì)靜止坐標(biāo)系旋轉(zhuǎn)的坐標(biāo)系固結(jié)在旋轉(zhuǎn)區(qū)域,螺旋槳遠(yuǎn)場(chǎng)的固定區(qū)域仍采用絕對(duì)靜止坐標(biāo)系.這樣,旋轉(zhuǎn)區(qū)域的流場(chǎng)在旋轉(zhuǎn)坐標(biāo)系下是定常的,而遠(yuǎn)場(chǎng)的流域在絕對(duì)靜止坐標(biāo)系下可近似認(rèn)為是定常的,因而簡(jiǎn)化了螺旋槳流場(chǎng)的求解.

圖1 螺旋槳表面網(wǎng)格

圖2 旋轉(zhuǎn)內(nèi)域網(wǎng)格

圖3 螺旋槳計(jì)算流域網(wǎng)格

圖4 重疊網(wǎng)格橫剖面(注意內(nèi)外域交界面處)

3.1.2內(nèi)外域交界面為滑移網(wǎng)格的螺旋槳流場(chǎng)網(wǎng)格

與重疊網(wǎng)格類似,內(nèi)外域分開,但是滑移網(wǎng)格的方法需要在外域提前人工挖去與內(nèi)域重合的部分,在內(nèi)外域交界面建立滑移網(wǎng)格交界面.

圖5 滑移網(wǎng)格橫剖面(注意內(nèi)外域交界面處)

3.2吊艙推進(jìn)器流場(chǎng)重疊網(wǎng)格

吊艙推進(jìn)器模型為加拿大海洋技術(shù)研究所設(shè)計(jì)的某拖式吊艙推進(jìn)器的模型,Mohammed F Islam、Liu Pengfei等對(duì)其直航及偏轉(zhuǎn)工況下的水動(dòng)力性能進(jìn)行了系列試驗(yàn)研究,試驗(yàn)數(shù)據(jù)參考文獻(xiàn)[8].吊艙槳和吊艙的主要參數(shù)見表1~2.

表1 吊艙槳的主要參數(shù)

表2 吊艙的主要參數(shù)

吊艙推進(jìn)器的幾何形狀見圖6,吊艙推進(jìn)器網(wǎng)格橫剖面見圖7,面網(wǎng)格與流場(chǎng)網(wǎng)格見圖8.

圖6 吊艙推進(jìn)器的外形輪廓

圖7 吊艙推進(jìn)器網(wǎng)格橫剖面

圖8 吊艙推進(jìn)器面網(wǎng)格與流場(chǎng)網(wǎng)格

重疊網(wǎng)格在吊艙偏轉(zhuǎn)后,只需再次建立吊艙子域網(wǎng)格與背景固定子域網(wǎng)格的插值關(guān)系,即可開始計(jì)算,無需重新生成,效率大幅提升.

4螺旋槳敞水性能預(yù)報(bào)

4381型螺旋槳計(jì)算流域分別采用重疊網(wǎng)格和滑移網(wǎng)格進(jìn)行了數(shù)值模擬,4119型螺旋槳計(jì)算流域采用重疊網(wǎng)格進(jìn)行了數(shù)值模擬,螺旋槳敞水性能的預(yù)報(bào)結(jié)果均與試驗(yàn)結(jié)果進(jìn)行了對(duì)比分析.其中,以4381型螺旋槳為例,將重疊網(wǎng)格和滑移網(wǎng)格獲得的計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行相互比較分析.計(jì)算過程中螺旋槳轉(zhuǎn)速保持n=1 200 r/min,改變進(jìn)速Va以得到不同的進(jìn)速系數(shù)J.

圖9 采用重疊網(wǎng)格計(jì)算的4381槳敞水性能對(duì)比

圖10 采用滑移網(wǎng)格計(jì)算的4381槳敞水性能對(duì)比

由表1可知,采用重疊網(wǎng)格計(jì)算得到的4381型螺旋槳的推力系數(shù)kt,在設(shè)計(jì)工況下(J=0.899)與試驗(yàn)值的相對(duì)誤差在3%以內(nèi),且所有計(jì)算點(diǎn)的相對(duì)誤差控制在3%以內(nèi).而10kq在設(shè)計(jì)工況下計(jì)算值與試驗(yàn)值的誤差控制在1%以內(nèi),且所有計(jì)算點(diǎn)的相對(duì)誤差控制在2%以內(nèi).

表3 4381型螺旋槳兩種網(wǎng)格計(jì)算值與試驗(yàn)值誤差的對(duì)比

圖11 采用重疊網(wǎng)格計(jì)算的4119槳敞水性能對(duì)比

從表3得出結(jié)論:重疊網(wǎng)格劃分的螺旋槳計(jì)算流域在0.3~1.0進(jìn)速系數(shù)系列中(除進(jìn)速系數(shù)為0.6時(shí)),推力系數(shù)計(jì)算值較滑移網(wǎng)格更接近試驗(yàn)值,且該優(yōu)勢(shì)在低進(jìn)速系數(shù)(J=0.3~0.5)時(shí)更明顯.

從圖11可以看出,采用重疊網(wǎng)格計(jì)算得到的4119型螺旋槳的推力系數(shù)kt在0.5~1.1進(jìn)速系數(shù)系列中,計(jì)算值與試驗(yàn)值吻合較好,最大相對(duì)誤差在6%以內(nèi);在設(shè)計(jì)工況(J=0.833)下計(jì)算值與試驗(yàn)值的相對(duì)誤差小于1%.10kq在設(shè)計(jì)工況下計(jì)算值與試驗(yàn)值的誤差小于2%,且所有計(jì)算點(diǎn)的相對(duì)誤差控制在6%以內(nèi).因此應(yīng)用重疊網(wǎng)格能夠模擬螺旋槳旋轉(zhuǎn)流場(chǎng)并獲得較高的螺旋槳敞水性能計(jì)算精度.

5吊艙推進(jìn)器偏轉(zhuǎn)工況水動(dòng)力性能預(yù)報(bào)

船舶在回轉(zhuǎn)和Z字操縱時(shí),吊艙推進(jìn)器工作在遠(yuǎn)離設(shè)計(jì)點(diǎn)的工況下,船體和吊艙推進(jìn)器因而承受巨大荷載,嚴(yán)重可致結(jié)構(gòu)破壞,因此偏轉(zhuǎn)工況下吊艙推進(jìn)器的水動(dòng)力性能預(yù)報(bào)是十分重要的.本節(jié)結(jié)合RANS方法和Realizable k-ε湍流模型對(duì)某吊艙推進(jìn)器在偏轉(zhuǎn)工況下的水動(dòng)力性能進(jìn)行預(yù)報(bào),主要包括吊艙槳推力、轉(zhuǎn)矩系數(shù),吊艙單元推力系數(shù)和橫向力系數(shù).偏轉(zhuǎn)角、推力和轉(zhuǎn)矩的定義見圖12.從船尾向船首看,吊艙推進(jìn)器向左舷偏轉(zhuǎn)偏轉(zhuǎn)角正,向右舷偏轉(zhuǎn)為負(fù),偏轉(zhuǎn)角以θ表示.

圖12 偏轉(zhuǎn)角、推力和轉(zhuǎn)矩的定義

進(jìn)速系數(shù),吊艙槳推力系數(shù)和轉(zhuǎn)矩系數(shù),吊艙單元推力和橫向力系數(shù)定義如下.

(11)

(12)

(13)

(14)

(15)

式中:Tp為吊艙槳的推力,N;Qp為吊艙槳的轉(zhuǎn)矩,N·m;Fx和Fy分別為吊艙單元沿x軸和y軸的推力和橫向力,N;D為吊艙槳直徑,m;n為吊艙槳轉(zhuǎn)速,r/min.

J=0.6時(shí),吊艙推進(jìn)器在±5°,±10°,±15°偏轉(zhuǎn)角下的水動(dòng)力性能試驗(yàn)數(shù)據(jù)和數(shù)值計(jì)算結(jié)果見表4,兩者之間的偏差見表5,吊艙槳推力系數(shù)ktp、吊艙槳10kq、吊艙單元推力系數(shù)ktx和吊艙單元橫向力系數(shù)kty的計(jì)算值與試驗(yàn)值的對(duì)比見圖13.

表4 吊艙推進(jìn)器偏轉(zhuǎn)工況下水動(dòng)力性能的試驗(yàn)數(shù)據(jù)

表5 吊艙推進(jìn)器偏轉(zhuǎn)工況下水動(dòng)力性能的數(shù)值計(jì)算結(jié)果

表6 計(jì)算值與試驗(yàn)值的偏差

由表4可見,數(shù)值計(jì)算值與試驗(yàn)數(shù)據(jù)整體上吻合較好.吊艙槳推力系數(shù)和轉(zhuǎn)矩系數(shù)的計(jì)算值與試驗(yàn)數(shù)據(jù)非常接近,誤差基本上在±3%以內(nèi);吊艙單元推力系數(shù)數(shù)值計(jì)算值整體上較試驗(yàn)數(shù)據(jù)低;吊艙單元橫向力系數(shù)的計(jì)算值在右舵角時(shí)較試驗(yàn)數(shù)據(jù)大,而左舵角時(shí)較試驗(yàn)數(shù)據(jù)小.

從圖13可得出以下結(jié)論:(1) 吊艙槳推力系數(shù)、轉(zhuǎn)矩系數(shù)的值在0°舵角時(shí)最小,并隨偏轉(zhuǎn)角的增大而增大;吊艙推進(jìn)器左右舵角大小相同時(shí),吊艙槳推力系數(shù)、轉(zhuǎn)矩系數(shù)大小基本相同;(2) 吊艙單元推力系數(shù)值在0°舵角時(shí)最大,其值隨偏轉(zhuǎn)角增大而減小;(3) 吊艙單元橫向力系數(shù)在0°舵角時(shí)最小,并隨偏轉(zhuǎn)角增大而增大;圖13c)中其數(shù)值曲線關(guān)于(0,0)點(diǎn)并不旋轉(zhuǎn)對(duì)稱,主要是吊艙槳產(chǎn)生橫向力及其后方尾流旋轉(zhuǎn)所致.

圖13 偏轉(zhuǎn)工況下吊艙推進(jìn)器水動(dòng)力性能計(jì)算值與試驗(yàn)值對(duì)比

6結(jié)論

1) 對(duì)敞水工況下4381型、4119型螺旋槳的流場(chǎng)進(jìn)行數(shù)值模擬,發(fā)現(xiàn)兩型螺旋槳敞水性能預(yù)報(bào)結(jié)果與試驗(yàn)值吻合較好,因此應(yīng)用重疊網(wǎng)格能夠獲得較高的螺旋槳敞水性能計(jì)算精度;

2) 分別使用重疊網(wǎng)格和滑移網(wǎng)格劃分流場(chǎng),在對(duì)4381型螺旋槳敞水性能預(yù)報(bào)結(jié)果的對(duì)比中,顯示出重疊網(wǎng)格在低進(jìn)速系數(shù)時(shí)對(duì)螺旋槳敞水性能預(yù)報(bào)的優(yōu)越性.

3) 使用重疊網(wǎng)格劃分流場(chǎng),預(yù)報(bào)了吊艙推進(jìn)器偏轉(zhuǎn)工況下的水動(dòng)力性能,計(jì)算值與試驗(yàn)值吻合較好.所用網(wǎng)格在吊艙偏轉(zhuǎn)后,只需再次建立吊艙子域網(wǎng)格與背景固定子域網(wǎng)格的插值關(guān)系,即可開始計(jì)算,無需重新生成,計(jì)算效率得以提升.

參 考 文 獻(xiàn)

[1]STEGER J L, DOUGHERTY F C, BENEK J A. A chimera grid scheme[C]. Mini Symposium on Advances in Grid Generation,Houston: ASME,1982.

[2]CAI Jinsheng, TSAI Hermann, LIU Feng. A parallel viscous flow solver on multi-block overset grids[J]. Computers Fluids,2006,35:1290-1301.

[3]XU Heyong, YE Zhengyin. Numerical simulation of rotor-airframe aerodynamic interaction based on unstructured dynamic overset grid[J]. Sci China Tech Sci, 2012,55(10):2798-2807.

[4]伍貽兆,田書玲,夏健.基于非結(jié)構(gòu)動(dòng)網(wǎng)格的非定常流數(shù)值模擬方法[J].航空學(xué)報(bào),2011,32(1):15-26.

[5]LI Yuwei, PIAK Kwangjun, TAO Xing, et al. Dynamic overset CFD simulation of wind turbine aerodynamics[J]. Renewable Energy,2012,37:285-298.

[6]趙發(fā)明.重疊網(wǎng)格在船舶CFD中的應(yīng)用研究[J].船舶力學(xué),2011,15(4):332-341.

[7]PABLO M C, ALEJANDRO M C, FREDRICK S. Self-propulsion computations using a speed controller and a discretized propeller with dynamic overset grids[J]. J Mar Sci Technol,2010,15:316-330.

[8] ISLAM M F, VEITCH B, AKINTURK A, et al. Experiments with podded propulsors in static azimuthing conditions[C]. 8th Canadian Marine Hydromechanics and Structures Conference,2007.

Numerical Simulation on Hydrodynamic Performance of Podded Propulsor in Azimuthing Conditions with Overset Mesh

XU Jiaqi1)XIONG Ying1)SHI Lipan2)

(DepartmentofShipEngineering,NavalUniversityofEngineering,Wuhan430033,China)1)

(NavalAeronauticalEngineeringInstitute-QingdaoCampus,Qingdao266000,China)2)

Abstract:In order to numerically simulate the hydrodynamic performances of podded propulsors in azimuthing conditions with more flexibility and higher efficiency, the computational flow domain of DTMB4381, DTMB4119 and a podded propulsor are partitioned by using overset mesh. The hydrodynamic performances of DTMB4381 and DTMB4119 are numerically predicted, and a comparison is made between the calculation results using both overset and sliding meshes. In addition, the feasibility of applying the overset mesh approach to the calculation of propellers' hydrodynamic performance is proved. The hydrodynamic performances of a podded propulsor in azimuthing conditions are predicted and analyzed. The results show that: the calculation results by using overset mesh have higher accuracy, and the overset mesh can more efficiently deal with problems of multi-body relative motions such as podded propulsors in azimuthing conditions.

Key words:overset mesh; propeller; podded propulsor; azimuthing condition; hydrodynamic performance; numerical simulation

doi:10.3963/j.issn.2095-3844.2016.02.023

中圖法分類號(hào):U661.31

收稿日期:2016-02-14

徐嘉啟(1991- ):男,碩士生,主要研究領(lǐng)域艦船流體動(dòng)力性能

*國家自然科學(xué)基金項(xiàng)目資助(51179198)

主站蜘蛛池模板: 中文精品久久久久国产网址 | 欧美区一区二区三| 亚洲国产精品不卡在线| 国产手机在线小视频免费观看| 亚洲一区国色天香| 日韩最新中文字幕| 狠狠色婷婷丁香综合久久韩国 | 伊人久久久久久久| 日本一区高清| 天堂亚洲网| 亚洲福利视频一区二区| 伊人欧美在线| 国产永久免费视频m3u8| 国产成人综合网| 日韩午夜片| 一区二区三区四区日韩| 成人亚洲天堂| 亚洲精品福利视频| 美女视频黄又黄又免费高清| 97国产在线播放| 伊人无码视屏| 国产国拍精品视频免费看| 亚洲国产系列| 国产波多野结衣中文在线播放| 精品久久久久久久久久久| 免费xxxxx在线观看网站| 亚洲男人的天堂网| 天堂在线视频精品| 伊人91在线| 中文字幕丝袜一区二区| 99国产精品免费观看视频| 在线视频亚洲色图| 99这里只有精品在线| 国产真实乱子伦视频播放| 久久精品中文无码资源站| 欧洲熟妇精品视频| 制服丝袜亚洲| 日本精品视频一区二区| 青青青国产精品国产精品美女| 亚洲人精品亚洲人成在线| 欧美综合区自拍亚洲综合绿色| 91精品啪在线观看国产91| 色综合婷婷| 亚洲av片在线免费观看| 欧美一级高清视频在线播放| 中文字幕乱码中文乱码51精品| 午夜精品国产自在| 青青热久麻豆精品视频在线观看| 久综合日韩| 无码人妻热线精品视频| 毛片视频网址| 欧美日本在线播放| 国产主播一区二区三区| 国产99精品视频| 亚洲日韩久久综合中文字幕| 日本伊人色综合网| 国产一级α片| 亚洲av中文无码乱人伦在线r| 91在线无码精品秘九色APP| 欧美精品1区2区| 亚洲中文无码av永久伊人| 91免费国产在线观看尤物| 99视频国产精品| 在线精品亚洲国产| 亚洲成人动漫在线观看| 国产va欧美va在线观看| 国产激情无码一区二区免费| 亚洲精品自拍区在线观看| 露脸国产精品自产在线播| 在线观看亚洲人成网站| 久久精品电影| 欧美狠狠干| 精品少妇三级亚洲| a天堂视频| 国产成人91精品| 国产精品无码作爱| 国产91色| 欧美一级大片在线观看| 丁香五月婷婷激情基地| 国产成在线观看免费视频| 波多野结衣无码视频在线观看| 亚洲精品日产AⅤ|