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

    擴(kuò)散模型和凝聚模型耦合作用下膠體凝聚動(dòng)力學(xué)的Monte Carlo模擬研究

    2014-10-18 05:27:48熊海靈楊志敏
    物理化學(xué)學(xué)報(bào) 2014年3期
    關(guān)鍵詞:體系模型

    熊海靈 楊志敏 李 航

    (1西南大學(xué)計(jì)算機(jī)與信息科學(xué)學(xué)院,重慶 400715;2西南大學(xué)土壤多尺度界面過程與調(diào)控重慶市重點(diǎn)實(shí)驗(yàn)室,重慶 400715;3西南大學(xué)資源環(huán)境學(xué)院,重慶 400715;4西南大學(xué)三峽庫(kù)區(qū)生態(tài)環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,重慶 400715)

    1 引言

    膠體穩(wěn)定性的研究在膠體科學(xué)中一直占有很重要的地位,Smoluchowski方程對(duì)膠體凝聚動(dòng)力學(xué)的描述是建立在平均場(chǎng)近似的基礎(chǔ)上的.實(shí)驗(yàn)上可測(cè)量膠體凝聚的速率常數(shù),但無(wú)法據(jù)此分析團(tuán)簇大小對(duì)速率常數(shù)的影響.膠體穩(wěn)定性理論(Derjaguin-Landau-Verwey-Overbeek(DLVO))雖可對(duì)膠體凝聚速率常數(shù)作定性解釋,但仍難以就團(tuán)簇結(jié)構(gòu)對(duì)速率常數(shù)的影響規(guī)律作出定量處理.1

    Monte Carlo模擬方法是近年來發(fā)展起來的一種可用于研究膠體顆粒相互作用的計(jì)算機(jī)模擬方法.2,3目前膠粒凝聚的Monte Carlo模擬最活躍的領(lǐng)域主要在四個(gè)方面,4即:膠體理論的進(jìn)展及計(jì)算方法的改進(jìn),熱力學(xué)性質(zhì)的模擬,動(dòng)力學(xué)性質(zhì)的模擬,團(tuán)簇結(jié)構(gòu)性質(zhì)和形貌特征的模擬.動(dòng)力學(xué)性質(zhì)的模擬包括使用Smoluchowski方程對(duì)體系中每一個(gè)顆粒與其它顆??赡苓M(jìn)行的“反應(yīng)”都作微觀模擬,5也包括了其它形式的模擬.6-11由于體系的復(fù)雜性,基于Smoluchowski方程的模型雖然有微觀機(jī)理方面的正確性,但其準(zhǔn)確性不足;其它形式的模型雖然準(zhǔn)確性較好,但普遍適應(yīng)性欠佳,模型的物理意義也不明確.目前還沒有一個(gè)模型既能明確地表達(dá)體系的微觀機(jī)理,又有宏觀上的普遍適應(yīng)性.4分析原因主要有以下兩點(diǎn):(1)基于Smoluchowski方程的模型需要建立在6條假設(shè)1的基礎(chǔ)上,這種基于微觀表象強(qiáng)加于模型上的約束條件,嚴(yán)重影響了它在真實(shí)系統(tǒng)中的應(yīng)用,其中對(duì)所有碰撞都導(dǎo)致粘結(jié),所有顆粒在碰撞前后均為實(shí)心球體的假設(shè)是該模型與實(shí)際相差太大的最根本原因;(2)其它形式的模型2,6-8,12-15大多還是從Smoluchowski方程的模型出發(fā),孤立地解決那6條假設(shè)所帶來的問題.如有的模型考慮到凝聚過程中形成的團(tuán)簇不是一個(gè)實(shí)心的、規(guī)則的球體,而人為地設(shè)定一個(gè)固定的形狀因子,顯然這樣一個(gè)形狀因子難以描繪團(tuán)簇的實(shí)際特征,也就更難以確定團(tuán)簇的擴(kuò)散機(jī)理.有的模型為了解決一碰撞就粘結(jié)的問題,也是人為地設(shè)定一個(gè)較小的凝聚概率,顯然這也無(wú)法表達(dá)不斷發(fā)展變化的各種大小的團(tuán)簇的凝聚機(jī)理.

    實(shí)際上,膠體的凝聚動(dòng)力學(xué)行為及團(tuán)簇的形貌特征是直接受團(tuán)簇的擴(kuò)散機(jī)理以及凝聚機(jī)理影響的,只有確定了正確的擴(kuò)散模型和凝聚模型,模擬的結(jié)果才可能準(zhǔn)確,模型也才能有較好的適應(yīng)性.

    目前在眾多模型中,7,8,15,16團(tuán)簇-團(tuán)簇凝聚(CCA)模型可以較好地用來描述膠粒凝聚的過程.17即該模型可以用來描述膠粒先聚合成小的結(jié)構(gòu)(團(tuán)簇),這些小的結(jié)構(gòu)再經(jīng)多次聚合而形成越來越大的結(jié)構(gòu)體的過程.以CCA模型為基礎(chǔ),膠粒凝聚過程一般認(rèn)為可分為兩類:擴(kuò)散控制過程的擴(kuò)散置限團(tuán)簇凝聚(DLCA)與反應(yīng)控制過程的反應(yīng)置限團(tuán)簇凝聚(RLCA).18-20在擴(kuò)散控制的凝聚中,兩個(gè)粒子(膠體顆粒或團(tuán)簇)一旦相遇即發(fā)生凝聚.所以,膠體顆粒的凝聚速度由碰撞概率決定,一切影響顆粒移動(dòng)速度的因素都影響碰撞概率,從而影響凝聚速度.而在反應(yīng)控制過程中相碰撞的兩個(gè)膠體顆?;驁F(tuán)簇不是一碰撞就凝聚,而是以一定的、概率值小于1的概率發(fā)生凝聚,這個(gè)概率可被稱為有效碰撞概率.此時(shí)膠體凝聚速度不僅受碰撞概率所制約,而且受有效碰撞概率所制約.由此可以看出,由于RLCA過程的凝聚概率等于碰撞概率與有效碰撞概率的乘積,所以相對(duì)來說RLCA過程是一種慢速進(jìn)行的凝聚過程.不同體系及同一體系不同條件下的碰撞概率和有效碰撞概率具有較大的差異.21

    然而,問題的關(guān)鍵是如何確定碰撞概率和有效碰撞概率,也就是說在利用CCA模型進(jìn)行計(jì)算機(jī)模擬時(shí)如何確定擴(kuò)散模型和凝聚模型.這兩個(gè)模型都要求模擬過程既能聯(lián)系諸多影響膠粒凝聚的熱力學(xué)條件,又要能體現(xiàn)膠體科學(xué)中已經(jīng)被人們普遍認(rèn)同的慢速凝聚和快速凝聚之間的聯(lián)系,這樣才能實(shí)現(xiàn)模型不但準(zhǔn)確性,普適性好,而且模型的物理意義也明確.

    多年來,人們?yōu)榇_定這兩個(gè)模型作了艱苦的努力,6-8直到分形理論5,9,13,22的誕生及其在膠體科學(xué)中的應(yīng)用,才使問題有了轉(zhuǎn)機(jī).現(xiàn)在許多實(shí)驗(yàn)表明,23-28膠體顆粒的凝聚體是分形結(jié)構(gòu),所以,傳統(tǒng)的模型應(yīng)該就分形特性的加入進(jìn)行調(diào)整.Ziff等1證明,如果考慮到團(tuán)簇的分形特性,Smoluchowski方程可以比較準(zhǔn)確地表達(dá)膠體的凝聚動(dòng)力學(xué).根據(jù)Kirkwood-Riseman理論,7,8,12對(duì)于不規(guī)則形狀的顆?;驁F(tuán)簇,其移動(dòng)的布朗擴(kuò)散系數(shù)Ds與其質(zhì)量s(實(shí)際可以用團(tuán)簇中的粒子數(shù)表達(dá))的關(guān)系可以用下式描述:

    其中D0是單個(gè)粒子的擴(kuò)散系數(shù),γ則稱為擴(kuò)散指數(shù)(diffusivity exponent).從公式(1)可以看出,對(duì)于γ=0的情況,實(shí)際就是大小團(tuán)簇的擴(kuò)散系數(shù)相同,與團(tuán)簇的質(zhì)量無(wú)關(guān);對(duì)于γ<0的情況,則對(duì)應(yīng)小團(tuán)簇比大團(tuán)簇移動(dòng)快;對(duì)于γ>0的情況,則對(duì)應(yīng)小團(tuán)簇比大團(tuán)簇移動(dòng)慢.已有文獻(xiàn)5,6,8,15表明,當(dāng)γ足夠大時(shí),在CCA模型下會(huì)得到與DLA模型類似的凝聚體.DLA模型所對(duì)應(yīng)的凝聚體是在有“凝聚中心”的前提下產(chǎn)生的,5,6,8,15但為什么會(huì)產(chǎn)生“凝聚中心”?為什么會(huì)出現(xiàn)小團(tuán)簇比大團(tuán)簇移動(dòng)慢的現(xiàn)象?也就是說γ>0的物理意義何在?在此以前還未見文獻(xiàn)報(bào)道.對(duì)于這些問題,本研究試圖給一個(gè)合理的解釋.

    根據(jù)Meakin的研究,8,15碰撞在一起的兩個(gè)顆?;驁F(tuán)簇粘結(jié)在一起的概率Pij可用下式描述:

    其中i、j分別代表碰撞在一起的兩個(gè)顆?;驁F(tuán)簇的質(zhì)量,σ稱為凝聚概率指數(shù)(sticking probability exponent),P0表示兩個(gè)初級(jí)粒子由于布朗運(yùn)動(dòng)碰撞而粘結(jié)在一起的概率,顯然P0≤1,當(dāng)Pij>1時(shí)則被設(shè)置為1.從公式(2)可以看出,對(duì)于σ=0的情況,實(shí)際就是大小團(tuán)簇碰撞的有效碰撞概率相同,與團(tuán)簇的質(zhì)量無(wú)關(guān);對(duì)于σ>0的情況,團(tuán)簇間的有效碰撞概率隨團(tuán)簇質(zhì)量的增加而增大;對(duì)于σ<0的情況,團(tuán)簇間的有效碰撞概率隨團(tuán)簇質(zhì)量的增加而減小.以前的研究2,8,15大多集中在σ≥0的情況,關(guān)于σ<0情況是否有明確的物理意義的探討,至今也還未見報(bào)道.

    Meakin等6分別利用公式(1)、(2),從團(tuán)簇的移動(dòng)性和團(tuán)簇的凝聚概率的角度,單獨(dú)考察了膠粒凝聚過程中團(tuán)簇大小的分布,發(fā)現(xiàn)團(tuán)簇大小的分布具有標(biāo)度性質(zhì).但Ziff等1指出僅從團(tuán)簇的大小分布來證明模型的正確性是不夠的,由此來研究膠體的凝聚動(dòng)力學(xué)行為也是不全面的.因此必須比較系統(tǒng)地考慮團(tuán)簇的移動(dòng)性和團(tuán)簇碰撞的凝聚概率對(duì)膠體凝聚動(dòng)力學(xué)的影響.

    先期的研究已表明雖然各種膠體的特性不一樣,但其凝聚的機(jī)理與動(dòng)力學(xué)卻有相似之處.20,21,29,30因此本研究將凝聚過程動(dòng)力學(xué)的各個(gè)影響因素的作用歸結(jié)為擴(kuò)散指數(shù)和凝聚概率指數(shù)的作用,通過調(diào)整擴(kuò)散指數(shù)和凝聚概率指數(shù)大小來探討膠體凝聚的動(dòng)力學(xué)過程.

    2 模擬研究方法

    模擬在邊長(zhǎng)為100個(gè)單位膠粒直徑長(zhǎng)度的立方體中進(jìn)行,立方體被劃分成邊長(zhǎng)為單位膠粒直徑長(zhǎng)度的小立方體格子(總共1003個(gè)),初始時(shí)N0=10000個(gè)粒子被無(wú)重疊地隨機(jī)置放在三維網(wǎng)格中,粒子在X、Y、Z軸的6個(gè)方向上做隨機(jī)運(yùn)動(dòng).模擬中膠?;驁F(tuán)簇移動(dòng)、碰撞、粘結(jié)的規(guī)則與文獻(xiàn)6,17相似,模擬時(shí)間按式(3)增長(zhǎng):12

    其中,Dmax是當(dāng)前體系中團(tuán)簇(含單粒)的最大擴(kuò)散系數(shù),是在模擬過程中進(jìn)行統(tǒng)計(jì)比較得出的,Nc是當(dāng)前體系中的團(tuán)簇?cái)?shù)(包括單粒),由此得到模擬實(shí)驗(yàn)的時(shí)間的標(biāo)度是在相同環(huán)境條件下單個(gè)膠體粒子移動(dòng)單位膠粒直徑所需的時(shí)間.模擬中考察不同擴(kuò)散指數(shù)和凝聚概率指數(shù)及其組合對(duì)凝聚過程的影響.

    對(duì)公式(1)、(2)而言,以前的研究結(jié)果8,12大多只是考慮了γ≤0和σ>0時(shí)的情況,因?yàn)檫@兩種情況與實(shí)際的膠體體系在常規(guī)情況下的凝聚特征是基本相符的.本模擬實(shí)驗(yàn)為了比較全面地研究這兩個(gè)因子的影響,以及有效地表征前人提出的三個(gè)凝聚模型(DLA、RLCA和DLCA)13,31之間的聯(lián)系,故意設(shè)定了擴(kuò)散指數(shù)和凝聚概率指數(shù)較寬的取值范圍.根據(jù)Einstein-Stokes定律,γ>0似乎是很不合理的.但后面的對(duì)模擬結(jié)果討論將表明,γ>0仍然可能有確定的物理意義.擴(kuò)散指數(shù)γ的取值為:2.0,1.0,0.5,0.1,0,-0.5至-3.0(步長(zhǎng)為-0.1);凝聚概率指數(shù)σ的取值為:3.0,2.0,1.0,0.5,0.1,0,-0.5至-2.0(步長(zhǎng)為-0.1).

    在模擬研究中,隨著時(shí)間t的延伸,團(tuán)簇中的粒子數(shù)不斷增多,團(tuán)簇的數(shù)目不斷減少,我們可以對(duì)含粒子數(shù)s的團(tuán)簇?cái)?shù)Ns(t)進(jìn)行統(tǒng)計(jì).在三維情形下,設(shè)邊長(zhǎng)為L(zhǎng)的立方體點(diǎn)陣中的總粒子數(shù)為N0,則粒子密度ρ為N0/L3,對(duì)團(tuán)簇?cái)?shù)Ns(t)除以L3進(jìn)行約化以消除點(diǎn)陣尺寸的影響,得到:

    這里的ns(t)是t時(shí)刻單位體積內(nèi)含s個(gè)粒子的團(tuán)簇的數(shù)目.研究表明,7,23-25團(tuán)簇大小分布可用式(5)描述:

    式中,t為聚集時(shí)間,τ的取值一般介于-2和-2.5之間,文獻(xiàn)7已對(duì)其取值作了比較詳細(xì)的研究,根據(jù)文獻(xiàn)7的研究結(jié)果,本文約定τ的取值為-2,S(t)為團(tuán)簇大小的質(zhì)量權(quán)重平均:

    公式(5)中的f(x)為動(dòng)態(tài)標(biāo)度函數(shù).已有的實(shí)驗(yàn)結(jié)果7,23-25表明,凝聚過程的早期S(t)隨時(shí)間t按指數(shù)規(guī)律增長(zhǎng),即有S(t)~eat,這里a為指數(shù)增長(zhǎng)因子;而在凝聚過程的后期S(t)隨時(shí)間t按冪函數(shù)規(guī)律增長(zhǎng),即S(t)~tb.與前人研究16相同,本模擬實(shí)驗(yàn)中對(duì)團(tuán)簇分布的標(biāo)度關(guān)系通過Ns(t)S(t)2對(duì)s/S(t)作雙對(duì)數(shù)圖來驗(yàn)證,而對(duì)團(tuán)簇大小的質(zhì)量權(quán)重平均隨時(shí)間的增長(zhǎng)關(guān)系可以通過對(duì)S(t)取對(duì)數(shù)對(duì)t作圖和通過S(t)對(duì)t作雙對(duì)數(shù)圖來觀察.

    3 結(jié)果與討論

    3.1 膠體體系隨時(shí)間的演變

    圖1表示了不同γ和σ及其組合作用下膠體體系在凝聚過程中體系團(tuán)簇總濃度及含有不同粒子數(shù)的團(tuán)簇的濃度隨時(shí)間的變化關(guān)系.圖中的藍(lán)虛線代表體系團(tuán)簇總濃度的變化趨勢(shì),另外的5條曲線分別代表了單粒子、2粒子團(tuán)簇、3粒子團(tuán)簇、4粒子團(tuán)簇和5粒子團(tuán)簇濃度的變化趨勢(shì),其它大小團(tuán)簇濃度的變化趨勢(shì)與5粒子團(tuán)簇類似,其曲線沒有畫出.團(tuán)簇濃度用體積分?jǐn)?shù)表達(dá),即由體系中這種大小的團(tuán)簇的總數(shù)除以立方體的格點(diǎn)數(shù)而得到.體系聚集時(shí)間的標(biāo)度是單粒子在體系中移動(dòng)單位粒子直徑距離所需要的時(shí)間.

    可以看出,膠體的凝聚過程是團(tuán)簇總濃度逐漸減小的過程,對(duì)于單粒子來說,其濃度都是單調(diào)下降,其它團(tuán)簇的濃度則都是先上升,而后下降,且團(tuán)簇濃度曲線上都有一最高點(diǎn).但對(duì)于不同的體系(即不同的擴(kuò)散指數(shù)和凝聚概率指數(shù)取值),膠體團(tuán)簇總濃度及單粒子濃度變化的快慢是不同的,其它團(tuán)簇濃度曲線的峰高及峰寬也是有較大差異.

    圖1 體系的團(tuán)簇濃度隨時(shí)間的變化關(guān)系Fig.1 Time evolution of concentrations of colloidal clusters in the system

    模擬結(jié)果表明,在σ取定值的前提下,當(dāng)γ<0時(shí),隨著γ的增大,峰的高度(團(tuán)簇的最大濃度)逐漸變低,而峰寬逐漸增大.原因可能是在γ較小的情形中,團(tuán)簇的擴(kuò)散系數(shù)受團(tuán)簇大小的影響大,團(tuán)簇越大,擴(kuò)散系數(shù)越小.在γ=0時(shí),團(tuán)簇的擴(kuò)散系數(shù)為一常數(shù),這就意味著,大小凝聚體移動(dòng)的速率相同,濃度變化的速度也就相同.因此,在γ=0時(shí),單粒子濃度下降線不與其它團(tuán)簇的變化曲線相交,實(shí)際上,其它各種團(tuán)簇濃度變化曲線也是不相交的,即理論上在某一時(shí)刻點(diǎn)始終存在各種團(tuán)簇共存的現(xiàn)象(圖1(a)).在γ<0時(shí),小團(tuán)簇總會(huì)逐漸消失殆盡,大團(tuán)簇的濃度逐漸增大并超過小團(tuán)簇的濃度(圖1(b)).在γ>0時(shí),表現(xiàn)出大團(tuán)簇比小團(tuán)簇運(yùn)動(dòng)快,當(dāng)γ很大時(shí),則類似于有單一凝聚核的DLA模型(圖1(e)).9

    圖2反映了γ對(duì)團(tuán)簇均勻性的影響.從圖2可以看出,γ很小時(shí),形成的團(tuán)簇非常均勻(圖2(a)).當(dāng)γ很大時(shí),則體系中基本上只有最大的那個(gè)團(tuán)簇在移動(dòng),而且這個(gè)團(tuán)簇非常類似于DLA模型形成的產(chǎn)物(圖2(b)).該結(jié)果表明,γ>0描述了體系有確定凝聚中心(可能是一個(gè)或多個(gè))時(shí)的膠體凝聚情形,而當(dāng)γ>>0時(shí),體系只有一個(gè)凝聚中心.那么在實(shí)際問題中,這種凝聚中心是如何產(chǎn)生的呢?毫無(wú)疑問,凝聚中心存在的前提是:大團(tuán)簇是一個(gè)“吸引中心”,即該“吸引中心”對(duì)體系中其它小團(tuán)簇或單粒的吸引力大于小團(tuán)簇與小團(tuán)簇、小團(tuán)簇與單粒及單粒與單粒間的吸引力.在這種情形下,膠體顆?;蛐F(tuán)簇的移動(dòng)不受布朗運(yùn)動(dòng)的制約,而可能是受一種由大團(tuán)簇產(chǎn)生的作用力場(chǎng)(如強(qiáng)長(zhǎng)程范德華力,32電場(chǎng)力等)或體系邊界處的外力場(chǎng)的作用.由此可以看出:γ<0適合于描述顆粒的碰撞概率由布朗運(yùn)動(dòng)決定的情況;而γ>0適合于描述顆粒的碰撞概率由顆粒的“定向移動(dòng)”決定的情況,而該“定向移動(dòng)”的推動(dòng)力可來自于大團(tuán)簇產(chǎn)生的強(qiáng)“長(zhǎng)程范德華力”,“電場(chǎng)力”等,或來自于體系邊界處的外力場(chǎng)作用(如重力場(chǎng)33中的沉降過程).

    在γ取定值的前提下,σ對(duì)膠體體系的凝聚過程的影響類似于γ對(duì)膠體凝聚過程的影響.團(tuán)簇總濃度隨著時(shí)間的延長(zhǎng)而逐漸降低是必然的趨勢(shì),但當(dāng)σ<<0時(shí),后期的凝聚幾乎停滯,因此計(jì)算機(jī)模擬時(shí)也很難達(dá)到最后形成一個(gè)團(tuán)簇的結(jié)果(圖1(c));當(dāng)σ>>0時(shí),則逐漸演變成為快速凝聚,直接受初級(jí)粒子有效碰撞概率的影響,即受公式(2)中P0的影響,凝聚過程與圖1(a)類似.σ從-2.0變化到3.0的總體趨勢(shì)是體系向終態(tài)(形成一個(gè)團(tuán)簇)發(fā)展所需的時(shí)間逐漸縮短.但σ增加到一定程度后,體系凝聚成一個(gè)團(tuán)簇的總時(shí)間就不會(huì)再縮短了,即保持快速凝聚所需的時(shí)間.σ具體增加到什么程度而使慢速凝聚變?yōu)榭焖倌?依初始粒子的有效碰撞概率P0而定.通過公式(2)可以推算出,當(dāng)滿足σ≥-ln P0/ln(ij)時(shí),大小為i、j的兩個(gè)團(tuán)簇碰撞必然發(fā)生凝聚;當(dāng)σ為一固定的正值,i、j滿足 ij≥時(shí),體系則由慢速凝聚過渡到快速凝聚.

    圖2 擴(kuò)散指數(shù)對(duì)團(tuán)簇大小均勻性的影響Fig.2 Effect of diffusivity exponent on clusters′uniformity of size grading

    在各種情形中單粒子的濃度都是單調(diào)減少的,在σ較小的情況下,大團(tuán)簇的濃度逐漸超過小團(tuán)簇的濃度(圖1(c)),在σ較大的情況下,小團(tuán)簇的濃度始終大于大團(tuán)簇的濃度,與γ的影響相同,最終這兩種情況下體系中團(tuán)簇的濃度都趨向一樣,與圖1(a)表達(dá)的情況類似.與γ不同的是在σ接近-2.0時(shí),峰的高度逐漸變高,峰寬也逐漸增大,在σ趨向0的過程中峰高逐漸變低,峰寬也逐漸減小,到σ>0后,峰高,峰寬則維持基本不變,這說明在γ=0的情況下,σ在這個(gè)范圍內(nèi)對(duì)體系團(tuán)簇的大小分布影響不太明顯,因?yàn)榈侥鄣暮笃趯?shí)際上已基本過渡到了快速凝聚階段,σ的大小只是影響過渡的早晚,σ越大過渡得越快,σ越小過渡得越晚.σ<0和γ<0時(shí),隨著凝聚過程的發(fā)生,體系團(tuán)簇都體現(xiàn)出均勻分布的特性.σ<0所對(duì)應(yīng)的均勻性來自于團(tuán)簇的有效碰撞概率隨團(tuán)簇質(zhì)量的增大而減小,碰撞主要發(fā)生在兩個(gè)小團(tuán)簇之間,因此大團(tuán)簇的繼續(xù)擴(kuò)張得到遏制,凝聚過程到后期也就很難繼續(xù)進(jìn)行;γ<0所對(duì)應(yīng)的均勻性來自于單?;蛐F(tuán)簇的移動(dòng)幾率大,從而發(fā)生碰撞的可能性也就大,而且這種碰撞主要是因?yàn)樾F(tuán)簇的移動(dòng)造成的,但這并不妨礙大團(tuán)簇的繼續(xù)擴(kuò)張,凝聚過程能繼續(xù)進(jìn)行.顯然,γ<0與實(shí)際膠體體系在常規(guī)情況下的凝聚非常相似,但對(duì)σ<0能否給出合適的物理解釋,還期待進(jìn)一步的研究.

    通過大量的γ和σ的組合模擬實(shí)驗(yàn),發(fā)現(xiàn)γ對(duì)凝聚過程中體系團(tuán)簇均勻程度的影響占主導(dǎo)作用(圖1(d,e,f)),而σ則主要影響團(tuán)簇的密實(shí)程度,這在后邊的研究中可以得到證實(shí).根據(jù)已有的研究,7,12圖1(d)能較好地反映膠體體系在常規(guī)情況下的凝聚過程,即體現(xiàn)了小團(tuán)簇比大團(tuán)簇移動(dòng)快,有效碰撞概率隨團(tuán)簇的增大而增大的特性.

    3.2 團(tuán)簇大小分布及其標(biāo)度關(guān)系

    圖3為體系團(tuán)簇的大小分布.每個(gè)子圖中的10條曲線分別代表不同時(shí)刻(體系團(tuán)簇?cái)?shù)減少至分別為 0.95N0、0.90N0、0.80N0、0.70N0、0.60N0、0.50N0、0.40 N0、0.30N0、0.20N0、0.10 N0),粒子數(shù)為s的團(tuán)簇個(gè)數(shù)Ns(t)對(duì)s的雙對(duì)數(shù)關(guān)系,其中N0為初始模擬體系的粒子數(shù).這也是Meakin等6從團(tuán)簇的移動(dòng)性和Family等7從團(tuán)簇凝聚概率的角度,單獨(dú)考察膠粒凝聚過程中團(tuán)簇大小分布的主要工作.他們認(rèn)為膠體體系團(tuán)簇大小分布存在轉(zhuǎn)折點(diǎn).在主要考慮γ的影響時(shí)(令σ=0),本模擬實(shí)驗(yàn)在γ<-0.7時(shí),團(tuán)簇大小呈現(xiàn)鐘罩形分布,而當(dāng)γ>-0.7,團(tuán)簇分布呈現(xiàn)單調(diào)下降;在主要考慮σ的影響時(shí)(令γ=0),當(dāng)σ<-0.9時(shí),團(tuán)簇大小呈現(xiàn)鐘罩形分布,而當(dāng)σ>-0.9,團(tuán)簇分布呈現(xiàn)單調(diào)下降.圖3(a,b)顯示了這個(gè)轉(zhuǎn)折點(diǎn)時(shí)團(tuán)簇的大小分布,這個(gè)轉(zhuǎn)折點(diǎn)的意義在于當(dāng)γ>-0.7時(shí)或σ>-0.9時(shí),只要s1<s2,Ns1(t)恒大于Ns2(t),即所含粒子數(shù)少的團(tuán)簇的個(gè)數(shù)始終大于所含粒子數(shù)多的團(tuán)簇的個(gè)數(shù),直觀表現(xiàn)是團(tuán)簇分布高度分散.

    本模擬實(shí)驗(yàn)還探討了γ和σ的組合影響,以及團(tuán)簇大小分布曲線的包跡(即各個(gè)曲線的公切線)的特征.

    模擬實(shí)驗(yàn)發(fā)現(xiàn)γ和σ對(duì)團(tuán)簇大小分布的組合影響體現(xiàn)一個(gè)正負(fù)效應(yīng)的關(guān)系,即γ和σ的取值如果都有使團(tuán)簇大小分布呈現(xiàn)鐘罩形分布的趨勢(shì),則組合后鐘罩形分布的趨勢(shì)更明顯;γ和σ的取值如果都有使團(tuán)簇大小分布呈單調(diào)下降的趨勢(shì),則組合后單調(diào)下降的趨勢(shì)更陡,如圖3(e)的曲線比圖3(f)的曲線更陡;如果γ和σ的取值效應(yīng)相反,則組合后彼此削弱,對(duì)比圖3(e)和3(f),σ=0.5大于轉(zhuǎn)折點(diǎn)(-0.9),其效應(yīng)是使團(tuán)簇分布呈現(xiàn)單調(diào)下降,但在圖3(f)中γ=-2,其效應(yīng)是使曲線變得彎曲,由于其值遠(yuǎn)小于轉(zhuǎn)折點(diǎn)(-0.7),所以曲線還是呈現(xiàn)了鐘罩形.

    圖3最重要的特征是各時(shí)間點(diǎn)團(tuán)簇大小分布曲線的包跡是一條直線,這說明團(tuán)簇大小分布有動(dòng)態(tài)標(biāo)度性質(zhì),7,8,12γ和σ組合作用時(shí)會(huì)對(duì)此斜率有一些改變(圖3(c,e)),但團(tuán)簇大小分布曲線的包跡仍然是一條直線,這說明團(tuán)簇大小分布的標(biāo)度性質(zhì)是依然存在的.另外一種考察團(tuán)簇大小分布是否具有標(biāo)度性質(zhì)的方法是通過Ns(t)S(t)2對(duì)s/S(t)作雙對(duì)數(shù)圖,發(fā)現(xiàn)不同時(shí)刻體系的團(tuán)簇分布基本落在同一條標(biāo)度曲線上(圖4),這說明團(tuán)簇的大小分布可被動(dòng)態(tài)標(biāo)度,7,25即各階段膠體的凝聚過程具有同一性.

    3.3 團(tuán)簇重均大小S(t)的變化

    在不同的P0、γ和σ作用下,一般都體現(xiàn)出在凝聚的初始階段S(t)隨t緩慢增加而后迅速增加的規(guī)律.這說明膠體的凝聚動(dòng)力學(xué)過程是一個(gè)存在正反饋機(jī)制的非線性動(dòng)力學(xué)過程.凝聚開始時(shí),體系形成了一些小的團(tuán)簇,但因小的團(tuán)簇的范德華力遠(yuǎn)遠(yuǎn)大于單粒,所以小的團(tuán)簇將加速凝聚速度.而凝聚速度的增加使這些小的團(tuán)簇以更快的速度凝聚成更大的團(tuán)簇,這反過來又使凝聚速度進(jìn)一步加快而形成更大的團(tuán)簇,這種過程可使體系在瞬間迅速完成聚合過程,圖5中團(tuán)簇重均大小S(t)隨時(shí)間的變化曲線可以清楚地反映出這種正反饋現(xiàn)象.

    P0、γ和σ對(duì)這種正反饋機(jī)制有不同的影響,在γ<0或σ>0,當(dāng)P0取定值時(shí),σ和γ對(duì)S(t)的影響并不明顯,即對(duì)S(t)隨t由緩慢增加何時(shí)向迅速增加過渡的影響似乎并不明顯;相反在γ>0或σ<0,當(dāng)P0取定值時(shí),γ的增大將推動(dòng)S(t)隨t由緩慢增加向迅速增加的過渡,而且這種過渡更迅速,而σ的減小將延遲S(t)隨t由緩慢增加向迅速增加的過渡,而且這種過渡也變得緩慢.當(dāng)σ<<0時(shí),凝聚動(dòng)力學(xué)過程就不再是正反饋機(jī)制的動(dòng)力學(xué)過程,相反,而變成負(fù)反饋機(jī)制的動(dòng)力學(xué)過程,由于在這種情況下團(tuán)簇之間的有效碰撞概率隨著團(tuán)簇質(zhì)量的增大而減小,所以后期的凝聚過程將變得越來越慢.σ<<0可能對(duì)應(yīng)了膠體體系中單?;驁F(tuán)簇之間存在某種排斥力場(chǎng),而且這種排斥力場(chǎng)隨著團(tuán)簇質(zhì)量的增大而逐漸增強(qiáng)的現(xiàn)象.

    圖3 體系團(tuán)簇的大小分布Fig.3 Cluster-size distribution in the system

    在γ和σ取定值時(shí),P0對(duì)S(t)隨t由緩慢增加何時(shí)向迅速增加過渡的影響非常明顯,圖5表示了γ=-0.5、σ=0.25,P0對(duì)這種過渡的影響.從圖5中可以看出,隨著P0的減小,這種過渡將大大延遲,而且這種過渡也逐漸變得緩慢.從模擬實(shí)驗(yàn)看出,P0和σ對(duì)S(t)隨t變化的影響最為明顯.

    圖4 團(tuán)簇大小分布時(shí)間標(biāo)度關(guān)系Fig.4 Scaling relations of the time-dependent cluster-size distribution

    較多的文獻(xiàn)討論了團(tuán)簇重均大小的標(biāo)度關(guān)系,7,8,12,18,19普遍認(rèn)為團(tuán)簇重均大小與體系團(tuán)簇大小分布的聯(lián)系可用公式(5)、(6)表達(dá).本模擬實(shí)驗(yàn)主要研究了快速凝聚和慢速凝聚情況下團(tuán)簇重均大小增長(zhǎng)規(guī)律的聯(lián)系和區(qū)別.圖6(a,b)主要表達(dá)了在初級(jí)粒子有效碰撞概率P0=1的情況下,γ和σ單獨(dú)作用對(duì)團(tuán)簇重均大小增長(zhǎng)規(guī)律的影響,其中最右綠色的實(shí)線則代表了在P0=0.01、γ=-0.5、σ=0.25時(shí)團(tuán)簇重均大小隨時(shí)間增長(zhǎng)的雙對(duì)數(shù)關(guān)系.

    圖5 團(tuán)簇重均大小S(t)隨時(shí)間的變化Fig.5 Time dependence of the weight-averaged cluster size S(t)

    從圖6(a)可以看出,在σ=0時(shí),當(dāng)γ>0時(shí),隨著γ的減少,團(tuán)簇的平均增長(zhǎng)冪函數(shù)指數(shù)也相應(yīng)減小,表現(xiàn)為團(tuán)簇的結(jié)構(gòu)變得疏松,實(shí)際上γ>>0時(shí),則對(duì)應(yīng)DLA模型.當(dāng)γ<0時(shí),各條S(t)隨t的增長(zhǎng)曲線重合在一起,這說明在擴(kuò)散控制的凝聚中,粒子的移動(dòng)與其大小和形狀有關(guān),兩個(gè)粒子一旦相遇即發(fā)生凝聚,且一個(gè)重要特征是不同體系的凝聚具有同一性,其變化規(guī)律是一致的.

    從圖6(b)可以看出,在γ=0時(shí),當(dāng)σ<0時(shí),隨著σ的減少,團(tuán)簇的平均增長(zhǎng)冪函數(shù)指數(shù)也相應(yīng)減小,表現(xiàn)為團(tuán)簇的結(jié)構(gòu)變得疏松,實(shí)際上σ<<0時(shí),則對(duì)應(yīng)慢速凝聚模型.當(dāng)σ>0時(shí),團(tuán)簇間的粘結(jié)與其大小和形狀有關(guān),團(tuán)簇間的有效碰撞概率隨σ的增大而增大,但到一定程度,實(shí)際已發(fā)展成為快速凝聚,即有效碰撞概率為1,圖6中σ>0時(shí)S(t)的各條增長(zhǎng)曲線幾乎完全重合在一起的現(xiàn)象正說明了這一點(diǎn).

    關(guān)于團(tuán)簇重均大小隨時(shí)間增長(zhǎng)規(guī)律的問題,理論和實(shí)際之間還有較大的差距.大多數(shù)計(jì)算機(jī)模擬的結(jié)果顯示S(t)隨t按冪函數(shù)規(guī)律增長(zhǎng),2,6但實(shí)驗(yàn)研究則出現(xiàn)S(t)先隨t按指數(shù)規(guī)律增長(zhǎng),而后按冪函數(shù)規(guī)律增長(zhǎng)的現(xiàn)象.23-25本模擬實(shí)驗(yàn)對(duì)P0=0.01、γ=-0.5、σ=0.25的情況,早期按對(duì)S(t)取對(duì)數(shù)后對(duì)t作圖,后期則作S(t)與t的雙對(duì)數(shù)圖,觀察發(fā)現(xiàn)圖6(c,d)確實(shí)都呈現(xiàn)直線關(guān)系,這在一定程度上證實(shí)了實(shí)驗(yàn)研究的結(jié)果,同時(shí)也說明本實(shí)驗(yàn)使用擴(kuò)散模型和凝聚模型相結(jié)合的方法來模擬膠體體系的凝聚動(dòng)力學(xué)行為是正確的.

    圖6中P0=0.01、γ=-0.5、σ=0.25的曲線表明,慢速凝聚在初期與快速凝聚有較大的差別,但普通的慢速凝聚進(jìn)行到一定程度后會(huì)發(fā)展為快速凝聚,即到后期團(tuán)簇重均大小的增長(zhǎng)規(guī)律彼此是一樣的.值得注意的是P0=0.01所代表的慢速凝聚與圖6(b)中σ<<0時(shí)所代表的慢速凝聚有本質(zhì)的區(qū)別,這個(gè)問題在3.4節(jié)有詳細(xì)的論述.

    3.4 擴(kuò)散指數(shù)和凝聚概率指數(shù)對(duì)團(tuán)簇分形維數(shù)的影響

    采用盒維數(shù)方法25計(jì)算體系團(tuán)簇的分形維數(shù),如表1所示.每個(gè)分形維數(shù)對(duì)應(yīng)3次模擬實(shí)驗(yàn),對(duì)結(jié)果求算術(shù)平均,表中除σ=-1.0,-1.5,-2.0外,其它都是最終凝聚形成了一個(gè)凝聚體.

    在初始條件相同的情況下,當(dāng)γ<0時(shí),凝聚體的分形維數(shù)基本上是相同的,都為1.93左右,與前面團(tuán)簇的平均大小S(t)的變化規(guī)律是一致的,這也進(jìn)一步說明不同體系快速凝聚的過程具有同一性.然而,當(dāng)γ>>0,在粒子濃度很低時(shí),團(tuán)簇凝聚變成與DLA相當(dāng),因?yàn)樵谶@種情況下實(shí)際上只有一個(gè)單獨(dú)的團(tuán)簇(即最大的團(tuán)簇)在運(yùn)動(dòng),并且在其擴(kuò)散運(yùn)動(dòng)時(shí),它將余下的單個(gè)粒子都收集于自身,表現(xiàn)為分形維數(shù)明顯增加.

    圖6 團(tuán)簇重均大小的標(biāo)度關(guān)系Fig.6 Scaling relations of the weight-averaged cluster size

    表1 擴(kuò)散指數(shù)和凝聚概率指數(shù)對(duì)分形維數(shù)的影響Table 1 Effects of diffusivity exponent and stickingprobability exponent on fractal dimension

    同樣,當(dāng)σ>0時(shí),凝聚體的分形維數(shù)為1.93左右,與σ=0、γ<0時(shí)的情況非常類似,這時(shí)膠體體系過渡到后期實(shí)際已完全變成擴(kuò)散控制凝聚,因?yàn)榕鲎苍谝黄鸬膬蓚€(gè)團(tuán)簇粘結(jié)在一起的概率隨著團(tuán)簇的質(zhì)量增大而增大,到一定程度有效碰撞概率Pij(σ)實(shí)際已等于1,即碰撞就導(dǎo)致粘結(jié),這說明反應(yīng)控制凝聚(RLCA)會(huì)最終過渡到擴(kuò)散控制凝聚(DLCA),這也是最終凝聚體分形維數(shù)相差不大的根本原因.

    σ對(duì)團(tuán)簇分形維數(shù)的影響最有趣的結(jié)果在于σ<0的情況.如前文所述,σ<0對(duì)應(yīng)了慢速凝聚過程,大多文獻(xiàn)8,12,13都認(rèn)為慢速凝聚得到的團(tuán)簇的分形維數(shù)要比快速凝聚的大,但為什么本研究結(jié)果得到的分形維數(shù)有隨著σ的減小而逐漸減小的趨勢(shì)呢?我們認(rèn)為其原因在于σ<0時(shí),碰撞在一起的兩個(gè)團(tuán)簇粘結(jié)在一起的概率隨著團(tuán)簇的質(zhì)量增大而減小,在凝聚過程中則體現(xiàn)為兩個(gè)級(jí)別小的團(tuán)簇更容易粘結(jié)在一起,也就是說大團(tuán)簇的產(chǎn)生是兩個(gè)小團(tuán)簇聚結(jié)而成的,而不是一個(gè)大團(tuán)簇與一個(gè)小團(tuán)簇聚結(jié)而成的(γ>>0的情況大多如此),在這種情況下形成的團(tuán)簇必然有較小的分形維數(shù).大多文獻(xiàn)中所說的慢速凝聚是因?yàn)樵诒3至税l(fā)生碰撞的兩個(gè)團(tuán)簇粘結(jié)在一起的概率隨著團(tuán)簇的質(zhì)量增大而增大的規(guī)律的前提下,初級(jí)粒子之間的有效碰撞概率P0就很低,在這種情況下,大團(tuán)簇大多是大團(tuán)簇與大團(tuán)簇或大團(tuán)簇與小團(tuán)簇聚結(jié)而成的,其分形維數(shù)必然要比σ<0時(shí)形成的團(tuán)簇的分形維數(shù)大.當(dāng)σ<0時(shí),體系成為先快后慢的慢速凝聚,即顆粒越大,有效碰撞概率越小,這可能對(duì)應(yīng)大團(tuán)簇為一排斥中心,即膠體顆粒存在“排斥力場(chǎng)”的現(xiàn)象.但這只是從凝聚模型得到的一個(gè)有趣現(xiàn)象,在什么情況下會(huì)存在“排斥中心”?是否有明確的物理意義?還需作進(jìn)一步的探索.

    總之,本研究結(jié)果表明,在初級(jí)粒子之間的有效碰撞概率就很低的情況下,所謂慢速凝聚從理論上可以分為先慢后快和先快后慢兩種情況,先慢后快對(duì)應(yīng)了σ>0的情況,先快后慢對(duì)應(yīng)了σ<0的情況.注意,這里提出的先快后慢或先慢后快是對(duì)各自的整個(gè)凝聚過程而言的,當(dāng)P0<<1時(shí),σ>0對(duì)應(yīng)的慢速凝聚過程可以用“越來越快”描述其凝聚過程,而σ<0對(duì)應(yīng)的慢速凝聚過程可以用“越來越慢”表達(dá).

    4 結(jié)論

    Meakin和Family各自獨(dú)立地考慮團(tuán)簇的移動(dòng)性和團(tuán)簇碰撞的凝聚概率,對(duì)膠體體系凝聚過程中團(tuán)簇大小分布進(jìn)行了研究,本模擬實(shí)驗(yàn)在此基礎(chǔ)上,以公式(1)為團(tuán)簇的移動(dòng)模型,以公式(2)為團(tuán)簇的凝聚模型,并將凝聚過程動(dòng)力學(xué)的各個(gè)影響因素的作用歸結(jié)為擴(kuò)散指數(shù)γ和凝聚概率指數(shù)σ的作用,系統(tǒng)地考慮團(tuán)簇的移動(dòng)性和團(tuán)簇碰撞的凝聚概率對(duì)膠體凝聚動(dòng)力學(xué)行為的影響,得到了如下幾個(gè)重要結(jié)論.

    (1)模擬實(shí)驗(yàn)證實(shí)了在凝聚過程中,早期的團(tuán)簇重均大小按指數(shù)規(guī)律增長(zhǎng),而后期按冪函數(shù)規(guī)律增長(zhǎng)的實(shí)驗(yàn)現(xiàn)象,這也進(jìn)一步說明了膠體的凝結(jié)動(dòng)力學(xué)過程在常規(guī)情況下,是一個(gè)存在正反饋機(jī)制的非線性動(dòng)力學(xué)過程,而在σ<<0時(shí),模擬研究顯示該凝聚動(dòng)力學(xué)過程具有負(fù)反饋特性,不過該情況是否具有明確的物理意義,目前尚不清楚.

    (2)以CCA模型為紐帶,考慮γ和σ的作用,將DLCA、RLCA與經(jīng)典的DLA模型有機(jī)地聯(lián)系在一起:

    根據(jù)γ,可以將快速凝聚分為DLCA式的快速凝聚和DLA式的快速凝聚.前者對(duì)應(yīng)了無(wú)固定凝結(jié)核的擴(kuò)散控制凝聚過程,后者對(duì)應(yīng)了有固定凝結(jié)核的擴(kuò)散控制凝聚過程.已有的研究表明,γ>0時(shí),前者可逐漸過渡到后者.本文的分析進(jìn)一步表明,γ>>0實(shí)際代表了體系中團(tuán)簇或單粒做“定向運(yùn)動(dòng)”而非無(wú)規(guī)則的布朗運(yùn)動(dòng)的情況.這種“定向運(yùn)動(dòng)”的推動(dòng)力可能來自于大團(tuán)簇產(chǎn)生的強(qiáng)“長(zhǎng)程范德華力”,“電場(chǎng)力”等,或來自于體系邊界處的外力場(chǎng)的作用.

    根據(jù)σ,可以將慢速凝聚分為先慢后快,即“越來越快”的慢速凝聚和先快后慢,即“越來越慢”的慢速凝聚.前者對(duì)應(yīng)了P0<<1、σ>0的情況,即常規(guī)的膠體凝聚,后者對(duì)應(yīng)了P0<<1、σ<<0的情況,即可能大團(tuán)簇為一排斥中心,膠體顆粒存在“排斥力場(chǎng)”的現(xiàn)象;前者為研究慢速凝聚向快速凝聚的過渡提供了依據(jù),后者對(duì)膠體的穩(wěn)定性研究可能有一定的參考價(jià)值.

    (3)以凝聚模型為基礎(chǔ),對(duì)慢速凝聚向快速凝聚的過渡進(jìn)行了定量描述:當(dāng)i×j固定,σ≥-ln P0/ln(ij)時(shí),大小為i、j的兩個(gè)團(tuán)簇碰撞就要發(fā)生凝聚;當(dāng)σ為一固定的正值,i、j滿足ij≥P0-1/σ時(shí),體系則由慢速凝聚過渡到快速凝聚.本研究發(fā)現(xiàn)擴(kuò)散控制凝聚中,不同體系的凝聚具有同一性,其變化規(guī)律是一致的.

    總之,模擬研究表明,膠體體系的凝聚實(shí)際是一種分形凝聚,σ>0時(shí),其凝聚動(dòng)力學(xué)過程是一個(gè)存在正反饋機(jī)制的非線性動(dòng)力學(xué)過程,σ<0時(shí)則具有負(fù)反饋的特性.顯然在進(jìn)行計(jì)算機(jī)模擬時(shí)必須綜合考慮團(tuán)簇的移動(dòng)特性和碰撞特性,才能較好地反映膠體凝聚的動(dòng)力學(xué)過程.

    (1)Ziff,R.M.;McGrady,E.;Meakin,P.J.Chem.P hys.1985,82(11),5269.doi:10.1063/1.448600

    (2)Rottereau,M.;Gimel,J.C.;Nicolai,T.;Durand,D.Eur.Phys.J.E 2004,15(2),133.doi:10.1140/epje/i2004-10044-x

    (3)Zhang,H.L.;Zhou,G.Journal of Fuel Chemistry and Technology 1998,26(1),79.[張慧寧,周 革.燃料化學(xué)學(xué)報(bào),1998,26(1),79.]

    (4)Li,J.S.;Bian X.X.Ch emistry 1998,61(11),9.[李敬生,邊選霞.化學(xué)通報(bào),1998,61(11),9.]

    (5)Evilevitch,A.;Rescic,J.;J?nsson,B.;Olsson,U.J.P hys.Chem.B 2002,106(45),11746.doi:10.1021/jp020467r

    (6)Meakin,P.;Vicsek,T.;Family,F.Phys.Rev.B 1985,31(1),564.doi:10.1103/PhysRevB.31.564

    (7)Family,F.;Meakin,P.;Vicsek,T.J.Chem.Phys.1985,83(8),4144.doi:10.1063/1.449079

    (8)Meakin,P.;Family,F.P hys.Rev.A 1988,38(4),2110.doi:10.1103/PhysRevA.38.2110

    (10)Zheng,D.M.;Zhu,S.G.;Jiang,W.T.;Dai,J.Journal of Nanchan g University(E ngineering&Technology)2011,33(2),129.[鄭典模,朱升干,蔣文天,戴 健.南昌大學(xué)學(xué)報(bào):工科版,2011,33(2),129.]

    (11)Liu,D.J.;Qiu,Z.M.Computers and Applied Chemistry 2012,29(4),451.[劉東京,邱祖民.計(jì)算機(jī)與應(yīng)用化學(xué),2012,29(4),451.]

    (12)Gonzalez,A.E.Phys.Rev.Lett.1993,71(14),2248.doi:10.1103/PhysRevLett.71.2248

    (13)Jullien,R.Conpemp.Ph ys.1987,28(5),477.doi:10.1080/00107518708213736

    (14)Haw,M.;Poon,W.;Pusey,P.P hys.Rev.E 1997,56(2),1918.doi:10.1103/PhysRevE.56.1918

    (15)Meakin,P.;Chen,Z.Y.;Deutch,J.J.Chem.P hys.1985,82(8),3786.doi:10.1063/1.448890

    (16)Meakin,P.J.Sol-Gel Sci.Technol.1999,15(2),97.doi:10.1023/A:1008731904082

    (17)Xiong,H.;Li,H.;Chen,W.;Xu,J.;Wu,L.J.Colloid Interface S ci.2010,344(1),37.doi:10.1016/j.jcis.2009.12.051

    (18)Wang,H.;Pan,Y.M.;Cai,Z.S.;Wang,G.C.;Zhao,X.Z.Chem.J.Chin.Univ.1999,20(9),1419.[王 惠,潘蔭明,蔡遵生,王貴昌,趙學(xué)莊.高等學(xué)?;瘜W(xué)學(xué)報(bào),1999,20(9),1419.]

    (19)Wang,H.;Wang,G.C.;Cai,Z.S.;Pan,Y.M.;Zhao,X.Z.Chem.J.Chin.Univ.2000,21(6),922.[王 惠,王貴昌,蔡遵生,潘蔭明,趙學(xué)莊.高等學(xué)校化學(xué)學(xué)報(bào),2000,21(6),922.]

    (20)Sandkühler,P.;Lattuada,M.;Wu,H.;Sefcik,J.;Morbidelli,M.Adv.Colloid Interface S ci.2005,113(2),65.

    (21)Klein,R.;Meakin,P.Nature 1989,339(6),360.

    (22)Loh,N.D.;Hampton,C.Y.;Martin,A.V.;et al.Nature 2012,486(7404),513.doi:10.1038/nature11222

    (23)Robinson,D.;Earnshaw,J.Phys.Rev.A 1992,46(4),2045.doi:10.1103/PhysRevA.46.2045

    (24)Robinson,D.;Earnshaw,J.Phys.Rev.A 1992,46(4),2055.doi:10.1103/PhysRevA.46.2055

    (25)Robinson,D.;Earnshaw,J.Phys.Rev.A 1992,46(4),2065.doi:10.1103/PhysRevA.46.2065

    (26)Zhu,H.L.;Li,H.;Jiao,M.Y.;Tian,R.A cta P edolo gica Sinica 2012,49(3),409.[朱華玲,李 航,賈明云,田 銳.土壤學(xué)報(bào),2012,49(3),409.]

    (27)Zhu,H.L.;Li,B.;Xiong,H.L.;Li,H.;Jia,M.Y.A cta Phys.-Chim.Sin.2009,25(6),1225.[朱華玲,李 兵,熊海靈,李 航,賈明云.物理化學(xué)學(xué)報(bào),2009,25(6),1225.]doi:10.3866/PKU.WHXB20090631

    (28)Wang,H.X.;Ma,J.;Song,Z.B.;Liu,D.Z.Acta Optica Sinica 2011,31(3),284.[王紅霞,馬 進(jìn),宋仔標(biāo),劉代志.光學(xué)學(xué)報(bào),2011,31(3),284.]

    (29)Lin,M.;Lindsay,H.;Weitz,D.;Ball,R.;Klein,R.;Meakin,P.Phys.Rev.A 1990,41(4),2005.doi:10.1103/PhysRevA.41.2005

    (30)Lin,M.;Lindsay,H.;Weitz,D.;Klein,R.;Ball,R.;Meakin,P.J.P hys.:Cond es.Matter 1990,2(13),3093.

    (31)Whitesides,G.M.;Grzybowski,B.Science 2002,295(5564),2418.doi:10.1126/science.1070821

    (32)Xiong,H.L.;Yuan,Y.Z.;Li,H.;Zhu,H.L.;Jiang,X.J.A cta Phys.-Chim.Sin.2007,23(8),1241.[熊海靈,袁勇智,李航,朱華玲,蔣先軍.物理化學(xué)學(xué)報(bào),2007,23(8),1241.]doi:10.3866/PKU.WHXB20070819

    (33)Yuan,Y.Z.;Xiong,H.L.;Li,H.;Zhu,H.L.;Jiang,X.J.A cta Phys.-Chim.Sin.2007,23(5),688.[袁勇智,熊海靈,李航,朱華玲,蔣先軍.物理化學(xué)學(xué)報(bào),2007,23(5),688.]doi:10.1016/S1872-1508(07)60044-8

    猜你喜歡
    體系模型
    一半模型
    構(gòu)建體系,舉一反三
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    探索自由貿(mào)易賬戶體系創(chuàng)新應(yīng)用
    3D打印中的模型分割與打包
    FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
    如何建立長(zhǎng)期有效的培訓(xùn)體系
    “曲線運(yùn)動(dòng)”知識(shí)體系和方法指導(dǎo)
    “三位一體”德育教育體系評(píng)說
    成人午夜精彩视频在线观看| 在线观看美女被高潮喷水网站| 人妻人人澡人人爽人人| 一边亲一边摸免费视频| 国产男人的电影天堂91| 国产精品三级大全| 黄片无遮挡物在线观看| 免费久久久久久久精品成人欧美视频 | 五月玫瑰六月丁香| 桃花免费在线播放| 国产福利在线免费观看视频| 在线精品无人区一区二区三| 成年人免费黄色播放视频| 欧美日韩一区二区视频在线观看视频在线| 成人亚洲精品一区在线观看| 久久亚洲国产成人精品v| 丝袜在线中文字幕| 亚洲五月色婷婷综合| 久久国产亚洲av麻豆专区| 欧美精品高潮呻吟av久久| 免费在线观看黄色视频的| 黑人巨大精品欧美一区二区蜜桃 | 国产精品久久久久久久久免| 天堂俺去俺来也www色官网| 老司机亚洲免费影院| 亚洲综合色惰| 免费黄色在线免费观看| 亚洲精品456在线播放app| 9191精品国产免费久久| 精品久久久精品久久久| 国语对白做爰xxxⅹ性视频网站| 好男人视频免费观看在线| 一区二区三区乱码不卡18| 九色亚洲精品在线播放| 97人妻天天添夜夜摸| videossex国产| 国产欧美日韩综合在线一区二区| 亚洲天堂av无毛| www.av在线官网国产| 国产免费一级a男人的天堂| 伦理电影大哥的女人| 亚洲,欧美,日韩| 欧美成人精品欧美一级黄| 18禁观看日本| 亚洲国产毛片av蜜桃av| 青春草国产在线视频| 日本欧美国产在线视频| 成人漫画全彩无遮挡| 国产片内射在线| 一级毛片 在线播放| 如日韩欧美国产精品一区二区三区| 男的添女的下面高潮视频| 黄片无遮挡物在线观看| av免费观看日本| 成人国产麻豆网| 1024视频免费在线观看| 亚洲av.av天堂| av片东京热男人的天堂| 一级片'在线观看视频| 久久久精品94久久精品| 草草在线视频免费看| 各种免费的搞黄视频| 国产激情久久老熟女| 国产欧美另类精品又又久久亚洲欧美| 国产 一区精品| 日本午夜av视频| 色婷婷av一区二区三区视频| 少妇被粗大猛烈的视频| 午夜av观看不卡| 成人手机av| 免费在线观看黄色视频的| 国产亚洲欧美精品永久| 久久人人爽av亚洲精品天堂| 热re99久久国产66热| 亚洲欧美中文字幕日韩二区| a级毛片黄视频| 国产乱来视频区| 亚洲图色成人| 久久国产精品男人的天堂亚洲 | 免费高清在线观看日韩| 亚洲在久久综合| 久久狼人影院| 久久毛片免费看一区二区三区| 人人妻人人爽人人添夜夜欢视频| 亚洲精华国产精华液的使用体验| 成年人免费黄色播放视频| 国产精品久久久久久久电影| 全区人妻精品视频| 国产精品国产三级专区第一集| 欧美精品一区二区免费开放| 五月伊人婷婷丁香| 免费播放大片免费观看视频在线观看| 中国三级夫妇交换| 国产综合精华液| 丰满迷人的少妇在线观看| 91在线精品国自产拍蜜月| 国产激情久久老熟女| 久久热在线av| 日韩欧美精品免费久久| 2021少妇久久久久久久久久久| 午夜激情av网站| 中文字幕免费在线视频6| 国产精品蜜桃在线观看| 男女午夜视频在线观看 | 高清黄色对白视频在线免费看| 伦理电影大哥的女人| 午夜福利网站1000一区二区三区| 国产女主播在线喷水免费视频网站| 日产精品乱码卡一卡2卡三| 日韩伦理黄色片| 黑人猛操日本美女一级片| 精品少妇久久久久久888优播| 在线免费观看不下载黄p国产| 纯流量卡能插随身wifi吗| 国产色爽女视频免费观看| 男人添女人高潮全过程视频| 精品人妻熟女毛片av久久网站| 寂寞人妻少妇视频99o| 免费人妻精品一区二区三区视频| 亚洲av中文av极速乱| 超色免费av| 精品一区二区免费观看| 亚洲欧洲精品一区二区精品久久久 | 亚洲欧洲精品一区二区精品久久久 | 中文字幕最新亚洲高清| 最后的刺客免费高清国语| 亚洲精品日本国产第一区| 国产精品久久久久久精品古装| 亚洲av在线观看美女高潮| a级片在线免费高清观看视频| 亚洲经典国产精华液单| 午夜福利网站1000一区二区三区| 国产精品偷伦视频观看了| 国产精品成人在线| 国产成人免费观看mmmm| 久久精品久久精品一区二区三区| 色婷婷久久久亚洲欧美| 天天躁夜夜躁狠狠躁躁| 大香蕉久久网| 天天躁夜夜躁狠狠久久av| 亚洲天堂av无毛| 一边亲一边摸免费视频| 人妻人人澡人人爽人人| 99热6这里只有精品| 男男h啪啪无遮挡| 亚洲五月色婷婷综合| 亚洲少妇的诱惑av| 国产成人一区二区在线| 插逼视频在线观看| 这个男人来自地球电影免费观看 | 97精品久久久久久久久久精品| 国产深夜福利视频在线观看| 国产成人精品福利久久| 2022亚洲国产成人精品| 爱豆传媒免费全集在线观看| 一二三四中文在线观看免费高清| 亚洲精品国产av成人精品| 超碰97精品在线观看| 久热这里只有精品99| 成人免费观看视频高清| 一区二区三区精品91| 黑人高潮一二区| 久久97久久精品| 999精品在线视频| 两个人免费观看高清视频| 午夜福利在线观看免费完整高清在| 看非洲黑人一级黄片| 男女高潮啪啪啪动态图| 国产在视频线精品| 中文字幕免费在线视频6| 久久久欧美国产精品| 久久女婷五月综合色啪小说| 日韩中字成人| 搡老乐熟女国产| 亚洲国产精品成人久久小说| av在线老鸭窝| 你懂的网址亚洲精品在线观看| 99香蕉大伊视频| 国产亚洲最大av| 精品久久久精品久久久| 韩国av在线不卡| 高清毛片免费看| 日韩电影二区| 女性生殖器流出的白浆| 如日韩欧美国产精品一区二区三区| 日韩中文字幕视频在线看片| 最后的刺客免费高清国语| 国产精品久久久久久久久免| 国产精品国产av在线观看| 男女下面插进去视频免费观看 | 亚洲五月色婷婷综合| 精品国产一区二区三区四区第35| 国产高清不卡午夜福利| 婷婷色综合大香蕉| 精品第一国产精品| 建设人人有责人人尽责人人享有的| 天天躁夜夜躁狠狠躁躁| 国产一区二区三区综合在线观看 | 欧美日韩av久久| 国产亚洲精品第一综合不卡 | 97在线视频观看| 99久久精品国产国产毛片| 亚洲av综合色区一区| 久久国产精品男人的天堂亚洲 | 嫩草影院入口| 97在线人人人人妻| 欧美激情极品国产一区二区三区 | 在现免费观看毛片| 精品99又大又爽又粗少妇毛片| 国产精品不卡视频一区二区| 成年美女黄网站色视频大全免费| 亚洲美女视频黄频| 青青草视频在线视频观看| 国产日韩欧美亚洲二区| 亚洲一区二区三区欧美精品| 制服人妻中文乱码| 日韩免费高清中文字幕av| 人人妻人人澡人人看| 男男h啪啪无遮挡| 涩涩av久久男人的天堂| 免费av不卡在线播放| 欧美日韩视频精品一区| 国产高清国产精品国产三级| 国产成人aa在线观看| 日韩精品免费视频一区二区三区 | 欧美精品一区二区大全| 制服人妻中文乱码| 国产高清三级在线| 如日韩欧美国产精品一区二区三区| 大片免费播放器 马上看| 少妇的逼水好多| 久久这里只有精品19| 一级黄片播放器| 大香蕉久久网| 夫妻午夜视频| 久久精品国产a三级三级三级| 男的添女的下面高潮视频| av卡一久久| 男女无遮挡免费网站观看| 久久精品久久久久久噜噜老黄| 久久久久网色| 夫妻性生交免费视频一级片| av在线播放精品| 日韩视频在线欧美| 十分钟在线观看高清视频www| 国产日韩欧美视频二区| 久久99一区二区三区| 亚洲精品久久午夜乱码| 免费看不卡的av| 一本色道久久久久久精品综合| 久久人人爽人人片av| 人人妻人人澡人人爽人人夜夜| 色哟哟·www| 精品99又大又爽又粗少妇毛片| 久久人妻熟女aⅴ| 日本爱情动作片www.在线观看| 丝袜喷水一区| 国产精品成人在线| 亚洲欧美中文字幕日韩二区| 成人亚洲欧美一区二区av| 国产亚洲最大av| 欧美精品高潮呻吟av久久| 少妇被粗大猛烈的视频| 麻豆精品久久久久久蜜桃| 国产男女超爽视频在线观看| 啦啦啦中文免费视频观看日本| 欧美97在线视频| 一区二区av电影网| 天天躁夜夜躁狠狠躁躁| 欧美人与性动交α欧美精品济南到 | 国产国拍精品亚洲av在线观看| 最后的刺客免费高清国语| 亚洲av在线观看美女高潮| videosex国产| 精品国产一区二区久久| 咕卡用的链子| 少妇精品久久久久久久| 亚洲欧美精品自产自拍| 精品久久蜜臀av无| 精品一区二区三区视频在线| 亚洲国产日韩一区二区| 欧美日韩视频高清一区二区三区二| 国产高清三级在线| av免费在线看不卡| 久久久久久久大尺度免费视频| 欧美精品一区二区大全| 在线 av 中文字幕| 成年美女黄网站色视频大全免费| 丝袜美足系列| 久久精品人人爽人人爽视色| 国产av国产精品国产| 国产av一区二区精品久久| 国产精品国产三级国产av玫瑰| 欧美3d第一页| 一区二区三区精品91| 国产精品久久久久久精品电影小说| 欧美激情极品国产一区二区三区 | 精品少妇黑人巨大在线播放| 9色porny在线观看| av视频免费观看在线观看| 久久久久视频综合| 乱码一卡2卡4卡精品| av在线观看视频网站免费| 99久国产av精品国产电影| 国产在线视频一区二区| 久久综合国产亚洲精品| 久热这里只有精品99| xxx大片免费视频| 久久人人97超碰香蕉20202| 黄片无遮挡物在线观看| 欧美3d第一页| 午夜av观看不卡| 秋霞在线观看毛片| av网站免费在线观看视频| 国语对白做爰xxxⅹ性视频网站| 国产探花极品一区二区| 日韩av免费高清视频| 少妇精品久久久久久久| 日韩中文字幕视频在线看片| 一级爰片在线观看| 中文精品一卡2卡3卡4更新| 国产乱人偷精品视频| av在线观看视频网站免费| 在线 av 中文字幕| 纵有疾风起免费观看全集完整版| 国语对白做爰xxxⅹ性视频网站| 十八禁网站网址无遮挡| 最近中文字幕高清免费大全6| 亚洲精品久久成人aⅴ小说| 亚洲一级一片aⅴ在线观看| 亚洲国产成人一精品久久久| 最黄视频免费看| 妹子高潮喷水视频| 免费观看a级毛片全部| 建设人人有责人人尽责人人享有的| 国产免费一级a男人的天堂| 国产免费一区二区三区四区乱码| 亚洲中文av在线| 啦啦啦中文免费视频观看日本| 热99久久久久精品小说推荐| 久久影院123| 国产av一区二区精品久久| 国产男女内射视频| 日韩欧美一区视频在线观看| 久久久精品94久久精品| 日韩视频在线欧美| 精品久久国产蜜桃| 老女人水多毛片| 最黄视频免费看| 香蕉精品网在线| 男人操女人黄网站| 欧美日韩亚洲高清精品| 又大又黄又爽视频免费| 亚洲国产精品国产精品| 视频区图区小说| 熟妇人妻不卡中文字幕| 最新中文字幕久久久久| 国产精品久久久久久精品古装| a级毛片在线看网站| 国产精品一区www在线观看| 国产av码专区亚洲av| 精品卡一卡二卡四卡免费| 少妇高潮的动态图| 黄片播放在线免费| 91精品三级在线观看| 人人妻人人添人人爽欧美一区卜| 一区二区日韩欧美中文字幕 | 久久久久久人人人人人| 99热国产这里只有精品6| 国产精品一二三区在线看| 日本黄大片高清| 99久久精品国产国产毛片| 黄色一级大片看看| 欧美日韩视频高清一区二区三区二| 日韩一区二区视频免费看| 各种免费的搞黄视频| 欧美精品高潮呻吟av久久| 免费不卡的大黄色大毛片视频在线观看| 国产成人欧美| 看非洲黑人一级黄片| 在线免费观看不下载黄p国产| 日本黄大片高清| 女性被躁到高潮视频| 男女下面插进去视频免费观看 | 男女下面插进去视频免费观看 | 亚洲av欧美aⅴ国产| 久久久久网色| 亚洲一区二区三区欧美精品| 另类精品久久| 看非洲黑人一级黄片| 青春草国产在线视频| 免费高清在线观看视频在线观看| 最近的中文字幕免费完整| 免费女性裸体啪啪无遮挡网站| 午夜福利视频精品| 免费av不卡在线播放| 美国免费a级毛片| 欧美日韩综合久久久久久| 中文字幕av电影在线播放| av有码第一页| 亚洲精品国产av成人精品| 免费女性裸体啪啪无遮挡网站| 99久久综合免费| 久久国产精品大桥未久av| 久久精品国产综合久久久 | 精品少妇内射三级| videosex国产| 大香蕉97超碰在线| 母亲3免费完整高清在线观看 | 婷婷色综合www| 永久网站在线| 又黄又粗又硬又大视频| 满18在线观看网站| 久久久久国产精品人妻一区二区| 国产成人精品无人区| 最近2019中文字幕mv第一页| 国产精品 国内视频| 欧美成人午夜精品| 99香蕉大伊视频| 国产深夜福利视频在线观看| 黑人高潮一二区| 乱码一卡2卡4卡精品| 亚洲伊人久久精品综合| 免费黄频网站在线观看国产| 亚洲内射少妇av| 国产色婷婷99| 国产成人精品婷婷| 亚洲精品一二三| 国产综合精华液| 九草在线视频观看| 一本久久精品| 国产永久视频网站| 久久99热6这里只有精品| 午夜福利网站1000一区二区三区| 成年人午夜在线观看视频| 最近最新中文字幕大全免费视频 | 自线自在国产av| 免费人妻精品一区二区三区视频| 成年人午夜在线观看视频| 日韩av不卡免费在线播放| 欧美精品亚洲一区二区| 欧美日韩成人在线一区二区| 国产爽快片一区二区三区| 亚洲av.av天堂| 精品亚洲乱码少妇综合久久| 又粗又硬又长又爽又黄的视频| 男人爽女人下面视频在线观看| 成年美女黄网站色视频大全免费| 美女国产高潮福利片在线看| 亚洲av中文av极速乱| 亚洲精品成人av观看孕妇| 街头女战士在线观看网站| 久久久久久久大尺度免费视频| 国产 一区精品| 亚洲精品中文字幕在线视频| 草草在线视频免费看| 美女内射精品一级片tv| 欧美日本中文国产一区发布| 丝袜美足系列| 国产女主播在线喷水免费视频网站| 热99国产精品久久久久久7| 伊人亚洲综合成人网| 国产乱来视频区| 又大又黄又爽视频免费| 天天影视国产精品| 欧美性感艳星| 亚洲成av片中文字幕在线观看 | 精品卡一卡二卡四卡免费| 亚洲第一区二区三区不卡| 老司机影院成人| 成年美女黄网站色视频大全免费| www日本在线高清视频| 丝袜人妻中文字幕| 天天躁夜夜躁狠狠躁躁| 在线观看三级黄色| 国产成人一区二区在线| xxxhd国产人妻xxx| 国产熟女欧美一区二区| av不卡在线播放| a级毛片在线看网站| 国产无遮挡羞羞视频在线观看| 超色免费av| 午夜精品国产一区二区电影| 国产一区亚洲一区在线观看| 亚洲国产毛片av蜜桃av| 久久ye,这里只有精品| 九草在线视频观看| www日本在线高清视频| 国产精品国产三级专区第一集| 青春草国产在线视频| 国产精品秋霞免费鲁丝片| 在线观看人妻少妇| 成人影院久久| videosex国产| 国产精品久久久久久精品电影小说| 亚洲av欧美aⅴ国产| 国产无遮挡羞羞视频在线观看| 日本色播在线视频| 成年av动漫网址| 亚洲av成人精品一二三区| 美女中出高潮动态图| 黄片播放在线免费| 日本免费在线观看一区| 欧美成人午夜精品| 精品久久久久久电影网| 午夜视频国产福利| 久久午夜综合久久蜜桃| 欧美日韩综合久久久久久| 国产精品久久久久成人av| 亚洲精品日韩在线中文字幕| 久久久久久久大尺度免费视频| 免费播放大片免费观看视频在线观看| 在线观看人妻少妇| 久久ye,这里只有精品| 人人妻人人澡人人爽人人夜夜| 看免费av毛片| 国产av码专区亚洲av| 一二三四在线观看免费中文在 | 欧美人与性动交α欧美软件 | 国产69精品久久久久777片| 国产精品一国产av| 在线观看人妻少妇| av天堂久久9| 亚洲中文av在线| 9热在线视频观看99| 欧美人与性动交α欧美精品济南到 | av女优亚洲男人天堂| 婷婷色麻豆天堂久久| 精品99又大又爽又粗少妇毛片| 99热网站在线观看| 久久久a久久爽久久v久久| 精品午夜福利在线看| 免费人成在线观看视频色| 久久午夜综合久久蜜桃| 精品国产露脸久久av麻豆| 少妇 在线观看| 日本wwww免费看| 成年动漫av网址| 久久97久久精品| 国产成人精品久久久久久| 欧美精品亚洲一区二区| 国产欧美日韩综合在线一区二区| 自线自在国产av| 日本猛色少妇xxxxx猛交久久| 国产精品国产av在线观看| 啦啦啦视频在线资源免费观看| 国产精品一区二区在线观看99| 精品卡一卡二卡四卡免费| a级毛片在线看网站| 黑人巨大精品欧美一区二区蜜桃 | 97超碰精品成人国产| av在线播放精品| 伦理电影大哥的女人| 日本欧美国产在线视频| 亚洲精品视频女| 亚洲伊人色综图| 黑人猛操日本美女一级片| 亚洲成色77777| 麻豆乱淫一区二区| 26uuu在线亚洲综合色| 人妻 亚洲 视频| 大片免费播放器 马上看| 亚洲一区二区三区欧美精品| 日日啪夜夜爽| 亚洲欧洲精品一区二区精品久久久 | 国产免费视频播放在线视频| 22中文网久久字幕| 韩国高清视频一区二区三区| 亚洲成国产人片在线观看| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 成人无遮挡网站| 只有这里有精品99| 91成人精品电影| 日本av手机在线免费观看| 午夜福利网站1000一区二区三区| 建设人人有责人人尽责人人享有的| 久久精品国产自在天天线| videos熟女内射| 免费大片18禁| 亚洲欧美一区二区三区黑人 | 久久婷婷青草| 热re99久久国产66热| 这个男人来自地球电影免费观看 | 欧美成人午夜精品| 亚洲精品乱久久久久久| 久久久久久久久久成人| 精品卡一卡二卡四卡免费| 婷婷色麻豆天堂久久| 涩涩av久久男人的天堂| 亚洲精品久久午夜乱码| av网站免费在线观看视频| 高清毛片免费看| 国产成人精品婷婷| 国产成人精品在线电影| 麻豆乱淫一区二区| 亚洲国产av影院在线观看| 欧美激情 高清一区二区三区| 18禁在线无遮挡免费观看视频| 精品福利永久在线观看| 久久人人97超碰香蕉20202| 亚洲av电影在线观看一区二区三区| 波野结衣二区三区在线| 日本欧美国产在线视频| 一区二区三区精品91| 久久狼人影院| 国产免费一区二区三区四区乱码| 中文天堂在线官网| 国产日韩欧美亚洲二区| 天美传媒精品一区二区| 国产有黄有色有爽视频| 亚洲精品456在线播放app| 99久久中文字幕三级久久日本| 亚洲欧美成人精品一区二区| 国产成人a∨麻豆精品| 精品熟女少妇av免费看|