(1.蘭州理工大學 材料科學與工程學院，甘肅 蘭州 730050； 2.蘭州理工大學 理學院，甘肅 蘭州 730050)

1 引 言

熱壓罐成型工藝是制造連續(xù)纖維增強熱固性樹脂基復合材料構件的主要方法。由于影響構件成型質(zhì)量的因素較多，而且實驗周期長、成本高，難以保證重復性。在這種情況下，采用數(shù)值模擬法建立一種比較完整的固化過程分析和預測模型尤為重要。

近二十幾年來，國內(nèi)外眾多學者針對復合材料固化過程進行了數(shù)值模擬研究。初期，對層合板固化的數(shù)值模擬多采用簡單的一維和二維模型[1-3]。隨著計算機技術的不斷進步，一些學者建立了復合材料固化的三維有限元分析模型[4-6]，實現(xiàn)了對復雜結構件固化的模擬研究。對復合材料固化過程中溫度場變化情況的研究較多[7]，少數(shù)學者還實現(xiàn)了復合材料固化的多物理場耦合計算[8-9]。

能否準確預測固化過程殘余應力的關鍵在于選擇合適的本構模型。在早期的研究中，彈性模型的運用較多，將其與經(jīng)典層合板理論結合可以較好地預測薄板的殘余應力[10-12]。但是，彈性模型僅僅考慮的是冷卻階段殘余應力的發(fā)展狀況，忽略了冷卻前產(chǎn)生的應力。眾所周知，樹脂在固化過程中會發(fā)生明顯的粘彈性性能的變化，尤其在升溫和保溫階段，這也引發(fā)了許多學者對粘彈性模型的青睞。由于粘彈性模型自身的復雜性，在研究初期，人們只是利用粘彈性模型預測冷卻階段的固化殘余應力[13]。后來，引入時間-溫度疊加原理和等效原理，將材料的熱膨脹、化學收縮和應力松弛等因素考慮在內(nèi)，使得對固化殘余應力發(fā)展變化的預測研究更加深入[14]。

本研究建立了層合板固化的三維有限元分析模型，考慮了材料性能在固化過程中的時變特性，并采用COMSOL多物理場耦合分析軟件對層合板固化過程中的溫度場、固化度場及殘余應力場進行耦合計算分析。

2 數(shù)學模型

2.1 熱化學模型

熱化學模型由熱傳導方程和固化動力學方程構成，可以用來計算復合材料固化過程中的溫度場和固化度場。復合材料固化過程中的熱量主要來源于兩個方面，一是經(jīng)模具輔助材料傳遞給復合材料結構的外界熱量，另一部分是樹脂基體固化過程中化學反應產(chǎn)生的非線性內(nèi)熱源。對于內(nèi)熱源的熱傳導問題，本研究采用Fourier熱傳導定律和能量平衡關系來建立數(shù)學模型。三維瞬態(tài)熱傳導控制方程為[2]：

(1)

式中：ρ、Cp、T和t分別為復合材料的密度、比熱、溫度和時間；kx、ky和kz分別為復合材料在x、y和z方向上的導熱系數(shù)；Q為內(nèi)部熱源，由下式確定：

(2)

樹脂基體的固化反應是一個熱激活的復雜的化學交聯(lián)反應，因此大多數(shù)有關固化反應動力學的方程建立在經(jīng)驗模型基礎上。本研究采用唯象動力學模型，該模型不考慮整個反應過程中的動力學機理，只用簡單的公式來表示反應中各參數(shù)的相互作用關系。本研究主要采用3501-6環(huán)氧樹脂對所建模型進行驗證，其固化動力學方程為[1-2]：

(3)

其中，Ki(i=1,2,3)為反應速率常數(shù)，由下式確定：

(4)

式中：Ai、Ei和R分別為頻率因子、活化能和普適氣體常數(shù)。在實際計算中，由于不能計算零的負整數(shù)次冪，因此初始條件可設為α|t=0=10-10。3501-6環(huán)氧樹脂的固化動力學參數(shù)如表1所示。

表1 3501-6 環(huán)氧樹脂的固化動力學參數(shù)[1-2]Table 1 Cure kinetic parameters of 3501-6 epoxy resin [1-2]

2.2 粘彈性力學模型

盡管粘彈性模型存在一些局限性，例如：對材料實驗數(shù)據(jù)需求量大，占用計算機內(nèi)存，計算時間長等，但是它是描述熱固性樹脂基復合材料固化過程中粘彈性行為最精準的數(shù)值計算方法之一。假設層合板在固化發(fā)生前無應變產(chǎn)生，則粘彈性材料的應力可以由下式確定[15]：

(5)

其中，

εtc=εth+εch=φΔT+φΔα

(6)

式中，C為松弛剛度矩陣，其分量的計算方法參見文獻[3]；εtot和εtc分別為總應變和熱化學應變，熱化學應變包括熱膨脹應變εth和化學收縮應變εch；φ和φ分別為熱膨脹系數(shù)和化學收縮系數(shù)；ξ(t)和ξ(τ) 分別表示當前和過去相對時間，由式(7)確定：

(7)

式中：αref為參考固化度；aT為移位因子。根據(jù)n單元并聯(lián)的廣義麥克斯韋模型，粘彈性復合材料的松弛模量可以表達為[2]：

(8)

(9)

2.3 材料熱物理性能

在固化過程中，復合材料及其組份的熱物理參數(shù)(包括密度、比熱容和導熱系數(shù))是隨溫度和固化度的變化而變化的。

AS4碳纖維的密度為常量1790kg·m-3，3501-6環(huán)氧樹脂的密度由式(10)計算得到[2]：

(10)

表3 AS4/3501-6復合材料組份的彈性力學性能[3]Table 3 Elastic mechanical properties of AS4/3501-6 composite constituents [3]

AS4/3501-6復合材料的密度可以用混合率公式計算如下[15]：

ρ=vfρf+(1-vf)ρr

(11)

式中，ρf為纖維的密度。

AS4碳纖維的比熱容與溫度呈線性關系[2]：

Cpf=1390+4.50T

(12)

3501-6環(huán)氧樹脂的比熱容是溫度和固化度的線性函數(shù)[2]：

Cpr=4184[0.468+5.975×10-4T-0.141α]

(13)

AS4/3501-6復合材料的比熱容用混合率公式計算[15]：

(14)

AS4碳纖維的導熱系數(shù)是溫度的線性函數(shù)[2]：

kf=0.742+9.02×10-4T

(15)

3501-6環(huán)氧樹脂的導熱系數(shù)是溫度和固化度的線性函數(shù)[2]：

kr=0.04184[3.85+(0.035T-0.141)α]

(16)

復合材料的導熱系數(shù)具有各向異性，復合材料平行于纖維方向和垂直于纖維方向的導熱系數(shù)可以用混合率公式計算得到[15]：

kL=vfkf+(1-vf)kr

(17)

(18)

其中，

(19)

3 基于COMSOL的數(shù)值模擬算法

在第2節(jié)中提出的理論模型基礎上，運用COMSOL 多物理場耦合有限元分析軟件(4.3b版)對纖維增強樹脂基復合材料層合板的固化過程進行數(shù)值模擬分析。理論模型包含兩部分：熱化學模型和粘彈性力學模型。熱化學模型由熱傳導和固化反應動力學組成。模型之間存在著一定的耦合關系，它們之間的耦合計算可通過COMSOL軟件中相應的應用模塊來實現(xiàn)，具體的模擬流程圖見圖1。

圖1 復合材料固化過程數(shù)值模擬程序流程圖Fig.1 Flow diagram of numerical simulation procedure during cure cycle for composite

首先，利用軟件中的“熱傳遞”應用模塊模擬熱傳導模型。在模塊中輸入材料的初始溫度、固化工藝參數(shù)及熱物理性能參數(shù)，同時調(diào)用第二個模塊(固化反應動力學模塊)中計算得到的固化度值，通過計算式(1)得到固化溫度分布。

其次，利用軟件中的“系數(shù)型偏微分方程”應用模塊模擬固化動力學模型。系數(shù)型偏微分方程的形式為：

β·u+au

(20)

式中，u和t分別為因變量和瞬態(tài)時間；c，a，ea，da，α，β和γ均為待定系數(shù)；f是源項。經(jīng)過設定系數(shù)和源項，將式(20)轉化為式(3)的形式，輸入初始固化度和固化動力學參數(shù)(見表1)，并調(diào)用由熱傳遞模塊計算的溫度分布，得到材料的固化度分布。

最后，利用軟件“結構力學”應用模塊中的“粘彈性材料”項模擬粘彈性力學模型。在該項中輸入邊界載荷、松弛因子、松弛時間(見表2)及彈性力學性能參數(shù)(見表3)，并調(diào)用由“熱傳遞”模塊和“系數(shù)型偏微分方程”模塊計算的溫度分布和固化度分布，得到層合板的彈性模量變化情況和固化殘余應力分布。

整個耦合過程在每一個時間增量內(nèi)反復迭代求解，直到達到所要求的精度，整個循環(huán)過程一直作用到固化完成。上述過程中計算得到的新的狀態(tài)參數(shù)有材料密度、比熱容、導熱系數(shù)、熱膨脹系數(shù)和模量等，在下一個時間增量內(nèi)進行耦合求解時，需要更新這些狀態(tài)參數(shù)。

4 算例驗證

為了驗證所構建模型的有效性和準確性，構造了一個正交各向異性四層復合材料層合板，鋪層順序為[0°/90°/90°/0°]，尺寸為10.16×10.16× 2.54cm，如圖2所示。材料選為碳纖維增強環(huán)氧樹脂預浸料AS4/3501-6，AS4碳纖維被視為橫向同性材料，因此，橫向各方向性能參數(shù)相同，3501-6環(huán)氧樹脂和AS4碳纖維的力學性能參數(shù)如表3所示。網(wǎng)格的劃分直接影響著數(shù)值模擬的收斂性和準確度，本研究通過有效網(wǎng)格驗證，選取了最佳網(wǎng)格，即采用六面體單元劃分網(wǎng)格，其離散的單元總數(shù)為9216，如圖3所示。

圖2 復合材料層合板幾何模型Fig.2 Schematic of the composite laminates

實驗采用的固化工藝如圖4中的短橫虛線所示，包括五個階段，在第一加熱階段，以2.5℃/min的速率從25℃升溫至 116℃，保溫1h后再次以2.5℃/min的速率加熱至177℃，隨后進入第二保溫階段，保溫2h，以-2.5℃/min的速率冷卻至25℃。在模具與構件的界面上設置了兩個邊界條件：熱量和壓力(見圖3)，在層合板的外表面設置模具加熱的固化工藝(見圖4)，在側面和底面施加壓力為689kPa。

5 分析與討論

5.1 熱化學分析

圖4 復合材料層合板中心點(5.08, 5.08, 1.27)的溫度變化情況Fig.4 Development of temperature in central point (5.08, 5.08, 1.27) of the composite laminate under manufacture recommended cure cycle

圖5 復合材料層合板中心點(5.08, 5.08, 1.27)的固化度變化情況Fig.5 Development of DOC in central point (5.08, 5.08, 1.27) of the composite laminate under manufacture recommended cure cycle

圖6 復合材料層合板表面點(5.08, 5.08, 0)到中心點(5.08, 5.08, 1.27)直線段上溫度變化情況Fig.6 Development of temperature on a line from surface point (5.08, 5.08, 0) to central point (5.08, 5.08, 1.27) of the composite laminate

圖4和圖5為固化過程中層合板中心點(5.08, 5.08, 1.27)的溫度和固化度變化情況。結合圖6中關于層合板厚度方向上溫度變化情況的數(shù)值計算結果進行固化的熱化學分析。在第一升溫階段，由于外部先受熱，而內(nèi)部樹脂尚未發(fā)生劇烈反應，所以，層合板的中心溫度低于環(huán)境溫度。在接下來的保溫階段，由于外部熱量向層合板內(nèi)部傳遞，再加上化學反應放熱的不斷積累，中心溫度迅速上升，并在43min時達到第一個峰值121℃。在第二升溫階段，由于復合材料較低的導熱系數(shù)和熱傳導滯后，層合板中心溫度再次低于環(huán)境溫度，但交聯(lián)反應更加劇烈。進入第二保溫階段后，隨著熱傳導和內(nèi)部固化放熱的積累，高溫區(qū)域再次向內(nèi)部轉移，并在129min時達到第二峰值191℃。最后，在冷卻階段，隨著固化反應基本完成，固化放熱減少，中心溫度漸漸接近室溫。從圖4和圖5可見，本研究模型關于材料固化溫度和固化度的計算結果與Kim 和 White[3]的研究結果基本一致。由于本研究模型考慮了材料熱物理性能隨固化溫度和固化度變化的特性，計算出的固化過程中層合板中心溫度在第一和第二保溫階段的峰值分別比文獻計算結果小2.3%和2.5%，固化反應開始時間提前了10min左右，而最終固化度值降低了約0.01。

5.2 力學分析

圖7和圖8分別為層合板中心點(5.08, 5.08, 1.27)厚度方向彈性模量和剪切模量隨時間的變化情況。從圖可見，彈性模量和剪切模量均隨固化度產(chǎn)生非線性變化。在第一升溫和保溫階段，樹脂還未發(fā)生劇烈反應，處于粘流態(tài)，模量很小，幾乎不發(fā)生變化；到了第二升溫階段，隨著固化反應的劇烈進行，模量迅速增大；約從180min開始，隨著固化反應減弱，模量的變化速率明顯降低，直至冷卻結束。

圖7 復合材料層合板中心點(5.08, 5.08, 1.27)厚度方向彈性模量變化情況Fig.7 Development of elastic modulus in through-thickness direction in central point (5.08, 5.08, 1.27) of the composite laminate under manufacture recommended cure cycle

圖8 復合材料層合板中心點(5.08, 5.08, 1.27)剪切模量變化情況Fig.8 Development of shear modulus in central point (5.08, 5.08, 1.27) of the composite laminate under manufacture recommended cure cycle

圖9和圖10分別為層合板在點(5.08，0，1.27)處的層間正應力σ3和0°層中點(5.08，5.08，2.22)處垂直于纖維的橫應力σ2隨時間的變化情況。從兩圖均可看出，在開始階段，樹脂還是液態(tài)，不產(chǎn)生殘余應力。之后，隨著固化反應迅速進行，由于溫度場和固化度場的不均勻性而產(chǎn)生殘余應力。在冷卻階段，殘余應力迅速增大。從圖9和圖10還可以看出，本研究模型關于應力的預測結果與White和Kim[10]的研究結果基本一致。冷卻結束時，模型計算所得的正應力為23.5MPa，比文獻計算結果小6.1%，而由本模型計算所得的橫應力為31.4MPa，比文獻計算結果小5.4%。

圖9 復合材料層合板間點(5.08，0，1.27)處正應力σ3的變化情況Fig.9 Development of interlaminar normal stress σ3 at the point (5.08, 0, 1.27) of the composite laminate under manufacture recommended cure cycle

圖10 復合材料層合板在0°層中點(5.08，5.08，2.22)處垂直于纖維的橫應力σ2的變化情況Fig.10 Development of transverse stress σ2 being perpendicular to the fiber in the 0° ply at point (5.08，5.08，2.22) of the composite laminate under manufacture recommended cure cycle

圖11為冷卻結束時，層合板的Von Minses 應力切面分布圖。如圖所示，由于層合板內(nèi)溫度梯度較大，內(nèi)部應力相比周圍區(qū)域相差較大。

圖11 復合材料層合板在固化結束時Von Minses 應力切面分布圖Fig.11 Sectional drawing of Von Minses stresses at the end of cure cycle from the two models

本研究建立了熱固性樹脂基復合材料層合板固化過程三維有限元分析模型，并考慮了材料性能參數(shù)的時變特性，實現(xiàn)了熱化學模型和粘彈性力學模型的多物理場耦合計算。

預測了AS4/3501-6復合材料層合板固化過程中溫度場、固化度場及殘余應力場的分布情況，與White 和 Kim所做研究相比，計算結果基本一致，但固化過程中層合板中心溫度峰值減小2.3%～2.5%，固化反應開始時間提前約10min，最終固化程度下降0.01，冷卻結束后的正應力減少1.4MPa，橫應力減少1.7MPa。

本研究的研究方法具有較強的適用性，對優(yōu)化固化工藝參數(shù)，提高固化質(zhì)量具有一定的指導意義。