鏈段剛性對(duì)非溶劑致相分離成膜過(guò)程影響的耗散粒子動(dòng)力學(xué)模擬

唐元暉，李春玉，林亞凱，張春暉，劉澤，余立新，王海輝，王曉琳

（1.中國(guó)礦業(yè)大學(xué)(北京)化學(xué)與環(huán)境工程學(xué)院,北京 100083；2.清華大學(xué)化學(xué)工程系,膜科學(xué)與工程北京市重點(diǎn)實(shí)驗(yàn)室,北京 100084）

近年來(lái)，以膜作為分離過(guò)程核心組件的膜分離技術(shù)在海水淡化、工業(yè)廢水及生活污水處理等水處理領(lǐng)域得到越來(lái)越廣泛的應(yīng)用.非溶劑致相分離（NIPS）法，通常也稱(chēng)為浸沒(méi)沉淀法，是20世紀(jì)60年代就被提出的多孔膜制備方法，其過(guò)程是在常溫下將聚合物溶液浸沒(méi)在聚合物的非溶劑（凝固?。┲?，使聚合物溶液內(nèi)的溶劑與非溶劑發(fā)生跨界面的雙擴(kuò)散傳質(zhì)過(guò)程，誘導(dǎo)聚合物溶液相分離并形成非對(duì)稱(chēng)膜結(jié)構(gòu).通過(guò)NIPS法制備的聚合物膜，微觀上通常是由疏松通透的多孔斷面主體和致密聚合物表層組成的非對(duì)稱(chēng)結(jié)構(gòu)，是目前應(yīng)用最廣泛的超濾/微濾膜材料［1］.其中，非溶劑與溶劑間傳質(zhì)誘導(dǎo)的相分離對(duì)最終聚合物的孔隙率、孔徑分布等膜的微觀結(jié)構(gòu)特征以及膜的強(qiáng)度和滲透性等性能有著決定性的作用［2］.目前對(duì)NIPS過(guò)程得到的膜結(jié)構(gòu)的優(yōu)化和膜性能的提升主要依賴(lài)于實(shí)驗(yàn)的反復(fù)嘗試，缺乏對(duì)相分離過(guò)程本身復(fù)雜的熱力學(xué)和動(dòng)力學(xué)因素的深入理解.為了更好地掌握操作參數(shù)與膜微觀結(jié)構(gòu)之間的聯(lián)系，從20世紀(jì)80年代開(kāi)始，Wang等［3］采用模型和模擬的方法來(lái)描述和研究NIPS過(guò)程，多數(shù)是基于Flory-Huggins熱力學(xué)理論和傳遞過(guò)程原理來(lái)對(duì)膜形成過(guò)程中的濃度場(chǎng)和溫度場(chǎng)進(jìn)行數(shù)值分析，并結(jié)合體系熱力學(xué)相圖得到系統(tǒng)在相分離中的相變路徑從而實(shí)現(xiàn)對(duì)膜結(jié)構(gòu)的預(yù)測(cè)［4］.不過(guò)傳質(zhì)計(jì)算無(wú)法直觀地展示相分離過(guò)程和膜微觀結(jié)構(gòu)特征.近20年來(lái)，隨著高性能計(jì)算在硬件和軟件上的快速發(fā)展，推動(dòng)了計(jì)算機(jī)模擬方法，包括分子動(dòng)力學(xué)（MD）、蒙特卡洛（MC）及耗散粒子動(dòng)力學(xué)（DPD）等，在流體相分離和流動(dòng)領(lǐng)域的廣泛應(yīng)用，使得可以基于分子或者粒子體系和合適的相互作用力場(chǎng)，直接模擬聚合物溶液的相分離過(guò)程.

DPD是1992年由荷蘭的Hoogerbrugge和Koelma在分子動(dòng)力學(xué)模擬的基礎(chǔ)上提出、1997年由英國(guó)的Groot和Warren完善的一種基于粗粒化分子的流體動(dòng)力學(xué)模擬方法［5］.由于其模擬體系的時(shí)空尺度能達(dá)到微米和毫秒的介觀尺度［4］，自提出后就受到了復(fù)雜流體動(dòng)力學(xué)領(lǐng)域研究者的關(guān)注，在聚合物體系相分離［5，6］、表面活性劑［7］及微通道流動(dòng)［8］等領(lǐng)域均有廣泛應(yīng)用.近年來(lái)，DPD也逐漸被用于研究聚合物溶液體系的NIPS動(dòng)力學(xué)過(guò)程和相應(yīng)的成膜機(jī)理中.

2008年，Lu等［9］率先使用DPD構(gòu)建了虛擬聚合物體系的二維模型來(lái)模擬NIPS過(guò)程.2018年，本課題組［10］針對(duì)NIPS法聚合物/溶劑/非溶劑體系的成膜過(guò)程的早期階段構(gòu)建了DPD三維模型，得到了體系相分離形貌演變及相區(qū)尺寸變化等重要結(jié)果，說(shuō)明DPD能夠成功應(yīng)用于聚合物溶液的NIPS過(guò)程的動(dòng)力學(xué)研究.本課題組［11］針對(duì)聚醚砜（PES）/N-甲基-2-吡咯烷酮（NMP）/水的實(shí)際體系，研究了NMP和水的跨界面?zhèn)髻|(zhì)和聚合物濃度等對(duì)PES相分離早期過(guò)程和膜孔徑結(jié)構(gòu)的影響.

在上述模擬NIPS過(guò)程的研究中，均采用了Groot和Warren在1997年提出的彈簧力來(lái)鏈接聚合物相鄰的重復(fù)單元［5］，其中彈簧力（）的表達(dá)式如下：

式中：K為彈簧常數(shù)，K=4；為聚合物相鄰重復(fù)單元間的距離，DPD模擬中涉及的單位通常按照一定的規(guī)則歸一化使用無(wú)因次量綱［5］，因此本文中DPD力場(chǎng)相關(guān)物理量的單位按照引用文獻(xiàn)的規(guī)則歸一化為1.

在此力場(chǎng)中聚合單元間僅存在距離限制，而缺少角度等限制，導(dǎo)致彈簧力場(chǎng)控制下的聚合物鏈段是很柔軟的.

目前常用的有機(jī)膜原料按聚合物鏈段柔性可分為柔性聚合物（如聚偏氟乙烯及聚丙烯等）和剛性聚合物（如聚砜和聚醚砜等），彈簧力場(chǎng)則主要適用于鏈接柔性聚合物的重復(fù)單元.由圖1可見(jiàn)，PES鏈段中相鄰重復(fù)單元存在空間位阻較大的砜基和苯基結(jié)構(gòu)，這意味著結(jié)構(gòu)單元間的角度是無(wú)法隨意彎曲的，呈現(xiàn)出一定的剛性.然而彈簧力場(chǎng)控制下的PES鏈段是很柔軟的［式（1）］，同時(shí)Dingeman等［12，13］研究發(fā)現(xiàn)，描述聚合物鏈段的力場(chǎng)會(huì)影響聚合物鏈段的徑向分布函數(shù)，進(jìn)而可能改變聚合物的相行為.

Fig.1 Atomic structure of PES(A)and schematics of the MD simulation strategy for the PES/NMP solution formation process(B)

為了更加準(zhǔn)確地模擬PES/NMP/水體系的NIPS過(guò)程并考察PES剛性對(duì)相分離過(guò)程的影響，本文基于量子力學(xué)計(jì)算構(gòu)建了PES/NMP溶液，通過(guò)對(duì)PES單體間距離和角度分布映射的方法將分子動(dòng)力學(xué)模擬所得的原子尺度信息與DPD粒子間的相互作用聯(lián)系起來(lái)，構(gòu)建了符合PES鏈段剛性的簡(jiǎn)諧力場(chǎng)；然后在前期已構(gòu)建的PES/NMP/水體系的DPD模擬平臺(tái)基礎(chǔ)上［11］，對(duì)比簡(jiǎn)諧力場(chǎng)與式（1）的彈簧力場(chǎng)對(duì)聚醚砜溶劑相分離過(guò)程的影響，并考察了不同聚合物濃度下該體系相分離的變化，最后結(jié)合文獻(xiàn)中的實(shí)驗(yàn)結(jié)果得到綜合評(píng)價(jià).

1 模擬方法及模型構(gòu)建

1.1 模擬方法

分子動(dòng)力學(xué)是一種基于牛頓力學(xué)來(lái)模擬由原子構(gòu)成的分子體系的動(dòng)態(tài)行為或宏觀性質(zhì)的方法.在分子動(dòng)力學(xué)模擬中，原子和分子間的相互作用分為成鍵相互作用（如鍵長(zhǎng)勢(shì)能、鍵角勢(shì)能等）和非鍵相互作用（如范德華相互作用、氫鍵等）.通常將這些相互作用以數(shù)學(xué)形式的力場(chǎng)進(jìn)行表達(dá)，常用的力場(chǎng)有GAFF力場(chǎng)［14］、AMBER力場(chǎng)［15］及CHARMM力場(chǎng)［16］等.本文根據(jù)量子力學(xué)采用Gaussian軟件［17］構(gòu)建PES/NMP溶液的全原子體系，使用Amber軟件［18～21］開(kāi)展PES溶液的MD模擬，選擇的力場(chǎng)是GAFF.由于本文的重點(diǎn)在于MD與DPD之間的映射和DPD模擬，關(guān)于MD的模擬細(xì)節(jié)請(qǐng)參考文獻(xiàn)［22，23］.

DPD是一種粗?；姆肿觿?dòng)力學(xué)，該方法將微觀粒子粗粒化為一定量原子或分子的集合.為了便于運(yùn)算，DPD模擬通常采用約化單位［5］，即用于描述DPD力場(chǎng)的物理量的單位均為1.在DPD體系中，每個(gè)粒子具有相同的質(zhì)量，粒子質(zhì)量中心的運(yùn)動(dòng)符合牛頓運(yùn)動(dòng)方程，且所有的DPD粒子都會(huì)受到3種力的作用［24］：保守力（FC）、耗散力（FD）和隨機(jī)力（FR），它們都是粒子對(duì)間在截?cái)喟霃剑╮c）距離內(nèi)的相互作用，其中保守力是用于表征體系特性的力，其表達(dá)式如下：

式中：αij為粒子間的排斥力參數(shù)；耗散力體現(xiàn)了體系黏度，隨機(jī)力反映了熱漲落的影響，表達(dá)式分別為：

式中：ωD和ωR是依賴(lài)于r的權(quán)重函數(shù)，并且當(dāng)rij＞rc時(shí)其值均為是一個(gè)符合高斯正態(tài)分布的隨機(jī)變量；參數(shù)γ和σ分別是耗散力和隨機(jī)力的系數(shù).耗散力和隨機(jī)力也是作用于兩個(gè)粒子質(zhì)心的連線，同時(shí)保證動(dòng)量守恒，耗散力使體系粒子的溫度降低，隨機(jī)力給體系注入能量，耗散力和隨機(jī)力分別起到了摩擦作用和噪聲作用，兩者結(jié)合起來(lái)成為一個(gè)熱浴.

1.2 MD到DPD的映射

為了構(gòu)建更加準(zhǔn)確的DPD力場(chǎng)來(lái)描述聚合物鏈段，近年來(lái)研究者們開(kāi)展了很多研究來(lái)更精準(zhǔn)地構(gòu)建從全原子到DPD粒子體系的粗?；成?Lu等［25］從MD模擬中獲得聚乙烯鏈段的徑向分布函數(shù)，進(jìn)而得到聚乙烯重復(fù)單元間的距離和角度分布，而后在DPD體系中依據(jù)距離和角度分布以及Rosenbluth規(guī)則進(jìn)行原子再生長(zhǎng)，從而得到結(jié)構(gòu)更加合理的聚乙烯鏈段；Ortiz等［26］從MD模擬中獲得了環(huán)氧乙烷-聚乙二烯（PEO-PEE）嵌段共聚物的均方根偏差（RMSD），而后將RMSD映射到DPD體系，進(jìn)而得到了符合PEO-PEE嵌段共聚物結(jié)構(gòu)特征的力場(chǎng).根據(jù)上述研究結(jié)果，本文將通過(guò)MD模擬獲取PES單體質(zhì)心間距離和角度的分布，并將其映射到DPD體系，從而構(gòu)建符合PES鏈段剛性的簡(jiǎn)諧力場(chǎng).

1.2.1 MD體系構(gòu)建 通過(guò)MD對(duì)PES/NMP在室溫下形成的溶液開(kāi)展模擬，得到熱力學(xué)平衡態(tài)，從而獲得PES鏈段在NMP中PES重復(fù)單元間的距離和角度的平衡分布結(jié)構(gòu)特征.為了構(gòu)建PES/NMP的全原子體系，本文基于量子化學(xué)方法，使用Gaussian軟件構(gòu)建PES單體以及NMP分子并對(duì)其進(jìn)行優(yōu)化，其中結(jié)構(gòu)優(yōu)化時(shí)選擇了MP2方法，基組為6-31g（d），靜電荷優(yōu)化時(shí)采用Merz-Kollman方法得到了擬合靜電勢(shì)電荷.基于此優(yōu)化的結(jié)構(gòu)采用Amber軟件構(gòu)建出聚合度為30的PES鏈，并在真空中進(jìn)行能量極小化和分子動(dòng)力學(xué)平衡，從而得到結(jié)構(gòu)合理且松弛的PES聚合物鏈段.最后利用Packmol軟件［27］構(gòu)建聚合物質(zhì)量分?jǐn)?shù)為16%的PES/NMP溶液盒子，盒子采用周期性邊界條件，由于MD體系時(shí)空尺度的限制，將溶液設(shè)置為僅含有5條重復(fù)單元數(shù)為30的PES鏈，溶液構(gòu)建過(guò)程如圖1所示，即將PES鏈和一定數(shù)量的NMP分子放入模擬盒子中，通過(guò)Amber軟件進(jìn)行能量極小化和分子動(dòng)力學(xué)過(guò)程，并在溫度為298 K條件下運(yùn)行100 ns的NPT系綜以保證體系達(dá)到穩(wěn)定狀態(tài)，從而得到穩(wěn)定的PES/NMP溶液，最終該體系總能量、溫度、壓力及密度均達(dá)到穩(wěn)定狀態(tài)（圖2）.

Fig.2 Evolutions of the density and pressure(A),energy and temperature(B)of PES/NMP solution in the NPT ensemble during 100 ns

本文中分子動(dòng)力學(xué)過(guò)程所采用的力場(chǎng)為Amber中的GAFF力場(chǎng)，其表達(dá)式如下［14］：

式中：Vpair為體系總勢(shì)能；r和θ分別為粒子間距離和角度；req和θeq分別為平衡距離和平衡角度；Kr和Kθ分別為鍵勢(shì)能常數(shù)和角勢(shì)能常數(shù)；Rij為粒子i和粒子j質(zhì)心間的距離；Aij，Bij和ε均為勢(shì)能參數(shù)；qi和qj分別為粒子i和粒子j所帶電荷數(shù).由式（5）可以看出，GAFF力場(chǎng)主要由鍵長(zhǎng)相互作用、鍵角相互作用以及非鍵相互作用三部分組成.

1.2.2 DPD體系構(gòu)建 從圖3（A）可以發(fā)現(xiàn)，PES單體或NMP分子整體呈平面結(jié)構(gòu)，與DPD粒子的球形結(jié)構(gòu)偏離較大，但從參考文獻(xiàn)［5，7，26］等中均可以發(fā)現(xiàn)，將分子粗?；癁镈PD粒子并沒(méi)有嚴(yán)格遵循分子本身的結(jié)構(gòu)，而是根據(jù)模擬體系本身的物性確定的粗?；桨?因此，本文基于本課題組已發(fā)表的工作［11］，將PES單體和NMP分子的實(shí)際體積映射到DPD模型中：分別將1個(gè)PES單體或2個(gè)NMP分子粗粒化為1個(gè)DPD粒子，1個(gè)DPD粒子所占據(jù)的體積為0.2 nm3.其映射模型分別如圖3（A）所示.在DPD模擬中，長(zhǎng)度單位為截?cái)喟霃剑╮c′，nm），質(zhì)量單位為粒子質(zhì)量（m′，g/mol），能量單位為［（kBT）′，J］（上標(biāo)′為模擬中帶單位的量）.為了便于計(jì)算，通常對(duì)各個(gè)單位進(jìn)行無(wú)量綱化處理［5］，由此得到長(zhǎng)度單位為rc，質(zhì)量單位為m，能量單位為kBT.本文中每個(gè)DPD粒子所占據(jù)的體積均為0.2 nm3，數(shù)密度為3，即單位體積的空間區(qū)域可容納3個(gè)粒子，進(jìn)而得到時(shí)間單位為τ′=其中玻爾茲曼常數(shù)kB=1.38×10-23J/K，參考溫度Tref=298 K，M是所有粒子的平均分子量，根據(jù)本課題組之前的研究結(jié)果［11］，本文中M=203，所以時(shí)間單位τ′=7.9×10-12s，并且在DPD模擬中所使用的時(shí)間步長(zhǎng)Δt=0.01τ′.

在粗?；Ｐ秃蜁r(shí)空尺度確定后，開(kāi)始構(gòu)建PES/NMP溶液的DPD模型體系.PES鏈的聚合度與MD體系均為30，溶液的聚合物質(zhì)量分?jǐn)?shù)也均為16%，由于后續(xù)模擬相分離過(guò)程需要用水作為凝固浴，而NMP的黏度大于水，為了將兩者黏度進(jìn)行區(qū)分，本文將NMP模型設(shè)置為兩個(gè)DPD粒子相連，而水模型為單獨(dú)的DPD粒子（圖3）.然后采用GALAMOST軟件［28］構(gòu)建尺寸為20rc×20rc×20rc的聚合物溶液盒子，盒子采用周期性邊界條件.然后在298 K條件下，對(duì)溶液進(jìn)行30000τ的均一化［圖3（B）］.

Fig.3 Schematics of the coarse-graining from the molecular level to the DPD particles(A)and a snapshot of the DPD equilibrium system of the PES/NMP solution(B)

在PES/NMP溶液的DPD體系中，粒子間的相互作用通過(guò)排斥參數(shù)αij實(shí)現(xiàn)，而排斥參數(shù)αij與Flory-Huggins相互作用參數(shù)χij有關(guān)，其關(guān)系式如下：

式中：當(dāng)溫度為298K時(shí)，kBT=1，相同粒子間的排斥參數(shù)αij=75kBT/ρ=25［5］.PES與NMP間的相互作用參數(shù)χij=0.16，排斥參數(shù)αij=25.56.

此外，對(duì)于構(gòu)成PES鏈的粒子之間的相互作用方式，本文引入了簡(jiǎn)諧鍵勢(shì)能（Harmonic bond potential）和簡(jiǎn)諧角勢(shì)能（Harmonic angle potential），從而取代Warren所提出的彈簧勢(shì)能，其中簡(jiǎn)諧鍵勢(shì)能［Vbond(r)］表達(dá)式為［28］：

式中：kb表示彈簧常數(shù)；rb表示平衡距離.

簡(jiǎn)諧角勢(shì)能［Vangle(θ)］表達(dá)式為［28］：

式中：ka表示勢(shì)能常數(shù)；θa表示平衡角度.

1.2.3 MD到DPD的映射通過(guò)映射的方法將MD模擬所得的原子尺度信息與DPD粒子間的相互作用聯(lián)系起來(lái)，具體方法如圖4所示.首先利用Amber中的cpptraj模塊對(duì)MD體系中PES相鄰重復(fù)單元的質(zhì)心間的距離和角度分布數(shù)據(jù)進(jìn)行統(tǒng)計(jì)及高斯擬合，得到的PES相鄰重復(fù)單元間的距離和角度示意圖如圖4（A）所示；然后在DPD體系中，通過(guò)預(yù)估簡(jiǎn)諧力場(chǎng)中［即式（7）和式（8）］的參數(shù)kb，rb，ka和θa，并按1.2.2節(jié)的方法得到均一穩(wěn)定的PES/NMP溶液，對(duì)構(gòu)成PES相鄰粒子質(zhì)心間的距離和角度分布數(shù)據(jù)進(jìn)行統(tǒng)計(jì)和高斯擬合，這些距離和角度可以近似地理解為DPD粒子間的鍵長(zhǎng)和鍵角；接下來(lái)將擬合得到的均值和方差與分子動(dòng)力學(xué)所得結(jié)果對(duì)比，不斷調(diào)整簡(jiǎn)諧力場(chǎng)參數(shù)，并重復(fù)上述步驟，直到DPD體系所得的鍵長(zhǎng)和鍵角分布與MD體系趨于一致［圖4（B）］，即在該力場(chǎng)參數(shù)條件下DPD體系的PES結(jié)構(gòu)最接近于MD體系.最終擬合得到的兩體系的均值和方差如表1所示，鍵長(zhǎng)和鍵角分布高斯曲線如圖5（A）和（B）所示，本文所構(gòu)建的DPD簡(jiǎn)諧力場(chǎng)表達(dá)式為：

Fig.4 Schematic diagram of the mapping of bond length and angle from MD(A)to DPD(B)

Table 1 Mean and variance of bond lengths and angles in MD and DPD simulation

此外，由表1和圖5可以發(fā)現(xiàn)，DPD體系的鍵角分布均值和方差與MD體系基本一致，這說(shuō)明對(duì)于聚合物鏈鍵角的優(yōu)化效果是較好的.雖然DPD體系的鍵長(zhǎng)分布的均值與MD體系基本一致，但兩者分布的方差間差距較大，這是因?yàn)樵贒PD體系中，鍵長(zhǎng)方差的大小取決于簡(jiǎn)諧鍵勢(shì)能表達(dá)式中參數(shù)kb的大小，kb越大，方差越小.但我們并沒(méi)有過(guò)度提高kb的值，否則會(huì)使簡(jiǎn)諧力場(chǎng)的剛性過(guò)強(qiáng)，降低模擬效率：導(dǎo)致DPD需要更小的步長(zhǎng)來(lái)避免體系發(fā)生可能的崩潰，降低DPD體系可達(dá)到的時(shí)空尺度.

Fig.5 Distribution of bond length(A)and angle(B)between consecutive PES repeat units

1.3 DPD模型構(gòu)建和模擬結(jié)果表征

采用與1.2.2節(jié)相同的方法得到PES/NMP/水體系的粗?；Ｐ停罱K該體系的粗?；瘏?shù)如表2所示，粗?；Ｐ腿鐖D6所示，粒子間的相互作用參數(shù)χij及排斥參數(shù)αij如表3所示.

Table 2 Molecular volume,molar mass,and the number of species per particle for PES,NMP and water

Fig.6 Schematics of the coarse-graining from the molecular level to the DPD particles

Table 3 Flory-Huggins interaction parameter χij and the corresponding DPD repulsio n parameter aij between different species

在DPD模擬階段，本文分別構(gòu)建凝固浴盒子（尺寸為90rc×50rc×90rc）和聚合物盒子（尺寸為90rc×90rc×90rc），并且兩盒子在x和z方向均為周期性邊界條件，y方向設(shè)置了中性基板限制.中性基板是由兩層凍結(jié)的DPD粒子構(gòu)成的，并且近基板區(qū)域設(shè)置了反彈邊界條件，本文中的中性基板與聚合物、溶劑的相互作用參數(shù)均等于10.8，該值已在文獻(xiàn)［5］中被證明能正確描述體系在邊界處的密度分布.

分別對(duì)聚合物和凝固浴進(jìn)行30000τ的均一化，以確保粒子均勻分布在各自的空間內(nèi).將平衡后的兩體系按照?qǐng)D7（A）所示的方式拼接在一起，并移除兩體系中間部分的中性基板［如圖7（A）中的箭頭所示］，此時(shí)該體系的尺寸為90rc×140rc×90rc.圖7（B）顯示了初始狀態(tài)下PES濃度為12%（體積分?jǐn)?shù)）時(shí)各個(gè)粒子的體積分?jǐn)?shù)分布以及總粒子數(shù)密度分布.值得注意的是，體系的總粒子數(shù)密度在y=50rc處略有降低［圖7（B）洋紅色虛線由3.0降至2.8］，這是移除兩體系中間部分的中性基板導(dǎo)致的，但總粒子數(shù)密度的變化幅度較小，因而其對(duì)整體的影響可以忽略.將聚合物溶液與凝固浴間的臨時(shí)壁移除后，使PES溶液與非溶劑（水）充分接觸混合，并對(duì)相分離過(guò)程進(jìn)行10000τ的DPD模擬.

Fig.7 Schematics of the DPD simulation strategy for the membrane formation process via NIPS

采用與前期工作［11］相同的表征方法對(duì)DPD模擬結(jié)果進(jìn)行表征.首先采用回轉(zhuǎn)半徑（Rg）對(duì)PES鏈段剛性進(jìn)行表征，其計(jì)算方法如下：

式中：ri是鏈段i的位置；rcm是聚合物鏈質(zhì)心的位置.由式（11）可知聚合物鏈段剛性越強(qiáng)，其重復(fù)單元距離鏈段質(zhì)心間的距離越大，進(jìn)而導(dǎo)致其均方回轉(zhuǎn)半徑越大.

相分離過(guò)程的形貌表征是通過(guò)程序計(jì)算得到體系中不同區(qū)域的密度分布，并基于Matlab軟件繪制了局部聚合物體積分?jǐn)?shù)大于或等于45%的等值面.本文中DPD體系數(shù)密度為3，而模擬的聚合物溶液的體積分?jǐn)?shù)為8%～16%，均小于33%，這意味著聚合物溶液均一時(shí)單位體積的空間區(qū)域內(nèi)存在0～1個(gè)聚合物粒子，當(dāng)聚合物溶液發(fā)生相分離時(shí)，聚合物富相區(qū)域單位體積空間內(nèi)可聚集2～3個(gè)聚合物粒子，因此通常認(rèn)為局部聚合物體積分?jǐn)?shù)達(dá)到40%～50%即為富相區(qū)域.本文選擇局部聚合物體積分?jǐn)?shù)45%作為聚合物富相與貧相的臨界點(diǎn)，前期工作中［11］同樣采用的是該值.由于體系中存在跨越相界面的粒子傳遞，必然會(huì)造成體系中沿流動(dòng)方向（即y軸）上的局部密度發(fā)生明顯的變化，本文計(jì)算了各個(gè)粒子沿y軸方向的密度分布以及每個(gè)區(qū)域內(nèi)的平均粒子數(shù).聚合物的相區(qū)尺寸是公認(rèn)的相分離過(guò)程的常規(guī)表征手段，本文采用與前期工作［11］相同的計(jì)算方法，有別于常規(guī)均勻體系的相區(qū)尺寸計(jì)算［29］，本文的相區(qū)尺寸去除了如圖7（A）所示的凝固浴部分，僅計(jì)算了鑄膜液一側(cè)的聚合物相區(qū)尺寸隨時(shí)間的變化.

除了相分離過(guò)程的表征，本文以相分離形貌較完善的體系認(rèn)定為可能形成的膜結(jié)構(gòu)，計(jì)算了此時(shí)模擬所得聚合物膜的二維和三維孔徑分布，其中二維孔徑的計(jì)算方法為：首先通過(guò)程序計(jì)算了體系中不同局部區(qū)域的密度場(chǎng)體系的密度場(chǎng)，然后沿y軸掃描，對(duì)于xz二維截面中PES粒子體積分?jǐn)?shù)小于45%的區(qū)域（聚合物貧相區(qū)域）即認(rèn)定為孔徑區(qū)域，通過(guò)統(tǒng)計(jì)沿y軸孔徑區(qū)域所占的面積即得到了聚合物膜的二維孔徑分布，這種二維孔徑分布可以直接地反映體系沿粒子傳質(zhì)方向上體系的不對(duì)稱(chēng)性變化；對(duì)于三維孔徑的獲取，首先去除了凝固浴所在部分的區(qū)域，而后使用Zeo++軟件進(jìn)行分析［30］.

2 結(jié)果與討論

首先在聚合物體積分?jǐn)?shù)為12%的條件下，探究聚合物鏈段剛性對(duì)PES/NMP/水體系相分離過(guò)程的影響；然后在簡(jiǎn)諧力場(chǎng)條件下，探究聚合物體積分?jǐn)?shù)從8%增加至16%時(shí)PES/NMP/水體系相分離過(guò)程的變化.

2.1 聚合物鏈段剛性對(duì)NIPS法成膜過(guò)程的影響

將PES鏈長(zhǎng)設(shè)定為90，聚合物初始體積分?jǐn)?shù)為12%，對(duì)比了兩種力場(chǎng)條件下體系相分離過(guò)程隨時(shí)間的變化.圖8為在4個(gè)時(shí)間步長(zhǎng)（2500τ，5000τ，7500τ和10000τ）下的相分離形貌圖，其中藍(lán)色區(qū)域?yàn)榫植烤酆衔矬w積分?jǐn)?shù)大于或等于45%的等值面.在PES/NMP溶液與水之間的界面上，水和NMP快速進(jìn)行物質(zhì)交換，導(dǎo)致聚合物在界面區(qū)域快速聚集，形成薄而致密的聚合物表層.從圖8可以觀察到，簡(jiǎn)諧力場(chǎng)的成膜速度更快，并且形成的膜也更加致密，在2500τ時(shí)［圖8（A）］已初步形成了聚合物表層.隨著聚合物表層的逐漸形成，界面處的NMP與水的交換速率逐漸變慢，并且由于非溶劑水的推動(dòng)，聚合物表層逐漸向盒子的右側(cè)移動(dòng)，并且隨著時(shí)間的推移，水逐漸擴(kuò)散到溶液內(nèi)部區(qū)域，聚合物層逐漸生長(zhǎng)，且變得更加致密.同時(shí)PES溶液內(nèi)部逐漸發(fā)生聚合物聚結(jié)，由此產(chǎn)生了表層下的多孔亞層結(jié)構(gòu)，這是典型的聚合物富相成核和生長(zhǎng)的特點(diǎn).

為了更加直觀地描述以上現(xiàn)象，并探究力場(chǎng)對(duì)膜結(jié)構(gòu)的影響，計(jì)算了體系中各粒子沿y軸的體積分?jǐn)?shù)分布以及總數(shù)密度分布（圖9），圖9中洋紅色虛線表明該體系是總數(shù)密度為3的粒子均勻分布體系，這與DPD初始設(shè)置的條件相吻合；通過(guò)比較同一力場(chǎng)的不同時(shí)刻N(yùn)MP和水的粒子分布（圖9中紅色和藍(lán)色虛線）可以發(fā)現(xiàn)，初始時(shí)水位于凝固浴側(cè)，NMP位于鑄膜液一側(cè)，隨著時(shí)間的推移兩者逐漸相互擴(kuò)散，并且在PES體積分?jǐn)?shù)分布（圖9中黑色虛線）的峰值處兩者的擴(kuò)散趨勢(shì)明顯減緩，這是由于在界面附近形成了致密的聚合物層導(dǎo)致的.

Fig.8 Morphology evolution of the systems with initial polymer volume fractions 12% in spring force field(A—D)and harmonic force field(E—H)at 2500τ(A,E),5000τ(B,F),7500τ(C,G)and 10000τ(D,H)during the DPD simulations

Fig.9 Profiles of volume fraction of PES(black),NMP(red)and H2O(blue)and total particle density(magenta)along the y axis in spring force field(A—D)and harmonic force field(E—H)at 2500τ(A,E),5000τ(B,F),7500τ(C,G)and 10000τ(D,H)during the DPD simulations

圖10給出不同力場(chǎng)下PES在不同時(shí)刻的體積分?jǐn)?shù)分布.初始的PES濃度在溶液側(cè)的分布是均勻的，由于非溶劑水和溶劑NMP在界面上的傳質(zhì)，使得PES在相界面上逐漸形成了致密層，正好對(duì)應(yīng)于PES體積分?jǐn)?shù)上的峰值，此峰的高度逐漸增加并且向右移動(dòng)，且其峰寬度保持相對(duì)恒定，進(jìn)一步說(shuō)明了PES表層的厚度不變但更為致密.在兩種力場(chǎng)條件下，聚合物形態(tài)隨時(shí)間變化（圖8）以及聚合物體積分?jǐn)?shù)分布演變（圖10）非常相似，這說(shuō)明兩者相分離機(jī)制是相似的.盡管在兩種力場(chǎng)條件下都能觀察到相似的行為，但是各組分的體積分?jǐn)?shù)變化速率存在明顯的差異，簡(jiǎn)諧力場(chǎng)條件下的NMP和水的傳質(zhì)速率以及成膜速度更快.此外，簡(jiǎn)諧力場(chǎng)中PES峰明顯高于彈簧力場(chǎng)，但是峰寬卻略小于彈簧力場(chǎng)，這表明簡(jiǎn)諧力場(chǎng)形成的膜表層更加致密且略薄.

圖11（A）給出不同力場(chǎng)條件下Rg隨時(shí)間的變化.首先，隨時(shí)間的延長(zhǎng)，Rg值迅速減小，簡(jiǎn)諧力場(chǎng)從6.59rc降至6.14rc，彈簧力場(chǎng)從6.59rc降至3.99rc，在1500τ后減小趨勢(shì)變得非常緩慢.這可能是因?yàn)槌跏紩r(shí)NMP與水發(fā)生快速的物質(zhì)交換，導(dǎo)致界面處聚合物快速析出，因而Rg迅速降低；隨后由于膜表層的形成，導(dǎo)致NMP和水傳質(zhì)速率降低，聚合物析出速率同樣降低，所以Rg減小趨勢(shì)也逐漸平緩.此外，還發(fā)現(xiàn)簡(jiǎn)諧力場(chǎng)的Rg值始終大于彈簧力場(chǎng)，這說(shuō)明彈簧力場(chǎng)下的PES鏈段更柔軟，而簡(jiǎn)諧力場(chǎng)在其基礎(chǔ)上增加了對(duì)聚合物結(jié)構(gòu)單元間鍵角的約束，使聚合物剛性更大，這與文獻(xiàn)［31］中的結(jié)論一致.

Fig.10 Profiles of volume fractions of PES along the y axis at spring force field(A)and harmonic force field(B)of different times

Fig.11 Evolution of radius of gyration Rg(A)and PES domain size(B)as functions of time for the systems in harmonic force field and spring force field

聚合物的相區(qū)尺寸變化是衡量聚合物體系相分離程度的常用手段，相區(qū)尺寸增大表明從鑄膜液中析出的聚合物富相在逐漸增長(zhǎng)，即相分離的程度增大.圖11（B）給出PES的相區(qū)尺寸隨時(shí)間變化.可見(jiàn)，隨著相分離的進(jìn)行，相區(qū)尺寸先快速增長(zhǎng)，而后緩慢增加.初始時(shí)增加較快可能是由界面處聚合物迅速成膜引起的，隨后成膜速度降低，導(dǎo)致相區(qū)尺寸增長(zhǎng)也趨于平緩.此外，在整個(gè)過(guò)程中，簡(jiǎn)諧力場(chǎng)的相區(qū)尺寸一直大于彈簧力場(chǎng)，這也進(jìn)一步印證了簡(jiǎn)諧力場(chǎng)條件下聚合物剛性較大，鏈段不容易被壓縮這一結(jié)論.

對(duì)于不對(duì)稱(chēng)的形態(tài)結(jié)構(gòu)，二維（2D）孔徑分布和三維（3D）孔徑分布都是有重要意義的.二維孔徑表示在二維截面中孔徑區(qū)域所占面積的大小，其可以間接地反映孔隙間相互連通的程度；而三維孔徑分布表示的是整體孔徑大小的分布，其與膜的滲透性有直接關(guān)系.由圖10可見(jiàn)，聚合物濃度在y軸方向上達(dá)到一定峰值是存在一個(gè)過(guò)度過(guò)程的，這是由于在界面處部分聚合物鏈段延伸到凝固浴中導(dǎo)致的，該部分在除去NMP和水后不能形成具有顯著分離效果的固體膜結(jié)構(gòu)，因此在聚合物濃度峰值處設(shè)置了一個(gè)截止值［圖12（A）］，僅計(jì)算截止值右側(cè)區(qū)域的孔徑.由圖12（B）可以看出，簡(jiǎn)諧力場(chǎng)和彈簧力場(chǎng)在對(duì)二維孔徑分布的影響上具有相似性，膜表面結(jié)構(gòu)相對(duì)于內(nèi)部更加致密.與彈簧力場(chǎng)相比，簡(jiǎn)諧力場(chǎng)作用下膜表層孔徑略小，這是由于聚合物剛性增強(qiáng)，導(dǎo)致水難以穿透膜表層與NMP進(jìn)行質(zhì)量交換所致，二維孔徑分布表明聚合物鏈段剛性增強(qiáng)有助于形成更加致密的膜表層.圖12（C）為不同力場(chǎng)下的膜的三維孔徑分布圖.與彈簧力場(chǎng)相比，簡(jiǎn)諧力場(chǎng)下的膜的孔徑分布略寬，但兩者的最大孔徑是相近的，這表明力場(chǎng)主要影響的是膜表面孔徑分布，對(duì)內(nèi)部疏松層的影響并不明顯.

Fig.12 Effects of harmonic force field and spring force field on the 2D and 3D pore size distributions

2.2 聚醚砜濃度對(duì)NIPS法成膜過(guò)程的影響

目前，人們普遍認(rèn)為NIPS法成膜過(guò)程中溶劑和非溶劑的質(zhì)量交換作用對(duì)最終膜結(jié)構(gòu)和性能有著決定性的作用［2］，其中溶劑和非溶劑間的傳質(zhì)與溶劑/非溶劑相互作用、聚合物濃度、聚合物鏈長(zhǎng)以及非溶劑溫度等有重要聯(lián)系.本課題組前期的模擬研究發(fā)現(xiàn)［10，11］，在聚合物溶液發(fā)生NIPS過(guò)程中，聚合物濃度對(duì)膜結(jié)構(gòu)的影響更大.因此，本文在簡(jiǎn)諧力場(chǎng)的條件下探究聚醚砜濃度對(duì)相分離過(guò)程的影響，并結(jié)合前期的工作［11］，對(duì)比在簡(jiǎn)諧力場(chǎng)和彈簧力場(chǎng)條件下聚醚砜濃度對(duì)相分離過(guò)程的影響.

由圖13可以看出，同一濃度下，隨著時(shí)間的推移膜形貌的變化與圖8～圖12結(jié)果相似.當(dāng)聚合物濃度較低（8%，體積分?jǐn)?shù)）時(shí)，初始時(shí)聚合物的高濃度區(qū)域較少，且表層更加多孔；當(dāng)聚合物濃度較高（16%，體積分?jǐn)?shù)）時(shí)，由于聚合物鏈更容易發(fā)生聚集，會(huì)增加聚合物的局部濃度，從而可以觀察到更多的高濃度聚合物區(qū)域，且其表層的孔也會(huì)更少.這些觀察結(jié)果與實(shí)驗(yàn)結(jié)果是一致的［32～34］，即膜表層為致密層，內(nèi)部為多孔亞層結(jié)構(gòu).

Fig.13 Morphology evolution of the systems with chain length n=90 and initial polymer volume fractions of 8%(A—D),12%(E—H),and 16%(I—L)at 2500τ(A,E,I),5000τ(B,F,J),7500τ(C,G,K)and 10000τ(D,H,L)during the DPD simulations

圖14表明，隨著時(shí)間的推移，峰的高度逐漸增加并向右移動(dòng)，隨著聚合物初始濃度的增加，峰的寬度也緩慢增長(zhǎng)，這說(shuō)明PES濃度的增加會(huì)導(dǎo)致在界面附近形成了厚且致密的聚合物層，這與實(shí)驗(yàn)觀察結(jié)果一致［33，35］.此外，隨著聚合物濃度的升高，聚合物溶液黏度變得更大，并且所形成的膜表層也會(huì)更加致密，進(jìn)而會(huì)導(dǎo)致溶劑與非溶劑的傳質(zhì)減慢，從而延遲相分離，這一現(xiàn)象在圖13（D），（H）和（L）的相分離形貌中得以體現(xiàn).濃度為16%的膜表層明顯更加致密，但其厚度卻低于濃度為12%的膜表層厚度；這一現(xiàn)象同樣可以在圖14中觀察到，濃度為16%的PES分布密度峰高明顯高于12%的，但其峰寬卻明顯小于12%濃度下的.

Fig.14 Profiles of volume fractions of PES along y axis at 0,2500τ,5000τ,7500τ and 10000τ

圖15（A）給出3種不同的PES濃度條件下，PES的Rg隨時(shí)間的變化.在初始階段，Rg下降得很快，而后又緩慢地減小.這可能是因?yàn)槌跏紩r(shí)聚合物鏈段較為舒展，導(dǎo)致Rg的初始值較大；而隨著非溶劑的接觸，導(dǎo)致了聚合物鏈快速聚集，因而Rg快速降低最后趨于穩(wěn)定.但不同濃度PES溶液的初始Rg也是不同的，濃度為8%，12%和16%時(shí)的初始Rg值分別為5.97rc，5.84rc和5.82rc，即當(dāng)聚合物濃度增加時(shí)，初始Rg值會(huì)略有降低，這與聚合物溶液理論是一致的［36，37］.但是當(dāng)體系逐漸穩(wěn)定后，Rg值卻隨濃度的增大而減小了.這可能是由于低濃度時(shí)，聚合物被溶劑和非溶劑壓縮導(dǎo)致其所占空間相對(duì)較小，因而分子鏈更加團(tuán)聚，所以穩(wěn)定后濃度低時(shí)Rg反而較小.

Fig.15 Radius of gyration Rg(A)and domain size(B)as functions of time for the systems with polymer volume fractions of 8%,12%and 16%

圖15（B）給出不同濃度下聚合物相區(qū)尺寸隨時(shí)間的變化.由圖15（B）可見(jiàn)，隨著時(shí)間的推移，相區(qū)尺寸首先是快速增加，而后增長(zhǎng)速度逐漸減緩.初始時(shí)相區(qū)尺寸快速增加可能是由界面處聚合物層的形成導(dǎo)致的，而后由于膜表層形成導(dǎo)致了物質(zhì)交換阻力變大，因而相區(qū)尺寸增長(zhǎng)幅度減小.此外，PES濃度越高，相區(qū)尺寸越大，這是因?yàn)榫酆衔餄舛仍酱笤饺菀自贜MP與水物質(zhì)交換的時(shí)候析出成膜.

圖16（A）給出不同濃度PES的二維分布孔徑.可見(jiàn)，不同濃度的孔徑分布總體趨勢(shì)是相同的，膜表層的二維孔徑均明顯小于內(nèi)部孔徑，這與實(shí)驗(yàn)結(jié)論一致［34，38］.隨著聚合物濃度從8%增加到16%，膜表層孔徑明顯減小，由7.1rc減小至0.4rc，這與文獻(xiàn)中的實(shí)驗(yàn)結(jié)果符合得很好［33，34］.與前期彈簧力場(chǎng)的模擬結(jié)果相比［11］，簡(jiǎn)諧力場(chǎng)條件下模擬得到的表面孔徑更加接近于實(shí)驗(yàn)所得到的納米級(jí)孔徑［34］，并且在低濃度時(shí)（8%）表現(xiàn)尤為明顯，表明聚醚砜鏈段剛性的增強(qiáng)使得模擬體系更加接近實(shí)際；但內(nèi)部疏松層的二維孔徑卻隨著濃度的升高而減小，這是因?yàn)榫酆衔餄舛壬咴谝欢ǔ潭壬献璧K了溶劑與非溶劑的傳質(zhì)作用，導(dǎo)致了內(nèi)部存在部分溶劑保留的現(xiàn)象，這一現(xiàn)象可以從實(shí)驗(yàn)中觀察到［34］.圖16（B）為不同濃度下膜的三維孔徑分布圖.隨著PES濃度從8%增加到16%，孔徑分布會(huì)明顯變窄，且隨著濃度的增加，三維孔徑分布峰由7.1rc降至3.9rc，與文獻(xiàn)［34］中實(shí)驗(yàn)數(shù)據(jù)的吻合性更好.

Fig.16 Effects of initial polymer volume fraction(8%,12%,16%)on the 2D(x-z plane)pore size profile(A)and 3D pore size distributions(B)

3 結(jié) 論

通過(guò)從MD模擬映射的方法構(gòu)建了符合PES剛性結(jié)構(gòu)的DPD簡(jiǎn)諧力場(chǎng)，證明了不同力場(chǎng)施加下的PES鏈段剛性對(duì)PES/NMP/水體系NIPS過(guò)程存在影響.首先，施加不同力場(chǎng)時(shí)，在整體上均可得到明顯的非對(duì)稱(chēng)結(jié)構(gòu)；其次，PES鏈段剛性的提升能夠明顯提升體系的相分離速度，同時(shí)會(huì)導(dǎo)致相界面處的PES薄層形成得更加快速而致密，孔徑更小，而對(duì)內(nèi)部的疏松結(jié)構(gòu)影響較??；此外，結(jié)合不同力場(chǎng)下聚合物濃度對(duì)相分離過(guò)程的影響發(fā)現(xiàn)，不同PES濃度下，鏈段剛性的提升對(duì)相分離過(guò)程的特征和演變趨勢(shì)沒(méi)有造成根本性的變化，與經(jīng)典的彈簧力場(chǎng)的模擬結(jié)果在整體趨勢(shì)上有相似性，這一點(diǎn)也說(shuō)明了彈簧力場(chǎng)的模擬結(jié)果也有相當(dāng)?shù)膮⒖純r(jià)值.本文中簡(jiǎn)諧力場(chǎng)的構(gòu)建從原理而言并不復(fù)雜，但實(shí)際操作上較為繁瑣且耗時(shí)，因此對(duì)于一些精確度要求不高的體系而言，經(jīng)典的彈簧力場(chǎng)依然是簡(jiǎn)便且有效的方法.