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

跨聲速面積律的近場機理研究

2016-07-05 12:53:35王鋼林
實驗流體力學 2016年4期
關鍵詞:設計

王鋼林,鄭 遂

跨聲速面積律的近場機理研究

王鋼林*,鄭 遂

(中國航空研究院飛行物理研究中心,北京 100012)

面積律過于定性的描述給實際的飛機設計工作帶來了一定的困惑和問題,其理論推導采用的小擾動線化假設也不適應未來空氣動力學設計越來越精細化的發展方向。針對具有典型高速飛行器外形特征的AGARD-B標模,結合CFD和優化方法,探討了實現最優減阻效果的機身修形形式,得出了較經典跨聲速面積律減阻效果更好的結果,給出了比經典面積律更為細致的減阻修形原則。以此為基礎,通過對各部件的減阻貢獻情況的分析,通過修形前后機體表面阻力、壓強及等壓線分布的對比,發現面積律減阻的實質是飛行器外形所造成的相鄰部件之間的壓力傳遞而形成的有利干擾。應用這一結論,研究并驗證了機身收縮剖面形狀對于減阻效果的影響。最后經過不同升力系數條件的對比,證明對于不同升力、不同迎角的飛行條件,面積律減阻的效果是相同的。

面積律;流動機理;飛機設計;空氣動力學;計算流體力學;干擾阻力

0 引 言

面積律的發現源于R.T.惠特科姆于1952年所進行的風洞實驗,根據實驗的現象,惠特科姆提出了跨聲速面積律,而后在理論上經過旋成體的小擾動線化理論分析得到證明[1-2]。

面積律作為上世紀50年代空氣動力學領域最偉大的發現之一,有力地推動了航空科學技術的發展,成功地使飛機的飛行跨越了聲速,此后各種以跨聲速或者超聲速飛行的飛行器的設計、研制,無一不遵循面積律的原則[3-4]。

但是跨聲速面積律的理論證明是基于小擾動線性化假設條件下完成的,隱藏在面積律背后的流動機理卻一直沒有得到很明確的研究和闡述,造成面積律的描述過于定性,在某些情況下給實際的飛機設計工作帶來了一定的困惑和問題。

另一方面,現在及未來飛機設計的要求在不斷提高,氣動力設計手段也從小擾動理論提升為CFD方法,空氣動力學設計朝著越來越精細化的方向發展,對于流場和流動現象的分析越來越細致[5-7]。相比較而言,建立在小擾動理論基礎上的面積律在設計空氣動力學迅速發展的今天則略顯“粗糙”[8],在對未來飛行器氣動布局及外形設計的研究中逐漸顯露其指導性略顯不足的問題[9-12]。為此一些研究者轉而試圖采用優化等方法解決跨聲速、超聲速設計方案的減阻問題[13-19];同時另外一些研究者則試圖從理論角度探討波阻的減小方法[20-26]。

本文的研究從跨聲速面積律出發,結合CFD方法和優化方法,分析面積律減阻的流動機理,尋找面積律的本質,探討面積律在更深和更精細的層次上對飛機設計的指導原則。同時探索利用CFD方法分析和解釋風洞試驗中發現的試驗現象,結合試驗和CFD方法總結相關規律和法則,并應用于指導飛機的設計與優化。

1 經典面積律及其應用

跨聲速面積律描述為:當飛行馬赫數接近于1時,為了使波阻最小,飛行器所有部件的橫截面積疊加在一起的分布應該相當于一個當量旋成體橫截面積的分布,或者分布曲線比較光滑而無不規則的變化。如圖1所示,(1)為實際飛行器,(2)為當量旋成體,(3)為實際飛行器和當量旋成體的橫截面積分布曲線,(a)為不滿足跨聲速面積律的面積分布形式,(b)為滿足跨聲速面積律的面積分布形式。

圖1 跨聲速面積律Fig.1 Transonic area law

面積律的描述簡單,很快即被應用于飛機設計實踐。美國的YF-102戰斗機在1954年試飛時由于跨聲速波阻過大而未超過聲速,后來通過其改型YF-102A的機身蜂腰設計達到了跨聲速面積律所需要的橫截面積分布,同年試飛時順利地突破音障,成為第一架采用面積律設計的飛行器。其后世界上各種以接近聲速或者超聲速飛行的飛行器無一不應用面積律設計。

2 面積律設計中的問題

面積律自提出后很快被成功應用于飛機設計,但是經由實驗現象總結、僅在小擾動線性化范圍內推導而得的面積律只有定性的闡述和粗略的定量描述,因此在飛機設計實踐中也給設計人員造成了一些困惑,帶來了一些困難。

首先,面積律未能給出其減阻的流動近場機理。近場機理的缺失使得應用面積律進行氣動外形設計的過程顯得過于“粗獷”,不利于精細化的氣動力設計,不易獲得最優的減阻效果。

其次,根據實驗結果總結和小擾動線化理論推導的面積律僅適用于細長體,而對于非細長體如何降低其波阻則沒有直接指導意義。

第三,面積律要求橫截面積分布曲線光滑,在實際工程中常用的外形修形方式是形如圖1(b)所示的收縮機身形成蜂腰。蜂腰減小中機身截面積將嚴重影響中機身的容積,而中機身又是容納飛機燃油、有效載荷等裝載的最佳場所。這樣對于其他一些飛機功能需求而言則是非常不利的。

3 研究方法

為了獲得面積律減阻的近場機理,本文選用了具有明確的參數化外形描述的AGARD-B標模[27],如圖2所示。該標模外形簡單,其主要參數如表1所示,AGARG-B標模具備典型的高速飛行器外形特征,用于研究面積律非常合適。

實際采用的計算數模如圖3所示,模型尾部增加了收縮尾錐(但在結果積分時不考慮這一段的貢獻)。同時為了簡化分析過程,主要針對馬赫數等于1.0、迎角為0°的狀態進行分析研究,最后結合跨聲速飛行的典型升力系數狀態進行驗證。

圖2 AGARD-B標準模型參數化外形定義Fig.2 Parameter definitions of standard model AGARD-B

表1 AGARD-B標準模型參數Table 1 Parameters of standard model AGARD-B

圖3 CFD計算采用的幾何數模Fig.3 Geometry model for CFD computation

由于面積律減阻的對象是波阻,因此在算法方面本文選擇了歐拉方法,在獲得本文分析所需的足夠計算精度的同時,避免RANS方法需要的巨大計算資源和時間。馬赫數1.0時,AGARD-B標模試驗數據的與本文所采用的CFD方法計算結果的升力、阻力對比情況分別如圖4和5所示。升力的計算結果與試驗完全吻合,阻力的計算結果較試驗略小,二者之間的差量為摩擦阻力,不影響本文對于面積律減小激波阻力的分析。

圖4 本文方法計算的升力與試驗結果的對比Fig.4 Lift comparison between test and computation

迎角0°時,計算數模的橫截面積分布如圖6所示。

按照經典的跨聲速面積律的要求,欲使波阻最小,需要橫截面積分布如同圖中黑色虛線所示,而由于機翼的存在,面積分布出現一個鼓包。因此在保持機翼不變的情況下,必須縮小機翼所在位置機身的直徑,使得此處單獨機身的橫截面積分布出現如同圖中紅色虛線所示的“凹坑”。為了簡化問題的復雜程度,“凹坑”深度初值取機翼最大橫截面積,在“凹坑”最低點與機翼外露根弦前后緣之間,機身橫截面積線性變化,形成圖6中虛線圍成的機身修形三角區,對應的CFD分析數模如圖7所示。

圖5 本文方法計算的阻力與試驗結果的對比Fig.5 Drag comparison between test and computation

圖6 計算數模橫截面積分布Fig.6 Cross section distribution of geometry model

圖7 面積律修形后的計算數模Fig.7 Geometry model modified with transonic area law

基于前文的設定,采用優化設計的自動程序,首先固定修形三角區下頂點高度坐標,搜索阻力最小的軸向位置(以下簡稱第1種搜索),然后將下頂點固定在最佳軸向位置,搜索阻力最小的法向位置(以下簡稱第2種搜索),從而得到機身修形的最佳方案,與跨聲速面積律結果進行對比分析,探討跨聲速面積律減阻的近場流動機理。

4 面積律的近場流動分析

4.1機身修形方式對減阻效果的影響

為了研究機身收縮程度以及最大收縮位置對于減阻效果的影響,分析過程中改變了一系列機身最大收縮直徑及其縱向位置(沿機身軸線方向)進行了優化計算。

首先根據AGARD-B標模的最大橫截面積,按照圖6黑色虛線所示的面積律要求確定機身軸對稱最大收縮直徑(最大收縮面積),采用第1種搜索方式,經過CFD計算之后得出相對應的一系列阻力值,擬合之后得到阻力隨機身軸對稱最大收縮位置之間的變化關系,如圖8所示。

圖8 阻力隨機身最大收縮直徑所在位置的變化Fig.8 Drag varies with fuselage location of maximum contraction diameter

由此可見,隨著收縮機身最小直徑所在機身軸線位置從機翼安裝根弦前緣向后緣的移動,機身修形之后的阻力值逐漸減小,直到某個位置阻力達到最小值,與修形前的等直機身相比,阻力減小約0.0041;而后隨著最大收縮直徑位置的繼續后移,阻力值則不降反升。

需要注意的是,由圖8的阻力變化規律可以發現一個經典面積律未能給出的現象,即:機身按照面積律原則進行橫截面積收縮時,收縮最大的最優位置并不位于機翼導致的橫截面積增量最大的位置(如圖8中紅色豎線所示位置),而在于該位置之后的某處。相對于在橫截面積增量最大位置進行的最大收縮構型,圖8給出的最優收縮構型還能獲得額外的約33%(約0.001)的減阻收益。由此可以得出結論,在應用面積律對飛行器進行跨聲速減阻設計時,如果將機身橫截面積最大收縮位置對應至原始橫截面積分布最突出位置的話,此時未必能獲得最佳的減阻效果。為了達到最佳的減阻效果,需要通過優化設計的方法尋找合適的機身最大收縮位置。

接下來使用第2種搜索方式研究橫截面積收縮程度對減阻效果的影響。如圖9所示,可以獲得相應的減阻效果如圖10所示。

圖10中的紅色豎線表示之前按照面積律原則確定的機身最大收縮直徑,由圖可見,在此直徑附近的阻力變化趨勢趨平,遠離此直徑時則阻力急劇增大。因此面積律所確定的機身最大修形面積可以獲得最好的減阻效果,同時在該面積值附近,減阻效果對于橫截面積的變化并不敏感。

綜合圖8和10的結果不難看出,采用優化設計的方法可以獲得與跨聲速面積律基本一致的結論。但是結合CFD方法的優化方法給出了減阻效果更好的結果,同時還揭示了經典面積律未能描述的現象。

圖9 改變機身最大收縮直徑的計算數模橫截面積分布Fig.9 Cross section distribution of geometry model for changing fuselage maximum reducing diameter

圖10 不同機身最大收縮直徑對減阻效果的影響Fig.10 Drag reduction with various minimal fuselage diameter after contraction

4.2飛行器各部件對面積律減阻的貢獻

面積律修形可以使全機獲得波阻減小的收益,但經典面積律并未指出收益從何而來。分析上一節優化得到的最優減阻構型,將其阻力按部件進行分解,可以得到如圖11所示的面積律減阻前后各部件阻力貢獻對比情況。

圖11 面積律減阻前后各部件阻力貢獻對比Fig.11 Comparison of drag contributions of different components with/without area rule

按照面積律原則對機身進行收縮修形之后,對于機身而言,所受到的阻力略有增加;對于機翼,則阻力得到大幅度降低;二者疊加,最終獲得了全機的阻力減小。由此可見,此例面積律減阻是通過修改機身的外形,在付出機身阻力的輕微代價的情況下,依靠機翼獲得的大幅度減阻收益,實現了全機減阻收益。

4.3面積律減阻的近場流動機理分析

根據優化的結果,利用CFD計算所得到的流場數據,可以繪制面積律減阻前后模型表面的阻力分布分別如圖12和13所示。

圖12 等直機身的計算數模阻力分布Fig.12 Drag distribution without area rule

圖13 收縮機身的計算數模阻力分布Fig.13 Drag distribution with area rule

對比圖12和13可以發現,未采取面積律修形之前,機翼后緣附近存在一個較大的阻力貢獻區域;經過機身修形之后,在其余部分基本保持一致的情況下,該區域的面積大為縮小,同時強度也得到抑制,從而在宏觀上獲得了減阻效果。

進一步再繪制面積律減阻前后計算模型表面的壓力和等壓線分布分別如圖14和15所示。

圖14 等直機身的計算數模壓力及等壓線分布Fig.14 Pressure distribution and isobar without area rule

圖15 收縮機身的計算數模壓力及等壓線分布Fig.15 Pressure distribution and isobar with area rule

如圖14所示,未修形時,流經機翼表面的流場正常地在機翼前半部迎風區形成高壓區,在機翼后半部背風區形成低壓區,綜合形成了較大的壓差阻力。

修形之后,機身的收縮和恢復至原直徑的變化分別使流場出現一道膨脹波和一道激波,并且都作用到了機翼表面,干擾了機翼的正常流場。機身收縮段在前,因此膨脹波的低壓區作用在機翼迎風的前半部;機身恢復段在后,故而激波的高壓區作用在了機翼背風的后半部。機身對機翼的干擾結果形成了附加的向前的力,在一定程度上抵消了機翼正常的壓差阻力,從而實現了減阻。

由此,前文中得到的機身最大收縮位置應設計在機翼(及其他部件)所導致的橫截面積增量最大的位置之后的結論也有了合理的解釋,即:當激波產生的高壓區能夠最大面積地覆蓋機翼最大厚度之后的背風區時,面積律減阻的效果最好。為了確保這一點的實現,自然需要機身最大收縮位置位于機翼導致的橫截面積增量最大位置之后。

由此可見,通過CFD的詳細分析,可以得出造成減阻效果的面積律流動機理,從而解釋面積率減阻的根本原因;同時進一步可以據此得出能夠指導飛機設計的相關設計準則,有利于跨聲速設計空氣動力學的精細化發展。

5 面積律近場流動機理的應用及驗證

5.1機身收縮剖面形狀的設計

根據第4節的分析,面積律減阻的本質是收縮機身產生的膨脹波和激波作用在機翼上所導致的有利干擾,因此收縮的目標為產生有利的膨脹波和激波,即收縮的截面形狀對于飛行器而言是可設計的。這對于飛行器的總體布置而言,無疑在一定程度上可以緩解前文第2節中提到的最后一個問題。

為了驗證這一結論,按照前文機身最大收縮面積位置沿軸線變化的規律構造了一系列兩側收縮(保持機身高度不變)和上下收縮(保持機身寬度不變)的機身,再次計算了阻力隨機身最大收縮面積所在位置之間的變化關系如圖16所示,同時也將軸對稱收縮的結果繪制在一起,以便進行比較。

圖16 3種收縮機身導致的阻力隨機身最大收縮面積所在位置的變化Fig.16 Drag comparison for three kinds of fuselage with location of maximum contraction area

由此可見,兩側收縮和上下收縮的機身修形方式對于減阻都是有效的,兩側收縮的效果與軸對稱收縮效果相當,上下收縮效果略遜一籌。這一結論無疑對于飛行器總體設計和氣動設計而言都是非常有用的。

5.2升力和迎角對于面積律減阻的影響

前文的計算、分析和結論都是在迎角為0°(此時升力系數也為0)時得到的。為了驗證飛行器正常飛行時以上各項結論是否仍然成立,本文選擇了超聲速飛行普遍需要的0.2的升力系數再次進行了計算,得出不同升力系數下阻力隨機身最大收縮面積所在位置的變化規律如圖17所示,不同升力系數下阻力隨機身最大收縮面積位置的變化規律如圖18所示。

圖17 不同升力系數下阻力隨機身最大收縮面積所在位置的變化Fig.17 Drag comparison for varying lift with fuselage location of maximum contraction area

圖18 不同升力系數下阻力隨機身最大收縮面積位置的變化Fig.18 Variational drag for varying lift with various minimal fuselage diameter after contracion

由圖可見,不同的升力系數僅造成了面積律減阻之前的基本阻力的變化,而未影響前文研究得到的各種規律和現象。應用面積律的設計可以在0°迎角狀態進行,所得到的減阻效果在各升力系數條件下均成立,因此可以使飛機跨聲速減阻設計的工作量和復雜程度得到一定程度的緩解。

6 結 論

經過本文的研究,基本明確了飛行器面積律減阻的近場流動機理,歸納起來可以得到以下結論:

(1)經典跨聲速面積律所確定的機身橫截面積最大收縮量是減阻效果最優的,但其位置位于原始橫截面積增量最大的位置之后的某處時,才能獲得最好的減阻效果;

(2)面積律的本質是相鄰部件的有利干擾,通過適當的機身外形修形,實現機身激波后的高壓區傳遞到機翼背風面,機身膨脹波后的低壓區傳遞到機翼迎風面,形成前低后高的壓力分布,產生向前的附加力,抵消了一部分正常阻力,從而達到減阻的效果;

(3)CFD方法可以成為風洞試驗的補充,結合試驗與CFD計算分析,可以有效地對試驗中得到的一些現象進行定性和定量的分析,并得到合理的解釋。

[1]Cole J D,Cook L P.Transonic aerodynamics[M].Amsterdam:Elsevier Science Publishers,2012.

[2]Chattot J,et al.Theoretical and applied aerodynamics[M].New York:Springer Science+Business Media,2015.

[3]Gu Y Q,Fan T X,Mou J G,et al.Drag reduction technology of jet-a review[J].International Journal of Engineering Research in Africa,2015,17(7):30-42.

[4]Jameson A,Ou K.50years of transonic aircraft design[J].Progress in Aerospace Sciences,2011,47(5):308-318.

[5]David A.Computational sensitivity analysis for the aerodynamic design of supersonic and hypersonic air vehicles[R].ADA622380,2015.

[6]Michael G.Computational analysis and characterization of RC-135external aerodynamics[R].ADA561659,2012.

[7]Forbes A,Patel A,Cone C,et al.Drag prediction for supersonic hydrogen-fueled airliners[R].AIAA-2011-3968,2011.

[8]Palaniappan K,Jameson A.Bodies having minimum pressure drag in supersonic flow:Investigating nonlinear effects[J].Journal of Aircraft,2010,47(4):1451-1454.

[9]Dehpanah P,Nejat A.The aerodynamic design evaluation of a blended-wing-body configuration[J].Aerospace Science and Technology,2015,43(6):96-110.

[10]Lyu Z,Joaquim R,Martins R.Aerodynamic design optimization studies of a blended-wing-body aircraft[J].Journal of Aircraft,2014,51(5):1604-1617.

[11]Roysdon P,Khalid M.Lateral-directional stability investigation of a blended-wing-body[R].AIAA-2010-9617,2010.

[12]Utsumi Y,Obayashi S.Design of supersonic biplane aircraft concerning sonic boom minimization[R].AIAA-2010-4962,2010.

[13]Ku Y C,Rho J H,et al.Optimal cross-sectional area distribution of a high-speed train nose to minimize the tunnel micropressure wave[J].Structural and Multidisciplinary Optimization,2010,42(6):965-976.

[14]Alyanak E.Modeling for conceptual design:an aeroelastic approach[R].AIAA-2012-1425,2012.

[15]Christopher L O,Choi S.Design of optimum equivalent-area target for high-fidelity low-boom aircraft design[R].AIAA-2015-2580,2015.

[16]Khalid A,Kumar P.Aerodynamic optimization of box wing-a case study[J].International Journal of Aviation,Aeronautics,and Aerospace,2014,1(4):1-46.

[17]Berger C,Carmona K,Espinal D,et al.Supersonic bi-directional flying wing configuration with low sonic boom and high aerodynamic efficiency[R].AIAA-2011-3663,2011.

[18]Zha G,Im H,Espinal D.Toward zero sonic-boom and high efficiency supersonic flight,part i:a novel concept of supersonic bi-directional flying wing[R].AIAA-2010-1013,2010.

[19]Wintzer M,Kroo I.Optimization and adjoint-based CFD for the conceptual design of low sonic boom aircraft[R].AIAA-2012-963,2012.

[20]Waddington M,McDonald R.Development of an interactive wave drag capability for the open VSP parametric geometry tool[R].AIAA-2015-2548,2015.

[21]Guan X H.Supersonic wing-body two-level wave drag optimization using extended far-field composite-element methodology[J].AIAA Journal,2014,52(5):981-990.

[22]Feng X Q,Li Z K,et al.Research of low boom and low drag supersonic aircraft design[J].Chinese Journal of Aeronautics,2014,27(3):531-541.

[23]Berci M,Vigevano L.Sonic boom propagation with a non linear geometrical acoustic model[R].AIAA-2011-2856,2011.

[24]Berci1M,Igevano L.Sonic boom propagation revisited:a nonlinear geometrical acoustic model[J].Aerospace Science and Technology,2012,23(1):280-295.

[25]Haas A,Kroo I.A multi-shock inverse design method for lowboom supersonic aircraft[R].AIAA-2010-843,2010.

[26]Jung T,Starkey R,et al.Design methods to create frozen sonic booms via lobe balancing using aircraft components[R].AIAA-2011-2857,2011.

[27]Damljanovic D,Vitic A,et al.Testing of AGARD-B calibration model in the T-38transonic wind tunnel[J].Scientific Technical Review,2006,56(2):52-62.

Research on mechanism of transonic area rule in near field

Wang Ganglin*,Zheng Sui
(Flight Physics Research Center,Chinese Aeronautical Establishment,Beijing 100012,China)

The qualitative descriptions of the area rule bring some confusion and problems to the actual aircraft design work.The linear perturbation assumption in conventional theoretical derivations does not suit the development for more and more refined aerodynamic design in the future.For AGARD-B standard model which has typical shape characteristics of high speed aircraft,we combined the CFD with optimization methods to probe the body modification form for optimal drag reduction.From that,a better drag reduction result and more detailed modification principles of drag reduction are obtained compared to those obtained from the traditional area rule method.Based on the present principles,through the analysis of drag force felt by each component and comparison of the drag forces on the body surface before and after modification,it is found that the essence of area rule drag reduction is the advantageous interference produced among the adjacent components of the aircraft configuration.Finally,the drag reduction effects of fuselage shrinkage cross-sectional shape are studied and verified.The comparison among different lift coefficient conditions validates that the drag reduction effect of area rule is the same under various lift coefficients and angles of attack condition.

area rule;flow mechanism;aircraft design;aerodynamics;computational fluid mechanics;interference drag

V211.4

:A

(編輯:李金勇)

1672-9897(2016)04-0001-07

10.11729/syltlx20160024

2016-02-01;

2016-04-25

*通信作者E-mail:wglxy@china.com

Wang G L,Zheng S.Research on mechanism of transonic area rule in near field.Journal of Experiments in Fluid Mechanics,2016,30(4):1-6,25.王鋼林,鄭 遂.跨聲速面積律的近場機理研究.實驗流體力學,2016,30(4):1-6,25.

王鋼林(1975-),男,博士,高級工程師。研究方向:飛機總體設計、新概念飛行器、飛機設計方法及相關領域。通信地址:北京市朝陽區安外北苑2號院中國航空研究院(100012)。E-mail:wglxy@china.com

猜你喜歡
設計
二十四節氣在平面廣告設計中的應用
河北畫報(2020年8期)2020-10-27 02:54:06
何為設計的守護之道?
現代裝飾(2020年7期)2020-07-27 01:27:42
《豐收的喜悅展示設計》
流行色(2020年1期)2020-04-28 11:16:38
基于PWM的伺服控制系統設計
電子制作(2019年19期)2019-11-23 08:41:36
基于89C52的32只三色LED搖搖棒設計
電子制作(2019年15期)2019-08-27 01:11:50
基于ICL8038的波形發生器仿真設計
電子制作(2019年7期)2019-04-25 13:18:16
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
從平面設計到“設計健康”
商周刊(2017年26期)2017-04-25 08:13:04
主站蜘蛛池模板: 亚洲午夜福利在线| 看国产一级毛片| 欧美亚洲一区二区三区在线| 国产91视频免费| 欧美日韩午夜视频在线观看| 国产日韩精品欧美一区喷| 国产屁屁影院| 亚洲一级无毛片无码在线免费视频 | 国产午夜在线观看视频| 亚洲一区二区三区中文字幕5566| 亚洲国产欧美目韩成人综合| 黄色三级毛片网站| 国产精品欧美在线观看| 97综合久久| 亚洲伦理一区二区| 99在线视频免费| 国产网站一区二区三区| 永久免费av网站可以直接看的| 精品国产www| 成人一级黄色毛片| 成人福利在线视频免费观看| 国产精品视频观看裸模| 亚洲欧洲日产国产无码AV| 亚洲日韩高清在线亚洲专区| 国产亚洲精久久久久久久91| 国产免费久久精品99re不卡| 婷婷六月激情综合一区| 思思热精品在线8| 亚洲天堂伊人| 亚洲av无码人妻| 国产96在线 | 精品在线免费播放| 日韩 欧美 国产 精品 综合| 五月婷婷伊人网| 亚洲三级成人| 91在线一9|永久视频在线| 69av免费视频| 国产成人免费视频精品一区二区| 日本精品一在线观看视频| 美女无遮挡免费网站| 好久久免费视频高清| 免费毛片a| 在线免费a视频| 国产极品粉嫩小泬免费看| 欧美狠狠干| 亚洲不卡影院| AV老司机AV天堂| 亚洲免费三区| 成人看片欧美一区二区| 福利在线不卡一区| 色综合中文字幕| 久草中文网| 57pao国产成视频免费播放| 欧美精品v欧洲精品| 国产精品美女自慰喷水| 日韩精品一区二区三区中文无码| 亚洲av无码人妻| 亚洲国产中文欧美在线人成大黄瓜| 国产av一码二码三码无码| 国产96在线 | 毛片免费在线视频| 国产精品无码久久久久AV| 色综合综合网| 欧美日韩亚洲国产主播第一区| 九九精品在线观看| 亚洲午夜综合网| 国产男女XX00免费观看| 亚洲人成在线精品| 九色91在线视频| 四虎永久免费地址在线网站 | av在线人妻熟妇| 无码综合天天久久综合网| 国产喷水视频| 亚洲日本韩在线观看| 99999久久久久久亚洲| 尤物亚洲最大AV无码网站| 国产欧美一区二区三区视频在线观看| 国产人成在线视频| 园内精品自拍视频在线播放| 国产无码制服丝袜| 国产网站免费| 午夜免费视频网站|