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

電推進(jìn)系統(tǒng)Xe物理特性計(jì)算方法

2016-02-13 07:43:58陳濤劉國西宋飛武蔥蘢
中國空間科學(xué)技術(shù) 2016年1期
關(guān)鍵詞:物理

陳濤,劉國西,宋飛,武蔥蘢

北京控制工程研究所,北京100190

電推進(jìn)系統(tǒng)Xe物理特性計(jì)算方法

陳濤*,劉國西,宋飛,武蔥蘢

北京控制工程研究所,北京100190

空間電推進(jìn)系統(tǒng)的工質(zhì)為Xe,其工作溫度范圍為-30~45℃,該范圍覆蓋Xe的臨界點(diǎn)。在臨界點(diǎn)附近,Xe可能呈現(xiàn)出多種形態(tài),且對溫度和壓力變化十分敏感,采用傳統(tǒng)狀態(tài)方程在上述范圍內(nèi)計(jì)算其物理參數(shù)偏差最大可達(dá)30%。為解決這一問題,以CH4作為參考流體,建立了一種基于對比態(tài)原理的Xe物理特性計(jì)算方法。該方法能夠?qū)Π庀唷⒁合唷⒊R界區(qū)域的所有狀態(tài)Xe物理性質(zhì)進(jìn)行準(zhǔn)確計(jì)算。試驗(yàn)研究與國外數(shù)據(jù)對比的結(jié)果表明,在整個壓力-溫度范圍內(nèi),計(jì)算誤差小于0.5%。

電推進(jìn);Xe;超臨界;對比態(tài)原理;物理特性;狀態(tài)方程

電推進(jìn)技術(shù)是一種利用電能加速推進(jìn)工質(zhì),從而實(shí)現(xiàn)高比沖的先進(jìn)航天器推進(jìn)技術(shù),由于其高比沖的優(yōu)勢,在先進(jìn)國家的大型衛(wèi)星平臺及微小衛(wèi)星平臺上得到日漸廣泛的應(yīng)用[1-8]。

Xe作為電推進(jìn)系統(tǒng)的工質(zhì),具有分子量大、熔點(diǎn)高、沸點(diǎn)高及臨界點(diǎn)高等特點(diǎn)。考慮到系統(tǒng)的穩(wěn)態(tài)工作、充氣加壓等過程,其工作溫度范圍為-30~45℃,壓力范圍為1~30 MPa。整個工作范圍覆蓋氙的臨界點(diǎn)。在臨界點(diǎn)附近,Xe狀態(tài)對溫度和壓力變化十分敏感,分子間的作用力不可忽略,表現(xiàn)出強(qiáng)烈的真實(shí)流體特性。除此之外,在上述范圍內(nèi),Xe可能呈現(xiàn)出液態(tài)、氣態(tài)、氣液混合態(tài)以及超臨界態(tài)等多種形態(tài)。

目前,比較常用的描述真實(shí)流體性質(zhì)的狀態(tài)方程可以分為兩大類:解析形式和非解析形式。解析形式的狀態(tài)方程主要特點(diǎn)是當(dāng)溫度T、壓強(qiáng)p給定時,體積V可以解析地求出,即方程中V的次數(shù)不超過4次,例如Redlich-Kwong方程[9]、Peng-Robinson方程[10]、Soave-Redlich-Kwong(SRK)方程[9]等;非解析形式的狀態(tài)方程由于V次數(shù)較高,只能采用數(shù)值方法進(jìn)行求解,例如Benedict-Webb-Rubin(BWR)/Modified Benedict-Webb-Rubin(MBWR)方程、Wagner方程等。

通過比較發(fā)現(xiàn),以上幾種傳統(tǒng)狀態(tài)方程無法同時準(zhǔn)確計(jì)算各種形態(tài)下氙單質(zhì)的p-V-T關(guān)系,計(jì)算偏差最大可達(dá)30%,精度無法滿足應(yīng)用需求。而精確獲得上述關(guān)系對于Xe的充氣加壓、發(fā)射場加注及在軌監(jiān)測都極為重要[11-12]。

因此,深入了解Xe物理特性,對電推進(jìn)系統(tǒng)設(shè)計(jì)具有重要意義[13]。目前,國內(nèi)外沒有公開發(fā)表的計(jì)算方法。為了解決這一問題,本文主要針對電推進(jìn)系統(tǒng)工質(zhì)Xe的物理特性開展理論分析及計(jì)算方法研究。

1 理論模型及計(jì)算方法

典型物質(zhì)的三維熱力學(xué)曲面如圖1所示[9]。按照流體所處的不同狀態(tài)可以分為固相、液相、氣相三種不同相存在。在工程應(yīng)用中,常常需要對不同狀態(tài)下物質(zhì)的物理特性進(jìn)行計(jì)算。

在大量的實(shí)際應(yīng)用中,研究人員提出了各種各樣的狀態(tài)方程來描述p-V-T之間的關(guān)系。狀態(tài)方程中最簡單、最常見的是理想氣體狀態(tài)方程,主要適用于氣體密度較低(即壓力不太高、溫度不太低)分子間相互作用可以忽略不計(jì)的情形。當(dāng)壓力較高或溫度接近飽和氣體溫度時,分子間的作用力不可忽略,為了能描述真實(shí)工質(zhì)的熱物理屬性,必須采用在低溫和高壓下都能對真實(shí)流體性質(zhì)進(jìn)行刻畫的狀態(tài)方程。

由于傳統(tǒng)真實(shí)流體模型在高壓區(qū)、液態(tài)區(qū)和臨界點(diǎn)附近計(jì)算結(jié)果偏差較大,本研究中采用對比態(tài)原理結(jié)合BWR方程[14]來計(jì)算Xe的p-V-T關(guān)系。

1.1 對比態(tài)原理

對比態(tài)原理(Corresponding State Principle, CSP)是指:對所有物質(zhì)來說,其狀態(tài)變量經(jīng)過恰當(dāng)?shù)臒o量綱化處理后,都遵從普遍的變化規(guī)律。例如,對p、V、T都以臨界點(diǎn)參數(shù)進(jìn)行歸一化處理后(即pr=p/pc,Tr=T/Tc,Vr= V/Vc),pr-Vr-Tr關(guān)系對于所有同類極性物質(zhì)都一樣。

對比態(tài)原理是目前應(yīng)用最廣泛的物性估算方法。也就是說,在獲得參考流體的物性試驗(yàn)數(shù)據(jù)后,就可以采用對應(yīng)的對比態(tài)理論來計(jì)算目標(biāo)流體的物性。對于具有不同分子結(jié)構(gòu)而顯現(xiàn)出不同極性的流體,按照對比態(tài)理論中使用參數(shù)的個數(shù),可以分為2參數(shù)CSP、3參數(shù)CSP和4參數(shù)CSP。不同的模型適用于不同的物質(zhì)類型。

由于Xe屬于單原子球形分子,本文選取同極性的具有球形分子結(jié)構(gòu)的CH4作為參考流體,采用2參數(shù)對比態(tài)原理結(jié)合BWR方程來計(jì)算Xe的物理特性。

2參數(shù)CSP是指,僅用兩個特征性質(zhì)來使?fàn)顟B(tài)條件無量綱化(不需引入新的無量綱特征參數(shù)),其特征參數(shù)為Tc和pc或者Tc和Vc。只有單原子物質(zhì)如Ar、Kr、Xe或簡單流體才精確地遵循2參數(shù)對比態(tài)原理,其他物質(zhì)都有一些偏離。

1.2 BWR狀態(tài)方程

BWR方程[14]能同時對氣相、液相區(qū)域的物性進(jìn)行準(zhǔn)確的估算。方程將溫度多項(xiàng)式與密度的冪級數(shù)和指數(shù)相結(jié)合,其具體表達(dá)形式為

其中,15個系數(shù)An的具體形式如表1所示,除了參數(shù)r外,與溫度有關(guān)的函數(shù)共涉及32個擬合參數(shù)。

表1 BWR方程中系數(shù)項(xiàng)Table 1 Coefficient items in BWR equation

BWR方程中共有32個擬合參數(shù),文獻(xiàn)基于CH4大范圍的試驗(yàn)數(shù)據(jù)給出了BWR模型中對應(yīng)的32個參數(shù)取值,見表2。文獻(xiàn)[15]指出,采用該組參數(shù)取值對CH4物性進(jìn)行計(jì)算,最大誤差不超過0.45%。

1.3 氣液狀態(tài)辨識方法

由于計(jì)算過程中涉及氣、液、氣液混合等多種狀態(tài),為了對流體狀態(tài)進(jìn)行辨識,本文基于美國國家標(biāo)準(zhǔn)技術(shù)研究院(NIST)數(shù)據(jù)庫中CH4的飽和點(diǎn)數(shù)據(jù),來獲得CH4飽和狀態(tài)下壓力和溫度的對應(yīng)關(guān)系。擬合得到的氣化線為

式中:p1=1.383×10-8,p2=-3.04×10-6,p3=0.0002654,p4=-0.01427,p5=0.4669。對于給定的溫度,可以得到該點(diǎn)下的飽和蒸氣壓。若流體壓力高于該數(shù)值,則流體處于液態(tài),反之處于氣態(tài)。

在此基礎(chǔ)上,Xe的氣液狀態(tài)辨識通過式(1)結(jié)合對比態(tài)理論實(shí)現(xiàn)。

1.4 模型求解流程

采用對比態(tài)原理對Xe密度進(jìn)行求解時,首先根據(jù)Xe的溫度、壓力參數(shù),無量綱化處理后,采用對比態(tài)理論獲得參考流體CH4對應(yīng)的參數(shù),在此基礎(chǔ)上,反求BWR方程得到CH4的密度,再根據(jù)對比態(tài)理論求得Xe的狀態(tài)參數(shù)。

2 試驗(yàn)系統(tǒng)及試驗(yàn)方案

為了對理論模型的正確性進(jìn)行驗(yàn)證,同時獲得Xe第一手試驗(yàn)數(shù)據(jù),開展Xe物理特性試驗(yàn)研究。

2.1 系統(tǒng)組成

Xe物理特性試驗(yàn)設(shè)備組成及工作原理如圖2所示,設(shè)備主要包括高純置換氣體、氣體凈化器、回收氣瓶、Xe運(yùn)輸罐、稱重裝置、高低溫箱、真空泵、高壓氣瓶、高精度測量與控制系統(tǒng)。其中高壓氣瓶是本試驗(yàn)中的主要試驗(yàn)裝置,通過充入不同質(zhì)量的Xe實(shí)現(xiàn)不同密度條件下的參數(shù)測量;高精度測量系統(tǒng)用于監(jiān)測Xe運(yùn)輸罐質(zhì)量,高低溫箱中Xe的質(zhì)量、記錄壓力容器內(nèi)溫度、壓力變化數(shù)據(jù)采集頻率不低于50 Hz;高低溫箱用于控制試驗(yàn)中高壓氣瓶所需恒定溫度,恒溫精度不低于0.5℃。

試驗(yàn)介質(zhì)選用高純Xe,純度≥99.995%。

圖2 Xe物理特性設(shè)備原理Fig.2 Schematic of the xenon physical property experiment equipment

2.2 試驗(yàn)方法

首先對高壓氣瓶容積進(jìn)行標(biāo)定。獲得其精確容積之后,采用定容法測定Xe的p-V-T關(guān)系。在依次完成抽真空、置換、高純置換之后,對高壓氣瓶吸入一定量的高純度Xe,后啟動高低溫箱對高壓氣瓶進(jìn)行升溫、降溫的過程,依次進(jìn)行該密度下Xe壓力隨溫度的變化規(guī)律測量;之后對高壓氣瓶再充入定量Xe,再重復(fù)上述降溫回溫過程,實(shí)現(xiàn)另一密度條件下的數(shù)據(jù)測量。以此類推,依次完成0.2 kg/L、0.4 kg/L、0.6 kg/L、0.8 kg/L、1.0 kg/L、1.2 kg/L、1.4 kg/L和1.6 kg/L密度條件下的數(shù)據(jù)測量。

3 計(jì)算結(jié)果及試驗(yàn)對比分析

根據(jù)BWR模型和對比態(tài)理論開展Xe物理特性計(jì)算,并與NIST數(shù)據(jù)[16]及試驗(yàn)結(jié)果進(jìn)行對比以驗(yàn)證計(jì)算方法的正確性。

Xe的臨界點(diǎn)參數(shù)如表3所示。當(dāng)溫度高于Tc時,Xe只處于氣體狀態(tài)。當(dāng)溫度壓力都高于臨界點(diǎn)參數(shù)時,此時稱為超臨界狀態(tài)。這一狀態(tài)下,流體的密度表現(xiàn)得像液態(tài),但輸運(yùn)特性表現(xiàn)得像氣態(tài)。

參考流體CH4的臨界點(diǎn)參數(shù)見表3。

表3 Xe及參考流體CH4的臨界點(diǎn)參數(shù)Table 3 Critical point parameters of Xe and CH4

3.1 飽和液線及飽和氣線

當(dāng)流體從液態(tài)變?yōu)闅鈶B(tài)時,會經(jīng)歷如下幾個狀態(tài):液態(tài)、飽和液體態(tài)、氣液混合態(tài)、飽和蒸氣態(tài)和氣態(tài)。采用對比態(tài)原理計(jì)算得到的Xe氣化線如圖3所示,計(jì)算溫度范圍為170~289.74 K,覆蓋從三相點(diǎn)到臨界點(diǎn)所有工況。可以看到Xe在飽和態(tài)或者氣液混合態(tài)壓力隨溫度的變化關(guān)系。對應(yīng)的飽和氣線和飽和液線分別如圖4和圖5所示。結(jié)果表明,采用本文提出的方法與NIST數(shù)據(jù)符合較好,飽和線上密度計(jì)算最大誤差小于0.1%,能較好地分辨Xe給定壓力和溫度下的物理狀態(tài)。

3.2 熱力學(xué)曲線

為了驗(yàn)證程序在大范圍內(nèi)計(jì)算結(jié)果的有效性,對橫跨液相、氣液混合、氣相區(qū)的溫度和壓力組合進(jìn)行密度計(jì)算,來獲得Xe的熱力學(xué)曲線,并與NIST數(shù)據(jù)[16]進(jìn)行對比。溫度計(jì)算范圍170~600 K,壓力計(jì)算范圍1~30 MPa,計(jì)算結(jié)果如圖6所示。可以看到,計(jì)算結(jié)果與NIST數(shù)據(jù)完全吻合,能準(zhǔn)確辨識流體狀態(tài)。證明計(jì)算結(jié)果準(zhǔn)確可信,可以用于給定溫度T、壓力p條件下電推進(jìn)系統(tǒng)工質(zhì)Xe的密度計(jì)算。

圖3 Xe的氣化線Fig.3 Evaporating line of xenon

圖4 飽和液線Fig.4 Saturation liquid line of xenon

圖5 飽和氣線Fig.5 Saturation gas line of xenon

圖6 1~30 MPa范圍內(nèi)Xe的T-ρ關(guān)系Fig.6T-ρrelationship in 1~30 MPa for Xe

對于工程應(yīng)用中的典型密度值,將計(jì)算數(shù)據(jù)(ρth為理論計(jì)算結(jié)果)匯總,并與試驗(yàn)結(jié)果進(jìn)行對比,如圖7所示,平均誤差小于0.5%。可以看出,各密度條件下,計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)符合良好,證明了計(jì)算方法可以完全滿足電推進(jìn)系統(tǒng)設(shè)計(jì)需求。

圖7 計(jì)算數(shù)據(jù)與試驗(yàn)數(shù)據(jù)對比Fig.7 Comparison between calculation and the experiment

由圖7中可以看到,在不同密度條件下,當(dāng)溫度逐步降低時,壓力都呈下降趨勢。同時,高壓氣瓶中Xe密度在臨界密度(ρc=1.1029 kg/L)附近和遠(yuǎn)離臨界點(diǎn)時,兩者的變化趨勢存在顯著不同。例如,ρ=2.0 kg/L時,圖7中p-T關(guān)系存在明顯的轉(zhuǎn)折點(diǎn),真實(shí)流體的特性表現(xiàn)明顯;而ρ=1.0 kg/L,壓力隨溫度變化平緩,曲線較為光滑。

結(jié)合Xe三相圖在密度壓力平面的投影,即圖8中的ρ-P圖可以發(fā)現(xiàn),在不同密度條件下,隨著溫度(或壓力)的降低,Xe經(jīng)歷的相變并不相同。當(dāng)密度大于臨界密度1.102 9 kg/L時,Xe從40℃降低到-30℃過程中,會經(jīng)歷從超臨界態(tài)變?yōu)橐簯B(tài),再到氣液混合態(tài)的情形;而當(dāng)密度小于1.102 9 kg/L時,會經(jīng)歷從超臨界態(tài)變?yōu)闅鈶B(tài),再到氣液混合態(tài)的情形。這是圖7中p-T關(guān)系在不同密度條件下變化趨勢存在區(qū)別的原因。

圖8 Xeρ-p關(guān)系Fig.8 Relationship betweenρ-pfor Xe

4 結(jié)束語

本文對電推進(jìn)系統(tǒng)工質(zhì)Xe的物理特性進(jìn)行了理論分析,建立了對應(yīng)的真實(shí)流體狀態(tài)方程,通過開展試驗(yàn)研究、并國外NIST數(shù)據(jù)庫對比分析,結(jié)果表明:

1)建立的計(jì)算方法原理簡單,計(jì)算方便,基于現(xiàn)有數(shù)據(jù),無需提供額外參數(shù)。

2)在整個壓力-溫度范圍內(nèi),該方法能夠準(zhǔn)確對氣相、液相、超臨界區(qū)域的Xe物理性質(zhì)進(jìn)行計(jì)算,誤差小于0.5%,能夠完全滿足電推進(jìn)貯供系統(tǒng)設(shè)計(jì)的應(yīng)用需求。

3)研究過程中采用試驗(yàn)方法合理有效,能用于Xe高精度物性測量。

4)由于臨界溫度較高,在電推進(jìn)系統(tǒng)的使用溫度范圍內(nèi),氙的真實(shí)流體的特性表現(xiàn)明顯,在不同密度條件下,p-T關(guān)系存在典型的差異。

References)

[1] ZHANG T P,WANG L,WU X M.Electric propulsion development for DFH-4SP satellite platform[C]∥Proceedings of 65th IAC,Toronto,Canada,IAC-14-C4.4.2,2014.

[2] VALERIE C T,JOSEPH M M,GBROWN J M.The dawn spacecraft[J].Space Sci.Rev.,2011,163(1): 175-249.

[3] ADAM S,ALEC D G,PETER Y P.Performance of a helicon Hall thruster operating with xenon,argon and nitrogen[J].Journal of Propulsion and Power, 2014,30(3):664-671.

[4] HANI K,WEN S H,THOMAS H.Overview of the development of the solar electric propulsion technology demonstration mission 12.5-kW Hall thruster[C]∥50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, Cleveland,OH,AIAA 2014-3898,2014.

[5] DAN L,RAANAN E,GAL A.CAM200 Hall thrusterdevelopment overview[C]∥Proceedings of 66th IAC, Jerusalem,Israel,IAC-15-C4.4.4,2015.

[6] CHIEN K R.L-3 communications ETI electric propulsion overview[C]∥The 29th International Electric PropulsionConference,Princeton University,October 31-November 4,IEPC2005-315,2005.

[7] SEMENKIN A V.Overview of electric propulsion activity in Russia[C]∥The 30th International Electric Propulsion Conference,Florence,Italy,September 17-20,IEPC 2007-275,2007.

[8] KOPPEL C R,MARCHANDISE F,PRIOUL M.The SMART-1 electric propulsion subsystem around the Moon:in flight experience[C]∥41st AIAA/ASME/ SAE/ASEE Joint Propulsion Conference&Exhibit,Joint Propulsion Conference,Tucson,Arizona,July 10-13, AIAA 2005-3671,2005.

[9] BRUCE E P,JOHN M P,JOHN P O.Properties of gases and liquids[M].New York:McGraw-Hill, 2001:111-163.

[10] PENG D Y,ROBINSON B D.A new two constant equation of state[J].Industry Engineering Chemical Fundamentals,1976,15(1):59-64.

[11] JOHN W D,JOSEPH C,ANTHONY D.Advanced Xenon feed system(AXFS)development and hot-fire testing[C]∥45th AIAA/ASME/SAE/ASEE Joint Propulsion Conference&Exhibit,Denver,Colorado, August 2-5,2009,AIAA 2009-4910.

[12] GANAPATHI G B,ENGELBRECHT C S.Performance of the xenon feed system on Deep Space One[J].Journal of Spacecraft and Rockets,2000,37(3):392-298.

[13] OSBORN M F,NETWALL C J.Higher performance xenon flow system(XFS)optimized for low mass, volume and cost[C]∥45th AIAA Joint Propulsion Conference,Denver,Colorado,August 2-5,2009: AIAA 2009-4909.

[14] JACOBSEN R T,STEWART R B.Thermodynamic properties of nitrogen including liquid and vapor phases from 63 K to 2 000 K with pressure to 10000 bar[J]. Journal of Physical and Chemical Reference Data,1973, 2(4):757.

[15] ZONG N.Modeling and simulation of cryogenic fluid injection and mixing dynamics under supercritical conditions[D].Pennsylvania:The Pennsylvania State University,2005.

[16] NIST.NIST standard reference database[EB/OL].http:∥webbok.nist.gov.

(編輯:車曉玲)

A theoretical method for xenon's physical property used in electrical propulsion system

CHEN Tao*,LIU Guoxi,SONG Fei,WU Conglong
Beijing Institute of Control Engineering,Beijing 100190,China

The possible operation temperature of propellant xenon in electrical propulsion system is between-30℃and 45℃,which covers the xenon's critical temperature.At the critical temperature,xenon's state is sensitive to the operating pressure and the ambient temperature.The propellant may exhibit different forms.These characteristics make the traditional equation of state ineffective to calculate thep-V-Trelationship in the above temperature range(max error may rise to 30%).To solve this problem,a new calculation method based on the corresponding state principle was put forward for physical property estimation.The theoretical outcomes were compared with experiment data and available database.The results show that among the whole temperature and pressure range,the method is capable to predict the xenon's physical property in gas,liquid,supercritical state with an error less than 0.5%.

electrical propulsion;xenon;supercritical;corresponding state principle; physical properties;equation of state

V43

:A

10.3780/j.issn.1000-758X.2016.0018

2015-11-09;

:2015-12-25;錄用日期:2016-01-18;< class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間:2016-02-24 13:42:48

http:∥www.cnki.net/kcms/detail/11.1859.V.20160224.1342.015.html

*

:陳濤(1986-),男,博士,工程師,364993478@163.com,主要研究方向?yàn)殡娡七M(jìn)系統(tǒng)設(shè)計(jì)

陳濤,劉國西,宋飛,等.電推進(jìn)系統(tǒng)Xe物理特性計(jì)算方法[J].中國空間科學(xué)技術(shù),2016,36(1):113-119.

CHEN T,LIU G X,SONG F,et al.A theoretical method for xenon′s physical property used in electrical propulsion system[J].Chinese Space Science and Technology,2016,36(1):113-119(in Chinese).

http:∥zgkj.cast.cn

猜你喜歡
物理
物理中的影和像
只因是物理
井岡教育(2022年2期)2022-10-14 03:11:44
高考物理模擬試題(五)
高考物理模擬試題(二)
高考物理模擬試題(四)
高考物理模擬試題(三)
留言板
如何打造高效物理復(fù)習(xí)課——以“壓強(qiáng)”復(fù)習(xí)課為例
處處留心皆物理
我心中的物理
主站蜘蛛池模板: 中文字幕人妻无码系列第三区| 天天色综网| 色成人综合| 香蕉eeww99国产在线观看| 91九色最新地址| 欧美一级视频免费| 亚洲a级毛片| 亚洲第一网站男人都懂| 成人毛片免费在线观看| 国产高清在线观看| 亚洲精品桃花岛av在线| 综合久久五月天| 亚洲中文无码av永久伊人| 国产一区成人| 国产成人高清在线精品| 亚洲无卡视频| 成人午夜网址| 91在线国内在线播放老师| 欧美三级自拍| 男人的天堂久久精品激情| 国产综合在线观看视频| 在线观看无码a∨| 中文字幕第4页| 国产精品青青| 色久综合在线| www亚洲天堂| 亚洲第一成网站| 亚洲国产精品无码AV| 天堂av综合网| 亚洲国产成人麻豆精品| 亚洲黄色成人| 毛片网站在线播放| 婷婷久久综合九色综合88| 欧美一区国产| 欧美一道本| 精品国产自在在线在线观看| 伊人久久青草青青综合| 欧美精品另类| 在线免费不卡视频| 亚洲日韩图片专区第1页| 国产91丝袜在线播放动漫 | 91无码网站| 天天色天天综合| 制服丝袜亚洲| 91麻豆国产在线| 欧美a在线| 欧类av怡春院| 99人体免费视频| 97在线碰| 丝袜久久剧情精品国产| 久久久受www免费人成| 亚洲青涩在线| 免费全部高H视频无码无遮掩| 欧洲亚洲一区| 亚洲毛片网站| 久久无码av三级| 中文字幕一区二区人妻电影| 日本www色视频| 四虎精品黑人视频| 综1合AV在线播放| 国产第一页屁屁影院| 国产精品无码久久久久久| 幺女国产一级毛片| 国产成人亚洲精品无码电影| 一区二区日韩国产精久久| 伊人丁香五月天久久综合| 日韩大片免费观看视频播放| 中文字幕在线看| 好吊日免费视频| 国产福利一区在线| 国产一区二区精品高清在线观看| 狂欢视频在线观看不卡| 沈阳少妇高潮在线| 欧美日韩专区| 91精品日韩人妻无码久久| 精品国产一区91在线| 嫩草国产在线| 91无码人妻精品一区| 丝袜久久剧情精品国产| 91年精品国产福利线观看久久 | 自拍欧美亚洲| 国产成a人片在线播放|