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

多機(jī)并聯(lián)火箭羽流流場及其底部熱環(huán)境分析

2022-01-20 07:16:18丁兆波楊進(jìn)慧
宇航學(xué)報 2021年11期
關(guān)鍵詞:發(fā)動機(jī)環(huán)境

楊 巖,田 原,丁兆波,楊進(jìn)慧

(北京航天動力研究所,北京 100076)

0 引 言

多機(jī)并聯(lián)火箭飛行時箭體底部流場分布極其復(fù)雜,且隨飛行高度變化而變化。隨之帶來的箭體底部熱環(huán)境極為惡劣,對箭體底部熱環(huán)境預(yù)估的準(zhǔn)確性將直接影響火箭底部防熱結(jié)構(gòu)設(shè)計的正確性。如果防熱設(shè)計考慮不足,飛行過程中艙外結(jié)構(gòu)可能發(fā)生燒蝕,影響飛行成敗;如果設(shè)計過于保守,又會使得結(jié)構(gòu)偏重,不利于火箭運(yùn)載性能提升[1-2]。

目前國外已有大量學(xué)者對多機(jī)并聯(lián)火箭羽流及其底部熱環(huán)境進(jìn)行了研究。早在1961年,Norman等[3]設(shè)計了小型四機(jī)并聯(lián)火箭并進(jìn)行了風(fēng)洞試驗,得到了高空飛行時箭體底部流場特征及其熱環(huán)境的影響因素。Otis等[4]對Titan III固體助推發(fā)動機(jī)羽流熱輻射進(jìn)行了試驗測試,結(jié)果表明剛起飛1 s內(nèi)輻射熱流是總熱流的主要來源,隨著飛行高度上升,輻射熱流減少,對流熱流升高。Terry等[5]對航天飛機(jī)一二級飛行時底部熱環(huán)境進(jìn)行了詳細(xì)分析,設(shè)計值與飛行實測值吻合較好,并指出航天飛機(jī)底部熱環(huán)境主要由發(fā)動機(jī)羽流輻射、自由氣流冷卻和回流加熱等構(gòu)成。上述分析大多采用理論計算或風(fēng)洞試驗對火箭底部熱環(huán)境進(jìn)行分析,隨著計算機(jī)技術(shù)的發(fā)展,越來越多的學(xué)者采用計算流體力學(xué)對火箭底部熱環(huán)境進(jìn)行研究,同時也得到了更為詳細(xì)的羽流流場結(jié)構(gòu)分布。如:Philippe等[6]對Ariane 5火箭飛行流場進(jìn)行了數(shù)值模擬,并開展了縮比模型風(fēng)洞試驗,得出影響對流熱流的主要因素是飛行高度、噴管數(shù)量以及噴管排列方式。Hideyo等[7]對H-IIA火箭羽流熱環(huán)境開展了數(shù)值計算,結(jié)果表明不同發(fā)動機(jī)羽流間的相互作用及其回流是火箭底部加熱的主要來源。Aaron等[8]則利用縮比模型風(fēng)洞試驗對空間發(fā)射系統(tǒng)(Space Launch System,SLS)飛行熱環(huán)境進(jìn)行了研究,并與數(shù)值計算結(jié)果進(jìn)行了對比分析。

中國常規(guī)系列運(yùn)載火箭研制初期,受限于研究基礎(chǔ)和計算機(jī)水平,對火箭羽流流場及其熱環(huán)境方面研究相對較少,大多通過飛行試驗進(jìn)行考核。近年來,伴隨著新一代運(yùn)載火箭的研制,國內(nèi)學(xué)者通過計算流體力學(xué)等手段對火箭底部羽流流場及其熱環(huán)境問題開展了一系列研究分析。羅天培等[9]對某型運(yùn)載火箭一子級動力系統(tǒng)試驗尾流輻射場進(jìn)行了數(shù)值模擬,計算結(jié)果表明箭體底部的換熱方式以輻射換熱為主。范瑞祥等[10]針對某新型運(yùn)載火箭二級尾艙防熱問題,采用CFD/DSMC相結(jié)合的辦法對四機(jī)并聯(lián)的高空羽流及其回流流場進(jìn)行了仿真分析,提出了可靠有效的防熱措施,改善了火箭二級尾艙的熱環(huán)境。李國良等[11]基于某單發(fā)動機(jī)飛行器,研究了發(fā)動機(jī)噴流對底部流動的影響。周志壇等[12]對某單發(fā)動機(jī)火箭底部熱環(huán)境進(jìn)行了研究,研究結(jié)果表明火箭底部輻射熱流在起飛階段最大,隨飛行高度上升,輻射熱流逐漸降低,對流熱流則呈現(xiàn)先升高后降低的趨勢。楊學(xué)軍等[13]采用理論預(yù)測和數(shù)值計算方法,對某固體火箭尾艙熱環(huán)境進(jìn)行了研究,結(jié)果表明尾艙熱環(huán)境存在天地差異,設(shè)計時應(yīng)充分考慮真實飛行狀態(tài)下的熱環(huán)境。

上述研究大多基于理論分析或地面試驗結(jié)果,且多針對單個火箭發(fā)動機(jī)噴流進(jìn)行分析,針對多機(jī)并聯(lián)火箭真實飛行過程的熱環(huán)境變化情況研究相對較少。本文針對某新型多機(jī)并聯(lián)火箭實際飛行工況,采用數(shù)值模擬方法對其羽流流場及箭體底部熱環(huán)境進(jìn)行計算,并與熱環(huán)境實測值進(jìn)行了對比分析,計算結(jié)果對后續(xù)熱防護(hù)結(jié)構(gòu)改進(jìn)優(yōu)化設(shè)計具有一定的參考意義。

1 計算方法

1.1 物理模型

某新一代運(yùn)載火箭由芯一級捆綁4臺助推器構(gòu)成,其芯一級及助推發(fā)動機(jī)布局如圖1所示。其中芯一級由兩臺氫氧發(fā)動機(jī)提供推力,四個助推由兩臺液氧煤油發(fā)動機(jī)提供推力。為研究實際飛行過程中發(fā)動機(jī)羽流對底部熱環(huán)境的影響,在芯一級防熱底上布置有3個熱流傳感器Q1~Q3,測量精度±3%,具體位置如圖1所示。

圖1 某新一代運(yùn)載火箭結(jié)構(gòu)示意圖Fig.1 Schematic configuration of a new generation rocket

1.2 數(shù)值計算

由于計算的流場結(jié)構(gòu)復(fù)雜,流場中包含了激波、膨脹波、激波干擾、邊界層和尾跡區(qū)等復(fù)雜流場結(jié)構(gòu),本文采用成熟的商業(yè)CFD計算軟件FLUENT進(jìn)行計算。湍流模型采用的k-ε雙方程模型,計算中考慮化學(xué)反應(yīng)及輻射散熱。采用壓力基求解器,利用SIMPLEC算法,梯度項采用最小二乘格式,壓力項采用標(biāo)準(zhǔn)格式,其余各項均采用一階迎風(fēng)格式進(jìn)行離散求解。

1)計算網(wǎng)格

為精確捕捉多個發(fā)動機(jī)羽流相互干擾的流場結(jié)構(gòu),采用商業(yè)網(wǎng)格劃分軟件ICEM進(jìn)行了全結(jié)構(gòu)化網(wǎng)格劃分。考慮計算能力,采用對稱面進(jìn)行簡化,計算域僅為實際的1/4。為準(zhǔn)確描述芯級發(fā)動機(jī)底部流場特征及其熱環(huán)境,對芯級發(fā)動機(jī)底部以及不同發(fā)動機(jī)羽流交匯處進(jìn)行了網(wǎng)格加密處理,模型網(wǎng)格數(shù)量在1100萬左右,網(wǎng)格劃分結(jié)果如圖2所示。

圖2 網(wǎng)格示意圖Fig.2 Mesh generation of the model

2)湍流模型

計算采用Standardk-ε模型,該模型為半經(jīng)驗?zāi)P?假設(shè)整個流場為湍流,忽略分子黏性的影響,在較大的湍流流動范圍內(nèi)表現(xiàn)出較好的穩(wěn)定性。

Gb-ρε-YM+Sk

(1)

(2)

式中:t為時間;x,y為坐標(biāo)軸;ρ為流體密度;μ為流體層流黏度;μt為流體湍流黏度;ui為速度在i方向上分量;Gk是由平均速度梯度而產(chǎn)生的湍流動能;Gb是由浮力而產(chǎn)生的湍流動能;YM是由在可壓縮湍流過渡到全部擴(kuò)散速率而引起的波動擴(kuò)張;σk和σε分別為k和ε的湍流普朗特數(shù);C1ε,C2ε和C3ε均為常數(shù);Sk和Sε為自定義源相。

3)化學(xué)反應(yīng)

為研究化學(xué)反應(yīng)對箭體底部熱環(huán)境的影響,分別對有無化學(xué)反應(yīng)進(jìn)行了計算,將計算結(jié)果與實測值進(jìn)行對比分析。考慮現(xiàn)有計算能力,本文采用FLUENT自帶的一步反應(yīng)。反應(yīng)機(jī)理如表1所示,均采用一步總包反應(yīng)機(jī)理,有限速率模型進(jìn)行計算。

表1 化學(xué)反應(yīng)機(jī)理Table 1 Chemical reaction mechanism

4)輻射模型

液體火箭飛行時底部熱輻射[14-15]主要包括兩個方面:

① 發(fā)動機(jī)羽流輻射作用:在發(fā)動機(jī)羽流中,存在著大量的H2O和CO2,對箭體底部艙外結(jié)構(gòu)有著較強(qiáng)的輻射作用。

② 各發(fā)動機(jī)噴管外壁面的輻射作用:液氧煤油發(fā)動機(jī)為再生冷卻噴管,外壁溫較低,其輻射可以忽略。氫氧發(fā)動機(jī)外壁面溫度變化較大,噴管出口位置溫度高達(dá)900 K,對渦輪排氣管有著一定的輻射作用,其外壁溫通過UDF(User defined dunctions)給定。

在FLUENT的輻射模型中,DO(Discrete ordinates radiation model)模型既能求解S2S(Surface to surface)的無介質(zhì)封閉區(qū)域問題,也能求解介質(zhì)參與的輻射問題。綜合考慮計算時間、計算精度以及各輻射模型的應(yīng)用范圍,本次計算過程中采用DO輻射模型對羽流的熱輻射情況進(jìn)行仿真,其中吸收系數(shù)由H2O和CO2吸收系數(shù)通過質(zhì)量平均獲得(即FLUENT中WSGGM模型)[16]。

計算過程中,先不考慮輻射模型,對羽流流場及箭體底部熱環(huán)境進(jìn)行計算,此時計算得到的熱流為對流熱流;計算收斂后,開啟輻射模型重新計算至收斂,此時得到的即為對流熱和輻射熱的總和。

5)邊界條件

計算時邊界條件設(shè)置如下:

① 噴管入口條件:采用壓力入口條件,芯級發(fā)動機(jī)入口壓力10.1 MPa,助推發(fā)動機(jī)入口壓力18 MPa,溫度及組分根據(jù)熱力計算結(jié)果給定。

② 出口條件:采用壓力出口條件,壓力值為所在高度大氣環(huán)境值。

③ 壁面條件:助推發(fā)動機(jī)為再生冷卻噴管,其外壁面輻射可以忽略,給定常溫邊界條件;芯級發(fā)動機(jī)噴管外壁面溫度參考一維傳熱計算結(jié)果,外壁溫通過UDF給定;箭體底部設(shè)為恒溫邊界條件,模擬熱流傳感器冷端,具體溫度根據(jù)所在高度大氣參數(shù)給定,實際該處溫度可能隨飛行變高,因此熱流計算結(jié)果可能存在一定偏差。

④ 遠(yuǎn)場條件:由于箭體頭部激波距離芯一級底部較遠(yuǎn),經(jīng)分析認(rèn)為來流不會對發(fā)動機(jī)羽流及其回流造成影響,計算中不考慮來流影響,壓力和溫度都設(shè)為所在高度大氣環(huán)境參數(shù)。

火箭飛行中不同高度下的大氣環(huán)境參數(shù)見表2。

表2 大氣環(huán)境參數(shù)Table 2 Atmospheric parameters

2 分析與討論

2.1 羽流流場分布

從芯級發(fā)動機(jī)截面流場參數(shù)分布(如圖3所示)可以看出,低空飛行時,各發(fā)動機(jī)流場互不干擾;隨著飛行高度增加,羽流逐漸擴(kuò)張,大約在88 s(16 km)時,各發(fā)動機(jī)羽流開始出現(xiàn)干擾。由于液氧煤油發(fā)動機(jī)的推力室室壓和噴管出口壓力均高于氫氧發(fā)動機(jī),芯級氫氧發(fā)動機(jī)流場受到助推液氧煤油發(fā)動機(jī)羽流的擠壓;飛行至約107 s(25 km)時,芯級底部開始出現(xiàn)回流,從而對箭體底部產(chǎn)生對流加熱。

從距離發(fā)動機(jī)噴口2.5 m位置處的流場參數(shù)分布(如圖3所示)同樣可以看出,低空飛行時,各發(fā)動機(jī)羽流互不干擾;隨飛行高度增加,環(huán)境壓力不斷減小,發(fā)動機(jī)羽流不斷擴(kuò)張,不同發(fā)動機(jī)羽流間出現(xiàn)嚴(yán)重干擾,在兩臺芯級發(fā)動機(jī)和助推發(fā)動機(jī)之間形成較明顯的高溫區(qū)和高速區(qū),從而在箭體底部形成極其復(fù)雜的流場環(huán)境,有可能對艙外結(jié)構(gòu)熱防護(hù)造成一定威脅。

圖3 不同飛行高度下的流場分布Fig.3 Flow field distribution at different fight altitudes

為進(jìn)一步分析箭體底部流場參數(shù)對熱環(huán)境的影響,將芯級底部兩臺氫氧發(fā)動機(jī)間燃?xì)鉁囟取毫Α⒚芏燃案鹘M分濃度參數(shù)進(jìn)行統(tǒng)計結(jié)果如圖4和圖5所示。從圖中可以看出:火箭起飛后燃?xì)鉁囟扔行》陆第厔?大約在100 s之后出現(xiàn)明顯上升。這是因為低空飛行時在羽流引射作用下,進(jìn)入箭體底部的空氣對相關(guān)結(jié)構(gòu)具有一定的冷卻作用,當(dāng)箭體底部出現(xiàn)回流后,燃?xì)鉁囟让黠@上升;而燃?xì)鈮毫兔芏入S飛行高度增加呈現(xiàn)不斷降低的趨勢;箭體底部水蒸氣和二氧化碳組分濃度在低空時基本不變,隨飛行高度增加,二者組分濃度也出現(xiàn)明顯上升,也進(jìn)一步說明約100 s之后,箭體底部出現(xiàn)了明顯的燃?xì)饣亓鳌?/p>

圖4 箭體底部溫度及壓力變化趨勢Fig.4 Temperature and pressure trends at bottom of the rocket

圖5 箭體底部組分濃度變化情況Fig.5 Component concentration trends at bottom of the rocket

2.2 熱環(huán)境分析

為分析箭體底部熱環(huán)境,對不同高度下芯級底部熱流計算結(jié)果進(jìn)行統(tǒng)計,并與飛行實測結(jié)果進(jìn)行對比分析,如圖6所示。從圖6可以看出,計算結(jié)果與實測結(jié)果趨勢基本吻合,火箭剛起飛時,箭體底部受發(fā)動機(jī)羽流的熱輻射影響最大,最高約550 kW/m2。此后隨飛行高度不斷增加,環(huán)境壓力降低,燃?xì)獠粩鄶U(kuò)張,輻射能力減弱,箭體底部所受輻射熱開始減小;約107 s左右,箭體底部出現(xiàn)回流后,由此產(chǎn)生的對流熱開始增加,直到助推分離時達(dá)到最大值。

圖6 熱流仿真值與飛行實測值對比Fig.6 Heat flux in numerical calculation compared with test data in flight

火箭低空飛行時,箭體外部的來流會對箭體底部發(fā)動機(jī)形成一定的冷卻作用,而高空飛行時,箭體底部形成明顯回流,對底部結(jié)構(gòu)件形成對流加熱。因本計算中未考慮來流的影響,低空飛行時熱流計算結(jié)果略高于實測值,高空飛行時對底部回流預(yù)估可能存在不足,導(dǎo)致計算熱流低于實測值。

為分析箭體底部熱流組成,將計算得到的總熱流和對流熱流進(jìn)行統(tǒng)計如圖7所示。從圖7可以看出,約107 s之后箭體底部對流熱流開始增加,到約145 s時對流熱流達(dá)到最大值,此后基本保持不變。此外隨飛行高度增加,羽流不斷擴(kuò)張并產(chǎn)生回流,箭體底部羽流的輻射角系數(shù)逐漸增大,因此100 s之后輻射熱也不斷增加,總熱流大約在助推分離前再次達(dá)到峰值。

圖7 箭體底部熱流組成Fig.7 Heat flux composed in bottom of the rocket

根據(jù)氣體輻射理論[17],當(dāng)同時存在水蒸氣和二氧化碳兩種組分時,氣體輻射率由下式計算:即氣體的輻射能力受氣體溫度、組分和沿途吸收性氣體分子數(shù)目有關(guān),沿途氣體分子數(shù)顯然與氣體分壓P和平均射線程長s的乘積Ps成正比。

εg=f(Tg,Ps)

(3)

式中:εg為氣體發(fā)射率,Tg為氣體溫度。

根據(jù)上述流場計算結(jié)果,火箭起飛后,隨著飛行高度增加,燃?xì)鈮毫懊芏炔粩嘟档?輻射能力減弱,輻射熱流也不斷減小;約100 s之后,箭體底部出現(xiàn)明顯回流,水蒸氣和二氧化碳濃度明顯增加,一方面燃?xì)饣亓鲗?dǎo)致對流熱流增加,同時由于組分濃度及羽流輻射角系數(shù)的變化,燃?xì)廨椛淠芰υ俅卧鰪?qiáng),進(jìn)而導(dǎo)致箭體底部熱流再次增加,直到助推分離時達(dá)到最大值。

如不考慮化學(xué)反應(yīng)即羽流與空氣的補(bǔ)燃,計算得到起飛時刻箭體底部熱流值約為126 kW/m2,與實測值偏差過大;從助推發(fā)動機(jī)溫度場分布也可以明顯看出各發(fā)動機(jī)噴流與周圍的空氣發(fā)生了補(bǔ)燃,導(dǎo)致局部溫度升高。綜合分析認(rèn)為發(fā)動機(jī)羽流與空氣發(fā)生了補(bǔ)燃,且對箭體底部熱輻射影響較大。

2.3 艙外結(jié)構(gòu)熱防護(hù)

為研究飛行中發(fā)動機(jī)艙外結(jié)構(gòu)件熱防護(hù)的可靠性,基于上述分析得到的熱環(huán)境,對某關(guān)鍵結(jié)構(gòu)開展了瞬態(tài)熱分析。計算相關(guān)邊界條件如下:

1)結(jié)構(gòu)初始溫度設(shè)為295 K,與燃?xì)饨佑|面設(shè)為恒溫邊界條件,溫度1100 K。

2)外壁面通過輻射方式向外散熱。

3)根據(jù)地面試驗與仿真校算結(jié)果,結(jié)構(gòu)發(fā)射率取為0.2。

4)結(jié)構(gòu)所處熱環(huán)境根據(jù)上述計算中得到的熱流給定。

為準(zhǔn)確評估發(fā)動機(jī)羽流對結(jié)構(gòu)熱防護(hù)的影響,在兩次飛行任務(wù)中對某關(guān)鍵結(jié)構(gòu)壁面溫度進(jìn)行了測量,與計算結(jié)果對比情況如圖8所示。從圖8可以看出,計算結(jié)果與實測溫度曲線基本吻合。受實際飛行參數(shù)、計算所用發(fā)射率等因素綜合影響,計算所得最高溫度1109 K與實測最高溫度1189 K存在一定誤差。實測最高溫度下所用材料仍具有一定的結(jié)構(gòu)強(qiáng)度,可以可靠工作。

圖8 助推發(fā)動機(jī)溫度場分布Fig.8 Temperature distribution of the rocket booster engine

相比飛行實測結(jié)果,計算得到的溫度曲線在低空飛行時上升較快,結(jié)合羽流流場分析,該差異主要是由于計算中未考慮來流對羽流流場及其熱環(huán)境的影響。而低空飛行時,在羽流引射作用下,來流對箭體底部具有一定的冷卻作用。

圖9 某結(jié)構(gòu)表面溫度計算值與飛行實測值對比Fig.9 Comparison of measured and calculated temperature values

3 結(jié) 論

通過對某火箭助推及一級發(fā)動機(jī)的羽流流場及底部熱環(huán)境進(jìn)行計算分析,可以得出以下結(jié)論:

1)低空飛行時各發(fā)動機(jī)羽流互不干擾,隨飛行高度增加,羽流逐漸擴(kuò)張,大約飛行至約16 km時,各發(fā)動機(jī)羽流開始互相干擾,且氫氧發(fā)動機(jī)的羽流明顯受到液氧煤油發(fā)動機(jī)羽流的擠壓;約25 km時,箭體底部開始出現(xiàn)回流,隨之產(chǎn)生對流熱,在此之前,箭體底部僅受羽流的輻射熱影響。

2)發(fā)動機(jī)羽流與周圍空氣發(fā)生了補(bǔ)燃,起飛時刻箭體底部熱環(huán)境最為惡劣,約550 kW/m2,助推分離時刻,對流熱占總熱流比例達(dá)到最大,約為30%,但對流熱流最大值僅約20 kW/m2,計算結(jié)果與熱流實測趨勢基本吻合。發(fā)動機(jī)艙外相關(guān)結(jié)構(gòu)件進(jìn)行熱防護(hù)時,應(yīng)主要采用輻射熱防護(hù)方法進(jìn)行熱防護(hù)。

3)發(fā)動機(jī)艙外結(jié)構(gòu)件瞬態(tài)熱分析結(jié)果與實測值基本吻合,相關(guān)結(jié)構(gòu)熱防護(hù)條件滿足使用需求,可以可靠工作。

猜你喜歡
發(fā)動機(jī)環(huán)境
長期鍛煉創(chuàng)造體內(nèi)抑癌環(huán)境
一種用于自主學(xué)習(xí)的虛擬仿真環(huán)境
元征X-431實測:奔馳發(fā)動機(jī)編程
2015款寶馬525Li行駛中發(fā)動機(jī)熄火
孕期遠(yuǎn)離容易致畸的環(huán)境
不能改變環(huán)境,那就改變心境
環(huán)境
孕期遠(yuǎn)離容易致畸的環(huán)境
新一代MTU2000發(fā)動機(jī)系列
發(fā)動機(jī)的怠速停止技術(shù)i-stop
主站蜘蛛池模板: 91久久天天躁狠狠躁夜夜| 国内熟女少妇一线天| 国产精品妖精视频| 91av成人日本不卡三区| 亚洲欧美国产视频| 国内精品小视频福利网址| 热伊人99re久久精品最新地| 91破解版在线亚洲| 亚洲最大综合网| 99热亚洲精品6码| 白浆视频在线观看| 自拍亚洲欧美精品| 成人字幕网视频在线观看| 欧美亚洲国产日韩电影在线| 成人在线观看不卡| 欧美激情首页| 四虎永久在线精品国产免费| 国产在线观看一区精品| 国产视频你懂得| 国产精品三级av及在线观看| 狂欢视频在线观看不卡| 欧美一区二区丝袜高跟鞋| 国产在线欧美| 亚洲欧美精品日韩欧美| 一本久道久综合久久鬼色| 国产18在线| 91免费观看视频| av午夜福利一片免费看| 国产99久久亚洲综合精品西瓜tv| 国内老司机精品视频在线播出| P尤物久久99国产综合精品| 日本高清在线看免费观看| 久久精品中文字幕免费| 国产成a人片在线播放| 久久久久夜色精品波多野结衣| 欧美在线免费| 国模粉嫩小泬视频在线观看 | 亚洲国产日韩在线成人蜜芽| 青青草国产免费国产| 美美女高清毛片视频免费观看| 国产成人精品在线1区| 国产 日韩 欧美 第二页| 国产自产视频一区二区三区| 一级爱做片免费观看久久| 亚洲欧洲日本在线| 成人免费一级片| 国产午夜福利在线小视频| 成人精品区| 日韩美女福利视频| 国产色婷婷视频在线观看| 热久久这里是精品6免费观看| 啪啪永久免费av| 色成人亚洲| 国产呦视频免费视频在线观看| 国产性生大片免费观看性欧美| 国产午夜一级毛片| 在线观看的黄网| 婷婷久久综合九色综合88| 国产成人综合久久精品下载| 国产高清不卡| 国产SUV精品一区二区6| 欧美在线视频a| 第九色区aⅴ天堂久久香| 有专无码视频| 自拍亚洲欧美精品| 亚洲 成人国产| 在线观看欧美国产| 噜噜噜久久| 亚洲欧美日韩动漫| 亚洲Av综合日韩精品久久久| 国产九九精品视频| 亚洲国产精品国自产拍A| 国产性生交xxxxx免费| 亚洲精品国产首次亮相| 国产原创自拍不卡第一页| 999在线免费视频| 好吊色妇女免费视频免费| 亚洲视频免费在线看| 麻豆国产在线观看一区二区 | 成人福利在线视频| 成人无码区免费视频网站蜜臀| 久久久久久高潮白浆|