橡膠非等溫硫化工藝的數值模擬

張韜杰，施 飛

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

橡膠硫化是決定橡膠制品質量的關鍵工序。在硫化過程中，經過高溫加熱，橡膠流體在硫化劑的作用下發生交聯反應，分子結構由鏈式結構變成網狀結構，經一段時間的保壓，制品定型，開啟模具，取出制品。

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

1 基本數學模型

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

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

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

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

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

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

在第三類邊界條件處：

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

橡膠硫化的反應方程為：

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

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

2 離散微分方程

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

圖1 四面體網格的控制體單元

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

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

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

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

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

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

其中：

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

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

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

3 計算案例

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

圖2 計算流程

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

表1 橡膠制品材料的熱力學參數

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

圖3 零件圖及劃分的網格

圖4 熱熔系數隨溫度的變化曲線

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

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

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

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

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

4 結論

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