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

HFO-1234yf/HFC-32氣-液相平衡特性研究

2016-04-25 06:10:51祁影霞
制冷技術(shù) 2016年1期
關(guān)鍵詞:方法

吳 東,祁影霞,楊 喜

(上海理工大學(xué),上海 200093)

?

HFO-1234yf/HFC-32氣-液相平衡特性研究

吳東*,祁影霞,楊喜

(上海理工大學(xué),上海 200093)

[摘 要]本文簡(jiǎn)單介紹了吉布斯蒙特卡洛(GEMC)模擬法原理并采用GEMC模擬了混合制冷劑HFO-1234yf/HFC-32質(zhì)量分?jǐn)?shù)分別為50%/50%和80%/20%,在溫度(273~333) K的氣液相平衡特性。本文計(jì)算結(jié)果與PR狀態(tài)方程計(jì)算值的比較表明,飽和氣液相密度相對(duì)偏差分別為-3%~+2%、-4.5%~+3.5%。結(jié)果表明,GEMC模擬方法誤差小,為進(jìn)一步研究混合制冷劑HFO-1234yf/HFC-32氣液相平衡特性奠定了基礎(chǔ)。

[關(guān)鍵字]模擬;HFO-1234yf/HFC-32;氣液相平衡;狀態(tài)方程

*吳東(1988-),男,碩士在讀,主要從事制冷劑熱物性研究。通信作者:祁影霞,女,副教授,地址:上海市楊浦區(qū)軍工路516號(hào)上海理工大學(xué)。E-mail:qipeggy@126.com。

0 引言

隨著CFCs和HCFCs類制冷劑逐步的禁止使用,國(guó)際上形成了兩種不同的制冷劑替代研究路線[1]:一條以美日為代表,主要采用合成HFCs類制冷劑,如R134a、R410A和R407C等;另外一條以德國(guó)和北歐一些國(guó)家為代表,主張使用天然制冷劑,如CO2、H2O和HCs等。當(dāng)前,汽車空調(diào)依然廣泛采用R134a作為制冷劑,主要是基于其ODP為0,安全性能高;根據(jù)歐盟F-gas法規(guī),自2011年起禁止使用GWP大于150的制冷劑,而其GWP為1,430[2-4]。陸岷山等[5]通過(guò)實(shí)驗(yàn)的方法對(duì)混合制冷劑R134a/R600a的PVT性質(zhì)進(jìn)行了實(shí)驗(yàn)研究。

從上可以看出,由于理想制冷劑應(yīng)符合熱力學(xué)性質(zhì)、遷移性質(zhì)、物理化學(xué)性質(zhì)和其他性質(zhì)等諸多方面的要求,因此新型混合制冷劑的研究將有助于找到更為理想的替代制冷劑[6-8]。近年,由于計(jì)算機(jī)技術(shù)的巨大進(jìn)步,建立在量子力學(xué)和統(tǒng)計(jì)熱力學(xué)基礎(chǔ)上的分子模擬方法已成為預(yù)測(cè)新型混合制冷劑的一項(xiàng)突破性技術(shù)[9]。當(dāng)前分子模擬方法主要有:蒙特卡洛(Monte Carlo,MC)、分子動(dòng)力學(xué)(Molecular Dynamics,MD)、分子力學(xué)(Molecular Mechanics,MM)和量子力學(xué)(Quantum Mechanics,QM)[10]。其中,蒙特卡洛法是隨機(jī)性的分子模擬方法,執(zhí)行了Gibbs的統(tǒng)計(jì)力學(xué),還將分子動(dòng)力學(xué)方法建立在經(jīng)典力學(xué)算符理論的基礎(chǔ)上,使分子模擬有了一個(gè)嚴(yán)格的理論依據(jù)。重慶大學(xué)蔣國(guó)柱等[11]應(yīng)用該方法對(duì)L-J流體的氣液相平衡開展了模擬研究。

混合制冷劑氣液相平衡研究的方法主要有模擬仿真、狀態(tài)方程和實(shí)驗(yàn)。目前應(yīng)用最廣泛的氣液相平衡模擬方法[12]就是吉布斯系統(tǒng)蒙特卡洛方法(GibbsEnsemble Monte Carlo,GEMC),通過(guò)分別模擬氣相和液相避免了對(duì)相界面的直接處理,對(duì)純工質(zhì)及混合工質(zhì)的氣液相平衡模擬結(jié)果與實(shí)驗(yàn)具有較高的重合度。2009年,RAABE等[13]開發(fā)了一種適用于預(yù)測(cè)R1234yf熱力學(xué)性質(zhì)的新型力場(chǎng),并用GEMC方法模擬了其氣液相平衡特性。目前,由于R32和R1234yf混合氣液相實(shí)驗(yàn)數(shù)據(jù)較少,本文分別采用GEMC法和PR狀態(tài)方程法對(duì)混合制冷劑R32/R1234yf進(jìn)行氣液相平衡模擬,并將MC模擬和狀態(tài)方程計(jì)算的結(jié)果進(jìn)行比較分析。

1 GEMC模擬

1.1 Monte Carlo基本原理

MC方法又稱為隨機(jī)取樣法,它以概率論中的大數(shù)定理和中心極限定理為基礎(chǔ),通過(guò)不斷產(chǎn)生隨機(jī)數(shù)序列來(lái)模擬物體的運(yùn)動(dòng)過(guò)程。MC模擬方法通過(guò)粒子的隨機(jī)移動(dòng)使構(gòu)型發(fā)生改變[14],進(jìn)行系綜上的平均,當(dāng)給初始構(gòu)型賦予體系后,對(duì)構(gòu)型內(nèi)每個(gè)粒子進(jìn)行隨機(jī)移動(dòng),比較移動(dòng)前后體系狀態(tài)的變化,按照一定的接受準(zhǔn)則對(duì)移動(dòng)進(jìn)行判斷,若移動(dòng)滿足這一準(zhǔn)則就接受,若不滿足就拒絕,重復(fù)移動(dòng)多次,就能得到體系的真實(shí)構(gòu)型,通過(guò)對(duì)系綜上的平均求出體系的宏觀性質(zhì)。圖1為MC模擬流程圖。

圖1 MC模擬流程圖

1.1.1 系綜

系綜是許多同樣結(jié)構(gòu)并處在給定宏觀條件之下的系統(tǒng)集合[15],其中每個(gè)系統(tǒng)各處于某一個(gè)微觀狀態(tài),并且彼此獨(dú)立;在相空間內(nèi)微觀運(yùn)動(dòng)狀態(tài)組成一個(gè)連續(xù)的區(qū)域,宏觀量與微觀量相對(duì)應(yīng),在一定的宏觀條件下是所有可能的運(yùn)動(dòng)狀態(tài)的平均值。在傳統(tǒng)MC模擬方法的應(yīng)用中,存在著兩種系統(tǒng)狀態(tài)的模擬:平衡態(tài)的模擬和非平衡態(tài)的模擬。在氣液相平衡的計(jì)算中,主要采用平衡態(tài)的模擬。平衡態(tài)系綜主要有以下幾個(gè)類型:微正則系綜(N、V、E)、正則系綜(N、V、T)、等溫等壓系綜(N、P、T)、巨正則系綜和吉布斯(Gibbs)系綜等。

1.1.2 力場(chǎng)

力場(chǎng)方法可以看成是勢(shì)能函數(shù)的經(jīng)驗(yàn)表達(dá)式,是MC模擬內(nèi)容的基礎(chǔ)。分子力場(chǎng)大致上可分為兩類:全原子力場(chǎng)和聯(lián)合原子力場(chǎng)。此外,也可根據(jù)選擇的函數(shù)和力場(chǎng)參數(shù)分為以下3類:傳統(tǒng)力場(chǎng),如AMBER、CHARMM、CVFF和MM等系列力場(chǎng);第二代力場(chǎng),如CFF91、PCFF、CFF95和MMF94等系列力場(chǎng);通用力場(chǎng),如ESFF、UFF和DREDING等系列力場(chǎng)。

1.2 GEMC方法計(jì)算原理

GEMC方法能夠同時(shí)在氣相和液相兩個(gè)模擬盒子中進(jìn)行MC模擬,并且兩者保持相對(duì)獨(dú)立,但是滿足相平衡條件(溫度T、壓力P和化學(xué)勢(shì)μ相等),總粒子數(shù)N、總體積V和溫度保持不變;為了能夠滿足相平衡條件,在模擬過(guò)程中需要進(jìn)行3種不同的MC移動(dòng),以不同的接受概率進(jìn)行判斷。

1)氣相和液相盒子中的粒子分別在盒子內(nèi)自由移動(dòng),如轉(zhuǎn)動(dòng)、平動(dòng)等,使盒子達(dá)到平衡狀態(tài),粒子移動(dòng)被接受的概率為:

2)在總體積不變的情況下,氣相和液相的體積發(fā)生變化,使兩個(gè)模擬盒子的壓力平衡,粒子移動(dòng)被接受的概率為:

3)在粒子總數(shù)N保持不變的情況下,在氣相和液相兩個(gè)模擬盒子中進(jìn)行交換,使兩個(gè)盒子的化學(xué)勢(shì)達(dá)到平衡狀態(tài),粒子移動(dòng)被接受的概率為:

式中:

E——盒子的能量;

V——體積;

N——粒子數(shù);

k——玻爾茲曼常數(shù)。

1.3 模擬方法

本文采用NPT+GEMC方法對(duì)R32/R1234yf進(jìn)行氣液相平衡的模擬,模擬建立氣相和液相兩個(gè)盒子,兩個(gè)盒子的初始結(jié)構(gòu)為面心立方(FCC)結(jié)構(gòu)。根據(jù)混合制冷劑的飽和蒸氣壓數(shù)據(jù),可以確定在質(zhì)量比20%/80%時(shí),液相分子數(shù)分別為20/44,氣相分子數(shù)分為為83/182;在質(zhì)量比25/75時(shí),液相分子數(shù)分別為23/40,氣相分子數(shù)分別為82/150。所有的Monte Carlo計(jì)算均采用MAPS(Materials and Process Simulation)的Towhee模塊進(jìn)行。GEMC模擬采用通用UFF力場(chǎng),電荷計(jì)算采用Gasteiger方法,系統(tǒng)總的電荷量為0,圖2為模擬所建立的氣相與液相盒子。

圖2 氣相與液相模擬盒子

模擬溫度范圍255 K~320 K,溫度間隔為3 K~4 K,初始設(shè)定的壓力P為溫度T對(duì)應(yīng)的飽和蒸汽壓,R32/R1234yf的飽和蒸氣壓數(shù)據(jù)采用KAMIAKA[16]的實(shí)驗(yàn)值。在模擬開始前需要對(duì)兩個(gè)盒子進(jìn)行優(yōu)化,優(yōu)化方法為最速下降法,運(yùn)行步數(shù)為1,000步,精度為0.01,van de Waals截?cái)嘁蜃訛?.5,花鍵移位因子為1.8。模擬循環(huán)的次數(shù)為10萬(wàn)次,前5萬(wàn)次循環(huán)用于使模擬體系達(dá)到氣液相平衡,后5萬(wàn)次循環(huán)用于統(tǒng)計(jì)氣液相平衡數(shù)據(jù)并取平均值。構(gòu)型偏倚公式采用Martin和Frischknecht(2006)形式,構(gòu)型偏倚設(shè)置按照Martin和Frischknecht方法,加入靜電作用項(xiàng),電介質(zhì)常數(shù)為1 Ang,靜電作用項(xiàng)的計(jì)算采用Ewald求和法。Van de Waals截?cái)喟霃饺橐合嗪虚L(zhǎng)的一半,截?cái)嗪蟛捎梦膊啃拚椒ā?/p>

模擬中Monte Carlo的移動(dòng)形式為:分子質(zhì)量中心的平移、質(zhì)量中心的轉(zhuǎn)動(dòng)、分子內(nèi)單原子的平動(dòng)、兩個(gè)盒子分子的構(gòu)型偏倚平動(dòng)、兩個(gè)盒子分子的構(gòu)型偏倚轉(zhuǎn)動(dòng)、各向異性體積移動(dòng)、各向同性體積移動(dòng)。移動(dòng)概率見表1。

表1 R32/R1234yf的MC移動(dòng)概率

1.4 模擬結(jié)果

不同配比的混合制冷劑模擬氣液相平衡密度如表2和表3所示。

由于缺乏R32/R1234yf的氣液相平衡密度的直接實(shí)驗(yàn)數(shù)據(jù),下節(jié)將與狀態(tài)方程在相同狀態(tài)下的計(jì)算結(jié)果進(jìn)行比較驗(yàn)證。

表2 R32/R1234yf(20/80,wt%)模擬氣液相平衡密度

表3 R32/R1234yf(50/50,wt%)模擬氣液相平衡密度

2 氣液相平衡特性計(jì)算及結(jié)果分析

目前,在NIST(National Institute of Standards and Technology)所發(fā)布的軟件Refprop9.0中還沒(méi)有R32/R1234yf的實(shí)驗(yàn)數(shù)據(jù),并且由于混合工質(zhì)的種類和配比繁多,實(shí)驗(yàn)測(cè)定工作量大,實(shí)驗(yàn)時(shí)間長(zhǎng),不可能一一加以實(shí)驗(yàn)研究。另外,根據(jù)KIM等[17]的研究可知,PR狀態(tài)方程計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)有很好的一致性。因此本文將用參數(shù)較少的PR方程對(duì)R1234yf/R32的氣液相平衡性質(zhì)進(jìn)行研究。

2.1 PR狀態(tài)方程介紹

PR方程形式如下所示[18]:

式(4)可寫成式(5)的形式:

式(5)根據(jù)系統(tǒng)中的相數(shù)會(huì)產(chǎn)生一個(gè)或兩個(gè)根,在兩相區(qū),最大的根是氣體的壓縮因子,最小的正根是液體的壓縮因子。式(4)應(yīng)用到臨界點(diǎn),有:

在非臨界點(diǎn)有:

式中:

Tr——對(duì)比溫度,T/Tc,

ω——偏心因子。

α(Tr,ω)的形式如下:

2.2 混合法則

一般來(lái)說(shuō),狀態(tài)方程只適用于純工質(zhì)計(jì)算,混合工質(zhì)則必須采用一定的混合法則以保證計(jì)算的精度。常用的混合法則有Van De Waals(VDW)、Huron-Vidal(HV)、Stryjek-Vera(SV)和MHV1等。其中,SV[19]混合法則是VDW的改進(jìn)形式,比VDW具有較高精確度。但是當(dāng)二元工質(zhì)混合時(shí),SV法則即可簡(jiǎn)化為VDW法則。因此,本文采用VDW混合法則計(jì)算混合工質(zhì)R32/R1234yf的氣液相平衡特性。

VDW混合法則的形式如下:

式中kij為二元相互作用系數(shù),xi為組分i的摩爾分?jǐn)?shù),ai和bi分別為組分i狀態(tài)方程參數(shù)。

在混合法則中,二元相互作用系數(shù)kij有著重要的影響。kij表示混合物的兩種分子間的非理想作用的特性和相互作用性質(zhì)[20],kij的準(zhǔn)確度對(duì)混合物氣液相平衡的計(jì)算至關(guān)重要,kij一般是由氣液相平衡數(shù)據(jù)優(yōu)化回歸得到。這里,R32/R1234yf的二元相互作用系數(shù)采用文獻(xiàn)[16]的實(shí)驗(yàn)數(shù)據(jù)擬合值0.037。

2.3 計(jì)算結(jié)果及分析

不同配比的混合二元制冷劑的飽和氣液相密度如表4和表5。由圖3和4可知,用GEMC模擬與PR狀態(tài)方程計(jì)算的R32/R1234yf飽和氣相密度相對(duì)偏差在-3.0%~+2.0%范圍內(nèi),相對(duì)偏差絕對(duì)值的平均偏差1.44%。飽和液相密度相對(duì)偏差在-4.5%~+3.5%范圍內(nèi),相對(duì)偏差絕對(duì)值的平均偏差2.25%。

表4 R32/R1234yf(20/80,wt%)飽和氣液相密度

表5 R32/R1234yf(50/50,wt%)飽和氣液相密度

圖3 R32/R1234yf的MC飽和氣相密度與PR數(shù)據(jù)的偏差

圖4 R32/R1234yf的MC飽和液相密度與PR數(shù)據(jù)的偏差

3 結(jié)論

GEMC模擬和PR狀態(tài)方程對(duì)于氣液相平衡性質(zhì)的計(jì)算都具有相當(dāng)高的精度,本文選擇了適當(dāng)?shù)膮?shù)對(duì)二元混合制冷劑R32/R1234yf分別進(jìn)行了氣液相平衡計(jì)算。結(jié)果表明,GEMC模擬方法誤差小,為進(jìn)一步研究混合制冷劑R32/R1234yf氣液相平衡特性奠定了基礎(chǔ)。

參考文獻(xiàn):

[1]陳斌,陳光明.R407C、R410A制冷系統(tǒng)相關(guān)特性研究進(jìn)展[J].制冷與空調(diào),2003,31(2):42-45.

[2]汪軍,李美玲,周文濤.核電站空調(diào)系統(tǒng)HCFC22冷水機(jī)組替代工質(zhì)性質(zhì)的研究[J].上海理工大學(xué)學(xué)報(bào),2001,23(1):8-13.

[3]葛志松,吳獻(xiàn)忠,崔曉鈺,等.三元非共沸混合工質(zhì)變組成容量控制空調(diào)系統(tǒng)可行性分析[J].上海理工大學(xué)學(xué)報(bào),2005,27(4):336-340.

[4]LONGO G A,ZILIO C,RIGHETTI G.Condensation of the low GWP refrigerant HFC152a inside a brazed plate heat exchanger[J].Experimental Thermal and Fluid Science,2015,68:509-515.

[5]陸岷山,祁影霞,楊喜,等.混合工質(zhì)R134a/R600a 氣相PVTx性質(zhì)的實(shí)驗(yàn)研究[J].制冷技術(shù),2014,34(6):27-31,36.

[6]鄭賢德.制冷原理與裝置[M].北京:機(jī)械工業(yè)出版社,2010.

[7]劉圣春,饒志明,楊旭凱,等.新型制冷劑R1234yf的性能分析[J].制冷技術(shù),2013,33(1):56-59.

[8]肖學(xué)智,周曉芳,徐浩陽(yáng),等.低GWP制冷劑研究現(xiàn)狀綜述[J].制冷技術(shù),2014,34(6):37-42.

[9]GLASS C W,REISER S,RUTKAI G,et al.ms2:A molecular simulation tool for thermodynamic properties,new version release[J].Computer Physics Communications,2014,185(12):3302-3306.

[10]莊昌青,岳紅,張慧軍.分子模擬方法及模擬軟件Material Studio在高分子材料中的應(yīng)用[J].塑料,2010,39(4):81-84.

[11]蔣國(guó)柱,薛榮書,魏順安.吉布斯系綜Monte Carlo技術(shù)模擬L-J流體汽液相平衡[J].重慶大學(xué)學(xué)報(bào),2003,26(8):64-67.

[12]李曉峰,趙立峰,孫淮.GEMC和GDI方法計(jì)算流體氣液相平衡的比較[J].物理化學(xué)學(xué)報(bào),2008,24(10):1824-1830.

[13]RAABE G,MAGINN E J.Molecular modeling of the vapor-liquid equilibrium properties of the alternative refrigerant 2,3,3,3-tetrafluoro-1-propene (HFO-1234yf)[J].Journal of Physical Chemistry Letters,2010,1(1):93-96.

[14]胡慧娟.吉布斯系綜模擬氮?dú)?水和氧氣-水的氣液相平衡[D].南京:南京工業(yè)大學(xué),2005.

[15]陳舜麟.計(jì)算材料科學(xué)[M].北京:化學(xué)工業(yè)出版社,2005.

[16]KAMIAKA T,DANG C,HIHARA E.Vapor-liquid equilibrium measurements for binary mixtures of R1234yf with R32,R125,and R134a[J].International Journal of Refrigeration,2013,36(3):965-971.

[17]JU H K,MIN S K,KIM Y.Vapor-liquid equilibria for pentafluoroethane + propane and difluoromethane + propane systems over a temperature range from 253.15 to 323.15 K[J].Fluid Phase Equilibria,2003,211(2):273-287.

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

[19]STRYJEK R,VERA J H.PRSV2:A cubic equation of state for accurate vapor-liquid equilibria calculations[J].The Canadian Journal of Chemical Engineering,1986,64:820-826.

[20]胡芃,高靜軒,陳龍祥,等.PR方程結(jié)合vdw混合法則推算二元及三元HFC/HC混合工質(zhì)氣液相平衡性質(zhì)[J].化工學(xué)報(bào),2012,63(2):350-355.

Investigation on Vapor-liquid Equilibrium Characteristics of HFO-1234yf/HFC-32

WU Dong*,QI Ying-xia,YANG Xi
(University of Shanghai for Science and Technology,Shanghai,200093)

[Abstract]Gibbs Ensemble Monte Carlo (GEMC) simulation is briefly introduced and adopted to model the vapor-liquid equilibrium characteristics of binary refrigerant HFO-1234yf/HFC-32 with mass fractions of 50%/50% and 80%/20% in the range of (273~333) K.The comparison between the simulated values and those calculated by PR state equations shows that,the relative deviations of the saturated vapor-liquid densites are –3%~+2% and –4.5%~+3.5%,respectively.It shows that,the deviation of the GEMC simulation method is small,which provides the basis for further researching the liquid-vapor equilibrium characteristics for the mixed refrigerant of HFO-1234yf/HFC-32.

[Keywords]Simulation; HFO-1234yf/HFC-32; Vapor-liquid equilibrium; Equation of state

基金項(xiàng)目:上海市教育委員會(huì)科研創(chuàng)新項(xiàng)目(No.11YZ119)。

doi:10.3969/j.issn.2095-4468.2016.01.106

猜你喜歡
方法
中醫(yī)特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數(shù)學(xué)教學(xué)改革的方法
化學(xué)反應(yīng)多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學(xué)習(xí)方法
用對(duì)方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡(jiǎn)單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 亚洲三级片在线看| 国产极品嫩模在线观看91| 午夜一级做a爰片久久毛片| 国产欧美视频在线观看| 粗大猛烈进出高潮视频无码| 波多野结衣一二三| 久久福利片| 日韩一二三区视频精品| 国产簧片免费在线播放| 日本黄色不卡视频| 久久黄色一级片| 久久久精品久久久久三级| 欧美激情伊人| 国产又爽又黄无遮挡免费观看| 国产性精品| 伊人中文网| 国产精品视频第一专区| 国产亚洲精品资源在线26u| 亚洲无码免费黄色网址| 国产网友愉拍精品视频| 亚洲欧美日韩成人高清在线一区| 伊人色在线视频| 亚洲天堂啪啪| 国产在线精品99一区不卡| 极品国产一区二区三区| 日本色综合网| 99精品热视频这里只有精品7| 成人亚洲视频| 中文字幕 91| 99久久精品国产麻豆婷婷| 久久99热66这里只有精品一| 久草性视频| 亚洲色图欧美在线| 国产玖玖玖精品视频| 国产成人高清精品免费5388| 亚洲国产综合精品中文第一| 精品人妻一区二区三区蜜桃AⅤ | 欧美福利在线| 久久性视频| 呦女亚洲一区精品| 大学生久久香蕉国产线观看 | 福利国产在线| 国产免费网址| 国产成人无码播放| 四虎综合网| 欧美成人看片一区二区三区| 极品国产一区二区三区| 为你提供最新久久精品久久综合| 久久久亚洲国产美女国产盗摄| 国产精品对白刺激| 天堂av高清一区二区三区| 草逼视频国产| 69视频国产| 亚洲综合欧美在线一区在线播放| 乱人伦99久久| 伊人精品视频免费在线| 亚洲最黄视频| 亚洲黄网视频| 国产麻豆精品手机在线观看| 91精品国产一区自在线拍| 无码日韩精品91超碰| 久久这里只精品国产99热8| 亚洲一区二区日韩欧美gif| 亚洲a免费| 22sihu国产精品视频影视资讯| 中文字幕乱妇无码AV在线| 国内精品一区二区在线观看 | 四虎AV麻豆| 99草精品视频| 老司机aⅴ在线精品导航| 国产成人三级在线观看视频| 国产成a人片在线播放| 亚洲欧美成人综合| 亚洲二三区| 香蕉蕉亚亚洲aav综合| 二级特黄绝大片免费视频大片| 国产91小视频在线观看| 日本道中文字幕久久一区| 久久精品最新免费国产成人| 久久96热在精品国产高清| 久草国产在线观看| 中国一级特黄大片在线观看|