彭亮，毛偉，趙美英

（西北工業(yè)大學航空學院，陜西西安 710072）

復合材料層合板結(jié)構(gòu)固化過程數(shù)值分析方法研究

彭亮，毛偉，趙美英

（西北工業(yè)大學航空學院，陜西西安 710072）

為了研究復合材料層板結(jié)構(gòu)在固化過程中的溫度場與固化度場分布情況，對ABAQUS子程序進行了二次開發(fā)，構(gòu)建了不同固化動力學方程的子程序模塊，實現(xiàn)了層合板固化過程的三維數(shù)值模擬。通過算例驗證，所建立的數(shù)值分析方法能精確求解固化過程中的溫度場和固化度場。基于該方法，進一步討論了層合板邊界條件、厚度、升溫速率對溫度場和固化度場的影響。上述研究方法具有較強的適用性，對改善固化工藝參數(shù)，控制結(jié)構(gòu)固化變形具有一定的指導意義。

固化過程；邊界條件；ABAQUS；有限元法

復合材料層合板在固化階段，不僅層合板外部環(huán)境溫度變化，而且層合板內(nèi)部固化反應產(chǎn)生化學放熱，兩者相互影響，是一個復雜的耦合過程，導致層合板內(nèi)部產(chǎn)生復雜的溫度梯度。這種不均勻分布的溫度場和固化度場，將在復合材料內(nèi)產(chǎn)生熱應力和變形，是復合材料制品中產(chǎn)生諸如尺寸誤差、翹曲等缺陷的根本原因。因此，研究復合材料層合板在固化工藝過程中的溫度和固化度分布及其變化規(guī)律，對提高復合材料工藝質(zhì)量有重要意義［1］。

近年來對樹脂基復合材料固化過程的有限元模擬取得了一定進展。Loos等［2］采用模塊化方法研究層合板固化過程，針對AS4/3501-6單向?qū)雍习褰⒘艘痪S數(shù)值模型；Bogetti等［3］采用二維有限元方法數(shù)值模擬了任意截面形狀和邊界條件層合板的固化過程；楊正林等［4］采用有限元與有限差分相結(jié)合的方法分析了二維層合板模型的固化變形；Cheung、張紀奎等［5-6］用有限元方法研究了三維條件下的固化過程。

由于固化過程中溫度場和固化度場耦合問題的復雜性，研究者基本上都是自行編制程序來實現(xiàn)數(shù)值仿真，增加了研究的難度和周期，而且程序的適用性也受到限制。本文從正交各向異性材料的三維熱傳導方程和固化動力學方程出發(fā)，基于有限元軟件ABAQUS進行二次開發(fā)，實現(xiàn)了樹脂基復合材料固化過程的三維數(shù)值模擬，通過對2種材料體系的固化過程進行了算例驗證，證明了程序的正確性，并討論了層合板邊界條件、厚度、升溫速率對溫度場和固化度場的影響。

本文研究所涉及的玻璃纖維/聚酯樹脂的固化動力學方程屬于n級反應模型［7］，AS4/3501-6的固化動力學方程為自催化反應模型［8］，當采用不同固化動力學方程時，只需依據(jù)利用Fortran語言在ABAQUS子程序的基礎上進行編輯即可完成，因而具有更強的適用性，縮短了程序的編寫周期?？捎糜诜治鰪碗s形狀和邊界條件下結(jié)構(gòu)固化過程的溫度場和固化度場，具有較強的適用性，對改善復合材料結(jié)構(gòu)固化工藝參數(shù)具有指導意義。

1　熱傳導控制方程和固化動力學方程

熱固性樹脂基復合材料的固化過程是一個熱與化學反應相互耦合的過程。復合材料構(gòu)件內(nèi)部的溫度分布由向復合材料的傳熱速率和固化反應產(chǎn)熱速率共同決定，復合材料固化溫度場的分析本質(zhì)上是一個具有非線性內(nèi)熱源的熱傳導問題，其中內(nèi)熱源是樹脂基體固化反應放出的熱量。通過對該熱傳導問題的求解，可以得到復合材料在固化過程中任意時刻、位置的溫度及固化度。由于固化階段樹脂基本不發(fā)生流動，可忽略對流傳熱影響，則根據(jù)Fourier熱傳導定律和能量平衡原理建立該問題的熱傳導控制方程：

式中：ρc、c、kii（i=x，y，z）分別為復合材料的密度、比熱、各向異性的熱傳導系數(shù)，內(nèi)部熱源項q□為樹脂發(fā)生化學反應放出的熱量，可以表示為：

式中：ρr為樹脂密度，Hr為固化反應完成時單位質(zhì)量樹脂放出的總熱量，α為樹脂固化度，t為時間。

樹脂的固化反應決定了熱傳導控制方程中內(nèi)熱源的熱量大小，綜合各類文獻的研究方法，對固化動力學模型的表征方法主要有2種：微觀水平（力學的）和宏觀水平（唯象的）。其中唯象模型以化學反應動力學為主要特征，忽略各組分之間相互作用的細節(jié)，也不關(guān)注樹脂的組成或配方，是最被廣泛使用的方法。

大多數(shù)唯象模型以方程（3）、（4）為基礎：

式中，f（α）為固化機理函數(shù)，由實驗數(shù)據(jù)確定；K（T）為固化速率系數(shù)，用阿累尼烏斯方程表示：

式中：A為頻率因子；E為活化能；R為普適氣體常數(shù)；T為溫度。

根據(jù)樹脂反應機理函數(shù)f（α）的形式，可以將現(xiàn)有模型分成n級反應模型和自催化模型2類，方程（5）和（6）的化學動力學模型分別屬于n級反應模型和自催化反應模型。

式中：k1和k2可用阿累尼烏斯方程表示為：

式中，A1，A2，ΔE1，ΔE2為實驗確定的常數(shù)。

2　基于ABAQUS用戶子程序的溫度場和固化度場模擬

ABAQUS具有強大的非線性分析能力，能很好地解決熱-力耦合問題。本文主要運用ABAQUS/ Standard提供的HETVAL、USDFLD用戶子程序，而在對流換熱條件和溫度邊界條件比較復雜的模型進行定義時，使用了FILM和DISP子程序。

固化過程中樹脂固化放熱的模擬是通過子程序HETVAL定義熱源項來實現(xiàn)的。該子程序用來在傳熱分析中定義內(nèi)部生熱與由于產(chǎn)熱導致的熱流，允許內(nèi)部生熱依賴于狀態(tài)變量；子程序USDFLD對STATEV的任何更新都可傳遞到HETVAL中。

表征化學反應程度的固化度場通過用戶定義場在材料屬性中進行定義，定義中用使用了用戶子程序USDFLD。狀態(tài)變量可以在USDFLD中更新，然后傳遞到子程序HETVAL中。

溫度場和固化度場的耦合求解如圖1所示。

圖1　溫度場與固化度場耦合求解流程圖

在子程序HETVAL中，通過固化動力學方程得到固化速率和固化度，并存儲在狀態(tài)數(shù)組STATEV中，更新后的STATEV數(shù)組傳遞到子程序USDFLD以得到固化度場；固化反應放熱通過FLUX模塊計算，產(chǎn)生的熱量傳遞到ABAQUS熱分析模塊中得到新的溫度場，并利用新溫度計算得到固化速率和固化度，如此循環(huán)，直到固化過程結(jié)束。

3　算例及程序驗證

根據(jù)上述溫度場和固化度場的耦合求解流程，基于用戶子程序?qū)BAQUS進行二次開發(fā)，編寫了用于計算樹脂基復合材料固化過程的子程序。本文采用了2個算例對不同固化動力學模型進行了驗證。

3.1算例1

［5］，采用［0/90］鋪層的玻璃纖維/聚酯樹脂層合板，長寬均為15.24 cm，厚為2.54 cm，共42層。玻璃纖維/聚酯樹脂固化動力學經(jīng)驗公式為：

該固化動力學方程屬于n級動力學模型，熱力學參數(shù)如表1所示，固化動力學參數(shù)如表2所示，其中R為普適氣體常數(shù)，取8.314。

表1　玻璃纖維/聚酯樹脂的熱力學參數(shù)

表2　玻璃纖維/聚酯樹脂的固化動力學參數(shù)

邊界條件用對流換熱系數(shù)h與導熱系數(shù)k11的比值表示。上表面h/k11=87 m-1，下表面h/k11= 125 m-1，四周絕熱。計算得到該層合板中心點溫度和固化度歷程，結(jié)果如圖2a）和圖2b）所示。由結(jié)果可知，溫度和固化度均與實驗結(jié)果符合良好，從而驗證了程序針對n級動力學模型的正確性。

3.2算例2

參考文獻［5］，算例2采用［0/90］鋪層的AS4/ 3501-6層合板，幾何模型和算例1相同，長寬均為15.24 cm，厚為2.54 cm，共42層。AS4/3501-6固化動力學經(jīng)驗公式如下：

圖2　玻璃纖維/聚酯樹脂層合板中心點溫度、固化度歷程

式中，K1、K2及K3分別為3501-6樹脂體系的反應速率常數(shù)，A1、A2及A3分別為其頻率因子，ΔE1、ΔE2及ΔE3分別為其活化能?？梢钥闯?，AS4/3501-6復合材料屬于自催化動力學模型，以分段函數(shù)給出，以固化度達到0.3時作為分界點，0.3前后采用不同的固化動力學模型。

表3　AS4/3501-6的熱力學參數(shù)

表4　AS4/3501-6的固化動力學參數(shù)

圖3　AS4/3501-6層合板中心點溫度、固化度歷程

給定邊界條件與算例1中相同。圖3a）和圖3b）分別為該層合板中心點溫度和固化度歷程，與文獻實驗結(jié)果進行比較，可以看出符合性較好，驗證了程序針對自催化模型時的正確性。

4　分析與討論

由計算結(jié)果可知，邊界條件、厚度、升溫速率對層合板溫度場和固化度場均有影響。以下討論中，除了研究厚度對溫度場和固化度場的影響時，試件厚度不同外，在研究邊界條件、升溫速率對溫度場和固化度場的影響時，層合板厚度均為2.54 cm；中心點指的是中間層（z=t/2）處。

4.1邊界條件的影響

熱傳導問題有3類邊界條件：①Dirichlet邊界條件，即已知溫度的邊界條件；②Neuman邊界條件，即已知熱流密度的邊界條件；③Robin邊界條件，即已知對流換熱系數(shù)的邊界條件。第3類邊界條件中，當對流換熱系數(shù)h→∞時，邊界條件就轉(zhuǎn)化為第1類邊界條件；當h→0時，邊界條件就轉(zhuǎn)化為第2類邊界條件。故可通過改變h的大小來改變邊界條件和特征。本文用h/k的比值變化來表示邊界條件的改變，其中k為邊界法向方向的熱傳導數(shù)，分別取h/k=100 m-1和500 m-1作為層合板上、下表面的邊界條件，側(cè)面仍為絕熱邊界。

圖4　邊界條件對AS4/3501-6層合板溫度、固化度的影響

圖4a）給出了層合板表層（z=0）中心點和中間層（z=t/2）的中心點溫度曲線，由計算結(jié)果可知：隨著h/k比值的增大，層合板升溫速率加快，中心點溫度峰值逐漸減小。圖4b）是對應條件下的固化度曲線，結(jié)果表明h/k越小，放熱期間固化度梯度將會增加。

4.2 厚度的影響

對于前述AS4/3501-6層合板，取1.27 cm、2.54 cm、3.81 cm 3種不同厚度，上下表面為第1類邊界條件，側(cè)面為絕熱邊界。圖5a）和圖5b）為這3種不同厚度層合板固化過程中心點的溫度和固化度變化歷程。由圖5a）的溫度曲線可知，在中心點溫度超過保溫平臺溫度之前，1.27 cm層合板的溫度最高，2.54 cm板次之，3.81 cm板溫度最小，超過保溫平臺溫度最晚，但溫度峰值明顯高于其他2個。由圖5b）的固化度曲線可知，1.27 cm層合板固化反應啟動時間最早，進行時間最長，而3.81 cm層合板固化反應啟動最晚，但最先結(jié)束。所以厚度越大，中心點溫度峰值就出現(xiàn)得越晚，峰值越大，中心點固化開始得就越晚，而完成固化所需時間就越短。

圖5　厚度對AS4/3501-6層合板溫度、固化度的影響

4.3升溫速率的影響

對于前面所述AS4/3501-6層合板，上下表面采用第1類邊界條件，側(cè)面為絕熱邊界，層合板在389 K溫度下固化60 min后，再分別施以1.5 K/min，3 K/min和6 K/min的加熱速率達到449 K。

圖6a）和圖6b）分別給出了不同升溫速率下，層板中心點溫度和固化度隨時間的變化曲線。由圖6a）可知，升溫速率越慢，中心點溫度升溫越慢，峰值越小，峰值過后的降溫過程也越慢；同時由圖6b）可知，升溫速率越慢，固化速率越慢，固化完成所需時間也越長。

圖6　升溫速率對AS4/3501-6層合板溫度、固化度的影響

5　結(jié) 論

1）根據(jù)熱傳導和唯象固化動力學理論，建立了正交各向異性復合材料層合板固化過程的三維有限元模型；研究了基于ABAQUS子程序耦合求解溫度場和固化度場的方法，通過對比算例與文獻實驗結(jié)果，驗證了本文方法的正確性。

2）由邊界條件對溫度場和固化度場的影響表明，h/k越大，層合板升溫越快，中心點溫度峰值越小，放熱期間固化度梯度減小，因而h/k的增大有利于減輕溫度場和固化度場的不均勻性。

3）由厚度對溫度場和固化度場的影響表明，層板越厚，固化開始得越晚，中心點溫度超過保溫平臺溫度越晚，溫度峰值越大，完成固化所需時間越短。

4）升溫速率對溫度場和固化度場的影響表明，升溫速率越慢，中心點溫度升溫越慢，峰值越小，降溫過程也越慢，從而降低了固化速率，使得固化完成所需時間延長，因而適當降低升溫速率有利于減輕溫度場和固化度場的不均勻性，但同時也會延長固化完成時間，二者之間要根據(jù)實際情況權(quán)衡。

5）本文基于ABAQUS子程序所建立的復合材料三維有限元固化數(shù)值模擬方法，適用于模擬具有不同類型動力學方程的復合材料固化過程，可分析復雜形狀和邊界條件下結(jié)構(gòu)固化過程溫度場和固化度場，因此對改善固化工藝參數(shù)具有一定的指導意義。

參考文獻：

［1］ Wang J，Kelly D.Finite Element Analysis of Temperature Induced Stresses and Deformations of Polymer Composite Components ［J］.Journal of Composite Materials，2000，34（17）：1456-1471

［2］ Loos A C，Springer G S.Curing of Epoxy Matrix Composites［J］.Journal of Composite Materials，1983，17（2）：135-169

［3］ Bogetti T A，Gillespie J W.Two Dimensional Cure Simulation of Thick Thermosetting Composite［J］.Journal of Composite Materials，1991，25（3）：239-250

［4］ 楊正林，陳浩然.層合板在固化過程中瞬態(tài)溫度場及固化度的有限元分析［J］.玻璃鋼/復合材料，1997（3）：3-7

Yang Zhenglin，Chen Haoran.The Finite Element Analysis of Transient Fields of the Temperature and Degree of Cure of Laminates during the Whole Curing Process［J］.Fiber Reinforced Plastics/Composites，1997（3）：3-7（in Chinese）

［5］ Cheung A，Yu Y，Pochiraju K.Three-Dimensional Finite Element Simulation of Curing of Polymer Composites［J］.Finite Elements in Analysis&Design，2004，40（8）：895-912

［6］ 張紀奎，酈正能，關(guān)志東，等.熱固性復合材料固化過程三維有限元模擬和變形預測［J］.復合材料學報，2009，26（1）：174-178

Zhang Jikui，Li Zhengneng，Guan Zhidong，et al.Three-Dimensional Finite Element Simulation and Prediction for Process-Induced Deformation of Thermoset Composites［J］.Acta Material Compositae Sinice，2009，26（1）：174-178（in Chinese）

［7］ Boey Fyc，Song Xl，Yue Cy，et al.Modeling the Curing Kinetics for a Modified Bismaleimide Resin［J］.Journal of Polymer Science Part A：Polymer Chemistry，2000，38（5）：907-913

［8］ Park H C，Lee S W.Cure Simulation of Thick Composite Structures Using the Finite Element Method［J］.Journal of Composite Materials，2001，35（3）：188-201

Finite Element Analysis for Curing Process of Composite Laminates

Peng Liang，Mao Wei，Zhao Meiying
（College of Aeronautics，Northwestern Polytechnical University，Xi′an 710072 China）

In this paper，three-dimensional numerical simulation of the curing process of composite laminates is realized with ABAQUS subroutine secondary development method.Being verified by two examples，this method can predict coupled temperature field and degree of cure field relative precisely in the curing process.Moreover，the influence of boundary conditions，thickness and heating rate on the temperature field and the degree of cure field are discussed with this method.When curing kinetic equations are different for different types of composites，just program with the Fortran language on the basis of ABAQUS subroutine according to the corresponding equations；this shortens the program cycle.The applicability of this method is universal and it provides meaningful guidelines to refine the curing process parameters.

curing process；boundary condition，ABAQUS，finite element method，subroutine，heating rate，temperature distribution

V258

A

1000-2758（2015）06-0900-06

2015-04-28基金項目：國家自然科學基金（11502205）資助

彭亮（1982—）西北工業(yè)大學講師、博士研究生，主要從事飛行器設計與適航研究。