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

民機風洞試驗半模墊板高度對氣動特性的影響

2017-11-20 03:32:51王繼明劉亦鵬
航空學報 2017年5期

王繼明, 劉亦鵬

上海飛機設計研究院, 上海 201210

民機風洞試驗半模墊板高度對氣動特性的影響

王繼明*, 劉亦鵬

上海飛機設計研究院, 上海 201210

半模作為提高大型商用飛機風洞試驗雷諾數的一種模擬手段而被廣泛應用。首先回顧了半模試驗的模擬方式及其優劣,進而選取當前發展趨勢的附面層墊板作為研究對象,采用數值模擬研究了墊板高度變化對氣動特性影響的內在機理。數值模擬結果和試驗吻合較好,數值計算采用速度分布入口可以較好模擬風洞核心段邊界層厚度,計算值和試驗值更加接近;墊板高度的增加使得升力系數增加、阻力系數減小及俯仰力矩系數增加;墊板在機翼上游區引起的上洗使得機翼沿展向各剖面當地迎角增加5%、動壓增加1%,從而使得機翼上翼面壓力分布朝負值方向移動。有別于以往認為墊板的洗流只影響內側機翼,結果表明墊板影響范圍擴展至全翼展,當地迎角的增加是主要影響因素,墊板對機翼展向各剖面影響量值不一致,對內側機翼影響較大。所得結論可更好用于民機半模風洞試驗的開展,具有一定的工程實用性。

半模; 風洞試驗; 墊板; 邊界層; 氣動特性

長期以來風洞試驗被用作獲取大型商用飛機[1]氣動數據及研究飛機氣動特性的重要工具,而風洞試驗能夠達到的雷諾數遠遠小于飛行雷諾數。雷諾數是影響飛機氣動特性尤其是失速及分離特性的重要參數。通常雷諾數的影響都有一個敏感范圍,即當雷諾數在該范圍內時氣動特性隨雷諾數變化劇烈,當雷諾數較高時氣動特性隨雷諾數變化較為緩和。因此增加試驗模擬的雷諾數或盡可能地使試驗模擬的雷諾數超過敏感范圍能獲取更加可靠的試驗數據,同時也降低了將試驗雷諾數外推到飛行雷諾數進行數據修正的風險。半模試驗是一種較為經濟的獲取較高雷諾數的手段,在相同尺寸的風洞中半模試驗的雷諾數較全模大一倍。對于目前中國4 m×3 m量級的風洞,全模雷諾數一般為1.4×106,而對于目前大型商用飛機其雷諾數敏感范圍一般在6×106以下,在雷諾數為3×106時,氣動特性變化尤為劇烈,可見半模試驗可有效避開該敏感范圍。

半模試驗可以在不增加模型加工及試驗成本的基礎上增加雷諾數,可以模擬更加細節的部位,無支架干擾;缺點是半模的機身對稱面流動和實際不符,半模機身及內側機翼與風洞洞壁的干擾較全模大,因此如何使半模對稱面、機身及內側機翼的流動更加接近實際情況是國內外研究的核心內容。

半模通常采用附面層墊板使得模型遠離風洞壁面來減小洞壁對模型的干擾,有研究[2]表明墊板厚度約在2δ*(δ*為空風洞核心段附面層位移厚度)模擬效果較好,且升力系數隨著墊板厚度的增加而增加。NASA亦有研究表明[3]墊板厚度在模型半展長的3%模擬效果較好。從各墊板高度下機翼各剖面的壓力分布和全模對比來看,墊板高度的增加使得整個上翼面流速增加,從而對應的壓力系數Cp減小,升力系數CL增加。

因此半模試驗應著重選擇合適的墊板高度,因為較高或較低的墊板高度會使得墊板影響較大或機身進入附面層影響區從而帶來模擬的差異。同時盡可能模擬對稱面及內側機翼的流動,消除或減弱機頭前方馬蹄渦[4]對下游區的影響。

1 半模試驗研究回顧

半模試驗模擬主要包括減弱邊界層對模型的影響,如反射板研究[5-6]、半模墊板高度研究[7-11]、半模墊板前緣外形研究[3,8]及減小模型對稱面處洞壁邊界層厚度如邊界層吹吸(主動控制)研究[12]。國外對半模試驗模擬較系統的研究報導可追溯到20世紀90年代初,國內報導[13-14]稍晚,研究內容與方法和國外類似。Milholen等[2]的研究結果表明墊板厚度從δ*~15δ*的變化過程中,升力系數單調增加且在墊板厚度2δ*處與全模數據吻合較好。分析其原因為墊板厚度增加使得整個上翼面流速增加,從而上翼面吸力增加,而墊板厚度小于2δ*時,上翼面流速降低,相應吸力減小。Eliasson[15]認為機翼展向交叉流是半模數據不模擬的一個重要因素。對稱面機身壓力分布[2-3]也是半模模擬的一項參考因素,墊板高度會使得機身對稱面上表面加速,吸力增加;對于沿流向剖面采用機身對稱面外形的墊板,其和洞壁邊界層相互干擾使得在墊板前緣上游處形成馬蹄渦[4],對于全模則不會有,有研究[3,8]表明,對墊板機頭處向內側倒圓角可以減弱或消除馬蹄渦的不利影響,從而使得半模和全模數據吻合較好。但也有學者[3]通過選擇合適的2D墊板高度可以獲得較好的模擬結果。前述減弱邊界層與機身干擾的方法是增加高度使得機身遠離洞壁,另外一種方法[16]是在上游吹除下游吸附以減弱邊界的影響,結果表明升力系數及俯仰力矩系數都有不同程度的增加,尤其是俯仰力矩系數增加較為明顯。

鑒于附面層墊板模擬簡單有效,當前大型商用飛機制造商波音和空客公司的主流機型都通過半模墊板方式獲取飛行雷諾數試驗數據。如波音系列Boeing737NG[17]、Boeing777[18]及Boeing787[19]等機型在美國低溫增壓風洞NTF[20]通過附面層墊板的半模模擬方式獲取高雷諾數試驗數據。空客A320等機型[21-22]在歐洲跨聲速低溫增壓風洞ETW通過附面層墊板半模試驗獲取飛行雷諾數試驗數據。

縱觀以上各種方法,其目的都是減弱邊界層對模型的影響以獲取更為準確的模擬。從發展趨勢來看,采用附面層墊板使得模型遠離洞壁的方法因簡單有效而得到廣泛應用。

2 試驗模型及風洞

本次試驗模型為下單翼翼吊常規布局民用飛機,機翼采用新一代超臨界翼型。試驗中不帶垂尾及平尾,模型比例為1∶7,總長為5 546 mm,半翼展為2 557 mm,墊板研究高度為60、80、100、120、140 mm。模型在風洞中安裝如圖1所示。墊板迷宮槽常用于減弱機身對稱面與墊板之間的串流,其尺寸如圖2所示。

承試風洞橫截面尺寸為4.5 m×3.5 m,試驗來流馬赫數Ma=0.2,常壓下以平均氣動弦長為參考尺寸的雷諾數Re=2.9×106。

圖1 風洞試驗模型 Fig.1 Test model in wind tunnel

圖2 墊板迷宮槽間隙尺寸示意圖 Fig.2 Schematic diagram of peniche labyrinth slot dimensions

3 計算網格

選取合適的邊界層網格高度對于邊界層模擬至關重要。為研究合適的邊界層網格參數,選取了第1層網格高度Δy為0.2,0.1,0.05,0.02,0.01 mm 這5個算例。增長因子相同為1.4,邊界層Prism網格總高度相同,風洞壁面采用20層Prism,模型固壁附近采用10層Prism(如圖3所示),計算域中設置最大網格尺寸為80 mm。相應網格數量為900萬、1 000萬、1 100萬、1 200萬及1 300萬。

圖4為第1層網格高度Δy對線性段升力系數的影響,Δy在0.05 mm以上,第1層網格高度對升力系數影響較大,0.2 mm(y+~40)高度較0.05 mm(y+~10)高度的升力系數大0.01。而從0.05 mm減小到0.01 mm(y+~2),升力系數變化在0.002左右,接近試驗的重復性,故認為在該網格尺寸以下,網格的影響可忽略。且隨著第1層網格尺寸的減小,CFD計算結果與試驗結果更加接近。

因此網格第1層厚度最終采用0.05 mm,網格單元數為1 100萬。

圖3 固壁面邊界層網格 Fig.3 Boundary layer mesh of solid surface

圖4 網格第1層高度對升力系數的影響 Fig.4 Effect of first node height on lift coefficient

4 洞壁邊界層模擬

由于半模接近風洞壁面,因此洞壁邊界層對模型流場的影響不可忽略。模擬風洞核心段邊界層厚度須設置合理的入口邊界條件,以確保邊界層發展到核心段和實驗值一致。湍流模型采用k-ω剪切應力輸運(SST),入口界面速度Vin考慮洞壁邊界層效應(如圖5所示),模擬入口及風洞核心段的邊界層厚度,表達式為

(1)

式中:γ=1.4為空氣比熱比;空氣氣體常數Rg=287 J/(kg·K);T為試驗時風洞來流溫度;hw為距離洞壁的高度;hin為空風洞入口邊界層厚度。

通常給定入口邊界條件: 恒定速度入口V(圖6(a))及速度分布入口Vin(圖6(b))。明顯恒定速度分布入口條件邊界層厚度是從0開始發展的,而速度分布入口則先有邊界層厚度再發展至核心段。模擬風洞核心段邊界層厚度有采用增加試驗段長度的方法,但該方法勢必增加較多網格數量。

圖5 考慮邊界層效應的入口速度設置 Fig.5 Inlet velocity setting considering boundary layer effects

圖6 風洞試驗核心段邊界層厚度模擬 Fig.6 Boundary layer thickness simulation of middle wind tunnel section

圖7為上述兩種入口邊界條件的風洞核心段邊界層內速度u分布對比,可見采用恒定速度入口條件其核心段邊界層厚度約為60 mm,與實驗值140 mm有較大差距,而采用速度入口分布條件可以較好模擬風洞核心段的邊界層厚度。

如圖8所示,從兩種入口邊界條件氣動力系數差量來看,恒定速度入口較入口速度分布Vin條件升力系數大。兩種入口邊界條件升力系數的差量為0.009。入口速度分布Vin條件更加接近于試驗值。

圖7 兩種入口邊界條件的風洞試驗核心段邊界層內速度分布對比 Fig.7 Comparison of middle wind tunnel test section velocity distribution between two inlet boundary conditions

圖8 升力系數在兩種入口邊界條件與試驗值對比 Fig.8 Comparison of lift coefficients of two inlet boundary conditions with test results

綜上,采用恒定速度入口邊界,其邊界層厚度從0開始發展,其風洞核心段邊界層厚度計算值約為60 mm,較試驗值小80 mm,得到的升力系數偏大。而采用速度分布入口Vin條件,相當于給定入口邊界層厚度使得發展至風洞核心段的邊界層厚度和試驗值一致,速度分布入口條件較好地模擬了試驗段的邊界層厚度,且氣動力系數的計算值更加接近于試驗值。本文將以速度分布作為入口條件分析墊板高度變化對大型商用民機氣動特性的影響。

5 半模墊板高度結果分析

民機的失速特性、升阻特性及力矩特性影響飛機的進場特性、爬升特性、安定性及尾翼配平和載荷特性。研究半模獲取的數據和全模的差異及民機氣動特性隨半模墊板高度的變化對于研究民機各項氣動特性至關重要。圖9為墊板高度變化對升力系數CL、阻力系數CD及俯仰力矩系數Cm影響的試驗結果。墊板高度h從60 mm增加到140 mm,升力系數增加0.02,阻力系數減小0.004 3,俯仰力矩系數增加0.03。與全模相比,半模的升力線斜率大4%,且半模的升力系數較全模的大。

墊板高度越小,半模值越接近于試驗值。但從趨勢來看,即使墊板高度減小到0,半模的結果和全模也不完全一致,但從流動分離與發展來看,半模和全模吻合較好,在該雷諾數下外翼分離在前,內翼分離在后。而由后續分析可知,當模型越靠近壁面,受到壁面影響,內翼分離會提前。故對于內翼先分離并發展引起失速的飛機,墊板高度不能太小,否則失速特性模擬不準,雖然墊板高度越低其線性段與全模越接近。

CFD模擬結果顯示,墊板高度從60 mm增加到100 mm,升力系數增加0.010,較風洞試驗結果略大。墊板高度從60 mm增加到140 mm,升力系數增加0.020,該差量和試驗結果吻合較好,如表1所示。CFD模擬結果較好地反映了墊板高度對氣動特性影響的趨勢,通過CFD流場分析可進一步獲取墊板高度變化對流場影響的機理。

為分析墊板高度增加引起飛機氣動特性變化的原因,圖10及圖11對比了墊板高度對機翼及機身壓力分布影響的數值模擬結果。從壓力分布來看,墊板高度的增加對整個翼展壓力分布都有影響,主要影響的是上翼面壓力分布,墊板越高上翼面吸力越大,這也是升力系數隨著墊板高度增加的主要原因。對影響區域進行分析,展向內側機翼壓力分布隨墊板高度變化更為顯著,弦向則30%c(c為弦長)之前變化更明顯。隨著墊板高度增加內側機翼吸力增加,相應增加了俯仰力矩。機翼壓力分布變化是全機升力系數變化的主要原因,而機身受邊界層干擾更為明顯。隨著墊板高度的增加,機身壓力分布朝負值方向移動且頭部壓力變化較中后部劇烈,因此使得阻力減小及俯仰力矩增加。

圖9 不同墊板高度氣動系數和全模對比 Fig.9 Comparison of aerodynamic coefficients with different peniche heights and full model

表1 墊板高度對線性段升力系數的影響

Table1Effectsofpenicheheightsonlinearliftcoefficients

Penicheheight/mmCLdifferenceWindtunnelCFD60?1000.0070.01060?1400.0200.020

圖10 墊板高度對機翼壓力分布的影響 Fig.10 Effect of peniche heights on wing pressure distribution

圖11 墊板高度對機身壓力分布的影響 Fig.11 Effect of peniche heights on fuselage pressure distribution

分析近壁面處(距洞壁1 mm)及機身對稱面處的流場可以發現,在近壁面處存在馬蹄渦(如圖12 所示),但該渦僅存在于墊板高度范圍內。馬蹄渦是風洞壁面邊界層與模型相互作用的產物,在全模試驗中并非存在。隨著墊板高度的增加,馬蹄渦的影響呈現先減小后增大的趨勢。馬蹄渦改變了墊板附近的流場,使得頭部上表面上洗增加,這種流場的變化也使得機身對稱面處的速度場發生變化(如圖13所示),隨著墊板高度的增加,機身對稱面處流速增加,進而影響機身及機翼的壓力分布。

墊板高度的增加使得上翼面壓力分布朝負值方向移動,升力隨之增加,而影響升力的因素主要有迎角及動壓。墊板本身在上翼面誘導的上洗流在增加當地來流速度Vlocal的同時也使得機翼當地迎角αlocal增加。通過對比分析各個墊板高度下機翼前方來流的當地迎角及流速的變化(如圖14及圖15所示),可以發現內側機翼受墊板影響最大,當墊板高度從60 mm增加到140 mm,當地迎角增加約0.6°,增加約5%;速度增加約0.4 m/s,相應動壓增加約1%。從數值來看,升力系數隨著墊板高度的增加而增加,其中當地迎角的增加是主導因素。

圖12 風洞近壁面(距洞壁1 mm)處馬蹄渦 Fig.12 Horse shoe vortex near wind tunnel wall (1 mm to tunnel wall)

圖13 機身對稱面處速度分布對比 Fig.13 Comparison of velocity distribution in symmetry plane

圖14 墊板高度對翼展剖面當地迎角的影響 Fig.14 Effect of peniche height on local angle of attack of different spanwise location

圖15 墊板高度對機翼剖面當地流速的影響 Fig.15 Effect of peniche height on local velocity of different spanwise location

6 結 論

1) 第1層網格高度采用0.05 mm,y+~10,網格數量1 100萬與試驗結果吻合較好。采用速度分布入口邊界條件可以更好地模擬風洞核心段的邊界層厚度,其計算值更加接近于試驗值。

2) 對于4.5 m×3.5 m量級風洞,超臨界機翼翼吊布局大型商用民機,墊板高度從60 mm增加到140 mm:升力系數增加0.02,阻力系數減小約0.004 3,俯仰力矩系數增加0.03。

3) 數值模擬結果表明墊板高度對升力影響的差量與試驗結果吻合較好,通過數值模擬研究半模的流場特性可知:近壁面存在馬蹄渦,隨著墊板高度的增加,馬蹄渦的影響呈現先減小后增加的趨勢,且馬蹄渦僅存在于墊板高度范圍的流場內。

4) 有別于以往認為墊板僅影響內側機翼的流動,研究表明墊板高度的增加使得整個翼展范圍的上翼面壓力分布朝負值方向移動,對內側機翼影響更大。

5) 墊板高度的增加誘導的上洗流使得機翼各剖面當地迎角及來流速度增加,其中內翼當地迎角增加近5%,內翼段來流動壓增加約1%。

6) 墊板高度的增加使得機身對稱面處流速增加,機身壓力分布朝負值方向移動,且頭部影響更大,從而使阻力系數減小及俯仰力矩系數增大。

[1] WHITE P J, GIBSON T M. Exploitation by airbus of high Reynolds number test capabilities in the European transonic wind tunnel: AIAA-2004-0768[R]. Reston: AIAA, 2004.

[2] MILHOLEN II W E, CHOKANI N, MCGHEE R J. Development of semi-span model test techniques: AIAA-1996-2412[R]. Reston: AIAA, 1996.

[3] GATLIN G M, PARKER P A, OWENS L R, Jr. Development of a semi-span test capability at the national transonic facility: AIAA-2001-0759[R]. Reston: AIAA, 2001.

[4] MILHOLEN II W E, CHOKANI N. Computational analysis of semi-span test techniques: AIAA-1995-2290[R]. Reston: AIAA, 1995.

[5] GATLIN G M, MCGHEE R J. Study of semi-span model testing techniques: AIAA-1996-2386[R]. Reston: AIAA, 1996.

[6] ZHANG F, KHALIDT M, XU H, et al. Numerical prediction of the reflection plate boundary layer effects on half plane model in wind tunnel: AIAA-2000-2376[R]. Reston: AIAA, 2000.

[7] VIEHWEGER G, EWALD B. Half model testing in the Cologne Cryogenic Tunnel (KKK): AIAA-1994-2511[R]. Reston: AIAA, 1994.

[8] MILHOLEN II W E. A design methodology for semi-span model mounting geometries: AIAA-1998-0758[R]. Reston: AIAA, 1998.

[9] GATLIN G M, MCGHEE R J. Experimental investigation of semispan model testing techniques[J]. Journal of Aircraft, 1997, 34(4): 500-505.

[10] EARNSHAW P B, GREEN A R, HARDY B C, et al. A study of the use of half-models in high-lift wind-tunnel testing[C]//High-Lift System Aerodynamics. Paris: AGARD, 1992: 20-1-20-9.

[11] TAKALLU M A. Reynolds-averaged Navier-Stokes computations of a high-lift transport model with and without semi-span standoff: AIAA-2000-4222[R]. Reston: AIAA, 2000.

[12] PETZ R, NITSCHE W. Active control of flow separation on a swept constant chord half model in a high-lift configuration: AIAA-2006-3505[R]. Reston: AIAA, 2006.

[13] 李艷亮, 董軍, 楊希明. 半模試驗翼身組合體墊塊的數值模擬分析[J]. 空氣動力學學報, 2009, 27(1): 73-77.

LI Y L, DONG J, YANG X M. Navier-Stokes simulation of the stand-off on the flow over a fuselage wing combination[J]. Acta Aerodynamica Sinica, 2009, 27(1): 73-77 (in Chinese).

[14] 李艷亮, 董軍, 楊希明. 墊塊法半模型風洞試驗的數值模擬研究[J]. 航空計算技術, 2006, 36(6): 36-39.

LI Y L, DONG J, YANG X M. Navier-Stokes simulation of influence of stand-off on flow over a semi-span wing[J]. Aeronautical Computing Technique, 2006, 36(6): 36-39 (in Chinese).

[15] ELIASSON P. Numerical validation of a half model high lift configuration in a wind tunnel: AIAA-2007-262[R]. Reston: AIAA, 2007.

[16] MALIK A, RENDER P M. Use of wall suction in half model wind tunnel testing: AIAA-2010-4828[R]. Reston: AIAA, 2010.

[17] SAUNDERS M. High Reynolds number testing of a conventional high lift model in a mild cryogenic environment: AIAA-2008-0837[R]. Reston: AIAA, 2008.

[18] GATLIN G M, TOMEK W G, PAYNE F M, et al. Recent improvements in semi-span testing at the national transonic facility: AIAA-2006-0508[R]. Reston: AIAA, 2006.

[19] PAYNE F, BOSETTI C, GATLIN G, et al. Progress in flaps down flight reynolds number testing techniques at the NTF: AIAA-2007-0751[R]. Reston: AIAA, 2007.

[20] GATLIN G M, PARKER P A, OWENS L R. Advancement of semispan testing at the national transonic facility[J]. Journal of Aircraft, 2002, 39(2): 339-353.

[21] GROSS N, QUEST J. The ETW wall interference assessment for full and half models: AIAA-2004-0769[R]. Reston: AIAA, 2004.

[22] QUEST J, WRIGHT M, HANSEN H, et al. First measurements on an airbus high lift configuration at ETW up to flight Reynolds number: AIAA-2002-0423[R]. Reston: AIAA, 2002.

(責任編輯: 李明敏)

URL:www.cnki.net/kcms/detail/11.1929.V.20160816.0859.002.html

Effectsofhalfmodelpenicheheightoncivilaircraftaerodynamiccharacteristicsinwindtunneltest

WANGJiming*,LIUYipeng

ShanghaiAircraftDesignandResearchInstitute,Shanghai201210,China

Halfmodelsimulation,asamethodtogethighertestReynoldsnumber,iswidelyusedinthedesignoflargecommercialtransportaircrafts.Thispaperreviewstheprosandconsofthehalfmodelsimulation,andthenstudiesthepenichesimulationwhichiswidelyaccepted.Themechanismoftheeffectofthepenicheheightonaerodynamiccharacteristicsisstudied.CFDsimulationisfoundtoagreewellwiththeexperimentalresult.Theboundarylayerthicknessofthemiddleofthewindtunneltestsectioncanbebettersimulatedbyusingvelocitydistributioninletcondition,andthesimulationresultsaremoreclosetotheexperimentalresults.Withtheincreaseofpenicheheight,liftcoefficientincreases,dragcoefficientdecreasesandpitchingmomentcoefficientincreases.Theupwashinducedbypenicheinthecomingflowinfrontofthewingincreasesthelocalangleofattackby5%anddynamicpressureby1%alongthefullspan,thusmakingthepressuredistributionmorenegative.Differentfromthetraditionalconceptsthattheupwashinducedbypenichecanonlyaffecttheinboardwing,resultsshowthatthepenicheeffectsextendtothewholespan.Theprimefactoristheincreaseofthelocalangleofattack.Theeffectsofpenicheheightsvarywiththespanwiselocation,havingmoreimpactsontheinboardwing.Theresultscanbebetterusedinthehalfmodelwindtunneltestwithcertainengineeringpracticability.

halfmodel;windtunneltest;peniche;boundarylayer;aerodynamiccharacteristics

2016-05-11;Revised2016-06-01;Accepted2016-08-05;Publishedonline2016-08-160859

s:AeronauticalScienceFoundationofChina(20153240003);CivilAircraftProjectResearch(MJ-2014-F-04-01)

.E-mailwangjiming@comac.cc

2016-05-11;退修日期2016-06-01;錄用日期2016-08-05; < class="emphasis_bold">網絡出版時間

時間:2016-08-160859

www.cnki.net/kcms/detail/11.1929.V.20160816.0859.002.html

航空科學基金 (20153240003); 民用飛機專項科研 (MJ-2014-F-04-01)

.E-mailwangjiming@comac.cc

王繼明, 劉亦鵬. 民機風洞試驗半模墊板高度對氣動特性的影響J. 航空學報,2017,38(5):120429.WANGJM,LIUYP.EffectsofhalfmodelpenicheheightsoncivilaircraftaerodynamiccharacteristicsinwindtunneltestJ.ActaAeronauticaetAstronauticaSinica,2017,38(5):120429.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0229

V211.753

A

1000-6893(2017)05-120429-09

主站蜘蛛池模板: 青青草91视频| 91青青草视频| AV不卡国产在线观看| 97一区二区在线播放| 国产欧美又粗又猛又爽老| 伊人精品成人久久综合| 色吊丝av中文字幕| 国产色伊人| 欧美一区日韩一区中文字幕页| 在线视频97| 国产精品yjizz视频网一二区| 亚洲人成网址| 伊人久久精品无码麻豆精品| 女同国产精品一区二区| 午夜视频日本| 人妻91无码色偷偷色噜噜噜| 久久伊人久久亚洲综合| 亚洲视频四区| 久久综合AV免费观看| 色首页AV在线| 日韩精品毛片人妻AV不卡| 亚洲精品在线91| 3344在线观看无码| 日韩 欧美 小说 综合网 另类| 精品综合久久久久久97| 美女免费黄网站| 国产成人精品高清在线| 亚洲国产精品无码AV| 国产在线八区| 精品91视频| 亚洲一级无毛片无码在线免费视频 | 国产欧美日韩综合在线第一| 国产欧美精品专区一区二区| 99国产在线视频| 久久婷婷国产综合尤物精品| 黑色丝袜高跟国产在线91| 精品夜恋影院亚洲欧洲| 国产小视频a在线观看| 久久久波多野结衣av一区二区| 国产女人在线视频| 欧美、日韩、国产综合一区| 一级香蕉人体视频| 精品亚洲欧美中文字幕在线看| 欧美另类视频一区二区三区| 9久久伊人精品综合| 污网站在线观看视频| 亚洲综合片| 18禁不卡免费网站| 亚洲国产精品一区二区第一页免 | 无码一区18禁| 亚洲精品色AV无码看| 狠狠亚洲五月天| 好紧太爽了视频免费无码| 中文字幕亚洲精品2页| 国产精品亚洲欧美日韩久久| 免费高清毛片| 亚洲精品欧美日本中文字幕 | 亚洲无码视频图片| 国产成人在线无码免费视频| 网友自拍视频精品区| 欧美性猛交xxxx乱大交极品| 91精品久久久无码中文字幕vr| 国产成人免费| 色老二精品视频在线观看| 国产精品免费电影| 无码久看视频| 国产精品成人AⅤ在线一二三四| 亚洲成人动漫在线观看| 欧美啪啪网| 成人无码一区二区三区视频在线观看 | 亚洲一区二区三区麻豆| 三级国产在线观看| 久久一级电影| 伊人精品视频免费在线| 女人18毛片久久| 国产精女同一区二区三区久| 欧美一级片在线| 久久精品午夜视频| 久久久久免费精品国产| 国产特级毛片aaaaaa| 亚洲一区色| 日韩美毛片|