连续纤维热塑性复合材料热成型仿真研究

毛晋轩 刘 东 史鹏程 颜 春 祝颖丹

（1 江西理工大学，赣州 341000）

（2 浙江省机器人与智能制造装备技术重点实验室，中国科学院宁波材料技术与工程研究所，宁波 315201）

文 摘 连续纤维热塑性复合材料（CFRTP）拥有优异的力学性能、抗疲劳性和设计灵活性，在航空航天、汽车、能源和海洋工程等军民领域具有广泛应用。为了满足大规模高效低成本制造的需求，可在CFRTP 产品开发阶段使用仿真模拟手段代替昂贵的“试错”实验。本文详细概述了CFRTP 热成型模拟仿真的各类仿真方法研究进展，总结了各类方法的优势及其局限性。最后展望了CFRTP热成型模拟仿真的国内未来研究重点及发展趋势。

0 引言

根据树脂基体类型的不同，连续纤维增强树脂基复合材料可分为热固性（CFRTS）和热塑性（CFRTP）［1］两种。相比前者，CFRTP 具有成型周期短、预浸料存储条件简单且保质期长、易于回收等优点［2-3］。当前CFRTP 成型工艺类型主要有热压、树脂传递模塑成型（T-RTM）、拉挤、缠绕和自动铺放等，其中热压成型是以热塑复合材料预浸料或层合板为原材料，经加热软化后在模具中快速热压成型，而TRTM 需要将干纤维预先铺覆到模具中，然后进行注胶、加热，树脂固结成型。上述过程中软化材料的热变形-固结行为以及纤维铺覆预成型的质量会对最终构件的力学性能、尺寸精度和缺陷产生决定性影响。与传统短纤、长纤增强热塑性复合材料的成型工艺相比，由于连续纤维变形延展性不足、树脂流变行为复杂等因素，CFRTP 的铺覆预成型和热压热成型阶段会涉及复杂的变形机理，并导致一些典型缺陷，如翘曲、褶皱，甚至纤维断裂等。而国内外大多仍然采用昂贵的“试错”实验进行产品开发，缺乏有效的理论模型开展工艺及性能预测研究，严重影响产品质量、研发周期和推广应用［4］，因此开展有关CFRTP 成型过程仿真研究，准确预测纤维变形、缺陷产生和残余变形等行为，对于优化CFRTP 成型工艺，提高制品质量具有重要意义。

当前用于CFRTP 成型仿真的方法大体上可以分为两类：运动学法与力学法。本文将评述各种模拟仿真方法的发展历程及近年来最新研究成果。

1 运动学法

运动学法也被称为映射法、渔网算法或销钉连接法（pin-jointed net），基本原理是利用几何投影的方法，将二维平面的织物投影至三维曲面表面［5］。C.Mack［6］最早建立了映射法的雏形，他忽略了纱线伸长和滑移带来的影响，推导出了可预测可展曲面铺覆效果的微分方程，并通过对球面进行铺覆实验，验证了该方程。后续学者对于运动学法铺覆变形的完善与拓展，普遍基于该模型。Van Der WEEËN［7］认为织物节点间的纱线投影到三维曲面上是测地线，研究了曲面上测地线的解法。BOROUCHAKI 等［8］提出了一种新的运动学模拟算法，该算法考虑了绘制在表面上的织物网格单元的真实几何形状。这样的织物网格单元由曲线四边形定义，该四边形的边缘是绘制在待铺覆表面上具有相同长度的测地线。杨波等［9］提出了一种基于几何信息的映射算法，该算法利用了相邻节点处的曲率及切向量，确定该节点的位置，避免了高强度的迭代计算，提高了仿真效率。KAUFMANN 等［10］利用运动学法与其建立的铺覆数据库结合，从纤维角度偏差、余料和材料剪切角度等方面评估仿真种子点和参考角的组合，试图在部件的结构性能和制造成本之间找到平衡。

也有研究者另辟蹊径，将运动学法与力学法相结合。SHARMA 等［11］提出了一种简化的有限元模型来模拟铺覆过程，其单胞模型如图1所示。

该单胞由销钉连接的刚性桁架网络组成，通过对角单元引入剪切刚度，材料属性可从拉伸试验获取。该单元模型能够代表平面内力对剪切变形和纤维滑移的影响，同时作者预测该模型也适用于预浸料，因为该模型可以考虑结合预浸料的剪切响应进行模拟。

运动学法极大地提高了计算效率，然而这种模型做出的假设条件过于理想，未考虑载荷、边界条件以及材料与模具之间的摩擦，缺乏预测褶皱的能力。因此，运动学法主要应用于预浸料的手工铺层与纯纤维织物铺覆性研究，或者在初步的设计阶段，对最终的成型效果做简单预估。

2 力学法

与运动学法相比，力学法有限元模拟仿真考虑了复合材料的力学性能，并且可以与成型的物理过程相耦合。用于模拟CFRTP 热成型的力学模拟方法可分为三种：连续法、离散法以及半离散法。每种力学法耗费的计算成本差别很大，连续法成本最低，而离散法最为昂贵，半离散法介于两者之间［12］。

连续法适用于宏观现象的模拟和全结构问题的计算，离散法适用于模拟介观尺度与微观尺度的现象，而半离散法因其特性对所有尺度均可适用。图2显示了织物复合材料的三种不同的等级尺度。

2.1 连续法

连续法是在复合材料连续介质力学的大框架下，假设复合材料是一个均质化的连续体，这是目前为止研究最为广泛的方法。因为该模型仅在材料行为方面与经典连续结构有所不同，所以对各类结构具有良好的适用性［13］。连续法模拟仿真分析本质上是金属成形仿真过程中深拉伸的一种拓展，可采用经典有限元软件实现［14］。

早期的连续法建模，采用了理想化的纤维增强体模型，假设单向或织物预浸料中的纤维为刚体，因此早期研究局限于成型主要的变形模式—面内剪切。如HSIAO 等［15］基于均质化假设，假设在黏性树脂流体中的是刚性纤维，对热成型过程进行了数值模拟，所提出的各向异性材料参数与温度有关。该模型既能计算纤维取向，又能计算冷却阶段产生的残余应力。然而，理想化的纤维增强模型不仅给后续的二次开发带来困难，也难以提高仿真精度。此后随着研究者进一步改良与有限元仿真技术发展，逐渐解决了有限元编程困难的问题。YU 等［16］提出了一种新的描述大变形下非正交材料行为的本构模型，旨在更好地表征独立于网格剪切的面内正应变。该模型能准确捕捉材料在不同加载路径下的响应，真实地反映复合材料中纤维重分布和重定向。Guzman-Maldonado 等［17］基于麦克斯韦流变模型的推广建立了一个非线性粘-超弹性模型，该模型是通过迭代的热分析和成型分析来实现的，保证了冲压变形和温度场之间的耦合。计算得到的剪切角与实验结果吻合较好，其模型的仿真结果如图3所示。结果表明，在成型过程中温度场发生了明显的变化，单胞变形和模具接触改变了局部热性能和温度，而由于温度变化和高应变率，面内剪切刚度增加，导致了频繁起皱。在实际的热成型过程中，温度对预浸料的机械性能有着极大的影响，由此可见热力耦合模型对于提高CFRTP的热成型仿真精度是不可或缺的。

随着近年来研究者们证明了弯曲刚度对缺陷褶皱产生有着不可忽视的影响，尤其是在确定褶皱的形状和大小方面之后［18-20］，不少研究者都将弯曲刚度纳入了本构模型。考虑了弯曲刚度的本构模型，仿真将更贴合实际情况，如图4所示［21］。

DOMINIK 等［22］基于沃伊特-开尔文模型和广义麦克斯韦模型提出了一种模拟单向或织物增强体的应变速率弯曲行为模型，该模型考虑了平面剪切、拉伸、弯曲刚度等力学行为。根据热塑性复合材料单向带在加工条件下的热表征结果，成功地对所提出的黏弹性模型进行了参数化。研究表明，广义麦克斯韦方法结合非线性粘弹性行为最适合预测褶皱的产生及发展过程。

在研究某些特定问题时，除了考虑面内剪切、拉伸、弯曲刚度等常用的力学行为之外，还会考虑其他力学因素，如SOULAT 等［23］研究成型过程中的孔隙率问题时，为了避免在实际成型仿真过程中对层合板每层进行单独计算，开发了一种具有横向应力的壳单元，该壳单元增加了一个自由度，允许单元模拟厚度的变化，通过解耦弯曲和收缩应变来避免体积锁定现象。结果表明，所获得的横向正应力与去除孔隙的结果非常吻合。

总体而言，考虑的力学行为越多，仿真模拟得到的结果会越接近实际情况，然而现实中我们却不得不考虑计算成本，所以根据研究问题的实际需要，对力学行为进行取舍，可以忽略对最终结果影响小的力学行为以提高计算效率。连续法建模时由于假设连续介质代替了连续纤维，在这种情况下，建模的困难之处在于宏观模型必须考虑到纤维材料的异质性，特别是在纤维大变形情况下，否则仿真结果与实验结果将出现较大的偏差。因此，连续法要用相对复杂的计算策略来设置单元中的纤维方向。

2.2 离散法

离散法将纤维增强相在介观尺度下看作一系列的离散单元，该方法下的有限元模型由桁架、梁、壳或膜的网格单元组成。相比于连续法，这种模型的主要优点是直接考虑了材料的异质性，不需要复杂的计算策略来设置单元中的纤维方向，即可很清晰地描述纤维增强体内部结构［24］。

一般情况下，使用离散模型来模拟复合材料成型过程的研究目的是预测无褶皱成型结构，最基本的模型仅考虑了变形过程中纤维方向的计算。例如SIDHU 等［25］提出了一种棋盘模型，该模型使用三维桁架单元和三维壳单元来模拟织物复合材料预制体的成型，其模型单胞如图5所示。

该模型考虑了纤维间摩擦、丝束干扰和纤维滑移。结果表明，大变形过程中随着剪切角的变化，纱线间出现明显的滑移。JAUFFRÈS 等［26］提出了一种基于亚弹性描述的显式有限元公式离散模型，纱线的拉伸行为由桁架、梁等一维单元模拟，织物的剪切行为由壳或薄膜单元模拟，该模型有限元分析结果与实验数据基本一致。SHERWOOD 等［27］在该模型的框架下，考虑了复合材料热冲压过程中模具与织物之间摩擦，使用恒定和可变摩擦因数做对比实验，进行半球冲压模拟。结果表明提高冲压速度会导致模具与织物界面的摩擦力增加，从而增加了纤维的拉伸应力，且冲压力与冲压速度无关。HARRISON等［28］提出了一个多尺度能量模型，使用桁架和壳单元来模拟黏性织物复合材料的热成型，通过将多尺度能量模型与宏观尺度模型相结合，使得该模型能够预测复合材料成型过程中的剪切应力-剪切角-剪切速率行为，适用于在不同剪切速率下的成型模拟。

此外，离散法的特性允许模型更容易预测一些细观现象，如BOISSE 等［29-30］基于纤维增强体离散模型的宏观成型，预测了纤维之间的大滑移现象。在该模型中每个织物单元被模拟为一组非常简化的壳单元，在保证自由度很低的情况下，还能描述纤维之间的摩擦接触，如图6所示。

离散法所耗费的计算资源巨大，计算效率低，RAMGULAM 等［31］认为，虽然详细的有限元分析对最终方案必不可少，但近似的力分析有利于在工艺优化过程中快速筛选大量的铺层方案。因此，他提出了一种快速计算织物褶皱的近似力学分析方法。利用一种基于微分几何的铺覆算法，如公式（1）、（2）所示，获得了织物在曲面上变形后的纱线取向。

式中，K是高斯曲率，α、β为织物内角。通过调整外部载荷和给定织物的剪切性能，可以快速确定织物铺覆时出现问题的区域。此外，压力分布可用于计算层间摩擦力。这种计算效率高的力/褶皱方法还可用于扩展运动学褶皱模型的能力，以补充全尺度有限元分析。

离散法主要的局限性在于因为将每根纤维看作独立单元，所以必须考虑大量纤维复杂的机械性能以及它们之间的相互作用，且平面内剪切和拉伸之间是相互依赖的，很难考虑刚度之间的耦合。离散法通常只用在介观尺度上分析小数量的单元，并且这个级别一般不适合分析整个成型模拟过程［32］。因此，离散法建模时需要在单胞模型的准确性和自由度总数之间进行折中。单胞的建模必须足够准确才能获得正确的整体宏观力学行为，且每个单胞的自由度数要尽量小，以便计算成千上万个单胞的成型过程［33］。

2.3 半离散法

半离散法是在介观尺度上构建的特定有限元。一方面，与离散法一样，材料的力学行为被拆分，仅考虑主要刚度（如拉伸、面内剪切或面外弯曲）；另一方面，用与连续法相同的方式，将所有刚度都视为在半离散单元内部，可以很容易地将刚度参数传递到宏观模型的半离散单元中［34］。

HAMILA 等［35］提出了一种考虑拉伸和面内剪切行为的三结点单元。经纱和纬纱方向相对于单元侧面是任意的，材料参数是通过标准拉伸和偏轴拉伸试验确定。而后续HAMILA 等［36］改良了这种半离散三角形单元，除了考虑拉伸刚度和平面内剪切刚度之外，还考虑了弯曲刚度。弯曲行为由悬臂弯曲试验给出，半离散单元的弯曲曲率由相邻单元的位移得到，就铺覆过程中褶皱的出现和发展分析了所考虑的三种刚度的影响。CHEN 等［37］集中研究了热成型过程中CFRTP 的层间剪切行为，提出了一种层合板粘弹性模型，其模型示意如图7所示。该模型应用了具有八个节点的大位移三维内聚单元，研究织物热塑性复合材料的层间剪切机理。模型考虑了不同加工温度和不同应变速率下的拉伸、压缩和面内剪切行为。通过该模型，可以模拟温度对应变率的影响规律。WANG 等［38］在该模型的基础上考虑了熔点附近不同温度下铺层的张力、面内剪切和弯曲刚度，模拟了多层CFRTP 热成型过程中的接触摩擦，并使用拉格朗日乘子法进行计算，仿真得到的成型剪切角和起皱角与热成型实验结果吻合较好。

半离散法旨在避免使用应力张量，仅通过张力以及面内剪切和弯曲刚度直接定义单元上的载荷。这些量简单地定义在一个单胞上，而且这些材料参数可以通过复合材料的标准测试直接获得。半离散法也可以与经典有限元技术一起使用，并且通过对特定有限元单元进行定义考虑纤维的异质性，进而提高计算精度。

3 结语

目前，在CFRTP 成型模拟仿真领域仍存在一些重要的科学问题有待解决，例如没有一种较为系统的流固热力耦合分析模型，缺少针对成型缺陷成因机理的探索，只是简单的定性分析。同时，在介观与宏观尺度建模的相关文献中可以看出，开发一个同步的多尺度模型来模拟仿真成型过程仍然是一个巨大挑战。

总体而言，我国CFRTP 成型的模拟仿真仍停留在实验室阶段，尚不能够满足实际生产的需求。为了解决高性能CFRTP 大批量生产所需求的高效低成本、稳定可靠的制造工艺，对CFRTP 成型的模拟仿真依然是未来树脂基复合材料领域的重要研究方向之一，值得进一步深入研究。