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

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

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

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

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

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

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

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

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

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

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

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

1.1 模擬方法

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

DPD是一種粗?；姆肿觿恿W(xué)，該方法將微觀粒子粗?；癁橐欢吭踊蚍肿拥募?為了便于運算，DPD模擬通常采用約化單位［5］，即用于描述DPD力場的物理量的單位均為1.在DPD體系中，每個粒子具有相同的質(zhì)量，粒子質(zhì)量中心的運動符合牛頓運動方程，且所有的DPD粒子都會受到3種力的作用［24］：保守力（FC）、耗散力（FD）和隨機力（FR），它們都是粒子對間在截斷半徑（rc）距離內(nèi)的相互作用，其中保守力是用于表征體系特性的力，其表達式如下：

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

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

1.2 MD到DPD的映射

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

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

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

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

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

在粗粒化模型和時空尺度確定后，開始構(gòu)建PES/NMP溶液的DPD模型體系.PES鏈的聚合度與MD體系均為30，溶液的聚合物質(zhì)量分數(shù)也均為16%，由于后續(xù)模擬相分離過程需要用水作為凝固浴，而NMP的黏度大于水，為了將兩者黏度進行區(qū)分，本文將NMP模型設(shè)置為兩個DPD粒子相連，而水模型為單獨的DPD粒子（圖3）.然后采用GALAMOST軟件［28］構(gòu)建尺寸為20rc×20rc×20rc的聚合物溶液盒子，盒子采用周期性邊界條件.然后在298 K條件下，對溶液進行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體系中，粒子間的相互作用通過排斥參數(shù)αij實現(xiàn)，而排斥參數(shù)αij與Flory-Huggins相互作用參數(shù)χij有關(guān)，其關(guān)系式如下：

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

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

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

簡諧角勢能［Vangle(θ)］表達式為［28］：

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

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

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體系基本一致，這說明對于聚合物鏈鍵角的優(yōu)化效果是較好的.雖然DPD體系的鍵長分布的均值與MD體系基本一致，但兩者分布的方差間差距較大，這是因為在DPD體系中，鍵長方差的大小取決于簡諧鍵勢能表達式中參數(shù)kb的大小，kb越大，方差越小.但我們并沒有過度提高kb的值，否則會使簡諧力場的剛性過強，降低模擬效率：導(dǎo)致DPD需要更小的步長來避免體系發(fā)生可能的崩潰，降低DPD體系可達到的時空尺度.

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/水體系的粗?；Ｐ?，最終該體系的粗?；瘏?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，該值已在文獻［5］中被證明能正確描述體系在邊界處的密度分布.

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

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

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

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

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

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

2 結(jié)果與討論

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

2.1 聚合物鏈段剛性對NIPS法成膜過程的影響

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

為了更加直觀地描述以上現(xiàn)象，并探究力場對膜結(jié)構(gòu)的影響，計算了體系中各粒子沿y軸的體積分數(shù)分布以及總數(shù)密度分布（圖9），圖9中洋紅色虛線表明該體系是總數(shù)密度為3的粒子均勻分布體系，這與DPD初始設(shè)置的條件相吻合；通過比較同一力場的不同時刻NMP和水的粒子分布（圖9中紅色和藍色虛線）可以發(fā)現(xiàn)，初始時水位于凝固浴側(cè)，NMP位于鑄膜液一側(cè)，隨著時間的推移兩者逐漸相互擴散，并且在PES體積分數(shù)分布（圖9中黑色虛線）的峰值處兩者的擴散趨勢明顯減緩，這是由于在界面附近形成了致密的聚合物層導(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給出不同力場下PES在不同時刻的體積分數(shù)分布.初始的PES濃度在溶液側(cè)的分布是均勻的，由于非溶劑水和溶劑NMP在界面上的傳質(zhì)，使得PES在相界面上逐漸形成了致密層，正好對應(yīng)于PES體積分數(shù)上的峰值，此峰的高度逐漸增加并且向右移動，且其峰寬度保持相對恒定，進一步說明了PES表層的厚度不變但更為致密.在兩種力場條件下，聚合物形態(tài)隨時間變化（圖8）以及聚合物體積分數(shù)分布演變（圖10）非常相似，這說明兩者相分離機制是相似的.盡管在兩種力場條件下都能觀察到相似的行為，但是各組分的體積分數(shù)變化速率存在明顯的差異，簡諧力場條件下的NMP和水的傳質(zhì)速率以及成膜速度更快.此外，簡諧力場中PES峰明顯高于彈簧力場，但是峰寬卻略小于彈簧力場，這表明簡諧力場形成的膜表層更加致密且略薄.

圖11（A）給出不同力場條件下Rg隨時間的變化.首先，隨時間的延長，Rg值迅速減小，簡諧力場從6.59rc降至6.14rc，彈簧力場從6.59rc降至3.99rc，在1500τ后減小趨勢變得非常緩慢.這可能是因為初始時NMP與水發(fā)生快速的物質(zhì)交換，導(dǎo)致界面處聚合物快速析出，因而Rg迅速降低；隨后由于膜表層的形成，導(dǎo)致NMP和水傳質(zhì)速率降低，聚合物析出速率同樣降低，所以Rg減小趨勢也逐漸平緩.此外，還發(fā)現(xiàn)簡諧力場的Rg值始終大于彈簧力場，這說明彈簧力場下的PES鏈段更柔軟，而簡諧力場在其基礎(chǔ)上增加了對聚合物結(jié)構(gòu)單元間鍵角的約束，使聚合物剛性更大，這與文獻［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ū)尺寸增大表明從鑄膜液中析出的聚合物富相在逐漸增長，即相分離的程度增大.圖11（B）給出PES的相區(qū)尺寸隨時間變化.可見，隨著相分離的進行，相區(qū)尺寸先快速增長，而后緩慢增加.初始時增加較快可能是由界面處聚合物迅速成膜引起的，隨后成膜速度降低，導(dǎo)致相區(qū)尺寸增長也趨于平緩.此外，在整個過程中，簡諧力場的相區(qū)尺寸一直大于彈簧力場，這也進一步印證了簡諧力場條件下聚合物剛性較大，鏈段不容易被壓縮這一結(jié)論.

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

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

2.2 聚醚砜濃度對NIPS法成膜過程的影響

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

由圖13可以看出，同一濃度下，隨著時間的推移膜形貌的變化與圖8～圖12結(jié)果相似.當聚合物濃度較低（8%，體積分數(shù)）時，初始時聚合物的高濃度區(qū)域較少，且表層更加多孔；當聚合物濃度較高（16%，體積分數(shù)）時，由于聚合物鏈更容易發(fā)生聚集，會增加聚合物的局部濃度，從而可以觀察到更多的高濃度聚合物區(qū)域，且其表層的孔也會更少.這些觀察結(jié)果與實驗結(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表明，隨著時間的推移，峰的高度逐漸增加并向右移動，隨著聚合物初始濃度的增加，峰的寬度也緩慢增長，這說明PES濃度的增加會導(dǎo)致在界面附近形成了厚且致密的聚合物層，這與實驗觀察結(jié)果一致［33，35］.此外，隨著聚合物濃度的升高，聚合物溶液黏度變得更大，并且所形成的膜表層也會更加致密，進而會導(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隨時間的變化.在初始階段，Rg下降得很快，而后又緩慢地減小.這可能是因為初始時聚合物鏈段較為舒展，導(dǎo)致Rg的初始值較大；而隨著非溶劑的接觸，導(dǎo)致了聚合物鏈快速聚集，因而Rg快速降低最后趨于穩(wěn)定.但不同濃度PES溶液的初始Rg也是不同的，濃度為8%，12%和16%時的初始Rg值分別為5.97rc，5.84rc和5.82rc，即當聚合物濃度增加時，初始Rg值會略有降低，這與聚合物溶液理論是一致的［36，37］.但是當體系逐漸穩(wěn)定后，Rg值卻隨濃度的增大而減小了.這可能是由于低濃度時，聚合物被溶劑和非溶劑壓縮導(dǎo)致其所占空間相對較小，因而分子鏈更加團聚，所以穩(wěn)定后濃度低時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ū)尺寸隨時間的變化.由圖15（B）可見，隨著時間的推移，相區(qū)尺寸首先是快速增加，而后增長速度逐漸減緩.初始時相區(qū)尺寸快速增加可能是由界面處聚合物層的形成導(dǎo)致的，而后由于膜表層形成導(dǎo)致了物質(zhì)交換阻力變大，因而相區(qū)尺寸增長幅度減小.此外，PES濃度越高，相區(qū)尺寸越大，這是因為聚合物濃度越大越容易在NMP與水物質(zhì)交換的時候析出成膜.

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

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