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

“鉆石背”前后翼高度差對(duì)氣動(dòng)特性影響的數(shù)值模擬研究*

2015-05-08 07:42:32陳軍葵王志軍張連博張娟婷吳國(guó)東
關(guān)鍵詞:模型

陳軍葵,王志軍,張連博,張娟婷,吳國(guó)東

(1 中北大學(xué)機(jī)電工程學(xué)院, 太原 030051; 2 93705部隊(duì), 河北唐山 064200;

3 北京航空航天大學(xué)宇航學(xué)院, 北京 100191; 4 中北大學(xué)電子測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 太原 030051)

“鉆石背”前后翼高度差對(duì)氣動(dòng)特性影響的數(shù)值模擬研究*

陳軍葵1,2,王志軍1,張連博3,張娟婷4,吳國(guó)東1

(1 中北大學(xué)機(jī)電工程學(xué)院, 太原 030051; 2 93705部隊(duì), 河北唐山 064200;

3 北京航空航天大學(xué)宇航學(xué)院, 北京 100191; 4 中北大學(xué)電子測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 太原 030051)

在“鉆石背”形翼幾何外形及設(shè)計(jì)參數(shù)的基礎(chǔ)上,選用NACA6411翼型,建立二維計(jì)算模型進(jìn)行數(shù)值模擬,研究前后翼垂直高度差ΔH對(duì)低雷諾數(shù)下“鉆石背”形翼氣動(dòng)特性的影響。結(jié)果表明,采用Transition SST湍流模型所得結(jié)果具有高的準(zhǔn)確性和可信性;ΔH對(duì)氣動(dòng)特性的影響,在前后翼前緣間距約為前翼弦長(zhǎng)的1.5倍附近,出現(xiàn)變化趨勢(shì)反轉(zhuǎn)的現(xiàn)象;優(yōu)化氣動(dòng)布局的有效方法之一是合理選取ΔH和前后翼梢弦前緣間距d的設(shè)計(jì)值,在計(jì)算范圍內(nèi),取ΔH=22 mm和d=100 mm,既能達(dá)到優(yōu)化效果,又便于工程實(shí)際應(yīng)用。

“鉆石背”形翼;湍流模型;垂直高度差;氣動(dòng)特性;數(shù)值模擬

0 引言

具有“鉆石背”氣動(dòng)布局的折疊方式,因其有收攏展開機(jī)構(gòu)簡(jiǎn)單、牢固、可靠;收攏狀態(tài)緊湊,占用空間小;升阻比大,滑翔增程能力強(qiáng);能改善大展弦比折疊翼的強(qiáng)度和剛度,提高顫振臨界速度等優(yōu)點(diǎn)[1],常常在設(shè)計(jì)彈翼或機(jī)翼需要折疊的飛行器時(shí)被考慮。目前,對(duì)于“鉆石背”氣動(dòng)布局的研究,實(shí)驗(yàn)與數(shù)值模擬同步進(jìn)行,多以加裝在高速巡飛彈或制導(dǎo)炸彈上定義為滑翔增程的“特殊”彈翼,以整體彈翼的氣動(dòng)特性研究為主。雷娟棉、吳小勝等人的研究成果表明[2-5],前翼?xiàng)l高、后翼?xiàng)l低配置的“鉆石背”彈翼的升力、升阻比大于前翼?xiàng)l低、后翼?xiàng)l高配置的“鉆石背”彈翼的升力和升阻比;“鉆石背”彈翼的外形滾轉(zhuǎn)阻尼力矩系數(shù)較大;具有“鉆石背”彈翼的飛行器不宜在大攻角下飛行。

低速和小尺度共同決定了小型飛行器的飛行雷諾數(shù)很低(1.5×104~5×105),而在低雷諾數(shù)條件下,黏性效應(yīng)和非定常效應(yīng)顯著,固定翼流場(chǎng)結(jié)構(gòu)和氣動(dòng)特性與高Re數(shù)顯著不同[6]。因此文中在已有高速和整體氣動(dòng)特性研究的基礎(chǔ)上,借助FLUENT平臺(tái),從二維角度入手,以不同垂直高度差和不同展向位置處的前后翼剖面作為仿真模型,采用Transition SST湍流模型進(jìn)行數(shù)值模擬,研究低雷諾數(shù)下前后翼垂直高度差對(duì)“鉆石背”形翼氣動(dòng)特性的影響,旨在為“鉆石背”形翼在小型無人機(jī)上的應(yīng)用和設(shè)計(jì)提供參考。

1 “鉆石背”形翼的物理及計(jì)算模型

1.1 物理模型

文中研究所采用的“鉆石背”形翼的平面形狀及幾何參數(shù)如圖1所示,其中L為左前后翼連接處和右前后翼連接處之間的距離,b為前翼根弦前緣到后翼根弦前緣的距離,d為前翼梢弦前緣到后翼梢弦前緣的距離,x為前翼任意前緣處到對(duì)應(yīng)后翼前緣處的距離(d≤x≤b),φ1為前翼后掠的角度,φ2為后翼前掠的角度,C1為前翼的有效弦長(zhǎng),C2為后翼的有效弦長(zhǎng),C0為前翼前緣到前翼后緣的垂直距離(考慮折疊后的布局問題,后翼前緣到后翼后緣的垂直距離也為C0)。此外,文中將前后翼垂直高度差ΔH定義為前翼下表面最低點(diǎn)到后翼上表面最高點(diǎn)之間的垂直距離。

圖1 “鉆石背”形翼的平面形狀及幾何參數(shù)

總體設(shè)計(jì)時(shí)首先確定b、d和C0,其他參數(shù)隨φ1、φ2的變化關(guān)系為:

C1=C0/cosφ1

(1)

C2=C0/cosφ2

(2)

(3)

值得注意的是,“鉆石背”形翼的展向長(zhǎng)度與L相關(guān),但并不由L表征,該值的大小根據(jù)翼梢的具體設(shè)計(jì)而定。

1.2 計(jì)算模型

為了研究前后翼垂直高度差ΔH對(duì)低雷諾數(shù)下“鉆石背”形翼氣動(dòng)特性的影響,選取NACA6411作為數(shù)值模擬翼型,該翼型在低雷諾數(shù)下具有良好的氣動(dòng)特性?!般@石背”形翼的幾何參數(shù)分別取為C0=130 mm,φ1=37°,φ2=11°,b=450 mm,d=0 mm,則C1=162 mm,C2=132 mm,L=950 mm。此外以x和ΔH作為主要變化量,根據(jù)前后翼的弦長(zhǎng)和前后翼前緣間距的變化及翼型厚度的不同,x以50 mm為間距遞增值,并在考慮個(gè)別特殊點(diǎn)的情況下,分別取值0 mm,15 mm,30 mm,50 mm,100 mm,150 mm,162 mm,200 mm,250 mm,350 mm,450 mm,共11個(gè)模擬工況,ΔH分別取值為0 mm,18 mm(后翼的厚度),22 mm(前翼的厚度),40 mm(前后翼厚度的和)共4類模擬工況。

計(jì)算網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格,利用Gambit軟件進(jìn)行劃分,翼型的上下及來流上游方向取25倍前翼弦長(zhǎng),翼型下游方向在x≤150 mm時(shí),取25倍前翼弦長(zhǎng),x≥162 mm時(shí),取28倍前翼弦長(zhǎng)作為計(jì)算域邊界。物面法向第一層網(wǎng)格高度為3.7×10-5C1~4.3×10-5C1,滿足y+<1。在保證計(jì)算精度的前提下,盡可能減少網(wǎng)格數(shù)量。圖2為ΔH=0 mm時(shí),x=250 mm處的“鉆石背”形翼前后翼剖面計(jì)算網(wǎng)格及翼型周圍網(wǎng)格局部圖。

圖2 ΔH=0 mm,x=250 mm時(shí),“鉆石背”前后翼剖面的計(jì)算網(wǎng)格

2 數(shù)值方法

2.1 計(jì)算方法

文中使用Fluent14.0軟件求解定常不可壓縮流動(dòng),控制方程為質(zhì)量加權(quán)的N-S方程,采用有限體積法進(jìn)行離散,選擇基于壓力的求解器,方程的求解采用SIMPLE算法,對(duì)流項(xiàng)采用二階迎風(fēng)格式,擴(kuò)散項(xiàng)采用中心差分格式。設(shè)定邊界條件為速度進(jìn)口和壓力出口,入口來流速度為34 m/s,攻角為4°,以162 mm為基準(zhǔn)弦長(zhǎng)和參考值,則雷諾數(shù)約為3.5×105。

2.2 湍流模型

Transition SST湍流模型的具體表達(dá)式及相關(guān)常數(shù)詳見參考文獻(xiàn)[11]。

3 計(jì)算結(jié)果及分析

3.1 計(jì)算方法及湍流模型驗(yàn)證

為了驗(yàn)證所采用計(jì)算方法及湍流模型的準(zhǔn)確性和可信性,針對(duì)伊利諾伊大學(xué)在低速亞音速風(fēng)洞中,關(guān)于幾種翼型進(jìn)行的低雷諾數(shù)實(shí)驗(yàn)[12-13],選取E387(E)(括號(hào)里的E是指實(shí)驗(yàn)數(shù)據(jù)采集模型中的第五種)翼型進(jìn)行數(shù)值計(jì)算和對(duì)比分析驗(yàn)證。

設(shè)定入口速度為20.25 m/s,攻角在-6°~12°之間以2°的增量進(jìn)行變化,基于弦長(zhǎng)的雷諾數(shù)Re=3.0×105。分別采用SSTk-ω湍流模型和Transition SST模型,計(jì)算所得的升力、阻力系數(shù),與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比如圖3、圖4所示。

圖3 E387(E)翼型升力系數(shù)計(jì)算值與實(shí)驗(yàn)值對(duì)比

圖4 E387(E)翼型阻力系數(shù)計(jì)算值與實(shí)驗(yàn)值對(duì)比

從圖3可以看出,采用兩種湍流模型所得的升力系數(shù)計(jì)算值與實(shí)驗(yàn)值,除了在失速攻角10°~12°之間有最大為3.4%的誤差外,其他各點(diǎn)吻合的均較好。

從圖4可以看出,采用Transition SST模型計(jì)算所得的阻力系數(shù)較SSTk-ω湍流模型更為接近實(shí)驗(yàn)值,尤其是在0°~12°之間,這種特性表現(xiàn)的更為明顯。就Transition SST模型計(jì)算所得的阻力系數(shù)與實(shí)驗(yàn)值相比較而言,在失速攻角12°附近時(shí)誤差較大,誤差為26%,其他點(diǎn)與實(shí)驗(yàn)值吻合的較好,整體呈相同的變化趨勢(shì)。

通過對(duì)比分析發(fā)現(xiàn),在-6°~10°之間的模擬點(diǎn)處,采用Transition SST模型所得的計(jì)算值較SSTk-ω湍流模型更為接近實(shí)驗(yàn)值,具有相對(duì)較高的可信度。但值得注意的是,翼型在失速位置,流動(dòng)很早于翼型后緣發(fā)生分離,且不會(huì)再附著,如圖5所示,分離后的湍流處于強(qiáng)烈的非平衡狀態(tài),而SSTk-ω湍流模型是建立在簡(jiǎn)單邊界平衡湍流基礎(chǔ)上的,它不能準(zhǔn)確預(yù)測(cè)湍流非平衡輸運(yùn)特性,這導(dǎo)致模型模擬翼型失速特性效果不佳[14]。此外,流動(dòng)進(jìn)入失速狀態(tài)后,轉(zhuǎn)捩并不是影響翼型氣動(dòng)特性模擬不準(zhǔn)的關(guān)鍵因素[15-16]。所以,兩種湍流模型在失速攻角12°附近的模擬結(jié)果均存在較大誤差。

圖5 攻角為12°時(shí),采用Transition SST模型得到的速度云圖和流線圖

因此,數(shù)值模擬求解對(duì)失速特性模擬要求不高的低雷諾數(shù)下翼型繞流問題時(shí),在上述模型建立方法的基礎(chǔ)上采用Transition SST湍流模型,所得結(jié)果具有相對(duì)較高的準(zhǔn)確性和可信性,滿足工程研究的要求。

3.2 結(jié)果分析

圖6~圖8給出了ΔH分別取0 mm、18 mm、22 mm和40 mm時(shí),升力系數(shù)、阻力系數(shù)和升阻比隨x的變化趨勢(shì)及對(duì)比圖。值得注意的是ΔH=0,x=0時(shí),該處的速度云圖和流線圖如圖9所示。從圖中可以看出,受氣流阻擋的影響,后翼很早發(fā)生大的分離,且不會(huì)再附著,這與翼型失速時(shí)的氣流狀態(tài)較像,從前面的分析可知,這種情況下所得的數(shù)據(jù),升力系數(shù)較實(shí)驗(yàn)值偏大,阻力系數(shù)較實(shí)驗(yàn)值偏小,從而使得x=0時(shí)的數(shù)據(jù)在整個(gè)變化趨勢(shì)中較為異常。考慮到其誤差較大,在后續(xù)的分析中,不再使用該點(diǎn)的數(shù)據(jù)。

3.2.1 ΔH對(duì)升力系數(shù)影響的分析

從圖6可以看出,當(dāng)ΔH取不同值時(shí),升力系數(shù)隨x的變化趨勢(shì)基本一致。當(dāng)0≤x≤250 mm時(shí),隨著ΔH的增大,同一x點(diǎn)的升力系數(shù)增大,增幅不單調(diào)變化;當(dāng)350 mm

圖6 ΔH分別取0、18、22和40 mm時(shí)的升力系數(shù)變化趨勢(shì)

3.2.2 ΔH對(duì)阻力系數(shù)影響的分析

從圖7可以看出,當(dāng)ΔH=18、22和40 mm時(shí),隨著x的增大,在0≤x≤50 mm這一段,阻力系數(shù)減小,在50 mm150 mm時(shí),隨著x的增加,4類模擬工況下的阻力系數(shù)均呈減小趨勢(shì),減幅不同。隨ΔH的增大,同一x點(diǎn)處的阻力系數(shù)在0≤x≤200 mm這一段整體呈減小趨勢(shì);當(dāng)x>200 mm時(shí),阻力系數(shù)呈增大趨勢(shì),增幅較小;阻力系數(shù)隨ΔH的變化趨勢(shì)轉(zhuǎn)折點(diǎn)位于x=200 mm處。

圖7 ΔH分別取0、18、22和40 mm時(shí)的阻力系數(shù)變化趨勢(shì)

3.2.3 ΔH對(duì)升阻比影響的分析

從圖8可以看出,受升阻力系數(shù)變化的影響,隨著ΔH的增大,在0≤x≤200 mm這一段,升阻比增大;當(dāng)x>250 mm時(shí),升阻比減小;升阻比隨ΔH的變化趨勢(shì)轉(zhuǎn)折點(diǎn)位于200 mm

圖8 ΔH分別取0、18、22和40 mm時(shí)的升阻比變化趨勢(shì)

根據(jù)以上分析,可以發(fā)現(xiàn),ΔH對(duì)氣動(dòng)參數(shù)的影響,在展向不呈現(xiàn)單調(diào)增加或減小的特點(diǎn),而是在200 mm≤x≤300 mm的區(qū)間內(nèi),亦即約為x=1.5C1附近,分別出現(xiàn)變化趨勢(shì)反轉(zhuǎn)的現(xiàn)象??傮w而言,轉(zhuǎn)折點(diǎn)之前,隨著ΔH的增大,升力系數(shù)增大,阻力系數(shù)減小,升阻比增大;轉(zhuǎn)折點(diǎn)之后,隨著ΔH的增大,升力系數(shù)減小,阻力系數(shù)增大,升阻比減小。轉(zhuǎn)折點(diǎn)之前遞增或者遞減的幅度要比轉(zhuǎn)折點(diǎn)之后遞減或者遞增的幅度大的多。

圖9 ΔH=0,x=0時(shí)的速度云圖和流線圖

此外,從圖6~圖8綜合來看,前后翼間距x和ΔH對(duì)氣動(dòng)特性的影響主要體現(xiàn)在0≤x≤1.5C1這一段,也就是說,對(duì)于低雷諾數(shù)下“鉆石背”氣動(dòng)布局的優(yōu)化,主要應(yīng)在該段做相應(yīng)的研究。合理選擇ΔH和d的設(shè)計(jì)值,即為優(yōu)化設(shè)計(jì)的有效方法之一。在文中幾何參數(shù)和計(jì)算模型的基礎(chǔ)上,選取ΔH=22 mm和d=100 mm的設(shè)計(jì)值,既能有效減小不利氣動(dòng)段在整個(gè)“鉆石背”形翼中所占的比重,達(dá)到氣動(dòng)布局優(yōu)化的目的,又能考慮到機(jī)翼折疊后所占空間問題和翼梢部位鉸接問題,滿足總體設(shè)計(jì)要求,便于工程實(shí)際應(yīng)用,較為理想。

4 結(jié)論

1)數(shù)值模擬求解低雷諾數(shù)下翼型繞流問題時(shí),如果對(duì)失速特性模擬要求不高,那么設(shè)定該流動(dòng)為定常不可壓縮,邊界條件為速度進(jìn)口、壓力出口,湍流模型采用Transition SST,所得結(jié)果具有較高的準(zhǔn)確性和可信性,可用于工程研究。

2)ΔH對(duì)低雷諾數(shù)下“鉆石背”形翼氣動(dòng)參數(shù)的影響,大約在x=1.5C1處發(fā)生變化趨勢(shì)反轉(zhuǎn)的現(xiàn)象。轉(zhuǎn)折點(diǎn)之前,隨著ΔH的增大,升力系數(shù)增大,阻力系數(shù)減小,升阻比增大;轉(zhuǎn)折點(diǎn)之后,隨著ΔH的增大,升力系數(shù)減小,阻力系數(shù)增大,升阻比減小。轉(zhuǎn)折點(diǎn)之前的增減幅度要比之后大。

3)優(yōu)化“鉆石背”形翼氣動(dòng)布局的有效方法之一是合理選取前后翼垂直高度差ΔH和前后翼梢弦前緣間距d的設(shè)計(jì)值。在文中,根據(jù)總體設(shè)計(jì)要求和工程實(shí)際應(yīng)用,選擇ΔH=22 mm,d=100 mm較為理想。

[1] 雷娟棉, 胡俊, 吳甲生. 小直徑炸彈(SDB)氣動(dòng)力問題 [J]. 彈箭與制導(dǎo)學(xué)報(bào), 2004, 24(3): 169-170.

[2] 雷娟棉, 吳甲生. 鉆石背彈翼外形參數(shù)對(duì)氣動(dòng)特性的影響 [J]. 北京理工大學(xué)學(xué)報(bào), 2006, 26(11): 945-948.

[3] 雷娟棉, 吳甲生. 鉆石背彈翼氣動(dòng)特性風(fēng)洞實(shí)驗(yàn)研究 [J]. 兵工學(xué)報(bào), 2007, 28(7): 893-896.

[4] 吳小勝, 雷娟棉, 吳甲生. “鉆石背”彈翼氣動(dòng)特性數(shù)值模擬研究 [J]. 兵工學(xué)報(bào), 2010, 31(8): 1048-1052.

[5] 吳小勝. 折疊式主彈翼氣動(dòng)特性研究 [J]. 北京理工大學(xué)學(xué)報(bào), 2010, 30(9): 1024-1027.

[6] 李鋒, 白鵬, 石文, 等. 微型飛行器低雷諾數(shù)空氣動(dòng)力學(xué) [J]. 力學(xué)發(fā)展, 2007, 37(2): 257-268.

[7] Langtry R B, Menter F R. Transition modeling for general CFD applications in aeronautics. AIAA Paper, 2005-0522 [R]. 2005.

[8] Menter F R, Langtry R B, Likki S R, et al. Acorrelation-based transition model using local variables part I-Model formulation, GT2004-53452 [R]. 2004.

[9] Menter F R, Langtry R B, Likki S R, et al. Acorrelation-based transition model using local variables part II-Test cases and industrial applications, GT2004-53454 [R]. 2004

[10] 劉沛清, 馬利川, 屈秋林, 等. 低雷諾數(shù)下翼型層流分離泡及吹吸氣控制數(shù)值研究 [J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2013, 31(4): 518-524.

[11] Fluent Incorporated. Fluent(V14.0) Theory Guides [CP]. 2011.

[12] Selig M S, Deters R W, Williamson G A. Wind tunnel testing airfoil sat Low Reynolds Numbers, AIAA 2011-875 [R]. 2011.

[13] Selig M S, McGranahan B D. Wind tunnel aerodynamic tests of six airfoils for use on small wind turbines, NREL/SR-500-34515 [R]. 2004.

[14] 文曉慶, 柳陽威, 方樂, 等. 提高k-ωSST模型對(duì)翼型失速特性的模擬能力 [J]. 北京航空航天大學(xué)學(xué)報(bào), 2013, 39(8): 1127-1132.

[15] Langtry R B, Gola J, Mente F R. Predicting 2D airfoil and 3Dwind turbine rotor performance using a transition model for general CFD codes, AIAA 2006-395 [R]. 2006.

[16] SheltonA, Abras J, Jurenko R, et al. Improving the CFD predictions of airfoils install, AIAA2005-1227 [R]. 2005.

Numerical Simulation for Influence of Vertical Height Difference to Aerodynamic Characteristics of “Diamond Back” Wing

CHEN Junkui1,2,WANG Zhijun1,ZHANG Lianbo3,ZHANG Juanting4,WU Guodong1

(1 School of Mechatronics Engineering, North University of China, Taiyuan 030051, China; 2 No.93705 Unit, Hebei Tangshan 064200, China; 3 School of Astronautics, Beijing University of Aeronautics and Astronautics, Beijing 100191, China; 4 National Key Laboratory for Electronic Measurement Technology, North University of China, Taiyuan 030051, China)

Based on geometry and design parameters of “diamond back” wing, influence of vertical height difference of the anterior and posterior wing aerodynamic characteristics at low Reynolds number were simulated by two-dimensional computational model which use NACA6411 airfoil. The simulation results indicate that the Transition SST turbulence model high accuracy and creditability in solving flowing around an airfoil at low Reynolds number; Near the station of leading edge distance between the anterior and posterior wing is about 1.5 times that of chord length of anterior airfoil, the influence of vertical height difference aerodynamic characteristics appear the phenomenon that change trend reversal; One of effective optimization methods for “diamond back” aerodynamic configuration is designing the vertical height difference and the leading edge distance between the tip chord of the anterior and posterior wing rationally, if these values are 22 mm and 100 mm, it both can achieve optimization and is convenient for engineering application.

“diamond back” wing; turbulence model; vertical height difference; aerodynamic characteristics; numerical simulation

2014-04-30

陳軍葵(1989-),男,甘肅秦安人,碩士研究生,研究方向:彈箭飛行仿真技術(shù)。

V211.3

A

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产精品对白刺激| 亚洲成a∧人片在线观看无码| 亚洲aaa视频| 亚洲无码视频一区二区三区| 午夜视频免费一区二区在线看| 久久精品国产亚洲AV忘忧草18| 亚洲中文字幕日产无码2021| 国产亚洲精久久久久久无码AV| 国产精品偷伦在线观看| 中文字幕第4页| 少妇露出福利视频| 成人伊人色一区二区三区| 2018日日摸夜夜添狠狠躁| 亚洲无线一二三四区男男| 亚洲欧美日韩精品专区| 国产在线一二三区| 2021国产乱人伦在线播放| 亚洲国产av无码综合原创国产| 欧美精品不卡| 特级做a爰片毛片免费69| 国产精品xxx| 国产成人免费| 亚洲人成网站色7777| 日韩欧美91| 午夜性爽视频男人的天堂| 激情无码字幕综合| 亚洲熟女中文字幕男人总站 | 中文字幕久久亚洲一区| 国产成人乱码一区二区三区在线| 国产偷倩视频| 国产午夜精品一区二区三区软件| 孕妇高潮太爽了在线观看免费| 1024国产在线| 成人福利在线免费观看| 欧美一区二区自偷自拍视频| 99久久国产综合精品女同| 这里只有精品在线播放| 夜夜操国产| 国产欧美精品一区二区| 精品国产网站| 天天干伊人| 午夜电影在线观看国产1区| 久久久久亚洲av成人网人人软件 | 欧美日本一区二区三区免费| 日本在线国产| 午夜国产大片免费观看| 国产在线视频导航| 国产亚洲一区二区三区在线| 午夜福利无码一区二区| 88国产经典欧美一区二区三区| 久草青青在线视频| 亚洲精品麻豆| 国产精品护士| 青青草原国产一区二区| 99热免费在线| 欧美国产日韩在线播放| 精品无码日韩国产不卡av| 欧美精品亚洲精品日韩专区va| 婷婷开心中文字幕| 色老二精品视频在线观看| 亚洲手机在线| 有专无码视频| 亚洲男人的天堂在线观看| 伊人久久婷婷| 日韩欧美色综合| 欧美日韩成人| 91亚洲影院| 一级成人a毛片免费播放| 国产大片黄在线观看| 日本国产在线| 一区二区午夜| 午夜一级做a爰片久久毛片| 婷婷综合在线观看丁香| 亚洲全网成人资源在线观看| 亚洲第一成年人网站| 正在播放久久| 欧美亚洲欧美区| 久久综合激情网| 亚洲一区网站| 3p叠罗汉国产精品久久| av无码久久精品| 日本午夜视频在线观看|