|
复合材料结构固化过程中，结构内部会产生固化残余应力，构件脱模后固化残余应力释放，结构必然产生固化变形，尤其是对于大型复合材料整体结构，固化后的变形将导致一系列的装配问题和使用问题。因此，如何能够有效地预测复杂结构固化后的变形是目前复合材料结构制造过程中面临的重要的问题。 本期文章特邀复合材料力学公众号专家库中长期从事复合材料固化过程仿真的老师来讲述固化温度场及固化变形预测方面的分析方法。 欢迎感兴趣的读者在留言区或者后台沟通交流。同时，也欢迎关注近期即将举办的2023第一期复合材料固化仿真分析培训。 培训课程推荐 目前热固性树脂基复合材料固化仿真分析方向，研究的人已经有所减少，这从去年珠海的复材大会中，该领域极少的PPT报告以及近些年的文章发表量都可以看出。但我认为还是有很多问题没有研究明白。从我自己了解的情况来看，几大飞机制造主机厂中西飞、成飞以及昌飞，在复材构件的固化变形这一问题上，实际上并没有很好的解决方法，其他诸如沈飞、哈飞，笔者没有去考察过，不敢妄下结论。从现有的研究文献来看，似乎我们对简单的诸如平板或者L型、U型构件，已经掌握了这些构件的固化变形规律与研究方法。但是关于实际使用中的构件，例如机翼壁板或者长桁，研究文献并不多，而且有些文献的计算结果或者仿真分析结果很难复现（也可能是笔者自身的学术水平不足）。所以这也是我想投稿的原因，希望可以抛砖引玉，和大家一起探讨该领域目前的一些问题。 热固性树脂基复合材料（以下简称复合材料）的固化仿真，现有的主流研究是把该问题转化为热—力问题，通过研究该热—力问题，研究复合材料内部温度、固化度、残余应力等参量的演化规律，继而求得固化变形。热方面，主要是要考虑树脂的固化反应放热，以及材料热力学参数，诸如比热、导热系数随温度、固化度的变化。再深入一些，要考虑树脂流动的影响，这对于厚截面复合材料的温度场计算，影响还是很大的。力方面，主要是材料本构模型的问题，主要是针对树脂，由最初的线弹性发展到今天的粘弹性。可参考之前公众号中忽忽沛幺同志做过的一个题为“基于Abaqus的复合材料固化变形及参与应力仿真简介。  基于Abaqus的复合材料固化变形及残余应力仿真简介 1 复合材料固化过程温度场研究 复合材料固化过程内部温度场分布，一方面影响构件成型后内部固化残余应力的分布情况，另一方面过高的温度峰值以及固化不均匀对构件成型后的力学性能也要很大影响。针对固化温度场的研究，现阶段主要使用的子程序有HETVAL、FILM、DISP以及USDFLD等，在庄茁先生的著作中有部分源代码，忽忽沛幺同志也进行过介绍，在此不再赘述。 我想要强调的是，HETVAL程序最大的问题，就是无法定义随温度和固化度变化的导热系数和比热容。实际上，复合材料在固化过程中，其材料性能参数并不是固定不变的，大多数情况下，与内部温度、固化度的变化密切相关。此时，UMATHT子程序便可以解决导热系数以及比热容随温度、固化度的变化问题，这是HETVAL程序所完成不了的。 如果大家对UMATHT子程序感兴趣，回头我可以再专门写一篇HETVAL和UMATHT子程序的比较分析，二者在处理树脂固化反应放热上是一致的，但UMATHT可以定义导热以及比热随温度、固化度的变化。 1.1 仿真模型 结构模型来源于文献[1]，如图1-1所示，材料为AS4/8552，固化后单层厚度0.2mm，构件铺层为[02/902]s，文献[1]中并没有进行相应温度场和固化度场的计算，在此我查了资料，找到了相应材料的热力学参数，进行温度场的仿真分析。这里我是利用HETVAL子程序进行树脂固化反应放热计算的。复合材料传热的问题，本质上可以看作是一个具有非线性内热源的热传导问题。 考虑到树脂基体的化学反应放热，并忽略树脂流动过程产生的热量传递，则复合材料制件固化过程中三维热传导方程可以用式（1-1）和（1-2）表示。在HETVAL程序中，FLUX(1)用于定义当前节点上的热流密度，即树脂的内部反应放热Q，见式（1-1）。STATEV一方面用于定义固化率（固化度的变化率，见式（1-2），以表征材料的固化动力学方程；另一方面用作子程序间传递的媒介。譬如常见的STATEV(1)表示固化度，即树脂固化程度；STATEV(2)表示固化速率。对于材料AS4/8552，其固化速率可以用公式（1-3）进行表示[2]。在温度场计算过程中，可以进一步计算出树脂玻璃化温度变化曲线，如公式（1-4）所示[2]，该曲线对于后续进行固化残余应力分析尤其重要。
式中：λx、λy、λz分别为全局坐标系下沿x，y，z三个方向复合材料制件的热传导系数；T为温度；Q为热生成率；ρc为复合材料密度；Cc为复合材料比热容；t为时间。详细参数的含义在此不细说，上述提到的论文中大家都可以找到。具体的取值见表1-1所示。温度场计算的边界条件为构件外表面均采用70W/m2·K的对流换热系数[2]，以此来表征热压罐内高温气体与复合材料构件之间的换热。固化过程中，初始温度和固化最终温度均为20℃，第一个保温平台温度为110℃，保温时长60min，第二个保温平台温度为180℃，保温时长120min，升降温速率均为2℃/min。固化工艺曲线和对流换热系数可以通过FILM子程序进行定义，也可以利用Abaqus的“Interaction”模块直接定义，如图1-2所示。
图1-1 结构模型[1]
图1-2 定义对流换热系数 表1 -1AS4/8552的热力学参数[2]
1.2 计算结果在这里需要说明两个问题。第一个问题是，我们在用实体单元做复合材料分析的时候，通常是采用Abaqus里面的复材单元模块进行铺层以及区域的选择（“Property”—“Creat Composite Layup”），而帮助文档里面明确说明“Continuum shell and solid composite layups are expected to have a single element through the entire thickness across a combination of all the regions specified in the composite layup. Each single element through the thickness contains the multiple plies that you defined in the ply table. If the region to which you assign your continuum shell or solid composite layup contains multiple elements through the thickness, each element will contain all of the plies that you defined in the ply table and the analysis results will not be as expected. If your model contains multiple continuum shell or solid elements through the thickness of a region, you can obtain correct results by defining a separate composite layup for each layer of elements. You can define a composite layup for each layer by selecting the elements of a native Abaqus/CAE mesh or the orphan elements when you specify the region of the layup. You must create a layup for each layer of elements.”这里大家都能看懂，我就不详细解释了。简单来说，就是你在厚度方向上定义一层网格和多层网格时，在复合材料建模的过程中需要注意厚度尺寸上网格数量的问题。用一层网格可能完全足够做力学分析，但是在传热过程中，更希望厚度上有更多层的网格，这样才能更好地观察到厚度方向上是否存在温度梯度或者固化不均匀。 第二个问题是，在使用Abaqus的复材模块时，是不能用来进行瞬时温度场分析的，这也是帮助文档里面提过的（The use of composite solids is limited to three-dimensional brick elements that have only displacement degrees of freedom.）。（不知道最新版本的Abaqus是否可以，如果可以了，请大家给我留言，指正我的错误，避免我误导他人。）在这里进行瞬时温度场分析时，我是采用厚度上四层网格，考虑到铺层为[02/902]s，因此我是两个连续相同的铺层作为一层网格，然后在“Property”模块人为的赋予每层网格材料方向。对于这个模型，四层网格应该足够了，但是如果构件厚度增加，厚度方向上可能需要更多的网格。（本来我想做一个4、8、12层网格的比较，但是到8层的时候，我的电脑表示它负担不了这些网格的计算了。。。）。然后在力学分析的时候，可以使用传热分析时的网格，然后采用复材模块进行铺层，将传热计算结果导入应力分析模型，进行固化变形计算。 固化过程温度和固化度计算结果如图1-3所示，可以看到，在第一个保温平台结束前，树脂实际上是几乎没有固化的，这一过程树脂黏度逐渐降低，以树脂流动为主，树脂固化主要发生在第二个升温阶段。最终固化完成后，树脂的固化度为0.9767，可以认为固化完全。以图1-3(c)中的A点，进行整个固化阶段温度、固化度以及Tg曲线的演化规律研究，如图1-4所示，可以看到由于树脂固化反应放热的原因，在第二个保温平台的初始时刻是存在温度过热现象的，但是并不明显，并且该温度峰值很快消退。因此可以认为构件在厚度方向上几乎是没有温度梯度发生的，在实际的固化残余应力计算过程中，对于该模型，可以将热-化学-力问题，转换为化学-力问题，从而减少建模和计算工作量。
图1-3 复合材料结构固化过程温度和固化度计算结果
图1-4温度、固化度以及Tg曲线的演化规律 2 复合材料固化变形研究 谈到固化变形，就很难回避固化残余应力。本质上来讲，固化过程导致的结果是复合材料构件内部产生固化残余应力，而工程上反映出来的，就是构件脱模后固化残余应力释放，产生固化变形，在此我没有把二者进行严格的区分。但需要指出的是，固化残余应力在宏细观不同尺度上，是有所不同的。宏观层面上，固化残余应力引起固化变形，带来的最主要问题是构件成型精度降低，引起装配问题；而细观尺度上的固化残余应力，也就是纤维、树脂以及界面相上的残余应力，还会进一步对构件的力学性能，例如横向拉压性能产生影响。我现在的想法是把二者放在一个统一的框架内进行计算求解，目前宏观到细观已经搞定，但是细观尺度上的残余应力如何影响宏观，还没有完全解决。 言归正传，固化变形可以说是连接复合材料结构设计、制造以及装配过程的一个核心问题之一。需要指出的是，目前在结构设计过程中，可能对于固化变形的关注并不明显，但相信随着复合材料结构设计制造一体化的进一步推动，在结构设计过程中考虑固化残余应力和固化变形，是不可避免的。 固化变形或者固化残余应力的影响因素很多，例如复合材料自身的材料各向异性，树脂的化学收缩，模具作用、固化不均匀等等。关于固化变形的研究，有两个不可回避的问题就是建模方法和材料的本构模型。建模方法中，如何高效地反映构件结构特征以及固化变形的诱因（例如平板构件和L型构件的固化变形主要诱因就有所不同，对称铺层与非对称铺层固化变形的主要诱因也不一样，等等），本构模型中如何反映材料的实际应力应变关系（例如线弹性与粘弹性模型的区别，或者更严格讲是材料参数的时变特性是否考虑），是能否准确计算固化变形的重点，这些内容现有文献已阐述很多，在此不做赘述。同样以上述模型为例，进行固化残余应力和固化变形的计算。 2.1 材料本构模型构件结构模型与铺层形式同样来自上一章的结构。鉴于上一章的计算结果，在此将热－化学－力问题，转换为化学－力问题。在固化残余应力分析中，关键的一点是确定材料的固化力学本构方程。本构方程表示的应力－应变－时间的关系，计算过程中需要确定当前计算步的雅克比矩阵和当前增量步中应力大小。对于固化过程中复合材料力学本构关系，早期采用线弹性本构，但是复合材料的性能与固化度、温度等多个因素有关，线弹性模型显然有很大的局限性。目前常用的固化力学本构模型有CHILE（cure hardening instantaneous linear-elastic）和黏弹性模型。还有一种是对CHILE模型的进一步简化，我称之为阶段跳跃模型（目前没有看到非常合适的称呼方法），即将整个固化阶段分为液态或者黏流态（viscous state）、橡胶态（rubbery state）和玻璃态（glassy state）三个状态，如图2-1所示，每个状态有各自的材料性能参数。对于AS4/8552而言，黏流态与橡胶态的分界点在固化度达到0.31，而橡胶态和玻璃态的分界点则发生在Tg大于固化工艺温度。并假设在黏流态没有残余应力产生，因此只需要计算橡胶态和玻璃态产生的残余应力即可。这种简化方式可以基本反映复合材料固化过程中的力学性能变化，同时计算效率又很高，对于处理厚度较薄的复合材料L型或者U型构件的固化变形问题是非常有效的[3-5]。通过上一章的温度和固化度计算结果，我们可以找出三个阶段的分割点，从而可以进一步将化学-力问题，转化为单纯的力问题，从而进一步简化模型。模型材料力学参数见表2-1。
图2-1 复合材料固化过程典型三阶段[5] 表2-1 模型中所用到的材料参数[3-5]
2.2 模型计算结果采用UMAT和UEXPAN两个子程序进行橡胶态和玻璃态两个阶段的材料参数整合，UMAT子程序负责处理模量和泊松比，UEXPAN则负责热膨胀系数和树脂化学收缩系数的编写。橡胶态和玻璃态时将构件内外表面约束，在回弹步，则整体自由放开，结果如下图2-2所示，可以看到，由于结构长度方向上存在一定的曲率，因此固化之后不仅弯边存在回弹，整体也存在扭曲现象，这种情况下修模的难度还是比较大的。这里需要指出的是，采用不同的回弹约束条件，固化变形情况是有所不同的，这还要根据实际操作过程和工程经验进行判断。
图2-2 固化变形（回弹时完全放开） 3 总结 可以看到，上述的计算模型与结构其实都相对简单，并且所采用的本构模型也很容易编入UMAT。但实际工程应用中的结构要比文中的结构复杂的多，图3-1是最近几年做过的一些实际结构。在此有几个问题，是我在一直思考的，希望可以集思广益。 （1）树脂的流动主要发生在液态（黏流态）阶段，这个时候树脂实际上并没有开始固化，或者固化度非常小，那这个阶段究竟对固化残余应力有多少影响？真的可以忽略吗？ （2）对于结构复杂的大尺寸构件，例如加筋壁板，采用实体单元进行计算固然准确，但计算效率如何保证？而采用壳单元，是否可以将固化变形的诱导因素都考虑在内？ （3）现有的模型多是分别针对不同类型结构的，例如L型或者U型主要考虑圆角处的回弹，平板则是考虑模具作用或者铺层角度的变化，但如果对于厚度较厚的复杂形状构件，则可能发生纤维体积分数分布不均、固化不均等其他因素，如何能在复杂模型中考虑尽可能多的因素，也是目前建模需要考虑的问题。 （4）虽然现有一些软件，例如ESI-RTM，Compro等软件，也可以做复合材料固化仿真，但如何突破国外软件的封锁，在复材领域也任重道远。 参考文献 [1] Wucher B, Lani F, Pardoen T, et al. Tooling geometry optimization for compensation of cure-induced distortions of a curved carbon/epoxy C-spar[J]. Composites Part A, 2014, 56:27-35. [2] Bellini C, Sorrentino L. Analysis of cure induced deformation of CFRP U-shaped laminates[J]. Composite Structures, 2018, 197:1-9. [3] Ersoy N, Garstka T, Potter K, et al. Modelling of the spring-in phenomenon in curved parts made of a thermosetting composite[J]. Composites Part A, 2010, 41(3):410-418. [4] Cinar K, Ersoy N . Effect of fibre wrinkling to the spring-in behaviour of L-shaped composite materials[J]. Composites Part, 2015, 69:105-114. [5] Cinar K, Oztrk UE, Ersoy N, et al. Modelling manufacturing deformations in corner sections made of composite materials. J Compos Mater, 2014, 48(7):799–813 |
||||||||||||
|
|
