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

        ?

        M-MOF-74吸附分離H2/He混合物的分子模擬研究

        2022-11-13 07:39:56王玉杰李申輝趙之平
        化工學(xué)報 2022年10期
        關(guān)鍵詞:環(huán)境壓力混合物選擇性

        王玉杰,李申輝,趙之平

        (1中國石化北京化工研究院,北京 100013; 2 北京理工大學(xué)化學(xué)與化工學(xué)院,北京 102488)

        引 言

        He 由于其化學(xué)惰性、分子量小、低沸點以及高導(dǎo)熱性等特點,廣泛用于航空航天、低溫超導(dǎo)、醫(yī)療等 領(lǐng) 域[1]。從 液 化 天 然 氣(liquefied natural gas,LNG)的不凝尾氣中提取He,是He 資源化利用的新興領(lǐng)域[2]。LNG 工廠的泡點氣因富集了He、Ar 等稀有氣體,是優(yōu)良的稀有氣體原料,通過對泡點氣進行分離來生產(chǎn)He 具有較好的發(fā)展前景[3]。此外,泡點氣中還含有大量H2、N2、CH4等雜質(zhì)氣體,對于N2、CH4等雜質(zhì)氣體可通過高壓冷凝等辦法來除去[4],但H2由于難以液化被濃縮在粗He 中, 其中H2的摩爾分數(shù)在7%~40%不等[2],并需要在粗He 精制前將H2除去。工業(yè)上一般采用催化氧化脫氫法來去除H2[5],但化學(xué)方法往往伴隨著催化劑的失活、反應(yīng)劑的消耗和H2的不可再生等問題。因此,采用物理的膜分離及吸附分離等方法來分離H2是一種更為綠色環(huán)保的手段。

        對于膜分離而言,Antunes 等[6]在壓力為7 bar(1 bar=105Pa)的條件下,研究了MFI-ZSM 沸石分子篩膜對H2∕He 混合物的分離能力,結(jié)果表明純H2及純He 的通量比僅為1.88 ± 0.03。Simplicio 等[7]通過MFI 型沸石分子篩中空纖維膜分離He 摩爾分數(shù)約為99.9%的H2∕He 混合物,依據(jù)環(huán)境溫度的差異,純H2與純He 的通量比最高可達2.31。Favvas 等[8]研究了活性炭中空纖維膜對H2、He、CO2、O2和CH45種氣體兩兩組合所形成混合物的分離能力,并指出隨著氣體分子分子動力學(xué)直徑的上升,氣體通量有下降趨勢,且純H2及純He 的通量比接近于1,這表明H2及He由于其相近的分子動力學(xué)直徑,較難依靠尺寸篩分的方法進行膜分離。

        然而,依據(jù)H2及He 在極性上的差異,理論上能夠通過兩者與吸附材料相互作用強度的差異來進行吸附分離。鑒于金屬有機框架材料(metal organic frameworks,MOFs)相較于無機分子篩與氣體分子具有更高的親和性[9],MOFs 可作為吸附分離材料的選擇對象。MOFs 材料上的有機配體通常以π-電子的弱范德華相互作用來捕獲H2[10],這使得MOFs 對H2的負載容量較低,達不到美國能源部(US Department of Energy,DOE)2020 年設(shè)置的目標(biāo)(5.5%(質(zhì)量)或者40 g·L-1)[11]。而Pham 等[12]的研究工作表明,含有金屬空配位點(open metal site,OMS)的MOFs由于H2與金屬的配位作用能夠顯著增加其對H2的吸附容量。因此,為了提高MOFs 對H2的親和性,可通過對MOFs 的去溶劑化處理,來暴露MOFs 的OMS。在眾多的包含OMS 的MOF 材料中,MOF-74 由于易于合成、孔徑及孔隙率大以及金屬中心離子易于調(diào)控等優(yōu)點,成為研究吸附分離混合氣體的熱點材料[13-16]。但通過實驗的方法研究H2∕He 混合物在MOF-74 上的吸附性能往往伴隨著人力和資源的消耗。隨著信息技術(shù)的發(fā)展,分子模擬技術(shù)已成為一種甄別MOFs 材料對特定體系是否具有分離性能,并能夠在分子層面的空間尺度解釋吸附機理的手段[17-19]。對于H2∕He 混合物的難分離問題,通過物理吸附法對H2∕He 進行分離的報道較少,這使得設(shè)計吸附H2∕He的MOFs材料缺乏先驗知識。因此,采用分子模擬的方法能夠通過極小的成本來獲取較多的信息,從而初步判斷某種MOFs 材料是否具有H2∕He分離的潛力。

        本文通過分子模擬的方法,首先,在不同的環(huán)境壓力下計算H2及He 在M-MOF-74(M= Mg、Co、Ni、Cu、Zn)的吸附量,來考察M-MOF-74 是否具有工業(yè)應(yīng)用的潛力。隨后,研究不同濃度的H2∕He 混合物在M-MOF-74 上H2和He 的競爭吸附行為。最后,通過探究H2和He 在M-MOF-74 上的吸附位點、吸附熱,以及不同金屬離子對H2的吸附強度,從而解析H2和He在M-MOF-74上的吸附機理。

        1 模擬方法

        M-MOF-74 的結(jié)構(gòu)文件從劍橋大學(xué)晶體數(shù)據(jù)庫(The Cambridge Crystallographic Data Centre,CCDC)下載獲得[20],并通過Chung 等[21]提供的Python腳本去除掉M-MOF-74 中吸附的溶劑分子。MMOF-74的結(jié)構(gòu)特點如表1所示,其結(jié)構(gòu)性質(zhì)由多孔材料分析軟件Zeo++計算得到[22]。通過GCMC 的方法來模擬M-MOF-74 對H2、He 以及H2∕He 混合物的吸附性能,且所有的GCMC 計算皆由RASPA(2.0 版本)軟件包來完成[23]。Monte Carlo 總的計算步數(shù)為5×107步,前2×107步用于平衡,后3×107步用于采樣。在模擬過程中保持MOFs 框架為剛性,且在每一步MC的模擬過程中,客體分子進入∕離開MOFs的概率以及在MOFs 中平移∕旋轉(zhuǎn)的概率皆為0.5。非鍵相互作用的截斷半徑為1.2 nm,且為了避免由周期性邊界條件(periodic boundary conditions,PBC)引起的最近鄰鏡像相互作用,將M-MOF-74 的單胞擴展為1×1×4 的超胞。采用TraPPE 力場描述H2和He 的非鍵相互作用勢[24],并采用五點模型來描述H2的分子結(jié)構(gòu),從而增加計算結(jié)果的精確性[25]。采用UFF 力場來描述M-MOF-74[26],原子電荷來自Pham 等[27]的工作。

        表1 M-MOF-74的結(jié)構(gòu)特點Table 1 Structural features of M-MOF-74

        M-MOF-74 的結(jié)構(gòu)及孔尺寸如圖1(a)所示,H2及He 的三維尺寸如圖1(b)所示,H2的五點模型如圖1(c)所示,H2與He 的物化性質(zhì)如表2 所示。MMOF-74 的力場參數(shù)及原子電荷如表3 所示。計算過程中,環(huán)境溫度維持為300 K。表4為本文計算得到的H2在M-MOF-74上吸附熱與實驗值的對比,結(jié)果表明本文使用的力場對Mg-MOF-74和Co-MOF-74具有較好的描述,分別高估和低估了Zn-MOF-74和Ni-MOF-74 與H2的吸附熱,但總體來說與實驗值具有相似的趨勢。此外,圖1(d)為M-MOF-74 在77 K 下對H2的吸附等溫線,其中的實驗數(shù)據(jù)來自文獻[36-38]。結(jié)果表明,模擬與實驗在不同壓力下的數(shù)值較為吻合,也說明了力場參數(shù)的合理性。綜合上述吸附熱及吸附等溫線的結(jié)果,能夠通過本文的模型定性對比不同金屬中心MOF-74 對H2∕He 混合物吸附分離能力。

        表2 H2與He的物化性質(zhì)Table 2 Physicochemical properties of H2 and He

        表3 M-MOF-74的力場參數(shù)及原子電荷Table 3 Force field parameters and atomic charges of M-MOF-74

        表4 本文計算得到的吸附熱與實驗值的對比Table 4 Comparison of heat of adsorption calculated in this work with experimental values

        圖1 (a)M-MOF-74的結(jié)構(gòu);(b)H2和He的分子尺寸;(c)H2的五點模型;(d)M-MOF-74對H2的吸附等溫線(其中Exp代表實驗測得的吸附等溫線,數(shù)據(jù)來自文獻[36-38];Sim代表本文模擬得到的數(shù)據(jù);1 ?=0.1 nm)Fig.1 (a)Structure of M-MOF-74;(b)Molecular dimensions of hydrogen and helium;(c)Five-point model of hydrogen;(d)Adsorption isotherms of H2 in M-MOF-74(where Exp represents the experimentally measured adsorption isotherm,and the data are from Refs.[36-38];Sim represents the data get from the simulation of this work)

        相較于分子力場方法,基于密度泛函理論(density functional theory, DFT)的相關(guān)方法在描述分子間相互作用時更為準(zhǔn)確。故本文將H2和He 在M-MOF-74上由GCMC方法得到的能量最低的吸附構(gòu)象作為DFT 的輸入構(gòu)象。為了研究OMS 位點對H2和He 的吸附作用,截取了包括H2和He 分子在內(nèi)的M-MOF-74 的一個金屬團簇作為研究對象。采用Gaussian(版本09)軟件,在泛函為M062X,基組為6-31++G(d, p)的級別下對所選團簇進行結(jié)構(gòu)優(yōu)化,且在優(yōu)化過程中,凍結(jié)M-MOF-74 的原子坐標(biāo),并依據(jù)式(1)計算H2和He與金屬中心的結(jié)合能:

        式中,EBindingEnergy為H2和He與M-MOF-74的結(jié)合能,kJ·mol-1;E(MOF-74+gas),EMOF-74及Egas分別為M-MOF-74 和H2或者He 同時存在時的能量,M-MOF-74 單獨存在時的能量及H2或者He 單獨存在時的能量,kJ·mol-1。隨后采用約化密度梯度(reduced density gradient, RDG)的方法[39]可視化了H2與M-MOF-74的弱相互作用。波函數(shù)分析軟件和可視化軟件分別為Multifwn[40](3.8版本)及VMD[41](1.9.3版本)。

        2 結(jié)果與討論

        2.1 環(huán)境壓力對吸附量的影響

        在給定體積與溫度的條件下,環(huán)境壓力越小氣相流體物質(zhì)的量越小,即壓力越低,氣相流體與吸附劑接觸的概率越低,處于動態(tài)平衡的吸附量也越低。因此,本文研究了環(huán)境壓力分別為1、5、10、15及20 bar 時純H2和純He 在M-MOF-74 上的吸附量,來考察M-MOF-74 是否在工業(yè)上具有應(yīng)用前景??紤]到不同金屬中心在摩爾質(zhì)量上有較大差異,吸附量的計算公式如式(2)所示

        式中,n為吸附H2或者He 的吸附量,mmol;V為M-MOF-74的體積,cm-3。由于不同金屬中心MOF-74 的體積接近,因此吸附量完全由金屬中心的差異來控制。如圖2(a)~(e)所示,隨著環(huán)境壓力的上升,H2和He 在M-MOF-74 上的吸附量皆上升,且H2的吸附量遠高于He 的吸附量。其中,H2和He 在Mg-MOF-74 上的吸附量最大,在20 bar 的環(huán)境壓力下Mg-MOF-74 對H2和He 的吸附量分別為3.33 和0.583 mmol·cm-3。而H2在Cu-MOF-74 上的吸附量最小,He 在Ni-MOF-74 上的吸附量最小。為了進一步對比H2和He 在不同金屬離子MOF-74 上吸附選擇性的差異,計算了不同環(huán)境壓力下純H2對純He 在M-MOF-74 上的吸附選擇性,計算式如式(3)所示

        圖2 環(huán)境壓力對H2及He在M-MOF-74上吸附量和選擇性的影響Fig.2 Effects of ambient pressure on the adsorption of H2 and He on M-MOF-74

        式中,S為選擇性;UH2、UHe分別為H2、He 在MMOF-74 上的吸附量。如圖2(f)所示,隨著環(huán)境壓力上升,純H2對純He 在M-MOF-74 上的吸附選擇性下降,即吸附量和選擇性存在著此消彼長的相互制約(trade-off)現(xiàn)象。其中,H2對He 在Ni-MOF-74上的吸附選擇性在環(huán)境壓力為1 bar時最高,大小為6.58,其次是Mg-MOF-74 的6.47,最小為Cu-MOF-74 的5.32。此外,對于不同的金屬中心離子的MOF-74,吸附選擇性的波動范圍皆在1.5 以內(nèi),說明MOF-74 金屬中心離子的變化不能顯著影響H2和He在M-MOF-74上的吸附性能。

        2.2 H2/He混合濃度對吸附性能的影響

        當(dāng)環(huán)境溫度、壓力及體積恒定時,氣相流體混合物某一組分的物質(zhì)的量與其密度相關(guān),且多組分在吸附劑表面存在競爭關(guān)系,而某組分濃度越大,其與吸附劑接觸的概率越大,使得該組分能夠在競爭吸附中取得優(yōu)勢。因此,通過計算H2摩爾分數(shù)分別為1%、5%、10%、20%以及40%的H2∕He混合物中H2對He 的分離因子,來修正He 由于較高的濃度所帶來的吸附競爭優(yōu)勢。分離因子的定義式如下

        式中,α為分離因子;UH2、UHe分別為H2、He 在M-MOF-74 上的吸附量;XH2為H2的摩爾分數(shù)。如圖3(a)~(e)所示,隨著H2濃度的上升,H2在MMOF-74 上的吸附量上升,而He 在M-MOF-74 上的吸附量下降,且H2吸附量上升的速度大于He 下降的速度。在H2摩爾分數(shù)為40%時,H2在Mg-MOF-74 上的吸附量最大,在Cu-MOF-74 上的吸附量最小。圖3(f)為H2摩爾分數(shù)對分離因子的影響,當(dāng)H2摩爾分數(shù)為40%時,H2∕He 混合物在Ni-MOF-74 上吸附時取得最高的分離因子,其值為6.78,而在Cu-MOF-74 上取得最低的分離因子,其值為5.48。此外,隨著H2摩爾分數(shù)的上升,不同金屬中心離子的MOF-74的分離因子皆無明顯變化,這表明H2∕He混合物中H2濃度的下降并不會減少H2與吸附劑相互作用的能力。

        圖3 H2∕He混合濃度對兩者在M-MOF-74上吸附量和分離因子的影響Fig.3 The effect of H2∕He mixture concentration on the adsorption of H2 and He on M-MOF-74

        2.3 M-MOF-74中溶劑分子對吸附量的影響

        為了定性描述OMS 的存在對H2及He 在MMOF-74 上吸附量的貢獻,計算了在環(huán)境壓力分別為1 和20 bar 時,M-MOF-74 含有OMS 以及OMS 被溶劑分子占據(jù)的情況下,純H2和純He 在M-MOF-74上的吸附量。如圖4(a)所示,OMS位點的存在對H2吸附量的影響遠大于對He 吸附量的影響。這提示OMS 為H2的主要吸附位點,而不吸附He。其中,H2吸附量受到OMS 影響最大的是Co-MOF-74,去除溶劑前后的吸附量分別為0.037 和0.178 mmol·cm-3,Ni-MOF-74 與之相近。而H2吸附量受到OMS影響最小的是Cu-MOF-74,去除溶劑前后的吸附量分別為0.063 和0.154 mmol·cm-3。當(dāng)環(huán)境壓力為20 bar 時,OMS 存在也會顯著影響H2在M-MOF-74 上的吸附量,但對比圖4(a)和圖4(b)可知,高壓情況下OMS 的存在對H2吸附量的影響要小于低壓情況下的影響。這主要是由于壓力較低時,氣體的物質(zhì)的量也較低,因此氣體分子會優(yōu)先占據(jù)相互作用較強的吸附位點,即在較低壓力時,H2分子主要占據(jù)OMS 位點,而隨著壓力的上升,其他位點也會逐漸被H2分子所占據(jù),從而稀釋了OMS位點對吸附量的貢獻。

        圖4 M-MOF-74上OMS位點對H2和He在M-MOF-74上吸附量及H2∕He選擇性的影響Fig.4 Effects of OMS sites on M-MOF-74 on the adsorption capacity of H2 and He and ideal selectivity of H2∕He on M-MOF-74

        2.4 吸附機理的分析

        2.4.1 吸附位點和能量分析 為了探究H2和He 在M-MOF-74 上吸附性能的差異,固定M-MOF-74 僅能吸附一個H2分子或He分子,隨后在溫度300 K 下進行4×107步的GCMC計算,前2×107步用于平衡,后2×107步用于采樣。由于為單分子吸附,其相互作用能即為該溫度下的吸附熱[42]。通過自編程的方法每隔2000步記錄H2分子或He分子的質(zhì)心位置和吸附熱,從而生成10000 個能夠關(guān)聯(lián)質(zhì)心位置和吸附熱的樣本,通過VMD(1.9.3版本)軟件對吸附熱的分布進行可視化分析。圖5 為H2和He 在M-MOF-74 上的吸附熱分布,其中顏色越接近藍色,代表吸附能力越弱,而顏色越接近紅色,代表吸附能力越強。由圖5(a)可知,H2的吸附主要發(fā)生在兩個位置,其一為OMS 位點,其二為有機配體的苯環(huán)附近,且前者的吸附熱為紅色,而后者為綠色,這表明H2在OMS 位點處的吸附熱高于苯環(huán)位點。而He 僅能夠吸附在有機配體的苯環(huán)附近,且吸附熱遠小于H2。此外,從兩個吸附位點附近的樣本數(shù)量可知,OMS位點所能夠容納的H2分子數(shù)量要遠小于有機配體苯環(huán)位點。然而,當(dāng)OMS 位點被溶劑分子占據(jù)時,從圖5(b)可看出,當(dāng)H2聚集在溶劑水分子周圍時,其吸附熱與吸附在苯環(huán)上的吸附熱接近,表明此時OMS不再是H2的吸附位點。

        圖5 H2、He在M-MOF-74上的吸附熱分布Fig.5 Adsorption heat map of H2 and He in M-MOF-74

        進一步將上述的10000個關(guān)聯(lián)質(zhì)心位置和吸附熱的樣本點進行統(tǒng)計分析,得到如圖6 所示的H2和He 在M-MOF-74 上的吸附熱頻率直方圖。從圖6(a)~(e)可看出,OMS 是否存在對He 在M-MOF-74上的吸附熱分布沒有較大影響,但會使得H2在MMOF-74 上的分布更加向低能量方向遷移,這意味著OMS主要以高能量的方式對吸附熱做貢獻。

        圖6 H2、He在M-MOF-74上的吸附熱頻率直方圖Fig.6 Histograms of adsorption heat frequencies of H2 and He on M-MOF-74

        2.4.2 金屬中心對吸附性能的影響 通過吸附熱分布圖和吸附熱頻率直方圖可知,不同金屬中心對H2和He 的吸附性能存在一定的差異。進一步通過DFT的方法,在電子大小的時空層面考察了H2和He與金屬中心相互作用的差異。如圖7(a)所示的RDG 分析,H2與M-MOF-74 皆存在著弱相互作用,表明M-MOF-74 的金屬位點的確為H2的一個吸附位點,這與上述GCMC得到的結(jié)果是一致的,且依據(jù)RDG 的大小和顏色差異可知相互作用強度Ni2+>Co2+>Mg2+>Zn2+>Cu2+。如圖7(c)所示,H2到MOF-74 金屬中心離子的距離依次為Cu2+(0.277 nm)>Zn2+(0.262 nm)>Mg2+(0.243 nm)>Co2+(0.217 nm)>Ni2+(0.211 nm),而He到MOF-74金屬中心離子的距離依次為Ni2+(0.342 nm)>Zn2+(0.335 nm)>Cu2+(0.333 nm)>Co2+(0.331 nm)>Mg2+(0.320 nm)。H2到金屬中心的距離越近,越易通過色散產(chǎn)生的負電中心與金屬離子進行相互作用,從而通過弱配位作用增強H2在M-MOF-74 金屬中心上的吸附。而He到所有金屬中心的距離皆大于0.3 nm,這使得He和金屬中心很難產(chǎn)生弱相互作用,這是因為He具有明顯的化學(xué)惰性,色散偶極矩極其微弱,無法與金屬中心及圖7(b)所示的芳香環(huán)位點產(chǎn)生明顯的弱相互作用。此外,He 到五種金屬中心的距離差異較小,這意味著五種金屬中心對He的吸附能力無明顯差異。圖7(d)為結(jié)合能的計算結(jié)果,H2與Ni2+的結(jié)合能絕對值最大為14.23 kJ·mol-1,H2與Cu2+的結(jié)合能絕對值最小為10.46 kJ·mol-1;He與所有金屬中心的結(jié)合能絕對值皆在3.70 kJ·mol-1左右,該結(jié)果進一步驗證了上述觀點。

        圖7 (a)H2在M-MOF-74上的金屬離子吸附位點及RDG分析(其中藍色代表一定的配位作用,綠色代表弱的范德華相互作用),等值面為0.2;(b)H2在M-MOF-74上的苯環(huán)吸附位點;(c)H2和He到M-MOF-74的平衡距離;(d)H2和He與M-MOF-74的結(jié)合能Fig.7 (a)Metal ion adsorption sites and RDG analysis of H2 on M-MOF-74(Blue represents certain coordination and green represents weak van der Waals interactions),isosurface at 0.2;(b)H2 adsorption sites on the benzene ring on M-MOF-74;(c)Equilibrium distances of H2 and He to M-MOF-74;(d)Binding energy of H2 and He to M-MOF-74

        3 結(jié) 論

        通過分子模擬的方法,首先,在不同的環(huán)境壓力下計算了H2及He 在M-MOF-74 上的吸附量。然后,研究了不同濃度的H2∕He 混合物在M-MOF-74上H2和He 的競爭吸附行為。最后,解析了H2和He在M-MOF-74上的吸附機理,得到以下主要結(jié)論。

        (1)采用M-MOF-74 作為吸附劑分離H2∕He 混合物時,Ni-MOF-74 具有最高的選擇性,而Mg-MOF-74 對H2具有最大的吸附量,其在壓力為20 bar時值為3.33 mmol·cm-3,為同等條件下He吸附量的5.71 倍。這表明MOF-74 作為吸附材料分離H2∕He混合物時在選擇性方面有一定的優(yōu)勢,具有應(yīng)用于工業(yè)生產(chǎn)的潛力。

        (2)M-MOF-74上OMS位點的H2吸附熱分布在12~17 kJ·mol-1之間,而有機配體對H2的吸附熱分布在0~12 kJ·mol-1。OMS 位點的存在幾乎不影響He吸附熱的分布,He 吸附熱分布在0~2.7 kJ·mol-1之間。由于M-MOF-74 上OMS 吸附位點的存在,強化了M-MOF-74與H2的相互作用,最終提高了H2在M-MOF-74上的吸附量和選擇性。

        (3)對于不同金屬離子的MOF-74 而言,H2與M-MOF-74 上不同金屬中心結(jié)合順序為Ni2+>Co2+>Mg2+>Zn2+>Cu2+,而He 與不同金屬中心結(jié)合的能力皆弱于H2且相近,這得益于H2為雙原子分子,存在較弱的色散偶極矩,而He 具備明顯的化學(xué)惰性,使得H2能夠通過與金屬中心產(chǎn)生弱的配位作用,來增強H2在M-MOF-74上的吸附。

        猜你喜歡
        環(huán)境壓力混合物選擇性
        Keys
        多組分纖維混合物定量分析通用計算模型研制
        正丁醇和松節(jié)油混合物對組織脫水不良的補救應(yīng)用
        故障狀態(tài)下純電動汽車環(huán)境壓力及海拔高度估算方法
        北京汽車(2021年1期)2021-03-04 13:05:46
        可替換牙刷
        選擇性聽力
        選擇性應(yīng)用固定物治療浮膝損傷的療效分析
        選擇性執(zhí)法的成因及對策
        混合物按照歐盟CLP進行分類標(biāo)簽
        貴陽市經(jīng)濟發(fā)展與環(huán)境壓力實證分析
        婷婷开心五月综合基地| 国产真实伦在线观看| 东北无码熟妇人妻AV在线| a√无码在线观看| 免费在线国产不卡视频| 亚洲国产一二三精品无码| 国产无遮挡又黄又爽在线视频| 日韩欧美国产亚洲中文| 92自拍视频爽啪在线观看| 77777亚洲午夜久久多喷| 一级黄色一区二区三区视频| 亚洲综合精品一区二区| 亚洲av色影在线| 亚洲国产精品久久久久秋霞影院| 日本女优中文字幕看片| av人妻在线一区二区三区| 日韩视频在线观看| 久久棈精品久久久久久噜噜| 国产乱子伦视频一区二区三区| 极品少妇高潮在线观看| 国产精品无码dvd在线观看| 亚洲一本大道无码av天堂| 国产欧美亚洲精品第二区首页| 国产精品女同一区二区软件| 日本丰满熟妇videossex一| 人妻在卧室被老板疯狂进入国产| 五月激情狠狠开心五月| 上海熟女av黑人在线播放| 精品亚洲成a人片在线观看| 亚洲天堂中文| 日韩伦理av一区二区三区| 国产无遮挡aaa片爽爽| 亚洲精品综合一区二区三| 午夜在线观看有码无码| 一区二区在线观看精品在线观看| 一进一出一爽又粗又大| 国产小毛片| 国产麻豆国精精品久久毛片| 中文字幕亚洲综合久久| a在线观看免费网站大全| 性色av成人精品久久|