树脂基复合材料曲面结构件固化变形数值模拟

1、树脂基复合材料曲面结构件固化变形数值模拟马云荣*, 贺继林, 李 栋, 谭 耀, 徐 雷(中南大学 机电工程学院, 长沙 410083)摘要：研究分析了碳纤维增强树脂基复合材料在固化过程中密度、模量、热胀系数、比热、热传导系数等材料物性的变化，并将上述变化引入到数值模拟当中。针对复合材料复杂曲面件，提出了利用定常流动的流线方程构建曲线坐标系的新方法。根据建立的曲线坐标系，运用有限单元法计算了某轻型飞机机翼上蒙皮板在固化过程中的内部温度、固化度和内应力的分布情况以及材料物性随固化度的变化情况；计算了由于内部不均匀的温度场和固化度场、热膨胀系数的各向异性和固化引起的树脂体积收缩而导致的结构变形。研

2、究表明：引入材料物性变化使固化过程更加合理、结果更加精确，利用定常流动的流线方程构建的曲线坐标系完全适用于复合材料曲面件的有限元仿真。关键词：复合材料；曲面构件；曲线坐标系；有限元分析；内应力；固化变形中图分类号：TB330.1 文献标志码：ANumerical simulation of curing deformation of resin matrix composite curved structureMA yunrong*,HE jilin,LI dong,TAN yao,XU lei(College of mechanical and electrical engineering

3、,Central South University,Changsha 410083,China)Abstract: The evolutions of physical properties during the curing process of the carbon fibre reinforced resin composite, such as density, modulus, coefficient of thermal expansion, specific heat capacity and thermal conductivity were analyzed and thos

4、e evolutions were introduced into the numerical simulation. The new approach of constructing the curvilinear coordinate system by streamline equation of steady flow was proposed with the complex curved structure composite as the study object. The distribution of internal temperature, degree of cure,

5、 internal stresses and evolutions of physical properties during the curing process of the composite skin plate of a light aircraft wing were calculated by finite-element method and the curvilinear coordinate system, then the deformation of the skin plate caused by the uneven distributions of tempera

6、ture field and cure degree field, anisotropic thermal expansion coefficients and volumetric shrinkage of resin were calculated by the same method. From the numerical simulation results, it can be concluded that the curing process will be more reasonable and more accurate when the evolutions of physi

7、cal properties of composites are adopted and the curvilinear coordinate system constructed by streamline equation of steady flow is fully applicable to finite element analysis of composite curved structure.Keywords: composites; curved structure; curvilinear coordinate system ; finite element analysi

8、s (FEA); internal stress；curing deformation收稿日期：2014-06-16；录用日期：2014-09-18；网络出版时间：网络出版地址：DOI:基金项目：国家973计划（2014CB046500）通讯作者：马云荣，硕士研究生，研究方向为高性能复合材料的制造工艺及失效模式、通用航空器的开发等 E-mail：thrashead引用格式：马云荣，贺继林，李栋，等树脂基复合材料曲面结构件固化变形数值模拟J复合材料学报, 2015, 32(x): xxx-xxx. MA yunrong, HE jilin, LI dong, et al. Numerical s

9、imulation of curing deformation of resin matrix composite curved structure J. Acta Materiae Compositae Sinica, 2015, 32(x):xxx-xxx.树脂基复合材料固化阶段是树脂发生化学反应、不断固化成型的过程。随着交联反应的进行，树脂体系会不断放热、由低弹性模量的粘性流体到高弹性模量的粘弹性直至弹性固体。在这个过程中，树脂会产生固化收缩，由于复合材料内部温度场和固化度场的不均匀会导致不同程度的热膨胀和固化收缩。这就是内部残余应力的产生和固化变形的直接原因。严重的甚至可能导致基体开裂

10、和分层现象。而为了消除固化变形而补偿性加工模具型面的传统方法既费时又浪费材料，而且具有严重的工艺和材料依赖性。问题解决的关键就在于建立一套完整的固化变形的分析和预测方法，从本质上改变基于经验和工艺试验的传统工作模式。因此，研究树脂基复合材料在固化过程中的化学物理变化、内应力场和固化变形具有重要的意义。大多数学者采用傅里叶热传导控制方程描述复合材料固化过程中的热传导。Costa等1采用了二维平面模型描述了复杂轮廓型面复合材料结构件内部的热传导，得到了温度场分布。谭华等2采用三维模型模拟了形状复杂的大型复合材料结构件内部的热传导现象，得到了内部三维温度场分布。Michaud等3研究了411-C50

11、树脂体系在固化过程中的密度和比热变化。Chen等4研究了3501-6环氧树脂在固化过程中的比热和热传导率的变化。Bogetti和Gillespie5,6对厚板复合材料的制造过程进行了研究。Prasatya等7通过模拟计算得到环氧树脂的化学收缩率对残余应力的影响占30%以上。很多学者在以往的研究中对固化过程中复合材料的等效密度、比热、热膨胀系数和热传导系数等物性参数取常量或平均值，这是不符合实际的。本文的研究中引入了固化过程中以上参数随温度和固化度的变化函数。本文作者主要研究了树脂基复合材料固化过程中的温度场和度化度场的分布、固化过程中树脂关键参数演变规律及其对残余应力和固化变形的影响；研究了树

12、脂的收缩应变过程、复合材料体积应变过程和残余应力分布；提出了建立空间曲线坐标系的新方法，并利用此方法对某型轻型飞机机翼上蒙皮进行了固化残余应力和固化变形的数值建模计算。1 树脂基复合材料固化的理论模型1.1热-化学模型反应热的温度场问题本质上是具有非线性内热源的热传导问题，其中内热源是树脂基体的固化反应放出的热量。由傅里叶热传导定律和固化动力学方程可以得到8,9：,Tc,TTt=xikijTxj+,THuddt (i,j=1,2,3) (1)边界条件：keffTn+heffTs-Tm=0 (2)式中，T为绝对温度，、c和kij分别为复合材料的等效密度、比热和热传导系数。在(1)式中，最后一项为

13、固化反应热项，Hu为固化反应完成时树脂放出的总热量。d/dt为树脂的瞬时固化率。为树脂的固化度。Ts为复合材料结构件表面的温度，Tm为模具温度，keff和heff分别为复合材料表面的等效热传导系数和对流换热系数。1.2固化动力学模型本文中所研究的复合材料为AS4/3501-6，而3501-6环氧树脂符合一种特殊的固化动力学模型5,10,11：ddt=K1+K21-0.47- (0.3)K31- (>0.3) (3)ki=Aie-EiRT (i=1,2,3) (4)其中，ki为反应速率常数，Ai是频率因子， Ei为活化能，R为通用气体常数。1.3 残余应力模型残余应力模型基于温度场和固化度

14、场来计算复合材料构件内部的热应力和固化收缩应力，从而计算固化变形。1.3.1本构方程包含温度和固化度影响的本构方程为：ijt=-tQijkl,T,t-kl- kltcdt (i,j,k,l=1,2,3) (7)其中:ij为应力分量，kl为总应变，kltc为热化学应变，Qijkl为材料刚度矩阵。1.3.2 热-化学应变纯树脂的热应变：mt=mT (8)其中m为树脂的热膨胀系数，T为温度差。复合材料纤维的热应变：fit=fiT (i=1,2,3) (9)其中fi为复合材料纤维三个方向上等效的热膨胀系数。树脂的化学收缩应变：mc=31+v-1 (10)其中, v为树脂的体积收缩率。由于复合材料纤维在

15、固化过程中没有化学变化，故其体积收缩率为0。复合材料的热-化学应变是热应变与化学应变的简单求和：kitc=kit+kic (i=1,2,3) (11)其中，kit为复合材料各方向的热应变，kic为复合材料各方向的化学应变。2 固化过程中复合材料的热物理性能变化2.1 模量树脂弹性模量可由下式确定5：Em=Em0 (<get)Em=1-modEm0+modEm (get)(12)mod=-get1-get (13)式中，Em为树脂弹性模量，Em0为树脂未固化时的弹性模量，Em为树脂固化完全时的弹性模量，get为树脂凝胶点。复合材料的弹性模量可由细观力学公式得到，在此不再赘述。2.2 密度在

16、本文的模拟中，密度由下式计算12：m=0.09+1.232 0.451.272 >0.45 (g/cm3)(14)假设复合材料纤维密度不变，则可以用混合率公式计算复合材料密度：=Vff+(1-Vf)m (15)式中，m为树脂的密度，f为纤维的密度，为复合材料密度，Vf是纤维的体积分数。2.3 热膨胀系数对树脂而言，热膨胀系数可以根据树脂在固化过程中的三种不同状态来计算13：m=ml <get and TTg()mr get and TTg() mg T<Tg() (16)式中，ml、mr、mg分别为树脂在粘流态、橡胶态、玻璃态的热膨胀系数。其中，橡胶态的热膨胀系数可以取玻璃态

17、热膨胀系数的2.5倍13。由于复合材料纤维的加入，复合材料的热胀系数也为各向异性。单向复合材料沿纤维方向和垂直于纤维方向的热膨胀系数计算公式为5,14：1=1fEfVf+mEmVmEfVf+EmVm (17)2=2f+v12f1fVf+m+vmmVm- v12fVf+vmVm1fEfVf+mEmVmEfVf+EmVm (18)其中1、2分别是复合材料沿纤维和垂直于纤维方向的热膨胀系数，1f、2f分别是平行于纤维方向和垂直于纤维方向的纤维热胀系数，Vm是树脂的体积分数。2.4 比热对于3501-6树脂，比热是温度与固化度的线性函数12,15：Cm=4.1840.468+5.975×10

18、-4T-0.141 (J/ g) (19)对于碳纤维，比热与温度呈线性关系16：Cf=1.390+4.50×10-3T (20)复合材料比热可由混合率计算14：c=mVmCm+fVfCf (21)式中，c、Cm、Cf分别是复合材料、树脂、纤维的比热。2.5 热传导系数对于3501-6树脂，热传导系数是温度与固化度的函数12,16：Km=0.041843.85+0.035T-0.41 W/m (22)对于碳纤维，热传导系数是温度的函数16： Kf=0.742+9.02×10-4T (23)复合材料的热传导系数也具有各向异性，复合材料垂直于纤维方向和平行于纤维方向的热传导系数为

19、17：KTKm=1-2Vf+1B-41-B2Vftan-11-B2Vf1+BVf (24)B=2(KmKf-1) (25)KL=KmVm+KfVf (26)式中，KT、KL分别为复合材料垂直于纤维方向和平行于纤维方向的热传导系数。3 三维实体复杂曲面构件的曲线坐标系反映材料各向异性的介质参量为三维二阶张量，且物理规律本身常常要求这些张量为对称张量。通过选择合适的曲线坐标系，可使得它们对角化，从而只出现3个对角元素，由它们完全确定相应张量的性质。以复合材料为例，它在平行和垂直于纤维方向的物理性质截然不同，相应出现的热传导系数、热膨胀系数、模量等均为张量。在这种情况下，可直接将纤维选为一条坐标曲线

20、，并选择其余两个坐标，共同构成正交曲线坐标系（其3个基矢彼此正交）。该坐标系中，上述张量将转化为对角张量，仅含3个对角元素。假设在三维实体曲面构件中有理想流体在做定常流动，则其流线与迹线重合。在复合材料曲面构件中，完全可以认为纤维的主方向与定常流动的流线方向一致。则可利用流线方程来构造曲线坐标系。利用欧拉运动微分方程可以求得流线方程，用以下边界条件对流线方程进行稳态求解可以得到三维实体曲面构件内的稳态流线方程：流动条件：v=u|u| -pI+u+uT=0u=0 (27)壁： un=0 (28)入口： un=un (29)出口： p=0 u+uTn=0 (30)式中，v为流线方向向量，u为速度矢

21、量，p为压强，I为横截面积，n为曲面法向量。则，复合材料曲线坐标系的第一基矢e1可以写成：e1=v|v| (31)取a为曲面第二坐标曲线的基矢向量，由下式构造复合材料曲线坐标系的第二基矢e2：e2=a-(ae1)e1|a-(ae1)e1| (32)容易看出，复合材料曲线坐标系的第三基矢e3为：e3=e1×e2 (33)设(xi,yi,zi)为三维空间ei的直角坐标，则曲线坐标系可表示为：ei=ei(xi,yi,zi)， i=1,2,3 (34)很容易证明，e1、e2与e3两两正交。由基矢的正交性立即可推得上述张量为对角张量，即构造出了复合材料曲面实体构件的正交曲线坐标系。4 某轻型飞

22、机机翼蒙皮板固化过程数值模拟模型由四个模块耦合计算得出结果：热传导模块、热变形模块、偏微分模块和曲线坐标系求解模块。热传导模块计算了蒙皮板内的热传导与温度分布；偏微分模块根据温度分布计算出蒙皮板内的固化度分布；热变形模块根据上述两个模块的计算结果计算出应力应变分布，最后得出变形量；曲线坐标系求解模块根据流动方法和边界条件求解出曲线坐标系，定义上述三个模块中的参数矢量方向。上述模块的计算方程和边界条件在上文中均已给出，在此不再赘述。利用上述理论，本文中对某公司的某型轻型飞机机翼蒙皮板进行了固化变形研究。蒙皮板采用真空注射工艺，热压罐固化成型。材料为AS4/3501-6，板厚0.008m，铺层类型

23、为0/45/-45/90，纤维体积含量为50%。由于复合材料在模具中约束条件较复杂，这里将约束条件简化为自由约束情况。AS4碳纤维和3501-6树脂材料性能参数分别见表1和表2。表1 AS4碳纤维材料性能参数Table 1 Material properties of AS4 carbon fiber E1/MPaE2/MPav12v13v23G12/MPaG13/MPaG23/MPa1/(1/)2/(1/)f/g/cm32.068×1051.72×102.758×1042.758×1046.894×103-9.0×

24、;10-77.2×10-61.78Note: E1,E2Elastic modulus; v12,v13,v23Poissons ratio; G12,G13, G23Shear moudulus; 1, 2Thermal expansion coefficient; fDensity表2 3501-6树脂材料性能参数Table 2 Material properties of 3501-6 resinEm0/MPaEm/MPavmgetmr1/mg1/v3.4473.447×1030.350.571.44×10-45.76×10-50.05A1/min

25、-1A2/min-1A3/min-1E1/(J/mol)E2/(J/mol)E3/(J/mol)Hu/(J/kg)2.102×109-2.014×1091.96×1058.07×1047.78×1045.66×104198900Note: Em0,EmElastic modulus；vmPoissons ratio; getgel point; mr,mgThermal expansion coefficient; vvolumetric shrinkage; A1,A2,A3Frequency factor; E1,E2,E3Act

26、ivation energy; HuHeat of reaction固化工艺温度曲线如图1所示。稳态求解机翼蒙皮板中流线（即0°角复合材料纤维主方向）结果如图2所示。图1AS4/3501-6复合材料真空注射工艺固化温度曲线Fig.1 Temperature curve of vacuum injection process of AS4/3501-6 composite curing图2 AS4/3501-6复合材料蒙皮板流线（即纤维主方向）Fig.2 Streamline of AS4/3501-6 composite skin plate根据上述理论计算得出的机翼蒙皮板复合材料的

27、局部坐标系如图3、图4（以0度铺层和45度角铺层为例）所示。(a) 0o(b) 45o图3 AS4/3501-6复合材料蒙皮板0o和45o角铺层材料坐标系Fig.3 Material coordinate system of 0° and 45o ply of AS4/3501-6 composite skin plate蒙皮板中心点与外表面温度随时间变化曲线如图4所示。 图4 AS4/3501-6复合材料蒙皮板中心点与外表面温度曲线Fig.4 Center point and surface temperature curve of AS4/3501-6 composite ski

28、n plate由于蒙皮板厚度很小，中心点温度与外表面温度基本一致。但还是可以看出，在2000s和8000s左右，中心点温度略高于外表面温度，也就是固化工艺温度。这是树脂固化放热的结果。蒙皮板中心点固化度随时间变化曲线如图5所示。蒙皮板中心点树脂杨氏模量和剪切模量随时间变化曲线如图6所示。图5 AS4/3501-6复合材料蒙皮板中心点固化度曲线Fig.5 Degree of cure curve of AS4/3501-6 composite skin plate center point图6 蒙皮板中心点3501-6树脂杨氏模量和剪切模量随时间的变化曲线Fig.6 Young's mo

29、dulus and shear modulus of resin with time for AS4/3501-6 composite skin plate center point从图中可以看出，第一阶段，树脂处于粘流态，模量很小，几乎可以忽略不计；第二阶段，模量随固化时间增加而增加，并且增加的速率越来越小；第三阶段，模量趋于稳定。蒙皮中心点复合材料体积应变与体积应力（残余应力）随时间变化曲线如图7所示。(a)(b)图7 AS4/3501-6复合材料蒙皮板中心点体积应变和应力随时间的变化曲线Fig.7 Volumetric strain and volumetric stress curve

30、 with time for AS4/3501-6 composite skin plate center point由图7可知，在0s到8000s左右，树脂为液态，复合材料体积应变几乎为零，8000s到12000s左右，体积应变主要变现为固化收缩，在固化基本完成后，也就是12000s左右之后，随着温度的下降，体积应变主要是热收缩。同样的道理，在0s到8000s左右，树脂为液态，不会产生内应力。8000s左右到18600s，由于不均匀的温度场和固化度场产生内应力。蒙皮板末端点z轴方向位移与变形云图随时间变化曲线如图8、图9所示.图8 AS4/3501-6复合材料蒙皮板最末端点z轴方向位移曲线F

31、ig.8 Displacement of z-axis direction curve of AS4/3501-6 composite skin plate center point图9 AS4/3501-6复合材料蒙皮板Z轴方向位移云图（固化完成）Fig.9 Displacement nephogram of z-Axis direction of AS4/3501-6 composite skin plate center point when curing is complete经上述计算，得到蒙皮板末端最大翘曲位移为23.151mm。根据模具型面补偿技术报告，实际翘曲位移为22 mm。

32、数值计算误差为5.23 %。5 结 论1.在树脂基复合材料固化过程中引入材料物性参数变量的函数定义使固化过程更加合理、结果更加精确。2.利用理想流体定常流动的流线方程建立的曲线坐标系完全适用于复杂曲面复合材料构件的三维实体有限元分析中的材料方向定义。3.本文中提出的流线方程构造曲线坐标系的方法可以应用于其他各向异性的曲面件三维实体有限元仿真的曲线坐标系定义中。参考文献：1 Costa V, Sousa A. Modeling of flow and thermo-kinetics during the cure of thick laminated composites J.Internati

33、onal Journal of Thermal Sciences, 2003, 42 (1):15-22.2 谭华, 晏石林. 热固性树脂基复合材料固化过程的三维数值模拟J. 复合材料学报, 2004,21(6): 167-172. Tan Hua , Yan Shilin. Three-dimensional simulation of curing process for thermoset compositesJ. Acta Materiae Compositae Sinica, 2004,21(6):167-172.3 Michaud D, Beris A. Dhurjati P Th

34、ick-sectioned RTM composite manufacturing, part II. Robust cure cycle optimization and control J. Journal of Composite Materials, 2002, 36(10): 1201-1230.4 Chern BC, Moon TJ, Howell JR, et al. New experimental data for enthalpy of reaction and temperature-and degree-of-cure-dependent specific heat a

35、nd thermal conductivity of the Hercules 3501-6 epoxy system J. Journal of Composite Materials, 2002, 36 (17): 2061-2072.5 Bogetti T A, Gillespie J W. Process-induced stress and deformation in thick-section thermoset composite laminates J. Journal of Composite Materials, 1992, 26(5):626-660.6 Bogetti

36、 T A, Gillespie J W. Two-dimensional cure simulation of thick thermosetting composites J. Journal of Composite Materials, 1991, 25(3): 239-273.7 Prasatya R McKenna G B, Simon S L. A viscoelastic model for predicting isotropic residual stresses in thermosetting materials: effects of processing parame

37、ters J. Journal of Composite Materials, 2001, 35(10): 826-849.8 Yan X. Finite element modeling of curing of epoxy matrix composites J. Journal of Applied Polymer Science, 2007,103 (4): 2310-2319.9 Svanberg J M, Holmberg J A. Prediction of shape distortions Part I. FE-implementation of a path dependent constitutive modelJ. Composites Part A: Applied Science and Manufacturing, 2004,35(6): 711-721