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

發(fā)動機噴流對飛機阻力傘性能的影響

2021-05-06 03:06:30馮傳奇孫建紅喻東明許常悅寇澤普
南京航空航天大學學報 2021年2期
關鍵詞:發(fā)動機

孫 智,馮傳奇,孫建紅,3,喻東明,許常悅,寇澤普

(1.南京航空航天大學航空學院飛行器環(huán)境控制與生命保障工信部重點實驗室,南京210016;2.中國飛行試驗研究院,西安710089;3.南京航空航天大學民航學院,南京211106;4.航空工業(yè)航宇救生裝備有限公司航空防護救生技術航空科技重點實驗室,襄陽441003)

飛機阻力傘是一種由柔性織物制造的氣動減速設備,通常折疊放置在飛機尾部的傘倉內,在飛機著陸時被拋出并打開從而輔助機輪制動,阻力傘產(chǎn)生的氣動阻力能夠有效縮短飛機滑跑距離,增加飛機減速制動的安全性。研究表明,阻力傘可使飛機著陸滑跑距離縮短30%~40%,有效延長了機輪外胎的使用壽命,保障了飛機的著陸安全[1]。然而,隨著航空部隊訓練的需求以及航空技術的發(fā)展,對戰(zhàn)斗機的著陸安全提出了更高的要求,在戰(zhàn)斗機中斷起飛以及著陸過程中,飛機發(fā)動機往往會處于未完全停機狀態(tài),此時發(fā)動機噴流產(chǎn)生的高速射流會對阻力傘工作流場產(chǎn)生較大的影響,進而影響阻力傘的氣動性能。因此,亟需對發(fā)動機噴流作用下的阻力傘工作過程進行研究,探究發(fā)動機噴流影響下的阻力傘氣動性能,為噴流作用下阻力傘的性能評估提供理論基礎。

對于單獨阻力傘的研究起步較早。在20 世紀70 年代,國外就己經(jīng)開始對降落傘工作過程進行數(shù)值模擬研究。降落傘的工作過程涉及非常復雜的氣動彈性問題,人們從早期對降落傘充氣過程的機理研究,逐步發(fā)展到對阻力傘相關流-構耦合的研究。例如,郭亮等[2]針對無人機傘回收進行了動力學分析,對無人機的整個回收過程進行仿真研究。程涵等[3]采用流-構耦合方法對某型傘低速氣流下的充氣過程進行了數(shù)值模擬,并對充氣過程中傘衣應力、流場速度矢量、壓力以及傘衣半徑變化等結果進行了分析。陳猛等[4]對某五環(huán)錐傘進行了無限質量情況降落傘充氣過程的流-構耦合數(shù)值模擬并進行了相關飛機減速過程開傘試驗。包進進等[5]對降落傘傘包拉出過程進行了仿真研究,分析了拉出過程的載荷變化。孫建紅等[6-7]采用質量阻尼彈簧模型分析了阻力傘拉直過程的引導傘阻力面積、傘系統(tǒng)彈性模量以及線密度等影響因素,也采用流固耦合方法對重裝空投降落傘充氣過程進行了仿真分析。Kimata 等[8]采用浸沒邊界法(Immersed boundary method,IBM)結合質量彈簧阻尼對半球殼形傘衣進行了馬赫數(shù)為2 的二維和三維流固耦合數(shù)值模擬,指出二維模擬的阻力系數(shù)振蕩幅度與試驗相比較小,三維模擬的流場中包含分離流和重建流兩種不同的流動模式。Borke 等[9]采用嵌入邊界法(Embedded boundary method,EBM)對二維傘衣進行馬赫數(shù)為2 的流固耦合研究,其中流場采用基于歐拉頂點的有限體積方法進行模擬,并將流場網(wǎng)格用自適應網(wǎng)格細化方法進行實時加密,結果表明在傘衣變形過程中,傘衣前方的激波能夠被準確地捕捉到。薛曉鵬等[10-11]采用IBM 研究了傘繩對拖曳比為2.38 的超聲速降落傘系統(tǒng)的影響,并采用這種方法對半球殼形傘衣進行了馬赫數(shù)為2 的流固耦合數(shù)值模擬,研究了傘衣攻角和前體攻角對降落傘性能的影響。Huang等[12]采用EBM 對2 馬赫的二維降落傘進行了充氣過程模擬,傘繩采用主從運動學進行模擬,發(fā)現(xiàn)初始的折疊狀態(tài)會顯著改變傘衣阻力系數(shù)和最大主應力。高興龍等[13]采用多材料任意拉格朗日-歐拉方法研究了前體的存在對超聲速盤縫帶傘充氣過程流場的影響,得到了前體尾渦結構和傘衣周圍的激波分布。蔡志軍等[14]采用流固耦合方法對戰(zhàn)斗機阻力傘載荷進行了仿真計算。楊品等[15]采用流固耦合方法對飛機阻力傘放傘過程進行了數(shù)值仿真,并對阻力傘拉直力及開傘過載等進行了分析。

上述研究主要針對阻力傘,而并未考慮發(fā)動機尾流的影響。對于發(fā)動機尾流場的研究主要關注于發(fā)動機本身以及噴流高溫尾流影響區(qū)域。例如,姚金華等[16]研究了固體火箭助推器產(chǎn)生的高溫燃氣尾流場對地面發(fā)射架造成的沖刷影響。李茂等[17]針對氫氧火箭發(fā)動機地面水平試車時尾流燃氣對地面熱防護的影響進行了研究,并采用數(shù)值模擬方法研究了尾流場溫度變化的趨勢。在飛機發(fā)動機尾流場計算方面,目前的研究主要側重于溫度場模擬、測量與對流場內設備的影響及防護。如于芳芳等[18]對某型發(fā)動機尾流流場及溫度場進行了數(shù)值計算,分析了不同發(fā)動機狀態(tài)下的溫度場特性,并研制了一套用于測量高溫條件下的尾流測溫系統(tǒng)。吳沿慶等[19]采用CFD 方法對飛機發(fā)動機尾流流場進行了數(shù)值模擬,采用反向蒙特卡洛方法對氣體紅外輻射特性進行了計算。黃爍橋等[20]采用試樣的方法,研究了噴流對飛機尾流渦的影響。

綜上可見,盡管國內外已有不少關于阻力傘和發(fā)動機噴流的研究,但鮮有對兩者相互影響的相關研究,特別是針對發(fā)動機噴流作用下阻力傘的流固耦合研究。因此,本文針對不同發(fā)動機噴流速度對阻力傘的氣動性能影響進行了數(shù)值仿真研究,分析了噴流速度對阻力傘阻力特性、穩(wěn)定性以及流場特性的影響,為發(fā)動機噴流作用下的阻力傘性能評估和阻力傘設計提供參考。

1 數(shù)值方法

1.1 控制方程

針對發(fā)動機噴流工況下的飛機阻力傘開傘過程,本文采用了時-空守恒元解元方法(CE/SE)方法對其進行求解。其N-S 方程可以寫成如下形式

式中:ρs為織物材料密度,xsi為結構單元位移。

1.2 耦合方法

文中飛機阻力傘織物間的接觸以及織物結構與流場間的耦合作用基于罰函數(shù)來實現(xiàn),采用具有二階時間精度的中心差分格式對流場與結構體進行顯示耦合計算,并對流體節(jié)點和結構體節(jié)點的速度以及位移進行求解,具體公式如下

式中:u 為速度矢量;M 為質量對角矩陣;Fext為外力矢量,F(xiàn)int為內力矢量,它們與體力和邊界條件相關聯(lián);s 為位移矢量。

本文在流體/結構體耦合計算中使用準約束方法,對于流體部件,解算器基于歐拉方法進行求解,而對于結構部件,則基于拉格朗日方法求解。界面邊界位置和速度由拉格朗日結構決定,并將該信息作為流體求解時每個時間步的界面條件,通過追蹤結構和流體之間的相對位移,計算出界面力,并作為外部力的一部分對耦合區(qū)域的速度、位移進行迭代,從而實現(xiàn)耦合計算。

2 模型與網(wǎng)格

為了對發(fā)動機噴流尾流場影響下的阻力傘阻力性能進行研究,文中選取某型飛機阻力傘進行研究[14-15],阻力傘模型如圖1 所示。圖中,縱橫交錯的線條為傘衣上的加強帶,傘衣臂上分布有狹小的開縫結構,縫兩側使用連接繩相連接,傘繩與加強帶在傘衣底邊相連,4 根加強帶在傘衣頂孔處匯聚為一點。根據(jù)十字傘實際折疊過程,得到簡化的十字形傘折疊模型,如圖1 所示。傘衣采用SHELL 單元,傘繩及加強帶采用BEAM 單元。

圖1 折疊狀態(tài)的阻力傘模型Fig.1 Model of drag parachute in folded state

這里以阻力傘名義直徑D 為參照,選取17D(流向)×10D(展向)×10D(法向)的矩形計算域進行計算。為了準確捕捉噴流流場以及阻力傘阻力特性,對發(fā)動機噴流位置以及阻力傘位置進行局部加密,流體域網(wǎng)格總數(shù)約為500 萬,網(wǎng)格如圖2所示。

圖2 發(fā)動機噴流中心垂直截面網(wǎng)格(充氣狀態(tài))Fig.2 Grid of vertical section at engine center(Inflated state)

對于軍用發(fā)動機,以某典型軍用發(fā)動機為例,其工作狀態(tài)主要可以分為:最大狀態(tài)(全加力狀態(tài))、過渡狀態(tài)(不加力最大推力)、額定/最大連續(xù)狀態(tài)、慢車狀態(tài)等工作狀態(tài)。其中,額定/最大連續(xù)狀態(tài)下推力為最大推力的85%~90%,連續(xù)工作時間小于30~60 min,該工作狀態(tài)主要用于起飛/爬升及高速飛行;慢車狀態(tài)推力為最大推力的3%~5%,主要用于下滑/進場著陸、地面滑行、地面待機等情況。某型號發(fā)動機噴流參數(shù)如表1 所示。發(fā)動機噴流作用下的阻力傘開傘主要包括中斷起飛開傘和發(fā)動機未停機地面滑行開傘,這些工況下發(fā)動機一般處于慢車狀態(tài)至額定狀態(tài)之間。文中主要針對發(fā)動機慢車狀態(tài)至額定狀態(tài)下的阻力傘特性進行研究,選取對應發(fā)動機噴口的出口速度uj范圍為250~500 m/s。

為了模擬發(fā)動機飛機滑行過程中發(fā)動機噴流對阻力傘阻力特性的影響,文中空氣來流以及發(fā)動機噴流采用速度入口邊界條件,來流速度為78 m/s,發(fā)動機噴流速度uj分別為0、250、350、500 m/s。地面為滑移壁面邊界條件,滑移速度與來流速度保持一致。其余各邊界均為無反射邊界條件。傘繩末端節(jié)點約束在噴流入口上方位置,該位置由實際發(fā)動機噴口中心與傘系統(tǒng)連接點的距離決定。

表1 發(fā)動機噴流參數(shù)Table 1 Parameters of engine nozzle jet

3 計算結果與分析

3.1 發(fā)動機噴流對阻力傘阻力特性的影響

這里以飛機滑行速度78 m/s 為例,分別對發(fā)動機噴流速度分別為0、250、350、500 m/s 的噴流影響進行仿真研究。充氣過程傘衣外形變化如圖3 所示,來流氣體在傘衣頂部聚集,隨后帶動傘衣迅速膨脹。經(jīng)過一段時間的耦合作用后達到穩(wěn)定充滿狀態(tài),充滿后的外形不再有大幅度變化。阻力傘充滿狀態(tài)不同噴流速度下的傘衣表面應力分布云圖如圖4 所示。由圖4 可知,隨著噴流速度的增大,傘衣表面應力也隨之增大。

圖3 阻力傘充氣過程傘衣外形變化Fig.3 Changes in shape of drag parachute canopy dur-ing inflation process

為了進一步定量研究發(fā)動機噴流對阻力傘阻力特性的影響,圖5 給出了不同噴流速度下阻力傘開傘過程中阻力動載F 隨時間的變化曲線,表2 給出了不同噴流速度下阻力傘動載峰值Fmax以及充滿狀態(tài)穩(wěn)態(tài)載荷(穩(wěn)態(tài)平均值)F-的變化情況。通過分析可知,當無噴流時,阻力傘動載峰值出現(xiàn)在0.089 s,動載峰值約為126.4 kN;當噴流速度為250 m/s 時,動載峰值出現(xiàn)在t=0.080 s,動載峰值約為142.0 kN;當噴流速度為350 m/s 時,動載峰值出現(xiàn)在t=0.074 s,動載峰值約為143.4 kN;當噴流速度為500 m/s 時,動載峰值出現(xiàn)在t=0.070 s,動載峰值約為162.2 kN。在無噴流情況下,阻力傘充滿后,穩(wěn)態(tài)載荷約為64.1 kN,有發(fā)動機噴流影響時,阻力傘穩(wěn)態(tài)載荷有所增大,但當噴流速度為250 m/s 時,穩(wěn)態(tài)載荷約為77.8 kN,穩(wěn)態(tài)載荷增加了21%;當噴流速度為350 m/s 時,穩(wěn)態(tài)載荷約為97.0 kN,穩(wěn)態(tài)載荷增加了51%;當噴流速度為500 m/s 時,穩(wěn)態(tài)載荷約為114.8 kN,穩(wěn)態(tài)載荷增加了79%。可見,隨著發(fā)動機噴流速度的增大,阻力傘開傘動載峰值不斷增大,并且動載峰值出現(xiàn)的時間不斷提前,充氣時間減小,穩(wěn)定充氣階段的動載值也會隨之增大。同時在充氣初始階段(t=0.05 s之前),開傘動載變化曲線因傘繩拉直的原因出現(xiàn)小幅度波動,并且噴流速度越大,其動載曲線在峰值前后的波動變化更劇烈,這會對阻力傘充氣的穩(wěn)定性和安全性產(chǎn)生一定的影響。因此,在阻力傘設計過程中,需要充分考慮發(fā)動機噴流帶來的影響。

表2 發(fā)動機噴流影響下的阻力傘阻力特性Table 2 Drag characteristics of drag parachute under different jet velocities

3.2 發(fā)動機噴流對阻力傘穩(wěn)定性的影響

穩(wěn)定性是阻力傘的另一個重要性能參數(shù),如果工作過程中阻力傘的擺動幅度過大,可能會使得阻力傘拖地,進而出現(xiàn)破損。因此,文中選取阻力傘在垂直方向的擺動距離進行分析研究。圖6 給出了不同噴流速度下,十字形阻力傘頂點(即徑向加強帶在傘頂孔處的交匯點)在開傘充氣過程中垂直方向(Y 方向)上的位移隨開傘過程的變化曲線。從圖中可以看出,由于發(fā)動機噴流的高速吹襲作用,阻力傘會向下偏斜,當噴流速度為250 m/s 時,阻力傘頂點向下偏斜位移為0.41 m;當噴流速度為350 m/s 時,阻力傘頂點向下偏斜位移為0.65 m;當噴流速度為500 m/s 時,阻力傘頂點向下偏斜位移為0.73 m 。通過對比發(fā)現(xiàn),由于發(fā)動機噴流噴口中心位置高度H 為1.85 m,噴口位置在阻力傘下方(阻力傘初始中心位置高度H 為3 m),高速的發(fā)動機噴流氣流主要吹襲到十字形阻力傘靠近地面的傘衣區(qū)域,傘衣靠近地面內側壓力偏大,傘衣會產(chǎn)生向下的偏向力。同時,高速氣流流經(jīng)傘與地面中間時,會產(chǎn)生一定的“狹管效應”,使得靠近地面?zhèn)阋峦鈧鹊膲毫档停M一步增加傘衣向下的偏移。在阻力傘頂點向下位移到最大位置后,在傘繩向上分力的作用下又迅速向上偏移,形成一種上下往復的周期性振蕩,并且擺動的幅度隨著噴流速度的增大而增大,同時偏斜發(fā)生的時間也隨著噴流速度的增大而前移。所以發(fā)動機噴流速度對阻力傘的穩(wěn)定性有一定的影響,噴流速度越大,阻力傘擺動振幅越大,阻力傘穩(wěn)定性越差。

圖6 不同噴流速度下阻力傘垂直方向擺幅隨開傘過程的變化Fig.6 Change of vertical displacement of drag para-chute during inflation process under different jet velocities

3.3 發(fā)動機噴流作用下的阻力傘流場特性分析

為了進一步研究發(fā)動機噴流對阻力傘阻力特性以及擺動穩(wěn)定性的影響,文中對不同發(fā)動機噴流情況下的阻力傘充滿狀態(tài)流場特性進行研究。圖7 為充滿狀態(tài)t= 0.2 s 時刻不同噴流速度中心截面上的速度分布。由圖可知,當發(fā)動機尾噴存在時,阻力傘傘前流場的速度比無發(fā)動機尾噴時要大,并且隨著噴流速度的增加,阻力傘前的流場速度也逐漸增大。發(fā)動機噴流吹襲到阻力傘傘衣偏下位置,會導致阻力傘傘衣膨脹程度更大。同時,噴流氣流在流經(jīng)阻力傘傘衣與地面形成的狹小通道時被加速,產(chǎn)生一定的“狹管效應”,使得傘衣向地面偏移,這種氣流加速的“狹管效應”會隨著發(fā)動機噴流速度的增加而越發(fā)明顯。

圖7 不同噴流速度下阻力傘充滿狀態(tài)垂直中心截面速度云圖Fig.7 Velocity distribution of vertical center section of drag parachute in inflated state under different jet velocities

通過阻力傘充滿狀態(tài)垂直中心截面上的速度分布云圖可以對發(fā)動機噴流對阻力傘流場特性的影響有一個定性的認識。為了進一步定量分析阻力傘流場特性,文中選取了發(fā)動機噴流中心線(H=1.85 m),阻力傘初始中心線(H =3 m)以及阻力傘前x=3、6、9 m 三條垂直線進行分析,具體位置如圖8 所示。不同噴流速度下發(fā)動機噴流中心線上的速度uce分布如圖9 所示。由圖可知,高速的氣流從發(fā)動機噴出后快速衰減,經(jīng)過阻力傘后速度出現(xiàn)負值,產(chǎn)生一定的回流區(qū),隨后速度又逐步回升,最后與來流速度基本一致。通過不同噴流速度之間的對比可以發(fā)現(xiàn),噴流速度越大,阻力傘之前的來流速度也越大,這也是噴流狀態(tài)下阻力傘阻力增大的原因。

圖8 阻力傘充滿狀態(tài)垂直中心截面特征線位置示意圖Fig.8 Schematic diagram of characteristic line position of vertical center section in inflated state

圖9 阻力傘充滿狀態(tài)下發(fā)動機噴流中心線上的速度變化曲線(H=1.85 m)Fig.9 Speed curve on centerline of engine nozzle in inflated state(H=1.85 m)

圖10 阻力傘充滿狀態(tài)下傘中心線上的速度變化曲線(H =3 m)Fig.10 Speed curve on centerline of drag parachute in inflated state(H =3 m)

不同噴流速度下阻力傘中心線上的速度ucp分布如圖10 所示。由圖可知,起始速度為來流速度,隨后由于高速噴流的卷吸作用,速度逐步升高,隨后又快速衰減,經(jīng)過阻力傘后同樣會出現(xiàn)回流區(qū),最后氣流速度逐步回升至來流大小。此高度上的阻力傘前來流速度會略小于噴流中心線上的氣流速度,但總體上仍大于無噴流狀態(tài)下的傘前速度。

阻力傘傘前不同垂直線上(x=3、6、9 m)上的速度分布如圖11 所示。由圖可知,噴流速度沿噴流中心衰減,且隨著流向距離的增加,噴流的速度也逐步減小,噴流的影響區(qū)域主要集中在0~4 m的范圍內,阻力傘大部分傘衣都在該范圍之內,因此,噴流對阻力傘的性能會產(chǎn)生影響。通過x=9 m 處氣流速度隨高度的變化曲線可以發(fā)現(xiàn),在靠近地面(阻力傘與地面之間)部分的氣流速度有明顯提升,這進一步佐證了氣流使得傘衣向下偏轉。

圖11 阻力傘充滿狀態(tài)時不同垂直線上的速度變化曲線(x=3,6,9 m)Fig.11 Speed curve at different distances in inflated state (x=3,6,9 m)

圖12 給出了充滿狀態(tài)t=0.2 s 時刻不同噴流速度中心截面上的壓力分布。由圖可知,在阻力傘工作過程中,傘衣內側壓力升高,阻力傘后方會形成一定的低壓區(qū),當存在發(fā)動機噴流時,阻力傘內側的高壓區(qū)壓力變大,且傘衣內側下方的壓力偏高,隨著噴流速度的增大,這種效果更加明顯,這種上下不對稱的傘衣壓力分布使得阻力傘產(chǎn)生向下偏轉力,使得阻力傘出現(xiàn)上下偏轉的現(xiàn)象。圖13~14 分別給出了發(fā)動機噴流中心線和阻力傘初始中心線上的壓力(Pce,Pcp)變化曲線。由圖可知,氣流經(jīng)過阻力傘后,壓力急速降低,出現(xiàn)一定的低壓區(qū)。通過對比不同噴流速度下的壓力分布可知,在傘前以及傘內側的壓力會隨著噴流速度的增大而增大,從而使得傘衣的阻力更大,傘衣充滿程度更明顯。

圖12 不同噴流速度下阻力傘充滿狀態(tài)垂直中心截面壓力云圖Fig.12 Pressure distribution of vertical center section of drag parachute in inflated state under differ-ent jet velocities

圖13 阻力傘充滿狀態(tài)時發(fā)動機噴流中心線上的壓力變化曲線(H=1.85 m)Fig.13 Pressure curve on centerline of engine nozzle in in-flated state(H=1.85 m)

圖14 阻力傘充滿狀態(tài)時傘中心線上的壓力變化曲線(H=1.85 m)Fig.14 Pressure curve on centerline of drag parachute in in-flated state(H=1.85 m)

4 結 論

本文對不同發(fā)動機噴流速度下的阻力傘開傘過程進行數(shù)值仿真研究,分析了不同噴流速度對阻力傘阻力特性、阻力傘穩(wěn)定性以及流場特性的影響,主要結論如下:

(1)發(fā)動機噴流會使得阻力傘傘前的流體速度變大,導致阻力傘開傘動載峰值出現(xiàn)時刻前移,開傘動載峰值變大,充滿狀態(tài)的阻力傘穩(wěn)態(tài)載荷變大。

(2)噴流速度為250 m/s 時,穩(wěn)態(tài)載荷增加21%;當噴流速度為350 m/s 時,穩(wěn)態(tài)載荷增加51%;當噴流速度為500 m/s 時,穩(wěn)態(tài)載荷增加79%。

(3)發(fā)動機噴流會使得傘衣內側下方的壓力偏大,傘衣壓力分布不對稱,從而使得阻力傘發(fā)生上下擺動,且噴流速度越大,阻力傘擺動振幅越大,阻力傘穩(wěn)定性越差。

猜你喜歡
發(fā)動機
元征X-431實測:奔馳發(fā)動機編程
2015款寶馬525Li行駛中發(fā)動機熄火
2012年奔馳S600發(fā)動機故障燈偶爾點亮
發(fā)動機空中起動包線擴展試飛組織與實施
奔馳E200車發(fā)動機故障燈常亮
奔馳E260冷車時發(fā)動機抖動
新一代MTU2000發(fā)動機系列
2013年車用發(fā)動機排放控制回顧(下)
VM Motori公司新型R750發(fā)動機系列
發(fā)動機的怠速停止技術i-stop
主站蜘蛛池模板: 日韩免费成人| 国产97视频在线| 试看120秒男女啪啪免费| 国产精品99r8在线观看| 伊人久久大香线蕉综合影视| 国产原创第一页在线观看| 色妺妺在线视频喷水| 又猛又黄又爽无遮挡的视频网站| 欧美午夜在线视频| 国产午夜精品一区二区三区软件| 热思思久久免费视频| 国产精品lululu在线观看| 野花国产精品入口| 97国产在线视频| 免费不卡视频| 国产精品真实对白精彩久久| 九九热在线视频| 国产成人乱码一区二区三区在线| 免费毛片视频| 国产91丝袜在线播放动漫 | 蜜臀AV在线播放| 91精品小视频| av在线无码浏览| 亚洲性色永久网址| 99热国产在线精品99| 亚洲久悠悠色悠在线播放| 国产无码高清视频不卡| 亚洲国产欧美国产综合久久| 亚洲av无码久久无遮挡| 国产00高中生在线播放| 亚洲色图欧美一区| 中文字幕亚洲专区第19页| 91无码国产视频| 亚洲成肉网| 成人国产免费| 99久久免费精品特色大片| 99久久精品国产自免费| 亚洲精品国产成人7777| 爆乳熟妇一区二区三区| 无码免费试看| 亚洲成人网在线播放| 浮力影院国产第一页| 麻豆AV网站免费进入| 成人在线观看一区| 国产91视频免费观看| 久久久久久尹人网香蕉| 中日韩一区二区三区中文免费视频 | 久久综合亚洲鲁鲁九月天| 午夜人性色福利无码视频在线观看| 伊人福利视频| 亚洲色图在线观看| 国产成人无码Av在线播放无广告| 在线国产综合一区二区三区 | 久久久久久高潮白浆| 亚洲第一成年网| 天天色天天综合网| 精品久久777| 欧美精品v| 欧美日本在线观看| 99精品视频九九精品| igao国产精品| 精品少妇人妻一区二区| 精品无码视频在线观看| 99ri精品视频在线观看播放| 国产在线观看高清不卡| 99久久精品视香蕉蕉| 亚洲视频在线观看免费视频| 国产一区二区福利| 熟女日韩精品2区| 亚洲视频色图| 亚洲三级a| 久久久久无码精品国产免费| 无码AV高清毛片中国一级毛片 | 欧美三级自拍| 久久久久免费精品国产| 日韩大乳视频中文字幕| 91精品视频在线播放| 91精品国产一区自在线拍| 国产成人精品免费视频大全五级 | 亚洲无码熟妇人妻AV在线| 在线色国产| 色视频国产|