橡膠非等溫硫化工藝的數(shù)值模擬

張韜杰，施 飛

（寧波工程學院 機械工程學院，浙江 寧波 315016）

橡膠硫化是決定橡膠制品質(zhì)量的關鍵工序。在硫化過程中，經(jīng)過高溫加熱，橡膠流體在硫化劑的作用下發(fā)生交聯(lián)反應，分子結(jié)構(gòu)由鏈式結(jié)構(gòu)變成網(wǎng)狀結(jié)構(gòu)，經(jīng)一段時間的保壓，制品定型，開啟模具，取出制品。

橡膠硫化工藝自從1839年被Goodyear發(fā)現(xiàn)以后，多年來一直吸引了眾多研究人員的注意，特別是利用數(shù)值模擬來研究硫化過程。J.C.Ambelang等[1]和H.P.Schlanger[2]利用有限差分方法模擬了輪胎的硫化過程，而G.D.Hubbard等[3]采用有限元法進行數(shù)值模擬研究，但他們沒有考慮到輪胎從硫化機中取出后在空氣中冷卻期間的硫化。M.Javadi等[4]和M.H.R.Ghoreishy等[5]以 三 維 瞬時熱平衡方程及經(jīng)驗的硫化動力學模型建立了有限元模型來模擬橡膠硫化工藝，通過實驗來驗證該模型的準確性和可靠性，并討論了優(yōu)化硫化時間的一些措施。R.Darwis等[6]基于Gauss-Seidel Red Black和Gauss-Seidel Method的方法開發(fā)了數(shù)值模擬工具，研究了計算硫化過程中溫度場和硫化度場的變化。P.M.Abhilasha等[7]模擬了具有各向異性溫度歷史的橡膠制品平板硫化的過程，該模型最有意義之處在于能對硫化成型制品的性能進行預測。在國內(nèi)，劉莉莉等[8]以遺傳算法優(yōu)化了硫化工藝，動態(tài)地模擬了非線性、非穩(wěn)態(tài)、多相、強耦合的樹脂硫化過程。閆相橋等[9-10]模擬了非均質(zhì)性、橡膠基復合材料傳熱性能的各向異性、輪胎硫化時橡膠材料伴有的反應生熱及輪胎硫化的實際升溫歷程。本工作以控制體/有限元（CV/FEM）為基礎，對硫化過程的熱擴散方程和硫化反應方程進行離散，建立控制方程組，并編寫代碼，進行數(shù)值模擬。

1 基本數(shù)學模型

由于在硫化階段沒有橡膠流動就沒有對流傳熱，因此僅考慮橡膠及嵌件與模具之間的熱擴散以及橡膠硫化所產(chǎn)生的熱量。控制方程包括兩類，即：通過熱傳導和熱傳遞方式進行的熱擴散方程以及橡膠硫化方程。

對于三維瞬時熱擴散，有：

式中，方程左邊描述了橡膠的能量變化，右邊第1—3項是橡膠在3個方向上的熱傳導，右邊第4項描述了橡膠硫化而放出的熱能。其中，ρ為密度；cp為比熱容；T為溫度；t為時間；kLx，kLy，kLz為橡膠熱傳導系數(shù)；s˙為硫化反應產(chǎn)生的熱量，有：

式 中，ΔH為 橡 膠 反 應 熱，而G(α,T) =dαdt（Kamal模型[11-13）]。

與之相關的溫度邊界條件主要分為兩類，即第一類邊界條件和第三類邊界條件。

在第一類邊界條件處：T=Tm

在第三類邊界條件處：

式中，nx，ny，nz為邊界上單位法向量；h為熱導率；T∞為邊界處模具的溫度。

橡膠硫化的反應方程為：

式中，α為硫化度，m為常數(shù)，k為Arrhenius類型的溫度常數(shù)，可以表示如下：

式中，k0為常數(shù)，E為活化能，R為氣體常數(shù)。在求解硫化反應方程時，其初始條件硫化度不能為0，否則方程的解無法收斂。

2 離散微分方程

基于控制體/有限元理論，單元網(wǎng)格按照以下方法進行劃分（詳見圖1）。從圖1可以看出，首先將橡膠制品劃分為四面體單元，再對四面體單元進行細劃分并重新裝配組合成新單元，即控制體單元，其既具有有限體積法的一些特性，也具有有限元法的一些功能，可在計算時充分使用它們的特點。

圖1 四面體網(wǎng)格的控制體單元

圖1中O是四面體單元ABCD的中心，點1，2和3分別是與C相鄰的面ABC，BCD，ACD的面中心，點4，5和6分別是與C相鄰的邊BC，AC，DC的邊中點。這樣四面體OC26就是C節(jié)點控制體與四面體單元ABCD相重合的一部分，同樣在與C點相關的其他四面體單元里也具有該相重合的部分，依次將這些重合部分組合疊加起來就組成了節(jié)點C的控制體單元。在其他節(jié)點處也進行類似運算，就可以得到橡膠制品的三維控制體網(wǎng)格，控制體單元之間沒有孔隙，遵循質(zhì)量守恒定律。此種單元網(wǎng)格各控制體之間關系并非固定，故也稱為非結(jié)構(gòu)化控制體網(wǎng)格，對各種形狀的零件具有很強的適應性。需要注意的是，在網(wǎng)格劃分過程中，最長邊與最短邊的比例應當盡量接近于1.0，一般不能超過3.0，否則有可能造成計算過程的錯誤。

按照以上方法，可以得到：

式中，δ可表示為溫度、硫化度等場量；Ni為形函數(shù)。

對于三維瞬時熱擴散方程（1）的離散，在空間上離散是通過差分格式進行的，考慮到求解過程的穩(wěn)定性，此處采用迎風格式。而對于時間上的離散，此處主要采用全隱式方法進行。各項分別積分得到：

式中，Sc和Sp為s˙進行線性化后的常數(shù)；0表示前一次的迭代結(jié)果；ΔV為控制體體積；Δt為時間步長。

對積分后的式（4），（5）和（6）進行重新組合[14]，得到如下關系式：

其中：

式中，N為與當前控制體相鄰的控制體數(shù)目；nxi為控制體單元表面向外的單元法向n的各個分量，λxi為當前控制體和相鄰控制體之間連線向量λ的分量。

而溫度邊界條件的離散，對于第一類邊界條件可以直接使用，對于第三類邊界條件，可以在相關邊界處積分，如式（11）所示。

對于硫化的反應方程（2），可以離散如式（12）所示。

3 計算案例

得到離散方程（7）和離散方程（12）及相關的邊界條件后，設計了該數(shù)值分析的計算流程圖（見圖2），并編寫了代碼。

圖2 計算流程

本算例的橡膠制品由橡膠與金屬嵌件組成并放置在金屬模具中，上下面是第一類邊界條件，其余為第三類邊界條件。橡膠材料的各參數(shù)[5]分別 為：k01.4×1025，m2.91，ΔH4.54×106J·m-3，E260 622 J·mol-1，T0293 K。表1所示為橡膠制品材料的熱力學參數(shù)[5]。

表1 橡膠制品材料的熱力學參數(shù)

本算例的橡膠制品圖及所劃分的四面體網(wǎng)格見圖3，熱熔系數(shù)隨溫度變化曲線見圖4。在計算中選取兩點A和D進行比較，位置參數(shù)及邊界條件見文獻[5]。

圖3 零件圖及劃分的網(wǎng)格

圖4 熱熔系數(shù)隨溫度的變化曲線

圖5所示為本工作代碼計算模擬曲線與文獻[5]結(jié)果曲線比較。從圖5可以看出，本文的計算結(jié)果與已發(fā)表文獻的數(shù)據(jù)是一致的，從而說明了該計算方法是正確的。但預測數(shù)據(jù)與文獻數(shù)據(jù)尚有些差異，其主要原因在于數(shù)值計算方法的不同。

圖5 點A和D的硫化度和溫度變化曲線

從圖5可以看出，在1 800 s左右的時候，A和D兩點的硫化度開始迅速增加，此時溫度在150 ℃附近，即對于這種橡膠材料，硫化溫度在150 ℃左右是合適的。

圖6是預測的橡膠制品中面在Z向的硫化歷史。由圖6可知，硫化過程的發(fā)展是非常迅速的，硫化度從接近0到0.99的時間在500 s左右。從圖6還可以看出橡膠制品各部分硫化質(zhì)量的差異。因此，通過對硫化工藝的優(yōu)化可以很好地控制硫化過程，縮短硫化時間，提高硫化質(zhì)量。

圖6 預測的零件中面在Z向的硫化歷史

4 結(jié)論

研究發(fā)現(xiàn)，模擬橡膠硫化工藝，可以得到零件內(nèi)溫度和硫化度分布圖，而無需通過多次實驗來測定，減少時間和金錢的浪費，并可以據(jù)此優(yōu)化硫化工藝，降低零件不均勻的溫度場和硫化度場，有利于減少零件變形，提高零件的硫化質(zhì)量。