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

超聲速下盤縫帶傘不同收口方式的氣動特性

2024-05-09 10:16:38代雨柔李健薛曉鵬榮偉
航空學報 2024年7期

代雨柔,李健,*,薛曉鵬,榮偉

1.北京空間機電研究所,北京 100094

2.中國航天科技集團有限公司 航天進入減速與著陸技術實驗室,北京 100094

3.中南大學 自動化學院,長沙 410083

盤縫帶降落傘在超聲速低密度條件下具有優良的阻力性能和穩定性能,結構簡單可靠,常用于火星探測器進入減速、地球高空科學探測等用途的減速傘[1]。以火星探測為例,美國從1976—2020 年共成功進行了9 次著陸探測,均使用了盤縫帶降落傘[2-5]。在中國天問一號火星探測著陸巡視器進入減速與著陸(Entry,Descent and Landing,EDL)階段,超聲速盤縫帶降落傘也起到了至關重要的作用[6]。

隨著中國航天事業的發展,超聲速降落傘應用領域在不斷擴展,性能指標在不斷提升,主要表現為被減速物體的質量不斷增加,開傘速度不斷提高,對降落傘的開傘力、質量與重量的約束越來越嚴苛。這就要求降落傘在大幅增加面積和提高強度的同時,要盡可能降低開傘力,減小重量和體積。

收口技術是解決上述工程問題的有效途徑之一。該方法在充氣過程中可有效控制降落傘的阻力特征,使降落傘逐級打開并可以抑制降落傘的過度充氣。該技術已在亞聲速領域得到了很好的應用,神舟系列飛船的減速傘和主傘均采用1 次收口開傘方法,以控制各級開傘過載不超過規定的限值要求[7]。但是在超聲速領域,該技術還未得到工程應用,國內外僅做過少量的空投試驗。超聲速高空降落傘試驗室(SHAPE)在馬赫數>2.5 的情況下,使用細長體和鈍體前體對收口降落傘進行了飛行試驗。相比于細長體,鈍體后面的收口降落傘由于前體尾流的影響,出現了破損的情況[8-9]。數值模擬方面未見相關文獻公開發表。由此可見,現階段缺乏關于超聲速收口盤縫帶傘有效的數值模擬方法,和關于不同收口比下的收口盤縫帶傘氣動表現的系統性研究。所以,關于盤縫帶傘超聲速收口技術的應用需求是迫切的,現有研究工作不滿足未來工程應用需求。

對比亞聲速收口降落傘,超聲速收口降落傘所面對的主要問題是:流場更加復雜,前體與降落傘間的干擾加劇,使得傘衣出現強烈的“呼吸”振動現象,有可能造成傘衣塌陷、結構破壞,導致減速失敗。本文以不同收口方式的盤縫帶傘為研究對象,采用流固耦合方法進行數值模擬,對比不同收口方式的盤縫帶傘的開傘和呼吸過程,總結不同收口比下的氣動特性規律,為超聲速收口盤縫帶傘的工程應用提供技術支撐。

1 研究對象

1.1 收口盤縫帶傘

盤縫帶降落傘是開縫傘的一種,傘衣由一個中心有通氣孔的圓“盤”和圓筒型“帶”組成。降落傘收口的常用方法是使用收口繩進行收口。該方法通過將收口繩用收口環固定在傘衣周向位置,從而對傘衣的充氣過程起一個限制作用。同時,收口繩上也會放置切割器,用于在規定時間解除收口[1]。由于盤縫帶傘的傘衣結構特殊,收口繩可以安裝在“盤”或者“帶”的邊緣,即盤縫帶傘的2 種收口方式,“帶”收口和“盤”收口,如圖 1所示[10]。

1.2 降落傘參數設置

使用美國“海盜號”的4%縮比傘為研究對象,該傘已有公開發表的相關試驗數據,可用于對比分析。其具體尺寸和相應尺寸分別如表 1和圖 2 所示[11]。

表1 盤縫帶傘結構尺寸Table 1 Structural size of disk-gap-band parachute 單位:m

數值模擬過程中,降落傘的傘衣材料設置為薄殼單元,傘繩和加強帶選用離散的梁/索單元,傘衣與傘繩處于初始拉直狀態。返回艙外形參考EXOMars 計劃的鈍體航天器模型[12],由一個半錐角70°的鈍頭前體和一個半錐角為45°的圓錐后體組成,使用薄殼單元搭建。返回艙的單元空間固定。整體結構如圖3 所示。

1.3 降落傘收口參數設置

工程上常用收口直徑比(也稱收口比)來定量描述收口程度,該參數為收口繩長度與傘衣底邊長度之比,也可以表示為收口繩直徑和傘衣名義直徑的比[13]。使用收口繩直徑和傘衣名義直徑比定義收口比,對馬赫數分別為1.6和1.8 時不同收口方式的不同收口比進行數值模擬,其中以馬赫數1.8 環境下的數值模擬結果作為主要研究對象,以馬赫數1.6 環境下的結果作為驗證對象。具體研究算例如表 2 所示。

表2 研究算例Table 2 Research examples

2 來流條件及研究方法

2.1 來流條件

以未來火星低密度大氣減速為需求背景,為更好地貼近真實的火星大氣環境,參照火星科學試驗室(Mars Science Laboratory,MSL)的探測巡視器著陸火星過程中其速度與離火星地面高度的對應關系[14]和NASA 給出的根據火星全球勘探者(Mars Global Surveyor)對火星大氣進行測量得到的擬合公式[15],計算出本研究使用的火星大氣參數,如表 3 所示。

表3 仿真使用來流條件Table 3 Freestream condition for simulation

2.2 網格無關性驗證及流場域模型

選取單個降落傘進行網格無關性驗證,對單元長度為0.05、0.04、0.03、0.02 m 的流場網格進行數值模擬計算。計算得到的降落傘的阻力系數如表 4 所示,以0.02 m 為基準算例,得到各算例的誤差。給出不同網格長度與阻力系數的關系圖,如圖 4 所示。與網格大小為0.02 m 的算例對比,各阻力系數的誤差均在10%以下,根據圖 4 所示結果,在網格數為0.02 m時,計算結果已收斂。同時對比圖 5 所示0.04、0.02 m 網格的局部速度云圖,觀察得到,盡管2 種不同大小的網格都能夠很好地展示降落傘后的不定常寬尾流。但相比于0.04 m 的網格,0.02 m 能夠更好地捕捉到降落傘頂孔的高速射流,流場細節更加清晰。所以在綜合考慮計算精度與時間成本的情況下,網格長度選取0.02 m。

表4 網格無關性驗證Table 4 Grid independence verification

所取得的流場域沿降落傘軸向、法向、側向的長度分別為5 m×3 m×3 m。最終六面體的流場域劃分為約526 萬個網格,建立的流場計算域及仿真模型的整體網格結構如圖 6 所示。由于傘衣面積小,開傘速度快,仿真流場采用理想氣體模型且采用如表 3 所示的恒定大氣密度和壓強進行仿真。流場頂端和壁面添加無反射約束,用以模擬遠場狀態,使仿真結果與實際無限質量開傘情況更為接近。

2.3 研究方法

降落傘流固耦合的研究一直是一個挑戰,為了模擬柔性降落傘的開傘和開傘后傘衣呼吸的全部過程,需要對流場中的激波和尾流進行良好的模擬。時-空守恒元解元(CE/SE)方法是一種物理概念清晰、計算精度高、格式構造簡單的流體力學方法[16],具有廣闊的發展前景。該方法由Chang 提出,由劉海燾等進行改進,可求解間斷流場CE/SE 格式的三維Navier-Stokes(N-S)方程,具有較高的分辨率。求解方法往往分為以下4步[17-21]:

步驟1在完全氣體的一維情況下,給出量綱為1 的非定常N-S 方程,并定義不同空間位置的離散點,如圖 7 所示的A1、Q、A2[21]。之后,在一維離散的基礎上將時間作為額外的空間維度,構建二維歐拉空間域E2,并建立該歐拉空間的坐標,令x1=x,且x2=t,則

步驟2利用高斯散度定理,將時間和空間完全統一,微分方程為

式中:S(V)是E2中任意區域V 的邊界。

步驟3建立SE,沿圖 7 所標定的時間空間的內部邊界,沿SE 域對流體變量進行近似的泰勒級數展開式進行逼近,令

同理可得:

步驟4將步驟2 所得式(2)用式(7)逼近,并將式(3)~式(5)代入,得到式(8):

此時在每個網格點的變量只有um和umx,將圖 7 所示的2 個守恒元CE-和CE+的變量代入式(8),得出計算結果。

所使用的DUALCE/SE 方法(dual-mesh CE/SE,又可簡稱為雙網格CE/SE 方法),對傳統的CE/SE 算法進行了小范圍的改進。傳統CE/SE方法求解器的流體變量都在網格的中心點上,而DUALCE/SE 方法采用雙網格交替在流體單元中心和流體單元節點上求解流體狀態變量。在采用相同的網格數下,該方法更為精確[22]。

試驗中,開傘的動態過程采用DUALCE/SE方法計算流場,固體結構分析使用有限元方法,兩者交替進行計算。流固耦合部分使用浸入邊界法來處理。浸入邊界法是一種數學建模方法,該方法在生成網格時不需要生成復雜的貼體網格,只需要生成正常的笛卡爾網格[23]。而針對充氣后的穩定形態,采用計算流體力學方法進行分析。

3 仿真模型準確性驗證

此前,已有部分學者使用時-空守恒元解元方法對超聲速降落傘的開傘過程進行了仿真,并得到了較為準確的結果[24-26]。本節的研究對象是在馬赫數1.6 的火星大氣條件下,對上述降落傘在不收口情況下展開的數值模擬結果,從阻力系數、傘衣充滿過程的外形變化方面和已有試驗結論進行對比。另外,超聲速高空降落傘試驗室(SHAPE)在1971 年進行了超聲速收口降落傘的飛行試驗,本文將對比收口降落傘的阻力系數、投影面積及氣動外形的相關數據。

如圖 8 所示,在馬赫數1.6時,數值模擬計算出的阻力系數平均值為0.60。圖 9 所示為NASA 綜合風洞、高空開傘、飛行試驗結果,給出的盤縫帶傘阻力系數隨馬赫數的變化曲線,中心線為標稱值,上下為偏差范圍[27]。其中,在馬赫數1.6時,阻力系數的高值為0.82,中值為0.66,低值為0.50。選取中值作為參考,數值模擬所得阻力系數位于預測阻力系數之間,并接近中值。

圖10 所示是先進超聲速降落傘充氣研究試驗(ASPIRE)在2019 年所做的盤縫帶傘高空開傘試驗的降落傘的展開仰拍圖[28]。圖中從左至右依次為傘繩拉直降落傘開始展開的時刻、降落傘張開過程中的時刻和降落傘首次張開的時刻。從圖中可以看出,雖然展開時間有所不同,但是3 組展開圖都是盤的部分先打開,之后帶的部分再打開。圖 11 所示是數值模擬的結果,圖中從左到右所展示的降落傘外形的順序與圖 10 相同,數值模擬結果與實際情況一致。

圖1 盤縫帶傘的不同收口方式[10]Fig.1 Different reefing ways of disk-gap-band parachute[10]

圖2 盤縫帶傘結構圖[11]Fig.2 Disk-gap-band parachute structure[11]

圖3 盤縫帶傘模型Fig.3 Disk-gap-band parachute model

圖4 不同網格大小的阻力系數Fig.4 Drag coefficients of different mesh sizes

圖5 速度云圖Fig.5 Velocity contour

圖6 仿真模型網格結構Fig.6 Grid structure of simulation model

圖7 CE/SE 一維求解空間域[21]Fig.7 CE/SE one-dimensional solution space domain[21]

圖8 阻力系數仿真時間歷程曲線Fig.8 Simulation time history curve of drag coefficient

圖9 不同速度下阻力系數范圍Fig.9 Range of drag coefficient at different speeds

圖10 探空火箭降落傘展開圖[28]Fig.10 Inflation process of sounding rocket parachute[28]

圖11 數值模擬降落傘展開過程Fig.11 Inflation process of numerical simulated parachute

對于超聲速收口盤縫帶傘的準確性驗證,超聲速高空降落傘試驗室(SHAPE)的飛行試驗得到了在收口比約為29%的情況下,不同馬赫數下的阻力系數,如圖12 所示[8]。其中馬赫數為1.8時,阻力系數為0.211。另一試驗項目得到了同大小的盤縫帶傘不收口情況下,馬赫數為1.8 時的阻力系數為0.535[29]。將兩者對比,得到在收口情況下阻力系數占未收口情況下阻力系數的39.4%。該文章也提供了降落傘在馬赫數為1.8時的阻力面積比(傘衣收口狀態的阻力面積與完全張滿傘衣的阻力面積之比[13])為50%。將該數據與數值模擬結果對比,見表 5。試驗結果與數值模擬結果誤差小于10%。同時該試驗室也提到,在降落傘展開之后,傘盤部分基本保持穩定,呼吸作用主要為傘帶部分發生膨脹和收縮。本文數值模擬的結果與該文描述一致相同,具體氣動特征見本文4.1節。

圖12 飛行試驗所得阻力系數[8]Fig.12 Drag coefficient of flight test[8]

圖13 盤收口盤縫帶傘相關結果Fig.13 Relevant result of disk-gap-band parachute with mid-gore reefing

圖15 盤收口降落傘壓力云圖Fig.15 Pressure contour of disk-gap-band parachute with mid-gore reefing

圖16 帶收口盤縫帶傘相關參數Fig.16 Relevant parameter of disk-gap-band parachute with skirt reefing

圖17 帶收口降落傘不同部分阻力系數Fig.17 Drag coefficients of different parts of disk-gapband parachute with skirt reefing

圖18 帶收口降落傘壓力云圖Fig.18 Pressure contour of disk-gap-band parachute with skirt reefing

圖19 不同收口比下的參數比值Fig.19 Parameter ratio of different reefing ratio

圖20 降落傘點位圖Fig.20 Different positions on parachute

由不收口降落傘的開傘外形和阻力系數以及收口降落傘的阻力系數,投影面積和氣動外形的對比,可以證明:本研究采用的數值模擬方法是有效的,計算結果是可靠的。

表5 參數對比Table 5 Parameter comparison

4 不同收口方式的影響

4.1 不同收口方式開傘和呼吸方式對比

4.1.1 盤收口

本節總結了收口比為30%的情況下收口繩的受力與時間的關系圖、降落傘上不同點的投影面積與時間的關系圖以及降落傘在呼吸過程中不同時間的仰視圖,如圖 13 所示。

在盤收口的開傘和呼吸過程中,降落傘的盤和帶表現出不同的氣動特點。降落傘的收口繩在約0.06 s 時第1 次收緊,在此之后的0.02 s內,降落傘盤部分點位的投影面積達到最大值,在這之后,盤部分點位的投影面積幾乎保持不變。而帶部分在收口繩第1 次收緊之后,該部分點位的投影面積依舊持續增大,直到約0.11 s 達到最大值。之后在0.192 s時,帶部分點位的投影面積到達第1 個波谷位置。在這個位置,收口繩受力也有一定的減小,這說明帶的收縮運動也微小地帶動了盤邊緣的收縮。之后在約0.21 s 的位置收口繩的受力再次增加,隨后帶的部分也再次張開,投影面積在0.24 s 時達到了第2 次膨脹的最大值,并在0.26 s 時出現了小范圍的波動。從波峰波谷對應的傘衣結構仰視圖也可以看到,在盤收口降落傘開傘的時候,盤的部分先充滿,之后帶的部分再逐漸張開。在降落傘呼吸的時候,盤的部分基本保持充滿狀態,而帶的部分發生不規則的膨脹和收縮。

流固耦合計算對流場中的激波、壓力捕捉精度有限。為了分析前體尾流與傘前激波的相互作用過程,使用計算流體力學的方法,提取收口比為30%的降落傘穩定后的平均外形,并假設該降落傘為不透氣剛體進行計算。算法基于密度求解器,采用有限體積法進行空間離散,使用三維可壓縮理想氣體N-S 方程為控制方程,仿真采用層流進行計算。圖 14 所示為使用計算流體力學方法所得到的盤部分、帶部分和整體受到的阻力系數。圖 15 所示為時間取0.062~0.070 s 這一個壓力周期內降落傘周圍的壓強云圖。在0.062 s時,降落傘阻力系數在波動周期的波谷,傘內處于一個相對低壓狀態,隨著時間的推進,傘前激波順流向推進,導致降落傘傘內的壓力升高。當傘內壓強增高到一定值后,傘內高壓會推動傘前激波逆向移動,致使傘內壓強再次降低,繼而完成一個周期。從云圖中可明顯看到,在不考慮織物透氣性的情況下,收口繩和降落傘盤的部分形成了一個只有頂孔透氣的傘包,傘盤中積累的高壓只能通過降落傘盤和帶之間的縫隙流出,所以該結構特點導致降落傘的傘盤部分始終處于充氣狀態。另一方面,雖然降落傘“縫”的位置形成了一個泄壓口,但由于傘后的膨脹波位置位于帶的外邊緣,導致帶的部分始終存在較高的內外壓差,從而使得該部分在開傘時會有一個較大的波動。

4.1.2 帶收口

收口比為30%的帶收口的收口繩受力和阻力系數如圖 16 所示,為了主要關注收口降落傘的開傘和呼吸過程,本節只截取上述參數0~0.4 s時間的數值。可以觀察到,收口繩受力在0.074 s時達到了第1 個峰值,但是最大投影面積還處于上升階段,在0.088 s 左右投影面積上升到最大值。收口繩受力在0.174 s 時達到了第2 次最大值,而投影面積則在0.182 s 時達到了第2 次的最大值。投影面積的最大值和收口繩受力的最大值的時間較為統一,但是投影面積的峰值要比收口繩受力的峰值出現的稍晚,這主要是因為在帶收口的降落傘充氣的時候,盤的邊緣部分要先于帶的部分張開,然后通過傘繩帶動帶的部分向四周張開,在帶的部分張開至收口繩完全撐開的時候,傘盤邊緣繼續張開直到最大程度后開始回彈。隨著時間的推進,收口繩受力和投影面積也越來越趨于穩定。從不同點的投影面積可以看出,除第1 次投影面積達到最高值的時候外,其他情況下,30%帶收口降落傘的最大投影面積位于N3 的位置,N4和N2 位置的投影面積大小隨降落傘的呼吸作用來回交替變換,所以帶收口的開傘和呼吸過程是整個降落傘一起展開和收縮的過程。

同樣,使用計算流體力學的方法,對30%帶收口降落傘的平均氣動外形進行了分析。圖 17所示為使用計算流體力學方法所得到的盤、帶和整體受到的阻力系數,并對自0.048~0.058 s 的降落傘周圍壓強部分進行分析,如圖 18 所示。從壓力云圖中可以看到壓力增強與激波移動的周期。但是和盤收口不同的是,帶收口的膨脹波在盤的邊緣,而不是帶的邊緣。這種特殊的結構使得帶的部分被包裹進了高壓區域,并使得該部分并沒有起到相應減速的作用,而各部分的阻力系數圖也同樣證明了這點。并且,從阻力系數圖也可以看到,帶收口的時候,降落傘盤部分的阻力系數和傘衣整體的阻力系數幾乎相同,而帶部分受到的力在正向力與負向力之間來回波動,這也會導致帶部分的不穩定。

4.2 不同收口方式的氣動規律

首先,使用數值模擬方法計算出不同收口方式和收口比的阻力系數和投影面積,為了更好地對比不同數據間的關系,對阻力系數和投影面積做以下處理:將不同收口方式的不同收口比的阻力系數除以不收口降落傘穩定后的阻力系數,同理,將不同收口方式的不同收口比穩定后的平均投影面積除以不收口降落傘穩定后的平均投影面積,得到如圖 19 所示的關系。從圖中可以看出,盤收口方式的阻力系數比和投影面積比隨著收口比的增大而增加,其中在30%~35%的過程中,投影面積占比和阻力系數占比增加的幅度較大,但是35%~40%之間,投影面積并沒有明顯增大。這主要是因為,在35%的收口比時,數值模擬時降落傘發生了較為嚴重的自旋狀況,隨著自旋速度的增加,降落傘的傘衣向離心方向擴展,使得該傘的整個氣動外型變得更加扁平,投影面積也相應增大。從數據分析來看,這種現象也間接使得降落傘的阻力系數有所增大。但整體來看,盤收口投影面積的占比始終比阻力系數的占比要大,并且能夠保持近似同步的增長幅度。

降落傘帶收口方式的結果有所不同,無論是阻力系數的比值還是投影面積的比值,收口比在30%~50%時,降落傘的投影面積和阻力系數改變都不明顯。為了探究該現象的原因,收集了不收口和30%、35%、40%帶收口在降落傘不同位置點位的投影面積大小,計算其對應的近似圓的直徑長度,并與名義直徑做比,點位圖如圖 20 所示,結果如表 6 所示。表中數據表明,在這幾種收口比下,降落傘的最大投影面積的位置(表中加粗數據)由N2 位置逐漸轉移到N3 位置,表明降落傘的外形隨著收口比的增加而變化,但是最大投影面積的值并沒明顯差別,并且跟不收口降落傘的投影面積相差較大。結合圖 19 可以看到,這些投影面積的改變并沒有引起阻力系數的明顯增加。

為了解釋這一現象,從降落傘結構、傘內壓力和傘衣的結構孔隙率方面進行說明。首先,從降落傘的結構方面分析,帶收口時降落傘的最大投影面積的位置位于N3 或N2 點位。所以,如圖 21 所示,帶收口時降落傘的傘衣形成了一個翻折,傘頂到傘底這段長度形成了圖中藍線所示的曲線,這里稱為迎風面和背風面。于是,當降落傘底部的收口繩增長時,也就是降落傘的底部橫向周長增加時,因此而導致的縱向長度的增長會同時被分配到迎風面和背風面2 面上,導致收口繩的增加對降落傘的阻力面積的影響很小。其次,由上文所給的帶收口降落傘的壓力云圖及不同部分的阻力系數可以看到,在帶收口的情況下,傘帶的部分被包裹在了傘前高壓中,這個時候傘帶受到的平行于來流的力非常小,即相當于傘帶的部分沒有為降落傘的阻力提供任何幫助,導致阻力系數增加不明顯。最后,從降落傘的結構透氣性并結合上文所展示的傘內壓力變化來解釋這一現象。降落傘膨脹的本質就是在孔隙率低的地方形成高壓區域。如上文提到盤收口是將降落傘的盤部分圍成了小型的傘包,若將該傘包的結構透氣性等效為該結構的有效孔隙率,結合前文所給參數可以計算出該傘包的有效孔隙率僅為0.77%。而帶收口所形成的整體部分的有效孔隙率為11.98%。因而超聲速情況下,有效孔隙率更大的情況會使得來流更易流出,更難形成較高的內外壓差,也導致降落傘的阻力系數增加不明顯。

表6 不同點位降落傘投影面積及比值Table 6 Projected area and its ratio of parachute at different positions

5 結論

基于流固耦合算法,得到了在收口比相同的情況下不同收口方式開傘和呼吸過程的特點,并得到了不同收口方式不同收口比的阻力系數和投影面積的關系,及其普遍氣動特點。現總結如下:

1)在盤收口降落傘開傘時,盤部分首先充氣,之后降落傘的帶部分再逐漸張開。降落傘呼吸的時候,盤的部分基本保持充滿狀態,帶的部分發生不規則的膨脹和收縮。

2)在帶收口降落傘開傘的過程中,盤的邊緣部分要先于帶的部分張開,然后帶動帶的部分向四周張開。呼吸過程是整個降落傘一起展開和收縮的過程。

3)盤收口方式的投影面積的占比始終比阻力系數的占比要大一些,自30%~50% 的收口比,兩者皆隨著收口比的增大而增大,并保持近似同步的增長幅度。

4)帶收口方式的阻力系數的比值和投影面積的比值在收口比為30%~50%之間并沒有明顯變化,這主要是因為相比起盤收口,帶收口的帶部分沒有起到阻力的作用,且結構孔隙率較大,相比盤收口,傘內外沒有形成較大的壓強。

5)在馬赫數1.8下,收口比相同時,盤收口能夠提供更大的阻力系數。結合文章所給數據,初步建議:盤收口的收口比應在50%以下,因為該收口比以上的收口阻力系數已非常接近全展開情況的阻力系數,收口繩已沒有明顯作用。

主站蜘蛛池模板: 女人18毛片久久| 成人第一页| 999福利激情视频| 亚洲啪啪网| 成人午夜视频网站| 五月天久久综合国产一区二区| 香蕉久人久人青草青草| 亚洲日韩精品欧美中文字幕| 国产69囗曝护士吞精在线视频| 欧美乱妇高清无乱码免费| 欧美一区二区人人喊爽| 亚洲高清国产拍精品26u| 白丝美女办公室高潮喷水视频| 亚洲精选无码久久久| 国产精品偷伦在线观看| 国产麻豆另类AV| 91无码视频在线观看| 国产黄色视频综合| 国产亚洲美日韩AV中文字幕无码成人 | 重口调教一区二区视频| 国产又爽又黄无遮挡免费观看| 在线观看精品国产入口| 欧美国产视频| 国产一级小视频| 亚洲国产成人综合精品2020| 欧美色综合网站| 日韩欧美国产另类| 亚洲欧美国产高清va在线播放| 91po国产在线精品免费观看| 高清乱码精品福利在线视频| 欧美天堂在线| 中文字幕有乳无码| 无码福利视频| 天天躁夜夜躁狠狠躁躁88| 国产小视频a在线观看| 婷婷综合在线观看丁香| 国产微拍精品| 久久99国产精品成人欧美| 国产高清免费午夜在线视频| 夜精品a一区二区三区| 日韩欧美国产成人| 久久黄色免费电影| 天天色综网| 成人av专区精品无码国产| 女人一级毛片| 国产超碰一区二区三区| 91精品久久久久久无码人妻| 欧美中文字幕一区| 五月激情婷婷综合| 尤物特级无码毛片免费| 国产成人亚洲无吗淙合青草| 色妞www精品视频一级下载| 久久免费视频播放| 国内精品自在自线视频香蕉| 精品视频一区在线观看| 国产91视频免费观看| 国产精品无码久久久久久| 一级毛片基地| 麻豆精品久久久久久久99蜜桃| 尤物精品视频一区二区三区| 亚洲中文无码av永久伊人| 国产1区2区在线观看| 99视频精品在线观看| 青草视频久久| 日韩国产精品无码一区二区三区| 色婷婷亚洲综合五月| 一本大道无码日韩精品影视| 福利视频一区| 亚洲福利片无码最新在线播放| 欧美日韩免费| 久久9966精品国产免费| 韩国自拍偷自拍亚洲精品| 熟妇无码人妻| 亚洲欧美精品日韩欧美| 欧美成a人片在线观看| 999在线免费视频| 波多野结衣一区二区三区四区视频| 国产亚洲精久久久久久无码AV| 另类欧美日韩| 欧美精品亚洲精品日韩专区| 国精品91人妻无码一区二区三区| 97免费在线观看视频|