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

彭亮，毛偉，趙美英

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

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

彭亮，毛偉，趙美英

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

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

固化過(guò)程；邊界條件；ABAQUS；有限元法

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

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

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

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

1　熱傳導(dǎo)控制方程和固化動(dòng)力學(xué)方程

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

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

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

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

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

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

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

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

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

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

2　基于ABAQUS用戶子程序的溫度場(chǎng)和固化度場(chǎng)模擬

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

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

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

溫度場(chǎng)和固化度場(chǎng)的耦合求解如圖1所示。

圖1　溫度場(chǎng)與固化度場(chǎng)耦合求解流程圖

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

3　算例及程序驗(yàn)證

根據(jù)上述溫度場(chǎng)和固化度場(chǎng)的耦合求解流程，基于用戶子程序?qū)BAQUS進(jìn)行二次開(kāi)發(fā)，編寫(xiě)了用于計(jì)算樹(shù)脂基復(fù)合材料固化過(guò)程的子程序。本文采用了2個(gè)算例對(duì)不同固化動(dòng)力學(xué)模型進(jìn)行了驗(yàn)證。

3.1算例1

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

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

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

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

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

3.2算例2

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

圖2　玻璃纖維/聚酯樹(shù)脂層合板中心點(diǎn)溫度、固化度歷程

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

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

表4　AS4/3501-6的固化動(dòng)力學(xué)參數(shù)

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

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

4　分析與討論

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

4.1邊界條件的影響

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

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

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

4.2 厚度的影響

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

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

4.3升溫速率的影響

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

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

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

5　結(jié) 論

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

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

3）由厚度對(duì)溫度場(chǎng)和固化度場(chǎng)的影響表明，層板越厚，固化開(kāi)始得越晚，中心點(diǎn)溫度超過(guò)保溫平臺(tái)溫度越晚，溫度峰值越大，完成固化所需時(shí)間越短。

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

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

參考文獻(xiàn)：

［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］ 楊正林，陳浩然.層合板在固化過(guò)程中瞬態(tài)溫度場(chǎng)及固化度的有限元分析［J］.玻璃鋼/復(fù)合材料，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］ 張紀(jì)奎，酈正能，關(guān)志東，等.熱固性復(fù)合材料固化過(guò)程三維有限元模擬和變形預(yù)測(cè)［J］.復(fù)合材料學(xué)報(bào)，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基金項(xiàng)目：國(guó)家自然科學(xué)基金（11502205）資助

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