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

肼類單組元發(fā)動機冷起動過程溫度敏感性研究

2021-12-23 11:19:00雷凡培梁樹強肖明杰
關(guān)鍵詞:催化劑發(fā)動機模型

冀 鵬,雷凡培,梁樹強,肖明杰

(1. 西安航天動力研究所液體火箭發(fā)動機重點實驗室,西安,710100;2. 中國船舶集團有限公司,北京,100044; 3. 航天推進技術(shù)研究院,西安,710100)

0 引 言

肼類單組元發(fā)動機是運載火箭和各類航天器姿態(tài)控制、軌道控制使用最廣泛的推進系統(tǒng),其性能以及穩(wěn)定性直接影響著運載火箭和航天器的控制精度和可靠性[1,2]。近年來隨著航天器飛行任務(wù)的增多,對單組元發(fā)動機的環(huán)境適應(yīng)性愈加嚴(yán)苛,單組元肼發(fā)動機在冷起動時,時常出現(xiàn)起動響應(yīng)延遲長、起動壓力峰大的現(xiàn)象。一方面較長的起動響應(yīng)延遲難以滿足總體對推進系統(tǒng)的響應(yīng)要求,另一方面較高的起動壓力峰也降低了發(fā)動機的結(jié)構(gòu)可靠性,可能導(dǎo)致發(fā)動機結(jié)構(gòu)破壞。亟待對肼類單組元催化分解發(fā)動機冷起動工作過程開展深入的研究工作。

采用催化分解產(chǎn)生熱燃?xì)獾碾骂悊谓M元發(fā)動機,推進劑進入推力室后涉及在催化床這種多孔介質(zhì)中的多相多尺度流動、均相熱分解和異相催化分解等復(fù)雜的物理化學(xué)過程,目前對該過程的研究主要依靠試驗。Hearn[3]在研究脈沖工作對催化床壽命影響時發(fā)現(xiàn),在脈沖工作模式下冷起動壓力峰不會發(fā)生在較短的脈沖過程中。Kagawa等[4,5]采用可視化技術(shù)研究催化床內(nèi)液態(tài)肼流動狀態(tài),發(fā)現(xiàn)相較于溫啟動和熱啟動,冷起動時液態(tài)肼的滲透面積更大、滲透距離更遠(yuǎn)。此外,火星著陸飛船[6]、阿里安5號運載火箭[7]、NASA空間飛行器[8,9]等所采用單組元發(fā)動機研制過程中,也主要采用試驗手段來考核和鑒定發(fā)動機的冷起動能力并評判其對發(fā)動機性能的影響,發(fā)現(xiàn)過低的起動邊界溫度會導(dǎo)致肼推力器的使用壽命大大降低。單組元發(fā)動機的冷起動仿真研究方面目前尚未有公開文獻(xiàn)發(fā)表。

一般認(rèn)為溫度較低時化學(xué)反應(yīng)速率會減慢,但在單組元肼類發(fā)動機冷起動過程中,會出現(xiàn)數(shù)倍于溫啟動或熱啟動的響應(yīng)延遲,并且延遲過程中室壓維持不變,然后突然產(chǎn)生明顯的壓力峰,該現(xiàn)象與低溫導(dǎo)致化學(xué)反應(yīng)變慢機理并不契合。因此考慮該現(xiàn)象與微觀催化分解過程有關(guān)。肼分解用催化劑并不是簡單的球形或柱形顆粒結(jié)構(gòu),而是由復(fù)雜的浸漬、焙燒等制備工藝得到的具有眾多微孔通道的層狀骨架結(jié)構(gòu)[10]。Sanglovanni等[11]通過試驗觀察到在低溫情況下,液態(tài)肼滴定到催化劑顆粒后迅速浸濕了疏松的催化劑顆粒,在試驗后發(fā)現(xiàn)催化劑顆粒發(fā)生了明顯的破碎。趙許群[12]在對低溫下催化劑的失活機理進行分析時提到,在低溫下肼會在催化劑微孔內(nèi)短時劇烈分解產(chǎn)生較大內(nèi)應(yīng)力導(dǎo)致催化劑顆粒破碎。

綜上,現(xiàn)在仍然缺乏對肼類單組元發(fā)動機冷起動過程的完整過程認(rèn)識和可靠的數(shù)值仿真方法,導(dǎo)致在發(fā)動機研制過程中主要依賴試驗來預(yù)判和評價發(fā)動機的低溫起動特性。因此結(jié)合催化劑微孔通道內(nèi)流動-反應(yīng)和經(jīng)典燃燒時滯理論,開展肼類單組元發(fā)動機冷起動過程數(shù)值仿真分析,研究初始溫度對冷起動壓力與延遲時間的影響,揭示發(fā)動機冷起動延遲長、壓力峰高的主導(dǎo)因素。

1 發(fā)動機典型組件的動力學(xué)模型

對于恒壓擠壓式的單組元姿控發(fā)動機,工作過程中貯箱壓力基本保持不變,因此假設(shè)貯箱內(nèi)氣體壓力恒定,主要建立起液路推進劑供應(yīng)系統(tǒng)以及催化床的動力學(xué)模型。

1.1 管路模型

在建立液體管路模型時不考慮管路的變形,認(rèn)為管道為剛性的;管路內(nèi)液體是一維的,管道橫截面上的速度分布是均勻的。本文主要關(guān)注的發(fā)動機起動過程屬于低頻動力學(xué)問題,因此采用分段集中參數(shù)的管路模型即可滿足精度要求。根據(jù)經(jīng)驗,一種比較保守的方法是將分段管長取波長的4%以內(nèi),即管路分段長度應(yīng)滿足:Lmax≤0.04a/fmax。此外,毛細(xì)管、閥門流道等在一定程度上也可以等效為液體管路。

對于特定的一段長度、直徑、截面積和體積分別為L、d、A、V的管路,考慮慣性、粘性時管路瞬變流:

式中qm,p1,p2分別為管路流量,管路入口和出口處壓力;j為管內(nèi)液體的慣性系數(shù),ALj/=;ξ為總流阻系數(shù),由管路沿程損失和局部損失組成,取決于流體管路的形狀和尺寸、管壁粗糙度及流體在管路的流動狀態(tài)等,其中流體流動狀態(tài)常采用雷諾數(shù)來表征。

考慮液體壓縮性的方程:

式中z為管路的流容系數(shù),z=V/a2,表征管路流體的壓縮性。

1.2 催化床模型

本文采用反映平均效應(yīng)的零維模型。在建立催化床模型時作如下假設(shè):分解產(chǎn)物為理想氣體;不考慮氣體的不均勻性和波動過程;不考慮催化床內(nèi)催化劑活性點處失活的影響。

1.2.1 冷起動延遲模型

在單組元發(fā)動機冷起動時,考慮到此時催化床和推進劑的溫度都比較低,液體推進劑噴入催化床內(nèi)后難以快速蒸發(fā)而以液態(tài)形式流入催化劑顆粒的間隙中,同時催化劑上微孔直徑較小(直徑范圍從10A到104A),在毛細(xì)作用下滲透到催化劑顆粒微孔的內(nèi)部。進入微孔后推進劑緩慢蒸發(fā),在催化劑微孔內(nèi)部形成氣液分界面,部分氣相推進劑通過質(zhì)量擴散到微孔壁面的活性位置分解得到氨氣、氫氣和氮氣等。這些分解氣體被催化劑外的液體推進劑堵塞在微孔內(nèi),隨著分解產(chǎn)物逐步增多微孔內(nèi)部壓力逐漸升高,推進劑分界面逐步向微孔外部移動,微孔內(nèi)部壓力最終會大于毛細(xì)壓力而使得分界面被擠壓到顆粒外部。大量的分解產(chǎn)物從微孔中溢出,氣相推進劑得以實現(xiàn)催化劑間隙到微孔活性位置的快速擴散分解,在這兩者共同的作用下催化床內(nèi)產(chǎn)生了較大的壓力峰。將冷起動延遲時間定義為從液體推進劑噴入催化床內(nèi)到氣體產(chǎn)物溢出催化劑微孔這段時間。

基于上述冷起動過程催化床內(nèi)部復(fù)雜過程的分析,下面將對催化劑顆粒內(nèi)部進行建模分析,計算氣液分界面在顆粒微孔內(nèi)部壓力的作用下逐漸推移到微孔外過程,得到單組元發(fā)動機冷起動延遲時間。

首先以催化劑顆粒內(nèi)部氣體為研究對象,分析催化劑顆粒微孔內(nèi)部壓升過程。根據(jù)微孔毛細(xì)理論[13],微孔內(nèi)部氣相產(chǎn)物沿顆粒徑向的平均速度通過泊肅葉方程描述如下:

式中vm為微孔內(nèi)部氣相的平均速度;R為大孔的特征半徑;μg為氣相產(chǎn)物的粘性系數(shù);δ為催化劑微孔彎曲因子;P為微孔內(nèi)壓力;r為流體徑向半徑。

式中J為氣相產(chǎn)物徑向單位面積上的質(zhì)量流量;γν為空隙分?jǐn)?shù),指單位體積顆粒內(nèi)氣相填充的比例;M為氣體產(chǎn)物的平均摩爾質(zhì)量;Ru為氣體常數(shù);T為氣體溫度,假設(shè)在較短的冷起動延遲時間內(nèi)氣相產(chǎn)物的溫度與初始溫度相差不大。

在微孔內(nèi)考慮氣體的連續(xù)方程:

聯(lián)立式(4)和式(5)可得催化劑顆粒內(nèi)部壓力:

在邊界條件方面,初始時刻微孔內(nèi)壓力為環(huán)境壓力;而在顆粒中心的壓力梯度為零,還有一個邊界條件是在微孔內(nèi)部移動的氣液分界面氣相質(zhì)量守恒,即微孔內(nèi)部氣相產(chǎn)物的質(zhì)量變化等于在分界面上蒸發(fā)和顆粒內(nèi)部擴散氣態(tài)肼的質(zhì)量變化。則初邊值條件如下:

式中DA為氣相推進劑在其分解產(chǎn)物中的擴散系數(shù),。

為了求解內(nèi)部壓力方程,還需要知道分界面隨時間的變化ai(t),以及在分界面上氣相肼的摩爾質(zhì)量沿微孔徑向分布。顆粒內(nèi)部氣液分界面的徑向位置可用描述液體在長管中流動規(guī)律的泊肅葉方程得到:

式中μl為液態(tài)推進劑粘性系數(shù);Pm為毛細(xì)壓力,;σ為液態(tài)推進劑表面張力,隨溫度變化;θ為接觸角。

根據(jù)在分界面后局部氣態(tài)肼考慮分解反應(yīng)的質(zhì)量守恒方程得到:

式中y為液體推進劑滲入催化劑顆微孔內(nèi)的長度;為液體推進劑的密度;ks(T)為推進劑的異相反應(yīng)速率常數(shù),ks=1010exp(-1389/T)。

這樣得到了完整的控制方程與初邊值條件,經(jīng)過無量綱化和小擾動法計算得到分界面位置ai(tlign)=0時,即氣體分解產(chǎn)物從微孔中溢出,此時得到了冷起動延遲時間tlign為[14]

式中 (ρA)vp為推進劑的飽和蒸汽密度,與溫度相關(guān)[1]。

1.2.2 催化床建壓過程模型

考慮到催化床內(nèi)復(fù)雜的物理化學(xué)過程,引入燃燒室時滯模型[15],假設(shè)在該轉(zhuǎn)化時間內(nèi)推進劑的轉(zhuǎn)化速率均勻,且分解后的氣體在任何瞬時分布均勻,即氨的解離過程也近似認(rèn)為與DT-3的分解過程同時進行。 則液體推進劑及氣體分解產(chǎn)物在催化床中的質(zhì)量變化分別為

式中mout是催化床出口的氣體質(zhì)量流量,;At為噴管的喉部截面積;ml、mg分別為催化床內(nèi)積存的液體推進劑質(zhì)量和氣體質(zhì)量;ql、分別為流入催化床內(nèi)的液體推進劑質(zhì)量流量和流出推力室的氣體質(zhì)量流量;τ為推進劑的轉(zhuǎn)化時間,可根據(jù)推進劑的分解延遲期試驗[16]得到。

由于催化床內(nèi)裝載著緊密的催化劑,因此推進劑通過催化床時需要克服阻力稱為床流阻,床流阻常采用經(jīng)驗公式計算[1]。

1.3 其他組件模型

用于單組元姿控發(fā)動機冷起動過程仿真的其他組件,如限流圈流阻模型,電磁閥考慮了閥芯動作過程的Γ-m模型,電磁閥后集液腔考慮氣體多變過程的充填排氣模型等[17]。限于文章篇幅不再詳述。

2 仿真模型的建立及校驗

某次單臺發(fā)動機地面試車考核了發(fā)動機在低溫下的工作特性,試驗系統(tǒng)如圖1所示。在該單機地面試驗系統(tǒng)中所采用的推進劑為液體單推-3[18](DT-3)。由于試驗系統(tǒng)中流量計的采樣頻率較低,不能反映發(fā)動機起動過程中的快速流量變化,故在此以推力室的室壓表征單組元發(fā)動機的冷起動特性。

圖1 單臺發(fā)動機地面試驗系統(tǒng)示意Fig.1 Engine Ground Test System Schematic

基于第1節(jié)所建立的典型組件數(shù)學(xué)模型,在MWorks軟件平臺上搭建該型單臺發(fā)動機的仿真模型,數(shù)值計算方法上采用4階Range-Kutta方法。計算結(jié)果和試驗測量結(jié)果如圖2所示。

圖2 單臺發(fā)動機冷起動過程仿真與試驗結(jié)果對比Fig.2 Simulation Results Compare with Ground Test Results of Engine Cold Starting Process

圖2中試驗值與計算值的比較結(jié)果表明,發(fā)動機室壓的計算值與試驗值變化趨勢吻合較好,但室壓計算值較試驗值偏高,引起這一誤差的主要原因是計算值是包含了催化床流阻部分的整個催化床平均壓力,而試驗測量值僅為推力室噴管入口處的壓力;并且發(fā)動機在冷起動時催化床分解效率是隨發(fā)動機工作時間逐漸升高的,而仿真中將分解效率設(shè)為恒定的穩(wěn)態(tài)分解效率,這2個因素導(dǎo)致了室壓計算值較高。其次,計算得到的室壓建壓曲線相比試驗建壓曲線更為陡峭,壓力峰的計算值稍快于試驗值,則是因為催化床模型采用了分解時滯模型反映推進劑的分解過程,認(rèn)為積存推進劑全部在瞬間發(fā)生快速分解,而發(fā)動機實際過程中積存的推進劑也是根據(jù)進入催化床順序先后逐步發(fā)生快速分解的,同時測量系統(tǒng)可能無法完全捕捉到室壓快速上升過程,考慮到這兩方面因素從而導(dǎo)致計算建壓過程快于試驗測量的建壓過程。

國外TRW公司對自發(fā)型催化劑Shell-405及其單組元肼發(fā)動機進行了充分的真空、地面點火試驗研究[19],得到了不同初始溫度下發(fā)動機的冷起動響應(yīng)時間,但文獻(xiàn)中未能給出試驗所用發(fā)動機系統(tǒng)的具體參數(shù),因此選取起動響應(yīng)時間作為關(guān)注量,起動響應(yīng)時間指推進劑進入催化床到推力室內(nèi)壓力達(dá)到穩(wěn)態(tài)室壓的5%所需時間。圖3給出了TRW公司地面試驗用發(fā)動機和本文所研究單組元發(fā)動機的冷起動響應(yīng)時間隨初始溫度變化的試驗值和仿真計算曲線。這里的溫度是指推進劑與催化床的初始溫度,考慮發(fā)動機地面試驗時第1次起動的催化劑與推進劑溫度相差不大,在此假設(shè)這兩者溫度一致。

圖3 響應(yīng)時間隨初始溫度變化對比Fig.3 Comparison of Test Data and Calculated Values of Response Time with Initial Temperature

本文的計算結(jié)果與文獻(xiàn)19中試驗數(shù)據(jù)的吻合情況較好。在初始溫度降至-20 ℃以下時,計算值與文獻(xiàn)19試驗值的誤差有增大的趨勢,分析認(rèn)為是該溫度范圍接近推進劑的冰點-30 ℃,數(shù)學(xué)模型中所采用的推進劑物性參數(shù)計算公式在該范圍內(nèi)適用性變差所導(dǎo)致。

3 仿真結(jié)果分析

3.1 單組元發(fā)動機冷起動過程研究

以第2節(jié)驗證過的仿真模型為基礎(chǔ),以發(fā)動機設(shè)計溫度下的狀態(tài)參數(shù)值為參考,討論和研究初始溫度對發(fā)動機冷起動過程中室壓和響應(yīng)時間的影響。圖4給出了不同初始溫度下發(fā)動機冷起動過程中無量綱室壓cP隨時間的變化。由圖4可知,催化床與推進劑的初始溫度對發(fā)動機的冷起動過程影響較大。初始溫度越低,起動響應(yīng)時間越長,冷起動壓力峰越高,發(fā)動機的冷起動性能越差,降低了發(fā)動機的結(jié)構(gòu)可靠性。還可以看到,在出現(xiàn)冷起動壓力峰后壓力曲線還會伴隨“凹坑”。這種“凹坑”現(xiàn)象的產(chǎn)生是由于高壓力峰的存在,使得整個系統(tǒng)壓差變化,發(fā)動機供應(yīng)系統(tǒng)出現(xiàn)短暫的推進劑供應(yīng)不足導(dǎo)致推力室產(chǎn)生了較大的壓力波動。冷起動壓力峰越高,所產(chǎn)生的壓力“凹坑”越突出。在初始溫度為10 ℃時冷起動壓力峰已經(jīng)高達(dá)設(shè)計值的4.5倍左右,并且伴隨了更為明顯的“凹坑”,壓力的最低值已經(jīng)接近0值而后緩慢爬升至設(shè)計穩(wěn)定值,導(dǎo)致發(fā)動機達(dá)到穩(wěn)態(tài)工作壓力的時間明顯延長,嚴(yán)重影響了發(fā)動機的工作可靠性。

圖4 初始溫度對冷起動過程的影響Fig.4 Influence of Initial Temperature on Cold Starting Process

進一步模擬不同初始溫度下發(fā)動機的冷起動過程,圖5為初始溫度在0~30 ℃推力裝置的無量綱冷起動壓力峰值及響應(yīng)時間隨初始溫度的變化。可以看到,冷起動壓力峰值和冷起動響應(yīng)時間呈現(xiàn)正相關(guān)的關(guān)系。這意味著冷起動響應(yīng)越快,冷起動壓力峰越低。

圖5 壓力峰值及峰值時刻隨初始溫度的變化Fig.5 Influence of Initial Temperature on Pressure Spike and Response Time

隨著初始溫度的升高,冷起動壓力峰和響應(yīng)時間隨溫度的變化近似呈指數(shù)變化。對于所研究地面試驗系統(tǒng)中的單組元發(fā)動機,當(dāng)催化床和推進劑的初始溫度從0 ℃升高到30 ℃時,壓力峰值降低了85.8%,響應(yīng)時間減小了67.5%。催化床和推進劑的初始溫度顯著影響著單組元發(fā)動機的冷起動過程,尤其是在初始溫度較低時,即便是較小程度的提高初始溫度也可以明顯地改善發(fā)動機冷起動延遲長、壓力峰高的現(xiàn)象。

3.2 初始溫度對冷起動延遲的作用機制分析

由冷起動延遲過程分析可知,初始溫度主要影響著包括催化分解反應(yīng)、液相推進劑蒸發(fā)以及氣相推進劑擴散這3種過程以及液相推進劑的表面張力,為進一步比較分析4種因素對冷起動延遲的影響程度,以20 ℃時各參數(shù)值為參照,對式(10)進行無量綱化后再取以為底的對數(shù)得到:

由于log函數(shù)不改變原函數(shù)的單調(diào)性,因此式(14)中右側(cè)分別代表著表面張力、初始溫度、飽和密度、氣相推進劑擴散系數(shù)以及分解反應(yīng)速率對冷起動延遲時間的影響因子。由于推進劑的冰點為-30 ℃,因此計算初始溫度在-20~30 ℃情況下各影響因子的變化如圖6所示。

圖6 不同初始溫度下各項的影響因子Fig.6 Influencing Factors at Different Initial Temperature

由圖6可知,隨著初始溫度的升高,飽和密度項和分解反應(yīng)項呈現(xiàn)一定的下降趨勢,而表面張力、擴散系數(shù)和溫度3項相應(yīng)的影響程度逐漸提高。其中表面張力的作用隨溫度上升提高明顯,在大約10 ℃其作用已經(jīng)超過了擴散系數(shù)和溫度項的影響,但是這3項的影響因子都小于5%。

在所研究的溫度范圍內(nèi),飽和密度項的影響因子一直處在較高值,并且超過其他4項之和的2倍以上。因此液體推進劑的蒸發(fā)過程是單組元肼發(fā)動機冷起動過程中導(dǎo)致延遲較長、壓力峰較高的主導(dǎo)因素;其次分解反應(yīng)項的影響因子超過了10%,表明分解反應(yīng)過程在冷起動過程中的影響也不容忽略。

3.3 不同推進劑的冷起動溫度敏感性

肼-70、無水肼和DT-3是3種常用肼類單組元推進劑[16,20],本節(jié)討論相同條件下采用不同推進劑的單組元發(fā)動機冷起動特性。圖7為采用3種推進劑的發(fā)動機在-15~30 ℃初始溫度條件下的冷起動延遲時間。由圖7可知,在所計算溫度范圍內(nèi),無水肼發(fā)動機冷起動延遲最長,且隨初始溫度的降低,冷起動延遲的差距愈加明顯,-15 ℃時無水肼和DT-3的差值已達(dá)到1倍以上。考慮3.1節(jié)得到的冷起動壓力峰與冷起動延遲時間正相關(guān),故認(rèn)為采用DT-3的單組元發(fā)動機在相同條件下較無水肼和肼-70具有更好的冷起動特性。

圖7 冷起動延遲時間隨初始溫度變化Fig.7 Influence of Initial Temperature on Response Time of Cold Starting

4 結(jié) 論

本文建立了反映冷起動延遲的單組元姿控發(fā)動機典型組件數(shù)學(xué)模型,分析了單組元姿控軌發(fā)動機在不同初始溫度下的冷起動特性,得到以下結(jié)論:a)單組元姿控發(fā)動機冷起動過程仿真模型的計算值與試驗值吻合較好,冷起動延遲時間與文獻(xiàn)中試驗值趨勢相同,驗證了所建立催化床冷起動延遲數(shù)學(xué)模型和相應(yīng)仿真模型,可為工程設(shè)計提供理論支持。b)肼類單組元發(fā)動機冷起動過程中壓力峰與起動延遲正相關(guān),并都與初始溫度近似指數(shù)變化的關(guān)系。c)相比于低溫下氣態(tài)推進劑的擴散過程和催化分解過程,液體推進劑進入催化床后的蒸發(fā)氣化過程是單組元發(fā)動機在較低溫度起動時延遲長、壓力峰高的主導(dǎo)因素。

相較于肼-70和無水肼,采用單推-3的單組元催化分解發(fā)動機具有更好的冷起動特性。

猜你喜歡
催化劑發(fā)動機模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
發(fā)動機空中起動包線擴展試飛組織與實施
直接轉(zhuǎn)化CO2和H2為甲醇的新催化劑
3D打印中的模型分割與打包
新型釩基催化劑催化降解氣相二噁英
掌握情欲催化劑
Coco薇(2016年2期)2016-03-22 02:45:06
V2O5-WO3/TiO2脫硝催化劑回收研究進展
新一代MTU2000發(fā)動機系列
主站蜘蛛池模板: 中文字幕 91| 国产一二三区视频| 国产91在线|日本| 欧美一级高清视频在线播放| 欧美成人精品在线| 欧美国产精品不卡在线观看| 亚洲va欧美ⅴa国产va影院| 国产永久免费视频m3u8| 69国产精品视频免费| 毛片网站观看| 亚洲欧美激情小说另类| 大学生久久香蕉国产线观看| 毛片在线看网站| 亚洲日本www| 呦视频在线一区二区三区| 亚洲天堂久久| yy6080理论大片一级久久| 亚洲中文字幕av无码区| 潮喷在线无码白浆| 亚洲欧美人成人让影院| 青青青视频蜜桃一区二区| 亚洲精品成人7777在线观看| 熟妇人妻无乱码中文字幕真矢织江 | 国产人成在线观看| 91丝袜美腿高跟国产极品老师| 永久免费无码成人网站| 六月婷婷激情综合| 欧美全免费aaaaaa特黄在线| 中文字幕色站| 四虎成人精品在永久免费| a毛片在线| 91精品啪在线观看国产91| 青青青伊人色综合久久| 国产在线八区| 国产欧美日韩视频怡春院| 成人精品亚洲| 欧美国产精品不卡在线观看| 全色黄大色大片免费久久老太| 亚洲人在线| 亚洲国产成人精品青青草原| 色欲综合久久中文字幕网| 在线观看的黄网| 精品国产欧美精品v| 国产亚洲高清视频| 人妻中文久热无码丝袜| 99热这里只有精品5| 综合天天色| 中文字幕在线永久在线视频2020| 婷婷久久综合九色综合88| 素人激情视频福利| 亚洲丝袜第一页| 午夜一区二区三区| 在线精品欧美日韩| 潮喷在线无码白浆| 99视频全部免费| 日韩欧美国产精品| 99一级毛片| 国产精欧美一区二区三区| www.亚洲天堂| 中文字幕久久亚洲一区| 超碰aⅴ人人做人人爽欧美| 四虎成人免费毛片| 伊大人香蕉久久网欧美| 久久一色本道亚洲| 亚洲中久无码永久在线观看软件| 九九九国产| 欧美第二区| 久久天天躁夜夜躁狠狠| 国产自在线播放| 色婷婷成人网| 日本精品影院| 欧美翘臀一区二区三区| 国产成人无码AV在线播放动漫| 一本色道久久88亚洲综合| 国产精品无码一二三视频| 亚洲码在线中文在线观看| 伊人色天堂| 国产三级毛片| 国产网站黄| 日韩AV无码免费一二三区| 国产视频只有无码精品| 亚洲天堂伊人|