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

功能型離子液體協同吸收NH3和CO2的密度泛函理論研究

2022-11-13 07:32:52朱先會王甫夏杰成袁金良
化工學報 2022年10期
關鍵詞:優化

朱先會,王甫,夏杰成,袁金良

(寧波大學海運學院,浙江 寧波 315832)

引 言

CO2捕集與封存(CCS)是我國應對氣候變化的重要戰略選擇,也是我國履約減排任務、應對未來挑戰的重要技術選擇之一[1]。以氨水為吸收劑的氨法碳捕集技術具備高CO2吸收容量和低再生能耗等優勢,且不存在設備腐蝕、氧化降解等問題,成為近年來國內外“低能耗碳捕集”的一個重要研究方向。然而,由于NH3的高揮發特性,NH3逃逸問題仍是氨法碳捕集工藝中的一大技術壁壘[2]。因此,開發合適的添加劑來抑制NH3的揮發顯得尤為重要。

離子液體(ILs)作為一種新型綠色溶劑,具有極低的揮發性、不會造成氣相污染、低損耗、解吸能耗低、常溫下可以溶解大多數氣體、穩定性好、性質可調等優勢,在氣體脫除領域獲得了廣泛的關注[3-5]。利用ILs 的低揮發性、低飽和蒸氣壓以及高穩定性進行NH3的吸收可以解決傳統水洗法脫除NH3存在的再生能耗高、廢水量大以及氨精餾困難等問題[6]。陳晏杰等[7]提出了以[C4mim][BF4]為吸收劑的NH3吸收與多級閃蒸回收工藝,其流程模擬結果顯示NH3的回收率高達93.3%。Zeng等[8]研究發現金屬ILs與硅膠結合的多孔材料可以大量吸附NH3。Xu 等[9]通過實驗發現ILs 添加劑可以降低氨法碳捕集中NH3的逃逸。Shi 等[10]通過分子模擬研究發現陽離子對NH3的吸收起主導作用,主要原因是[Emim]陽離子與NH3之間形成了較強的氫鍵。Shang 等[11]根據ILs的結構可調性設計合成了三種含有不同數量氫鍵的離子液體([Bmim][NTf2]、[Bim][NTf2]、[HOOC(CH2)3mim][NTf2]),實驗研究結果顯示,[Bim][NTf2]具有最高的NH3吸收容量(超過2.69 mol NH3·(mol ILs)-1)。曾少娟等[12]通過總結發現羥基功能化的ILs比普通的ILs可以吸收更多的NH3,主要在于羥基與NH3形成了較強的氫鍵。

進一步地,ILs 在吸收CO2方面也呈現出了較大的潛力。Blanchard 等的研究表明CO2在[Bmim][PF6]中具有很高的溶解度[13],隨后對多種常規ILs溶解度的研究證實了其對CO2具有較高的選擇性[14]。在大多數ILs 吸收CO2的過程中,陰離子發揮了重要作用[15],主要表現在CO2在ILs 中的物理溶解和ILs 與CO2弱相互作用的結合[16],以及ILs 在水中的溶解度和親和力。同時,一些含有特定化學反應的特殊基團被設計引入到傳統的ILs 中,以開發出具有低蒸氣壓、低再生能耗、抗氧化能力強和高吸收容量等特性的功能型ILs[17]。比如,一些研究學者借鑒有機胺吸收CO2的原理,將氨基引入到ILs 的陽離子上,從而極大地改善了ILs吸收CO2的能力[18-19]。

雖然ILs 在吸收CO2和NH3方面均表現出了良好的效果,但由于ILs 對CO2和NH3吸收作用方式的不同,目前研究的ILs 難以同時兼顧CO2和NH3的吸收效果[9]。同時,ILs 對CO2和NH3的吸收能力均受ILs中陰、陽離子的極性,酸堿性,密度,黏度,以及功能基團等因素的影響[20],需對ILs的結構和功能基團進行設計與優化,以達到對NH3和CO2的協同吸收效果。基于此,本文根據文獻中ILs 與NH3和CO2作用時最可能的結構以及ILs 陰、陽離子對CO2和NH3作用方式的不同,設計五種功能型ILs,利用密度泛函理論對其結構進行優化及頻率計算,探討ILs 單獨吸收NH3和CO2及協同吸收的作用機理以及吸收能力。

1 計算方法

密度泛函理論(density functional theory, DFT)是一種研究多電子體系電子結構的方法。DFT在物理和化學等領域都有廣泛的應用,特別是在研究分子和凝聚態的性質等方面,是計算化學領域最常用的方法之一。本研究采用DFT 中的B3LYP∕6-31’++G(d,p)基組研究ILs 構型及其與氣體分子之間的相互作用。電子密度拓撲分析采用Multiwfn 軟件包[21]來實現,其他所有的計算均依據Gaussian 09[22]程序完成。

1.1 相互作用能的計算

相互作用能是離子對與其相應的陰、陽離子之間的能量差,由式(1)進行計算得到。在相互作用能的計算中,通常會有基組重疊誤差(base set superposition error, BSSE)產生,導致計算出的相互作用能偏大。因而,在研究分子間的弱相互作用時,需要消除基組重疊誤差的影響。利用Boys 等[23]提出的理論,由式(2)進行基組疊加誤差計算。進一步地,通過式(1)和式(2)可推導出式(3)所示的ILs離子間相互作用能的計算方法。

式中,EILs是陰陽離子對的電子能量;Ecation是陽離子的電子能量;Eanion是陰離子的電子能量;EBSSE是基組疊加誤差。

在DFT 理論6-31 G(d,p)基組條件下,首先對五種ILs 的構型進行預優化,通過ILs 相互作用能的計算得到離子對較為穩定的構型。然后選用更復雜的6-31’++G(d, p)基組進行幾何構型優化,以獲得能量最低的穩定結構。

1.2 AIM理論分析

基于優化后的ILs 構型,通過Bader 分子中的原子理論(AIM)[24],在Multiwfn 軟件中進行電子密度拓撲理論分析。AIM 方法已被證明可用于ILs 體系中的氫鍵研究[25-26]。通過計算可得到電子密度、拉普拉斯電子密度值,并用式(4)得到氫鍵結合能(EHB,hydrogen bonding energy),以此來分析形成的氫鍵強弱。同時,使用Multiwfn 可以分析靜電勢(ESP)、原子偶極矩以及原子電荷,以預測親電和親核最可能的攻擊位置,為ILs 與氣體分子的相互作用提供理論基礎。

2 結果與分析

2.1 離子液體的構型優化及相互作用能

在B3LYP 雜化泛函6-31’++G(d,p)基組水平上對 設 計 的 五 種ILs([HEMim][Glu]、[HEMim][Asp]、[HEBim][Asp]、[HEBim][Ala]、[HEBim][His])進 行 結構優化。這五種ILs 由兩種羥基化的咪唑型陽離子和四種氨基酸陰離子組合而形成不同的空間結構。圖1 所顯示的是能量最低、最穩定的結構。經過結構優化和頻率分析證實各ILs 均無虛頻,可以穩定存在。其中咪唑環上的二面角N15-C1-N16-C2 和N16-C2-C3-N15 分別為0.1195°和0.1123°,表明咪唑環上的各原子基本處在同一個平面上;咪唑環上的原子與其周圍原子的二面角大都在179.6°左右,表明咪唑環與其周圍各個原子也基本處在同一個平面內,這為ILs自身的穩定提供了重要支持。

根據價鍵理論,當兩原子間距小于相應兩原子之間的范德華半徑之和而大于共價鍵鍵長,并且其相應角度大于90°時則認為形成了氫鍵。從圖1 可以看出,這五種ILs 內部形成了較多的氫鍵,這為離子液體的穩定性提供了基本保證。根據文獻中的結果:氫鍵應該有相對較大的電子密度(ρBCP)和拉普拉斯電子密度(?2ρBCP)。其中ρBCP用來表示原子之間成鍵的強弱,ρBCP越大形成的氫鍵就越強。?2ρBCP是ρBCP的二階導數,可以用來判斷原子之間的成鍵類型。當?2ρBCP<0,則表明鍵點(BCP)處的電荷是積聚的,相鄰的兩個原子之間以共價鍵的形式存在;相反,如果?2ρBCP>0,則表明BCP 處的電荷是發散的,相鄰兩個原子之間以氫鍵、離子鍵等閉殼層作用存在。一般來說,ρBCP的取值范圍為0.002~0.04 a.u.,?2ρBCP的取值范圍為0.02~0.15 a.u.。然而這只是傳統氫鍵的取值范圍,有一些氫鍵的ρBCP和?2ρBCP的取值范圍會比傳統氫鍵的取值范圍大。表1列出了五種ILs 的ρBCP和?2ρBCP以及氫鍵作用能(EHB),由此可判斷五種ILs 內部氫鍵的形成及強弱。例如:[HEMim][Asp]中形成氫鍵的?2ρBCP都為正值,H 原子和O 原子形成的化學鍵的鍵長為1.74~2.28 ?(1 ?=0.1 nm),ρBCP的范圍為0.013~0.045 a.u.,符合形成氫鍵的判據且形成了三個氫鍵,其中C1—H6…O28 處的氫鍵最強,其氫鍵結合能為38.45 kJ·mol-1。

表1 在B3LYP/6-31’++G(d,p)基組水平下計算得到的五種離子液體的電子密度性質Table 1 Properties of electron density of BCP of five ionic liquids calculated at B3LYP/6-31’++G(d,p)level

圖1 在B3LYP∕6-31’++G(d,p)基組水平下優化得到的五種離子液體較為穩定的構型Fig.1 Optimized structures of five ionic liquids at B3LYP∕6-31’++G(d,p)level

圖2 列出了所設計的五種ILs 的零點作用能以及經過BSSE 校正的相互作用能。從圖中可以看出,五種ILs 的零點作用能和經過BSSE 校正的相互作用能在-421~-390 與-416~-388 kJ·mol-1之間。與張營[27]研究的[C3mim] [Glu] 的相互作用能(-350.91~-306.42 kJ·mol-1)相比,本研究的ILs 相互作用能的絕對值更大。這是因為ILs 功能化之后,陰陽離子之間形成了更多的氫鍵,使得其穩定性更好。通過以上分析可以看出[HEBim][His]最穩定,其經過BSSE 校正之后的相互作用能達到-415.73 kJ·mol-1。

圖2 在B3LYP∕6-31’++G(d,p)基組水平下計算得到的五種離子液體的相互作用能Fig.2 Interaction energies for five ionic liquids at B3LYP∕6-31’++G(d,p)level

2.2 離子液體的靜電勢分析與原子電荷分析

靜電勢分析可以更好地理解ILs 與氣體分子之間的相互作用機理,找出ILs 與氣體相互作用的最佳位點。圖3 是五種ILs 的靜電勢能面示意圖,ILs表面的藍色和紅色分別對應的是正靜電勢和負靜電勢(ESP),其中顏色越深代表靜電勢越強。從圖中可以看出,正電勢大多集中在咪唑環附近,其中咪唑環上的H 原子正電勢最強,陽離子羥基中的H 原子正電勢次之。而負靜電勢大多集中在陽離子的O原子以及陰離子的O原子和N原子附近。上述結果說明咪唑環上的H 原子和羥基中的H 原子較易與NH3相 互 作 用,O 原 子 和N 原 子 較 易 與CO2相 互作用。

圖3 在B3LYP∕6-31’++G(d,p)基組水平下得到的五種離子液體的靜電勢能面圖Fig.3 Electrostatic potential surface of five ionic liquids at B3LYP∕6-31’++G(d,p)level

表2 給出了ILs 與氣體分子可能的相互作用點的ESP、Hirshfeld和ADCH電荷。Hirshfeld電荷通常比ADCH 電荷小,ESP 的再現性差,因為它忽略了原子偶極矩。與Hirshfeld 電荷相比,ADCH 電荷能很好再現ESP,并且可以準確地預測潛在的作用位點。例如:[HEMim][Glu]的電荷分析結果顯示H14 原子具有最大的正電荷,陽離子羥基中H14 原子是最有可能與NH3相互作用的原子,這與Li 等[28]和Yuan等[29]的結論相一致。[HEBim][Ala]的電荷分析得出H29 原子具有最高的正電荷,其可能是與NH3相互作用的最佳位點,而N39原子具有最高的負電荷,可能是與CO2相互作用的最佳位點。然而,H6 原子和H29 原子應該具有較高的正電荷,但靜電勢分析的結果卻與電荷分析結果不一致,最主要的原因是H6與H29原子和O原子之間形成氫鍵導致靜電勢分析結果出現誤差。

表2 在B3LYP/6-31’++G(d,p)基組下計算的離子液體的ESP、ADCH和Hirshfeld 電荷Table 2 ESP,ADCH and Hirshfeld charges for ionic liquids calculated at B3LYP/6-31’++G(d,p)level

2.3 離子液體與NH3的相互作用

在B3LYP雜化泛函6-31’++G(d,p)基組水平上對ILs與NH3的相互作用進行結構優化和頻率計算,以確定ILs 吸收NH3的作用方式。圖4 是優化后能量最低、最穩定的ILs與NH3相互作用的結構。可以看出,NH3主要與ILs陽離子中的羥基相互作用。通過鍵長與鍵角分析可以得出陽離子羥基中的H 原子與NH3中的N 原子形成了O—H…N 型氫鍵,這與前面的靜電和電荷分析中潛在的作用位點相一致。表3 給出了ILs 與NH3相互作用的電子密度拓撲分析結果。結果顯示:ILs 的陽離子羥基中的H 與NH3的N 形成的氫鍵的鍵長在1.780~1.796 ? 之間,ρBCP在0.043~0.045 a.u.之間,形成的氫鍵結合能在37.0~38.6 kJ·mol-1之間,其鍵長、ρBCP與?2ρBCP均符合形成氫鍵的條件,且也符合ρBCP越大氫鍵越強的規律。從以上分析可以得出五種ILs 吸收NH3的能力為[HEBim] [His] > [HEBim] [Ala] > [HEMim] [Asp] >[HEMim][Glu] >[HEBim][Asp]。其中,[HEBim][His]吸收NH3的能力最強,其形成氫鍵的鍵長為1.7883 ?,電子密度為0.0446 a.u.,形成的氫鍵結合能為38.52 kJ·mol-1。

表3 在B3LYP/6-31’++G(d,p)基組水平下計算得到的五種離子液體與氣體作用的鍵長以及電子密度性質Table 3 The bond lengths and electron density of five ionic liquids interacting with gases calculated at B3LYP/6-31’++G(d,p)level

圖4 在B3LYP∕6-31’++G(d,p)基組水平下優化得到的五種離子液體和氨相互作用的穩定構型Fig.4 Optimized structures of five ionic liquids and ammonia interaction at B3LYP∕6-31’++G(d,p)level

2.4 離子液體與CO2的相互作用

在B3LYP 雜化泛函6-31’++G(d,p)基組水平上進行ILs 與CO2相互作用的結構優化和頻率計算。圖5 是經過優化后的最穩定結構。可以看出CO2主要與ILs陰離子氨基中的N原子起作用,這也驗證了上述電荷分析和ILs 靜電勢分析的結果。經過鍵長與鍵角分析可以看出,CO2與氨基中的N 原子形成C—N…C 型氫鍵,鍵長大都在2.8~3.2 ? 之間。從表3 可以看出ILs 中的N 原子和CO2的C 原子的ρBCP在0.007~0.015 a.u.之間,形成的氫鍵結合能在3.0~11.0 kJ·mol-1之間,屬于較弱的氫鍵。根據形成氫鍵的強弱可以得出所設計的五種離子液體吸收CO2的強弱順序為[HEBim][Ala]>[HEBim][His]>[HEMim][Glu]>[HEMim][Asp]>[HEBim][Asp]。其中,[HEBim][Ala]吸收CO2的能力最強,其形成的氫鍵的鍵長為2.8215 ?,ρBCP為0.0142 a.u.,形成的氫鍵結合能為10.15 kJ·mol-1。

圖5 在B3LYP∕6-31’++G(d,p)基組水平下優化得到的五種離子液體和二氧化碳相互作用的穩定構型Fig.5 Optimized structures of five ionic liquids and carbon dioxide interaction at B3LYP∕6-31’++G(d,p)level

2.5 離子液體與NH3和CO2的相互作用

從ILs 分別與NH3和CO2的作用可以看出,雖然NH3和CO2與ILs 均是以氫鍵的方式進行結合的,但作用的分別是ILs 的陰離子和陽離子。考慮到氨法碳捕集中ILs 添加劑協同吸收NH3和CO2的效果,對其協同吸收時的相互影響需進一步分析。針對單獨作用時優選的[HEBim][His]與[HEBim][Ala],在B3LYP 雜化泛函6-31’++G(d,p)基組水平上對其同時與NH3和CO2的相互作用進行結構優化和頻率計算,其優化后的ILs與氣體相互作用的結構如圖6所示,結果表明:ILs同時與NH3和CO2相互作用時結構較穩定,無虛頻存在。NH3主要與陽離子中的羥基形成O—H…N 型氫鍵,CO2則主要與陰離子中的氨基形成C—N…C 型氫鍵,與NH3和CO2單獨與ILs 作用的結論相一致。兩種ILs 與氣體分子協同作用的電子密度拓撲分析如表4 所示,可以看出ILs 與NH3和CO2同時作用時的吸收能力均有不同程度的下降。例如:[HEBim][His]與NH3作用形成O17—H18…N48 型氫鍵,其鍵長為1.7886 ?、鍵角為169°,ρBCP為0.0445 a.u.,?2ρBCP為0.1034 a.u.,氫鍵結合能為38.43 kJ·mol-1。[HEBim][His]與CO2作用形成C41—N44…C53 型氫鍵,其鍵長為2.9830 ?,鍵角為110°,ρBCP為0.0112 a.u.,?2ρBCP為0.0345 a.u.,氫鍵結合能為7.35 kJ·mol-1。上述結果符合氫鍵的判據,同時也符合鍵長越短電子密度越大,氫鍵結合能越大的規律。通過分析可以看出對于吸收NH3的能力,[HEBIM][His]為38.43 kJ·mol-1,稍強于[HEBim][Ala]。而對于吸收CO2的能力,[HEBim][Ala]則稍強,為8.93 kJ·mol-1,兩者在抑制氨逃逸和吸收CO2方面的差異較小。為了進一步研究ILs 協同吸收NH3和CO2時吸收能力下降的原因,在B3LYP 雜化泛函6-31’++G(d,p)基組水平上對其進行結構優化和頻率計算,得到如圖7 所示的光譜圖。通過光譜分析可以看出:[HEBim][His]同時與NH3和CO2作用時羥基的吸收峰由3179左移到3115 cm-1,表明力常數變小,化學鍵變弱,使得吸收NH3的能力下降。

圖7 優化后的離子液體與氣體相互作用的光譜圖Fig.7 Spectra of interaction between ILs and gases after optimization

表4 在B3LYP/6-31’++G(d,p)基組水平下計算得到的兩種離子液體與氣體共同作用的鍵長、結合能以及電子密度等性質Table 4 The bond lengths,binding energies and electron densities of the two ILs interacting with gas calculated at the B3LYP/6-31++G(d,p)level

圖6 在B3LYP∕6-31’++G(d,p)基組水平下優化得到的兩種離子液體與氣體同時作用的穩定構型Fig.6 Two stale configurations of two ILs interacting with gases obtained at B3LYP∕6-31’++G(d,p)level

為了進一步研究ILs 同時與NH3和CO2的作用實質,采用了式(5)所示的約化密度梯度(RDG)方法來研究其弱相互作用[30],這種方法可以描述包括氫鍵在內的不同類型分子之間的弱相互作用。

因為電子密度只能描述相互作用的強度,而相互作用類型需要由第二個Hessian 矩陣的特征值符號來確定。例如sign(λ2)=1 和sign(λ2)=-1 分別表示無鍵和有鍵的相互作用。通過圖8 可以看出,在低密度和低梯度的區域內存在許多尖峰,表明ILs 與氣體之間存在弱相互作用。這些尖峰可以分為三類:強吸引相互作用如氫鍵[大且負的sign(λ2)ρ(r)]、范德華相互作用(接近0)以及強排斥相互作用如空間位阻效應[大且正的sign(λ2)ρ(r)]。圖8(a)的尖峰表示分子內和分子間的氫鍵。ILs 與氣體相互吸引作用的sign(λ2)ρ(r)值總結在表5 中。從表中可以看出[HEBim][His]與NH3和CO2之間都形成了氫鍵,其中[HEBim][His]與NH3的sign(λ2)ρ(r)值達到-0.0452 a.u.,表明O17—H18…N48 之間形成了較強的氫鍵,與CO2的sign(λ2)ρ(r)值 為-0.0092 a.u.,表 明C41—N44…C53之間形成了較弱的氫鍵,符合上述電子密度拓撲分析的結論。

圖8 通過第二個Hessian矩陣的特征值計算得出的離子液體與氣體相互作用的電子密度RDG圖Fig.8 Plot of RDG versus electron density multiplied by the sign of the second Hessian eigenvalue(λ2)for ILs-gas structures

表5 離子液體與氣體之間的氫鍵的sign(λ2)ρ(r)值Table 5 sign(λ2)ρ(r)values for H-bonds between ILs and gases

3 結 論

本文針對氨法碳捕集中NH3的高揮發性和CO2吸收慢的問題,設計了五種含有羥基及氨基的功能型ILs,即[HEMim] [Glu]、[HEMim] [Asp]、[HEBim][Asp]、[HEBim][Ala]、[HEBim][His],以探究其協同吸收NH3和CO2的作用方式。通過量子化學密度泛函理論研究了ILs 之間的構型及其與氣體分子之間的作用機理。相互作用能計算結果表明:[HEBim][His]的結構最為穩定,其內部形成了較多的氫鍵。電子密度拓撲分析表明:ILs的陽離子羥基中的H原子與NH3中的N 原子形成了O—H…N 型氫鍵,而CO2主要與陰離子中的氨基形成了較弱的C—C…N 型氫鍵。五種ILs 吸收NH3的能力為[HEBim][His]>[HEBim] [Ala] >[HEMim] [Asp] >[HEMim] [Glu] >[HEBim] [Asp],吸 收CO2的 能 力[HEBim] [Ala] >[HEBim] [His] > [HEMim] [Glu] >[HEMim] [Asp] >[HEBim][Asp]。在同時抑制NH3逃逸和促進CO2吸收性能方面,[HEBim][His]和[HEBim][Ala]均表現出了較好的效果,是比較適合的吸收劑。

猜你喜歡
優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
PEMFC流道的多目標優化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
圍繞“地、業、人”優化產業扶貧
今日農業(2020年16期)2020-12-14 15:04:59
事業單位中固定資產會計處理的優化
消費導刊(2018年8期)2018-05-25 13:20:08
4K HDR性能大幅度優化 JVC DLA-X8 18 BC
幾種常見的負載均衡算法的優化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 亚洲天堂免费观看| 中文字幕第4页| 五月天在线网站| 一区二区三区四区精品视频 | 国产日韩欧美成人| 天天干天天色综合网| 日韩一区二区三免费高清| 国产99欧美精品久久精品久久| 人妻一区二区三区无码精品一区 | AV不卡国产在线观看| 九色在线视频导航91| 国产精品开放后亚洲| 超薄丝袜足j国产在线视频| 97综合久久| 国产精品区视频中文字幕| 漂亮人妻被中出中文字幕久久 | 永久免费无码成人网站| 欧美专区在线观看| 国产粉嫩粉嫩的18在线播放91 | 在线毛片免费| 香蕉精品在线| 国产精品久久精品| 97se亚洲综合不卡 | 一本一本大道香蕉久在线播放| 在线观看热码亚洲av每日更新| 亚洲欧美一级一级a| 中文字幕波多野不卡一区| 国产美女在线观看| 国产精品黑色丝袜的老师| 欧美精品1区2区| 亚洲综合激情另类专区| 强奷白丝美女在线观看| 99视频全部免费| 九九九精品视频| 色香蕉影院| 乱码国产乱码精品精在线播放| 熟女成人国产精品视频| 91尤物国产尤物福利在线| 国产女人18水真多毛片18精品 | 欧美综合中文字幕久久| 国产精品视频白浆免费视频| 91人妻日韩人妻无码专区精品| 97国产成人无码精品久久久| 91成人在线观看视频| 国产制服丝袜91在线| 青草视频在线观看国产| 亚洲第一黄色网址| 玖玖精品视频在线观看| 国内自拍久第一页| 22sihu国产精品视频影视资讯| h网址在线观看| 无码丝袜人妻| 亚洲永久视频| 欧美精品H在线播放| 国产高清毛片| a级毛片一区二区免费视频| 福利姬国产精品一区在线| 国产性爱网站| 国产69精品久久久久孕妇大杂乱| 亚洲女同欧美在线| 亚洲性日韩精品一区二区| 国产一级二级在线观看| 国产成人91精品| 欧美怡红院视频一区二区三区| 中文字幕第1页在线播| 一级毛片在线播放| 亚洲欧美成人在线视频| 免费无码又爽又黄又刺激网站| 性网站在线观看| 欧美日韩久久综合| 欧美视频二区| 中文字幕1区2区| 免费av一区二区三区在线| 亚洲国产精品VA在线看黑人| 无码aaa视频| 国产精品理论片| 黄色一级视频欧美| 免费a在线观看播放| 色欲色欲久久综合网| 亚洲va精品中文字幕| 天天综合网色中文字幕| 亚洲最新地址|