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

類火星環(huán)境中火箭起降過程噴流影響研究

2023-12-14 05:11:22林曉輝李宗儒許常悅
宇航總體技術(shù) 2023年6期
關(guān)鍵詞:環(huán)境研究

林曉輝,秦 曈,李宗儒,杜 濤,許常悅

(1.南京航空航天大學(xué)飛行器環(huán)境控制與生命保障工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室,南京 210016;2.北京宇航系統(tǒng)工程研究所,北京 100076)

0 引言

SpaceX公司于2015年成功發(fā)射的可重復(fù)使用火箭引起了廣泛關(guān)注,與之相關(guān)的超聲速減速推進(jìn)[1-2]、氣動(dòng)面設(shè)計(jì)[3]、發(fā)動(dòng)機(jī)噴流特性[4-6]等研究方向也成為熱點(diǎn)。其中,與發(fā)動(dòng)機(jī)噴流特性相關(guān)的研究對(duì)火箭發(fā)動(dòng)機(jī)的設(shè)計(jì)及火箭回收等具有重要的參考價(jià)值。在過度膨脹的發(fā)動(dòng)機(jī)噴管中,流動(dòng)在噴管壁面出現(xiàn)分離,導(dǎo)致不對(duì)稱的側(cè)向載荷,這可能會(huì)損毀噴管,因此預(yù)測(cè)噴管內(nèi)部流動(dòng)分離點(diǎn)是試驗(yàn)和數(shù)值研究的重點(diǎn)[7]。Nguyen等[8-9]通過試驗(yàn)測(cè)量了處于過度膨脹狀態(tài)的鐘形噴管壁面壓力。他們觀察到自由沖擊分離(Free Shock Separation,F(xiàn)SS)和受限沖擊分離(Restricted Shock Separation,RSS)以及FSS-RSS轉(zhuǎn)換狀態(tài)。他們還觀察到在FSS狀態(tài)下分離線具有不穩(wěn)定性,在RSS狀態(tài)下,噴管的出口存在端部效應(yīng)。Deck等[10-11]對(duì)火箭噴管出口處不對(duì)稱和不穩(wěn)定分離流動(dòng)引起的側(cè)向載荷進(jìn)行了數(shù)值研究,研究發(fā)現(xiàn)出口的流動(dòng)不穩(wěn)定性誘發(fā)了大幅度的低頻側(cè)向載荷。在FSS狀態(tài),不穩(wěn)定性來自超聲速射流和環(huán)境空氣相互作用導(dǎo)致的絕對(duì)不穩(wěn)定性。為了研究RSS狀態(tài)的流動(dòng)特性,Deck[12]使用延遲分離渦模擬(Delayed Detached Eddy Simulation,DDES)方法對(duì)過度膨脹噴管中的流動(dòng)進(jìn)行了數(shù)值仿真。研究結(jié)果表明在RSS狀態(tài)及特定的壓力下,當(dāng)流動(dòng)再附著在出口的時(shí)候,流動(dòng)會(huì)呈現(xiàn)強(qiáng)烈的全局不穩(wěn)定性。

有關(guān)火箭噴管流動(dòng)特性的研究主要集中在對(duì)不同噴管壓力比(Nozzle Pressure Ratio,NPR)下FSS流態(tài)和RSS流態(tài)的機(jī)理分析,對(duì)于大壓力比下噴流的地面效應(yīng)的相關(guān)研究較少。火箭著陸過程中噴流沖擊地面會(huì)導(dǎo)致復(fù)雜的湍流現(xiàn)象,激波、剪切層和邊界層存在強(qiáng)烈的相互作用,隨著著陸距離縮短,這種湍流效應(yīng)更加明顯[12-13]。對(duì)火箭噴流沖擊地面這一過程開展試驗(yàn)的成本巨大且很難獲得準(zhǔn)確的數(shù)據(jù),因此在相關(guān)研究中計(jì)算流體力學(xué)(Computational Fluid Dynamics,CFD)占據(jù)了主導(dǎo)地位。如Deck等[10-11]的數(shù)值研究對(duì)Nguyen等的試驗(yàn)結(jié)果進(jìn)行了補(bǔ)充,進(jìn)一步闡釋了噴管內(nèi)部側(cè)向載荷的來源。Tsutsumi等[14]為了解第一級(jí)火箭發(fā)動(dòng)機(jī)布局對(duì)發(fā)射臺(tái)內(nèi)部流場(chǎng)的影響,對(duì)H3運(yùn)載火箭進(jìn)行了數(shù)值分析,研究發(fā)現(xiàn)來自相鄰發(fā)動(dòng)機(jī)的射流在撞擊發(fā)射臺(tái)后導(dǎo)致上游方向的側(cè)向射流反轉(zhuǎn)。

超聲速噴流沖擊地面的流動(dòng)特性對(duì)火箭著陸具有重要的工程指導(dǎo)意義,在作者已知的公開文獻(xiàn)中,相關(guān)的研究仍然較少。我們之前的研究基于穩(wěn)態(tài)數(shù)值仿真結(jié)果,對(duì)火箭著陸過程中噴管的布局和離地高度的影響展開了分析,發(fā)現(xiàn)火箭噴流沖擊地面具有明顯的非定常特性[15]。為了進(jìn)一步對(duì)火箭噴流地面效應(yīng)的作用機(jī)理進(jìn)行分析,本文采用尺度自適應(yīng)模擬(Scale Adaptive Simula-tion,SAS)方法對(duì)火箭噴流展開了非定常數(shù)值仿真。本文重點(diǎn)關(guān)注了在海平面環(huán)境和類火星環(huán)境中單個(gè)噴管內(nèi)部的流場(chǎng)結(jié)構(gòu)以及噴流和地面的相互作用。相關(guān)的研究成果對(duì)火星著陸具有一定的指導(dǎo)意義。

1 數(shù)值計(jì)算方法

1.1 控制方程

本文采用基于SST兩方程湍流模型的SAS方法對(duì)噴流沖擊地面的非定常特性進(jìn)行研究。SAS方法是一種改進(jìn)的瞬態(tài)雷諾應(yīng)力時(shí)均(Unsteady Reynolds Average Navier-Stokes,URANS)方法,可以動(dòng)態(tài)調(diào)整解析尺度,該方法已經(jīng)成功應(yīng)用于許多高速可壓縮流動(dòng)問題,適合在不穩(wěn)定流動(dòng)條件下解析湍流。SAS方法本質(zhì)上仍是一種URANS方法,因此其控制方程仍為Navier-Stokes(N-S)方程組,在笛卡爾坐標(biāo)系下可壓縮SAS方法的控制方程如下

(1)

(2)

(3)

這里建立的SST-SAS模型為

(4)

(5)

為了實(shí)現(xiàn)SAS的計(jì)算,需要在渦量密度方程(5)中添加SAS源項(xiàng)QSAS

QSAS=

(6)

以上各式中的參數(shù)信息詳見參考文獻(xiàn)[16]。

1.2 計(jì)算細(xì)節(jié)

計(jì)算物理模型采用德國(guó)航空航天中心設(shè)計(jì)的軸對(duì)稱噴管,如圖1所示。相關(guān)研究人員利用該噴管進(jìn)行大量的試驗(yàn),積累了較多的試驗(yàn)數(shù)據(jù),被廣泛用于數(shù)值計(jì)算的驗(yàn)證[8,10,12-13],詳細(xì)參數(shù)見表1。

圖1 火箭發(fā)動(dòng)機(jī)噴管輪廓Fig.1 The contours of the engine nozzle

表1 噴管設(shè)計(jì)參數(shù)

本文采用的計(jì)算域如圖2所示。外流場(chǎng)計(jì)算域直徑為10D(D為火箭噴管出口直徑),噴管出口到地面的距離為2D,地面傾斜角α=4°。當(dāng)噴流沖擊水平地面時(shí)α=0°,其他條件和噴流沖擊傾斜地面時(shí)保持一致。

圖2 計(jì)算域示意圖Fig.2 Diagram of computational region

如圖3所示,噴管內(nèi)部流動(dòng)采用結(jié)構(gòu)網(wǎng)格進(jìn)行計(jì)算,在噴管內(nèi)壁和噴流沖擊壁面區(qū)域局部加密。噴管內(nèi)壁面的第一層網(wǎng)格高度為10-6rt,底部噴流的第一層網(wǎng)格高度為10-4rt,計(jì)算網(wǎng)格總數(shù)為990萬。

圖3 噴管的O形網(wǎng)格Fig.3 O-type grid of the nozzle

1.3 算法驗(yàn)證

為了驗(yàn)證本文采用的SAS方法的正確性,將計(jì)算結(jié)果和文獻(xiàn)[8-9]的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。如圖4所示,噴管壁面的無量綱壓力pw/pa實(shí)驗(yàn)結(jié)果和SAS計(jì)算結(jié)果一致性較好。需要說明的是,由于缺乏流動(dòng)分離前的公開實(shí)驗(yàn)數(shù)據(jù),當(dāng)流向位置y/L<0.24的時(shí)候,本文的計(jì)算結(jié)果和Lijo等[17]的計(jì)算結(jié)果具有相似的流場(chǎng)結(jié)構(gòu)。URANS和大渦模擬(Large Eddy Simulation,LES)等CFD方法在低環(huán)境壓力的應(yīng)用在近年來得到了發(fā)展,其準(zhǔn)確性和適用性被一些公開文獻(xiàn)證實(shí)[18-20]。本文涉及的在火星環(huán)境下的噴流壓比較大,試驗(yàn)難度較高。通過對(duì)比文獻(xiàn)[21]、文獻(xiàn)[22]的數(shù)值結(jié)果發(fā)現(xiàn),SAS計(jì)算結(jié)果和通過差分求解N-S方程與直接模擬蒙特卡羅(DSMC)耦合方法的研究結(jié)果在當(dāng)前類火星環(huán)境下具有相似的流場(chǎng)結(jié)構(gòu)。綜上所述,SAS方法能夠滿足準(zhǔn)確性要求。

圖4 噴管壁面的無量綱壓力pw/pa對(duì)比Fig.4 Non-dimensional pressure (pw/pa)comparison of the nozzle wall

2 計(jì)算結(jié)果與分析

2.1 海平面環(huán)境下地面傾斜度的影響

火箭著陸時(shí)面臨復(fù)雜的地勢(shì)環(huán)境,當(dāng)火箭發(fā)動(dòng)機(jī)噴管和地面的距離較小時(shí),噴流和地面之間存在復(fù)雜的相互作用,在之前的研究已經(jīng)進(jìn)行了初步探索[15]。當(dāng)前工作研究了離地高度為2D時(shí),噴流沖擊地面的流場(chǎng)特性。當(dāng)NPR=19時(shí),噴管內(nèi)部為FSS狀態(tài),即流動(dòng)在噴管內(nèi)部出現(xiàn)分離,在噴管內(nèi)形成了馬赫盤。噴流形成馬赫盤后沿著流向向下游流動(dòng),直至噴流沖擊地面,沿著徑向流出。如圖5所示,噴流沖擊水平地面后在馬赫盤的正下方形成了一個(gè)回流區(qū),回流區(qū)的幾何結(jié)構(gòu)隨時(shí)間發(fā)生變化,這說明噴流沖擊地面的過程具有顯著的非定常效應(yīng)。如圖6所示,噴流沖擊傾斜地面后的流動(dòng)速度出現(xiàn)明顯的時(shí)間演化現(xiàn)象,噴流和地面形成的回流區(qū)出現(xiàn)更明顯的變化,這說明噴流沖擊傾斜地面后具有更強(qiáng)烈的流動(dòng)不穩(wěn)定性。

(a)t=0.15 s

(b)t=0.16 s

(c)t=0.17 s

(d)t=0.18 s

(a)t=0.15 s

(b)t=0.16 s

(c)t=0.17 s

(d)t=0.18 s

對(duì)噴流沖擊傾斜地面和水平地面的流動(dòng)機(jī)理展開分析,圖7和圖8給出了對(duì)稱截面上的壓力梯度幅值分布。在兩種工況下,噴管內(nèi)部均形成了交替分布的、螺旋狀的壓力梯度分布。在喉部位置和馬赫盤位置,壓力變化最為明顯。當(dāng)噴流沖擊水平地面時(shí),噴管出口的壓力變化基本一致;當(dāng)噴流沖擊傾斜地面時(shí),噴管出口及其下游的壓力受到了擾動(dòng),交替分布的壓力梯度被破壞。對(duì)比圖7和圖8,可以發(fā)現(xiàn)在圖7中噴流沖擊水平地面的區(qū)域形成了一個(gè)穩(wěn)定的半環(huán)狀高壓力梯度區(qū)域,該區(qū)域內(nèi)部的壓力變化較大。從圖8可以看出,噴流沖擊傾斜地面后,外圈的穩(wěn)定環(huán)狀結(jié)構(gòu)被破壞,引起部分回流,進(jìn)而導(dǎo)致噴管的側(cè)向載荷發(fā)生變化。

(a)t=0.15 s

(b)t=0.16 s

(c)t=0.17 s

(d)t=0.18 s

(a)t=0.15 s

(b)t=0.16 s

(c)t=0.17 s

(d)t=0.18 s

2.2 類火星環(huán)境下地面傾斜度的影響

當(dāng)前研究還關(guān)注了在類火星環(huán)境下地面效應(yīng)對(duì)噴流的影響。已有的研究表明,較低的NPR下,流動(dòng)分離出現(xiàn)在噴管的內(nèi)部[8-9],離地距離較大時(shí),地面效應(yīng)對(duì)馬赫盤的位置影響較小。在與火星類似的低溫低壓環(huán)境下,由于噴管靜壓室的壓力不變,而外界的環(huán)境壓力約為0.669%個(gè)地球大氣壓[17-19],此時(shí)噴管的NPR迅速增大。參考文獻(xiàn)[18]的火星環(huán)境參數(shù),本小節(jié)對(duì)類火星環(huán)境下的噴流開展研究。

從圖9可以看出,在類火星環(huán)境下,噴流在噴管出口處完全膨脹,與地面碰撞形成扇形膨脹界面,這與海平面環(huán)境下觀察到的現(xiàn)象明顯不同。在不同時(shí)刻,噴流沖擊水平地面時(shí),其速度均呈對(duì)稱分布。如圖10所示,噴流沖擊傾斜地面時(shí),其速度分布不再對(duì)稱,在噴管出口的下游形成了明顯的馬赫盤結(jié)構(gòu),且馬赫盤的位置隨時(shí)間出現(xiàn)上下移動(dòng)現(xiàn)象。

(a)t=0.15 s

(b)t=0.16 s

(c)t=0.17 s

(d)t=0.18 s

(a)t=0.15 s

(b)t=0.16 s

(c)t=0.17 s

(d)t=0.18 s

在類火星環(huán)境下,噴流與不同地面發(fā)生碰撞,流場(chǎng)中的壓力梯度幅值分布顯著不同,如圖11和圖12所示。當(dāng)噴流沖擊水平地面時(shí),噴管喉部、噴流與地面相互作用區(qū)域壓力變化較為明顯。當(dāng)噴流沖擊傾斜地面時(shí),在噴管喉部、噴管出口馬赫盤的位置以及底部回流區(qū)壓力變化更加顯著。如圖11所示,噴流和地面反彈氣流之間存在較大的壓力變化,說明噴流和地面反彈氣流之間存在強(qiáng)烈的相互作用。反彈氣流在膨脹噴流的作用下沿徑向流出,在遠(yuǎn)離噴流處則向上擴(kuò)散,呈現(xiàn)對(duì)稱分布。由于地面反彈氣流的阻礙作用,膨脹氣流的范圍受到限制。如圖12所示,當(dāng)噴流沖擊傾斜地面時(shí),噴流中心線兩側(cè)的壓力梯度出現(xiàn)不對(duì)稱分布,引起回流區(qū)的氣流失穩(wěn),進(jìn)而導(dǎo)致馬赫盤的位置出現(xiàn)上下運(yùn)動(dòng)現(xiàn)象。

(a)t=0.35 s

(b)t=0.36 s

(c)t=0.37 s

(d)t=0.38 s

(a)t=0.15 s

(b)t=0.16 s

(c)t=0.17 s

(d)t=0.18 s

2.3 類火星和海平面環(huán)境下噴管壁面載荷分析

如圖13(a)所示,當(dāng)噴流沖擊水平地面時(shí),噴管內(nèi)壁的壓力幾乎不發(fā)生變化。如圖13(b)所示,當(dāng)噴流沖擊傾斜地面時(shí),噴流從壁面分離后呈現(xiàn)明顯的壓力波動(dòng)。這種壓力波動(dòng)可能導(dǎo)致噴管結(jié)構(gòu)疲勞,從而影響噴管的性能。如圖14所示,由于噴流在噴管內(nèi)部完全膨脹,在類火星環(huán)境下噴流沖擊水平地面和傾斜地面其壁面壓力載荷均維持穩(wěn)定,其載荷大小幾乎一致。

如圖15所示,位于海平面環(huán)境時(shí),噴管的內(nèi)壁面溫度呈現(xiàn)明顯非定常特性。在噴流發(fā)生流動(dòng)分離前,噴管的壁面溫度非定常特性較弱,均沿流向逐漸減低,在分離點(diǎn)處達(dá)到最低溫度。當(dāng)噴流出現(xiàn)流動(dòng)分離后,噴管的壁面溫度均上升。如圖16所示,在類火星環(huán)境下,當(dāng)噴流沖擊水平地面時(shí),噴管壁面溫度較為穩(wěn)定,幾乎不隨時(shí)間發(fā)生變化。當(dāng)噴流沖擊傾斜地面時(shí),在靠近噴管出口處,噴管內(nèi)壁的溫度出現(xiàn)波動(dòng),這與回流區(qū)的氣流失穩(wěn)有關(guān)。

(a)水平地面

(b)傾斜地面

(a)水平地面

(b)傾斜地面

(a)水平地面

(b)傾斜地面

(a)水平地面

(b)傾斜地面

2.4 類火星環(huán)境下航天器發(fā)動(dòng)機(jī)布局分析

從前文的分析中可以發(fā)現(xiàn),研究采用的噴管在火星環(huán)境下噴流完全膨脹,導(dǎo)致其噴管側(cè)向壓力載荷明顯小于地球海平面環(huán)境,這從側(cè)面說明了高空發(fā)動(dòng)機(jī)在海平面試車,一般采用短噴管狀態(tài)[23]。對(duì)于用于類似于SpaceX公司星艦構(gòu)型的航天器來說,其既具有垂直返回地面的能力,又具有星際航行的需求,因此從總體設(shè)計(jì)考慮,其底部發(fā)動(dòng)機(jī)布局應(yīng)兼顧二者需求。其采用大噴管的發(fā)動(dòng)機(jī)用于提供火箭上升段的高比沖和大推力以及在月球或者火星表面進(jìn)行航天器垂直著陸;其采用短噴管的發(fā)動(dòng)機(jī)也可用于提供火箭上升段推力,但更為重要的是支持在地球表面返回。

因此,如果在星際往返過程中,航天器采用垂直起降的方式,則需要根據(jù)目標(biāo)行星的大氣環(huán)境開展有針對(duì)性的發(fā)動(dòng)機(jī)噴管設(shè)計(jì)。

3 結(jié)論

本文利用SAS方法研究了火箭位于海平面環(huán)境和火星環(huán)境下,發(fā)動(dòng)機(jī)噴流沖擊水平地面和傾斜地面的流場(chǎng)結(jié)構(gòu),分析了在不同工況下噴管內(nèi)壁的側(cè)向壓力載荷和溫度載荷。主要結(jié)論如下:

1)噴流沖擊傾斜地面時(shí)具有更強(qiáng)烈的流動(dòng)不穩(wěn)定性,噴流出現(xiàn)明顯的時(shí)間演化現(xiàn)象。

2)類火星環(huán)境下,噴流沖擊傾斜地面導(dǎo)致馬赫盤上下運(yùn)動(dòng)。噴流沖擊水平地面時(shí),噴流和地面反彈流相互作用,阻礙了噴流馬赫盤的形成。

3)海平面環(huán)境下,在當(dāng)前大面積比噴管內(nèi)部出現(xiàn)流動(dòng)分離,噴管內(nèi)壁的壓力先減小,在分離點(diǎn)處壓力達(dá)到最低點(diǎn),然后壓力增大至環(huán)境壓力;類火星環(huán)境下,噴管內(nèi)部流動(dòng)不發(fā)生分離,噴管內(nèi)壁的壓力沿著流向逐漸減小。

4)噴管的壁面溫度在類火星環(huán)境下沿著流向逐漸減低,相對(duì)水平地面,噴流沖擊傾斜地面時(shí)溫度較高。在海平面環(huán)境下,壁面溫度和壁面壓力呈現(xiàn)相同的趨勢(shì)。

本文對(duì)噴管在不同環(huán)境下的流動(dòng)結(jié)構(gòu)及噴管壁面載荷進(jìn)行了研究,研究結(jié)果具有一定的指導(dǎo)性。但是受限于當(dāng)前的試驗(yàn)數(shù)據(jù),僅對(duì)低NPR工況下進(jìn)行了方法驗(yàn)證。為了對(duì)該類問題進(jìn)行更深刻的研究,這要求在之后的研究中獲得更多的試驗(yàn)數(shù)據(jù)。同時(shí),也亟須發(fā)展和改進(jìn)適用于更高馬赫數(shù)的算法。由于大渦模擬和直接數(shù)值模擬等算法計(jì)算量巨大,為了兼顧計(jì)算效率和精確度,SAS和DDES等類大渦模擬方法為研究此類問題提供了一個(gè)可行的途徑。

猜你喜歡
環(huán)境研究
FMS與YBT相關(guān)性的實(shí)證研究
長(zhǎng)期鍛煉創(chuàng)造體內(nèi)抑癌環(huán)境
2020年國(guó)內(nèi)翻譯研究述評(píng)
遼代千人邑研究述論
一種用于自主學(xué)習(xí)的虛擬仿真環(huán)境
孕期遠(yuǎn)離容易致畸的環(huán)境
視錯(cuò)覺在平面設(shè)計(jì)中的應(yīng)用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
不能改變環(huán)境,那就改變心境
EMA伺服控制系統(tǒng)研究
環(huán)境
主站蜘蛛池模板: 久久香蕉国产线看观看式| 国产精品太粉嫩高中在线观看| 国产中文一区a级毛片视频| 欧美一区二区精品久久久| 在线日韩日本国产亚洲| 国产精品视频第一专区| 日韩无码真实干出血视频| 中文字幕永久在线看| 欧美高清国产| 日本亚洲欧美在线| 制服丝袜一区| 尤物国产在线| 第九色区aⅴ天堂久久香| 嫩草影院在线观看精品视频| 99在线观看免费视频| 国产欧美日韩另类| 亚洲动漫h| 中文字幕有乳无码| 亚洲AⅤ波多系列中文字幕 | 国产成人91精品免费网址在线| 一级香蕉人体视频| 欧美狠狠干| 久久亚洲国产最新网站| 人妻无码中文字幕第一区| 国产粉嫩粉嫩的18在线播放91| 一级毛片在线免费视频| 亚洲日韩AV无码一区二区三区人 | 国产美女精品一区二区| 国产精品三区四区| 国产精品专区第1页| 亚洲五月激情网| 99久久无色码中文字幕| 成年看免费观看视频拍拍| 久久综合久久鬼| 美女扒开下面流白浆在线试听| 四虎国产永久在线观看| 东京热av无码电影一区二区| 国产精品无码影视久久久久久久| 国产精品无码久久久久AV| 欧美激情视频一区| 青青草91视频| 欧美精品高清| 欧美一区中文字幕| www.精品国产| 强奷白丝美女在线观看| 一区二区日韩国产精久久| 亚洲综合色在线| 国产一区二区福利| 国产美女无遮挡免费视频| 国产手机在线ΑⅤ片无码观看| 亚洲精品自产拍在线观看APP| 91精品国产麻豆国产自产在线| 亚洲黄色高清| 亚洲日韩图片专区第1页| 91成人免费观看| 欧美精品亚洲精品日韩专| 91欧美在线| 无码日韩视频| 啪啪啪亚洲无码| 色婷婷成人网| 无码国产偷倩在线播放老年人| 国产成人午夜福利免费无码r| 欧美综合中文字幕久久| 亚洲最黄视频| 国产精品美人久久久久久AV| 久久久精品无码一二三区| 国产精品视频公开费视频| 日韩 欧美 小说 综合网 另类 | 天天躁狠狠躁| 中文毛片无遮挡播放免费| 国产欧美亚洲精品第3页在线| 国产欧美成人不卡视频| 国产成人亚洲毛片| 色婷婷在线影院| 91精品国产自产在线观看| 亚洲第一福利视频导航| 91在线无码精品秘九色APP| 香蕉色综合| 国产乱肥老妇精品视频| 一区二区三区成人| 精品成人一区二区| 91精品国产自产在线老师啪l|