编织复合材料预制体铺覆成型的数值模拟

王 波 李 昂 杨振宇

（1 中国航空工业集团有限公司防务工程部，北京 100022）

（2 北京航空航天大学固体力学所，北京 100083）

文 摘 针对平面编织复合材料预制体在曲面铺覆过程中的局部变形，建立了三维有限元模型，利用商业有限元软件Abaqus 模拟了预制体在铺覆成型过程中的纤维束变形规律。研究了0°和45°碳纤维织物在球面铺覆过程的变形过程和局部剪切变形。结果表明，纤维束之间的滑动和纤维束起皱是该预制体在球面铺覆过程中的典型变形模式。在0°织物的球形铺覆变形中，0°和90°纤维束的剪切变形角最小，45°方向纤维束的剪切变形角最大；在45°织物的铺覆变形中，0°和90°纤维束的剪切变形角最大，45°方向纤维束的剪切变形角最小。

0 引言

编织复合材料在航空航天、船舶和汽车等领域具有非常广泛的应用。随着结构的复杂程度不断提高，编织复合材料常常需要制备具有复杂曲率的三维结构，比如球面、锥面等[1]。平面编织复合材料在纤维铺覆过程中，需要根据结构的形状对预制体施加预先的变形，因此铺覆过程中纤维的位置、纤维取向和局部纤维体积分数会重新分布。复合材料内部的纤维分布和取向决定了各类载荷作用下的复合材料内部的传力路径，纱线束之间的相互作用形式也决定了复合材料的损伤和破坏机制[2]。因此，获取预制体铺覆成型后纤维束的变形规律，对于认识最终复合材料结构的整体力学性能至关重要。

关于纤维预制体的在铺覆过程中的变形研究，受到国内外学者的广泛关注。C.MACK 等[3]提出了渔网模型对织物的铺覆过程进行了预测，建立了可以适用于大部分曲面的通用微分方程，并且利用半球铺覆实验验证了其可行性。R.E.ROBERTSON 等[4]在渔网模型基础上提出了更接近实际的纤维假设，可以更加精确的预测连续纤维织物的无褶皱铺覆。此后，研究者对渔网模型的映射算法、假设和初始条件等进行了完善[5]，并将渔网法应用到了更加复杂的曲面上，如双曲率曲面[6-7]。为了准确的仿真织物铺覆过程中的变形过程，有限元方法开始在编织复合材料铺覆的仿真预测中被广泛的开发和采用[8]。M.A.KHAN等[9]建立了一种以3D 桁架单元代替单向纤维、以3D 膜单元代替树脂的预浸料模型。M.DUHOVIC等[10]通过模拟纤维的实际制造过程开始着重于模拟纤维中复杂的接触关系，并且以3D 梁单元来模拟单向纤维。彭雄奇等[11]基于各向异性超弹性本构模型，研究了碳纤维编织复合材料在冲压成型过程中的纤维重新排布和重新取向，并通过实验对数值模拟进行了验证。

本文通过建立平面编织复合材料球形铺覆的有限元模型，模拟了碳纤维平面编织预制体在成型过程中纤维束的变形规律，详细分析了各个部位的纤维位置和取向的变化。本文的研究为评估预制体质量以及后续预测结构的整体力学行为提供了有意义的参考。

1 有限元模型

1.1 平面编织复合材料预制体单胞结构

图1为单胞模型的几何尺寸示意图，其中纱线束的横截面几何尺寸a，b 是纱线的横截面长、短轴的长度。单胞模型中的单胞经纱/纬纱方向长度Wx与单层平纹织物厚度方向上的长度Wy决定了单元胞体的大小。单元胞体模型中几何参数之间的关系可由如下的公式计算而得。

图1 平面编织预制体的单元胞体模型Fig.1 Unit cell model for the plane woven composites preforms

纱线束曲率半径大小R 为

纱线束之间厚度方向的间隙大小f 为

根据纤维间的集合关系，纱线束直线段的长度大小L 为

纱线束横截面的面积大小S 为

纱线束中直段的倾角φ 的大小为

由于经纱和纬纱方向上，单胞模型各包含两条纱线束长度的纤维，因此单胞模型中单方向的完整纱线束长度L0为

纤维体积分数V 为

式中，η 为纤维敛集率。根据以上参数之间的联系[12]，即可在有限元软件中构建出平面编织复合材料的单元胞体模型。

1.2 预制体成型的有限元模型

平纹织物预制体的球形铺覆实验如图2所示，采用一个球形压头对纤维布作为一个整体进行冲压。

图2 铺覆成型的有限元模型Fig.2 Finite element model for the draping

为了保证模型网格的计算精度，纱线束均采用结构化网格，并选择Abaqus 中的线性减缩积分单元C3D8R。由于经纱和纬纱的边界处由于截面变化较大，需要对截面进行网格细分来控制边界处网格的质量。同时，为了减小求解器的计算量，我们将装配网格模型分别导出为经纱纬纱的孤立网格，这样在不进行质量缩放的情况下大大提高模型计算效率。

由于纤维预制体在受到冲压作用而变形的过程中涉及到大量的接触作用，本文中均采用罚函数的方法来处理接触问题。惩罚接触力可以表示为

式中，k 是惩罚刚度，d0是未解决过程的初始补偿距离，dcur是当前补偿距离。如果dcur减小至d0，则将d0重置为dcur。标准的库仑摩擦模型认为，当等效摩擦应力τeq小于与接触压力p 成比例的临界应力τcrit=μp 时不产生滑移。

因此，可以根据临界剪应力提出滑移条件如下

式中，τmax为用户自定义的最大剪应力值。在等效应力大于等于临界应力时，会发生滑移。如果此时摩擦是各向同性的，则滑移的方向和摩擦应力方向一致，则有

根据库仑模型，摩擦系数被定义为一个与等效滑移速率和接触压力有关的公式

因此，本文中摩擦系数只与滑移速率和接触压力有关。根据已有实验，本文模拟中采用的纤维束间的摩擦系数为0.22[13]，摩擦性质的方向性为各向同性。

2 结果与讨论

2.1 织物的铺覆变形

由图3可见，随着半球模具的冲压作用，0°纤维预制体不断向半球表面铺覆。同时由于压边圈的作用，大部分纤维呈拉伸状态。对于纤维取向为0°和90°的纤维来说，其在轴向的弹性模量远大于其横向上的弹性模量，因此纤维预制体的0°和90°方向上的纤维拉伸变形较为明显，边界中点附近产生了很大的向内缩进，边角处几乎无缩进。45°方向上相邻经纬纱之间的角度变化较为明显，可以看出这部分纤维主要靠互相之间的角度变化来适应整体结构的大变形过程。图4为45°纤维预制体在半球表面铺覆过程的变形情况，边角附近产生了很大的向内缩进，而边界中点缩进并不明显。0°和90°方向上相邻经/纬纱之间的角度变化较为显著，表明存在较大的局部剪切变形。

图3 0°织物的变形过程Fig.3 Deformation process of the 0° fabric

图4 45°织物的变形过程Fig.4 Deformation process of the 45° fabric

2.2 纤维束的局部变形

图5给出了0°平纹编织碳纤维预制体在半球曲面上的纤维局部剪切变形角分布情况，其中图5（a）横坐标为测量点处的半径Ri与半球半径R 之比。由图5（a）可以看出，纤维预制体在0°和90°两个方向上的剪切变形角十分接近于0，变化趋势趋于平缓。很显然这是由于0°织物的纤维取向为0°和90°，碳纤维在这个方向的铺覆成型过程中主要承受轴向拉伸作用，变形主要为轴向变形。同时碳纤维的轴向拉伸模量较大，因此相邻纱线间并没有产生明显的角度变化；相反，在±45°上的纤维随着数据点的位置Ri/R 的增大，剪切变形角不断增大，且变化趋势越来越快，最终在临近半球赤道处达到了剪切变形角的最大值40.66°，并产生了局部起皱现象。这主要是由于在±45°方向，纤维轴向和横向均受力，而碳纤维的横向弹性模量远低于其拉伸模量，因此在横向产生的变形远大于轴向，导致纱线间产生的角度变化较大，剪切变形明显。同时相邻纱线间的缝隙也相应增大，局部产生明显的间隙缺陷。在半球面上，纬度60°处，纤维预制体的纤维取向的变化较小，剪切变形不明显，如图5（b）所示。而在纬度为30°处，纤维预制体的纤维取向的变化较大，基本规律为0°、90°、180°经线处剪切变形角趋近于0；在45°和135°经线处剪切变形角分别达到最大值27.6°和27.8°，同时0°～90°区域和90°～180°经线区域的变化规律基本呈对称趋势。

同样，图6给出了45°平纹编织碳纤维预制体在半球曲面上的局部剪切变形角分布情况。由图6（a）可以看出，碳纤维预制体在45°和-45°两个方向上的剪切变形角十分接近于0，变化并不显著。相反，在0°和90°上纤维的剪切变形角随着半径比Ri/R 不断增大，最终在临近半球赤道处达到了剪切变形角的最大值41.29°。图6（b）可见，在纬度60°处，纤维预制体的纤维取向的变化较小，剪切变形不明显；而在纬度30°处，纤维预制体的纤维取向的变化较大。表现为45°、135°处纤维的剪切变形角趋近于0°，而在0°、180°和90°经线处剪切变形角分别达到最大值23.5°、22.4°和23.5°，同时0°～90°区域和90°～180°经线区域的变化规律呈对称趋势。

图5 0°织物半球曲面上不同方向的纤维局部剪切变形角分布Fig.5 Local shear deformation angle distribution of fibers in different directions on the hemispherical surface of 0° fabric

图6 45°织物半球曲面上不同方向的纤维局部剪切变形角分布Fig.6 Local shear deformation angle distribution of fibers in different directions on a hemispherical surface of 45° fabric

通过数值模拟的结果与已有实验之间的对比（图7），可以发现本文的有限元模拟很好的预测了球形曲面铺覆中平面编织纤维网络的变形。

图7 数值模拟结果与实验[14]的对比Fig.7 Comparison between the numerical simulation and experiments

从成型过程中纤维整体的变形规律来看，预制体中纱线都会受到的拉伸、压缩和剪切，其中拉伸导致了边界的缩进；压缩（主要由压边圈限制作用）导致了部分纤维的屈曲；剪切导致了纤维层间的大角度变化。这三者的共同作用之外还会产生层间滑动和起皱等缺陷，如图8所示。

图8 纤维束之间的滑动和弯曲变形Fig.8 Sliding between fiber bundles and bending deformation

起皱是由于纤维预制体在铺覆过程中需要由平面转变为三维曲面，因此在平面向曲面过渡的区域内预制体局部的几何形状会发生突变，单元会产生较大的变形，即纱线间产生更大的接触压力。此时在过渡区域会产生局部的较大的压缩变形，纤维会发生类似梁的“翘曲”，使得半球赤道处的过渡区域的经纱和纬纱发生起皱现象，与BOISSE 等[15]提出的由于纱线压缩变形导致的起皱现象吻合。

3 结论

本文采用有限元法分别对0°和45°碳纤维平面编织布在球形压头冲压成型的过程进行了仿真，获得了与实验一致的变形模式，详细分析了纤维束在成型过程中的剪切变形与分布情况。这些数据对不同取向编织物铺覆成型的结构设计和编织布的剪裁具有重要的意义，对于该碳纤维布的成型性能的准确预测，也具有很强的实际应用价值。通过对比分析，得到了以下结论:

（1）滑动和起皱是编织复合材料预制体在球面铺覆过程中的典型变形模式，纤维的局部剪切变形角可以用作碳纤维编织物铺覆成型性能的一种评判标准；

（2）铺覆过程中纤维的局部变形与织物中纤维的取向角密切相关，通过改变织物的铺层角度，可以调控成型后预制体在不同方向上纤维的变形程度；

（3）本文的有限元方法很好的预测了球形铺覆的成型过程，并且对平纹编织复合材料预制体的变形做出合理的预测，该方法还可望推广应用于其他复杂曲率的不可展曲面铺覆结构的预测。