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

基于激波裝配法的乘波體設計與分析

2017-07-03 16:08:48陳冰雁劉傳振紀楚群
空氣動力學學報 2017年3期
關鍵詞:設計

陳冰雁, 劉傳振, 紀楚群

(中國航天空氣動力技術研究院, 北京 100074)

?

基于激波裝配法的乘波體設計與分析

陳冰雁, 劉傳振*, 紀楚群

(中國航天空氣動力技術研究院, 北京 100074)

乘波體外形通過在給定的高超聲速激波流場中使用流線追蹤法設計得到,傳統的乘波體設計一般采用二維楔形或者軸對稱錐形流場生成,設計空間受限。為拓展乘波體設計空間,引入激波裝配法,數值計算帶激波的流場,建立了普適于一般三維流場的乘波體設計方法,突破了傳統乘波體設計方法的流場限制。引入“導波體”定義三維流場的生成外形,采用平切導波體流場激波及正向追蹤流線的方法設計乘波體外形。通過幾類導波體流場生成乘波體的算例,分析了導波體與乘波體之間的關系,包括縱向截面外形、橫截面外形以及前緣外形等。分析表明,雙錐體在后體產生的激波或膨脹波能有效地改變導波體速度沿縱向的分布特性,利用該特性可在俯視平面形狀受約束的情況下實現對乘波體構型縱向穩定性的按需設計;通過改變導波體橫截面外形,可以改變速度沿展向的分布特性,從而實現乘波體橫截面外形改變;乘波體前緣外形主要與導波體激波在俯視面上的外形相關,通過后體擴張的雙錐外形能夠生成雙前緣乘波體外形。本文指出乘波體與導波體外形之間存在定性對應關系,其關系可作為乘波體設計時導波體選擇的參考依據,為擴大設計空間、設計優良乘波體奠定了基礎。

激波裝配法;三維流場;乘波體;氣動外形;數值模擬

0 引 言

乘波體設計是當前國際上高超聲速飛行器氣動布局研究的重點和熱點之一。乘波體構型在設計狀態將激波附著于下表面前緣,阻止了氣流泄露,具有很高的升阻比。升阻比是反應高超聲速飛行器氣動特性非常重要的一個參數,升阻比很大程度上決定了高超聲速飛行器能夠實現哪一種飛行軌道。對于滑翔式和巡航式高超聲速飛行器的氣動設計,升阻比往往是非常重要的一個技術指標,因為升阻比直接關系到飛行器能夠達到的航程、橫向機動能力等關鍵戰技術指標。因此,乘波體往往被選擇為高超聲速飛行器氣動布局的基礎構型。

乘波體的設計要素包括兩點:帶激波的流場與沿流線追蹤生成型面。本文定義生成激波流場的外形為“導波體”,與在流場中追蹤流線得到的外形“乘波體”相對應。傳統的乘波體設計方法均可以據此定義導波體,例如楔形流場乘波體[1]的導波體為二維楔形,錐導乘波體[2-3]對應圓錐,密切錐方法[4-5]的流場是由已知激波形狀反推得到,基準流場為錐形流,導波體也可以認為是錐體的一種?;鶞柿鲌鰧Τ瞬w的性能有著根本性的影響,選擇導波體就是選擇基準流場,是乘波體設計過程的重要環節,也是保證乘波體滿足氣動、裝載等設計要求的基礎。

圓錐的激波流場一般通過求解Taylor-Maccoll方程數值[6],在乘波體設計中應用較多,但這也限制了導波體的選擇范圍,設計空間不大,同時流場模擬不夠準確[7]。為了規避傳統乘波體設計方法的局限性,相關學者提出使用三維流場進行乘波體設計[8-9]。求解三維流場時,激波曲面的準確分辨是一個難題,Marcu使用激波捕捉法計算了圓錐流場,并沿流線根據壓力梯度確定激波位置,這一方法計算量較大,激波曲面的光滑性難以保證,推廣到復雜外形存在困難。激波裝配法[10]可以準確定位激波曲面,為尋找性能更優的流場進行乘波體設計提供了良好的工具。前期工作中筆者使用三維流線追蹤技術[11],探索了三維流場乘波體快速設計方法,在擴大設計空間方面具有優勢。

在三維流場乘波體設計中,隨著設計空間的擴大,合理選擇流場變得越來越重要。如果能找到乘波體與導波體外形間的對應關系,就可以先根據乘波體的要求尋找導波體,再用所得的導波體設計外形,從而提高乘波體設計的靈活性。本文使用激波裝配法計算帶激波的流場,固定流線追蹤(FCT,Flow Capture Tube)初始線為水平直線,投影到激波作為流線追蹤的起始點,即“平切”激波,正向追蹤流線生成乘波體外形。研究幾類導波體生成乘波體的算例,指出乘波體與導波體之間存在一定程度上的對應關系,可作為未來設計中導波體外形選擇的參考依據。

1 激波裝配法

帶激波的超聲速流場通常采用激波捕捉法數值求解。激波捕捉法無需將激波作為未知邊界進行求解,因此所得到的激波邊界實際上是一個流動量連續急劇變化的狹小區域,給激波的分辨帶來了困難,難以用于乘波體設計。

與激波捕捉法不同,激波裝配法則將流場的外激波作為未知外邊界同流場一起求解(圖1)。為此引入三維非定常激波的求解方程,包括Rankine-Hugoniot激波關系式及描述激波邊界處特征波傳播特性的特征相容關系式。將上述方程聯立數值求解,可求得波后流動量及激波速度。由于外激波邊界在求解過程中是移動的,因此求解過程需在動網格中進行。當方程解收斂時,激波速度趨近于零,非定常激波收斂為定常激波,即得到流場外激波的準確外形,所得激波面可直接用于乘波體設計。

圖1 激波捕捉法及裝配法比較Fig.1 Comparison between shock-capturing method and shock-fitting method

采用超聲速Euler方程激波裝配法數值計算程序求解導波體激波及流場。選擇文獻[12]所給物形作為測試算例,其中Ma=6,α=15°,圖2為計算得到的激波,圖3和圖4分別為激波位置rs和激波斜率rsx與文獻[12]數據的比較,可以看出兩者非常吻合。

圖2 激波外形計算結果Fig.2 Shock simulation results using shock-fitting method

圖3 激波位置計算結果比較Fig.3 Comparison of shock locations

圖4 激波斜率計算結果比較Fig.4 Comparison of shock wave slope

2 導波體

2.1 外形基本特征描述方法

在描述乘波體和導波體外形時采用如圖5所示的坐標系,其中Ux、Uy、Uz表示流場速度分量。

圖5 坐標系定義Fig.5 Definition of coordinate system

在此坐標系下,乘波體和導波體由三個坐標面上的外形描述:1)xy方向上描述縱向截面外形;2)yz方向對應橫截面外形;3)xz方向描述前緣外形。

固定流線追蹤初始線時,決定乘波體外形基本特征的因素包括:縱向截面外形與導波體流場速度Uy/Ux沿縱向(x向)分布特性相關;橫截面與Uy/Ux、Uz/Ux沿展向(z向)的分布特性相關;前緣外形及長寬比與激波形狀和流線追蹤初始線切割激波的位置等因素相關。因此乘波體外形的基本特征可以通過改變導波體的激波外形及速度分布得到。

下面給出兩類可改變激波外形及流場速度沿縱向及展向分布的導波體——錐形導波體和雙錐導波體,作為重點研究對象。

2.2 錐形導波體

錐形導波體的外形為尖頭,母線為直線,橫截面可為任意外形。圖6給出了幾種典型的橫截面形狀,其中截面為圓的外形即為錐導乘波體中的圓形尖錐導波體。

圖6 錐形導波體的橫截面外形Fig.6 Cross sections of conical SGBs

錐形流的流場特性與長度無關,即沿縱向分布的各橫截面具有相同的流場特性。改變橫截面外形,可以改變錐形導波體流場的速度Uy/Ux和Uz/Ux沿z向的分布特性,進而改變乘波體的橫截面外形。

2.3 雙錐導波體

雙錐導波體的外形特征是尖頭錐形前體、后接擴張或收縮的后錐體,橫截面為任意外形。過前體的流場為錐形流場,當后體收縮時,存在沿縱向傳播的膨脹波,后體擴張時,存在激波。

通過后體出現的膨脹波或激波,可以改變流場速度Uy/Ux、Uz/Ux沿縱向的分布特性,從而改變乘波體縱向截面外形;俯視外形前后錐體的長度、角度決定激波外形,進而決定乘波體的前緣外形;改變橫截面外形,可以改變流場速度沿展向的分布特性,從而改變乘波體橫截面外形。

3 乘波體與導波體外形對應關系

3.1 縱向截面外形對應關系分析

乘波體外形通過追蹤流線生成,縱向截面流場對于乘波體性能影響很大,目前常用的錐導乘波體等使用圓錐流動,流場分布較簡單,不存在激波膨脹波等現象,但也限制了設計空間。而增加波系則可以使外形具有更好的特性,比如呂偵軍等[13]提出了多級壓縮錐導乘波體的概念,使用多波系充分發揮前體的預壓縮作用,為進氣道的正常工作提供均勻流場。本文使用準確的三維流場進行設計,可以方便地構造流場內激波或膨脹波,提高了設計靈活性及精度。

3.1.1 流場內波系對速度分布的影響

導波體流場內存在激波或膨脹波時,由于波前波后速度的明顯變化,流線穿過這些波時走向會發生明顯改變。分析三個典型外形:尖錐、后體收縮的雙錐、后體擴張的雙錐,考察波系對速度的分布。

圖7給出了尖錐和后體收縮、擴張的雙錐導波體流場內Uy/Ux的分布特性。圖8給出了取定導波體流場內的同一位置截面ab,查看速度分布的特性??梢钥吹绞湛s后體產生膨脹波,波后Ux大于波前值,而Uy小于波前值,因此膨脹波波后Uy/Ux小于波前值;擴張雙錐內存在激波,激波波后速度分布與膨脹波相反,Uy/Ux大于波前值。

(a) 尖錐

(b) 后體收縮膨脹波

(c) 后體擴張激波圖7 雙錐導波體流場的Uy/Ux分布特性Fig.7 Uy/Ux distribution characteristics of double-conical flow

根據波前波后速度分布Uy/Ux的變化特性,流線穿過內膨脹波及內激波的趨勢可以大體確定,如果雙錐導波體后體收縮,所生成的乘波體邊界后段也向上收縮;雙錐導波體后體擴張,所生成的乘波體邊界后段也向下擴張。

圖8 雙錐導波體流場ab截面的Uy/Ux分布特性Fig.8 Uy/Ux distribution characteristics at ab cross section

3.1.2 具體算例分析

考慮三個導波體生成流場進行乘波體設計,設計狀態為Ma=6,攻角α=0°,其中外形1為8.5°半錐角的尖錐外形;外形2為后體收縮的雙錐外形(10°半錐角轉7°);外形3為后體擴張的雙錐外形(7°半錐角轉10°)。橫截面均為圓,圖9所示為三個導波體外形流場的Uy/Ux云圖。

(a) 錐形

(b) 收縮雙錐

(c) 擴張雙錐圖9 導波體外形的流場速度Uy/Ux云圖Fig.9 Velocity Uy/Ux contour of SGB flowfield

構造乘波體外形,下表面使用追蹤流線得到,上表面則考慮擴充容積設計[14],截面形狀用三次函數表示,形狀為“凸”。在不同乘波體的構造過程中保持長寬比不變,得到的乘波體外形為W1、W2和W3,對應于圖9中的導波體外形1、2和3,其縱向截面、橫截面及前緣外形的比較如圖10所示??梢娙齻€乘波體的前緣外形和橫截面外形均相差不大,主要區別是在縱向截面外形。

圖10 乘波體W1、W2和W3外形特征比較Fig.10 Configurations of W1, W2 and W3

使用CFD方法計算乘波體的流場,圖11給出了三個乘波體外形縱向截面和迎風面表面的壓力分布云圖。三個外形的壓力分布差別很大:W1的壓力分布比較均勻;W2由于膨脹波的影響后體表面壓力明顯低于前體;W3則由于激波的影響后體表面壓力明顯高于前體。這種不同的壓力分布特性對乘波體的縱向穩定特性有著顯著的影響。

圖11 乘波體壓力云圖對比Fig.11 Pressure contour of the waveriders

圖12給出了三個乘波體的升阻及縱向壓心特性,即升力CL·S(CL為升力系數,S為參考面積),升阻比L/D和縱向壓心系數Xcp的比較。由圖11給出的乘波體迎面壓力云圖可以看出,內激波的存在導致W3的迎風面高壓區較大,迎風面壓力高于W1,而膨脹波流場得到的外形迎風面壓力低于W1,因此W3的升力最大,W2升力最小。與升力特性相對應,在設計點(α=0°)附近,W3的L/D最大而W2最小,但是隨著攻角增大,W3的阻力迅速增加,導致其升阻比低于W1和W2。此外,由于W3迎風面高壓區集中在后部,導致壓心靠后,W2迎風面低壓區集中在后部,導致壓心靠前。

(a) CL·S~α

(b) L/D~α

(c) Xcp~α圖12 乘波體氣動性能比較Fig.12 Aerodynamic properties via angle-of-attack

上述分析可知,雙錐體在后體產生的激波或膨脹波能有效地改變導波體流場內速度沿縱向分布的特性,利用該特性可以設計不同縱向截面外形的乘波體,從而實現在俯視平面形狀受約束的情況下對乘波體構型縱向穩定性的按需設計,而這種乘波體性能設計的靈活性正是常規乘波體設計方法所難以提供的。

本節算例導波體采用了雙錐軸對稱外形,實際應用中只要求物形在縱向截面的迎風面存在收縮或擴張即可。

3.2 橫截面外形對應關系分析

乘波體的橫截面外形與導波體流場內速度Uy/Ux、Uz/Ux沿展向的分布特性有關,此速度特性可由導波體外形的橫截面確定。分析三個導波體外形:尖錐A0、雙橢圓錐D1、雙橢圓錐D2,如圖13所示,這三個導波體外形的縱向截面外形完全一致,橫截面不同——A0橫截面為圓,D1和D2橫截面為兩個橢圓的組合,D2的下部份外形(迎風面)比D1分布較為平坦。

(a) A0

(b) D1

(c) D2圖13 雙橢圓錐導波體與尖錐導波體Fig.13 Double elliptical cone SGB

圖14給出了導波體A0、D1和D2的橫截面比較,并給出了迎風面處沿cd位置流場速度Uy/Ux、Uz/Ux沿展向分布的比較。導波體在相同迎風面錐角條件下,橫截面迎風面越平,Uy/Ux的負值越大,而且在整個cd段都保持這一特性。Uz/Ux的分布特性有所不同,在cd前段,隨著導波體迎風面逐漸變平,Uz/Ux逐漸變小,而在cd后端,Uz/Ux快速增大。上述特性導致乘波體橫截面外形及氣動特性差別明顯。

在要求長寬比相同的條件下,由導波體A0、D1和D2生成的乘波體W-A0、W-D1和W-D2的外形如圖15所示。三個乘波體的縱向截面外形和前緣外形基本不變,區別主要在橫截面外形。由前面給出的速度沿展向分布特性可以得出,導波體迎風面越平,

(a) 三個導波體橫截面及激波形狀

(b) 速度Uy/Ux沿展向cd分布圖

(c) 速度Uz/Ux沿展向cd分布圖圖14 導波體流場的橫截面速度分布Fig.14 Velocity on cross section of different SGBs

圖15 乘波體三視圖比較Fig.15 Waverider configurations with different cross section profiles

生成的乘波體越平,內側厚度越大,而展向外緣則較薄。這種橫截面特征的對應關系可以在對乘波體外形的迎風面進行擴容設計時獲得應用。

3.3 前緣外形對應關系分析

乘波體前緣由導波體流場的激波曲面決定,由于錐形流導波體的激波外形為錐形或近似錐形,采用FCT為直線平切激波曲面時,其生成的乘波體前緣形狀一般為近似拋物線,如圖16所示。

而我們擴大設計空間,則可以得到前緣線具有特殊性質的乘波體外形,比如雙后掠外形等,這為我們改善乘波體低速階段的性能缺陷提供了一種新的思路。如果需要得到雙后掠外形的乘波體,則需采用俯視圖為雙前緣的導波體。如圖17所示的兩個導波體外形,Y1橫截面為雙橢圓,Y2為橢圓,縱向截面外形均為錐形,錐角為8°;俯視平面外形均為后體擴張的雙錐,錐角分別為7°和18°。

圖16 典型類錐導流場生成的乘波體外形Fig.16 Conical flow waverider

圖17 雙前緣導波體外形Fig.17 Double swept conical SGB

從圖18看到,導波體Y1、Y2的激波的俯視圖為擴張的雙前緣形狀,可生成前緣為雙后掠的乘波體外形,同時由于Y1、Y2的橫截面不同,生成的乘波體WY1的橫截面比WY2更平(圖19所示)。

圖20為乘波體WY1和WY2氣動特性的比較。與其他乘波體類似,升力呈線性增長,升阻比在α=4°左右取得最大,同時壓心靠后,并隨攻角增大壓心逐步前移。雖然WY1和WY2的縱向對稱面外形及前緣外形相差不大,但由于WY1橫截面下表面更平,橫截面積比WY2大,導致兩者氣動特性不同。

(a) 側視圖 (b) 俯視圖

(e) WY1 (f) WY2圖18 雙錐導波體激波流場及乘波體外形Fig.18 Shock wave from double conical SGBs

圖19 雙后掠乘波體外形Fig.19 Views of the double swept waverider

(a) CL·S~α

(b) L/D~α

(c) Xcp·S~α圖20 雙后掠乘波體氣動特性對比Fig.20 Aerodynamic performances of the double swept waverider

4 結 論

通過引入激波裝配法建立了普適于一般三維流場的乘波體設計方法,突破了導波體外形的限制。采用平切導波體激波及正向跟蹤流線生成乘波體,通過錐形流導波體和雙錐導波體生成乘波體的算例分析了導波體與乘波體在縱向截面外形、橫截面外形、前緣外形方面的對應關系。分析表明,雙錐體在后體產生的激波或膨脹波能有效地改變導波體速度沿縱向分布特性,利用該特性可實現在俯視平面形狀受約束的情況下對乘波體構型縱向穩定性的按需設計;通過改變導波體橫截面外形,可以改變速度沿展向的分布特性,從而實現乘波體橫截面外形改變;乘波體前緣外形主要與導波體激波在俯視面上的外形相關,通過后體擴張的雙錐外形能夠生成雙前緣乘波體外形。

[1]Nonweiler T R.Aerodynamic problem of manned space vehicles[J].Journal of the Royal Aeronautical Society, 1959, 63: 521-530.

[2]Rasmussen M L.Waverider configurations derived from inclined circular and elliptic cones[J].Journal of Spacecraft and Rockets, 1980, 17 (6) : 537-545.

[3]Jones J G, Moore K C, Pike J, et al.A method for designing lifting configurations for high supersonic speeds using axisymmetric flow field[J].Archive of Applied Mechanics,

1968, 37(1): 56-72.

[4]Center K, Sobieczky H, Dougherty F.Interactive design of hypersonic waverider geometries[R].AIAA 1991-1697, 1991.

[5]Sobieczky H, Dougherty F C, Jones K.Hypersonic waverider design from given shock wave[C]//The First International Waverider Symposium.Maryland: University of Maryland, 1990.

[6]Anderson J D.Fundamentals of aerodynamics[M].3rd edition.McGraw-Hill Companies, 2001.

[7]Mark J L, Maire L C.Shock-based waverider design with pressure gradient corrections and computational simulations[J].Journal of Aircraft, 2005, 42(5): 1350-1352.

[8]Marcu A.L, Kojiro S.Experimental investigation of a Mach 3.5 waverider designed using computational fluid dynamics[J].AIAA Journal, 2015, 53(6): 1590-1601.

[9]Le G G, Ma D W, Li Z Y.Computation of hypersonic flowfields for elliptic-cone-derived waverider[J].Journal of Nanjing University of Science and Technology (Natural Science), 2006, 30(3): 257-260.(in Chinese)樂貴高, 馬大為, 李自勇.橢圓錐乘波體高超聲速流場數值計算[J].南京理工大學學報(自然科學版), 2006, 30(3): 257-260.

[10]Moretti G.Thirty-six years of shock fitting method[J].Computers & Fluids, 2002, 31: 719-723.

[11]Liu C Z, Bai P, Chen BY, et al.Rapid design and optimization of waverider from 3D flow[J].Journal of Astronautics, 2016, 37(5): 535-543.(in Chinese)劉傳振, 白鵬, 陳冰雁, 等.三維流場乘波體快速設計方法及多目標優化[J].宇航學報, 2016, 37(5): 535-543.

[12]Любимов А Н, Русанов В В .Течения газа около тупых тел[M].НАУКА.МОСКВА, 1970[13]Lyu Z J, Wang J F, Wu Y Z, et al.Design and analysis of multistage compression cone-derived waverider configuration[J].Journal of Astronautics, 2015, 36(5): 518-523.(in Chinese)呂偵軍, 王江峰, 伍貽兆, 等.多級壓縮錐導乘波體設計與分析[J].宇航學報, 2015, 36(5): 518-523.

[14]Bowcutt K G.Optimization of hypersonic waveriders derived from cone flows including viscous effects[D].University of Maryland, 1986.

Waverider design and analysis based on shock-fitting method

CHEN Bingyan, LIU Chuanzhen*, JI Chuqun

(ChinaAcademyofAerospaceAerodynamics,Beijing100074,China)

Most commonly, a waverider is defined as an aerodynamic configuration that is inversely created from a prescribed hypersonic flow field based on a planar or conical shock wave.In this work, the prescribed shock wave was numerically created from an arbitrary configuration using shock-fitting method.A Shock Generating Body (SGB) was defined as the object to generate 3D flow field with shock wave for waverider creation.The correlation of the longitudinal section, cross section, and leading edge feature between the waverider configuration and its corresponding SGB configuration were analyzed.The analysis results show that the shock wave or expansion wave generated from the aft-body of a double-cone SGB can effectively alter the velocity distribution characteristics of the SGB flow field.This feature can be used to design the longitudinal characteristics of a waverider configuration when its planar profile is restricted.The velocity distribution on cross section can be altered by changing the cross section profile of the SGB, thereby achieving the cross section profile alteration of the generated waverider.A waverider with double-sweep leading edge feature can be generated using a double-cone with expanding aft-body.Qualitative correlations exist between the waverider and the SGB, which provides valuable insight as the guideline to choose appropriate flow field for the design of high performance waveriders.

shock-fitting method; three-dimensional flow field; waverider; aerodynamic shape; numerical simulation

0258-1825(2017)03-0421-08

2017-03-06;

2017-03-27

自然科學基金(11672281)

陳冰雁(1977-),男,廣東惠州人,研究員,研究方向:飛行器氣動設計.E-mail: chen_binyan@hotmail.com

劉傳振*(1989-),男,山東德州人,博士,研究方向:氣動外形設計.E-mail: chuanzhenliu@126.com

陳冰雁, 劉傳振, 紀楚群.基于激波裝配法的乘波體設計與分析[J].空氣動力學學報, 2017, 35(3): 421-428.

10.7638/kqdlxxb-2017.0039 CHEN B Y, LIU C Z, JI C Q.Waverider design and analysis based on shock-fitting method[J].Acta Aerodynamica Sinica, 2017, 35(3): 421-428.

V411.4

A doi: 10.7638/kqdlxxb-2017.0039

猜你喜歡
設計
二十四節氣在平面廣告設計中的應用
河北畫報(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
瞞天過?!律O計萌到家
藝術啟蒙(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
主站蜘蛛池模板: 热久久这里是精品6免费观看| 久久人人97超碰人人澡爱香蕉| 欧美日韩国产在线播放| 日本道综合一本久久久88| 国产人成网线在线播放va| 亚洲国产成人久久精品软件| 国产精品网址在线观看你懂的| 人妻21p大胆| 亚洲手机在线| 97久久超碰极品视觉盛宴| 在线视频亚洲欧美| 91亚洲影院| 国产91麻豆免费观看| 中国黄色一级视频| 日本影院一区| 秋霞一区二区三区| 国产亚洲欧美在线中文bt天堂| 精品一区二区三区水蜜桃| 少妇露出福利视频| 波多野结衣一二三| 看国产一级毛片| 国产精品浪潮Av| 亚洲天堂自拍| 色成人综合| 香蕉伊思人视频| 日韩高清一区 | 亚洲色图狠狠干| 日韩黄色大片免费看| 欧美一区国产| 久久精品aⅴ无码中文字幕| 亚洲国产成人精品一二区| 精品国产女同疯狂摩擦2| 国产精品无码制服丝袜| 婷五月综合| 天天视频在线91频| 亚洲欧美另类色图| 亚洲一区二区三区麻豆| 为你提供最新久久精品久久综合| 日韩在线成年视频人网站观看| 无码在线激情片| 久热精品免费| 久草视频精品| 青草视频免费在线观看| 国产精品一线天| 国产裸舞福利在线视频合集| 91免费观看视频| 亚洲一级无毛片无码在线免费视频| 激情亚洲天堂| 午夜无码一区二区三区| 91在线丝袜| 亚洲一区黄色| 黑人巨大精品欧美一区二区区| 五月六月伊人狠狠丁香网| 五月婷婷亚洲综合| 国产精品99r8在线观看| 人妻免费无码不卡视频| 精品无码一区二区三区电影| 72种姿势欧美久久久大黄蕉| 九九免费观看全部免费视频| 国产精品lululu在线观看| 日韩av电影一区二区三区四区| 亚洲综合在线网| 超薄丝袜足j国产在线视频| 国产在线91在线电影| 国产午夜无码专区喷水| 无码中文字幕精品推荐| 亚洲av日韩综合一区尤物| 成AV人片一区二区三区久久| 亚洲av无码片一区二区三区| 91午夜福利在线观看精品| 日韩毛片免费视频| 福利在线不卡一区| 精品久久久久成人码免费动漫| 亚洲精品在线观看91| 久草性视频| 日日噜噜夜夜狠狠视频| 玖玖精品视频在线观看| 成人在线视频一区| 91破解版在线亚洲| 亚洲欧美成人综合| 伊人无码视屏| 亚洲第一色视频|