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

腐蝕介質在緩蝕劑膜中擴散行為的分子動力學模擬

2010-10-14 03:42:52喬貴民任振甲胡松青燕友果
物理化學學報 2010年11期
關鍵詞:擴散系數體系

喬貴民 任振甲 張 軍 胡松青 燕友果 提 陽

(中國石油大學物理科學與技術學院,山東東營 257061)

腐蝕介質在緩蝕劑膜中擴散行為的分子動力學模擬

喬貴民 任振甲 張 軍*胡松青 燕友果 提 陽

(中國石油大學物理科學與技術學院,山東東營 257061)

采用分子動力學模擬方法,從緩蝕劑膜阻礙腐蝕介質粒子(H2O、H3O+和HCO-3)向金屬表面擴散的角度,研究了4種1-R1-2-十一烷基-咪唑啉緩蝕劑(R1:羧甲基(A),羥乙基(B),氨乙基(C),氫(D))抑制碳鋼CO2腐蝕的緩蝕機理,并對其緩蝕性能進行了理論評價.腐蝕介質粒子在不同緩蝕劑膜中的擴散系數、粒子與膜的相互作用能以及膜的自擴散性能的計算結果表明:4種緩蝕劑均可形成穩定的緩蝕劑膜,能有效阻礙腐蝕介質粒子向金屬表面的擴散,達到抑制或延緩腐蝕的目的;隨親水支鏈(R1)極性的增加,緩蝕劑膜對腐蝕介質粒子擴散行為的抑制能力逐漸增強;同種緩蝕劑膜對正負離子H3O+和HCO-3比對中性的H2O分子具有更強的擴散抑制能力.綜合計算及分析結果,4種緩蝕劑緩蝕性能的理論評價結果為A>B>C>D,與文獻實驗結果吻合.

分子動力學模擬; 緩蝕劑膜; 腐蝕介質; 擴散系數

Abstract:The corrosion inhibition mechanism of four 1-R1-2-undecyl-imidazoline inhibitors(R1:CH2COOH(A),CH2CH2OH(B),CH2CH2NH2(C),H(D))for carbon steel against carbon dioxide corrosion was investigated by molecular dynamics simulation,from the aspect of corrosive medium particle (H2O,H3O+,and HCO-3)diffusion to the metal surface hindered by the corrosion inhibitor membrane.The corrosion inhibition performance of these compounds was also evaluated by the theoretical method.The diffusion coefficients in various corrosion inhibitor membranes,the interaction energies between particles and membranes,and the self-diffusion performance of the membranes indicated that the four imidazoline inhibitors could form stable membranes,which could effectively limited the diffusion of corrosive medium particles to the metal surface,in order to inhibit or delay corrosion.With an increase in the polarity of the hydrophilic chain(R1),the ability of the membrane to hinder particle diffusion enhanced.The membrane was better at limiting the diffusion of charged particles(H3O+and HCO-3)than that of a neutral particle(H2O).Based on the above analysis,we found that theoretically the corrosion inhibition performance of the four imidazoline inhibitors decreased as follows:A>B>C>D,which agrees with previous experimental results.

Key Words: Molecular dynamics simulation;Corrosion inhibitor membrane;Corrosive medium;Diffusion coefficient

腐蝕是工業生產過程中一個極為突出的問題,造成了巨大的經濟損失和安全隱患.在眾多的防腐蝕措施之中,添加緩蝕劑是廣泛采用的方法之一.常用的緩蝕劑主要為吸附型緩蝕劑,它可穩定吸附在金屬表面并在金屬表面形成致密的保護膜,阻礙腐蝕介質向金屬表面遷移擴散,以達到抑制或延緩腐蝕的目的[1].目前研究人員[2-6]主要通過腐蝕電流、電化學阻抗及緩蝕效率等宏觀參量的變化來評價緩蝕劑膜性能的優劣,由于實驗技術本身的限制,還不能揭示緩蝕劑膜阻礙腐蝕介質擴散的微觀機理,緩蝕機理還有待于進一步完善.

目前,分子動力學模擬(MD)方法已發展成為一種能從分子水平對復雜體系進行研究的有效手段,可跟蹤復雜體系隨時間的動態演化過程[7-9],描述微小時間尺度內粒子的動態擴散行為.國內外學者[10-16]曾采用MD方法對聚合物膜、反滲透膜等體系中小分子的擴散動力學過程進行研究,為深入分析小分子擴散行為的微觀機制提供了依據.但針對腐蝕介質在緩蝕劑膜中擴散行為的MD研究卻鮮有報道.鑒于此,本文以4種含有不同親水支鏈的1-R1-2-十一烷基-咪唑啉緩蝕劑為研究對象(見表1[17]),采用分子動力學模擬方法,研究CO2腐蝕環境中存在的3種腐蝕介質粒子在不同緩蝕劑膜中的擴散行為;通過計算腐蝕介質粒子的擴散系數、粒子與緩蝕劑膜的相互作用能和緩蝕劑膜的自擴散系數等來分析緩蝕劑膜抑制腐蝕介質粒子擴散的微觀機理,并在此基礎上對4種咪唑啉緩蝕劑抑制CO2腐蝕的緩蝕性能進行理論評價.

1 計算方法

本文所有模擬工作均由Accelrys公司開發的Materials Studio軟件包[18]完成,力場為COMPASS力場[19],它能夠準確給出孤立體系和凝聚態的結構與性質.首先利用visualizer模塊構建A、B、C和D的分子模型,通過discover模塊下的minimizer工具采用smart minimizer方法對其結構進行優化,然后利用amorphous cell模塊分別構建包含80種緩蝕劑分子的無定形組織結構,采用等溫等壓系綜(NPT)進行分子動力學模擬,平衡后統計體系密度的平均值,并把該值作為計算體系中緩蝕劑膜的密度.模擬體系由三層結構組成:第一層為Fe(001)面,其尺寸為3.153 nm×3.153 nm×1.576 nm;第二層為包含80個緩蝕劑分子和1個腐蝕介質粒子的無定形組織結構;第三層是厚度為2 nm的真空層;然后通過visualizer模塊下的“build layer”命令將三層結構組合在一起即可獲取本文模擬所需的計算模型,如圖1所示.

在室溫條件下,金屬表面原子的熱振動可以忽略,因此在計算過程中“凍結”表面體系中的所有Fe原子[20].腐蝕介質粒子在緩蝕劑膜中的擴散行為的模擬通過discover模塊的正則系綜(NVT)[21]來實現,模擬溫度為298 K,溫度和壓力分別采用Andersen[22]和Berendsen[23]方法控制,各分子起始速度由Maxwell-Boltzmann分布隨機產生,在周期性邊界條件和時間平均等效于系綜平均等假設基礎上,運用velocity verlet算法[24]求解牛頓運動方程.范德華和庫侖相互作用采用charge group方法[25]計算.截斷半徑選取1.5 nm(spline width:0.10 nm,buffer width:0.05 nm),截斷距離之外的分子間相互作用按平均密度近似方法進行校正.模擬時間為2000 ps,步長為1.0 fs,每1000步記錄一次體系的軌跡信息(幀),共輸出2000幀.

通過溫度和能量判據來判斷體系是否已達到平衡,圖2為H2O分子在D緩蝕劑膜中擴散時體系的能量和溫度隨時間演化曲線,可以看出500 ps后體系溫度和能量波動平緩,能量偏差僅為2.35%,這表明此后體系已達到平衡,分析其它模擬體系也可得出相同的結論.因此各體系均取500-1900 ps內的R1:hydrophilic group,R2:long chain alkyl—(CH2)10CH3;bThe inhibition efficiency(IE)is taken from Ref.[17],which was investigated through linear polarisation resistance method.數據(最后100 ps數據存在統計誤差故舍去)計算腐蝕介質粒子的擴散系數、粒子與膜的相互作用能及緩蝕劑分子的自擴散系數等參量.

表1 四種咪唑啉類緩蝕劑的分子結構和緩蝕效率實驗值Table 1 Conformations and experimental inhibition efficiencies(IE)for four imidazoline inhibitors

2 結果與討論

界面型緩蝕劑緩蝕性能的優劣,關鍵在于緩蝕劑膜能否有效抑制腐蝕介質向金屬表面擴散.腐蝕介質在緩蝕劑膜中的擴散行為是一個極其復雜的過程,它不僅取決于粒子自身的結構、電荷分布和遷移能力等,同時也與腐蝕介質和膜之間的相互作用以及緩蝕劑膜的自擴散性能等密切相關.

2.1 腐蝕介質粒子在緩蝕劑膜中的擴散系數

擴散系數可定量描述腐蝕介質粒子在緩蝕劑膜中擴散能力的強弱.腐蝕介質粒子的擴散系數越小,說明粒子的擴散遷移能力越弱,緩蝕劑膜的緩蝕性能越好;反之,緩蝕性能則越差.腐蝕介質粒子的擴散系數(D)可通過愛因斯坦關系式[26]求得,公式如下:

其中,Ri(t0)和Ri(t)分別為第i個分子或粒子在初始時刻和t時刻的位置,n為分子或粒子的總數.

腐蝕介質粒子(H2O、H3O+和HC)在純水體系和緩蝕劑膜中的擴散系數列于表2.由表2可見,純水體系中H2O分子自擴散系數的計算值為2.3658×10-9m2·s-1,與實驗室測得擴散系數的數值2.30×10-9-2.50×10-9m2·s-1相吻合[27-29],說明本文所采用的計算方法和模型是合理和可靠的.與純水體系相比,3種腐蝕介質粒子在緩蝕劑膜中的擴散系數均大幅度下降,說明4種咪唑啉緩蝕劑所形成的緩蝕劑膜均能有效抑制腐蝕介質向金屬表面擴散,都可以起到緩蝕的效果.對于同一緩蝕劑膜,分析3種腐蝕介質粒子在其中的擴散系數可以發現,H3O+和HC的擴散系數均比H2O的擴散系數小,這說明緩蝕劑膜對帶電粒子擴散的抑制能力更強,可有效提高腐蝕原電池的電阻,減緩腐蝕反應速率.

表2 腐蝕介質粒子在純水體系及緩蝕劑膜中的擴散系數Table 2 Diffusion coefficients of corrosive particles in aqueous phase and corrosion inhibitor membranes

從表2中還可以看出,3種腐蝕介質在緩蝕劑膜A-D中的擴散系數逐漸變大,說明緩蝕劑膜抑制粒子擴散的能力逐漸降低,4種咪唑啉緩蝕劑的緩蝕效率應滿足:A>B>C>D,這與實驗室測得緩蝕效率的變化規律相吻合[17].其原因可能與親水支鏈R1中官能團—COOH,—OH,—NH2和—H的極性有關.羧基(—COOH)是一個強極性基團,極性比羥基(—OH)和氨基(—NH2)都強,羧基中的羥基和羰基都可以形成氫鍵,提高緩蝕劑膜自身的穩定性,從而形成對腐蝕介質粒子緊密地“包裹”,因此緩蝕劑膜A對腐蝕介質粒子擴散的抑制能力應最強;羥基(—OH)和氨基(—NH2)也可形成分子間氫鍵,但是羥基的極性比氨基強,其分子間氫鍵也比較強[30],因此緩蝕劑膜B的緩蝕性能應比緩蝕劑膜C好,而緩蝕劑膜D的性能最差.

2.2 緩蝕劑膜與腐蝕介質粒子之間的相互作用

腐蝕介質在緩蝕劑膜中運移,將不可避免地與緩蝕劑膜發生相互作用,兩者之間相互作用強弱將直接影響腐蝕介質在膜中的擴散能力.為此,文中計算了緩蝕劑膜與腐蝕介質粒子之間的相互作用能(Eparticlefilm),公式如下:

其中,Efilm是不包含腐蝕介質粒子的膜體系能量,Eparticle是孤立腐蝕介質粒子的能量,Efilm+particle是包含腐蝕介質粒子的膜體系能量.

圖3為緩蝕劑膜與腐蝕介質粒子相互作用能的變化曲線.可見,同種腐蝕介質粒子與4種緩蝕劑膜的相互作用能數值均滿足:A>B>C>D,這說明親水支鏈官能團的極性越強,緩蝕劑膜對腐蝕介質粒子的“包裹”越緊密,越能有效抑制腐蝕介質向金屬表面擴散,緩蝕性能也越強,因此4種咪唑啉緩蝕劑的緩蝕性能應為A>B>C>D.

此外,從圖3可以看出,同種緩蝕劑膜與H3O+和HC的相互作用能遠大于與H2O分子的作用能,這可能與兩者之間相互作用的方式有關.從分子動力學模擬角度來說,緩蝕劑膜與腐蝕介質粒子之間的相互作用包括庫侖作用和范德華作用,而粒子本身的性質將顯著影響這兩種作用方式的強弱.表3給出了緩蝕劑膜與腐蝕介質粒子的庫侖相互作用能(ECoulomb)和范德華相互作用能(EvdW)數值.可知,3種腐蝕介質粒子與緩蝕劑膜之間的庫侖相互作用均占主導地位,占兩者相互作用能總量的80%以上;對于H3O+和HC來說,粒子均帶一個單位的電荷,與緩蝕劑膜之間的庫侖相互作用明顯高于H2O.因此,與中性粒子(H2O)相比,緩蝕劑膜對帶電粒子(H3O+和HC)向金屬表面擴散的抑制能力更強,有效阻礙帶電粒子在金屬表面聚集,減緩腐蝕反應的速率,從而起到緩蝕效果.

2.3 緩蝕劑膜的自擴散性能

緩蝕劑膜的自擴散行為對膜內孔洞結構的穩定性具有重要影響,可在一定程度上反映緩蝕劑膜抑制腐蝕介質擴散能力的強弱[31].如果緩蝕劑分子的自擴散能力越強,緩蝕劑膜的瞬時孔洞變化越劇烈,所包裹的腐蝕介質粒子的運動也會隨之加劇,緩蝕性能則越差.緩蝕劑膜的自擴散性能可通過均方位移(MSD)和自擴散系數曲線來衡量.MSD曲線[26]定義如下:

表3 腐蝕介質粒子與緩蝕劑膜的庫侖相互作用能(ECoulomb)和范德華相互作用能(EvdW)Table 3 Coulomb interaction energies(ECoulomb)and van der Waals interaction energies(EvdW)between corrosive particles and corrosion inhibitor membranes

根據MSD曲線利用公式(1)可計算出緩蝕劑分子的自擴散系數,緩蝕劑膜的MSD曲線越陡,自擴散系數越大,自擴散能力就越強.

4種緩蝕劑膜的MSD曲線和自擴散系數的變化規律見圖4.從圖4(a)中可以看出,4種緩蝕劑膜的MSD曲線均隨時間增加呈上升趨勢,這說明由于緩蝕劑分子自身的熱運動,膜內的孔洞是不斷變化的,可能導致腐蝕介質粒子在膜內的遷移.此外,緩蝕劑膜A的MSD曲線增加幅度最小,說明緩蝕劑膜A內的瞬時孔洞變化最小,膜內孔洞的連通性最差,膜結構最穩定,最不利于腐蝕介質粒子在膜內的擴散遷移,而B次之,再次之是C和D.

從圖4(b)可以看出,緩蝕劑膜A、B、C和D的自擴散系數分別為 4.9×10-12、6.2×10-12、7.7×10-12和15.2×10-12m2·s-1.4種緩蝕劑膜的自擴散系數大小為A<B<C<D,數值逐漸變大,表明膜內緩蝕劑分子的運動性能逐漸增強,膜結構穩定性逐漸降低.這主要是由于親水支鏈官能團的極性造成的,官能團極性越強,緩蝕劑分子之間交織作用越劇烈,形成緩蝕劑膜的致密性越好,因此膜內緩蝕劑分子彼此受到的牽制作用越明顯,運動遷移能力越弱,緩蝕性能應越好.從緩蝕劑膜的穩定性可以推斷出4種緩蝕劑的緩蝕效率順序為:A>B>C>D.

3 結 論

采用分子動力學模擬的方法,研究了4種1-R1-2-十一烷基-咪唑啉緩蝕劑膜(R1:羧甲基(A),羥乙基(B),氨乙基(C),氫(D))阻礙腐蝕介質粒子(H2O、H3O+和HC)向金屬表面擴散的機理,并對其抑制碳鋼CO2腐蝕的緩蝕性能進行了理論評價.擴散系數的計算結果顯示,緩蝕劑膜均可有效阻礙腐蝕介質向金屬表面擴散,從而達到緩蝕效果,且腐蝕介質粒子在不同膜內的擴散系數順序為A<B<C<D.相互作用能數據表明緩蝕劑膜與腐蝕介質粒子之間存在強烈的庫侖作用,同種緩蝕劑膜對帶電粒子(H3O+和HC)擴散的抑制能力明顯強于對中性粒子(H2O).緩蝕劑膜的自擴散系數說明,4種緩蝕劑均可形成穩定的膜結構,且穩定性大小為A>B>C>D.綜合上述計算結果,4種緩蝕劑緩蝕性能的理論評價結果為A>B>C>D,這與實驗結果完全吻合.

1 Zhang,T.S.Corrosion inhibitor.Beijing:Chemical Industry Press,2008:131-143 [張天勝.緩蝕劑.北京:化學工業出版社,2008:131-143]

2 Jovancicevic,V.;Ramachandran,S.;Prince,P.Corrosion,1998,55:449

3 Zhang,G.A.;Chen,C.F.;Lu,M.X.;Chai,C.W.;Wu,Y.S.Mater.Chem.Phys.,2007,105:331

4 Liu,X.Y.;Chen,S.H.;Tian,F.;Ma,H.Y.;Shen,L.X.;Zhai,H.Y.Surf.Interface Anal.,2007,39:317

5 Zhang,Z.;Chen,S.H.;Li,Y.H.;Li,S.H.;Wang,L.Corrosion Sci.,2009,51:291

6 Zhang,J.;Du,M.;Yu,H.H.;Wang,N.Acta Phys.-Chim.Sin.,2009,25:525 [張 靜,杜 敏,于會華,王 寧.物理化學學報,2009,25:525]

7 Zhang,J.;Zhao,W.M.;Guo,W.Y.;Wang,Y.;Li,Z.P.Acta Phys.-Chim.Sin.,2008,24:1239 [張 軍,趙衛民,郭文躍,王 勇,李中譜.物理化學學報,2008,24:1239]

8 Domiínguez,H.Langmuir,2009,25:9006

9 Kornherr,A.;Nauer,G.E.;Sokol,A.A.;French,S.A.;Catlow,C.R.A.;Zifferer,G.Langmuir,2006,22:8036

10 Huang,Y.;Liu,Q.L.Journal of Xiamen University:Nature Science,2006,45:664 [黃 宇,劉慶林.廈門大學學報:自然科學版,2006,45:664]

11 Zhou,J.H.;Zhu,R.X.;Zhou,J.M.Polymer,2006,47:5206

12 Pan,F.S.;Peng,F.B.;Jiang,Z.Y.Chem.Eng.Sci.,2007,62:703

13 Yang,J.Z.;Liu,Q.L.;Wang,H.T.J.MembraneSci.,2007,291:1

14 Liu,Q.Z.;Yang,D.F.;Hu,Y.D.Chem.J.Chin.Univ.,2009,30:568 [劉清芝,楊登峰,胡仰棟.高等學校化學學報,2009,30:568]

15 Tao,C.G.;Feng,H.J.;Zhou,J.;Lü,L.H.;Lu,X.H.Acta Phys.-Chim.Sin.,2009,25:1373 [陶長貴,馮海軍,周 健,呂玲紅,陸小華.物理化學學報,2009,25:1373]

16 Pan,F.S.;Ma,J.;Cui,L.;Jiang,Z.Y.Chem.Eng.Sci.,2009,64:5192

17 Liu,X.;Zheng,Y.G.Corros.Eng.Sci.Technol.,2008,43:87

18 Materials Studio.Revision 4.2W.San Diego,USA:Accelrys Inc.,2005

19 Sun,H.J.Phys.Chem.B,1998,102:7338

20 Hansal,W.E.G.;Besenhard,J.O.;Kronberger,H.;Nauer,G.E.;Zifferer,G.J.Chem.Phys.,2003,119:9719

21 Heermann,D.W.Computer simulation methods in theoretical physics.Trans.Qin,K.C.Beijing:Peking University Press,1996:42-48 [Heermann,D.W.理論物理學中的計算機模擬方法.秦克成,譯.北京:北京大學出版社,1996:42-48]

22 Andersen,H.C.J.Chem.Phys.,1980,72:2384

23 Berendsen,H.J.C.;Postma,J.P.M.;van Gunsteren,W.F.;DiNola,A.;Haak,J.R.J.Chem.Phys.,1984,81:3684

24 Wu,X.H.;Xiang,J.Z.Modern materials computation and design.Beijing:Electronic Industry Press,2002:18-21 [吳興惠,項金鐘.現代材料計算與設計教程.北京:電子工業出版社,2002:18-21]

25 Leach,A.P.Molecular modeling:principles and application.England:Person Education Limited,2001:324-334

26 Lin,Y.C.;Chen,X.Chem.Phys.Lett.,2005,412:322

27 McCall,D.W.;Douglass,C.D.J.Chem.Phys.,1965,69:2001

28 Trappeniters,N.J.;Gerritsma,C.J.;Oosting,P.H.Phys.Lett.,1965,18:256

29 Aouizerat-Elarby,A.;Dez,H.;Prevel,B.;Jal,J.F.;Bert,J.;Dupuy-Philon,J.J.Mol.Liq.,2000,84:289

30 Zhou,L.Organic chemistry.Beijing:Science Press,2009:355-362[周 樂.有機化學.北京:科學出版社,2009:355-362]

31 Bernard,M.C.;Duval,S.;Joiret,S.;Keddam,M.;Ropital,F.;Takenouti,H.Prog.Org.Coat.,2002,45:399

Molecular Dynamics Simulation of Corrosive Medium Diffusion in Corrosion Inhibitor Membrane

QIAO Gui-Min REN Zhen-Jia ZHANG Jun*HU Song-Qing YAN You-Guo TI Yang
(College of Physics Science and Technology,China University of Petroleum,Dongying 257061,Shandong Province,P.R.China)

O641

Received:May 6,2010;Revised:July 26,2010;Published on Web:September 9,2010.

*Corresponding author.Email:dynamic_zh@163.com;Tel:+86-546-8391592.

The project was supported by the China National Petroleum Corporation Innovation Foundation(07E1021,2008D-5006-02)and Natural Science Foundation of Shandong Province,China(Y2006B35).

中國石油中青年創新基金(07E1021,2008D-5006-02)及山東省自然科學基金(Y2006B35)資助項目

猜你喜歡
擴散系數體系
構建體系,舉一反三
探索自由貿易賬戶體系創新應用
中國外匯(2019年17期)2019-11-16 09:31:14
一類具有變擴散系數的非局部反應-擴散方程解的爆破分析
基于Sauer-Freise 方法的Co- Mn 體系fcc 相互擴散系數的研究
上海金屬(2015年5期)2015-11-29 01:13:59
FCC Ni-Cu 及Ni-Mn 合金互擴散系數測定
上海金屬(2015年6期)2015-11-29 01:09:09
如何建立長期有效的培訓體系
現代企業(2015年1期)2015-02-28 18:43:18
非時齊擴散模型中擴散系數的局部估計
“曲線運動”知識體系和方法指導
Ni-Te 系統的擴散激活能和擴散系數研究
上海金屬(2013年4期)2013-12-20 07:57:07
香蕉凍干加工的水蒸氣擴散系數
食品科學(2013年13期)2013-03-11 18:24:13
主站蜘蛛池模板: 高潮毛片免费观看| 五月天综合婷婷| 日韩精品一区二区三区swag| 国产日本欧美在线观看| 国产啪在线91| 国产不卡网| 欧美天堂在线| 伊人久久精品亚洲午夜| 日韩欧美在线观看| 日韩精品免费一线在线观看| 亚洲精品无码AV电影在线播放| 久久伊人色| 亚洲人成电影在线播放| 综合久久久久久久综合网| 在线观看免费国产| 国产一级毛片在线| 亚洲成人福利网站| 亚洲一区二区三区香蕉| 日韩无码真实干出血视频| 亚洲国产日韩一区| 天天综合网色中文字幕| 男女男免费视频网站国产| 久久大香伊蕉在人线观看热2| 亚洲最新地址| 午夜限制老子影院888| 97国产在线视频| 精品自拍视频在线观看| 亚洲欧洲综合| 亚洲综合网在线观看| 99一级毛片| 欧美中出一区二区| 全裸无码专区| 国产 日韩 欧美 第二页| 日韩在线视频网| 国产精品女熟高潮视频| 久青草国产高清在线视频| 国产乱子精品一区二区在线观看| 99久久国产综合精品女同| 国产XXXX做受性欧美88| 四虎精品国产AV二区| 午夜电影在线观看国产1区| 亚洲视频影院| 香蕉久久国产超碰青草| 国产呦精品一区二区三区下载 | 一级毛片免费观看久| 国产成人三级| 无码aaa视频| 亚洲综合色区在线播放2019 | 国产中文一区a级毛片视频| 色精品视频| 久久精品视频亚洲| 亚洲av成人无码网站在线观看| 精品91在线| 六月婷婷精品视频在线观看 | 久久熟女AV| 午夜无码一区二区三区在线app| 天天综合网色中文字幕| 欧美在线免费| 日本免费新一区视频| 国产精品亚洲天堂| 一级毛片免费高清视频| 国产精品va免费视频| 国产精品成人AⅤ在线一二三四 | 91精品免费高清在线| 人妻精品久久久无码区色视| 伊人久久婷婷| 国产成人艳妇AA视频在线| 国产午夜精品一区二区三区软件| 国产乱子伦一区二区=| 国产一级片网址| 2020国产在线视精品在| 91国内视频在线观看| 亚洲中文字幕在线精品一区| 波多野结衣一二三| 国产自在线拍| 老汉色老汉首页a亚洲| 无码专区第一页| 久久黄色一级片| 国产成人高清精品免费软件| 一级毛片免费的| 欧美精品xx| 国产精品亚洲欧美日韩久久|