整体穿刺预制体是一种三维立体织物，由于在厚度方向引入了均匀排列的增强纤维，与二维叠层织物相比，具有良好的层间性能和制备工艺简单的优点^{[ [1]  李彦璋, 贺辛亥, 程攀, 等. 碳纤维立体穿刺预制体成型技术发展现状与思考[J]. 纺织器材, 2024, 51(1): 65–68.LI Yanzhang, HE Xinhai, CHENG Pan, et al. Development status and consideration of carbon fiber stereo puncture prefabrication molding technology[J]. Textile Accessories, 2024, 51(1): 65–68.
[2]  官威, 李文晓, 戴瑛, 等. 纺织复合材料预制体变形研究综述[J]. 航空制造技术, 2021, 64(1/2): 22–37.GUAN Wei, LI Wenxiao, DAI Ying, et al. A review of study on deformation of textile composite preforms[J]. Aeronautical Manufacturing Technology, 2021, 64(1/2): 22–37.
1-2 ]}，是制作高性能防热隔热碳/碳复合材料的优良基材，广泛应用于航空航天、军工等领域^{[ [3]  YAO T L, LI D S, JIANG L. Parametric theoretical prediction of elastic properties and strength of fine weave pierced C/C composites[J]. Composite Structures, 2023, 325: 117614.
[4]  吴小军, 杨杰, 杨云鹏, 等. 细编穿刺C/C喉衬材料烧蚀微结构及损伤机理[J]. 固体火箭技术, 2022, 45(4): 594–600.WU Xiaojun, YANG Jie, YANG Yunpeng, et al. Ablative microstructure and damage mechanism of the fine weave pierced C/C throat material[J]. Journal of Solid Rocket Technology, 2022, 45(4): 594–600.
3-4 ]}。

加压密实工艺是指将预制体沿Z方向压缩到指定的高度，此工艺不仅可以提高预制体的纤维体积含量，还可以影响复合材料成型中的渗碳过程和孔隙结构，是决定复合材料力学性能的重要环节^{[ [5]  葛敬冉, 刘增飞, 乔健伟, 等. 航空复杂结构纤维预制体成型工艺与复合材料性能仿真研究进展[J]. 航空制造技术, 2022, 65(16): 14–30.GE Jingran, LIU Zengfei, QIAO Jianwei, et al. Research progress on molding processes of fiber preforms and performances simulation of composites for aeronautical complex structures[J]. Aeronautical Manufacturing Technology, 2022, 65(16): 14–30.
[6]  梅鸣, 周珺晗, 韦凯. 纤维增强复合材料自动化成型中织物变形研究进展[J]. 复合材料学报, 2023, 40(5): 2507–2524.MEI Ming, ZHOU Junhan, WEI Kai. Advance of fabric deformation in automated forming of fiber reinforced composites[J]. Acta Materiae Compositae Sinica, 2023, 40(5): 2507–2524.
5-6 ]}。为探究预制体压实过程中的变形机制，国内外学者采用数值仿真法揭示织物压缩过程的变形行为。Jaganathan等^{[ [7]  JAGANATHAN S, VAHEDI TAFRESHI H, SHIM E, et al. A study on compression-induced morphological changes of nonwoven fibrous materials[J]. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 2009, 337(1–3): 173–179.
7 ]}提出了一种基于图像的建模技术，研究了非织造纤维材料在压缩载荷下孔径分布的变化情况。Mahadik等^{[ [8]  MAHADIK Y, HALLETT S R. Finite element modelling of tow geometry in 3D woven fabrics[J]. Composites Part A: Applied Science and Manufacturing, 2010, 41(9): 1192–1200.
8 ]}基于动力学理论建立了三维织物压实过程有限元模型，研究了不同压实程度对织物几何结构的影响。Green等^{[ [9]  GREEN S D, LONG A C, EL SAID B S F, et al. Numerical modelling of 3D woven preform deformations[J]. Composite Structures, 2014, 108: 747–756.
9 ]}基于多链数字单元法构建了三维正交机织复合材料有限元模型，分析了压缩过程中纤维的形态变化。Liu等^{[ [10]  LIU Q, LU Z X, ZHU M, et al. Experimental and FEM analysis of the compressive behavior of 3D random fibrous materials with bonded networks[J]. Journal of Materials Science, 2014, 49(3): 1386–1398.
10 ]}采用试验和有限元结合的方法讨论了三维随机材料的压缩性能，揭示了其压缩破坏机理。Daelemans等^{[ [11]  DAELEMANS L, TOMME B, CAGLAR B, et al. Kinematic and mechanical response of dry woven fabrics in through-thickness compression: Virtual fiber modeling with mesh overlay technique and experimental validation[J]. Composites Science and Technology, 2021, 207: 108706.
11 ]}提出了一种结合有限元和虚拟纤维模型的混合虚拟纤维技术，通过ABAQUS软件预测干织物在厚度方向压缩下的力学行为。Mei等^{[ [12]  MEI B L, DONG J Z, REN H Q, et al. Research and modelling of fiber deformation mechanism of 3D four-direction preform under compression loading[J]. Textile Research Journal, 2023, 93(17–18): 4329–4341.
12 ]}在ABAQUS仿真平台上建立了反映三维四向预制体微观结构的数值模型，并构建了纤维数值生长模型用来描述压缩载荷下纤维束截面的变化。Fang等^{[ [13]  FANG J, WANG J J, XIE J B, et al. A microscale modeling method for predicting the compressive behavior of 3D needled nonwoven fiber preforms[J]. Materials & Design, 2024, 243: 113078.
13 ]}提出了一种微尺度建模方法，用于预测三维针刺非织造纤维预制体的压缩行为。

综上所述，国内外学者通过建立预制体数值模型分析二维机织物或非织造纤维的压缩性能，但整体穿刺碳纤维毡预制体较上述织物结构复杂，对其压实特性研究鲜有报道。本文利用梁单元对虚拟纤维进行建模，基于碳纤维毡微观结构参数化建模方法构建其数值模型，并在ABAQUS软件中赋予梁单元实际的碳纤维力学性能。进一步，基于整体穿刺预制体压实工艺建立了含Z方向钢针的碳纤维毡压实模型，对模型压实过程进行有限元仿真并观察预制体的结构变化。同时，采用万能试验机对预制体进行压实试验以验证压实模型的正确性。

整体穿刺碳纤维毡预制体制备过程包括钢针布放、整体穿刺、加压密实、钢针置换等阶段。将碳纤维毡裁剪叠层后铺放到钢针阵列上端，在穿刺模板的推动下沿钢针下移进行整体穿刺并加压密实，重复该过程至设定的高度后，再由碳纤维逐一替代钢针，制成整体穿刺预制体^{[ [14]  杨景朝, 蒋秀明, 董九志, 等. 基于机器学习的整体穿刺加压参数预测方法[J]. 纺织学报, 2019, 40(8): 157–163.YANG Jingzhao, JIANG Xiuming, DONG Jiuzhi, et al. Prediction method of integrated piercing pressure parameters based on machine learning[J]. Journal of Textile Research, 2019, 40(8): 157–163.
14 ]}。

加压密实是整体穿刺预制体制备的关键过程，在预制体制备过程中，每铺放5层碳纤维毡，就需要采用自制的压实工装对预制体进行一次压实，以提高预制体的纤维体积分数和层间密度，压实原理如图1所示。在穿刺模板的推动下，预制体厚度沿Z方向逐渐减小，直到预制体厚度不再发生变化，即完成加压密实工艺。

采用康本HT–1 3K系列碳纤维丝束来梳理面密度为25 g/m²的碳纤维毡，制备出4组不同尺寸的整体穿刺预制体，具体工艺参数见表1。

利用SHIMADZU AGS–X型号万能拉伸实验机和自制的整体穿刺预制体压实工装进行压实试验，试验装置如图2所示。压实试验具体过程为：将穿刺完成的碳纤维毡及穿刺工装夹具放置在万能拉伸机中间，采用等速法进行压实。压实阶段压头以1 mm/min的恒定速度对样件进行压实，直到预制体高度不再变化时进行载荷释放，压头以10 mm/min的速度回到原位，通过试验设备记录预制体在压实过程中的位移变化曲线。对1# ~ 4#试样分别进行试验，每组试验各进行3次后取平均值，并绘制试样的压实载荷–位移均值曲线。

碳纤维毡是由大量随机取向的碳纤维通过气流成网形成的一种轻质多孔材料，精确地控制预制体的几何尺寸和内部纤维结构是进行有限元分析的前提。从微观尺度看，碳纤维毡内部纤维相互穿插、随机分布，想要建立与实际完全相同的纤维几何模型非常困难。Xie等^{[ [15]  XIE J B, FANG J, CHEN L, et al. Micro-scale modeling of 3D needled nonwoven fiber preforms[J]. Composite Structures, 2022, 281: 114995.
15 ]}利用Micro-CT技术研究了纤维毡的微观结构，对纤维的直径、取向和聚集度等几何参数进行了统计分析，确定纤维的初始取向角分布区间为[0，π]。本文提出了基于碳纤维毡微观结构的参数化建模方法，假设碳纤维毡是由大量形状为圆柱的纤维个体组合而成，在不影响纤维形状和分布的基础上，定义纤维取向分布、空间分布和尺寸分布的微观结构几何参数，在初始取向角分布区间内选择参数，以建立真实反映织物形貌的碳纤维毡模型。

纤维的分布形式主要取决于纤维的位置和取向，纤维的位置可由纤维端点确定，纤维的取向以线段方位角表示^{[ [16]  李志男, 杨振宇, 卢子兴. 随机纤维网络材料力学性能研究进展[J]. 航空制造技术, 2023, 66(15): 71–95.LI Zhinan, YANG Zhenyu, LU Zixing. Research progress on mechanical properties of random fiber networks[J]. Aeronautical Manufacturing Technology, 2023, 66(15): 71–95.
16 ]}。任意取一段长度为l的纤维随机分布在l_x×l_y×l_z的三维空间中，如图3所示。纤维与z轴的夹角为θ，纤维在xy平面上的投影与x轴的夹角为φ，故任意纤维取向可以用 （θ，φ）来定义，其中θ、φ∈[0，π]。本文采用均匀分布函数描述纤维在指定范围内的随机分布^{[ [17]  LI Z Y, LIU Z, XUE Y B, et al. A novel algorithm for significantly increasing the fiber volume fraction in the reconstruction model with large fiber aspect ratio[J]. Journal of Industrial Textiles, 2022, 51(1): 506S–530S.
[18]  LI D T, XIA W, FANG Q Z, et al. Experimental and numerical investigations on the tensile behavior of 3D random fibrous materials at elevated temperature[J]. Composite Structures, 2017, 160: 292–299.
17-18 ]}，纤维取向分布函数Ω（θ，φ）在三维空间中满足

纤维的起始点P₀（x₀，y₀，z₀）和端点P₁（x₁，y₁，z₁）的坐标可由式 （2）和 （3）表示。

式中，Ω（0，l_x）、Ω（0，l_y）和Ω（0，l_z）分别是0到l_x、l_y和l_z范围内的随机分布函数，该函数由式 （1）确定。

将纤维毡中任意两根相邻的纤维AB、CD的分布形式简化为图4。其中，M、N分别为两根纤维上任意点，定义a = AM/AB（0 < a < 1），b = CN/CD（0 < b < 1），故M、N两点的向量可以通过a、b表示，如式 （4）所示。

相邻两根纤维的最短距离为

式中，d_min为两纤维最短距离，mm。纤维毡建模时要避免纤维与纤维之间相互交叉而产生干涉，故纤维在空间分布时相邻两根纤维的最短距离要大于纤维直径。

选取一块面密度为25 g/m²的碳纤维毡，提取碳纤维毡内部的100根纤维，对其长度和直径进行测量统计，得到碳纤维毡的纤维长度和直径，统计结果如图5所示。由图5（a）可知，纤维长度变化具有一定的随机性，长度主要集中在30~50 mm之间 （占47%），使用式 （6）中的Gamma分布来拟合纤维长度分布情况，并定义虚拟纤维长度取值范围为30 ~50 mm。由图5（b）可知，纤维直径范围主要集中在10 ~30 μm之间，纤维直径变异性对碳纤维毡整体性能影响不大^{[ [19]  XIE J B, CHEN X M, ZHANG Y F, et al. Experimental and numerical investigation of the needling process for quartz fibers[J]. Composites Science and Technology, 2018, 165: 115–123.
19 ]}。为减少纤维数量、提高仿真效率，本文假设纤维为直径均匀的圆柱体，对纤维直径统计结果取平均值，确定虚拟纤维的直径为24 μm。

式中，l为纤维长度，mm；α和β分别为形状参数和尺度参数。

基于前文中对碳纤维毡内部纤维各参数的定义和统计结果编写Python脚本程序，并在ABAQUS软件中执行生成纤维毡几何模型。具体步骤为：（1）定义随机纤维分布几何区域大小，生成第一根纤维确定起始点坐标；（2）设置纤维的长度和取向角，判断纤维端点是否超出区域边界；（3）随机定义纤维的移动位移，并生成下一根纤维，判断两根纤维的最短距离是否满足空间分布要求；（4）重复上述过程直到生成最后一根纤维；（5）通过式（7）计算纤维毡模型的面密度是否达到要求 （本文使用1500根虚拟纤维生成面密度为25 g/m²的纤维毡 （图6））。

式中，ρ_s为纤维毡模型的面密度，g/m²；N为虚拟纤维数量，个；R_f为虚拟纤维半径，mm；L_s为虚拟纤维的长度，mm；ρ_f为碳纤维的体积密度，g/cm³；A为纤维毡模型的面积，mm²。

基于ABAQUS仿真平台和Python程序设计语言，建立了整体穿刺碳纤维毡预制体压实模型，对预制体的压实过程进行精确描述。

根据整体穿刺预制体压实工艺和结构特点，本文采用代表性体积单元 （RVE）策略构建有限元模型以提高计算精度和效率。选取3×3钢针阵列和3层面积为10 mm×10 mm的碳纤维毡以及孔间距为2.3 mm的穿刺模板构建的含Z方向钢针阵列的整体穿刺碳纤维毡模型作为RVE单胞，该单胞尺寸的确定考虑了纤维毡内部纤维分布情况和层间的相互作用。

整体穿刺碳纤维毡压实过程仿真如下： （1）对虚拟纤维毡、钢针和穿刺模板进行材料参数赋值及网格划分，其中钢针和穿刺模板采用C3D8R实体单元进行网格划分，碳纤维则采用B31梁单元进行网格划分，通过对碳纤维网格密度的收敛性分析，多次运行模拟并确定最佳网格密度，具体材料类型和网格大小如表2所示；（2）采用通用接触来设置纤维之间及纤维与穿刺模板之间的接触，定义纤维之间摩擦系数μ = 0.3^{[ [19]  XIE J B, CHEN X M, ZHANG Y F, et al. Experimental and numerical investigation of the needling process for quartz fibers[J]. Composites Science and Technology, 2018, 165: 115–123.
19 ]}，设置法向接触为硬接触； （3）约束上压板的水平自由度，约束底板和钢针的所有自由度，赋予上穿刺模板向下的位移边界条件，使纤维毡施加压缩载荷； （4）通过ABAQUS软件中显示动力学算法模拟穿刺纤维毡压实过程，设置分析步时间为10 s； （5）在上穿刺模板的中心位置设置参考点RP，将该参考点与上穿刺模板进行耦合约束，并选择Amp振幅曲线加载位移载荷，即可得到压缩载荷与预制体压实高度曲线。在碳纤维毡压缩过程中，该参考点所记录的反作用力就是纤维毡所承受的压缩载荷，上穿刺模板的位移即为纤维毡受力的压缩位移。

基于上述穿刺碳纤维毡压实过程有限元仿真，记录不同时刻碳纤维毡的位移云图和应力云图如图7所示，可以观察预制体在压实过程中不同阶段的结构变化。压缩前，纤维在自身重力作用下分布松散，相邻纤维之间的距离较大，大部分纤维处于伸直状态，层间纤维并未完全接触且存在大量孔隙，如图7（a）所示；施加压缩载荷后 （图7（b）），纤维毡上部的纤维最先受到挤压，随着穿刺模板位移的增加，相邻的纤维相互接触，受挤压的纤维越来越多，甚至部分纤维发生弯曲变形，纤维毡也由蓬松状态逐渐变得密实，纤维毡高度沿载荷方向不断减小，纤维之间孔隙不断被填充；如图7（c）所示，当纤维毡高度不再随压缩载荷增加而发生变化时，即完成预制体压实工艺，此时预制体的纤维体积分数最大，从图7（c）的应力云图可以看出，纤维接触区域出现了明显的应力集中。

为了更接近真实碳纤维毡的几何结构，在数值模拟中通过调整纤维取向角的取值范围来调整压缩载荷与纤维形变量的关系，采用上述压实模型对1.3节中4组压实预制体的过程进行有限元模拟，预测其压缩载荷与压实高度映射关系，4组穿刺碳纤维毡预制体的仿真和试验压实载荷–位移曲线如图8所示。

从图8中可以看出，在压实载荷作用下，预制体压缩位移变化过程可分为3个阶段，第1阶段为预制体压实初始阶段，较小的加压载荷使预制体迅速压缩，纤维层间孔隙很快被纤维填充，从宏观上看碳纤维毡由松散状态逐渐压实，曲线ab段压实载荷与压缩位移近似呈线性阶段；第2阶段为非线性过渡阶段，当压实载荷小于纤维之间的摩擦力时，纤维单丝之间发生接触，当压实载荷大于纤维之间摩擦力时，纤维会发生滑移甚至产生弯曲变形进一步填充孔隙，曲线bc段压缩位移呈非线性趋势增加，随着压实载荷增加预制体越难被压缩，故曲线斜率逐渐减小；第3阶段为紧密压实阶段，该阶段压实载荷不断增加，预制体高度变化不大，表明此阶段已完成预制体压实工艺。

通过对比4组预制体压实结果可知，碳纤维毡层数越多，层间纤维分布越松散，压实初期预制体压缩位移变化越大，预制体1# ~ 4#的实际压实载荷分别为7078.069 N、8430.799 N、14988.466 N、17136.224 N，实际压缩位移分别为7.925 mm、13.731 mm、7.504 mm、14.897 mm。在碳纤维毡层数相同的条件下，预制体压实后高度近似相同，所需的压实载荷与碳纤维毡表面积成正比，故预制体压实高度与碳纤维毡铺层数有关，与表面积关系不大。

在数值仿真中，当纤维取向角θ，φ∈[π/8，π/2]时得到的压缩载荷与预制体高度变化关系更接近实际情况，且具有相同的压缩载荷–位移非线性变化趋势，两者误差不超过5.5%，验证了仿真模型的正确性。但实际试验中，碳纤维毡中的纤维存在交叉、偏转或弯曲现象，因在模型中忽略了碳纤维毡穿刺后纤维会发生绕针弯曲现象以及纤维直径的变异性，导致实际预制体的压实位移小于数值模拟预测值。在后续工作中，将针对碳纤维毡在穿刺过程中的损伤行为和绕针弯曲现象进行详细研究，进一步提高模型的准确度和适用性。

本文提出了碳纤维毡微观结构参数化建模方法，构造了碳纤维毡数值模型，并对穿刺碳纤维毡压实过程进行有限元模拟以观察预制体在压实过程中的微观形貌，并得出以下2个结论。

（1）考虑了碳纤维毡中纤维取向、长度、位置等具有一定的随机性，对纤维的几何参数进行定义以建立更加真实的碳纤维毡几何模型。基于整体穿刺预制体压实工艺构建了纤维随机分布的穿刺碳纤维毡压实模型。

（2）基于ABAQUS有限元数值分析平台，采用Abaqus/Explicit算法对整体穿刺纤维毡压实过程进行数值模拟，讨论了预制体施加压缩载荷前后的结构变化，从宏观上看，碳纤维毡在压实过程中由蓬松状态逐渐变得密实。此外，本文研究了压缩载荷与预制体压实高度的映射关系，当模拟仿真中的纤维取向角θ，φ∈[π/8，π/2]时与试验结果最吻合，两者之间最大误差<5.5%。证明了该建模方法可用于预测不同面密度穿刺碳纤维毡预制体的压缩载荷–位移响应，为后续预制体力学性能计算和预测提供了基础。