劉 睿, 陳 曦

聚合反應(yīng)過程中的蒙特卡洛方法

劉 睿, 陳 曦

(浙江大學(xué) 工業(yè)控制技術(shù)國家重點(diǎn)實驗室, 控制科學(xué)與工程學(xué)院, 浙江 杭州 310027)

蒙特卡洛(Monte Carlo)方法是一種統(tǒng)計性方法，以概率模型為基礎(chǔ)，通過多次模擬實驗推定未知特性量。將蒙特卡洛方法應(yīng)用于聚合反應(yīng)過程，可以完整地得到聚合物的微觀結(jié)構(gòu)信息，且建模簡單，適用于復(fù)雜聚合反應(yīng)機(jī)制的模擬，對產(chǎn)品性能優(yōu)化、產(chǎn)率提升等有幫助作用。然而，蒙特卡洛方法在聚合反應(yīng)過程的模擬以及聚合物微觀信息的預(yù)測中計算速度慢，計算效率較低，因此，模擬的加速技術(shù)對提升該方法的有效性和實用性具有重要作用。今基于國內(nèi)外已有研究成果和課題組的前期工作，對蒙特卡洛方法在聚合反應(yīng)中的應(yīng)用進(jìn)行了詳細(xì)闡述，并圍繞聚合物微觀質(zhì)量指標(biāo)預(yù)測進(jìn)行了蒙特卡洛模擬加速算法的討論。

蒙特卡洛方法；聚合反應(yīng)過程；微觀質(zhì)量指標(biāo)；模擬加速

1 前 言

高分子材料是現(xiàn)代工業(yè)和高新技術(shù)的重要基石，已成為國民經(jīng)濟(jì)的重要基礎(chǔ)性產(chǎn)業(yè)和國家戰(zhàn)略性新型產(chǎn)業(yè)[1-2]。我國高分子材料的生產(chǎn)主要以低端材料為主，同發(fā)達(dá)國家相比，在高端材料的研發(fā)上仍有一定距離[3]。加強(qiáng)產(chǎn)品質(zhì)量控制，實現(xiàn)技術(shù)獨(dú)立自主是突破高端裝備制造和戰(zhàn)略新興產(chǎn)業(yè)發(fā)展瓶頸的關(guān)鍵。

在體現(xiàn)高分子材料性能及質(zhì)量差異的諸多因素中，微觀質(zhì)量指標(biāo)具有關(guān)鍵作用——直接決定了聚合物產(chǎn)品的最終用途和特性[4-6]。聚合物的微觀質(zhì)量指標(biāo)包括鏈長分布(chain length distribution，CLD)、相對分子質(zhì)量分布(molecular weight distribution，MWD)、化學(xué)組成分布(chemical composition distribution，CCD)、序列分布(sequence length distribution，SLD)、支鏈分布(branching distribution)等[7-8]。獲取聚合物在聚合反應(yīng)過程中各個階段的微觀信息，提供分子層面的產(chǎn)品分析，有助于提高高分子材料的產(chǎn)率與性能。然而，準(zhǔn)確并快速地預(yù)測出聚合物的微觀質(zhì)量指標(biāo)是一大難點(diǎn)。通過實驗手段得到聚合物的微觀指標(biāo)，所需時間較長[9]；構(gòu)建機(jī)理模型表征微觀過程需要準(zhǔn)確地描述出反應(yīng)機(jī)理和速率，計算復(fù)雜性高。

計算聚合物微觀質(zhì)量指標(biāo)主要分為確定性方法和統(tǒng)計性方法(也稱隨機(jī)性方法)[10-12]。求解基于物料守恒的平衡方程獲得微觀信息是典型的確定性方法[13-17]。該方法具有求解快，結(jié)果確定性高的特點(diǎn)，但難以應(yīng)用于復(fù)雜聚合過程。蒙特卡洛(Monte Carlo)方法[18-21]隸屬于統(tǒng)計性方法，無需大量的數(shù)學(xué)求解，往往是對具體聚合反應(yīng)的直接模擬，模型化方法簡單且易于理解。但是，蒙特卡洛方法的隨機(jī)性要求樣本數(shù)必須足夠多，大大延長了模擬的時間[22]。根據(jù)蒙特卡洛模擬過程是否具有時間維度，又可分為穩(wěn)態(tài)蒙特卡洛模擬和動態(tài)蒙特卡洛模擬。為快速預(yù)測聚合物微觀信息，不同的蒙特卡洛模擬加速方法被相繼提出。在此基礎(chǔ)上，本文將對蒙特卡洛方法在聚合反應(yīng)工程領(lǐng)域中的研究和應(yīng)用進(jìn)行綜述及展望。

2 蒙特卡洛方法概述

蒙特卡洛方法是一種數(shù)值計算方法，通過生成隨機(jī)數(shù)，在所有可能發(fā)生的事件集合中根據(jù)概率選擇下一個要發(fā)生的事件。隨著實驗數(shù)據(jù)量的增大，蒙特卡洛方法得到的結(jié)果也越精確，經(jīng)典的蒲豐投針實驗就是蒙特卡洛方法的一個典型問題[23]。蒙特卡洛方法已廣泛應(yīng)用于金融[24]、生物醫(yī)學(xué)[25]、計算物理[26]、項目管理[27]、材料科學(xué)[28]、分子模擬[29]等諸多領(lǐng)域。隨著計算機(jī)技術(shù)的高速發(fā)展，以蒙特卡洛方法為代表的概率模型的建立和求解變得更加高效快捷。

以求解圓周率π為例，簡要說明蒙特卡洛方法的基本原理。如圖1所示，半徑為1的圓位于邊長為2的正方形中，向正方形內(nèi)均勻隨機(jī)地撒點(diǎn)，計算落在圓內(nèi)點(diǎn)與所有點(diǎn)個數(shù)的比值并乘以4，便可得到π值，計算過程如圖2所示。撒的點(diǎn)越多，得到的π值越精確。蒙特卡洛方法計算圓周率的模擬次數(shù)與計算結(jié)果見表1，計算得到的π值為10次不同實驗下的平均值。從表1可以看出，隨著模擬量的增加，計算得到的π值更接近真實值。

圖1 蒙特卡洛方法計算圓周率

表1 蒙特卡洛方法計算圓周率：模擬次數(shù)與計算結(jié)果

與分子動力學(xué)結(jié)合，蒙特卡洛方法可以模擬聚合反應(yīng)過程并獲得聚合物的微觀信息。在上述計算圓周率的例子中，若每落下一個點(diǎn)看作一個事件，那么在聚合反應(yīng)過程的模擬中，每條反應(yīng)的發(fā)生(如鏈引發(fā)反應(yīng)、鏈增長反應(yīng)、鏈終止反應(yīng)等)可視為一個事件。蒙特卡洛方法模擬聚合反應(yīng)過程的簡要流程如圖3所示。在當(dāng)前狀態(tài)下生成隨機(jī)數(shù)，根據(jù)各個事件在當(dāng)前系統(tǒng)狀態(tài)下的概率進(jìn)行選擇后更新系統(tǒng)。模擬不斷迭代直至滿足終止條件，最終獲得聚合物的微觀質(zhì)量指標(biāo)。結(jié)果準(zhǔn)確度和模擬時間均與迭代次數(shù)正相關(guān)[30]。

圖2 蒙特卡洛方法計算圓周率的流程示意圖

圖3 蒙特卡洛方法模擬聚合反應(yīng)過程

3 聚合反應(yīng)過程中的穩(wěn)態(tài)和動態(tài)蒙特卡洛模擬

根據(jù)聚合反應(yīng)系統(tǒng)狀態(tài)是否隨著時間改變而發(fā)生變化，聚合反應(yīng)過程中的蒙特卡洛方法可分為穩(wěn)態(tài)蒙特卡洛方法與動態(tài)蒙特卡洛方法[31]。

3.1 穩(wěn)態(tài)蒙特卡洛方法

穩(wěn)態(tài)蒙特卡洛模型(steady-state Monte Carlo model)用于模擬計算聚合物微觀質(zhì)量指標(biāo)時一般不具有反應(yīng)時間維度信息，常應(yīng)用于連續(xù)反應(yīng)下的聚合物微觀結(jié)構(gòu)信息的計算[32]。

表2 三元聚合反應(yīng)的化學(xué)反應(yīng)機(jī)制

圖4 穩(wěn)態(tài)蒙特卡洛方法模擬三元聚合反應(yīng)的流程圖

蒙特卡洛方法可以給出復(fù)雜聚合反應(yīng)體系下的分子靜態(tài)行為。Beigzadeh等[33]通過穩(wěn)態(tài)蒙特卡洛方法模擬聚合物的長鏈支化并得到了聚乙烯樹脂的支化結(jié)構(gòu)。Simon等[34-35]研究了鏈行走機(jī)理下的乙烯聚合，并提出用穩(wěn)態(tài)蒙特卡洛方法得到短支鏈分布結(jié)果。Tobita[36]通過蒙特卡洛方法計算得到加入交聯(lián)和降解反應(yīng)后的聚合物分子量分布。Alshaiban等[31]應(yīng)用蒙特卡洛方法建立了一種丙烯聚合反應(yīng)過程模型，描述了聚合物的MWD以及序列分布等。

3.2 動態(tài)蒙特卡洛方法

當(dāng)系統(tǒng)狀態(tài)隨著反應(yīng)時間而改變時，動態(tài)蒙特卡洛(dynamic Monte Carlo, DMC)模型可以描述聚合反應(yīng)過程的動態(tài)變化。區(qū)別于穩(wěn)態(tài)蒙特卡洛方法，動態(tài)模擬具有反應(yīng)時間維度信息，即除了要完成穩(wěn)態(tài)模型中所必需的反應(yīng)選擇，還需要在每次事件選擇后更新反應(yīng)系統(tǒng)狀態(tài)并重新計算各個反應(yīng)事件的概率。DMC方法可以記錄任意時刻下的反應(yīng)系統(tǒng)狀態(tài)，實時追蹤反應(yīng)過程中聚合物的結(jié)構(gòu)特征，對于復(fù)雜反應(yīng)的模擬和復(fù)雜微觀質(zhì)量指標(biāo)的預(yù)測具有優(yōu)勢。對比DMC方法，基于物料守恒的平衡方程求解[37]和矩方法難以完整地表征出聚合物的微觀細(xì)節(jié)信息[11,38]。

DMC應(yīng)用于聚合反應(yīng)過程，需要進(jìn)行宏觀量到微觀量的轉(zhuǎn)換[39]：

最終可以得到

式中：0為所有反應(yīng)傾向函數(shù)之和，即

SSA通過每次模擬一個反應(yīng)事件達(dá)到模擬聚合反應(yīng)過程的目的，根據(jù)計算參數(shù)對()方法的不同，主要分為直接計算法(the direct method)、第一反應(yīng)法(the first reaction method)以及下一反應(yīng)法(the next reaction method)等。

直接計算法[40]通過在每一次迭代中單獨(dú)生成兩個[0，1]間均勻分布的隨機(jī)數(shù)rnd1和rnd2來計算系統(tǒng)時間步進(jìn)以及下一個要發(fā)生的反應(yīng)Re：

(1) 系統(tǒng)初始化：輸入反應(yīng)系統(tǒng)的初始狀態(tài)，將計數(shù)器counter和當(dāng)前反應(yīng)時間置0，并設(shè)置隨機(jī)數(shù)生成器種子；

(2) 計算傾向函數(shù)：計算當(dāng)前系統(tǒng)狀態(tài)下每條反應(yīng)的傾向函數(shù)a, (=1,…,)以及所有傾向函數(shù)之和0；

(3) 計算時間步進(jìn)并選擇反應(yīng)：生成[0,1]間均勻分布的隨機(jī)數(shù)rnd1和rnd2，根據(jù)式(7)和 (8)計算系統(tǒng)時間步進(jìn)以及選擇反應(yīng)Re；

(4) 更新系統(tǒng)狀態(tài)：更新反應(yīng)時間()；更新計數(shù)器counter(counter=counter+1)；根據(jù)所選反應(yīng)，更新系統(tǒng)狀態(tài)，即所選反應(yīng)中反應(yīng)物的個數(shù)相應(yīng)減少，生成物的個數(shù)相應(yīng)增加。

(5) 判斷迭代結(jié)束：判斷模擬是否結(jié)束；如果結(jié)束，退出，否則返回(2)；

步驟(5)中的模擬結(jié)束條件可以是反應(yīng)時間或轉(zhuǎn)化率到達(dá)指定值。步驟(2)-步驟(3)-步驟(4)可認(rèn)為是DMC在直接計算法下的一次模擬迭代。

Chung[42]將直接計算法用于研究自由基調(diào)聚反應(yīng)過程，成功表征了調(diào)節(jié)聚合中聚合物的CLD指標(biāo)。Al-Harthi等[43]通過直接計算法，建立了半間歇操作下的原子轉(zhuǎn)移自由基聚合反應(yīng)DMC模型，預(yù)測得到了梯度共聚物的平均鏈重、多分散性指數(shù)以及MWD、CLD等微觀質(zhì)量指標(biāo)。Tongtummachat等[44]用直接計算法建立了與停留時間分布結(jié)合的DMC模型，模擬連續(xù)攪拌釜下的鏈穿梭聚合反應(yīng)并得到了烯烴嵌段共聚物的微觀結(jié)構(gòu)。Agarwal等[45]結(jié)合直接計算法提出了一種蒙特卡洛集成方法，用于模擬單中心催化劑存在下的溶液共聚反應(yīng)，成功描述了聚合物的流變和結(jié)晶行為。Wang等[46]通過建立基于動力學(xué)蒙特卡洛的通用計算框架來追蹤梯度共聚物中的每個分子鏈，預(yù)測了共聚物的顯性序列。He等[47]采用直接計算法對多功能引發(fā)劑存在下的超支化共聚物進(jìn)行研究，模擬得到超支化共聚物的MWD、支化點(diǎn)比率和支化度。

表3 均聚反應(yīng)的化學(xué)反應(yīng)機(jī)制

如表3所示為均聚反應(yīng)的化學(xué)機(jī)制，表中，0為活性位點(diǎn)；Mo為單體；P和P+1分別為鏈長為和+1的活鏈；D為鏈長的死鏈。以該反應(yīng)機(jī)制為例對直接計算法應(yīng)用于聚合反應(yīng)過程模擬做更具體地說明。其中的模擬參數(shù)值列于表4，系統(tǒng)中的氫氣濃度(H2)和單體濃度(Mo)為定值。此時傾向函數(shù)為各個反應(yīng)的微觀反應(yīng)速率：

表4 均聚反應(yīng)機(jī)制的初始參數(shù)值

時間步進(jìn)為

對表3中的反應(yīng)機(jī)制設(shè)置不同的反應(yīng)時間和反應(yīng)體積，進(jìn)行2次模擬，得到的聚合物CLD如圖5所示，其中圖5(a)表示反應(yīng)時間為0.5 s，初始的P0數(shù)量為106的模擬結(jié)果，計算模擬所需時長約34 min；圖5(b)表示反應(yīng)時間為30 min，初始的P0數(shù)量為104的模擬結(jié)果，計算模擬所需時長約17 h。DMC模擬耗時與反應(yīng)時間正相關(guān)。模擬的分子數(shù)量越多，DMC方法得到的結(jié)果越精確，但計算模擬時間也相應(yīng)延長。

區(qū)別于直接計算法，第一反應(yīng)法[48]在計算參數(shù)對()上有所不同。第一反應(yīng)法通過生成與反應(yīng)個數(shù)相同數(shù)量的[0，1]間均勻分布的隨機(jī)數(shù)來計算時間步進(jìn)

其中具有最小值的反應(yīng)Re為被選定的反應(yīng)事件。第一反應(yīng)法通過生成多個隨機(jī)數(shù)來確定參數(shù)對()的值，其算法可描述如下：

(1) 系統(tǒng)初始化：確定反應(yīng)系統(tǒng)的初始狀態(tài)，將計數(shù)器counter和當(dāng)前反應(yīng)時間置0，初始化隨機(jī)數(shù)生成器；

(2) 計算傾向函數(shù)：計算每條反應(yīng)的傾向函數(shù)a, (=1,…,)以及所有傾向函數(shù)之和0；

(3) 計算時間步進(jìn)并選擇反應(yīng)：生成個[0,1]間均勻分布的隨機(jī)數(shù)，計算τ, (=1)。根據(jù)式(13)和 (14)計算，并選擇反應(yīng)Re；

(4) 更新系統(tǒng)狀態(tài)：更新反應(yīng)時間()；更新計數(shù)器counter(counter=counter+1)；根據(jù)上步所選反應(yīng)，更新系統(tǒng)狀態(tài)；

(5) 判斷迭代結(jié)束：判斷模擬是否結(jié)束，如果結(jié)束，退出，否則返回(2)。

第一反應(yīng)法與直接計算法得到的模擬結(jié)果一致，但由于第一反應(yīng)法每步迭代中浪費(fèi)(-1)個隨機(jī)數(shù)，模擬效率較低。在第一反應(yīng)法的基礎(chǔ)上，Gibson等[49]改進(jìn)了數(shù)據(jù)處理策略并提出了下一反應(yīng)法，在每一步的反應(yīng)事件選擇和系統(tǒng)狀態(tài)更新中僅需要生成一個隨機(jī)數(shù)，模擬效率明顯提升。Rubenstein等[50]基于直接計算法和下一反應(yīng)法獲得了成核聚合反應(yīng)過程(nucleated polymerization)中的聚合物動態(tài)特性。

值得注意的是，下一反應(yīng)法在提高模擬效率的同時對計算機(jī)存儲，尤其是CPU存儲能力有較高要求。為進(jìn)一步提高計算效率，Cao等[51]提出了優(yōu)化直接法(optimized direct method)；McCollum等[52]提出了分類直接法(sorting direct method)。將這一系列隨機(jī)算法與聚合反應(yīng)過程結(jié)合，模擬預(yù)測聚合物性質(zhì)時，應(yīng)針對具體的反應(yīng)機(jī)制、反應(yīng)參數(shù)以及應(yīng)用需求等，權(quán)衡數(shù)據(jù)存儲與數(shù)據(jù)計算時間，選擇合適的方法用于反應(yīng)模擬[53]。

4 蒙特卡洛模擬加速方法

當(dāng)模擬復(fù)雜系統(tǒng)時，每步迭代僅模擬一個反應(yīng)事件的方法效率較低。即便很多學(xué)者提出了模擬效率高于直接計算法的改進(jìn)SSA方法，但受限于“一步一個反應(yīng)”，依舊無法完成快速的聚合反應(yīng)過程模擬。蒙特卡洛模擬的計算時間和結(jié)果精度嚴(yán)格依賴于系統(tǒng)中的分子數(shù)以及發(fā)生的反應(yīng)事件次數(shù)，Al-Harthi等[54]研究了動態(tài)蒙特卡洛方法模擬自由基聚合反應(yīng)過程，并比較了不同反應(yīng)體積下得到的MWD結(jié)果。當(dāng)模擬體積過小時，MWD結(jié)果曲線有明顯的鋸齒狀，隨機(jī)波動嚴(yán)重；當(dāng)模擬體積增加1 000倍后，得到的模擬結(jié)果精確的同時模擬時間增加17 308倍。很多學(xué)者提出了基于計算機(jī)算力提升和以犧牲結(jié)果精確度為代價換取模擬效率的加速方法，主要包括跳躍方法(-leaping method)及其改進(jìn)算法、縮放方法(scaling method)、混合方法(hybrid method)、并行方法等。

4.1 τ跳躍方法

傳統(tǒng)的SSA在時間步進(jìn)內(nèi)，有且僅有一個反應(yīng)發(fā)生。這種方法的優(yōu)點(diǎn)在于可以完整地記錄系統(tǒng)在時間維度上的所有狀態(tài)信息，但大量依次發(fā)生的反應(yīng)事件導(dǎo)致計算時間過長?？紤]到時間維度上完整的系統(tǒng)狀態(tài)信息并非完全必要，Gillespie[55]提出了跳躍方法，允許在時間間隔內(nèi)發(fā)生一定數(shù)量的多個反應(yīng)。如圖6所示，跳躍方法將時間維度劃分為多個子時間間隔[，+leap)，計算在各個子時間間隔內(nèi)反應(yīng)事件的發(fā)生次數(shù)以及事件類型。在模擬過程中，跳躍方法只需確定這些子時間間隔內(nèi)的系統(tǒng)狀態(tài)變化，這種“一步多個反應(yīng)”的方法，計算效率較高。

圖6 傳統(tǒng)隨機(jī)模擬算法與τ跳躍方法對比

對跳躍方法進(jìn)行數(shù)學(xué)描述，首先定義概率函數(shù)：

式中：1,2,…,K分別為第1,2,…,個反應(yīng)在時間間隔內(nèi)的發(fā)生次數(shù)。對于(1,2,…,K)這個整數(shù)隨機(jī)變量，定義：

為確定(1,2,…,K|leap;,)中的leap值，反應(yīng)系統(tǒng)需要滿足跳躍條件(leap condition)：leap(leap>0)值不能讓時間間隔[,+leap)內(nèi)傾向函數(shù)a()的變化過于明顯。當(dāng)滿足跳躍條件后，式(16)可表示為

leap值的確定對跳躍方法的發(fā)展起到關(guān)鍵作用，而以往靠經(jīng)驗選擇的leap值無法高效模擬反應(yīng)過程并得到系統(tǒng)狀態(tài)及分子微觀結(jié)構(gòu)。Gillespie等[56]對跳躍方法進(jìn)行擴(kuò)展，提出leap值的選擇策略。此外，由于各個反應(yīng)在反應(yīng)系統(tǒng)中的發(fā)生頻率可能會有很大的差別，出現(xiàn)頻率少的反應(yīng)在跳躍方法中很難被采樣到，從而導(dǎo)致最終結(jié)果發(fā)生偏差，這類偏差問題稱為“系統(tǒng)剛性問題”。Rathinam等[57]針對跳躍方法導(dǎo)致的系統(tǒng)剛性問題提出了隱性leap值計算方法。Cao等[58-59]在已有工作的基礎(chǔ)上對泊松跳躍方法中的負(fù)值問題進(jìn)行修正，并提出了最大leap值的選擇方法以加速模擬。Brand?o等[41]分別用直接計算法和跳躍方法模擬伴隨長支鏈形成的共聚反應(yīng)過程，同樣的反應(yīng)時間和系統(tǒng)體積下，跳躍方法的計算效率明顯高于直接計算法。

4.2 縮放方法(scaling method)

將蒙特卡洛方法應(yīng)用于模擬聚合反應(yīng)過程并計算聚合物的分子量分布、記錄各個分子鏈的單體序列等微觀結(jié)構(gòu)信息時，為了能準(zhǔn)確地表征出反應(yīng)過程，至少需要一定數(shù)量的分子數(shù)?？s放方法可以有效縮放模擬的分子數(shù)量和反應(yīng)速率常數(shù)，從而減少模擬耗時。Gao等[60]推導(dǎo)了共聚反應(yīng)的顯性表達(dá)式計算模擬所需的分子數(shù)下限以及縮放因數(shù)，通過對丙烯酸酯共聚機(jī)制進(jìn)行模擬得到了聚合物的MWD以及序列分布，模擬效率提升明顯。Lin等[61]提出了部分縮放法(partial scaling)，可以對反應(yīng)系統(tǒng)進(jìn)行動態(tài)縮放：考慮系統(tǒng)中的每一個單獨(dú)反應(yīng)，若放縮尺度小于提前設(shè)置的閾值，則對系統(tǒng)進(jìn)行縮放。該方法為每一個反應(yīng)計算出各自的縮放因數(shù)，對“剛性系統(tǒng)”的模擬準(zhǔn)確率和效率提升明顯。

4.3 并行策略

近幾年來計算集群、工作站的搭建，給并行算法提供了算力支持。將并行策略應(yīng)用于蒙特卡洛方法，為加速預(yù)測聚合物微觀質(zhì)量指標(biāo)提供了新思路。蒙特卡洛方法應(yīng)用于聚合反應(yīng)過程，模擬結(jié)果依賴于系統(tǒng)分子量。并行計算策略通過將各個模擬或迭代任務(wù)分布在不同的計算節(jié)點(diǎn)上，實現(xiàn)了加速模擬。Chaffey-Millar等[62]通過分布式策略以及CPU多核多線程并行策略，將計算數(shù)據(jù)分配到不同的計算機(jī)節(jié)點(diǎn)上并行處理，加速了聚合物微觀質(zhì)量指標(biāo)的模擬。Li等[63]提出了一種離散隨機(jī)模擬算法來預(yù)測無法通過確定性算法獲取到的分子動態(tài)行為，并結(jié)合基于圖像處理單元(graphics processing unit platform，GPU)的并行策略，得到了較好的加速效果。Weng等[64]在GPU上進(jìn)行多核并行編程，完成了穩(wěn)態(tài)蒙特卡洛預(yù)測共聚物微觀質(zhì)量指標(biāo)的加速。該方法結(jié)合GPU的高性能并行運(yùn)算，將每條分子鏈的模擬任務(wù)分配到不同的GPU線程中，實現(xiàn)了和GPU核數(shù)相同數(shù)量級的大規(guī)模并行，準(zhǔn)確得到了聚合物微觀信息的同時，模擬加速效果明顯。Ma等[65]結(jié)合GPU并行技術(shù)和蒙特卡洛方法預(yù)測了聚合物的CCD，并提出由自適應(yīng)模擬算法與連續(xù)邊界收縮算法組成的內(nèi)嵌蒙特卡洛模型的系統(tǒng)化優(yōu)化方法，該方法加快了優(yōu)化過程中的模擬計算。此外，Lemarchand等[66]將并行算法與分子動力學(xué)結(jié)合，對生成的長聚合物鏈進(jìn)行模擬；De Buyl等[67]將并行算法應(yīng)用于逐步聚合反應(yīng)過程與鏈增長聚合反應(yīng)過程的模擬；Zhang等[68]建立了一個并行多尺度模擬框架，克服了傳統(tǒng)方法模擬慢的問題，并成功應(yīng)用于具有超支化的復(fù)雜聚合反應(yīng)；Liu等[69]結(jié)合不同的聚合反應(yīng)機(jī)制和硬件架構(gòu)，提出了基于不同硬件平臺的動態(tài)蒙特卡洛模擬加速算法。

4.4 混合方法(hybrid method)

混合方法結(jié)合了確定性模型和統(tǒng)計性模型的優(yōu)點(diǎn)，通過確定性方法求解聚合物基本信息，再結(jié)合蒙特卡洛方法模擬得到具體的聚合物微觀性質(zhì)?；旌戏椒ㄔ讷@得聚合物微觀質(zhì)量指標(biāo)的同時模擬量大大減少。Schuette等[70]提出了混合伽遼金-蒙特卡洛算法(hybrid Galerkin-Monte Carlo method)，將伽遼金法確定性求解和蒙特卡洛統(tǒng)計性方法同時應(yīng)用于模擬計算分子鏈長分布；He等[71]用混合方法模擬了活性自由基聚合中的CLD，在交換反應(yīng)平衡的假設(shè)下，用解析表達(dá)表征活化及失活交換反應(yīng)，用蒙特卡洛方法模擬鏈增長反應(yīng)和不可逆的鏈終止反應(yīng)。Tripathi等[72]提出了一種應(yīng)用在自由基聚合反應(yīng)中的高效率混合算法，明顯縮短了聚合物微觀結(jié)構(gòu)信息的計算時間，同時還能保證結(jié)果的準(zhǔn)確性與完整性。在該方法中，原來在動態(tài)蒙特卡洛方法中占據(jù)大量模擬時間的鏈增長反應(yīng)由微分方程處理，而對于鏈終止反應(yīng)和鏈引發(fā)反應(yīng)則使用統(tǒng)計性方法處理。Ozcam等[73]在已有研究的基礎(chǔ)上，提出了新的混合方法——逐鏈蒙特卡洛方法(chain-by-chain Monte Carlo method)，先用矩方法獲得聚合物的整體信息，再將其與蒙特卡洛方法結(jié)合，得到分子鏈的細(xì)節(jié)信息。

5 結(jié) 論

以蒙特卡洛方法為代表的統(tǒng)計性模型(隨機(jī)性模型)具有建模簡單、適用于復(fù)雜聚合反應(yīng)過程模擬等優(yōu)點(diǎn)，在化學(xué)工程以及高分子領(lǐng)域的應(yīng)用愈發(fā)廣泛。針對聚合物微觀質(zhì)量指標(biāo)的預(yù)測，蒙特卡洛方法可以在時間及空間維度下得到聚合物分子內(nèi)和分子間完整的微觀結(jié)構(gòu)信息，對高端高分子材料產(chǎn)品的生產(chǎn)及質(zhì)量提升有著深刻影響。

蒙特卡洛方法應(yīng)用在以產(chǎn)品結(jié)構(gòu)優(yōu)化為技術(shù)導(dǎo)向的高端材料研究領(lǐng)域意義重大，聚合反應(yīng)過程中的蒙特卡洛方法在精細(xì)化建模、基于統(tǒng)計性模型的加速優(yōu)化、模擬過程中的數(shù)據(jù)存儲優(yōu)化、復(fù)雜聚合反應(yīng)過程模擬、多樣性微觀質(zhì)量指標(biāo)的預(yù)測等方向有很大的探索空間。另外，隨著機(jī)器學(xué)習(xí)的快速發(fā)展，蒙特卡洛方法還有望與強(qiáng)化學(xué)習(xí)、深度學(xué)習(xí)等技術(shù)整合，優(yōu)化生產(chǎn)流程以及聚合產(chǎn)品性能。

[1] 喬金樑. 我國高分子材料產(chǎn)業(yè)轉(zhuǎn)型發(fā)展的思考 [J]. 石油化工, 2015, 44(9): 1033-1037.

QIAO J L. Further development of polymer industry in china [J]. Petrochemical Technology, 2015, 44(9): 1033-1037.

[2] 崔毅杰, 宋盛菊, 屈小中, 等. 中國制造2025視野下高分子材料產(chǎn)業(yè)發(fā)展的再思考與展望 [J]. 工程研究-跨學(xué)科視野中的工程, 2017, 9(6): 568-576.

CUI Y J, SONG S J, QU X Z,. Rethinking and outlook on the development of polymeric materials industry under the view of “made in China 2025” [J]. Journal of Engineering Studies, 2017, 9(6): 568-576.

[3] 王瑜鑫. 我國高分子材料產(chǎn)業(yè)發(fā)展現(xiàn)狀與前景——以塑料橡膠為例 [J]. 環(huán)渤海經(jīng)濟(jì)瞭望, 2018, 286(7): 68-68.

WANG Y X. The development status and prospects of polymer materials industry in China--taking plastics and rubber as an example [J]. Economic Outlook the Bohai Sea, 2018, 286(7): 68-68.

[4] 張晨. 基于分子量分布的聚合過程建模、優(yōu)化與流程重構(gòu) [D]. 杭州: 浙江大學(xué), 2016.

ZHANG C. Modeling, optimization, and flowsheet reconfiguration of polymerization process with embedded molecular weight distribution [D]. Hangzhou: Zhejiang University, 2016

[5] 翁金祖. 基于微觀質(zhì)量指標(biāo)的聚合過程模擬與優(yōu)化方法研究 [D]. 杭州: 浙江大學(xué), 2016.

WENG J Z. Microstructure-based simulation and optimization methods for polymerization processes [D]. Hangzhou: Zhejiang University, 2016

[6] 馬彥楠. 面向微觀結(jié)構(gòu)質(zhì)量優(yōu)化的模型重構(gòu)與求解加速方法研究 [D]. 杭州: 浙江大學(xué), 2019.

MA Y N. Model reconstruction and optimization acceleration with microstructural quality indices [D]. Hangzhou: Zhejiang University, 2019.

[7] RIAHINEZHAD M, MCMANUS N, PENLIDIS A. Effect of monomer concentration and pH on reaction kinetics and copolymer microstructure of acrylamide/acrylic acid copolymer [J]. Macromolecular Reaction Engineering, 2015, 9(2): 100-113.

[8] LUTZ J F, JAHED N, MATYJASZEWSKI K. Preparation and characterization of graft terpolymers with controlled molecular structure [J]. Journal of Polymer Science Part A-Polymer Chemistry, 2004, 42(8): 1939-1952.

[9] SHAN C, SOARES J, PENLIDIS A. Ethylene/1-octene copolymerization studies with in situ supported metallocene catalysts: Effect of polymerization parameters on the catalyst activity and polymer microstructure [J]. Journal of Polymer Science Part A-Polymer Chemistry, 2002, 40(24): 4426-4451.

[10] 童偉達(dá), 陸建明, 孫猛, 等. 自由基聚合反應(yīng)的Monte Carlo模擬方法 [J]. 高等學(xué)?；瘜W(xué)學(xué)報, 1991, 12(7): 979-982.

TONG W D, LU J M, SUN M,. Monte Carlo simulation method of radical polymerization [J]. Chemical Journal of Chinese University, 1991, 12(7): 979-982.

[11] MASTAN E, ZHU S. Method of moments: A versatile tool for deterministic modeling of polymerization kinetics [J]. European Polymer Journal, 2015, 68:139-160.

[12] CHEN X, SHAO Z, GU X,. Process intensification of polymerization processes with embedded molecular weight distributions models: An advanced optimization approach [J]. Industrial & Engineering Chemistry Research, 2019, 58(15SI): 6133-6145.

[13] ZHANG C, SHAO Z, CHEN X,. Simulation and optimization of polymer molecular weight distribution with nonideal reactors [J]. Computers & Chemical Engineering, 2017, 106: 744-757.

[14] KANG J, SHAO Z, CHEN X,. Fast and reliable computational strategy for developing a rigorous model-driven soft sensor of dynamic molecular weight distribution [J]. Journal of Process Control, 2017, 56: 79-99.

[15] ZHANG C, SHAO Z, CHEN X,. Optimal flowsheet configuration of a polymerization process with embedded molecular weight distributions [J]. AIChE Journal, 2016, 62(1): 131-145.

[16] WENG J, CHEN X, BIEGLER L T. A multi-thread parallel computation method for dynamic simulation of molecular weight distribution of multisite polymerization [J]. Computers & Chemical Engineering, 2015, 82(2): 55-67.

[17] ZHANG C, SHAO Z, CHEN X,. Kinetic parameter estimation of HDPE slurry process from molecular weight distribution: estimability analysis and multistep methodology [J]. AIChE Journal, 2014, 60(10): 3442-3459.

[18] 陸建明, 孫猛, 楊玉良. 不等活性縮聚反應(yīng)動力學(xué)及其分子量分布Monte Carlo模擬 [J]. 高分子學(xué)報, 1994, 1(5): 565-572.

LU J M, SUN M, YANG Y L. Monte Carlo studies on the condensation polymerization of unequal functional group reactivity [J]. Acta Polymerica Sinica, 1994, 1(5): 565-572.

[19] 楊玉良, 張紅東. 高分子科學(xué)中的Monte Carlo方法 [M]. 上海: 復(fù)旦大學(xué)出版社, 1993.

YANG Y L, ZHANG H D. Monte Carlo method in polymer science [M]. Shanghai: Fudan University Press, 1993.

[20] LYU J, GAO Y, ZHANG Z,. Monte Carlo simulations of atom transfer radical (homo)polymerization of divinyl monomers: applicability of Flory-Stockmayer theory [J]. Macromolecules, 2018, 51(17): 6673-6681.

[21] ZHAO Y R, MCAULEY K B, IEDEMA P D,. Advanced Monte Carlo modeling using weight-based selection of arborescent polyisobutylene molecules in a batch reactor [J]. Macromolecular Theory and Simulations, 2016, 25(2): 134-154.

[22] 張勝, 駱德莘, 楊玉良, 等. 三元共聚反應(yīng)序列分布的Monte Carlo模擬 [J]. 功能高分子學(xué)報, 1993, 6(1): 41-45.

ZHANG S, LUO D X, YANG Y L. Monte Carlo simulation of sequence distribution of terpolymerization [J]. Functional Polymer, 1993, 6(1): 41-45.

[23] 黃朝霞. 蒲豐投針問題研究 [J]. 集美大學(xué)學(xué)報(自然科學(xué)版), 2005, 10(4): 381-384.

HUANG Z X. Study on Buffon throwing-pin problem [J]. Journal of Jimei University (Natural Science), 2005, 10(4): 381-384.

[24] AVRAMIDIS A N, L'ECUYER P. Efficient Monte Carlo and quasi-Monte Carlo option pricing under the variance gamma model [J]. Management Science, 2006, 52(12): 1930-1944.

[25] WANG L H, JACQUES S L, ZHENG L Q. MCML-Monte-Carlo modeling of light transport in multilayered tissues [J]. Computer Methods and Programs in Biomedicine, 1995, 47(2): 131-146.

[26] 司瑞芳, 杜永清. 蒙特卡洛方法在高能核物理研究中的一些應(yīng)用 [J]. 太原師范學(xué)院學(xué)報(自然科學(xué)版), 2010, 9(2): 84-86.

SI R F, DU Y Q. Some applications of Monte Carlo method in investigations of high energy nuclear physics [J]. Journal of Taiyuan Normal University (Natural Science Edition), 2010, 9(2): 84-86.

[27] VAN DORP J R, DUFFEY M R. Statistical dependence in risk analysis for project networks using Monte Carlo methods [J]. International Journal of Production Economics, 1999, 58(1): 17-29.

[28] SERGEI P, KRUZ K, YOUSRA N,. Large scale hybrid Monte Carlo simulations for structure and property prediction [J]. NPJ Computational Materials, 2018, 4: 80.

[29] MEIMAROGLOU D, KIPARISSIDES C. Review of Monte Carlo methods for the prediction of distributed molecular and morphological polymer properties [J]. Industrial & Engineering Chemistry Research, 2014, 53(22): 8963-8979.

[30] REGO A S C, BRANDAO A L T. General method for Speeding up kinetic Monte Carlo simulations [J]. Industrial & Engineering Chemistry Research, 2020, 59(19): 9034-9042.

[31] ALSHAIBAN A, SOARES J B P. Mathematical modeling of the microstructure of poly(propylene) made with Ziegler-Natta catalysts in the presence of electron donors [J]. Macromolecular Reaction Engineering, 2011, 5(2): 96-116.

[32] BRAND?O A L T, SOARES J B P, PINTO J C,. When polymer reaction engineers play dice: Applications of Monte Carlo models in PRE [J]. Macromolecular Reaction Engineering, 2015, 9(3): 141-185.

[33] BEIGZADEH D, SOARES J B P, DUEVER T A,. Analysis of branching structure in polyethylene resins synthesized with constrained-geometry catalyst systems, using monte carlo simulation [J]. Polymer Reaction Engineering, 1999, 7(2):195-205.

[34] SIMON L C, SOARES J, DE SOUZA R F. Monte-Carlo simulation of branching distribution in Ni-diimine catalyzed polyethylene [J]. AIChE Journal, 2000, 46(6): 1234-1240.

[35] SIMON L C, WILLIAMS C P, SOARES J,. Effect of polymerization temperature and pressure on the microstructure of Ni-diimine-catalyzed polyethylene: parameter identification for Monte-Carlo simulation [J]. Chemical Engineering Science, 2001, 56(13): 4181-4190.

[36] TOBITA H. Simulation model for the modification of polymers via crosslinking and degradation [J]. Polymer (Guilford), 1995, 36(13): 2585-2596.

[37] RIVERO P, HERRERA R. Modeling the kinetics of anionic polymerization in cyclohexane as a non-complexing solvent [J]. Journal of Polymer Research, 2011, 18(4): 519-526.

[38] AL-HARTHI M, SOARES J B P, SIMON L C. Dynamic Monte Carlo simulation of ATRP with bifunctional initiators [J]. Macromolecular Reaction Engineering, 2007, 1(1): 95-105.

[39] HE J, ZHANG H, YANG Y. Monte Carlo simulation of chain length distribution in radical polymerization with transfer reaction [J]. Macromolecular Theory and Simulations, 1995, 4(4): 811-819.

[40] GILLESPIE D T. Exact stochastic simulation of coupled chemical-reactions [J]. Journal of Physical Chemistry, 1977, 81(25): 2340-2361.

[41] BRAND?O A L T, SOARES J B P, PINTO J C,. Comparison of different dynamic Monte Carlo methods for the simulation of olefin polymerization [J]. Macromolecular Symposia. 2016, 360(1): 160-178.

[42] CHUNG I. Monte Carlo simulation of free radical telomerization [J]. Polymer, 2000, 41(15): 5643-5651.

[43] AL-HARTHI M, KHAN M J, ABBASI S H,. Gradient copolymers by ATRP in semibatch reactors: Dynamic Monte Carlo simulation [J]. Macromolecular Reaction Engineering, 2009, 3(4): 148-159.

[44] TONGTUMMACHAT T, ANANTAWARASKUL S, SOARES J B P. Dynamic Monte Carlo Simulation of olefin block copolymers (OBCs) produced via chain-shuttling polymerization: Effect of kinetic rate constants on chain microstructure [J]. Macromolecular Reaction Engineering, 2018, 12(4): 1800021.1-1800021.9.

[45] AGARWAL U, DEN DOELDER J. Integrated kinetic Monte Carlo–structure–rheology model for solution copolymerization of ethylene and α-olefins [J]. Macromolecules, 2018, 52(1): 292-310.

[46] WANG L, BROADBELT L J. Explicit sequence of styrene/methyl methacrylate gradient copolymers synthesized by forced gradient copolymerization with nitroxide-mediated controlled radical polymerization [J]. Macromolecules, 2009, 42(20): 7961-7968.

[47] HE X, LIANG H, PAN C. Monte Carlo simulation of hyperbranched copolymerizations in the presence of a multifunctional initiator [J]. Macromolecular Theory and Simulations, 2001, 10(3): 196-203.

[48] GILLESPIE D T. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions [J]. Journal of Computational Physics, 1976, 22(4): 403-434.

[49] GIBSON M A, BRUCK J. Efficient exact stochastic simulation of chemical systems with many species and many channels [J]. The Journal of Physical Chemistry A, 2000, 104(9): 1876-1889.

[50] RUBENSTEIN R, GRAY P C, CLELAND T J,. Dynamics of the nucleated polymerization model of prion replication [J]. Biophysical Chemistry, 2007, 125(2/3): 360-367.

[51] CAO Y, LI H, PETZOLD L. Efficient formulation of the stochastic simulation algorithm for chemically reacting systems [J]. Journal of Chemical Physics, 2004, 121(9): 4059-4067.

[52] MCCOLLUM J M, PETERSON G D, COX C D,. The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior [J]. Computational Biology and Chemistry, 2006, 30(1): 39-49.

[53] LECCA P. Stochastic chemical kinetics:A review of the modelling and simulation approaches [J]. Biophysical Reviews, 2013, 5(4): 323-345.

[54] AL-HARTHI M, SOARES J B P, SIMON L C. Dynamic Monte Carlo simulation of atom-transfer radical polymerization [J]. Macromolecular Materials and Engineering, 2006, 291(8): 993-1003.

[55] GILLESPIE D T. Approximate accelerated stochastic simulation of chemically reacting systems [J]. Journal of Chemical Physics, 2001, 115(4): 1716-1733.

[56] GILLESPIE D T, PETZOLD L R. Improved leap-size selection for accelerated stochastic simulation [J]. Journal of Chemical Physics, 2003, 119(16): 8229-8234.

[57] RATHINAM M, PETZOLD L R, CAO Y,. Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method [J]. Journal of Chemical Physics, 2003, 119(24): 12784-12794.

[58] CAO Y, GILLESPIE D T, PETZOLD L R. Avoiding negative populations in explicit Poisson tau-leaping [J]. Journal of Chemical Physics, 2005, 123(5): 54104.

[59] CAO Y, GILLESPIE D T, PETZOLD L R. Efficient step size selection for the tau-leaping simulation method [J]. Journal of Chemical Physics, 2006, 124(4): 044109.1-044109.11.

[60] GAO H, OAKLEY L H, KONSTANTINOV I A,. Acceleration of kinetic Monte Carlo method for the simulation of free radical copolymerization through scaling [J]. Industrial & Engineering Chemistry Research, 2015, 54(48): 11975-11985.

[61] LIN Y T, FENG S, HLAVACEK W S. Scaling methods for accelerating kinetic Monte Carlo simulations of chemical reaction networks [J]. Journal of Chemical Physics, 2019, 150(24): 244101.1-244101.17..

[62] CHAFFEY-MILLAR H, STEWART D, CHAKRAVARTY M M T,. A parallelised high performance Monte Carlo simulation approach for complex polymerisation kinetics [J]. Macromolecular Theory and Simulations, 2007, 16(6): 575-592.

[63] LI H, PETZOLD L. Efficient parallelization of the stochastic simulation algorithm for chemically reacting systems on the graphics processing unit [J]. International Journal of High Performance Computing Applications, 2010, 24(2): 107-116.

[64] WENG J, CHEN X, YAO Z,. Parallel Monte Carlo simulation of molecular weight distribution and chemical composition distribution for copolymerization on a graphics processing unit platform [J]. Macromolecular Theory and Simulations, 2015, 24(5): 521-536.

[65] MA Y, CHEN X, BIEGLER L T. Monte-Carlo-simulation-based optimization for copolymerization processes with embedded chemical composition distribution [J]. Computers & Chemical Engineering, 2018, 109: 261-275.

[66] LEMARCHAND C A, BOUSQUET D, SCHNELL B,. A parallel algorithm to produce long polymer chains in molecular dynamics [J]. Journal of Chemical Physics, 2019, 150(22): 224902.1-224902.16.

[67] DE BUYL P, NIES E. A parallel algorithm for step- and chain-growth polymerization in molecular dynamics [J]. Journal of Chemical Physics, 2015, 142(13): 134102.

[68] ZHANG Z, KRAJNIAK J, SAMAEY G,. A parallel multiscale simulation framework for complex polymerization: ab2‐type monomer hyperbranched polymerization as an example [J]. Advanced Theory and Simulations, 2018, 2(2): 1800102.

[69] LIU R, ARMAOU A, CHEN X. Adaptable parallel acceleration strategy for dynamic Monte Carlo simulations of polymerization with microscopic resolution [J]. Industrial & Engineering Chemistry Research, 2021, 60 (17): 6173-6187.

[70] SCHUETTE C, WULKOW M. A hybrid Galerkin-Monte-Carlo approach to higher-dimensional population balances in polymerization kinetics [J]. Macromolecular Reaction Engineering, 2010, 4(9/10): 562-577.

[71] HE J P, ZHANG H D, CHEN J M,. Monte Carlo simulation of kinetics and chain length distributions in living free-radical polymerization [J]. Macromolecules, 1997, 30(25): 8010-8018.

[72] TRIPATHI A K, SUNDBERG D C. A hybrid algorithm for accurate and efficient Monte Carlo simulations of free-radical polymerization reactions [J]. Macromolecular Theory and Simulations, 2015, 24(1): 52-64.

[73] OZCAM D D, TEYMOUR F. Chain-by-chain Monte Carlo simulation: A novel hybrid method for modeling polymerization. Part I. Linear controlled radical polymerization systems [J]. Macromolecular Reaction Engineering, 2017, 11(1): 1600042.1-160042.22.

Review on the Monte Carlo method in polymerization processes

LIU Rui, CHEN Xi

(State Key Laboratory of Industrial Control Technology, College of Control Science and Engineering,Zhejiang University, Hangzhou 310027, China)

The Monte Carlo method is a statistical method that generates random numbers to obtain the final result according to the probabilities of various events through multiple simulations. With the Monte Carlo method, it is possible to simulate a polymerization process with complex kinetics and predict microscopic properties of a polymer. However, the simulation of polymerization process and prediction of polymer microscopic information by the Monte Carlo method are usually time-consuming. Acceleration of the Monte Carlo simulations for polymerization processes is of great importance. This work is to review the Monte Carlo method and its applications in polymerization processes. The acceleration algorithms of Monte Carlo methods for predicting polymer microscopic qualities are discussed.

Monte Carlo method; polymerization; microscopic qualities; acceleration

TQ 316.3

A

10.3969/j.issn.1003-9015.2021.03.001

1003-9015(2021)03-0389-11

2020-12-17；

2021-02-08。

國家重點(diǎn)研發(fā)計劃(2017YFE0106700)。

劉睿(1997-)，女，江蘇徐州人，浙江大學(xué)博士生。

陳曦，E-mail：xi_chen@zju.edu.cn