• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    離子對CO2-鹽水系統(tǒng)界面特性的影響

    2017-01-19 01:08:08季佳圓趙伶玲李偲宇
    化工學報 2017年1期
    關(guān)鍵詞:原子質(zhì)量鹽濃度陽離子

    季佳圓,趙伶玲,李偲宇

    (東南大學能源與環(huán)境學院,能源熱轉(zhuǎn)換及其過程測控教育部重點實驗室,江蘇 南京 210096)

    離子對CO2-鹽水系統(tǒng)界面特性的影響

    季佳圓,趙伶玲,李偲宇

    (東南大學能源與環(huán)境學院,能源熱轉(zhuǎn)換及其過程測控教育部重點實驗室,江蘇 南京 210096)

    應用分子動力學模擬的方法對不同濃度及組分比例的 CO2-NaCl、CO2-CaCl2和 CO2-(NaCl+CaCl2)等 CO2-鹽水系統(tǒng)的界面張力進行了研究,從分子層面上分析了不同種類的鹽離子對界面張力的影響。結(jié)果表明,鹽離子的存在會增加CO2-純水的界面張力,界面張力增量主要由陽離子電荷數(shù)導致,其余則由陽離子蘭納-瓊斯勢參數(shù)、原子質(zhì)量等原子特性參數(shù)引起;3種系統(tǒng)的界面張力增量均與離子強度呈線性關(guān)系,且其曲線的斜率和鹽溶液種類無關(guān),為IFT預測提供理論指導。

    CO2-鹽水系統(tǒng);界面張力;分子模擬;原子特性參數(shù);離子強度

    引 言

    深部鹽水層分布廣泛、存儲容量大,因此被作為CO2地質(zhì)封存的儲存選址而備受推崇[1-2]。在CO2注入深部鹽水層時,流體之間的界面張力(interfacial tension,IFT)直接決定了封存安全性、注入能耗及最大存儲量,因此研究CO2與鹽水之間的界面相互作用及界面性質(zhì)具有重大意義。

    深部鹽水層由于種類及所在地的不同而存在離子成分差異,但大部分鹽水均以 Cl?為主要陰離子,Na+、Ca2+為主要的陽離子,附加少量的 K+、Mg2+等[3]。實驗發(fā)現(xiàn)鹽離子的存在會增加CO2-鹽水系統(tǒng)的界面張力,且IFT增量與鹽濃度呈線性增加關(guān)系[4-5]。此外,不同的陽離子對IFT增量具有不同的影響,如Ca2+的電荷數(shù)是Na+的2倍,而CaCl2對IFT增量的增加斜率(3.53 mN·L·m?1·mol?1)約為NaCl(1.43 mN·L·m?1·mol?1)的2.5倍[5],有實驗研究者認為,是由Ca2+與Na+的陽離子價態(tài)差異引起的[5]。然而,不同鹽成分導致不同CO2-鹽水系統(tǒng)IFT增量的內(nèi)在機理還尚未見報道。此外,由于離子電荷數(shù)、蘭納-瓊斯勢(Lennard-Jones (L-J) potential)參數(shù)、原子質(zhì)量等特性參數(shù)隨著離子同時變化,導致實驗中不能采用控制變量法進行單一原子特性參數(shù)變化的研究,所以難以實現(xiàn)離子各特性參數(shù)對 CO2-鹽水系統(tǒng)界面張力變化影響規(guī)律的單獨討論。

    除實驗研究外,分子動力學模擬作為研究CO2-鹽水體系界面張力的有效方法已被國內(nèi)外學者廣泛使用。da Rocha等[6]在 318 K和 20 MPa下采用SPC/E 水分子和 EPM2的 CO2分子模型研究了CO2-純水界面的 IFT,發(fā)現(xiàn)了分子熱波動引起的毛細波現(xiàn)象。Li等[7]采用EPM2的CO2模型以及SPC/E和TIP4P的水分子模型,在303~393 K and 2~50 MPa的范圍內(nèi)模擬了CO2-NaCl和CO2-CaCl2系統(tǒng),發(fā)現(xiàn)了離子電荷較離子半徑對IFT值的影響更大。Zhao等[8]采用柔性F3C水分子模型以及柔性EPM2 CO2分子模型得到了與實驗數(shù)據(jù)相一致的CO2-純水以及CO2-CaCl2系統(tǒng)的IFT值,同時分析出原子力對IFT具有一定的影響作用并運用周向分布函數(shù)推測出界面處分子的排布規(guī)律。然而,目前這些MD研究尚未對鹽溶質(zhì)的原子特性參數(shù)與IFT值的關(guān)系進行詳細研究。

    本文應用分子動力學模擬的方法,提出了虛擬系統(tǒng)與真實系統(tǒng)的比照模擬體系,研究了343 K和20 MPa 條 件 下 CO2-NaCl、 CO2-CaCl2和CO2-(NaCl+CaCl2)鹽水系統(tǒng)的界面張力等微觀性質(zhì),分析了不同陽離子種類及成分、陽離子電荷數(shù)、蘭納-瓊斯勢L-J勢能參數(shù)、原子質(zhì)量等原子特性參數(shù)對IFT的影響,探討了以上各參數(shù)對IFT的單獨作用機理。此外,本文還建立了IFT值與離子強度之間的直接聯(lián)系,為CO2在深部鹽水層中的存儲提供了理論依據(jù)。

    1 研究對象及方法

    本文的研究對象為343 K和20 MPa條件下,CO2-NaCl、CO2-CaCl2溶液和 CO2-(NaCl+CaCl2)混合鹽溶液界面系統(tǒng)。計算中考慮了分子間非鍵結(jié)作用及分子內(nèi)鍵結(jié)作用。其中分子間的非鍵結(jié)作用包括范德華力和庫侖靜電力,用 L-J勢能函數(shù)[9]描述范德華力,庫侖定律描述分子間庫侖靜電力,采用PME技術(shù)[10]模擬分子間長程庫侖作用力,范德華作用截距設定為0.9 nm[11]。而分子內(nèi)的鍵結(jié)作用則包括鍵拉伸和鍵角彎曲,分別采用諧波勢能函數(shù)[8]進行模擬。兩原子i和j之間的L-J勢能函數(shù)詳見式(1)。

    式中,εij=(εiiεjj)1/2,σij=(σiiσjj)1/2,rij為原子i和j之間的距離;ε為L-J勢能曲線中的勢能阱深度;勢能阱越深,兩點間引力越大,當勢能為零時,兩點間的距離等于σ。

    計算系統(tǒng)中CO2、水分子及離子的力場參數(shù)列于表1。其中,水分子采用了柔性F3C模型[12],CO2采用柔性 EPM2模型[13],Na+、K+和 Cl?采用Chandrasekhar等[14]開發(fā)的力場,Ca2+采用Aqvist[15]開發(fā)的力場。

    表1 CO2、H2O及鹽離子參數(shù)Table 1 Parameters of CO2, H2O and salt ions

    1.1 虛擬系統(tǒng)

    離子的電荷數(shù)(q)、L-J勢參數(shù)(ε、σ)、原子質(zhì)量(m)等均可表征鹽離子原子特性參數(shù)。本研究采用控制變量法依次研究上述參數(shù)對IFT的單獨作用。對給定離子類型的某種鹽而言,由于離子電荷和 L-J參數(shù)是相互依賴并共同作用于原子特性,從而在實驗角度上不可能實現(xiàn)改變單一變量的操作。然而MD模擬則可設置只改變某一參數(shù)而確保其他參數(shù)不變的鹽離子虛擬系統(tǒng),故而對于實現(xiàn)控制變量的研究思路具有其獨特的優(yōu)勢。為此,本文設置了CO2-NaCl和CO2-CaCl2溶液的虛擬系統(tǒng),列于表2。表2中系統(tǒng)1 (CO2-NaCl)、系統(tǒng)2(CO2- CaCl2)、系統(tǒng)11(CO2-KCl)為真實系統(tǒng),系統(tǒng)3*~10*(上標為“*”的系統(tǒng))為虛擬系統(tǒng),各虛擬系統(tǒng)的變化參數(shù)列于相應括號內(nèi)。由于氯離子在深部鹽水層中占據(jù)陰離子大部分份額[3],故本文僅選用Cl?作為陰離子且不探討陰離子種類及原子參數(shù)對 IFT的影響。

    表2 真實系統(tǒng)和虛擬系統(tǒng)的計算工況Table 2 Actual and fictitious simulation systems of CO2-electrolyte solution

    1.2 系統(tǒng)建模

    本文應用分子模擬軟件 Gromacs4.5[16]進行計算。計算中涉及的CO2-鹽水系統(tǒng)均采用三維周期性邊界條件,初始計算域(模擬盒子)的尺寸為4 nm× 4 nm×16 nm。其中,模型中部占z軸長度為8 nm區(qū)域內(nèi)為 4323個水分子和與鹽濃度一致的相應個數(shù)的鹽離子,例如鹽濃度為1.8 mol·L?1時的NaCl溶液,建模時設置139個Na+及Cl?。鹽水兩側(cè)分別為占z軸長度為4 nm的CO2,為使計算更快達到平衡,CO2分子個數(shù)根據(jù)計算工況溫度壓強下的CO2實驗密度值確定,343 K和20 MPa條件下CO2分子數(shù)為 1472。以 CO2-NaCl溶液(鹽濃度 1.8 mol·L?1)系統(tǒng)為例,建立的模型示于圖1。

    圖1 CO2-NaCl溶液系統(tǒng)分子模擬示意圖Fig.1 Snapshot of supercritical CO2-NaCl solution

    計算采用NPZT系綜條件。溫度及壓強的耦合分別采用基于 Berendsen耦合器的 Velocity Rescaling方法[17]和Semi-isotropic方法[18]。模擬選用Leap-Frog算法[19]以1 fs的時間步長求解運動方程。由于系統(tǒng)20 ns后達到平衡,故取25~30 ns內(nèi)的模擬結(jié)果作為有效數(shù)據(jù)進行分析,IFT(γ)的定義為界面上正向(z方向)和側(cè)向(x和y方向)壓強張量分量的差異[20]。因系統(tǒng)包括兩個界面,γ可表示為

    式中,Pxx、Pyy、Pzz分別為沿x、y、z方向壓強張量對z向的角分量。

    2 結(jié)果與討論

    2.1 模型驗證

    本文分別模擬了343 K和20 MPa條件下,鹽濃度為 0.9、1.8、2.7 mol·L?1的 CO2-CaCl2和CO2-NaCl鹽水系統(tǒng),IFT增量(Δγ)的計算結(jié)果示于圖2。其中,計算Δγ時所選取的基準為CO2-純水系統(tǒng),該系統(tǒng)在此溫度和壓強下的IFT值為31.8 mN·m?1[8]。由圖2可以看出,CO2-NaCl與CO2-CaCl2溶液系統(tǒng)的 IFT增量值均與鹽濃度呈線性關(guān)系;Ca2+的電荷數(shù)是Na+的2倍,而CO2-NaCl系統(tǒng)的IFT增量值與鹽濃度關(guān)系的實驗值曲線斜率為CO2-CaCl2系統(tǒng)的2.5倍[5];CO2-鹽水系統(tǒng)的IFT模擬結(jié)果與實驗數(shù)據(jù)吻合較好[5],說明本文的計算模型具有一定的準確性和可靠性。

    圖2 CO2-NaCl及CO2-CaCl2系統(tǒng)IFT增量與鹽濃度的關(guān)系Fig.2 Incremental IFTs for CO2-CaCl2and CO2-NaCl solution systems as a function of salinity

    2.2 陽離子原子特性參數(shù)的影響

    鹽離子的原子特性參數(shù)包括:陽離子電荷數(shù)(q)、L-J勢參數(shù)(ε、σ)、原子質(zhì)量(m)等。討論陽離子對IFT增量的影響必然要對以上各因素進行單獨分析。本節(jié)通過前文介紹的虛擬系統(tǒng)進行模擬,探討了各因素對IFT增量的影響及規(guī)律。

    2.2.1 陽離子電荷數(shù)(q)的影響 為將陽離子電荷數(shù)從諸多影響因素中獨立出來,得到僅因電荷數(shù)差異而導致的IFT變化規(guī)律,本文對表2中的虛擬系統(tǒng)3*(Na2+)與系統(tǒng)4*(Ca+Cl?)進行了模擬。在這些虛擬系統(tǒng)中,只改變了陽離子電荷數(shù)而保持其他陽離子特性參數(shù)不變。例如,將 NaCl的鈉離子由Na+變?yōu)镹a2+(系統(tǒng)3*),將CaCl2的鈣離子由Ca2+變?yōu)镃a+(系統(tǒng)4*)。虛擬系統(tǒng)3*、4*的鹽濃度均選為1.8 mol·L?1,計算結(jié)果示于圖3。

    圖3 虛擬及真實系統(tǒng)IFT增量計算結(jié)果Fig.3 Incremental IFTs for fictitious and real systems at salinity of 1.8 mol·L?1

    由圖3陽離子電荷數(shù)比照組可以看出,在鹽濃度為1.8 mol·L?1時,模擬所得到的Na+Cl?、Ca2+、Na2+和Ca+Cl?作為溶質(zhì)系統(tǒng)的Δγ分別為2.3、5.7、5.1和2.9 mN·m?1。由此可見,當鹽濃度一致時,在相同陽離子電荷數(shù)的情況下(系統(tǒng)1與4*、系統(tǒng)2與3*),Δγ值十分接近;二價陽離子Ca2+(系統(tǒng)2)和Na2+(系統(tǒng)3*)所引起的Δγ值約為一價陽離子Na+(系統(tǒng)1)與Ca+(系統(tǒng)4*)的2倍;由此看出,陽離子電荷數(shù)是導致Δγ值變化的主要原因。但是,兩個二價系統(tǒng)(Ca2+和Na2+)之間及兩個一價系統(tǒng)(Na+和Ca+)之間的Δγ值并不相等,均存在0.6 mN·m?1的差異,由此可預測仍有其他原子特性參數(shù)會對Δγ產(chǎn)生影響。

    2.2.2 L-J勢能參數(shù)(ε、σ)的影響 為得到L-J勢能參數(shù)ε和σ對Δγ值的影響,本文另增加了KCl溶液參照系統(tǒng)(表2中系統(tǒng)11)。Na+Cl?(系統(tǒng)1)、Ca2+(系統(tǒng)2)、K+Cl?(系統(tǒng)11)3個真實系統(tǒng)的L-J參數(shù)比較關(guān)系為εK+ < εNa+ < εCa2+和σK+ > σNa+> σCa2+(表1)。計算所得3個真實系統(tǒng)的IFT增量為:ΔγK+ < ΔγNa+ < ΔγCa2+,與實驗數(shù)據(jù)相一致[5,21-22]。但僅從這些真實系統(tǒng)的數(shù)據(jù)仍無法分析ε與σ對Δγ值的單獨影響。因此,本文繼續(xù)研究了系統(tǒng) 5*(Na+Cl?(with εCa2+))、6*(Na+Cl?(with σCa2+))、8*(Na+Cl?(with εK+))、9*(Na+Cl?(with σK+)),計算結(jié)果示于圖3。當保持系統(tǒng)1(Na+Cl?)Na+的電荷數(shù)及σ等原子特性參數(shù)不變,僅將參數(shù)ε改為εCa2+(系統(tǒng)5*)和εK+(系統(tǒng)8*)時,得到ΔγεK+<ΔγεNa+<ΔγεCa2+。改用 εK+時,Δγ相對于 εNa+下降了 0.4 mN·m?1,而用 εCa2+時 Δγ則相對 εNa+上升了 1.1 mN·m?1(見圖3中L-J勢參數(shù)ε比照組)。由此可推論,參數(shù)ε值越大(L-J勢的勢井越深),Δγ值越大。依據(jù) Lorentz-Berthelot混合定律可對此推論進行理論分析,在離子與水中的H原子或O原子之間的ε越大,離子與水分子間的范德華力作用就越強,進而含水電解質(zhì)溶液中分子的結(jié)合更加穩(wěn)固,最終導致了CO2與鹽水界面的IFT值越大。

    同樣,在系統(tǒng)1的基礎(chǔ)上,將Na+的參數(shù)σ改變?yōu)棣褻a2+(系統(tǒng)6*)和σK+(系統(tǒng)9*),可得結(jié)果:ΔγσCa2+< ΔγσNa+ <ΔγσK+。由圖3中L-J勢參數(shù)σ比照組(系統(tǒng)1、6*和9*)可知,改用σCa2+時Δγ相對于σNa+下降了0.6 mN·m?1,而用σK+時則相對σNa+上升了0.1 mN·m?1。由此可預測,σ越大,Δγ值越大,即增大離子半徑將導致IFT增大,但參數(shù)σ對Δγ的影響小于ε。總體來說,參數(shù)ε與σ對Δγ的影響程度均不及陽離子電荷數(shù)對Δγ的主導作用。

    2.2.3 原子質(zhì)量(m)的影響 與討論 L-J參數(shù)方法類似,本文還對虛擬系統(tǒng)7*[Na+Cl?(with mCa2+)]、系統(tǒng)10*[Na+Cl?(with mK+)]與真實系統(tǒng)1(Na+Cl?)(圖3原子質(zhì)量比照組)進行原子質(zhì)量的影響分析。已知原子質(zhì)量的大小關(guān)系為:mNa+ < mCa2+ ≈mK+。由圖 3 的計算結(jié)果可以發(fā)現(xiàn):ΔγmNa+ <ΔγmCa2+≈ΔγmK+,改用mCa2+(系統(tǒng)7*)和mK+(系統(tǒng)10*)的Δγ相對于mNa+(系統(tǒng)1)均下降了0.2 mN·m?1。由此可見,原子質(zhì)量越大,Δγ就越大,但原子質(zhì)量對于Δγ的影響程度相較電荷數(shù)及LJ參數(shù)而言較小。

    2.3 陽離子特性及離子濃度的綜合影響

    由上述討論可知,溶液電解質(zhì)不同時,陽離子的電荷數(shù)是導致 CO2-鹽水系統(tǒng) IFT變化的主導因素;而電解質(zhì)(或是陽離子電荷數(shù))相同時,界面處鹽離子的解吸能力會隨著鹽濃度升高而增強[8],進而導致 IFT增量隨著離子濃度呈現(xiàn)出線性關(guān)系(圖 2)。然而,目前綜合考慮陽離子及離子濃度兩大因素對IFT增量進行定量的相關(guān)性分析尚未見報道。圖2中CO2-CaCl2和 CO2-NaCl系統(tǒng)無論是模擬值亦或是實驗值,Δγ隨鹽濃度變化的曲線斜率比值均約為 2.5,與離子濃度比(1.5:1)或是陽離子電荷數(shù)比(2:1)均不一致。為此本文選取離子強度(I)這一物理量,綜合考慮離子電荷數(shù)及離子濃度因素,對IFT變化進行分析。離子強度的計算公式為

    式中,I為離子強度,ci為離子i的摩爾濃度,zi為離子i的電荷數(shù)。離子強度對IFT的影響研究可能揭示溶液電性強弱對界面結(jié)構(gòu)的影響。

    本文計算得到的CO2-NaCl和CO2-CaCl2溶液系統(tǒng)的Δγ值與離子強度的關(guān)系示于圖4(a)。由圖中可看到,這兩個系統(tǒng)的 Δγ與離子強度的關(guān)系曲線的斜率十分接近,其中CO2-CaCl2溶液系統(tǒng)的斜率略微小于CO2-NaCl系統(tǒng)的斜率。

    為進一步驗證所觀測到的Δγ和I之間的關(guān)系,本文模擬分析了CO2-混合鹽系統(tǒng),分別設置了總鹽濃度為0.9、1.8和2.7 mol·L?1的NaCl和CaCl2混合鹽溶液,并同時改變各溶液中 NaCl的摩爾分數(shù)為 20%、40%、60%和 80%依次進行模擬。不同NaCl/CaCl2(摩爾比)時IFT隨濃度的變化規(guī)律示于圖4(b)。由該圖可知,Δγ和混合鹽的總濃度呈線性關(guān)系,且混合鹽系統(tǒng)的 Δγ和鹽濃度的斜率介于純NaCl和純CaCl2系統(tǒng)之間。

    本文通過分析IFT斜率與兩種鹽摩爾比之間的關(guān)系,發(fā)現(xiàn)混合鹽引起的IFT斜率δmix可用插值法進行計算

    圖4 CO2-鹽水系統(tǒng)IFT增量與離子強度、鹽濃度的關(guān)系Fig.4 Incremental IFTs for CO2-salt solution systems as functions of ionic strength and salinities

    式中,xi為鹽類i占混合鹽的比例,δi為純鹽i引起的IFT增加斜率。式(4)可以為預測混合鹽的IFT提供依據(jù)。

    計算所得Δγ隨著NaCl組分變化的關(guān)系示于圖5(a)。由圖可知,Δγ與 NaCl組分的關(guān)系曲線斜率隨著混合鹽總濃度的上升而下降;在 NaCl組分固定時,Δγ隨著混合鹽總濃度呈現(xiàn)出與實驗數(shù)據(jù)[5]一致的線性變化。

    本文對圖5(a)的數(shù)據(jù)進一步處理,縱坐標改變?yōu)閱挝浑x子強度下的Δγ,結(jié)果示于圖5(b)。由圖中可以看出,在不同混合鹽總濃度情況下,單位離子強度的Δγ隨NaCl摩爾組分幾乎不變。由此可知,在強電解質(zhì)溶液中,將離子產(chǎn)生的電場影響進行量化并取單位值之后,Δγ則與溶液中鹽的成分(NaCl/CaCl2)無明顯關(guān)聯(lián),該現(xiàn)象為混合鹽系統(tǒng)IFT的預測提供了重要依據(jù)。

    根據(jù)以上研究,本文進一步預測:CO2-鹽水系統(tǒng)的單位離子強度的 Δγ的大小可能與鹽溶液電解質(zhì)的類型關(guān)系不大,在某一確定環(huán)境條件下,Δγ/I極可能對于通用型礦物電解質(zhì)(例如本文所考慮的Na+、Ca2+、K+和 Cl?所構(gòu)成的鹽)為一定值,可為工程上預測不同鹽溶液的類型、組分及濃度下的IFT值提供指導。

    3 結(jié) 論

    本文采用分子動力學模擬的方法對343 K和20 MPa條件下不同鹽濃度的CO2-NaCl、CO2-CaCl2和CO2-(NaCl+CaCl2)混合溶液系統(tǒng)進行了鹽離子對IFT增量的影響研究,得到以下結(jié)論。

    (1)IFT增量的變化主要來自于陽離子電荷數(shù),而L-J勢能力場參數(shù)ε、σ及原子質(zhì)量m也對IFT增量產(chǎn)生影響,但相對于電荷數(shù)影響較小。

    (2)CO2-NaCl和CO2-CaCl2界面的IFT增量隨離子強度成線性變化且斜率相近,溶液中離子強度對IFT的影響是具備量化條件的。

    (3)不同濃度及混合鹽摩爾比下單位離子強度IFT基本趨于定值,即CO2-鹽水系統(tǒng)中Δγ直接與I相關(guān),而與鹽溶液電解質(zhì)的類型和組分無關(guān),可為工程上預測IFT值提供指導。

    圖5 CO2-混合鹽系統(tǒng)IFT增量及單位離子強度IFT增量Fig.5 Incremental IFTs and incremental IFTs per ionic strength for CO2- mixed salts solution systems as functions of NaCl molar compositions

    符 號 說 明

    I ——離子強度,mol·L?1

    m ——原子質(zhì)量

    Pxx——x方向壓強張量對z向的角分量

    Pyy——y方向壓強張量對z向的角分量

    Pzz——z方向壓強張量對z向的角分量

    q ——離子的電荷數(shù),e

    rij——原子i和j之間的距離,nm

    xi——鹽類i占混合鹽的比例

    γ——界面張力,mN·m?1

    Δγ——界面張力增量,mN·m?1

    δi——純鹽i引起的IFT增加斜率混合鹽引起的IFT斜率

    ε——L-J勢能曲線中的勢能阱深度,nm

    σ——勢能為零時,兩點間的距離,nm

    [1] ECCLES J K, PRATSON L, NEWELL, R G, et al. Physical and economic potential of geological CO2storage in saline aquifers[J]. Environ. Sci. Technol., 2009, 43(6): 1962-1969.

    [2] SZULCZEWSKI M L, MACMINN C W, HERZOG H J, et al. Lifetime of carbon capture and storage as a climate-change mitigation technology[J]. Proc. Natl. Acad. Sci., 2012, 109(14): 5185-5189.

    [3]BACHU S, BENNION D B. Interfacial tension between CO2, freshwater, and brine in the range of pressure from (2 to 27) MPa, temperature from (20 to 125)℃, and water salinity from (0 to 334000) mg·L?1[J]. J. Chem. Eng. Data, 2009, 54(3): 765-775.

    [4] CHALBAUD C, ROBIN M, LOMBARD J M, et al. Interfacial tension measurements and wettability evaluation for geological CO2storage[J]. Adv.Water Resour., 2009, 32(1): 98-109.

    [5]AGGELOPOULOS C A, ROBIN M, PERFETTI E, et al. CO2/CaCl2solution interfacial tensions under CO2geological storage conditions: influence of cation valence on interfacial tension[J]. Adv. Water Resour., 2010, 33(6): 691-697.

    [6] DA ROCHA S R P, JOHNSTON K P, WESTACOTT R E, et al. Molecular structure of the water-supercritical CO2interface[J]. J. Phys. Chem. B, 2001, 105(48): 12092-12104.

    [7] LI X S, ROSS D A, TRUSLER J P, et al. Molecular dynamics simulations of CO2and brine interfacial tension at high temperatures and pressures[J]. J. Chem. Eng. Data, 2013, 117(18), 5647-5652.

    [8] ZHAO L L, LIN S C, MENDENHALL J D, et al. Molecular dynamics investigation of the various atomic force contributions to the interfacial tension at the supercritical CO2-water interface[J]. J. Phys. Chem. B, 2011, 115(19): 6076-6087.

    [9] 范康年. 物理化學[M]. 2版. 北京:高等教育出版社,2005:222. FAN K N. Physical Chemistry [M]. 2nd ed. Beijing: Higher Education Press, 2005: 222.

    [10]DARDEN T, YORK D, PEDERSEN L. Particle Mesh Ewald: an N·log (N) method for Ewald Sums in large systems[J]. J. Chem. Phys., 1993, 98(12): 10089-10092.

    [11]DE-LARA L S, MICHELON M F, MIRANDA C R. Molecular dynamics studies of fluid/oil interfaces for improved oil recovery processes[J]. J. Phys. Chem. B, 2012, 116(50): 14667-146676.

    [12]LEVITT M, HIRSHBERG M, SHARON R, et al. Calibration and testing of a water model for simulation of the molecular dynamics of proteins and nucleic acids in solution[J]. J. Phys. Chem. B, 1997, 101(25): 5051-5061.

    [13]NIETO-DRAGHI C, DE BRUIN T, PEREZ-PELLITERO J, et al. Thermodynamic and transport properties of carbon dioxide from molecular simulation[J]. J. Chem. Phys., 2007, 126(6): 064509.

    [14]CHANDRASEKHAR J, SPELLMEYER D C, JORGENSEN W L.Energy component analysis for dilute aqueous solutions of lithium (1+), sodium (1+), fluoride (1?), and chloride (1?) ions[J]. J. Am. Chem. Soc., 1984, 106(4): 903-910.

    [15]AQVIST J. Ion-water interaction potentials derived from free energy perturbation simulations[J]. J. Chem. Phys., 1990, 94(21): 8021-8024.

    [16]PRONK S, PALL S, SCHULZ R, et al. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit[J]. Bioinformatics, 2013, 29(7): 845-854.

    [17]BUSSI G, DONADIO D, PARRINELLO M. Canonical sampling through velocity rescaling[J]. J. Chem. Phys., 2007, 126(1): 014101.

    [18]BERENDSEN H J, POSTMA J P M, VAN GUNSTEREN W F, et al. Molecular dynamics with coupling to an external bath[J]. J. Chem. Phys., 1984, 81(8): 3684-3690.

    [19]HOCKNEY R, GOEL S, EASTWOOD J. Quiet high-resolution computer models of a plasma[J]. J. Comput. Phys., 1974, 14(2): 148-158.

    [20]ALEJANDRE J, TILDESLEY D J, CHAPELA G A. Molecular dynamics simulation of the orthobaric densities and surface tension of water[J]. J. Chem. Phys., 1995, 102(11): 4574-4583.

    [21]LI X S, BOEK E S, MAITLAND G C, et al. Interfacial tension of (brines + CO2): CaCl2(aq), MgCl2(aq), and Na2SO4(aq) at temperatures between (343 and 423) K, pressures between (2 and 50) MPa, and molalities of (0.5 to 5) mol·kg?1[J]. J. Chem. Eng. Data, 2012, 57(5): 1369-1375.

    [22]PEGRAM L M, RECORD M T J. Hofmeister salt effects on surface tension arise from partitioning of anions and cations between bulk water and the air-water interface[J]. J. Phys. Chem. B, 2007, 111(9): 5411-5417.

    Ion effect on interfacial properties of CO2-brine system

    JI Jiayuan, ZHAO Lingling, LI Siyu
    (Key Laboratory of Energy Thermal Conversion and Control of Ministry of Education, School of Energy & Environment, Southeast University, Nanjing 210096, Jiangsu, China)

    Interfacial tension of CO2-NaCl, CO2-CaCl2and CO2-(NaCl+CaCl2) solution systems was studied at different salt concentration and ratio by molecular dynamics simulation, in order to understand influence of different salt ions at molecular level. The results showed that the presence of salt ions increased interfacial tension between CO2and water, which majority was caused by number of cationic charges and the rest was by cation atomic properties such as Lennard Jones potential parameters and atomic mass. The interfacial tension increment of these three systems were all linearly dependent on ionic strength, but the slope had no relationship to types of salt solutions. This founding would provide theoretical guidance for IFT prediction.

    CO2-brine system; interfacial tension; molecular simulation; atomic characteristic parameter; ionic strength

    Prof. ZHAO Lingling, zhao_lingling@seu.edu.cn

    TQ 021.2

    :A

    :0438—1157(2017)01—0223—07

    10.11949/j.issn.0438-1157.20160771

    2016-06-02收到初稿,2016-09-30收到修改稿。

    聯(lián)系人:趙伶玲。

    :季佳圓(1993—),女,碩士研究生。

    國家自然科學基金項目(51106027)。

    Received date: 2016-06-02.

    Foundation item: supported by the National Natural Science Foundation of China (51106027).

    猜你喜歡
    原子質(zhì)量鹽濃度陽離子
    再識相對原子質(zhì)量
    2023年高考化學模擬試題A卷
    不同鹽分條件下硅對兩個高羊茅品種生物量分配和營養(yǎng)元素氮、磷、鉀吸收利用的影響
    高考化學模擬試題
    混合鹽堿脅迫對醉馬草種子萌發(fā)及幼苗生理特性的影響
    烷基胺插層蒙脫土的陽離子交換容量研究
    硅對不同抗性高羊茅耐鹽性的影響
    主鹽濃度對Au-Pt 合金催化劑性能的影響
    陽離子Gemini表面活性劑的應用研究進展
    仿生陽離子多烯環(huán)化
    eeuss影院久久| 中文亚洲av片在线观看爽| 精品99又大又爽又粗少妇毛片| 午夜精品在线福利| 青青草视频在线视频观看| 亚洲av成人av| 亚洲国产精品sss在线观看| 国产精品.久久久| 午夜视频国产福利| 亚洲精品影视一区二区三区av| 国产又色又爽无遮挡免| 欧美丝袜亚洲另类| 观看美女的网站| 国产伦精品一区二区三区视频9| ponron亚洲| h日本视频在线播放| 欧美97在线视频| av线在线观看网站| 乱系列少妇在线播放| 国产成人a∨麻豆精品| 噜噜噜噜噜久久久久久91| 99久久精品热视频| av在线观看视频网站免费| 女人久久www免费人成看片 | 女人久久www免费人成看片 | 日韩欧美三级三区| 亚洲自拍偷在线| 1000部很黄的大片| 三级经典国产精品| 麻豆一二三区av精品| 欧美成人午夜免费资源| 51国产日韩欧美| 两个人视频免费观看高清| 国产欧美另类精品又又久久亚洲欧美| 国产 一区精品| 国产精品久久久久久久电影| 国产伦一二天堂av在线观看| 十八禁国产超污无遮挡网站| 久久人人爽人人片av| 2021少妇久久久久久久久久久| 欧美激情在线99| 日韩强制内射视频| АⅤ资源中文在线天堂| 色噜噜av男人的天堂激情| 国产高潮美女av| 日本黄色视频三级网站网址| 看免费成人av毛片| 小蜜桃在线观看免费完整版高清| 日本三级黄在线观看| 亚洲欧美中文字幕日韩二区| 久久6这里有精品| 一级毛片电影观看 | 美女xxoo啪啪120秒动态图| 在线免费观看的www视频| 黄片wwwwww| 欧美又色又爽又黄视频| 你懂的网址亚洲精品在线观看 | 九色成人免费人妻av| 亚洲自拍偷在线| 久久久午夜欧美精品| av在线蜜桃| 美女脱内裤让男人舔精品视频| 中文资源天堂在线| 联通29元200g的流量卡| 久久久久久久亚洲中文字幕| 久久久久久久久久久丰满| 蜜桃久久精品国产亚洲av| 麻豆乱淫一区二区| 好男人在线观看高清免费视频| 一区二区三区乱码不卡18| 成人午夜高清在线视频| 国产成人福利小说| 日韩欧美三级三区| 亚洲18禁久久av| 亚洲精品乱久久久久久| 成人国产麻豆网| 国内揄拍国产精品人妻在线| av国产免费在线观看| 国产精品国产三级专区第一集| 纵有疾风起免费观看全集完整版 | 成人亚洲精品av一区二区| 99热这里只有是精品在线观看| 美女cb高潮喷水在线观看| 亚洲国产精品成人综合色| 美女xxoo啪啪120秒动态图| 成年女人永久免费观看视频| 免费av毛片视频| a级毛片免费高清观看在线播放| 精品久久久久久电影网 | 免费观看的影片在线观看| 亚洲精品国产av成人精品| 亚洲国产精品成人综合色| 国产成人a∨麻豆精品| 免费搜索国产男女视频| 一级二级三级毛片免费看| 男女边吃奶边做爰视频| 亚州av有码| 亚州av有码| 国产成年人精品一区二区| 久久久久精品久久久久真实原创| 人人妻人人澡欧美一区二区| 精品国产一区二区三区久久久樱花 | 99热精品在线国产| 久久久久久久亚洲中文字幕| 日本一本二区三区精品| 亚洲av不卡在线观看| 欧美极品一区二区三区四区| 久久精品久久久久久久性| 最近的中文字幕免费完整| 日本一二三区视频观看| 能在线免费看毛片的网站| 一二三四中文在线观看免费高清| 免费看美女性在线毛片视频| 欧美一区二区国产精品久久精品| 男女国产视频网站| 国产国拍精品亚洲av在线观看| 国产欧美另类精品又又久久亚洲欧美| 亚洲自偷自拍三级| 精品午夜福利在线看| 国产视频首页在线观看| 色哟哟·www| 亚洲内射少妇av| 亚洲欧洲日产国产| 久久久久久久午夜电影| 99在线人妻在线中文字幕| АⅤ资源中文在线天堂| 国国产精品蜜臀av免费| 亚洲国产欧美人成| 成人特级av手机在线观看| 一级毛片久久久久久久久女| 美女xxoo啪啪120秒动态图| 天美传媒精品一区二区| 天美传媒精品一区二区| 两个人视频免费观看高清| 亚洲高清免费不卡视频| 亚洲五月天丁香| 午夜免费激情av| 99在线视频只有这里精品首页| 少妇人妻精品综合一区二区| 国产美女午夜福利| 亚洲国产色片| 久久久久性生活片| 99久国产av精品| 欧美丝袜亚洲另类| 中文乱码字字幕精品一区二区三区 | 亚洲精品乱码久久久v下载方式| 18禁动态无遮挡网站| 99久久成人亚洲精品观看| 中文在线观看免费www的网站| 亚洲天堂国产精品一区在线| 久久久精品大字幕| 国产在线男女| 亚洲最大成人中文| 亚洲成色77777| 久久精品91蜜桃| 午夜福利网站1000一区二区三区| 国产精品无大码| 国内揄拍国产精品人妻在线| 99久久中文字幕三级久久日本| 美女黄网站色视频| 综合色av麻豆| 亚洲精品影视一区二区三区av| 亚洲精品影视一区二区三区av| 久久久色成人| 成人av在线播放网站| 99热6这里只有精品| 少妇熟女欧美另类| 成人二区视频| 国产精品永久免费网站| 成人美女网站在线观看视频| 国产伦精品一区二区三区四那| 搞女人的毛片| 欧美区成人在线视频| 三级国产精品欧美在线观看| 国产亚洲5aaaaa淫片| 久久久久精品久久久久真实原创| 国产精品人妻久久久影院| 中文字幕熟女人妻在线| 免费看美女性在线毛片视频| 老司机影院成人| 99热这里只有是精品50| 小说图片视频综合网站| 性色avwww在线观看| 亚洲婷婷狠狠爱综合网| 简卡轻食公司| 麻豆av噜噜一区二区三区| 精品久久久久久成人av| 真实男女啪啪啪动态图| 午夜免费男女啪啪视频观看| 岛国毛片在线播放| 日本黄大片高清| 人人妻人人澡欧美一区二区| www日本黄色视频网| 高清日韩中文字幕在线| 天堂√8在线中文| 一个人免费在线观看电影| 91精品伊人久久大香线蕉| 亚洲欧美日韩高清专用| 韩国高清视频一区二区三区| 国产91av在线免费观看| 国产 一区 欧美 日韩| 精品久久久久久久人妻蜜臀av| 色5月婷婷丁香| av天堂中文字幕网| 日本五十路高清| 午夜a级毛片| 五月玫瑰六月丁香| 久久久久久大精品| 欧美激情在线99| 久久国内精品自在自线图片| 在线免费十八禁| 99久国产av精品| 亚洲av成人精品一区久久| 精品久久久久久久人妻蜜臀av| 成年版毛片免费区| 永久免费av网站大全| 日韩中字成人| 亚洲国产欧美人成| 国产黄片视频在线免费观看| 能在线免费看毛片的网站| 最近2019中文字幕mv第一页| 晚上一个人看的免费电影| 亚洲在线观看片| 国产探花极品一区二区| 中文字幕av成人在线电影| 黄片wwwwww| 97在线视频观看| 卡戴珊不雅视频在线播放| 久久久久网色| 免费黄网站久久成人精品| 中文字幕人妻熟人妻熟丝袜美| 国产伦在线观看视频一区| av天堂中文字幕网| 汤姆久久久久久久影院中文字幕 | 联通29元200g的流量卡| 亚洲,欧美,日韩| www.av在线官网国产| 全区人妻精品视频| 国产精品福利在线免费观看| 七月丁香在线播放| 久久久久久国产a免费观看| 精品久久久久久久末码| 国产美女午夜福利| 99久久人妻综合| 日本熟妇午夜| 国产精品一区二区性色av| 美女大奶头视频| av又黄又爽大尺度在线免费看 | 国产精品一二三区在线看| 大又大粗又爽又黄少妇毛片口| 高清视频免费观看一区二区 | 国产精品.久久久| 成人毛片a级毛片在线播放| 久久久久精品久久久久真实原创| 床上黄色一级片| 亚洲aⅴ乱码一区二区在线播放| 免费黄色在线免费观看| 99热这里只有精品一区| 国产在线男女| 精品不卡国产一区二区三区| 亚洲国产精品sss在线观看| 色综合站精品国产| 亚洲五月天丁香| 青春草国产在线视频| 国产色婷婷99| 色综合色国产| 日韩av在线免费看完整版不卡| 国产伦精品一区二区三区视频9| 久久久精品94久久精品| 性插视频无遮挡在线免费观看| 非洲黑人性xxxx精品又粗又长| 尤物成人国产欧美一区二区三区| 日韩欧美在线乱码| 最近中文字幕2019免费版| 人人妻人人澡欧美一区二区| 国产一区亚洲一区在线观看| 美女内射精品一级片tv| 美女高潮的动态| a级一级毛片免费在线观看| 全区人妻精品视频| 91av网一区二区| 亚洲精品国产成人久久av| 午夜福利视频1000在线观看| 精品久久久久久成人av| 亚洲无线观看免费| 寂寞人妻少妇视频99o| 亚洲怡红院男人天堂| 嫩草影院新地址| 一卡2卡三卡四卡精品乱码亚洲| 熟妇人妻久久中文字幕3abv| 国产免费福利视频在线观看| 热99在线观看视频| 亚洲欧美日韩卡通动漫| 久久精品久久久久久噜噜老黄 | 免费黄网站久久成人精品| 日韩欧美在线乱码| 亚洲欧美精品自产自拍| 精品国产三级普通话版| 国国产精品蜜臀av免费| 午夜免费激情av| 建设人人有责人人尽责人人享有的 | 欧美日本亚洲视频在线播放| 亚洲av成人av| 毛片一级片免费看久久久久| 国产视频内射| 久久精品熟女亚洲av麻豆精品 | 三级国产精品片| 亚洲在线自拍视频| 欧美区成人在线视频| 精品久久久久久久久亚洲| 久久久久久久午夜电影| 在线免费十八禁| 成人综合一区亚洲| 国内精品美女久久久久久| 亚洲av成人精品一二三区| 搡老妇女老女人老熟妇| 少妇熟女欧美另类| 精品国内亚洲2022精品成人| 中文亚洲av片在线观看爽| 亚洲av成人精品一二三区| 免费观看人在逋| 三级国产精品片| 中文资源天堂在线| 久久欧美精品欧美久久欧美| 久久人人爽人人爽人人片va| 国产在视频线精品| 九九爱精品视频在线观看| 国产伦一二天堂av在线观看| 日韩av不卡免费在线播放| 天堂√8在线中文| 国产亚洲av嫩草精品影院| 国产乱人偷精品视频| 久久鲁丝午夜福利片| 亚洲自拍偷在线| 99九九线精品视频在线观看视频| av女优亚洲男人天堂| 欧美97在线视频| 亚洲欧美精品自产自拍| 亚洲综合色惰| 亚洲国产最新在线播放| 老女人水多毛片| 日本黄色视频三级网站网址| 美女国产视频在线观看| 久久久久九九精品影院| 国产综合懂色| 精品国产露脸久久av麻豆 | 日韩国内少妇激情av| 久久99蜜桃精品久久| 久久精品综合一区二区三区| 99在线视频只有这里精品首页| 国产 一区 欧美 日韩| 亚洲电影在线观看av| 国产精品一及| 亚洲婷婷狠狠爱综合网| 国产在视频线在精品| 全区人妻精品视频| 黑人高潮一二区| 日韩亚洲欧美综合| 亚洲丝袜综合中文字幕| av在线蜜桃| 亚洲国产欧美人成| 亚洲人与动物交配视频| 看黄色毛片网站| 久久精品人妻少妇| 久久久久久大精品| 久久精品夜色国产| 免费av观看视频| 日韩一本色道免费dvd| 精品久久国产蜜桃| 狂野欧美白嫩少妇大欣赏| 国产精品不卡视频一区二区| 国产av在哪里看| 中国国产av一级| 国产单亲对白刺激| 亚洲最大成人中文| 精品久久久久久久久亚洲| 日韩成人av中文字幕在线观看| 成人二区视频| 一个人看视频在线观看www免费| 中文欧美无线码| 日韩成人av中文字幕在线观看| 久久精品国产99精品国产亚洲性色| 色5月婷婷丁香| 日本五十路高清| 深爱激情五月婷婷| 亚洲av不卡在线观看| 国产精品一区二区三区四区免费观看| 久久婷婷人人爽人人干人人爱| 99久久精品国产国产毛片| av在线老鸭窝| 国产一级毛片七仙女欲春2| 精品午夜福利在线看| 99热6这里只有精品| 国产午夜精品论理片| 日韩中字成人| 69av精品久久久久久| 高清毛片免费看| 国产av码专区亚洲av| 精品午夜福利在线看| 好男人在线观看高清免费视频| 国产淫语在线视频| 亚洲乱码一区二区免费版| 视频中文字幕在线观看| 国产黄色视频一区二区在线观看 | 日本免费一区二区三区高清不卡| 色哟哟·www| 国产亚洲一区二区精品| 少妇熟女aⅴ在线视频| 一夜夜www| 欧美成人午夜免费资源| 一个人看的www免费观看视频| 亚洲怡红院男人天堂| 国产精品麻豆人妻色哟哟久久 | 中文字幕人妻熟人妻熟丝袜美| 国产极品精品免费视频能看的| 中国国产av一级| 91精品国产九色| 亚洲成色77777| 亚洲熟妇中文字幕五十中出| 免费人成在线观看视频色| 亚洲人成网站在线播| 国产伦在线观看视频一区| 2022亚洲国产成人精品| 国产午夜精品久久久久久一区二区三区| 九九久久精品国产亚洲av麻豆| 色哟哟·www| 一二三四中文在线观看免费高清| 最近2019中文字幕mv第一页| 十八禁国产超污无遮挡网站| 国产爱豆传媒在线观看| 久久婷婷人人爽人人干人人爱| 日本av手机在线免费观看| 韩国av在线不卡| 99热精品在线国产| 亚洲精品aⅴ在线观看| 五月玫瑰六月丁香| 一本久久精品| 18禁动态无遮挡网站| 中文字幕久久专区| 插逼视频在线观看| 国内揄拍国产精品人妻在线| 免费大片18禁| 又爽又黄a免费视频| 日韩av不卡免费在线播放| 国模一区二区三区四区视频| 日韩欧美国产在线观看| 国产精品一区二区性色av| 中文字幕熟女人妻在线| 久久人人爽人人片av| .国产精品久久| 少妇丰满av| 午夜爱爱视频在线播放| 欧美日韩综合久久久久久| 国产精品女同一区二区软件| 久久久a久久爽久久v久久| 亚洲精品,欧美精品| 99在线视频只有这里精品首页| 波野结衣二区三区在线| 国产精品三级大全| 久久久精品大字幕| 国内揄拍国产精品人妻在线| 久久精品夜色国产| 国产精品,欧美在线| 亚洲av免费高清在线观看| 国产精品一区二区三区四区免费观看| 亚洲精品国产成人久久av| 老女人水多毛片| 3wmmmm亚洲av在线观看| a级毛色黄片| 亚洲中文字幕一区二区三区有码在线看| 亚洲高清免费不卡视频| 国产免费男女视频| 91久久精品国产一区二区三区| 国产v大片淫在线免费观看| 中文天堂在线官网| 美女内射精品一级片tv| 精品酒店卫生间| 一级毛片久久久久久久久女| 大又大粗又爽又黄少妇毛片口| 高清日韩中文字幕在线| 午夜精品国产一区二区电影 | 国产真实伦视频高清在线观看| 国产淫片久久久久久久久| 一级黄色大片毛片| 综合色丁香网| 久久久欧美国产精品| 七月丁香在线播放| 免费看av在线观看网站| 中文天堂在线官网| 毛片一级片免费看久久久久| 国产免费男女视频| 成人特级av手机在线观看| 又黄又爽又刺激的免费视频.| 亚洲精品成人久久久久久| 非洲黑人性xxxx精品又粗又长| 午夜福利在线观看免费完整高清在| 亚洲国产精品成人综合色| 日韩国内少妇激情av| 淫秽高清视频在线观看| 欧美97在线视频| 国产精品一及| 99热这里只有精品一区| 成人三级黄色视频| 日本黄色片子视频| 边亲边吃奶的免费视频| 在现免费观看毛片| 99视频精品全部免费 在线| videossex国产| 波野结衣二区三区在线| 一边摸一边抽搐一进一小说| 成人毛片60女人毛片免费| 欧美丝袜亚洲另类| 亚洲欧美成人综合另类久久久 | 综合色av麻豆| 大香蕉久久网| 国产精华一区二区三区| av在线观看视频网站免费| 午夜精品国产一区二区电影 | 国产亚洲av嫩草精品影院| h日本视频在线播放| 我要看日韩黄色一级片| 午夜久久久久精精品| 欧美日韩在线观看h| 日本欧美国产在线视频| 天天一区二区日本电影三级| 欧美精品国产亚洲| 午夜福利成人在线免费观看| 麻豆久久精品国产亚洲av| 美女cb高潮喷水在线观看| 一个人看视频在线观看www免费| 欧美精品国产亚洲| 国产大屁股一区二区在线视频| 少妇人妻一区二区三区视频| 1024手机看黄色片| 国产免费福利视频在线观看| 国产免费视频播放在线视频 | kizo精华| 国产精品国产三级国产av玫瑰| 午夜福利成人在线免费观看| 国产三级中文精品| 国产片特级美女逼逼视频| 国产精品一区二区在线观看99 | www.色视频.com| 国产在线一区二区三区精 | 超碰av人人做人人爽久久| 成人午夜精彩视频在线观看| 日韩三级伦理在线观看| 久久久久久久久久久丰满| 亚洲av成人精品一二三区| 内地一区二区视频在线| 搡女人真爽免费视频火全软件| 中文字幕亚洲精品专区| 日本色播在线视频| 久久久国产成人精品二区| 亚洲欧美精品专区久久| 超碰97精品在线观看| 亚洲自拍偷在线| 91精品国产九色| 日韩中字成人| 国产精品久久久久久久久免| 男人舔女人下体高潮全视频| 亚洲av不卡在线观看| 亚洲欧洲国产日韩| 99久国产av精品| 国内揄拍国产精品人妻在线| 国产在视频线在精品| 日韩三级伦理在线观看| 美女被艹到高潮喷水动态| 我要搜黄色片| 久久精品91蜜桃| 色视频www国产| 日韩欧美三级三区| 内射极品少妇av片p| 好男人在线观看高清免费视频| 国产真实伦视频高清在线观看| 又黄又爽又刺激的免费视频.| 免费观看性生交大片5| 精品一区二区三区人妻视频| 中文字幕亚洲精品专区| 国产亚洲午夜精品一区二区久久 | 亚洲人成网站在线观看播放| 男女视频在线观看网站免费| 在线观看av片永久免费下载| 国产在线一区二区三区精 | 大又大粗又爽又黄少妇毛片口| 99国产精品一区二区蜜桃av| ponron亚洲| 国产91av在线免费观看| 国产高清国产精品国产三级 | 色视频www国产| 国产欧美另类精品又又久久亚洲欧美| 国产免费男女视频| 乱系列少妇在线播放| 少妇被粗大猛烈的视频| 国产精品电影一区二区三区| 亚洲av熟女| 男女啪啪激烈高潮av片| av黄色大香蕉| 免费看美女性在线毛片视频| 永久网站在线| 成人亚洲精品av一区二区| 色综合站精品国产| 网址你懂的国产日韩在线| 51国产日韩欧美| 国产精品国产三级国产av玫瑰| 日韩国内少妇激情av| 国产黄色视频一区二区在线观看 | 国产在线男女| 18禁裸乳无遮挡免费网站照片| 2021天堂中文幕一二区在线观| 啦啦啦观看免费观看视频高清| 欧美激情在线99| 男人舔奶头视频| 亚洲av成人精品一区久久| 69av精品久久久久久|