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

HTV2第二次飛行試驗氣動熱環境及失效模式分析

2017-09-04 02:29:07國義軍張昊元代光月王安齡周述光中國空氣動力研究與發展中心四川綿陽621000
空氣動力學學報 2017年4期
關鍵詞:環境

國義軍, 曾 磊, 張昊元, 代光月, 王安齡, 邱 波, 周述光, 劉 驍(中國空氣動力研究與發展中心, 四川 綿陽 621000)

HTV2第二次飛行試驗氣動熱環境及失效模式分析

國義軍*, 曾 磊, 張昊元, 代光月, 王安齡, 邱 波, 周述光, 劉 驍
(中國空氣動力研究與發展中心, 四川 綿陽 621000)

采用數值模擬和工程計算相結合的方法對HTV2第二次飛行試驗的熱環境進行了復現,發現在40 km以下,翼前緣駐點線會發生邊界層轉捩,引起前緣熱流比層流情況增加55%,最大熱流達到11 MW/m2,燒蝕量約為3 mm,前緣高熱流導致法向應力超過碳布層與層之間的粘接強度,而縱向應力小于碳布拉伸破壞極限。因此本文認為,HTV2第二次飛行試驗失利的原因主要是:燒蝕疊加應力破壞,即在翼前緣由于燒蝕導致多層碳布被燒破,從而在翼前緣沿展向駐點線出現較長的破損口,而法向應力導致碳布層與層之間的粘接失去作用,在氣動力作用下,可能從燒破的地方開始將碳布掀起,嚴重影響氣動性能,最終導致飛行器無法控制。

HTV2;熱環境;燒蝕;熱應力;湍流加熱

0 引 言

HTV-2(見圖1)是美國空軍和國防部預先研究計劃局(DARPA)聯合開展的Falcon計劃[1]的一部分,由洛克希德·馬丁公司制造,主要用于驗證高超聲速滑翔機動飛行器的氣動布局、氣動熱防護設計、材料、控制等關鍵技術,目標是確保美國在近遠期具備全球快速精確打擊能力。2010年4月22日進行了首次飛行試驗,結果以失敗告終[2],主要原因是控制要求超出了飛行器的能力范圍。

(a)

(b)

2011年8月11日,美國DARPA進行了HTV2第二次飛行試驗,仍然以失敗告終。調查飛行異常原因的美國獨立工程審查委員會(ERB)經過為期7個月的大量分析和額外地面測試,公布的第二次飛行失利的調查報告認為:作為實現在不到一小時內抵達全球任何地方實現快速打擊的能力的技術演示驗證和數據采集平臺,HTV2第二次試驗飛行證實了飛行器的氣動設計是有效的。此次飛行成功驗證了以高達Ma20的速度實現接近3 min的穩定氣動控制飛行。期間飛行器經歷了超出設計能夠承受的100倍的最初激波擾動,而飛行器能夠恢復并繼續可控飛行。在試驗飛行9 min時,飛行器異常地經歷了一系列強烈震動,自主飛行安全系統試圖利用飛行器氣動系統實現可控的降落并濺落于海洋。ERB總結到:“HTV2第二次飛行提前終止飛行最可能的原因是沒有預料到氣動殼體退化,產生了多個增加嚴重性的意外,最終激活了飛行安全系統”。 調查報告中提到,基于先進模型、高溫材料地面試驗,以及對其他已熟知的飛行機制的熱效應的認識,預計飛行器蒙皮在達到應力容忍限度時會產生一定的梯度性剝落。然而,飛行器蒙皮從氣動結構上剝落的部分遠遠大于預期的程度。當飛行器以每小時13000英里速度飛行時,因此而產生的縫隙將在飛行器周圍產生強烈的脈沖激波,因而導致飛行器突然滾轉。根據首次飛行試驗獲得并集成進第二次飛行的認識,飛行器的氣動穩定性使得其能夠在幾次激波導致的滾轉后成功糾正自身。盡管如此,連續擾動的嚴重程度最終超出了飛行器自我恢復的能力。HTV2第二次飛行試驗中采集的數據揭示了對熱防護材料特性的新認識,以及在大氣層內Ma20速度飛行的不確定性。第二次飛行的數據顯示,從已知的飛行機制的推斷和僅依靠先進的熱建模和地面測試是無法成功預測20馬赫大氣飛行下嚴酷現實的。

從以上公開的美國調查結論可以看出,第二次飛行失利的原因可能是氣動熱和防熱方面出現了問題,最終導致飛行器出現強烈振動從而無法控制,但并沒有給出具體細節的描述。

本文根據HTV2外形、彈道和防熱結構,就HTV2熱環境、燒蝕、溫度場和熱應力進行了深入計算分析,初步推測出導致飛行失敗的原因。

1 全機熱環境計算

為了分析確定防熱系統到底是哪個方面出現了問題,首先需要把熱環境搞準。HTV2飛行器布局看似簡單, 但受熱特征卻較為復雜,如果采用的方法不當,給出的熱環境數據會有很大差異。初步計算分析表明,其熱環境主要特點有:

1) 駐點屬于三維駐點,端頭縱向和橫向曲率半徑不同,俯仰平面R=17.6 mm,水平面R=24.3 mm; 2) 翼前緣后掠角很大(73°),前緣駐點線可能發生邊界層轉捩,駐點線上的湍流熱流顯著高于層流熱流;

3) 頭激波和前緣激波交匯對前緣和迎風面會形成一定的干擾,會引起水平翼前緣熱流局部增加,在迎風面會形成“條狀”干擾熱流。除此之外,激波干擾還會引起翼前緣邊界層提前轉捩,從而使干擾點后的整個翼前緣處于湍流狀態,導致熱環境大幅升高。

1.1 熱環境計算方法

為了搞準HTV2的熱環境,本文采用CARDC三套數值計算程序和兩套工程熱環境計算程序進行了對比計算。

數值計算[3-5]采用完全和非平衡兩種氣體模型,基于有限體積法,考慮層流和湍流兩種情況,采用兩套網格125萬和1000萬,35公里時網格雷諾數分別為32和16。

工程計算采用等價錐法和流線法兩套程序[6-8],也進行了層流和湍流、完全氣體和平衡氣體對比計算。

根據比對結果,最終選用經過數值計算校核的工程模型沿彈道計算全機熱環境??紤]到HTV2為扁平體外形,針對不同的展向截面,采用二維片條法計算每一截面的熱環境,之后再將所有片條組合起來,插值獲得全機熱環境。

計算分析表明,翼前緣的熱環境準確與否,對整個分析至關重要。對層流情況,采用下式計算后掠前緣駐點線熱流[9]:

對湍流情況,本文采用下式計算后掠圓柱前緣加熱[11]:

(2)

其中qcyl和qcyt分別為鈍頭體前緣駐點線上層流和湍流加熱率,kT/kL為湍流和層流加熱放大因子之比,rc為鈍頭體前緣曲率半徑。

1.2 熱環境計算結果

1.2.1 駐點熱流

圖2給出了HTV2彈道特性和駐點熱流計算結果,從頭體分離開始算起,第一次拉起最低點時刻為105.67 s,飛行高度為35.109 km,迎角 7.8°,馬赫數18.68。拉起最低點形成駐點熱流和壓力的峰值,工程計算給出的峰值熱流密度為24.2 MW/m2,峰值壓力為269 kPa。

圖2 駐點熱流沿彈道變化情況Fig.2 Heating rate at stagnation point along trajectory

1.2.2 全機熱環境

1) 數值計算結果

圖3給出了全層流和全湍流情況下典型時刻飛行器熱流分布計算結果。需要說明的是,實際上機身頭部x=0駐點附近的區域不可能是湍流狀態,因此駐點附近的湍流計算結果是不真實的,但除了x=0及其下游很小的區域外,其它地方都有可能出現湍流狀態,包括機翼前緣駐點線。從圖中可以看出,機身表面湍流熱流顯著高于層流熱流,特別是翼前緣,層流情況熱流為5 MW/m2左右,而湍流情況熱流高達8.2 MW/m2。

(a) 前緣駐點線熱流

(b) Z=500翼剖面熱流

(a) 機身下表面熱流

(b) 激波結構

(c) 激波干擾產生的高熱流區

數值計算表明,HTV2熱環境有一個顯著特點就是頭激波和機翼前緣激波交匯對前緣和迎風面會形成一定的干擾,會引起水平翼前緣熱流局部增加,在迎風面會形成“條狀”干擾熱流(見圖4),激波交匯干擾會引起翼前緣(Z=100~200)熱流局部增加20~30%左右。

表1給出了不同方法計算的駐點熱流比較,這里同時給出了數值方法和工程方法計算的熱流,可以看出,數值方法的結果與工程計算結果吻合很好,考慮到工程方法采用平衡氣體模型,其結果略高于完全氣體數值結果是合理的。

表1 不同方法計算的駐點熱流比較Table 1 Comparison of heating rate at stagnation point predicted by numerical to engineering methods

2) 工程計算結果

本文最終采用工程方法[6-12]的計算結果,這里考慮了翼前緣轉捩和湍流熱流,并考慮了頭激波與翼激波交匯干擾引起的熱增量。圖5給出了典型時刻計算得到的不同展向截面熱流分布,機身頭部駐點熱流為24 MW/m2左右,機翼前緣駐點線在頭激波與翼激波交匯前為層流狀態,最大熱流為6 MW/m2左右,交匯點最大熱流為11.2 MW/m2左右,交匯點后整個翼前緣都處于湍流狀態,最大熱流為9.4 MW/m2左右。為了考察邊界層轉捩的發展情況,圖6給出了對稱面不同時刻熱流分布,大概從45km起開始從尾部出現邊界層轉捩,到最低拉起點35.1 km時,身部邊界層轉捩起始點已經移到x=0.683 m處,考慮到激波交匯干擾會對邊界層轉捩產生影響,對于翼前緣,可以認為邊界層轉捩提前到z=200 mm處。工程計算結果與數值計算結果的對比情況見表1,考慮到翼前緣的湍流狀態和真實氣體效應,本文認為翼前緣大部分區域熱流應取9.4 MW/m2左右,比層流的6 MW/m2情況高出55%。

圖5 工程計算得到的不同展向截面熱流 (H=35 km,M=18.68,迎角7.8°)Fig.5 X direction distribution of predicted heating rate at Z cross sections parallel to pitching plane

圖6 對稱面不同時刻熱流分布Fig.6 Predicted heating rate distribution on body centerline at different times

2 典型部位的燒蝕情況

本文選擇了彈體上一些特征點進行了燒蝕防熱計算,表2給出了HTV2防熱材料和結構。

表2 HTV2防熱材料和結構Table 2 HTV2 TPS material and structure

圖7給出了展向Z=0 mm剖面特征點熱流、壁溫和燒蝕量沿彈道隨時間變化情況計算結果。為了便于分析和應用,熱流計算同時給出了冷壁熱流和熱壁熱流的計算結果,可以看出,沿再入彈道,受壁溫的影響,冷壁熱流與熱壁熱流的差別越來越大,后者遠遠小于前者。駐點最高外表面溫度達到3700 K,截止到161 s,端頭燒蝕量約為12.3 mm,身部大面積最高壁溫為2200 K左右,燒蝕量小于0.25 mm。

圖8給出了展向不同位置翼前緣駐點線熱流和燒蝕量沿彈道變化情況。Z=100 mm截面對應層流加熱,截止到161 s,燒蝕量為2 mm;Z=200 mm截面對應激波交匯點加熱,燒蝕量達到3.28 mm;Z=500 mm截面對應湍流加熱,燒蝕量為2.57 mm??紤]到每層碳布的厚度不到1 mm,所以沿翼前緣有2~3層碳布被燒破。

(a) 駐點熱流

(b) 其它特征點冷壁熱流

(c) 駐點內外壁溫

(d) 其它特征點壁溫

(e) 駐點燒蝕量

(f) 其它特征點燒蝕量

(a) 翼前緣駐點線熱流

(b) 翼前緣駐點線燒蝕量

3 熱應力計算結果

基于自研的三維溫度場和熱應力計算軟件[13-14],選取翼面(迎風+背風)一個條帶(Z=50 mm至Z=400 mm之間)作為計算分析對象,該區域包含了頭激波與翼前緣激波交匯的位置(Z≈200 mm),并認為在Z>200 mm之后加載的熱環境為湍流氣動熱。

HTV2殼體為二維碳布包裹結構,材料本身屬各向異性導熱材料,表 3給出了材料的有關物性參數。由于材料的具體編織和纏繞方式未知,這里選取K=4 W/m·K和K=45 W/m·K兩個導熱系數分析溫升歷程,實際情況應介于這兩種極限條件之間。

表3 2D純碳/碳熱物理及力學性能參數Table 3 Thermal physical and mechanical performance of 2D C/C material

圖9給出了兩種材料特性下翼前緣溫度變化歷程,最大溫升出現在105 s左右,超過3000 K。圖10給出了t=105 s時刻表面溫度分布,高溫區主要集中在前緣附近,前緣和翼面溫度差異巨大。

圖9 前緣點不同導熱系數下的溫升歷程對比Fig.9 Leading-edge temperature time histories

溫度變化劇烈的區域通常熱應力也較大,本文采用50 mm×50 mm×25 mm的平板模型模擬翼前緣各向異性材料的應力情況??紤]到不同邊界條件對結果的影響,這里分別使用X和Y方向約束和全部無約束兩種邊界條件(圖11、圖12)。

(a) 導熱系數K=4 W/m·K

(b) 導熱系數K=45 W/m·K

(a) 105 s等效應力云圖 (b) 105 s Z方向位移量云圖

1) 各向異性X、Y方向約束,Z方向無約束:計算結果表明,105 s等效應力最大值538 MPa。

2) 各向異性無約束:105 s等效應力最大值52 MPa。

(a) 105 s等效應力云圖 (b) 105 s Z方向位移量云圖

真實情況應該介于無約束和全約束之間,即最大應力應在52 MPa至538 MPa之間,都超過了層與層之間材料的連接強度(見表3),說明翼前緣附近碳布層與層之間粘接都已失效。

考慮到材料的燒蝕情況,飛行器前緣會因燒蝕出現沿展向的破裂縫。由于熱應力和氣動力共同作用,2D碳布防熱材料可能會有3層被掀起,從而導致飛行器外形發生較大變化,并最終導致飛行器失穩。

4 結 論

本文根據HTV2外形和第二次飛行彈道,采用數值模擬和工程計算相結合的方法確認熱環境,發現在40 km以下,翼前緣駐點線可能發生邊界層轉捩,引起前緣熱流相對于層流情況增加55%左右,由此引起燒蝕量顯著增大。計算結果表明:

頭部駐點最大熱流為24 MW/m2左右,燒蝕量為12.3 mm;

翼前緣50~200 mm為層流加熱,最大熱流為6 MW/m2左右,燒蝕量小于2 mm;

翼前緣Z=200 mm處由于激波交匯干擾,最大熱流為11.2 MW/m2,燒蝕量高達3.28 mm;

翼前緣Z≥200 mm在40 km以下出現湍流加熱,最大熱流為9.4 MW/m2左右,燒蝕量為2.57 mm;

身部大面積區域最大熱流為3 MW/m2左右,燒蝕量很小。

通過熱應力計算發現,飛行器前緣區域由于高熱流和大熱流梯度導致法向應力超過碳布層與層之間的粘接強度,使得粘接層失效,而縱向應力小于碳布拉伸破壞極限。

因此本文認為HTV2第二次飛行試驗失利的原因主要是:燒蝕疊加應力破壞。HTV2是由多層2D碳布包裹而成的,每層碳布厚度不足1mm,碳布的層與層之間采用粘接方式。計算表明,翼前緣燒蝕量達到2~3.3mm,導致2~3層碳布被燒破,從而在翼前緣沿展向駐點線出現較長的破損縫,而法向向外的拉應力導致碳布層與層之間的粘接失去作用,在氣動力作用下,可能從燒破的地方開始將碳布掀起,嚴重影響氣動性能,并最終導致飛行器無法控制。

本文認為不太可能是純應力拉伸破壞。盡管翼前緣存在高溫和大熱流梯度,可能會產生拉伸應力,但由于膨脹幅度不大,而碳布拉伸強度很高,不太可能會出現拉伸破壞。

本文認為也不太可能是擠壓破壞。翼前緣高溫膨脹會使碳布從前緣向飛行器中部擠壓,可能會使碳布隆起,但由于碳布層與層之間已經剝離,不太可能會將碳布折斷,而且碳布采用二維編織結構,在高溫情況下會發生結構變形使應力松弛掉。

[1]Walker S H, Sherk J. The DARPA/AF Falcon Program: The Hypersonic Technology Vehicle #2 (HTV-2) flight demonstration phase[R]. AIAA 2008-2539, 2008

[2]Li Jianlin. Research on development of hypersonic near space vehicle[M]. Beijing: China Astronautic Publishing House, 2012. (in Chinese)李建林. 臨近空間高超聲速飛行器發展研究[M]. 北京: 中國宇航出版社, 2012

[3]Zhang Haoyuan, Zong Wengang, Gui Yewei. Numerical investigation of flow in leading-edge gap of hypersonic vehicle[J]. Chinese Journal of Astronautics, 2014, 35(8): 893-900. (in Chinese)張昊元, 宗文剛, 桂業偉. 高超聲速飛行器前緣縫隙流動的數值模擬研究[J]. 宇航學報, 2014, 35(8): 893-900

[4]Li Zuowu. Study on the dissipative effect of approximate riemann solve on hypersonic heatflux simulation[J]. Chinese Journal of Theoretical and Applied Mechanics, 2008, 40(1): 19-25. 黎作武. 近似黎曼解對高超聲速氣動熱計算的影響研究[J]. 力學學報, 2008, 40(1): 19-25

[5]Dong Weizhong, Ding Mingsong, Gao Tiesuo, et al. The influence of thermo-chemical non-equilibrium model and surface temperature on heat transfer rate[J]. Acta Aerodynamica Sinica, 2013, 31(06): 692-698. (in Chinese)董維中, 丁明松, 高鐵鎖, 等. 熱化學非平衡模型和表面溫度對氣動熱計算影響分析[J]. 空氣動力學學報, 2013, 31(06): 692-698;

[6]Thermal Environment and Ablation/Erosion Analysis Software[AEROHEATS, V1.0]. China Aerodynamics Research and Development Center, Computer Software Copyright Registration Certificate(Registration Mark: 2013SR132872, Certificate No. 0638634), 2013高超聲速飛行器熱環境及燒蝕/侵蝕綜合分析軟件系統[簡稱AEROHEATS]V1.0版. 中國空氣動力研究與發展中心計算空氣動力研究所, 中華人民共和國計算機軟件著作權登記證書(登記號: 2013SR132872, 證書號: 0638634號), 2013

[7]Wang Anling, Guiyewei, Tang Wei, et al. Physical model establishment for the thermal torridor of a reusable launch vehicle[J]. Journal of Engineering Thermophysics, 2006, 27(5): 856-858. (in Chinese)王安齡, 桂業偉, 唐偉, 等. 可重復使用飛行器熱走廊物理建模研究[J]. 工程熱物理學報, 2006, 27(5): 856-858

[8]DeJarnette F R. Calculation of inviscid surface streamlines and heat transfer on shuttle type configurations[R]. NASA CR-111921, 1971

[9]Zhang Zhicheng, Pan Meilin, Liu Chuping. Hypersonic aerothermodynamics and thermal protection[M]. Beijing: National Defence Industry Press, 2003. (in Chinese)張志成, 潘梅林, 劉初平. 高超聲速氣動熱和熱防護[M]. 北京: 國防工業出版社, 2003

[10]Poll D. The effect of wing sweep back upon transtion in hypersonic flow[R]. AIAA-95-6090, 1995

[11]Engel C D. Miniver upgrade for the avid system, volume 1: lanmin user’s manual[R]. NASA CR-172212, 1983

[12]Guo Yijun, Dai Guangyue, Gui Yewei, et al. Engineering calculation of non-equilibrium effects on thermal environment of reentry vehicles[J]. Acta Aerodynamica Sinica, 2015, 33(5): 581-587. (in Chinese)國義軍, 桂業偉, 童福林, 等. 再入飛行器非平衡氣動加熱工程計算方法研究[J]. 空氣動力學學報, 2015, 33(5): 581-587

[13]Huang Qian, Gui Yewei, Geng Xiangren. Numerical simulation of thermal-stress in multi-layer plate[J]. Journal of Engineering Thermophysics, 2005, 26(5): 862-864. (in Chinese)黃謙, 桂業偉, 耿湘人. 層狀平板內熱應力的計算研究[J]. 工程熱物理學報, 2005, 26(5): 862-864

[14]Geng Xiangren, Gui Yewei, He Lixin, et al. Numerical study on heat transfer and thermal stress for infra-window with externally cooled and internally colled techniques[J]. Acta Aerodynamica Sinica, 2008, 26(3): 329-333. 耿湘人, 桂業偉, 賀立新, 張來平, 紅外窗口不同冷卻方式下的結構傳熱和熱應力特性計算研究[J]. 空氣動力學學報, 2008, 26(3): 329-333.

Investigation on aerothermodynamic environment and ablation which lead to HTV-2 second fight test failing

GUO Yijun*, ZENG Lei, ZHANG Haoyuan, DAI Guangyue, WANG Anling, QIU Bo, ZHOU Shuguang, LIU Xiao
(China Aerodynamics Research and Development Center, Mianyang 621000, China)

The second flight of the hypersonic technology vehicle 2 (HTV-2) failed on August 11, 2011. According to the engineering review board (ERB) analysis report, the most probable cause of the premature flight termination is unexpected aeroshell degradation, creating multiple upsets of increasing severity that ultimately activated the flight safety system. In order to investigate this failure, the aerothermodynamic environment, ablation, and stress are calculated and analyzed in this paper using numerical simulation and empirical calculation combined method, based on reconstructed HTV-2 configuration and the flight trajectory. It has been found that, at the height of 40 km, there is a possibility of boundary layer transition, leading to turbulent flow along the leading edge of the wing. Especially, due to the shock interaction, the transition moves forward, and the heating rates are 55% higher than those with laminar flow condition atZ=200 mm in spanwise direction at leading edge. The maximal value of cold wall heating rate reaches 11 MW/m2, and the ablation recession is approximately 3mm. Since the thickness of 2-D carbon cloth is only 1mm, there are two to three layers of carbon cloth are burned up at body leading edge. At the same time, the normal stress exceeds the bonding strength between carbon cloth layers. The most probable cause for the termination of the HTV2 second premature flight can be concluded according to the present study. The aeroshell degradation process can be described as follows. An unexpected great ablation at the leading edge breaks several layers of the carbon cloth, resulting in a long breakage, meanwhile the normal stress invalidates the bonding between carbon cloth layers. Under the influence of aerodynamic force, several layers of the carbon cloth can be lifted up from the breakage. This behaviour has a strong impact on the stable aerodynamically controlled flying, and finally activates the vehicles autonomous flight safety system to make a controlled descent and splashdown in the ocean.

HTV2; aerothermodynamic environment; ablation; thermal stress; turbulent heating

0258-1825(2017)04-0496-08

2016-10-08;

2016-12-20

國家重點基礎研究發展計劃資助(2014CB744100);國家重點研發計劃“大科學裝置前沿研究”重點專項資助(基金號2016YFA0401200)

國義軍*(1966-),男,博士,研究員,博導,主要從事高超聲速氣動熱和防熱研究.E-mail:13778169233@163.com

國義軍, 曾磊, 張昊元, 等. HTV2第二次飛行試驗氣動熱環境及失效模式分析[J]. 空氣動力學學報, 2017, 35(4): 496-503.

10.7638/kqdlxxb-2016.0114 GUO Y J, ZENG L, ZHANG H Y, et al. Investigation on aerothermodynamic environment and ablation which lead to HTV-2 second fight test failing[J]. Acta Aerodynamica Sinica, 2017, 35(4): 496-503.

V211.3

A doi: 10.7638/kqdlxxb-2016.0114

猜你喜歡
環境
長期鍛煉創造體內抑癌環境
一種用于自主學習的虛擬仿真環境
孕期遠離容易致畸的環境
不能改變環境,那就改變心境
環境與保護
環境
孕期遠離容易致畸的環境
高等院校環境類公選課的實踐和探討
掌握“三個三” 兜底環境信訪百分百
我國環境會計初探
中國商論(2016年33期)2016-03-01 01:59:38
主站蜘蛛池模板: 久久久精品无码一二三区| 欧美日韩资源| 国产在线小视频| 天天色综网| 久久九九热视频| 久久综合成人| 日本黄色不卡视频| 伊在人亞洲香蕉精品區| 又粗又硬又大又爽免费视频播放| 国外欧美一区另类中文字幕| 久久精品aⅴ无码中文字幕| 91麻豆精品国产91久久久久| 国产女人在线| 亚洲日本精品一区二区| 2021国产精品自产拍在线观看| 欧美久久网| 毛片在线播放a| 亚洲综合一区国产精品| 高h视频在线| 超薄丝袜足j国产在线视频| 一本综合久久| 欧美精品一二三区| 一级不卡毛片| 亚洲日韩精品伊甸| 97久久精品人人做人人爽| 午夜国产理论| 福利视频一区| 高潮爽到爆的喷水女主播视频| 澳门av无码| 亚洲精品中文字幕无乱码| 国产精品久久久久无码网站| 亚洲中文字幕国产av| 亚洲综合专区| 成人夜夜嗨| 啪啪永久免费av| 国产午夜一级淫片| 亚洲69视频| 免费一级毛片不卡在线播放| 精品国产免费人成在线观看| 亚洲中文字幕久久精品无码一区| 美女被操黄色视频网站| 日韩无码真实干出血视频| 国产成人资源| 日日噜噜夜夜狠狠视频| 国产美女无遮挡免费视频网站| 亚洲视频欧美不卡| 一级看片免费视频| 男女精品视频| 国产在线自乱拍播放| 色哟哟国产精品一区二区| 欧美日韩专区| 午夜爽爽视频| 高清国产va日韩亚洲免费午夜电影| 国产精品女人呻吟在线观看| 国产精品制服| 日本欧美视频在线观看| 福利小视频在线播放| 丁香六月激情综合| 久久久久人妻一区精品| 成人在线欧美| 一本色道久久88综合日韩精品| 一级全黄毛片| 免费不卡在线观看av| 国产视频入口| 91年精品国产福利线观看久久 | 日韩欧美中文字幕在线精品| 毛片三级在线观看| 免费三A级毛片视频| 亚洲综合香蕉| 亚洲一道AV无码午夜福利| 国产黄色爱视频| a级毛片免费看| 2020亚洲精品无码| 亚洲欧美在线综合图区| 成年人国产网站| 国产农村妇女精品一二区| 19国产精品麻豆免费观看| 欧美一级高清片欧美国产欧美| 日韩欧美色综合| 美女无遮挡免费网站| 亚洲精品自拍区在线观看| 国产精品自在拍首页视频8|