三维编织复合材料是将纤维束经一定的编织工艺技术形成预编织件，然后与基体材料固化成型得到的。三维编织技术起源于20世纪60年代，目前常用的制备工艺有：四步法编织、二步法编织和实体编织技术。其中四步法编织于1982年由Florentine^[1]发明出来，该项编织技术得到最广泛的应用。基于编织工艺规律以及编织物真实细观结构，一些学者针对编织物的细观结构和几何建模进行了大量分析和研究。Ko^[2]首次用一个四角交叉纤维束的正方体对矩形编织预制体的纤维结构进行描述。由于该模型描述过于简单，Wu^[3]和Wang等^[4-5]对三维编织预制件进行了拓扑分析，把结构定义为可区分的内部、表面和棱角处三种单胞组成。对于三单胞模型的选取，目前主流有两种划分方案^[6]：水平划分法和45°划分法，其中45°单胞划分方案中只适用于主纤维束行列数均为偶数情况，若其中有一个为奇数，该划分方法就不能使宏观结构完全地包含三种单胞，而是在角处出现多余的正方形部分。研究者^[7-9]利用有限元软件对三维编织过程进行模拟，并建立相关实体模型，或者建立起细观编织结构和编织参数之间的关系，并给出了相应的数学公式^{[6, 10]}。

目前学者获得三维编织物弹性常数的方法主要有三种：理论计算、数值分析以及实验测试。在理论分析方面，主要有刚度平均化方法^[11-12]和均匀化理论^[13]，这些方法主要基于纤维束和基体力学性能通过计算和转换对整个编织体的弹性常数进行预测。在数值分析方面，一些学者利用有限元软件ANSYS, ABAQUS等对单胞模型或预制体细观结构模型进行数值分析，预测整个编织体的弹性常数。在实验方面，文献[14-16]中主要开展三维编织物的拉伸、压缩、剪切等实验，获得三维编织物的弹性常数。

大部分学者仅从理论计算、数值仿真以及实验三种方法中三选二或三选一研究编织复合材料刚度特性。本工作不仅从理论计算、数值仿真以及实验三个方面研究三维编织复合材料弹性常数随编织角变化规律，还对现有文献中纤维束不同截面形状假设对纤维体积含量的影响进行对比与总结，同时指出了三维编织复合材料应力应变曲线随编织角的变化趋势。本工作理论计算采用刚度平均化方法，不仅将其应用于内胞模型，还应用于面胞模型；数值分析方法则利用ABAQUS软件不仅在内胞分析时施加了周期性边界条件，还在对面胞进行刚度分析时，利用编写的Python程序，在与相邻面胞接触的两个面亦施加了周期性边界条件。

1 刚度平均化方法

三维四向编织复合材料由纤维束和基体组成。在预测其弹性常数时，一般认为纤维束是线弹性、横观各向同性材料，树脂基体是线弹性、均匀各向同性材料。刚度平均化方法具体思路是将各个方向局部坐标系下纤维束的刚度矩阵通过应力转换矩阵转换为整体坐标系下的刚度矩阵，然后和基体的刚度矩阵基于单胞中各组分材料的体积含量进行刚度平均求和，得到整体坐标系下单胞的刚度矩阵，求逆得到整体坐标系下单胞的柔度矩阵，进而得到整体坐标系下单胞的弹性常数。图 1为刚度平均化方法预测内胞和面胞弹性性能的编程思路流程图。

如图 2所示，设定总体坐标系O-XYZ，纤维束局部坐标系O-123，总体坐标系中Z轴为编织方向轴，在局部坐标系中1轴为单向纤维束主轴。局部坐标系中1轴、2轴、3轴相对于总体坐标系X轴、Y轴、Z轴的方向余弦分别为l₁，m₁，n₁，l₂，m₂，n₂，l₃，m₃，n₃，相应关系和表达式见表 1。

表 1中：l₁=sin(ϕ)cos(φ), m₁=sin(φ), n₁=cos(ϕ)cos(φ), l₂=sin(ϕ)sin(φ), m₂=-cos(φ), n₂=cos(ϕ)sin(ϕ), l₃=cos(φ), m₃=0, n₃=-sin(φ), φ=45°。

可知应力转换矩阵为：

局部坐标系下刚度矩阵和整体坐标系下刚度矩阵之间关系为：

(1)

得到纤维束在整体坐标下的刚度矩阵后，根据内部单胞中纤维束和基体的体积含量以及刚度矩阵，可以得到内部单胞的刚度矩阵C_n，见式(2)：

(2)

式中：C_m为基体的刚度矩阵；V_nfi表示第i种方向纤维在内部单胞中体积含量; V_nf为所有纤维束在内部单胞中体积含量，且。

对于内胞和面胞中纤维束的分析，唯一区别在于，内胞中16根纤维束在空间走向都是直线，有4种空间角度，而面胞中采用文献[17]的处理结果，面部区域中纤维束在空间走向为折线，有8种空间角度，具体见表 2。

表 2中第一个坐标表示纤维束与编织方向(Z轴正方向)的夹角，即图 2中的ϕ，第二个坐标表示纤维束在XY的投影与X轴正方向的夹角，即图 2中的φ，γ为内部编织角，θ为面部编织角，它们与编织角α之间的关系为，可知编织角为20°，30°，45°的编织物另外两个角度见表 3。

得到内胞和面胞在整体坐标系下的刚度矩阵，根据实验件中内胞和面胞所含的体积含量，加权求和(见式(3))，可以预测实验件的弹性常数。

(3)

式中:V_n为编织物中内部单胞占实验件的体积分数；V_m为编织物中面部单胞占实验件的体积分数；C_m为面胞的刚度矩阵。

实验件材料选取天津工业大学研制的三维四向碳纤维/环氧编织复合材料，预成型件采用三维四向1×1四步法编织制造，通过RTM工艺固化成型，编织物中纤维体积分数61%~62%。增强纤维为12K-T700碳纤维，基体为TED-86环氧树脂。细观组分材料的性能参数见表 4。根据文献[18]以及厂家提供的m和n值可以计算出厚度为4mm，编织角为20°，30°，45°实验件的内胞体积分数和面胞体积分数，见表 5。

经过编程计算，得到编织角为20°，30°，45°三维四向编织复合材料内胞和面胞以及加权求和后预测的弹性常数，具体结果见表 6。

对于大厚度的编织复合材料，编织物内的内胞含量至少超过60%，因此可以用内胞的弹性常数反映整个编织体的弹性常数。图 3，4，5分别是根据内胞模型得到三维四向编织材料拉伸弹性模量、剪切模量以及泊松比随编织角的变化规律，纤维体积分数取62%。可以看出：(1)随着编织角α的增加，E_Z逐渐减小，E_X，E_Y，G_XY逐渐增加，主要原因是α增大，纤维束在Z轴方向刚度分量逐渐减小，在X，Y方向刚度分量逐渐增加；(2)随着α的增加，G_YZ，G_ZX先增大后减小，在32°附近取得最大值，这是因为随着编织角的增加，X，Y方向纤维束的增强作用一直增加，而Z方向的增强作用逐渐减少。当试件编织角达到一定角度时，沿着横向的变形达到最大，而随着编织角继续增加，试件的横向变形逐渐降低。(3)μ_YZ，μ_XZ先增大后减小，在23°附近取得最大值，这是因为随着编织角的增加，受纤维束空间位置的影响，试件三个方向的变形逐渐增大，但X，Y方向变形的增加趋势更大。在编织角达到一定角度时，X，Y方向的变形和Z方向变形的比值达到最大，此时纵向泊松比达到最大值。随着编织角度继续增加，试件X，Y方向变形受纤维束的约束而逐渐减小，泊松比逐渐降低，这使得纵向泊松比先增加后减小。同理，μ_XY先减小后增大，在30°附近取得最小值。

2 有限元分析

由于编织复合材料细观结构的复杂性，现有有限元分析工作通常忽略材料外部区域面胞、角胞的影响，假设材料横截面尺寸足够大，研究占材料主体的内胞的力学行为来表征整个结构的力学性能。

表 7为纤维束在不同截面形状假设下计算所得的纤维体积分数，从表中可以看出，圆形截面假设对应的纤维体积分数最小，六边形截面下的编织物纤维最大体积分数最大。根据天津工业大学提供的三维四向碳/环氧编织复合材料拉伸实验件，编织物中纤维体积分数61%~62%，可知八边形截面假设下，编织物纤维体积含量最接近实验件的纤维体积含量，再结合文献[19]介绍和分析，本工作将采用受挤压纤维束的截面为八边形的假设。此外八边形截面的单胞模型在有限元网格划分中要比圆形截面和椭圆截面要容易一些。

本工作基于内胞和面胞模型，利用HYPERMESH软件对其进行周期性网格划分(其中面胞模型仅在与相邻面胞接触两个面进行周期性网格划分)，然后在ABAQUS软件中施加周期性边界条件，分别对内胞和面胞进行弹性常数预测。

以20°编织角的内胞为例，在HYPERMESH软件中对内胞进行网格划分时，对于周期性网格，必须保证主从面上网格节点数目和位置分布完全一致，采用面网格复制的方法来实现。本工作对纤维束和基体采用四面体单元进行离散，内胞模型由28450个单元，5870个节点组成，如图 6所示。

在ABAQUS软件中利用基于Python软件编写的程序对单胞施加周期性边界条件^[20-22]，然后对纤维束和基体分别添加材料属性、截面属性，考虑到内胞中的纤维束有4种空间方向主轴，对相应的纤维束添加相应的材料轴局部坐标系。然后分别施加6种基于简单应变的方程边界条件(ε_X⁰, ε_Y⁰, ε_Z⁰, γ_XY⁰, γ_XZ⁰, γ_YZ⁰，设定其应变值为0.01)，基于单胞的应力和变形形态，确定其弹性常数。

由于三维四向编织复合材料纤维束空间结构的轮换特征，单胞模型在X拉伸和Y拉伸、XZ剪切和YZ剪切载荷下变形和应力分布具有相似性，以下仅给出X拉伸、Z拉伸、XY剪切及XZ剪切4种载荷下的分析结果。

从图 7单胞变形状态可以看出，相对两个面变形一致或保持常值，变形连续，单胞受载后保持了其在空间上的周期性。从图 8单胞应力状态可以看出，在各种载荷下，单胞模型中纤维束应力均明显大于基体，承受着主要载荷，这主要是由纤维束的刚度远大于基体，载荷按刚度分配决定的。但是，在不同应变载荷工况下，纤维束的承载机理并不相同，所以单胞模型的力学响应特征也不同。

表 8为数值计算得到20°，30°，45°编织角内胞和面胞及其加权求和的弹性常数。

3 实验验证

根据GB/T 1447-2005纤维增强塑料拉伸性能实验方法, 标准设计厚度为4mm、编织角为20°, 30°, 45°的拉伸实验件，拉伸实验件的设计方案如图 9所示，对于纵向拉伸实验件采用Ⅱ型试样，横向拉伸实验件采用Ⅰ型试样(厂家提供的编织物平板尺寸限制，故选用Ⅰ型)，鉴于横向拉伸实验件Ⅰ型尺寸较小，故在拉伸实验中应变率采用1×10^-5 s^-1，而纵向拉伸实验件Ⅱ型尺寸较大，故纵向拉伸应变率采用1×10^-4 s^-1。由于3D编织复合材料的各向异性，测试其面内性质时，分别进行纵向拉伸实验和横向拉伸实验，实验所用设备为MTS793拉压和疲劳试验机。

对实验后的数据进行处理，发现编织角20°编织复合材料试样的纵向拉伸应力应变曲线近似弹性，在最终断裂前显示出弱的屈服；编织角30°编织复合材料试样的纵向拉伸应力应变曲线，相对于20°，30°编织复合材料试样的纵向拉伸应力应变显示出明显非线性；编织角45°编织复合材料试样的纵向拉伸呈现双线性。

对于纵向拉伸，随着编织角增大，可以发现应力应变曲线由线弹性、弱屈服转换为非线性，进而转化为双线性；而对于横向拉伸，应力应变曲线基本呈现线弹性，45°时应力应变曲线略呈非线性。

随着编织角增大，实验件应力应变曲线发生变化的原因，分析如下：(1)实验件是在编织体中切割加工而成的，切割边缘导致实验件中部分纤维变成短纤维，编织角越大，实验件中纤维越短。(2)对于纵向拉伸过程，刚开始时，由纤维和基体共同承受作用力，由于纤维和基体在受力下变形的差异，随后部分纤维与基体界面发生脱胶，编织角较小时，界面脱胶并不剧烈，纤维和基体呈脆性断裂，纤维在纵向拉伸方向分量较大，纵向刚度较大；编织角较大时，短纤维与拉伸方向夹角较大，纤维在纵向拉伸方向分量较小，纵向刚度较小，由于大编织角的实验件的短纤维较多，脱胶的纤维较多，因而后续承受力主要是剩余基体和少量的纤维，因而应力应变关系逐渐呈现双线性；总体来说，随着编织角增大，编织复合材料的破坏从脆性破坏转向延性破坏。(3)对于横向拉伸过程，拉伸方向与编织方向垂直，且实验件宽度较纵向拉伸实验件宽度小，从而在拉伸过程中，纤维脱胶明显，主要是基体和少量的短纤维承受拉力，因此应力应变曲线呈现线弹性断裂，当编织角增大，纤维束的纵向刚度在横向分量增加，从而实验件横向刚度增加。

表 9和表 10分别是实验得到弹性模量和泊松比。表 9和表 10中实验件的代号意义：ST表示静态拉伸，ST后面数字为编织角，L表示纵向，T表示横向。

4 结果与分析

基于实验数据可以得到不同编织角下刚度参数，并与理论计算和数值分析结果进行比较，见表 9。由表 9可知，理论计算和数值分析在实验件编织方向(Z轴方向)刚度的预测精度较高，对实验件横向刚度预测精度比较差，预测值比实验值大，究其原因主要是无论是理论计算还是数值仿真都忽略了纤维与基体在实验件受拉伸过程中发生的脱胶现象，在横向拉伸中由于实验件宽度较小，在拉伸过程中，主要靠基体和短纤维承载，纤维与基体之间界面开裂非常明显。至于泊松比，实验测得泊松比受应变片粘贴位置的影响，可能贴在纤维上、基体上或纤维与基体交汇处，在拉伸过程中，贴片处的纤维出现抽拔或断裂以及基体开裂等现象影响，导致实验测得泊松比仅是实验件局部的泊松比而非整个实验件的泊松比，因此误差较大，实验测得泊松比普遍偏大。从表 11中也可以看出弹性常数随编织角变化趋势，随着编织角的增大，编织方向的刚度减小，横向刚度增大，与图 5中相应参数随编织角的变化趋势吻合。

5 结论

(1) 理论计算发现，随着编织角增大，三维四向编织复合材料纵向刚度减小，横向刚度增大，横向剪切模量增加，这是因为编织角增大，纤维束在纵向分量减小，在横向分量增大。三维四向编织复合材料纵向剪切模量先增大后减小，在32°附近取得最大值，横向泊松比先减小后增大，在30°附近取得最小值，纵向泊松比先增大后减小，在23°附近取得最大值。

(2) 在数值分析中，基于纤维束截面形状的分析，圆形截面假设对应的编织物纤维体积分数最小，六边形截面假设下的编织物纤维最大体积分数最大，其中八边形截面假设最适合用来预测本工作选取实验件的弹性常数。

(3) 基于内胞和面胞模型，刚度平均化方法和数值分析方法在预测三维编织复合材料纵向刚度方面精确度较高，但在横向刚度方面预测精度较差，究其原因是忽略了纤维与基体脱胶现象。

(4) 20°编织角纵向拉伸应力应变曲线基本呈线弹性，30°纵向拉伸应力应变曲线呈现非线性，45°纵向拉伸应力应变曲线呈现双线性。随着编织角增大，编织复合材料的破坏从脆性破坏转向延性破坏。20°, 30°, 45°横向拉伸应力应变曲线基本呈线性。