连续纤维增强复合材料（Continuous fiber reinforced composite，CFRC）具有出色的比强度、比刚度等力学性能，已在多个行业得到广泛应用。传统CFRC的加工工艺如模压成型、树脂转移成型（Resin transfer molding，RTM）等，需要借助预制模具并需要在热烘箱中进行固化^[1–3]，因此其几何复杂度和结构设计自由度受到了很大的限制。同时在加工含多孔复合材料结构时，传统工艺涉及成型、装配、蚀刻和后处理等多个步骤，制造过程复杂耗时，成本较高^[4–5]。相较之下，增材制造技术作为一种自由成型制造工艺，不仅能节省成本，还可提供更快速、更灵活的制造方式。因此，利用增材制造技术制备CFRC，能够充分发挥设计和制造的灵活性^[6–8]。

随着增材制造技术的不断进步，开发具有轻量化和先进功能性能的复合材料结构的需求愈发迫切，性能驱动设计逐渐成为一个新兴的研究领域^[9–12]。CFRC的轻量化设计需要充分利用增材制造技术中纤维的空间排列潜力。调控复合材料结构中纤维方向和纤维填充密度的工艺规划策略，对于从纤维微观性能到三维结构拓扑性能的驱动设计至关重要^[13]，是传统成型或层压工艺无法实现的。

具有原位浸渍功能的熔融沉积成型（Fused deposition modeling，FDM）是一种简便的增材制造工艺，可用于制造具有数字化编程路径的CFRC，实现性能驱动的纤维排列设计^[14–15]。三维模型的路径规划通常涉及切片、路径填充设计及工艺参数匹配，除了适当的切片算法外，还可通过设计独特的填充形态以响应外部载荷^[16–17]。

在路径填充设计方面，通用的方法是通过填充区域内轮廓的等高线生成打印路径，以保证纤维/基体按照结构的几何特征分布。Le Duigou等^[18]的研究表明，通过控制每层的纤维体积含量和厚度可以实现对复合材料力学性能的调控。还有一些研究集中在平面填充区域的轮廓和内部填充模式上，Sugiyama等^[19]提出了一种新的连续纤维填充规划策略，通过改变连续纤维的局部体积比可以打印出具有可变刚度的CFRC结构。Pedersen^[20–21]的研究从理论方面证实了正交各向异性的复合材料可通过与主应力方向对齐而达到最佳的力学性能。Li等^[22]开发了利用孔洞和拓扑结构填充工程组件的多种路径规划方法。增材制造工艺虽然可以很好地用于制造具有不同填充模式的CFRC结构，然而目前尚缺乏对填充路径进行数字化调整的参数化设计策略。

针对高效路径规划方法的研究仍处于起步阶段，这使得通过增材制造工艺加工的CFRC结构仅能具备简单的几何形状，同时被限制在有限区域内^[23]。纤维路径的填充方法不仅影响结构的几何形状，而且决定了具有不同纤维取向的打印结构件的机械性能^[24–25]。目前，最常用的3D打印连续纤维复合材料结构的路径填充方法有格栅路径^[26]、ZigZag路径^[27]、希尔伯特曲线路径、螺旋曲线路径、轮廓等距路径和混合路径^[28]等，上述路径填充方法具有简单和高适用性的特点。然而，CFRC的增材制造工艺在制造复杂的结构部件时需要连续且平滑的打印路径，避免纤维发生磨损、断裂等从而影响打印质量。

近期，有研究人员结合各向异性复合材料提出了一种以刚度为主导的结构拓扑优化框架。Wang等^[29]提出了一种应力矢量追踪算法，将固体正交各向异性材料惩罚法（Solid orthotropic material with penalization，SOMP）与拓扑优化方法相结合，创建了连续纤维的载荷依赖路径。Papapetrou^[30]和Fernandes ^[31]等的研究发现，拓扑结构的SOMP方法和水平集方法可以进一步与流线提取、轮廓偏移和等间距等形态学方法相结合，这些平面上的两步顺序优化法可为复合材料增材制造工艺生成理想的纤维打印路径。上述路径规划方法主要实现了紧凑的填充线来构建拓扑结构，并应用单一矢量场来引导纤维对准和定向。然而，这些方法并不适用于周期性或非周期性的多孔填充路径规划，因此需要开发一种具有随形和可调形态、密度的自适应填充路径设计方法。

本文提出了一种基于应力场的连续纤维路径规划方法，使用波投影函数（Wave projection function，WPF）将应力场映射为填充模式层，局部纤维间距和纤维取向分别与应力大小和方向对应。然后从波投影函数的峰值中提取出单个纤维轨迹，填充纤维路径的连接顺序建模为最小化纤维间断与路径总长度的旅行商（Traveling salesman problem，TSP）问题，生成连续打印路径并匹配合适的工艺参数，从而实现CFRC的逐层制造过程。

1　基于波投影的连续填充路径设计方法

在给定载荷条件下，通过有限元法提取主应力场，并利用相应的主应力场方向来控制纤维方向。采用局部纤维取向与局部应力主方向对齐，是一种提高CFRC增材制造结构力学性能的启发式策略。为获取各单元主应力方向并确定结构的传力路径，本文采用线弹性各向同性本构模型作为材料输入^[32]，借助现有商用有限元分析软件获得平滑的主应力方向场，作为路径投影的初始输入。以经典米歇尔梁为例，在有限元分析软件中对设计域进行网格单元划分，添加相应的载荷和约束，获得各单元中心位置与全局坐标系下的应力分量（图1）。其中，F为集中力；L为梁的长度；H为梁的高度。通过后处理可以计算结构主应力场f，提取得到拉伸为主的主应力场f₁和压缩为主的主应力场f₂。

图1　基于有限元法的主应力场提取

Fig.1　Principal stress field extraction based on finite element method

1.1　基于波投影函数的应力路径生成法

在获得图1所示的方向场后，基于此计算来投影函数。投影纤维路径的图形周期阵列形式选用平面余弦函数，该映射纤维路径的周期阵列投影方法名为波投影法。引入波投影法旨在建立起方向场和相场（亦称标量场）的关系，从而对矢量场进行映射，并表示出纤维轨迹。周期性波传播的方向与其振动方向垂直，纤维束将沿着余弦波的峰谱线排列，波长表示纤维间的距离。因此，可以利用波函数的局部周期性和全局周期性分别控制纤维间距和填充密度。Rumpf等^[33]提出的综合算法揭示了相函数（亦称映射关系函数）Φ（r）与空间变化方向向量场K（r）之间的关系。相函数Φ（r）可以表示为

令$r={[xy]}^{\text{T}}$，$K(r)=[\frac{1}{2}{\frac{1}{2}]}^{\text{T}}$，则由式（1）可知，$\Phi (r)=\Phi (x,y)=\frac{1}{2}$（x+y），向量场K（r）则为相函数Φ（r）的梯度，利用相函数可以构造波投影函数Γ（x，y），其解析式如式（2）所示。

则波投影函数Γ（x，y）顶峰区域的分布形式和相函数Φ（x，y）等值线的分布形式是一致的，并且波投影函数Γ（x，y）的分布也可以投射梯度向量场K（r）。以余弦函数为基础构造的波投影函数可以作为填充图像模式来代表纤维，如图2所示。

图2　基于余弦函数构造的波投影函数示意图

Fig.2　Schematic diagram of constructed wave projection function based on cosine function

为了得到可以映射设计域内主应力方向场的填充图案，波投影函数可以写作如下形式。

式中，K（r）为波投影函数的传播方向向量场，与纤维的排布方向向量场f（r）互相垂直，二者的定义如下。

同时，定义周期函数为

式中，P为周期函数；d为间距因子（Spacing factor, SF）；在不同的角度θ和间距因子d下，波投影函数可以控制纤维间距和填充密度。

然而，当向量场K（r）随着空间变化时，式（3）不再适用。为了在设计域内得到符合给定主应力方向分布的填充模式，需要通过映射相函数Φ（r）来匹配空间变化的向量场K（r），即

则波投影函数可以重写为

可以使用多种解析或数值方法来求解式（7），本文通过最小化二范数来获得最小二乘解。然而由于最小二乘效应，$\nabla \Phi (r)$与K（r）之间依然存在差异，这部分差异不可忽略，尤其是在高应力区域，将会降低结构整体的承载性能。为了减少这方面的影响，引入带约束的最小二乘求解式，通过基于应力的排序函数对设计空间进行划分，利用惩罚放松的数学处理方法缓解由最小二乘解带来的图像方向与应力方向不匹配的问题，最终在高应力区域实现图像方向与应力方向对齐，其表达式如下。

式中，S（r）为基于排序处理的应力场；Ω_u和Ω_r为对应分类的子区域；N为单元的总数量。

联合式（8）和（11）即可得到应力场波投影函数Γ（r），如图3所示。

图3　拉伸、压缩主应力场的波投影函数示意图

Fig.3　Schematic diagram of wave projection function for tensile and compressive principal stress field

1.2　连续路径生成方法

从纤维路径投影图像到打印路径生成还需要进一步的处理。在提取连续纤维路径轨迹的过程中，图像分割是很重要的一个步骤，基于图像二值化的阈值分割技术已经广泛应用于岩石显微结构图像识别等领域。阈值的选择通常有两种方法，一种是基于灰度直方图波谷的手动选择，另一种是传统的阈值自动选择技术，如Otsu方法、最大熵方法和模糊聚类方法^[34–36]。本文需要阈值过滤来实现二值化，考虑到波投影函数的简谐性，采取直接过滤是最有效的方法。通过阈值T_r来实现对波投影函数的二值化，即

阈值T_r取为0.85（通常为0.8~0.9），波投影函数的二值化处理过程如图4所示。在对图像进行二值化处理得到二值投影函数Γ_r（r）后，可以得到沿主应力分布的带状图形，每一个带状图形的骨架线（中轴线）是先前波投影函数波峰区域所代表的纤维轨迹。随后采用骨架线提取算法对二值投影函数进行骨架线提取，得到纤维轨迹。

图4　波投影函数二值化处理

Fig.4　Binarization processing of wave projection function

在实施具体打印时，采取的策略是将两个主应力主导层进行交替打印。为获得连续的纤维路径以保证结构的力学性能，需将两条不同填充路径的端点沿设计区域边界进行连接，从而减少纤维连续打印断点的产生，边界上的连接路径即为设计域的外轮廓。需要注意的是，边界上的某一段可能会被多次打印。在这种情况下，通过设置适当的路径偏移可以避免发生过度重叠。此时，虽然可以通过上述方法来消除打印间断，但是不同的打印路径顺序会产生不同的总长度。因此，本文选择使用模拟退火算法进行求解，该算法对打印总距离和连续打印路径迭代过程的优化如图5所示。最终优化生成的打印路径序列可以保证具有最短的路径总长度。

图5　连续打印路径的迭代优化过程

Fig.5　Iterative optimization process of continuous printing path

2　结果与讨论

2.1　试验材料与工艺参数

连续纤维选用型号为HTA40–E15–1K的碳纤维束，1束碳纤维中含1000根碳纤维（日本东邦公司）；基体和支撑材料分别为聚乳酸（PLA）丝材（深圳拓竹科技有限公司）和水溶性聚乙烯醇（PVA）丝材（珠海市三绿实业有限公司）。试验采用自主研发的多自由度机器人辅助增材制造系统，采用图6所示的原位浸渍逐层熔融沉积成型工艺进行材料的制备，喷嘴直径为1.2 mm。

图6　原位浸渍逐层熔融沉积成型工艺

Fig.6　Forming process of in-situ impregnation layer-by-layer fused deposition

在执行打印的过程中，连续纤维受到打印平台上已固化树脂与打印喷嘴之间的拉力作用被挤出。当固化树脂与打印喷嘴之间的纤维张力方向与已固化纤维的方向不平行时，纤维从尚未完全固化的树脂中被拉出，导致偏离其原始路径。为保证打印过程中纤维轨迹的准确性，喷嘴经过含有转弯的路径时，需要降低喷头运动速度为树脂固化提供足够的时间。但降低全局打印速度会降低制造效率，因此，结合实际打印时喷头运动速度、路径偏转角处曲率对纤维偏离路径的影响程度，通过计算生成路径的局部几何特征（路径中相邻两个方向向量/3个打印点之间的夹角）对工艺参数进行调整，使喷头运动速度v_p与路径局部偏转角α相匹配。偏转角范围与喷头运动速度的匹配规则如表1所示。

表1　偏转角范围与喷头运动速度的匹配规则

Table 1　Matching rules for deflection angle range and nozzle velocity

2.2　填充密度对承载性能的影响

在填充模式设计过程中，填充密度是影响复合材料结构力学性能的重要因素之一。波投影函数Γ（r）的填充形态与周期函数P和向量场K（r）有关，应力驱动投影的填充策略能够基于间距因子d调控填充密度。间距系数是由设计者人为调整的设计变量，在具体路径填充设计过程中，可根据设计需求，考虑最小许用线间距、承载能力、填充质量等因素选择合适的间距因子d，从而获得最佳的轻量化结构。为探究填充密度对结构承载性能的影响，选取典型悬臂梁模型为研究对象。考虑纤维路径的最小许用线间距为1.2 mm，分别设定间距因子d为4、6、8，以调控悬臂梁设计区域内的路径填充密度（图7）。CFRC增材制造系统分别根据3种设定的间距因子打印出3种不同填充率的悬臂梁填充结构。采用填充率ρ（打印结构占设计空间的体积比例）衡量基于波投影函数生成路径的填充密度，具体定义如下。

式中，w为路径成型宽度；t为打印层厚；n为所有轨迹数；L_i为第i个路径的长度；V为设计空间的总体积。

图7　典型悬臂梁的填充设计与制造

Fig.7　Design and manufacturing of infill for typical cantilever beam

可以看出，3种间距因子悬臂梁结构对应的填充率分别为61.14%、43.42%和32.30%，即随着间距因子的增大，相邻纤维路径之间的距离增大，结构的填充率相应减小。模拟图7（a）所示的悬臂梁工况对3种不同填充率结构进行准静态加载试验，试验结果以载荷–位移曲线表示，如图8所示。在初始加载阶段，3条载荷–位移曲线斜率差别不大，表明3种试件结构刚度相近，这是因为此时结构变形程度较小，主要由试件的矩形轮廓结构承受外部载荷。随着竖向载荷的增加，试件加载处变形逐渐增大，此时由矩形轮廓结构与内部填充结构共同承受外部载荷，设定不同的间距因子使得3种试件具有不同的填充密度，因而试件结构刚度表现出较大的差异。同时，由载荷–位移曲线可知，随着填充密度的增大，试件的承载极限逐渐增大。悬臂梁试件准静态加载试验表明，应力驱动的波投影法可以通过设定不同的间距因子设计不同的填充密度（以填充率表示），进而有效调控CFRC结构的承载性能。

图8　不同填充率悬臂梁的载荷–位移曲线

Fig.8　Force–displacement curves of cantilever beams with different filling rates

2.3　二维应用案例

2.3.1　非均匀面外载荷卫星舱板

选取承受面外非均匀载荷的卫星舱板作为研究对象，载荷分布在面板上的4个局部区域，面板的四角为固定约束，对面板进行加筋设计，如图9所示。对非均匀面外载荷卫星舱板（图9（a））进行有限元分析，提取其主应力场作为波投影法的输入，得到的波投影函数经过迹线提取、最小化迹线连接等步骤生成连续的加筋填充路径。选取间距因子d为8，匹配合适的工艺参数后作为连续CFRC增材制造系统的输入路径（图9（b））。路径成型宽度为1.2 mm，打印层厚度为0.4 mm，加筋厚度为4.8 mm，打印实物如图9（c）所示。由式（13）计算得到加筋的填充率为30.50%，最大和最小纤维间距分别为21 mm和6 mm。这种面板加筋的拓扑填充结构，是应力驱动的复合材料结构设计制造一体化应用于平面结构的一次深入探索，满足了CFRC结构的定制化设计制造需求。

图9　非均匀面外载荷的面板加筋设计与制造

Fig.9　Design and manufacturing of stiffened panel under non-uniform out-of-plane load

2.3.2　多载荷工况无人机翼段

选取具有典型曲面特征的多载荷工况无人机翼段结构作为研究对象，进行填充设计与制造（图10）。该结构的工况为升力方向的均布载荷和机翼中部后缘的集中力，翼根部分为固定约束，如图10（a）所示。然而，对于机翼这种曲面结构，无法直接使用基于平面余弦函数的波投影法进行主应力迹线生成，需要先进行曲面切片并将分割后的主应力场映射至二维平面，本文以上翼面作为保形面进行曲面切片。将映射至平面的应力场作为波投影法的输入，生成平面连续纤维路径，选取间距因子d为6。随后进行逆映射变换至原曲面，匹配合适的工艺参数，最终获得曲面连续纤维路径，如图10（b）所示。执行翼段结构的打印之前，以下翼面为参考生成对应的共形支撑结构，支撑部分选用水溶性聚乙烯醇丝材打印以方便后续的去除，最终得到填充率为49.43%的CFRC翼段填充结构（图10（c））。在考虑机翼多载荷工况的同时，须兼顾曲面切片、支撑约束、防喷头碰撞干涉、支撑结构与目标结构交替打印等技术难点。将应力驱动拓扑填充设计方法与增材制造技术相结合并应用于三维无人机翼段结构，为复合材料结构的轻质承载一体化设计制造提供了新思路。

图10　多载荷工况的无人机翼段填充设计与制造

Fig.10　Infill design and manufacturing of UAV wing segment under multi-load conditions

3　结论

（1）引入波投影函数和基于应力的惩罚约束，实现了连续纤维方向与局部主应力方向的对齐，同时将填充间距与局部主应力大小相映射。采用基于旅行商（Traveling salesman problem，TSP）的路径规划算法，实现了最少纤维间断与最短路径总长度的打印轨迹生成。通过控制填充密度设计参数（间距因子），可设计具有不同承载能力的填充结构。

（2）借助多自由度机器人辅助增材制造平台，完成了卫星舱板加筋设计与无人机翼段填充设计的实物打印，证实了该方法在自适应填充设计制造方面的应用潜力与优势。