复合材料层合板在固化阶段，不仅层合板外部环境温度变化，而且层合板内部固化反应产生化学放热，两者相互影响，是一个复杂的耦合过程，导致层合板内部产生复杂的温度梯度。这种不均匀分布的温度场和固化度场，将在复合材料内产生热应力和变形，是复合材料制品中产生诸如尺寸误差、翘曲等缺陷的根本原因。因此，研究复合材料层合板在固化工艺过程中的温度和固化度分布及其变化规律，对提高复合材料工艺质量有重要意义[1]。
近年来对树脂基复合材料固化过程的有限元模拟取得了一定进展。Loos等[2]采用模块化方法研究层合板固化过程，针对AS4/3501-6单向层合板建立了一维数值模型；Bogetti等[3]采用二维有限元方法数值模拟了任意截面形状和边界条件层合板的固化过程；杨正林等[4]采用有限元与有限差分相结合的方法分析了二维层合板模型的固化变形；Cheung、张纪奎等[5, 6]用有限元方法研究了三维条件下的固化过程。
由于固化过程中温度场和固化度场耦合问题的复杂性，研究者基本上都是自行编制程序来实现数值仿真，增加了研究的难度和周期，而且程序的适用性也受到限制。本文从正交各向异性材料的三维热传导方程和固化动力学方程出发，基于有限元软件ABAQUS进行二次开发，实现了树脂基复合材料固化过程的三维数值模拟，通过对2种材料体系的固化过程进行了算例验证，证明了程序的正确性，并讨论了层合板边界条件、厚度、升温速率对温度场和固化度场的影响。
本文研究所涉及的玻璃纤维/聚酯树脂的固化动力学方程属于n级反应模型[7]，AS4/3501-6的固化动力学方程为自催化反应模型[8]，当采用不同固化动力学方程时，只需依据利用Fortran语言在ABAQUS子程序的基础上进行编辑即可完成，因而具有更强的适用性，缩短了程序的编写周期。可用于分析复杂形状和边界条件下结构固化过程的温度场和固化度场，具有较强的适用性，对改善复合材料结构固化工艺参数具有指导意义。
1 热传导控制方程和固化动力学方程热固性树脂基复合材料的固化过程是一个热与化学反应相互耦合的过程。复合材料构件内部的温度分布由向复合材料的传热速率和固化反应产热速率共同决定，复合材料固化温度场的分析本质上是一个具有非线性内热源的热传导问题，其中内热源是树脂基体固化反应放出的热量。通过对该热传导问题的求解，可以得到复合材料在固化过程中任意时刻、位置的温度及固化度。由于固化阶段树脂基本不发生流动，可忽略对流传热影响，则根据Fourier热传导定律和能量平衡原理建立该问题的热传导控制方程：
树脂的固化反应决定了热传导控制方程中内热源的热量大小，综合各类文献的研究方法，对固化动力学模型的表征方法主要有2种：微观水平(力学的)和宏观水平(唯象的)。其中唯象模型以化学反应动力学为主要特征，忽略各组分之间相互作用的细节，也不关注树脂的组成或配方，是最被广泛使用的方法。
大多数唯象模型以方程(3)、(4)为基础：
根据树脂反应机理函数f(α)的形式，可以将现有模型分成n级反应模型和自催化模型2类，方程(5)和(6)的化学动力学模型分别属于n级反应模型和自催化反应模型。
ABAQUS具有强大的非线性分析能力，能很好地解决热-力耦合问题。本文主要运用ABAQUS/Standard提供的HETVAL、USDFLD用户子程序，而在对流换热条件和温度边界条件比较复杂的模型进行定义时，使用了FILM和DISP子程序。
固化过程中树脂固化放热的模拟是通过子程序HETVAL定义热源项来实现的。该子程序用来在传热分析中定义内部生热与由于产热导致的热流，允许内部生热依赖于状态变量；子程序USDFLD对STATEV的任何更新都可传递到HETVAL中。
表征化学反应程度的固化度场通过用户定义场在材料属性中进行定义，定义中用使用了用户子程序USDFLD。状态变量可以在USDFLD中更新，然后传递到子程序HETVAL中。
温度场和固化度场的耦合求解如图 1所示。
|
| 图 1 温度场与固化度场耦合求解流程图 |
在子程序HETVAL中，通过固化动力学方程得到固化速率和固化度，并存储在状态数组STATEV中，更新后的STATEV数组传递到子程序USDFLD以得到固化度场；固化反应放热通过FLUX模块计算，产生的热量传递到ABAQUS热分析模块中得到新的温度场，并利用新温度计算得到固化速率和固化度，如此循环，直到固化过程结束。
3 算例及程序验证根据上述温度场和固化度场的耦合求解流程，基于用户子程序对ABAQUS进行二次开发，编写了用于计算树脂基复合材料固化过程的子程序。本文采用了2个算例对不同固化动力学模型进行了验证。
3.1 算例1参考文献[5]，采用[0/90]铺层的玻璃纤维/聚酯树脂层合板，长宽均为15.24 cm，厚为2.54 cm，共42层。玻璃纤维/聚酯树脂固化动力学经验公式为：
| ρ/(kg·m-3) | C/(J·kg·℃-1) | k11/(W·m-1·℃) | k22=k33/(W· m-1·℃-1) |
| 1 890 | 1 260 | 0.432 6 | 0.216 3 |
边界条件用对流换热系数h与导热系数k11的比值表示。上表面h/k11=87 m-1，下表面h/ k11=125 m-1，四周绝热。计算得到该层合板中心点温度和固化度历程，结果如图 2a)和图 2b)所示。由结果可知，温度和固化度均与实验结果符合良好，从而验证了程序针对n级动力学模型的正确性。
|
| 图 2 玻璃纤维/聚酯树脂层合板中心点温度、固化度历程 |
参考文献[5]，算例2采用[0/90]铺层的AS4/3501-6层合板，几何模型和算例1相同，长宽均为15.24 cm，厚为2.54 cm，共42层。AS4/3501-6固化动力学经验公式如下：
| ρ/(kg·m-3) | C/(J·kg·℃-1) | k11/(W·m-1·℃) | k22=k33/(W·m-1·℃-1) |
| 1 578 | 862 | 12.83 | 0.413 5 |
| A1/(min-1) | A2/(min-1) | A3/(min-1) | ΔE1/(min-1) |
| 2.102×109 | -2.014×109 | 1.96×105 | 8.07×104 |
| ΔE2/(min-1) | ΔE3/(min-1) | R/(J·(mol·K-1) | Hr/(kJ·kg-1) |
| 7.78×104 | 5.66×104 | 8.314 | 198.9 |
给定边界条件与算例1中相同。图 3a)和图 3b)分别为该层合板中心点温度和固化度历程，与文献实验结果进行比较，可以看出符合性较好，验证了程序针对自催化模型时的正确性。
|
| 图 3 AS4/3501-6层合板中心点温度、固化度历程 |
由计算结果可知，边界条件、厚度、升温速率对层合板温度场和固化度场均有影响。以下讨论中，除了研究厚度对温度场和固化度场的影响时，试件厚度不同外，在研究边界条件、升温速率对温度场和固化度场的影响时，层合板厚度均为2.54 cm；中心点指的是中间层(z=t/2)处。
4.1 边界条件的影响热传导问题有3类边界条件：①Dirichlet边界条件，即已知温度的边界条件;②Neuman边界条件，即已知热流密度的边界条件；③Robin边界条件，即已知对流换热系数的边界条件。第3类边界条件中，当对流换热系数h→∞时，边界条件就转化为第1类边界条件；当h→0时，边界条件就转化为第2类边界条件。故可通过改变h的大小来改变边界条件和特征。本文用h/k的比值变化来表示边界条件的改变，其中k为边界法向方向的热传导数，分别取h/k=100 m-1和500 m-1作为层合板上、下表面的边界条件，侧面仍为绝热边界。
图 4a)给出了层合板表层(z=0)中心点和中间层(z=t/2)的中心点温度曲线，由计算结果可知：随着h/k比值的增大，层合板升温速率加快，中心点温度峰值逐渐减小。图 4b)是对应条件下的固化度曲线，结果表明h/k越小，放热期间固化度梯度将会增加。
|
| 图 4 边界条件对AS4/3501-6层合板温度、固化度的影响 |
对于前述AS4/3501-6层合板，取1.27 cm、2.54 cm、3.81 cm 3种不同厚度，上下表面为第1类边界条件，侧面为绝热边界。图 5a)和图 5b)为这3种不同厚度层合板固化过程中心点的温度和固化度变化历程。由图 5a)的温度曲线可知，在中心点温度超过保温平台温度之前，1.27 cm层合板的温度最高，2.54 cm板次之，3.81 cm板温度最小，超过保温平台温度最晚，但温度峰值明显高于其他2个。由图 5b)的固化度曲线可知，1.27 cm层合板固化反应启动时间最早，进行时间最长，而3.81 cm层合板固化反应启动最晚，但最先结束。所以厚度越大，中心点温度峰值就出现得越晚，峰值越大，中心点固化开始得就越晚，而完成固化所需时间就越短。
|
| 图 5 厚度对AS4/3501-6层合板温度、固化度的影响 |
对于前面所述AS4/3501-6层合板，上下表面采用第1类边界条件，侧面为绝热边界，层合板在389 K温度下固化60 min后，再分别施以1.5 K/min，3 K/min和6 K/min的加热速率达到449 K。
图 6a)和图 6b)分别给出了不同升温速率下，层板中心点温度和固化度随时间的变化曲线。由图 6a)可知，升温速率越慢，中心点温度升温越慢，峰值越小，峰值过后的降温过程也越慢；同时由图 6b)可知，升温速率越慢，固化速率越慢，固化完成所需时间也越长。
|
| 图 6 升温速率对AS4/3501-6层合板温度、固化度的影响 |
1)根据热传导和唯象固化动力学理论，建立了正交各向异性复合材料层合板固化过程的三维有限元模型；研究了基于ABAQUS子程序耦合求解温度场和固化度场的方法，通过对比算例与文献实验结果，验证了本文方法的正确性。
2)由边界条件对温度场和固化度场的影响表明，h/k越大，层合板升温越快，中心点温度峰值越小，放热期间固化度梯度减小，因而h/k的增大有利于减轻温度场和固化度场的不均匀性。
3)由厚度对温度场和固化度场的影响表明，层板越厚，固化开始得越晚，中心点温度超过保温平台温度越晚，温度峰值越大，完成固化所需时间越短。
4)升温速率对温度场和固化度场的影响表明，升温速率越慢，中心点温度升温越慢，峰值越小，降温过程也越慢，从而降低了固化速率，使得固化完成所需时间延长，因而适当降低升温速率有利于减轻温度场和固化度场的不均匀性，但同时也会延长固化完成时间，二者之间要根据实际情况权衡。
5)本文基于ABAQUS子程序所建立的复合材料三维有限元固化数值模拟方法，适用于模拟具有不同类型动力学方程的复合材料固化过程，可分析复杂形状和边界条件下结构固化过程温度场和固化度场，因此对改善固化工艺参数具有一定的指导意义。
| [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 |
| Click to display the text | |
| [2] | Loos A C, Springer G S. Curing of Epoxy Matrix Composites[J]. Journal of Composite Materials,1983,17(2): 135-169 |
| Click to display the text | |
| [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 |
| Click to display the text | |
| [4] | 杨正林,陈浩然. 层合板在固化过程中瞬态温度场及固化度的有限元分析[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) |
| Cited By in Cnki (21) | Click to display the text | |
| [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 |
| Click to display the text | |
| [6] | 张纪奎, 郦正能, 关志东,等. 热固性复合材料固化过程三维有限元模拟和变形预测[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) |
| Cited By in Cnki (31) | Click to display the text | |
| [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 |
| Click to display the text | |
| [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 |
| Click to display the text |
