阮春蕾，劉春太

（1河南科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，河南 洛陽(yáng) 471023；2橡塑模具國(guó)家工程研究中心，河南 鄭州 450001）

?

剪切流場(chǎng)中聚乙烯結(jié)晶過(guò)程的建模與模擬

阮春蕾1，2，劉春太2

（1河南科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，河南 洛陽(yáng) 471023；2橡塑模具國(guó)家工程研究中心，河南 鄭州 450001）

基于Eder模型推導(dǎo)了剪切流場(chǎng)中球晶、串晶形態(tài)演化的數(shù)學(xué)模型，將第一法向應(yīng)力差作為串晶成核的驅(qū)動(dòng)，并引入兩相懸浮模型描述體系，認(rèn)為其由無(wú)定形相和半結(jié)晶相組成，分別用FENE-P模型和剛性啞鈴模型描述?；谏鲜鰯?shù)學(xué)模型，分別構(gòu)造了捕捉球晶、串晶生長(zhǎng)的Monte Carlo法與體系控制方程求解的有限差分法，成功模擬了二維剪切流場(chǎng)中聚乙烯的結(jié)晶過(guò)程，給出了球晶、串晶的形態(tài)演化，分析了剪切時(shí)間、剪切速率對(duì)串晶成核密度、結(jié)晶速率、流體黏度等的影響。數(shù)值結(jié)果表明：所構(gòu)造的Monte Carlo法合理有效，不僅成功捕捉了晶體的生長(zhǎng)與碰撞，而且較為準(zhǔn)確地預(yù)測(cè)了結(jié)晶速率。此外，提高剪切時(shí)間或剪切速率，將增加串晶成核密度、提高結(jié)晶速率、使流體黏度突增的時(shí)間點(diǎn)提前。

流動(dòng)誘導(dǎo)結(jié)晶；形態(tài)學(xué)；動(dòng)力學(xué)；Monte Carlo模擬；結(jié)晶

DOI：10.11949/j.issn.0438-1157.20151472

引　言

聚合物的結(jié)晶過(guò)程及其形態(tài)演化是決定制品力學(xué)性能的關(guān)鍵因素［1］。廣義而言，聚合物必須經(jīng)過(guò)成型加工才能成為有用的制品。成型過(guò)程中，聚合物受到流場(chǎng)的作用，高分子鏈發(fā)生運(yùn)動(dòng)，在低溫條件下形成凝聚態(tài)結(jié)構(gòu)，并最終形成球晶、串晶等復(fù)雜結(jié)晶形態(tài)［1］。由于成型中以剪切流場(chǎng)居多，因此，本文將研究限定在簡(jiǎn)單剪切流場(chǎng)下，并以聚乙烯為代表，研究結(jié)晶形態(tài)與動(dòng)力學(xué)，流體性質(zhì)等的改變，揭示外場(chǎng)因素與結(jié)晶過(guò)程間的相關(guān)規(guī)律，并為復(fù)雜成型條件下的結(jié)晶過(guò)程提供指導(dǎo)。

流場(chǎng)中的聚合物結(jié)晶又被稱為流動(dòng)誘導(dǎo)結(jié)晶（flow induced crystallization，F(xiàn)IC）［1-2］。流動(dòng)誘導(dǎo)結(jié)晶的實(shí)驗(yàn)研究表明［1］：與靜態(tài)結(jié)晶相比，流場(chǎng)中的結(jié)晶發(fā)生了實(shí)質(zhì)性變化，不僅提高了結(jié)晶速率，而且改變了結(jié)晶形態(tài)（不同于靜態(tài)條件下的球晶，流場(chǎng)作用下可獲得高度異相的串晶結(jié)構(gòu)）?；谏鲜鰧?shí)驗(yàn)結(jié)果，許多學(xué)者提出了不少關(guān)于流動(dòng)誘導(dǎo)結(jié)晶的動(dòng)力學(xué)模型，其建立大多基于Nakamura方程［3］和Kolmogorov方程［4］。如，Doufas等［5］、Tanner等［6-7］通過(guò)修正Nakamura方程得到了流動(dòng)誘導(dǎo)結(jié)晶動(dòng)力學(xué)模型；Koscher等［8］、Zheng等［9］認(rèn)為流動(dòng)提高了晶體的成核數(shù)，并將其反映在Kolmogorov方程中，得到了流動(dòng)誘導(dǎo)結(jié)晶動(dòng)力學(xué)模型。該類方法的主要缺點(diǎn)在于不能提供較為詳細(xì)的微觀結(jié)構(gòu)信息。Eder等［10］首次提出了基于結(jié)晶形態(tài)的數(shù)學(xué)模型。他們把球晶看成是不斷長(zhǎng)大的球，串晶看成是不斷長(zhǎng)大的圓柱體，得到一系列微分方程組。Zuidema等［11］對(duì)其進(jìn)行了改進(jìn)，用可恢復(fù)形變替代Eder模型中流場(chǎng)應(yīng)變率作為串晶成核的驅(qū)動(dòng)。他們的模型較為成功地反映了串晶的微觀信息，但相關(guān)研究并沒(méi)有顯式地采用模擬方法給出其成核、生長(zhǎng)與碰撞的演化細(xì)節(jié)。因此，仍然擺脫不了使用結(jié)晶動(dòng)力學(xué)模型。Guo等［12］引入了大分子鏈形變因子，他們認(rèn)為當(dāng)形變因子大于某臨界值時(shí)，結(jié)晶按串晶形式生長(zhǎng)；反之，按球晶形式生長(zhǎng)。他們利用該方法成功預(yù)測(cè)了iPP注塑制件的皮層厚度，但未得到結(jié)晶形態(tài)細(xì)節(jié)。

國(guó)內(nèi)關(guān)于流動(dòng)誘導(dǎo)結(jié)晶數(shù)值模擬方面的研究也不多見(jiàn)。王錦燕等［13］、榮彥等［14］對(duì)剪切流場(chǎng)中的聚合物結(jié)晶進(jìn)行了模擬，定性分析了剪切對(duì)結(jié)晶過(guò)程及流體性質(zhì)等的影響。他們的研究并未涉及到結(jié)晶形態(tài)。Zhou等［15］采用Guo等［12］的方法對(duì)成型中的結(jié)晶形態(tài)進(jìn)行了顯示，區(qū)分了球晶和串晶，但相關(guān)研究并沒(méi)有將其與結(jié)晶動(dòng)力學(xué)建立聯(lián)系，而只是給出了結(jié)晶形態(tài)的可能表示。

本文基于Eder模型推導(dǎo)了剪切流場(chǎng)中球晶、串晶的形態(tài)演化模型，并采用Monte Carlo法［16］成功模擬了二維剪切流場(chǎng)中聚乙烯的結(jié)晶過(guò)程，給出了球晶、串晶的形態(tài)演化，分析了剪切時(shí)間、剪切速率對(duì)串晶成核密度、結(jié)晶速率、流體黏度、體系應(yīng)力的影響。

1　模型與算法

1.1結(jié)晶動(dòng)力學(xué)模型

Avrami方程［17］是聚合物結(jié)晶動(dòng)力學(xué)模型中應(yīng)用最廣的模型，它也可用于描述流場(chǎng)中的聚合物結(jié)晶，其具體表達(dá)式為

其中，0φ為球晶的“虛幻體積”，0ψ為串晶的“虛幻體積”。Eder等［10］給出了0φ和0ψ的表達(dá)式，其中，描述三維球晶生長(zhǎng)的方程組為

這里，Ns、Rtot、Stot、Vtot分別為球晶的晶核數(shù)、總體直徑長(zhǎng)度、總體表面積、總體體積；而a為成核速率相關(guān)的參數(shù)，Gs為球晶的生長(zhǎng)速率。三維串晶生長(zhǎng)的方程組為

這里，Ns-k、Ltot、S~tot、V~tot分別為串晶的晶核數(shù)、總體長(zhǎng)度、總體表面積、總體體積；而τn為與溫度相關(guān)的串晶成核率參數(shù)，為串晶成核的驅(qū)動(dòng)參數(shù)，γ.為應(yīng)變率，為擬合參數(shù)；τl為與溫度相關(guān)的串晶軸向生長(zhǎng)率參數(shù)，為串晶軸向生長(zhǎng)率的驅(qū)動(dòng)參數(shù)，為擬合參數(shù)，Gs-k，r為串晶徑向生長(zhǎng)速率。

Eder模型中將球晶視為不斷長(zhǎng)大的球，而串晶視為不斷長(zhǎng)大的圓柱體。本文研究二維情況下的晶體生長(zhǎng)，此時(shí)，認(rèn)為球晶按圓形生長(zhǎng)，而串晶按長(zhǎng)方形生長(zhǎng)。

由式（2）可導(dǎo)出三維球晶的等價(jià)微分方程組，將其降為二維，則有

從中可知，球晶的生長(zhǎng)由成核數(shù)sN及生長(zhǎng)速率sG所決定。采用如下的成核密度公式［2］

其中，sN是過(guò)冷度TΔ的函數(shù)，過(guò)冷度TΔ定義為為平衡熔點(diǎn)，0N和?是經(jīng)驗(yàn)參數(shù)。球晶的生長(zhǎng)速率用Hoffman-Lauritzen方程表示［18］

其中，G0、Kg為參考因子，U?為聚合物的分子活化能，Rg為氣體常數(shù)，T∞=Tg-30，Tg為玻璃化轉(zhuǎn)變溫度，

由式（3）可導(dǎo)出三維串晶的等價(jià)微分方程組，同樣，將其降為二維，有

這里，認(rèn)為參數(shù)nτ=∞，lτ=∞［11］，則有

由式（8）可知，串晶的生長(zhǎng)由串晶成核密度Ns-k、串晶軸向生長(zhǎng)速率Gs-k，l以及串晶徑向生長(zhǎng)速率Gs-k，r決定。串晶軸向生長(zhǎng)速率Gs-k，l，可記為［10］

一般情況下，認(rèn)為串晶徑向生長(zhǎng)速率Gs-k，r與球晶生長(zhǎng)速率Gs一致［11］，即

問(wèn)題的關(guān)鍵在于求解串晶成核密度。采用如下公式［8］

其中，C為參數(shù)，N1為體系的第一法向應(yīng)力差。

1.2無(wú)定形相與半結(jié)晶相的模型

式（11）中出現(xiàn)了第一法向應(yīng)力差，因此，必須對(duì)晶體所在體系做一下闡述。這里，采用Zheng等［9］的思想，認(rèn)為半結(jié)晶相懸浮于無(wú)定形相之中，無(wú)定形相采用FENE-P彈性啞鈴模型描述，而半結(jié)晶相采用棒狀啞鈴模型來(lái)描述。

無(wú)定形相的FENE-P彈性啞鈴模型為［19］

其中，C為大分子鏈構(gòu)型張量，λa（T）為分子松弛

其中，，0aλ為溫度為0T時(shí)的分子松弛時(shí)間，ag/ER為擬合常數(shù)。

無(wú)定形相對(duì)應(yīng)力的貢獻(xiàn)為［9，19］

其中，aτ為無(wú)定形相應(yīng)力，n為分子密度，k為Botzmann常數(shù)。

半結(jié)晶相中的分子不能拉伸，只能取向，故采用硬棒模型來(lái)描述［9］，其動(dòng)力學(xué)方程為

其中，為分子構(gòu)型張量，γ.為形變率張量，λsc（x，T）為棒狀啞鈴的松弛時(shí)間，采用如下公式［9］α為相對(duì)結(jié)晶度，A為經(jīng)驗(yàn)值。式（15）中的四階張量<>RRRR，采用如下的二次封閉近似［9］半結(jié)晶相對(duì)體系應(yīng)力的貢獻(xiàn)為［9］

這里，ηsc（x，T）為半結(jié)晶相的黏度，與無(wú)定形相的黏度間滿足如下關(guān)系［9］

至此，整個(gè)體系就完整描述出來(lái)。由式（14）、式（18）可得到體系的整體應(yīng)力表達(dá)式而式（11）中的第一法向應(yīng)力差N1=τ11-τ22由式（20）計(jì)算得到。

1.3算法分析

1.3.1Monte Carlo法本文采用Monte Carlo法［16］來(lái)捕捉結(jié)晶形態(tài)的演化。這里，研究對(duì)象設(shè)為空間中某塊特定區(qū)域上［0，1］mm×［0，1］mm的聚合物。給定溫度、剪切速率及剪切時(shí)間，并按式（5）、式（6）計(jì)算球晶成核密度Ns、球晶生長(zhǎng)速率Gs；按式（9）～式（11）計(jì)算串晶軸向生長(zhǎng)速率Gs-k，l、串晶徑向生長(zhǎng)速率Gs-k，r以及隨時(shí)間變化的串晶成核密度Ns-k。

實(shí)施時(shí)，先將區(qū)域剖分為Ntot個(gè)小單元，這里取Ntot=106，并對(duì)每個(gè)單元賦底色。在t=0時(shí)生成密度為Ns的球晶晶核和密度為Ns-k的串晶晶核，并對(duì)晶核賦顏色。在t時(shí)刻，計(jì)算球晶的半徑R=Gst，對(duì)于每一個(gè)球晶，生成足夠多的隨機(jī)點(diǎn)，若該點(diǎn)的顏色為底色，則判斷其坐標(biāo)是否落入該球晶范圍內(nèi)，若是則認(rèn)為該點(diǎn)被此球晶覆蓋，賦上該球晶的顏色；對(duì)于串晶，由于其成核與時(shí)間相關(guān)，先由式（11）計(jì)算串晶的成核密度，隨機(jī)產(chǎn)生新生核的坐標(biāo)并賦值，同時(shí)，計(jì)算串晶的長(zhǎng)度l=2Gs-k，l及半徑R=Gs-k，r，對(duì)于每一個(gè)串晶，生成足夠多的隨機(jī)點(diǎn)，判斷其是否落入某串晶長(zhǎng)度及半徑覆蓋的長(zhǎng)方形中，若是則賦上該串晶的顏色。如此循環(huán)，直到所有單元都被晶體覆蓋，則結(jié)束。需要指出的是，上述t~的計(jì)算與串晶的晶核生成時(shí)間相關(guān)，從晶核生成后開(kāi)始計(jì)時(shí)。此外，相對(duì)結(jié)晶度的計(jì)算，這里采用式（21）來(lái)表示。

其中，晶體單元數(shù)可根據(jù)單元顏色來(lái)統(tǒng)計(jì)，而總體單元數(shù)為totN。Monte Carlo法的優(yōu)勢(shì)在于可以避開(kāi)使用結(jié)晶動(dòng)力學(xué)模型。常見(jiàn)的結(jié)晶動(dòng)力學(xué)模型，如Avrami方程和Kolmogorov方程，在處理球晶、串晶共混體系時(shí)遇到困難，很難確定形狀因子mC和n。而Monte Carlo法直接從形態(tài)演化中統(tǒng)計(jì)獲得可信的相對(duì)結(jié)晶度，避免了由結(jié)晶動(dòng)力學(xué)模型中參數(shù)的不確定而導(dǎo)致的誤差。

1.3.2有限差分法無(wú)定形相與半結(jié)晶相方程的計(jì)算采用有限差分法。演化方程式（12）和式（15）采用時(shí)間一階向前進(jìn)行離散，即

2　結(jié)果與討論

2.1材料參數(shù)

本文以聚乙烯為代表，考察結(jié)晶過(guò)程及流體性質(zhì)。模擬所用參數(shù)為：N0=17.4×106，?=0.155， G0=2.83×102m/s，，Kg=5.5×105K2，Tg=269 K，=2.69×10-7， C=105Pa-1·s-1·m-1，λa，0=4.00×10-2s，T0=476.15 K，Ea/Rg=5.602×103K，b=5，n= 1.26×1026/m3，k=1.38×10-23，β=9.2，β1=0.05，A=0.44。其中，結(jié)晶形態(tài)參數(shù)見(jiàn)文獻(xiàn)［8，11］，無(wú)定形相及半結(jié)晶相參數(shù)見(jiàn)文獻(xiàn)［9］。

2.2算法有效性驗(yàn)證

為了驗(yàn)證算法的有效性，將模擬結(jié)果與Eder模型預(yù)測(cè)結(jié)果進(jìn)行比較。這里，假定球晶和串晶晶核瞬時(shí)產(chǎn)生，且有成核密度Ns=107m-2，Ns-k=107m-2，球晶按生長(zhǎng)速率Gs=10-6m·s-1進(jìn)行生長(zhǎng)，串晶按徑向生長(zhǎng)速率Gs=10-6m·s-1及軸向生長(zhǎng)速率Gs-k，l=10-5m·s-1進(jìn)行生長(zhǎng)。圖1給出了模擬結(jié)果與理論模型預(yù)測(cè)所得相對(duì)結(jié)晶度的比較。由圖可知，模擬結(jié)果與理論預(yù)測(cè)值吻合得較好。模擬結(jié)果主要在結(jié)晶后期發(fā)生偏離，這可能是由于Eder模型的建立是基于Avrami方程的，而該方程在建模時(shí)并未考慮到晶體間的碰撞。這也是Avrami方程的主要局限性。

圖1　數(shù)值解與Eder模型的比較Fig.1　Comparison of simulation result with Eder model

2.3剪切時(shí)間的影響本節(jié)討論剪切時(shí)間的影響。為了便于分析，這里固定剪切速率γ.=10 s-1，溫度T=137℃。

2.3.1結(jié)晶速率及形態(tài)的變化圖2為剪切時(shí)間t=0、1、5、10 s時(shí)所生成的串晶成核密度。其中，t=0代表不對(duì)流體進(jìn)行剪切，即靜態(tài)條件。施加剪切后，串晶的成核密度將有所增加；且剪切時(shí)間越大，串晶成核密度越大。當(dāng)剪切停止后，成核密度將不再隨時(shí)間的變化而變化。

圖2　不同剪切時(shí)間下的串晶成核密度Fig.2　Nucleation density of shish-kebabs with different shear time

圖3　不同剪切時(shí)間下的相對(duì)結(jié)晶度比較Fig.3　Relative crystallinity comparison with different shear time

圖3為不同剪切時(shí)間下的相對(duì)結(jié)晶度演化。與靜態(tài)條件（t=0）相比，剪切時(shí)間越大，結(jié)晶速率越快。結(jié)晶的加速作用主要是由剪切誘導(dǎo)的串晶所貢獻(xiàn)的。由于流場(chǎng)的剪切作用，促使串晶晶核的產(chǎn)生，并提供軸向的生長(zhǎng)速率，使串晶成長(zhǎng)，進(jìn)而加速結(jié)晶過(guò)程。在剪切時(shí)間較小的情況下，所提供的串晶晶核相對(duì)較少，其加速作用不甚明顯。文獻(xiàn)［9］也報(bào)道了大剪切時(shí)間對(duì)結(jié)晶速率的加速作用。一些學(xué)者［13-14］認(rèn)為這種加速不是無(wú)限制的，當(dāng)剪切時(shí)間達(dá)到某臨界值時(shí)，加速作用會(huì)變得緩慢。

圖4給出了剪切時(shí)間為0、5、10 s時(shí)的結(jié)晶形態(tài)比較。由圖可知，隨著剪切時(shí)間的增加，結(jié)晶形態(tài)發(fā)生較大變化：不僅提高了串晶成核密度，而且提高了串晶的各向異性。這種變化趨勢(shì)與Zhou等［15］的實(shí)驗(yàn)結(jié)果也是一致的。

2.3.2流體黏度及應(yīng)力的變化圖5給出了剪切時(shí)間t=0、1、5、10 s時(shí)流體黏度的變化。黏度隨著時(shí)間的增加在某臨界值時(shí)會(huì)發(fā)生突增現(xiàn)象。這是由相對(duì)結(jié)晶度引起的。由式（19）知，ηsc（x，T）=，當(dāng)α接近常數(shù)A時(shí)，ηsc將趨向于無(wú)窮；故體系黏度急劇變化。此外，剪切時(shí)間越大，流體黏度越早發(fā)生突增。這與Zheng等［9］的結(jié)果一致。

圖6為剪切時(shí)間t=10 s時(shí)的剪切應(yīng)力和第一法向應(yīng)力差。應(yīng)力值在前期變化較為緩慢，而到臨界值附近變化劇烈，該臨界值對(duì)應(yīng)于相對(duì)結(jié)晶度。當(dāng)結(jié)晶完成時(shí)，這些應(yīng)力將鎖在材料內(nèi)部，形成殘余應(yīng)力［20］。

2.4剪切速率的影響

本節(jié)討論剪切速率的影響。為了便于分析，這里固定剪切時(shí)間t=10 s，溫度T=137℃。

2.4.1結(jié)晶速率及形態(tài)的變化圖7給出了剪切速率0γ=.、1、5、10 s-1時(shí)串晶成核密度的比較。其中，γ=.0代表不對(duì)流體進(jìn)行剪切，即靜態(tài)條件。由圖可知，剪切速率越大，串晶成核密度越大。當(dāng)剪切停止后，成核密度將不再隨時(shí)間的變化而變化。

圖8為剪切速率γ=.0、1、5、10 s-1時(shí)的相對(duì)結(jié)晶度發(fā)展。由圖可知，由于剪切作用的存在，結(jié)晶速率將比靜態(tài)條件下（0γ=.）的結(jié)晶速率有所增加，但這種增加在小剪切速率下（γ=.1 s-1）較微弱。同樣，結(jié)晶的這種加速作用也是由剪切誘導(dǎo)的串晶所貢獻(xiàn)的。

結(jié)晶形態(tài)的演化受剪切速率的影響與受剪切時(shí)間的影響類似。隨著剪切速率的增加，串晶成核密度將有所增加，并且串晶的各向異性也更為明顯。這里將不再作圖顯示，圖形與圖4類似。

2.4.2流體黏度的變化圖9為不同剪切速率下流體黏度的變化。流體隨時(shí)間的變化首先較為緩慢，但到某臨界值時(shí)發(fā)生突增，這也是由結(jié)晶度所引起的。此外，剪切速率越大，流體黏度突增的時(shí)間點(diǎn)越早。這與Zheng等［9］的結(jié)果也是一致的。

圖4　不同剪切時(shí)間下結(jié)晶形態(tài)的比較Fig.4　Morphology comparison with different shear time

3　結(jié)　論

本文給出了一個(gè)二維簡(jiǎn)單剪切流場(chǎng)下聚乙烯結(jié)晶過(guò)程的建模與模擬?；贓der模型推導(dǎo)了球晶、串晶的形態(tài)演化模型；并用Zheng的兩相懸浮模型描述聚合物流體，將體系的第一法向應(yīng)力差作為串晶成核的驅(qū)動(dòng)；與此同時(shí)，構(gòu)建了模擬球晶、串晶形態(tài)演化的Monte Carlo法，分析了剪切時(shí)間、剪切速率對(duì)結(jié)晶速率、結(jié)晶形態(tài)、體系黏度等的影響。得出如下結(jié)論。

（1）本文所建立的Monte Carlo法合理有效，成功區(qū)分了晶體生長(zhǎng)模式，捕捉了晶體生長(zhǎng)前沿，較為直觀地給出了晶體的演化發(fā)展；同時(shí)，統(tǒng)計(jì)獲得了可信的相對(duì)結(jié)晶度。該方法的優(yōu)勢(shì)在于成功避開(kāi)了使用結(jié)晶動(dòng)力學(xué)模型，避免由結(jié)晶動(dòng)力學(xué)模型中參數(shù)的不確定而導(dǎo)致預(yù)測(cè)的失真。

（2）剪切時(shí)間對(duì)結(jié)晶速率、結(jié)晶形態(tài)影響重大；剪切時(shí)間越長(zhǎng)結(jié)晶速率越快，串晶對(duì)結(jié)晶形態(tài)的影響越明顯，各向異性更為顯著。剪切時(shí)間對(duì)流體黏度也有顯著影響，這種影響主要體現(xiàn)在對(duì)結(jié)晶度的依賴上，當(dāng)結(jié)晶度接近于某臨界值時(shí)，體系黏度將發(fā)生突增現(xiàn)象。

（3）剪切速率對(duì)結(jié)晶速率、結(jié)晶形態(tài)、體系黏度等的影響與剪切時(shí)間的影響類似；增加剪切速率能提高結(jié)晶速率、提高串晶的各向異性以及使體系黏度發(fā)生突增的時(shí)間點(diǎn)提前。

圖5　不同剪切時(shí)間下流體黏度的變化Fig.5　Evolution of melt viscosity with different shear time

圖6　剪切應(yīng)力和第一法向應(yīng)力差（剪切時(shí)間t=10 s）Fig.6　Shear stress and first normal stress （shear time t=10 s）

圖7　不同剪切速率下的串晶成核密度Fig.7　Nucleation density of shish-kebabs with different shear rate

圖8　不同剪切速率下的相對(duì)結(jié)晶度比較Fig.8　Relative crystallinity comparison with different shear rate

圖9　不同剪切速率下流體黏度的變化Fig.9　Evolution of melt viscosity with different shear rate

符號(hào)說(shuō)明

b——啞鈴的有限拉伸參數(shù)

C——彈性啞鈴構(gòu)型張量，m

Gs，Gs-k，l，Gs-k，r——分別為球晶、串晶軸向、串晶徑

向生長(zhǎng)速率，m·s-1

Ns，Ns-k——分別為球晶、串晶成核密度，m-2

N1——第一法向應(yīng)力差，N·m-2

R——?jiǎng)傂詥♀徣∠蛳蛄浚琺

T——熱力學(xué)溫度，K

t——時(shí)間，s

Δt——時(shí)間步長(zhǎng)，s

α——相對(duì)結(jié)晶度

ηa，ηsc——分別為無(wú)定形相、半結(jié)晶相黏度，Pa·sa

λ，scλ——分別為無(wú)定形相、半結(jié)晶相分子

松弛因子時(shí)間，s

τ，aτ，scτ——分別為總體、無(wú)定形相、半結(jié)晶

相應(yīng)力，N·m-2

?——梯度算子

〈〉——系綜平均

下角標(biāo)

a ——無(wú)定形相

s ——球晶

s-k ——串晶

sc ——半結(jié)晶相

References

［1］KENNEDY P， ZHENG R. Crystallization//Flow Analysis of Injection Molds ［M］. Munich： MunichCarl Hanser Verlag， 2013： 141-153.

［2］ZHE M， LUIGI B， GIUSEPPE P， et al. Flow induced crystallization in isotactic polypropylene during and after flow ［J］. Polym.， 2014， 55：6140-6151.

［3］NAKAMURA K， WATANABE T， KATAYAMA K， et al. Some aspects of nonisothermal crystallization of polymers （Ⅰ）：Relationship between crystallization temperature， crystallinity， and cooling conditions ［J］. J. Appl. Polym. Sci.， 1972， 16： 1077-1091.

［4］KOLMOGOROV A N. On the statistic of crystallization development in metals ［J］. Bull. Akad. Sci. USSR， Class Sci.， Math. Nat.， 1937， 1：355-359.

［5］DOUFAS A K， MCHUGH A J， MILLER C. Simulation of melt spinning including flow-induced crystallization （1）： Model development and predictions ［J］. J. Non-Newtonian Fluid Mech.， 2000， 92： 27-66.

［6］TANNER R I， QI F Z. A comparison of some models for describing polymer crystallization at low deformation rates ［J］. J. Non-Newtonian Fluid Mech.， 2005， 127： 131-141.

［7］WO D L， TANNER R I， FLETCHER D F. A numerical treatment of crystallization in tube flow ［J］. Prog. Polym. Sci.， 2012， 56： 1356-1366.

［8］KOSCHER E， FULCHIRON R. Influence of shear on polypropylene crystallization： morphology development and kinetics ［J］. Polym.，2002， 43： 6931-6942.

［9］ZHENG R， KENNEDY P. A model for post-flow induced crystallization： general equations and predictions ［J］. J. Rheol.， 2004， 48： 823-842.

［10］EDER G， JANESCHITZ-KRIEGL H， LIEDAUER S. Crystallization processer in quiescent and moving polymer melts under heat transfer conditions ［J］. Prog. Polym. Sci.， 1990， 15： 629-714.

［11］ZUIDEMA H， PETERS G W M， MEIJER H E H. Development and validation of a recoverable strain-based model for flow induced crystallization of polymers ［J］. Macromol. Theory Simul.， 2001， 10：447-460.

［12］GUO X， ISAYEV A I， GUO L. Crystallinity and microstructure in injection modelings of isotactic polypropylenes （Ⅰ）： a new approach to modeling and model parameters ［J］. Polym. Eng. Sci.， 1999， 39（10）： 2096-2114.

［13］王錦燕， 陳靜波， 劉春太， 等. 聚合物流動(dòng)誘導(dǎo)結(jié)晶數(shù)值模擬 ［J］.化工學(xué)報(bào)， 2011， 62 （4）： 1150-1156. WANG J Y， CHEN J B， LIU C T， et al. Numerical simulation of shear-induced crystallization of polymer ［J］. CIESC Journal， 2011， 62（4）： 1150-1156.

［14］榮彥， 賀惠萍， 曹偉， 等. 基于兩相模型的聚合物流動(dòng)誘導(dǎo)結(jié)晶數(shù)值模擬 ［J］. 化工學(xué)報(bào)， 2012， 63 （7）： 2252-2257. RONG Y， HE H P， CAO W， et al. Numercial simulation of flow-induced crystallization of polymer based on two-phase model ［J］. CIESC Journal， 2012， 63 （7）： 2252-2257.

［15］ZHOU Y G， TURNG L， SHEN C Y. Modeling and prediction of morphology and crystallinity for cylindrical-shaped crystals during polymer processsing ［J］. Polym. Eng. Sci.， 2010： 1-10.

［16］楊玉良， 張紅東. 高分子科學(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.

［17］AVRAMI M. Kinetics of phase change （Ⅰ）： General theory ［J］. J. Chem. Phys.， 1939， 7： 1103-1112.

［18］HOFFMAN J D， MILLER R L. Kinetics of crystallization from the melt and chain folding in polyethylene fractions revisited： theory and experiment ［J］. Polym.， 1997， 38： 3151-3212.

［19］NASIR A， MUHAMMAD A J， MUHAMMAD S. Theoretical analysis of the exiting thickness of sheets in the calendaring of FENE-P fluid ［J］. J. Non-Newtonian Fluid Mech.， 2015， 225： 28-36.

［20］ABOUHAMZEH M， SINKE J， JANSEN K， et al. Closed form expression for residual stresses and warpage during cure of composite laminates ［J］. Compos. Struct.， 2015， 133： 902-910.

Numerical simulation of morphology and kinetics of polyethylene in shear flow

RUAN Chunlei1，2， LIU Chuntai2
（1School of Mathematics and Statistics， Henan University of Science and Technology， Luoyang 471023， Henan， China；2National Engineering Research Center for Advanced Polymer Processing Technology， Zhengzhou 450001， Henan， China）

The mathematical model of morphology evolution of spherulites and shish-kebabs in the shear flow is deduced based on the Eder model. The model considers that the nucleation and growth of spherulites are determined by the static temperature， while the nucleation and growth of shish-kebabs are determined by the flow，which is depended on the first normal stress difference of the system and the shear rate respectively. In order to calculate the nucleation density of the shish-kebabs， the two-phase suspension model of Zheng is introduced. The model treats the stress as the combination of the amorphous phase and semi-crystalline phase. The amorphous phase is described by FENE-P model while the semi-crystalline phase is depicted by a rigid dumbbell model. Based on the mathematical model， the Monte Carlo method and the finite difference method are constructed，respectively. The former is to capture the crystal growth while the latter is to calculate the equation of the system. By using these methods， the 2D crystallization in the shear flow is simulated. The evolution of spherulites and shish-kebabs is given. Also， the effects of shear time and shear rate on the nucleation density of shish-kebabs，crystallization rate， the viscosity of the fluid and the system stress are discussed. Numerical results show that the Monte Carlo method is valid which not only captures the morphology evolution of crystals successfully， but alsopredicts the crystallization rate well. In addition， the increase of the shear time or the shear rate will increase the nucleation density of shish-kebabs and the crystallization rate.

date： 2015-09-21.

RUAN Chunlei， ruanchunlei622@mail.nwpu. edu.cn

supported by the National Natural Science Foundation in China （11402078， 51375148）.

flow induced crystallization； morphology； kinetics； Monte Carlo simulation； crystallization

O631.+3

A

0438—1157（2016）05—2144—08

2015-09-21收到初稿，2016-01-11收到修改稿。

聯(lián)系人及第一作者：阮春蕾（1983—），女，博士，副教授。

國(guó)家自然科學(xué)基金項(xiàng)目（11402078，51375148）。