亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        新型二元混合制冷劑(R290+R227ea)氣液相平衡研究

        2014-03-07 03:48:40陳秀萍祁影霞趙勝喜張華
        制冷學報 2014年5期
        關(guān)鍵詞:制冷劑物性電荷

        陳秀萍 祁影霞 趙勝喜 張華

        (上海理工大學能源與動力工程學院 上海 200093)

        新型二元混合制冷劑(R290+R227ea)氣液相平衡研究

        陳秀萍 祁影霞 趙勝喜 張華

        (上海理工大學能源與動力工程學院 上海 200093)

        混合制冷劑氣液相平衡數(shù)據(jù)是新型制冷劑的熱力學參數(shù)的重要組成部分。理論預測和實驗測定混合工質(zhì)的氣液相平衡數(shù)據(jù)成為研究混合工質(zhì)替代問題的一項迫切需要。本文基于量子化學理論,采用真實溶劑似導體模型(COSMO-RS)模擬了二元混合制冷劑(R290+R227ea)氣液相平衡性質(zhì),模擬的結(jié)果與實驗結(jié)果有很好的一致性。表明運用COSMO-RS模型預測混合制冷劑的氣液相平衡是可行的。

        氣液相平衡;量子化學;COSMO-RS;R290+R227ea

        近些年,在《聯(lián)合國氣候變化框架公約》和《京都議定書》的基礎(chǔ)上,加速推動溫室氣體減排、減緩全球變暖趨勢的呼聲日益高漲[1]。要開發(fā)一種新型的制冷劑,首先要從熱物性方面入手,氣液相平衡性質(zhì)是混合物最基本的熱力學性質(zhì),是其他熱力學性質(zhì)的研究基礎(chǔ)。國際上已經(jīng)發(fā)展多種氣液相平衡實驗方法,如循環(huán)法,靜態(tài)法,外推法等。

        我國清華大學侯樹鑫,段遠源研究員采用循環(huán)法建立高精度氣液相平衡實驗系統(tǒng)測定了R125+R290 在263.15 K,293.15 K,323.15 K下的性質(zhì)。實驗測定不僅實驗系統(tǒng)結(jié)構(gòu)復雜,實驗設(shè)備和操作環(huán)節(jié)較多,對于有毒性、易爆性或其他有害物質(zhì),用實驗方法測量熱力學性質(zhì)也很難行得通[2]。除實驗法外,通過狀態(tài)方程結(jié)合混合原則計算某些混合工質(zhì)熱物性。中國科學院理化技術(shù)研究所的公茂瓊等[3]采用nonrandom-two-liquid activity coefficient(NRTL)模型,Peng-Robinson(PR)狀態(tài)方程結(jié)合Huron-Vidal(HV)混合原則研究R1234ZE+R290的氣液相平衡性質(zhì)。張宇等[4]利用同樣的方法研究了二元混合制冷劑R170+R116,并將計算結(jié)果與實驗結(jié)果比對。而韓國首爾國立大學的Ju Hyok Kim等[5]將PR狀態(tài)方程與van-der-waols混合原則結(jié)合計算了環(huán)保型制冷劑CO2+R290在溫度253.15 K~323.15 K的氣液相平衡,法國的 Christophe Coquelet等[6]研究了 R32+R227ea混合制冷劑溫度從283.03 K~363.21 K的氣液相平衡。Hao Guo等[7]通過實驗測定與狀態(tài)方程結(jié)合相比較研究了混合制冷劑R152a+R134在溫度258.15 K~288.15 K的氣液相平衡。法國的F. Rivollet等[8]利用靜態(tài)分析法實驗測定CO2+R32在溫度283.12 K~343.25 K,壓力最高達7.46 MPa下的氣液相平衡。近年,因為層出不窮的新物質(zhì)使得狀態(tài)方程法越來越受限制,而計算機技術(shù)的發(fā)展,計算機模擬的優(yōu)勢也越來越明顯。趙勝喜等[9]利用吉布斯蒙特卡羅模擬了R410A的氣液相平衡并得到其臨界狀態(tài)點。意大利里雅斯特大學化學系的 Oliver Milocco等[10]列舉了部分制冷劑利用COSMO-RS與分子動力學模擬兩種方法預測熱物性。趙勝喜[11]等通過分子動力學模擬飽和液態(tài)制冷劑氨的熱物理性質(zhì)(比焓+密度)。相比較幾種模擬方法,COSMO-RS它繞開了勢能函數(shù)及其參數(shù)的選擇,僅依賴于單個分子的結(jié)構(gòu),因而具有較高的理論精度,是一種新穎的統(tǒng)計熱力學方法。

        本文基于以上現(xiàn)狀分析,結(jié)合前人已研究的成果,嘗試模擬分析了兩種新型二元混合制冷劑(R290 +R227ea)的氣液相平衡性質(zhì),采用由 Andreas Klamt[12]1995年提出的真實溶劑似導體屏蔽模型(Conductor like Screening Model for Realistic Solvents),簡稱COSMO-RS。是一種基于分子的量子化學計算信息結(jié)合統(tǒng)計力學理論進行物質(zhì)的熱物性計算的方法。該方法只需要物質(zhì)單個分子的結(jié)構(gòu)信息,就可以預測物質(zhì)的熱物性,是目前量子力學世界和工程熱力學聯(lián)系最有效的途徑。

        1 COSMO-RS模擬原理

        1.1 COSMO-RS概念

        COSMO是一種連續(xù)介質(zhì)溶劑化模式[13-15],模型中將連續(xù)介質(zhì)的介電常數(shù)設(shè)為無窮大(理想導體),這樣可以將屏蔽電荷限制在界面上,從而分子和溶劑間沒有電場,導體內(nèi)沒有電荷。在COSMO基礎(chǔ)上計算表面屏蔽電荷密度,再結(jié)合統(tǒng)計力學方法得到各組分的化學勢,最終計算獲得其他物理性質(zhì)如剩余焓、熵、蒸汽壓等。COSMO-RS通過對單個分子的量子化學計算來預測多元體系的相平衡[16]。下面介紹預測流體(R290+R227ea)氣液相平衡熱物性方面的應用。

        1.2 σ-profile

        分子表面的屏蔽電荷密度是采用COSMO-RS理論最重要的數(shù)據(jù)。因為我們要描述的真實體系都沒有真實的面,為了描述真實體系,我們必須引入一個虛擬的分子間相互作用的節(jié),并且這些節(jié)要有一定的表面積,每個節(jié)的表面屏蔽電荷能用一個合適的平均值進行表示,能滿足相互配對。我們可以采用量子化學方法得到這些節(jié)的表面電荷密度,得到一種類型分子的有效概率函數(shù)p(σ),稱為分子的“σ-profile”,它表示在一個相互接觸的分子節(jié)上面計算得到平均屏蔽電荷密度σ的概率,由量化計算產(chǎn)生,經(jīng)過一個基于下述公式的半徑γαν的面積上的系綜平均步驟:

        式中:dμυ為片段μ和υ之間的距離;γμ為根據(jù)片段面積計算得片段μ的平均半徑;γαυ為可調(diào)節(jié)參數(shù)。

        然后通過電荷平均得到σ-profile。分子i的σprofile定義為:

        式中:Ai(σ)為電荷密度為σ的所有片段的總表面積;Ai為總空穴表面積;ni(σ)為電荷密度為σ的片段數(shù)量;ni為總片段數(shù)。

        對于混合物,發(fā)現(xiàn)具有屏蔽電荷密度為σ的片段的幾率ps(σ),由體系中各種組分i的σ-profile,pi(σ)及其摩爾分數(shù)xi加權(quán)平均得到:

        2 COSMO-RS模擬步驟

        量化計算時可同時得到分子的表面積(molecular surface)Ai和總的空腔體積(cavity volume)Vi,上述信息可用于計算活度系數(shù)的組合部分;通過解自恰方程可得到片段活度系數(shù),混合物中物質(zhì)的活度系數(shù)的剩余部分可從片段活度系數(shù)得到。下面詳細介紹R290+R227ea的模擬步驟。

        2.1 σ-profile

        分子表面被分成若干面積為aeff的小塊,這些小塊相互之間是相互獨立的,每個小塊的表面屏蔽電荷是由量子化學中COSMO方法計算得到的。圖1,圖2分別是COSMO計算得到的R290,R227ea的σ-profile;圖3,圖4分別是直觀COSMO表面電荷密度。

        圖1 R290的σ-profileFig.1 σ-profile of R290

        2.2 COSMO-RS表面電荷的相互作用

        COSMO-RS理論的分子間相互作用力是根據(jù)分子表面片段的相互作用上得到的。對于理想的屏蔽分子,如果兩個節(jié)所帶的表面密度電荷一正一負,則它們的相互作用能肯定為零。實際上這種所有分子的理想配對是不可能的。對于實際流體,我們往往采用相互接觸的片段對的凈屏蔽電荷密度σ1和σ2來計量兩個節(jié)在真實體系和理想體系中的能量差別(misfit energy):

        圖2 R227ea的σ-profileFig.2 σ-profile of R227ea

        圖3 R290的COSMO表面電荷密度Fig.3 COSMO surface charge density of R290

        圖4 R227ea的COSMO表面電荷密度Fig.4 COSMO surface charge density of R227ea

        如果我們研究的體系是強極性物質(zhì),我們還應該考慮氫鍵的相互作用,Klamt和Eckert也提出了計算氫鍵相互作用的簡單表達式:

        式中:σdonor=min[σ1,σ2],σaccept=max[σ1,σ2]。CHB和σHB都是可調(diào)參數(shù)。

        另外,分子間的作用力,除了考慮Emisfit和EHB之外,還考慮了分子片段間的范德華相互作用力Evdw:

        上式中的τvdw和τ′vdw都是基本元素的可以調(diào)節(jié)的參數(shù)。范德華的作用能量與組成物質(zhì)的原子類型有關(guān)。

        總的分子相互作用能量可由下式得出:

        在得到分子間相互作用能量之后,經(jīng)過一系列嚴格的統(tǒng)計熱力學理論計算,可以得到系統(tǒng)中各組分的化學勢。在此基礎(chǔ)上,計算工程需要的熱物理性質(zhì)。

        2.3 二元混合物氣液相平衡數(shù)據(jù)

        由于大多數(shù)制冷工質(zhì)的壓力不高,我們可以把它們的氣相狀態(tài)看做理想狀態(tài),對于二元體系,由相平衡的熱力學關(guān)系可以得到:

        式中:psat1和psat2為純組分1和純組分2的飽和蒸汽壓;ptotal為系統(tǒng)的總壓力;x1和x2為組分1和組分2在液相的摩爾分數(shù);γ1和γ2為采用COSMO-RS方法預測的組分1和組分2的活度系數(shù)。

        根據(jù)道爾頓分壓力定律,組分1的分壓力可用下式表示:

        式中:y1為組分1在氣相中的摩爾分數(shù)。

        根據(jù)拉烏爾定律,組分1的分壓力也可用下式表示:

        聯(lián)立上面三式可得,氣液相平衡時組分1在汽相中的摩爾分數(shù):

        二元共沸體系中,共沸點處的壓力對組分的一階導數(shù)為0,即:

        定義COSMO-RS模擬值與實驗[15]數(shù)據(jù)相對誤

        差:

        式中:ps為COSMO-RS模擬壓力值,MPa;pe為實驗[17]壓力值,MPa;y1s為組分R290(1)氣相COSMO-RS模擬物質(zhì)的量,mol;y1e為組分R290(1)氣相實驗[17]物質(zhì)的量,mol;x1s為組分R290(1)液相COSMO-RS模擬物質(zhì)的量,mol;x1e為組分R290(1)液相實驗[17]物質(zhì)的量,mol。

        3 R290+R227ea氣液相平衡模擬結(jié)果

        HFO-227ea,1,1,1,2,3,3,3—·—七氟丙烷(R227ea)是國際上最近提出的一種新的不破壞臭氧層的替代物,已被作為滅火劑用于替代哈龍,并被認為是一種很有希望的制冷劑替代物,尤其是作為混合物的一種組分用于R502和R22的替代。R290,丙烷C3H8,分子量44.9,沸點-42.2℃,ODP=0,GWP=0.01,可用于替代R22,R502,是一種環(huán)保型制冷劑。基于國內(nèi)外行情,不僅對臭氧層空洞問題日益關(guān)注,溫室效應也越來越受到廣泛重視,丙烷以其極小的GWP以及零值ODP,低沸點的優(yōu)勢被研究者用來探索新型制冷劑。但是因為它易燃易爆,須添加阻燃劑降低其燃爆可能性,R227ea是理想的滅火劑,本文嘗試將這兩種制冷劑混合研究了定溫下(293.16 K,303.15 K,313.14 K,333.15 K,343.16 K,353.188 K)氣液相平衡時的總壓力p(MPa)隨組分R290(1)的變化,分析利用COSMO-RS模型預測制冷劑熱物性的可行性。用于替代現(xiàn)行個別的具有較高GWP的中溫制冷劑。模擬結(jié)果如下表1所示,由模擬結(jié)果可知,混合制冷劑在209.16 K~353.18 K下氣液相平衡時壓力的最大相對誤差在4.5%內(nèi),共沸點處的組分1摩爾分數(shù)最大相對誤差在4%??傮w看來,雖然有些偏差,但是從圖5,圖6可以看出與文獻[17]提供的實驗數(shù)據(jù)具有很好的一致性。另外,圖表說明該模型在低溫吻合的更好,如何改進高溫部分還值得進一步研究。

        實線:COSMO-RS模擬數(shù)值結(jié)果;散點:實驗[17]數(shù)據(jù)。根據(jù)圖表,找到不同溫度下,R290+R227ea的共沸點,結(jié)果如表2,圖7所示,與實驗[17]測得的值誤差在4%內(nèi),證明了該模擬方法在預測氣液相平衡具有一定的可行性。

        圖5 制冷劑(R290+R227ea)的氣液相平衡熱物理性質(zhì)Fig.5 Vapor-liquid equilibrium thermodynamic properties of refrigerant R290+R227ea

        圖6 制冷劑(R290+R227ea)的氣液相平衡壓力模擬與實驗值相對誤差Fig.6 Deviations of the pressure between the experiment data and the simulation of refrigerant(R290+R227ea)

        圖7 制冷劑(R290+R227ea)的共沸點模擬與實驗值相對誤差Fig.7 Deviations of the azeotropic composition x1and pressure p between the experiment data and the simulation of refrigerant(R290+R227ea)

        4 結(jié)論

        本文基于Kalmt等人提出量子力學理論,采用真實溶劑似導體屏蔽模型,模擬新型二元混合制冷劑(R290+R227ea)在293.16 K~353.18 K溫度下氣液相平衡的熱物理性質(zhì),模擬結(jié)果壓力及共沸點數(shù)據(jù)證實與實驗值相對誤差控制在4%以內(nèi),有很好的一致性,表明該模型可以計算混合制冷工質(zhì)的氣液相平衡性質(zhì)。

        表1 制冷劑(R290+R227ea)的氣液相平衡熱物理性質(zhì)Tab.1 Vapor-liquid equilibrium thermodynamic properties of refrigerant(R290+R227ea)

        333.15 K 0.5194 2.1312 2.1490 0.84 0.6214 0.6776 9.04 0.6331 2.2110 2.2141 0.14 0.6956 0.7268 4.49 0.7304 2.2450 2.2450 0.00 0.7596 0.7679 1.09 0.8283 2.2430 2.2520 0.40 0.8296 0.8178 -1.42 0.9316 2.1892 2.2100 0.95 0.9197 0.8987 -2.28 1.0000 2.1168 2.1168 0.00 1.0000 1.0000 0.00 343.16 K 0.0000 1.4874 1.4874 0.00 0.0000 0.0000 0.1637 1.9623 2.0331 3.61 0.2801 0.3792 35.38 0.3110 2.2891 2.3610 3.14 0.4327 0.5396 24.71 0.5108 2.5976 2.6220 0.94 0.5890 0.6649 12.89 0.6246 2.7060 2.7038 -0.08 0.6713 0.7170 6.81 0.7369 2.7579 2.7468 -0.40 0.7540 0.7671 1.74 0.8341 2.7505 2.7512 0.03 0.8303 0.8196 -1.29 0.9401 2.6709 2.6897 0.70 0.9285 0.9099 -2.00 1.0000 2.5868 2.5868 0.00 1.0000 1.0000 0.00 353.18 K 0.0000 1.8583 1.8583 0.00 0.0000 0.0000 0.1541 2.3607 2.4593 4.18 0.2374 0.3529 48.66 0.2343 2.5808 2.6925 4.33 0.3236 0.4552 40.67 0.3865 2.9276 3.0134 2.93 0.4565 0.5846 28.06 0.4511 3.0486 3.1089 1.98 0.5078 0.6244 22.96 0.5681 3.2264 3.2340 0.24 0.6011 0.6851 13.97 0.6538 3.3157 3.2934 -0.67 0.6729 0.7249 7.73 0.7653 3.3629 3.3317 -0.93 0.7690 0.7785 1.24 0.8119 3.3391 3.3390 0.00 0.8096 0.8048 -0.59 0.9450 3.2341 3.2430 0.28 0.9354 0.9167 -2.00 1.0000 3.1319 3.1319 0.00 1.0000 1.0000 0.00

        表2 組分1R290在共沸點的摩爾數(shù)及氣液相平衡壓力Tab.2 The azeotropic composition x1and pressure p for R290+R227ea at each temperature

        [1] 賈磊,史敏,張秀平,等.合成類HCFCs替代制冷劑的研究進展[J].制冷與空調(diào),2011,2(1):116-120.(Jia Lei,Shi Min,Zhang Xiuping,et al.The advanced developments of alternative refrigerants as anabolic HCFCs[J].Refrigeration and Air-condition,2011,2(1):116-120.)

        [2] 侯樹鑫,段遠源.高精度氣液相平衡實驗系統(tǒng)的建立和實驗測定[J].工程熱物理學報,2009,30(7):1093-1097.(Hou Shuxin,Duan Yuanyuan.The establishment and the experimental measurement of vapor-liquid equilibrium experimental system with high-precision[J].Journal of Engineering Thermophysics,2009,30(7):1093-1097.)

        [3] Xueqiang Dong,Maoqiong Gong,Jun Shen,et al.Experimental measurement of vapor-liquid equilibrium for(trans-1,3,3,3-tetrafluoropropene(R1234ze(E))+propane (R290))[J].International journal of refrigeration,2011,34:1238-1243.

        [4] Yu Zhang,MaoQiong Gong,HongBo Zhu,et al.Vaporliquid equilibrium measurements and correlations for an azeotropic system of ethane+hexafluoroethane[J].Fluid Phase Equilibria,2006,240:73-78.

        [5] Ju Hyok Kim,Min Soo Kim.Vapor-liquid equilibria for the carbon dioxide+propane system over a temperature range from 253.15 to 323.15 K[J].Fluid Phase Equilibria,2005,238:13-19.

        [6] Christophe Coquelet,Alain Valtz,Dominique Richon.Vapor-liquid equilibrium data for the difluoromethane(R32) +dimethylether(RE170)system at temperatures from 283.03 to 363.21 K and pressures up to 5.5 MPa[J]. Fluid Phase Equilibria,2005,232:44-49.

        [7] Hao Guo, Maoqiong Gong, XueqiangDong, etal. (Vapour+liquid)equilibrium data for the azeotropic{1,1-difluoroethane(R152a)+1,1,2,2-Tetrafluoroethane (R134)}system at various temperatures from(258.15 to 288.15)K[J].J.Chem.Thermodynamics,2012,54: 129-133.

        [8] F Rivollet,A Chapoy,C Coquelet,et al.Vapor-liquid equilibrium data for the carbon dioxide(CO2)+difluoromethane(R32)system at temperatures from283.12 to 343.25 K and pressures up to 7.46 MPa[J].Fluid Phase Equilibria,2004,218:95-101.

        [9] 趙勝喜,祁影霞.制冷劑R410a的氣液相平衡的吉布斯蒙特卡羅模擬[J].低溫與超導,2012,40(3):69-72. (Zhao Shengxi,Qi Yingxia.Gibbs Ensemble Monte Carlo simulation of vapor-liquid equilibrium properties of R410a [J].Cryogenics and Superconductivity,2012,40(3):69-72.)

        [10]Oliver Milocco,Maurizio Fermeglia,Sabrina Pricl.Prediction of thermophysical properties of alternative refrigerants by computational chemistry[J].Fluid Phase Equilibria,2002,199:15-21.

        [11]趙勝喜,祁影霞.飽和液態(tài)制冷劑氨的熱物理性質(zhì)的分子動力學模擬[J].制冷學報,2012,33(1):32-34. (Zhao Shengxi,Qi Yingxia.Molecular dynamics simulation of thermodynamic properties of the saturated liquid ammonia[J].Journal of Refrigeration,2012,33(1):32-34.)

        [12]Klamt A.Conductor-like screening model for real solvents: A new approach to the quantitative calculation of solvation phenomena[J].J.Chem.Phys.,1995,99:2224-2235.

        [13]Klanlt A,Jonas V,Burger T,et al.Refinement and Parameterization of COSMO-RS[J].The Journal of Physical Chemistry,1998,102(26):5074-5085.

        [14]Klamt A,Eckert F.COSMO-RS:a novel and efficient method for the a priori prediction of thermophysical data of liquids[J].Fluid Phase Equilibria,2000,172(1):43-72.

        [15]李勝迎.醇+(酮、離子液體)二元體系的過量焓測定、關(guān)聯(lián)和COSMO-type模型的應用[D].浙江大學理學院,2008.

        [16]Klamt A,Schuurmann G.COSMO:a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient[J].Journal of the Chemical Society,Perkin Transactions,1993,5:799-805.

        [17]A Valtz,C Coquelet,A Baba-Ahmed,et al.Vapor-liquid equilibrium data for the propane+1,1,1,2,3,3,3-heptafluoropropane(R227ea)system at temperatures from 293.16 to 353.18 K and pressures up to 3.4 MPa[J]. Fluid Phase Equilibria,2002,202:29-47.

        Vapor-liquid Equilibrium Properties of New Binary Mixture Refrigerant(R290+R227ea)

        Chen Xiuping Qi Yingxia Zhao Shengxi Zhang Hua

        (College of Energy and Power Engineering,University of Shanghai for Science and Technology,Shanghai,200093,China)

        The vapor-liquid equilibrium data is an important part of thermodynamic parameters of a new refrigerant.Theoretical predictions and experimental measurements of vapor-liquid equilibrium data of a new mixed working fluid become an urgent need for alternative refrigerant researches.The article is based on the theory of quantum chemistry.The vapor-liquid equilibrium properties of binary refrigerant(R290+R227ea)were modeled by COSMO-RS simulations.The simulation results could accord with the data of experiments well.It shows that COSMO-RS simulation method is feasible to predict the vapor-liquid equilibrium properties of mixture refrigerants.

        vapor-liquid equilibrium;quantum chemistry;COSMO-RS;R290+R227ea

        TB61+1;TK124

        A

        0253-4339(2014)05-0094-07

        10.3969/j.issn.0253-4339.2014.05.094

        陳秀萍,女(1989-),在讀碩士,上海理工大學能源與動力工程學院,18818251838,E-mail:shirley12141019@126.com。研究方向:COSMO-RS預測新型替代制冷工質(zhì)的熱物性。

        2013年11月22日

        About the author

        Chen Xiuping(1989-),female,Master Candidate,College of Energy and Power Engineering,University of Shanghai for Science and Technology,18818251838,E-mail:shirley12141019@126. com.Research fields:COSMO-RS studies of the thermodynamics properties of the new alternative refrigerant.

        猜你喜歡
        制冷劑物性電荷
        連續(xù)分布電荷體系電荷元的自能問題*
        物理通報(2024年4期)2024-04-09 12:41:28
        揚州市出臺“制冷劑六條”提升制冷劑安全管理水平
        R1234ze PVTx熱物性模擬計算
        能源工程(2022年1期)2022-03-29 01:06:26
        電荷知識知多少
        中韓天氣預報語篇的及物性分析
        LKP狀態(tài)方程在天然氣熱物性參數(shù)計算的應用
        煤氣與熱力(2021年6期)2021-07-28 07:21:30
        電荷守恒在化學解題中的應用
        低孔低滲儲層物性下限確定方法及其適用性
        近共沸制冷劑R134a/R1234yf的PVTx性質(zhì)的實驗研究
        靜電現(xiàn)象有什么用?
        在线观看高清视频一区二区三区| 国产亚洲精品资源在线26u| 午夜成人理论无码电影在线播放| 国产91 对白在线播放九色| 女人被躁到高潮嗷嗷叫免费软| 久久色悠悠综合网亚洲| 青青草精品在线视频观看| 亚洲视频在线观看| 国产精品va在线观看无码| 中文乱码字幕高清在线观看| 青青草最新在线视频观看| 91国产精品自拍在线观看| 日韩av高清在线观看| 麻豆一区二区99久久久久| 中国精品视频一区二区三区| 中文字幕久久熟女人妻av免费| 日韩精品一区二区免费 | 在线亚洲精品国产成人二区| 亚洲熟女熟妇另类中文| 久久婷婷国产综合精品| 精品欧洲av无码一区二区| 在线人妻无码一区二区 | 东方aⅴ免费观看久久av| 亚洲三级在线播放| 中文字幕一区二区三区6| 国产av剧情一区二区三区| 久久久久成人精品无码| 欧美黑人乱大交| 美女叉开双腿让男人插| 亚洲精品中字在线观看| 少妇愉情理伦片丰满丰满| 国产色a在线观看| 中文字幕大乳少妇| 亚洲精品综合中文字幕组合| 国产97色在线 | 国产| 台湾佬自拍偷区亚洲综合| 精品久久亚洲一级α| 男生自撸视频在线观看| 国产精品成人无码久久久久久| 日本理伦片午夜理伦片| 欧美成人高清手机在线视频|