3D打印技术在建筑、航空航天和生物医学等领域得到广泛应用^[1-3]. 传统3D打印采用三轴设备，在3D打印过程中须添加支撑结构，不仅会造成材料的浪费，还会增加支撑结构拆除的难度. 在金属3D打印中，拆除支撑结构还可能会损坏构件^[4]. 与三轴方法相比，多轴3D打印引入旋转打印平台，在打印过程中打印方向可动态调整，制造灵活性更高，有助于缓解甚至解决悬垂效应带来的支撑问题. 多轴3D打印领域的研究主要集中在具有5个或6个自由度(DOF)的打印系统上^[5]. Li等^[6]提出基于向量场的曲面层切片和打印路径规划方法，用来消除支撑结构、防止碰撞并提高打印构件的机械性能. Dai等^[7]采用两阶段分解方法克服多轴3D打印中的悬垂问题：将零件分解为曲面打印层后，进一步分解为曲线打印路径. Xu等^[8-9]也提出组件分解打印方法，将给定的结构分解为多个单独的部分，通过为各部分指定不同的打印方向来避免悬垂问题，但该方法未考虑打印碰撞问题. 上述研究主要关注打印路径规划以及减少使用支撑结构和打印时间等方面^[10-11]，研究多针对已知形状的结构，多轴3D打印与结构设计优化相结合的研究鲜见.

结构优化领域的理论基础可追溯至Michell的开创性工作^[12]. 该领域研究主要分为离散桁架优化与连续体拓扑优化. 对于连续体拓扑优化方法，Bendsøe等^[13]提出通过计算在均质同性材料中引入周期性分布小孔的极小值来构造各向异性材料在空间中的最优分布，推动了连续体拓扑优化方法的发展. 后来逐渐发展出各种优化分支，如双向渐进结构优化法(bidirectional evolutionary structural optimization, BESO)^[14]、水平集法^[15]、移动可变形组件法(moving morphable component, MMC)^[16-17]和可移动变形孔洞法(moving morphable voids, MMVs)^[18-19]等连续体拓扑优化方法. 对于离散桁架优化方法，Dorn^[20]提出基结构法并对桁架结构进行优化. Gilbert等^[21]在基结构法的基础上提出的构件添加法，显著降低了计算成本且得到最优解. He等^[22]引入几何优化方法进一步简化了结构布局. 在实际应用中，连续体优化方法更适于处理高体积比(优化结构与设计域的比值)问题，离散桁架优化更适于处理低体积比问题. 考虑到悬垂效应在低体积比问题中更为明显^[23]，本研究主要关注离散桁架优化方法.

在连续体拓扑优化中，基于密度解决悬垂问题是常用的方法. 根据附近元素的存在状况来评估单个元素的可打印性，Gaynor等^[24]提出楔形过滤器，用于设计由常规网格组成的自支撑结构，Garaigordobil等^[25]使用边缘检测方法识别和消除优化结构中的悬垂特征. 还有研究者采用前沿传播技术^[26-27]和水平集方法^[28-29]来解决悬垂问题. 在离散桁架优化领域，He等^[30]在传统布局和几何桁架优化方法中添加硬性角度约束，实现了桁架结构自支撑. 林晓阳等^[31-32]利用结构布局、几何优化算法，考虑力学性能与3D打印可制造型的优化方法，生成三维桁架自支撑优化结构，开展了实际模型打印测试. 上述方法虽可有效得到自支撑结构，但主要面向三轴3D打印，考虑悬垂约束还会明显增加优化结构的体积. 在结构优化与多轴3D打印相结合方面，Lu等^[33]将设计域分解为多个分区，再利用优化方法确定每个分区的局部打印方向. 该方法适用于二维桁架结构，在解决三维问题时会出现数值计算困难的情况. 本研究将问题扩展到三维桁架结构，针对多轴3D打印中三维自支撑桁架结构优化存在数值稳定性差、可打印性指标低的问题，提出基于多轴3D打印的两步优化方法.

为了克服多轴3D打印中的悬垂效应，优化生成自支撑结构. 从传统的布局优化中获得下界解，再考虑悬垂约束，步骤如下. 1)识别结构的打印曲面：使用新的优化方法来识别给定结构的打印曲面，使结构可打印比例(无悬垂体积/总体积)最大化. 步骤1）中结构保持不变，无法保证结构被完全打印. 当结构可打印比例低于设定阈值时，执行步骤2). 2)结构和打印曲面整体优化：以步骤1)的结果为起点，同时对结构和打印曲面进行优化. 基于几何优化的问题，进一步提高结构的可打印比率和结构性能.

通过布局优化方法，得到特定设计条件(如设计域、载荷和边界约束)下高效的结构布局. 如图1所示，该过程分3步：指定设计域、载荷和边界条件，生成节点和潜在构件以形成基结构，通过求解优化问题来确定最佳布局. 结构布局优化数学模型为

(1a)$ \mathop {\min }\limits_{a,q} {{ V = }}{{\boldsymbol{l}}^{\text{T}}}{\boldsymbol{a}} .$

(1b)$ {\text{s}}{\text{.t}}{\text{. }}{\boldsymbol{B}}{{\boldsymbol{q}}_\alpha }{\text{ = }}{{\boldsymbol{f}}_\alpha } , $

(1c)$ - {\sigma ^ - }{\boldsymbol{a}} \leqslant {{\boldsymbol{q}}_\alpha } \leqslant {\sigma ^+}{\boldsymbol{a}} , $

(1d)$ {{a}}_i \geqslant {\text{0 }} . $

式中：$ \boldsymbol{a}= [{a}_{1},{a}_{2},\cdots ,{a}_{i},\cdots ,{a}_{m}{]}^{{\mathrm{T}}} $为构件截面积向量，m为构件数量；$ {\boldsymbol{q}}_{\alpha }= [{q}_{1,\alpha },{q}_{2,\alpha },\cdots ,{q}_{m,\alpha }{]}^{{\mathrm{T}}} $为构件在第$ \alpha $种工况下的内力向量，受拉为正，受压为负；$ {{V}} $为结构总体积；$ \boldsymbol{l}= [{l}_{1},{l}_{2},\cdots ,{l}_{m}{]}^{{\mathrm{T}}} $为构件长度向量；$ \boldsymbol{B} $为平衡矩阵；$ {{\boldsymbol{f}}_\alpha } $为第$ \alpha $种工况下的节点载荷向量；$ {\sigma ^ - } $，$ {\sigma^+} $分别为材料在受压和受拉时的极限强度. 布局优化得到的优化结构通常比较复杂，为此将节点坐标视为额外设计变量的几何优化，进一步简化结构^[22]. 结构几何优化数学模型为

(2a)$ \min {{V = }}{\boldsymbol{l}}{({\boldsymbol{x}},{\boldsymbol{y}})^{\text{T}}}{\boldsymbol{a}} . $

(2b)$ {\text{s}}{\text{.t}}{\text{. }}{\boldsymbol{B}}{\text{(}}{\boldsymbol{x}},{\boldsymbol{y}}{\text{)}}{{\boldsymbol{q}}_\alpha }{\text{ = }}{{\boldsymbol{f}}_\alpha } , $

(2c)$ - {\sigma ^ - }{\boldsymbol{a}} \leqslant {{\boldsymbol{q}}_\alpha } \leqslant {\sigma ^+}{\boldsymbol{a}} , $

(2d)$ {{a}}_i \geqslant {\text{0}} . $

式中：$ \boldsymbol{x}= [{x}_{1},{x}_{2},\cdots ,{x}_{d}{]}^{{\mathrm{T}}} $和$ \boldsymbol{y}= [{y}_{1},{y}_{2},\cdots ,{y}_{d}{]}^{{\mathrm{T}}} $分别为节点的横纵坐标向量，$d$为节点数量. 几何优化利用布局优化结果作为初始解，通过在迭代过程中进行如合并节点和创建交叉点的几何操作实现结构布局的简化，更多细节见文献[22].

在二维问题中，将设计域划分为多个矩形分区，假定每个分区对应1个局部打印方向；根据分区内部构件方向确定相应分区的局部打印方向，识别原理如图2所示. 在分区(2,1)中，每根构件都有各自的可打印范围. 当不同构件的可打印范围存在重叠区域时，这些构件能够沿着重叠区域的共同打印方向顺利打印，定义重叠区域为共同打印范围. 若多根构件不存在共同打印范围，则选择可使结构无法打印体积最小化的方向作为最优局部打印方向. 打印方向识别方法的数学模型为

(3a)$ \mathop {\min }\limits_{\eta ,\phi } {\eta _{\text{t}}} = \sum\limits_\mu {\sum\limits_\nu {{\boldsymbol{\eta}} _{_{\mu ,\nu }}^{\text{T}}{{\boldsymbol{v}}_{_{\mu ,\nu }}}} } . $

(3b)$ \begin{split} {\text{s}}{\text{.t}}{\text{. }}{{\boldsymbol{\eta}} _{\mu ,v}} &= \left| {{\phi _{\mu ,v}}{\boldsymbol{e}} - {\theta _{\mu ,v}} - {\phi _{{\text{max}}}}{\boldsymbol{e}}} \right|+\left| {\phi _{\mu ,v}}{\boldsymbol{e}} -\right.\\ & \left.{\theta _{\mu ,v}}+{\phi _{{\text{max}}}}{\boldsymbol{e}} \right| -2{\phi _{{\text{max}}}}{\boldsymbol{e}},\; \mu \in {\bf{Z}_\mu }, \nu \in {\bf{Z}_v} ;\end{split} $

(3c)$\left({\phi _{i,j}} -{\phi _{i,j-1}} \right)\leqslant {\phi _{{\text{t,max}}}}. $

式中：$ {\eta _{\text{t}}} $为总悬垂违反值；${{\boldsymbol{\eta}} _{\mu ,v}}$,${{\boldsymbol{v}}_{\mu ,v}}$分别为每根构件的悬垂违反值和体积，$ \mu $和$ \nu $分别为分区的行、列索引；${\bf{Z}_\mu } = \{ 1,2,\cdots,{n_\mu }\} $和${\bf{Z}}_v = \{ 1,2,\cdots,{n_v}\} $分别为分区的行、列索引集；$ {{{n}}_\mu } $和$ {{{n}}_\nu } $分别为分区的行和列的数量；$ {\phi _{\mu ,v}} $为局部打印方向；$ {\phi _{{\text{max}}}} $为最大悬垂角；$ {\theta _{\mu ,v}} $为构件的方向角；e为0-1映射矩阵，用于定位构件与所属分区；${\phi _{{\text{t,max}}}}$为最大打印旋转角；${{\boldsymbol{\phi}} _{\mathrm{t}}} = [{\phi _{1,2}} - {\phi _{1,1}}, {\phi _{1,3}} - {\phi _{1,2}}, \cdots ,{\phi _{i,j}} - {\phi _{i,j-1}},\cdots, {\phi _{{{{n}}_\mu },{{{n}}_\nu }}} - {\phi _{{{{n}}_\mu },{{{n}}_\nu } - 1}}]^{\text{T}}$为转角向量.

与2D方法类似，三维问题中用切割面将设计域分割成若干立方体分区，如图3所示.

二维问题中的悬垂违反值利用局部打印方向和构件方向计算. 简单起见，在三维问题中，将局部打印方向和构件方向投影到2个相互垂直的平面$ xoz $和$ yoz $上，如图4所示. 约束条件1：投影面的悬垂违反值为

(4a)$ \beta = \left| {{\phi _x} - {\theta _x} - {\phi _{{\text{max}}}}} \right|+\left| {{\phi _x} - {\theta _x}+{\phi _{{\text{max}}}}} \right| - 2{\phi _{{\text{max}}}}, $

(4b)$ \gamma = \left| {{\phi _y} - {\theta _y} - {\phi _{{\text{max}}}}} \right|+\left| {{\phi _y} - {\theta _y}+{\phi _{{\text{max}}}}} \right| - 2{\phi _{{\text{max}}}} .$

式中：$\beta $、$ \gamma $分别为构件在$ xoz $和$ yoz $平面上的投影悬垂违反值；${\phi _x}$、$ {\phi _y} $分别为分区的局部打印方向在平面$ xoz $和$ yoz $上的投影；${\theta _x}$、$ {\theta _y} $分别为构件的方向角在平面$ xoz $和$ yoz $上的投影. 约束条件1量化了构件在平面$ xoz $和$ yoz $上的投影方向与理想打印方向之间的偏差，目的是确保构件的打印方向不会超过打印机的最大悬垂能力，即最大悬垂角. 当构件的方向角${\theta _x}$在${\phi _x}-{\phi _{{\text{max}}}}$与${\phi _x}$+${\phi _{{\text{max}}}}$之间时，$\beta $=0，表示没有悬垂违反；当${\theta _x}$超出这个范围时，$\beta $表示超出的最大悬垂限制的程度. 平面$ yoz $同理. 投影面的悬垂违反值如图5所示. 尽管上述约束条件导致优化结果保守，但可显著降低数值求解难度. 多轴打印具有灵活性，因此对角度采用投影处理带来的影响很小.

与二维问题不同，在三维问题中，相邻的局部打印面须完全连接，以确保整体打印面的连续性，如图6所示. 通过约束条件2来限制局部打印方向：

(5a)$ {\phi _{x,\left( {i,v,w} \right)}} = {\phi _{x,\left( {j,v,w} \right)}}, \;(i,j) \in {\bf{Z}_\mu }, \nu \in {\bf{Z}}_v, w \in {\bf{Z}} _w; $

(5b)$ {\phi _{y,\left( {\mu ,k,w} \right)}} = {\phi _{y,\left( {\mu ,l,w} \right)}}, \;(k,l) \in {\bf{Z}_\mu }, \nu \in {\bf{Z}}_v, w \in {\bf{Z}}_w. $

式中：${\phi _{x,(\mu ,v,w)}}、{\phi _{y,(\mu ,v,w)}}$分别为分区局部打印方向在平面$ xoz $、$ yoz $上的投影，$\mu ,v,w$分别为分区的列、行、层索引. 式（5a）表示在平面$ xoz $上，对于沿x轴相邻的列索引i和j，当位于同一y行和z层时，分区局部打印方向在x轴上的投影分量（沿x轴的投影）必须相等. 相邻分区的打印方向在x轴方向上的一致保持了打印面的连续性. 式（5b）表示在平面$ yoz $上，对于沿y轴相邻的行索引k和l当位于同一x列和z层时，分区局部打印方向在y轴上的投影分量（沿y轴的投影）必须相等. 打印面的连续约束条件确保了在3D打印过程中，局部打印方向在相邻分区之间是连续的. 平面$ xoz $、$ yoz $上投影的一致性可以保证打印面在空间中正确对齐，实现平滑和精确的打印效果.

曲面若有较大旋转角可能会导致打印喷嘴与已打印的部件发生碰撞，为此设置如图7所示的约束条件3，限制相邻分区之间的投影打印旋转角：

(6)$ {{\boldsymbol{\phi}} _{x,{\mathrm{t}}}} \leqslant {\phi _{{\text{t,max}}}}{\boldsymbol{e}},\;\;{{\boldsymbol{\phi}} _{y,{\mathrm{t}}}} \leqslant {\phi _{{\text{t,max}}}}{\boldsymbol{e}}. $

将3D打印面的优化问题表示为

(7a)$ \mathop {\min }\limits_{{\boldsymbol{\eta}} ,{\boldsymbol{\phi}} } {\eta _{\mathrm{t}}} = \sum\limits_{{\mu}} {\sum\limits_v {\sum\limits_w {{{\boldsymbol{\eta}} _{\mu ,v,w}}^{\mathrm{T}}{{\boldsymbol{v}}_{\mu ,v,w}}} } } . $

(7b)$ \begin{split} & {\text{s}}{\text{.t}}{\text{. }}{\overline {\boldsymbol{\eta}} _{\mu ,v,w}} = \left| {{{\boldsymbol{\phi}} _{\mu ,v,w}}{\boldsymbol{e}} - {{\boldsymbol{\theta }}_{\mu ,v,w}} - {\phi _{{\text{max}}}}{\boldsymbol{e}}} \right|+\\& \left| {{{\boldsymbol{\phi}} _{\mu ,v,w}}{\boldsymbol{e}} - {{\boldsymbol{\theta}} _{\mu ,v,w}}+{\phi _{{\text{max}}}}{\boldsymbol{e}}} \right| - 2{\phi _{{\text{max}}}},\;(\mu ,\nu ,w) \in \bf{Z} ;\end{split} $

(7c)$ \begin{split} & {\phi _{x,\left( {i,v,w} \right)}} = {\phi _{x,\left( {j,v,w} \right)}},{\phi _{y,\left( {\mu ,k,w} \right)}} = {\phi _{y,\left( {\mu ,l,w} \right)}},\\& (i,j) \in {\bf{Z}_\mu }, (k,l) \in {\bf{Z}_\nu }, (\mu ,\nu ,w) \in \bf{Z};\end{split} $

(7d)$ \left( {\phi _{i,j}}-{\phi _{i,j-1}} \right)\leqslant {\phi _{{\text{t,max}}}} . $

式中：${{\boldsymbol{\eta}} _{\mu ,v,w}} = {{\boldsymbol{\beta}} _{\mu ,v,w}}+{{\boldsymbol{\gamma}} _{\mu ,v,w}},$ ${\overline {\boldsymbol{\eta}} _{\mu ,v,w}} = \left[{{{\boldsymbol{\beta}} _{\mu ,v,w}}} ,\, {{{\boldsymbol{\gamma}} _{\mu ,v,w}}} \right]^{\mathrm{T}},$ ${{\boldsymbol{\beta}} _{\mu ,v,w}}$和${{\boldsymbol{\gamma}} _{\mu ,v,w}}$分别为构件在平面$ xoz $和$ yoz $上的投影悬垂违反值向量；$ {{\boldsymbol{\phi}} _{\mu ,v,w}} = \left[ {{{\boldsymbol{\phi}} _{x,(\mu ,v,w)}}} ,\, {{{\boldsymbol{\phi}} _{y,(\mu ,v,w)}}} \right] ^{\mathrm{T}}$，$ {{\boldsymbol{\phi}} _{x,(\mu ,v,w)}} $和$ {{\boldsymbol{\phi}} _{y,(\mu ,v,w)}} $分别为构件2个方向的局部打印方向投影向量；$ {{\boldsymbol{\theta}} _{\mu ,v,w}} = \left[{{{\boldsymbol{\theta }}_{x,(\mu ,v,w)}}} ,\, {{{\boldsymbol{\theta }}_{y,(\mu ,v,w)}}} \right] ^{\mathrm{T}}$、$ {{\boldsymbol{\theta }}_{x,(\mu ,v,w)}} $,$ {{\boldsymbol{\theta}} _{y,(\mu ,v,w)}} $分别为构件在平面$ xoz $和$ yoz $上的方向投影向量；${\bf{Z}} = {\bf{Z}_\mu } \times {\bf{Z}_\nu } \times {\bf{Z}}_w$为${\bf{Z}_\mu }$，${\bf{Z}_\nu }$、${\bf{Z}}_w$的笛卡尔积；$ {{\boldsymbol{\phi}} _{x,{\mathrm{t}}}} $、$ {{\boldsymbol{\phi}} _{y,{\mathrm{t}}}} $分别为x、y轴方向上的旋转角向量. 式(7)为线性规划问题，可用内点法求解器^[34]求解，保证优化结果为全局最优.

优化结构是确定的，由式（7）得到的打印面不能保证该结构完全可打印. 为了减少悬垂构件，进行优化问题整合并调整结构和打印面：结合式(2)和式(7)并调整1）由因子$ p $惩罚结构悬垂违反值，获得违反悬垂约束程度最小的结构(约束条件4)：

(8)$ \min \boldsymbol{l}^{\mathrm{T}} {\bf{a}}+p\left(\sum_{{\mu}} \sum_{{v}} \sum_{{w}} \boldsymbol{\eta}_{\mu, v, w}{ }^{\mathrm{T}} \boldsymbol{v}_{\mu, v, w}\right). $

2）每根构件的悬垂违反值根据其最初的分区计算得到，因此添加节点坐标以限制每个节点在设计域最初的分区内移动. 3）将分区边界视为额外的变量，增加整体问题的设计自由度. 如图8所示，在优化过程中，未填充区域会导致规划打印路径时无法完全填充打印层，因此同一行、列排列的分区须共用边界位置变量$ {{\boldsymbol{w}}^x} $和$ {{\boldsymbol{w}}^y} $，排列在同一层的分区共用位置变量$ {{\boldsymbol{w}}^z} $. 4）为了防止分区无限缩小，在调整3）中添加分区边界变量约束条件. 将式(2)、(7)、(8)以及调整3）、4）结合，整体优化模型表示为

(9a)$ \mathop {\min }\limits_{a,q,x,y,z,{\boldsymbol{\eta}} ,{\boldsymbol{\phi}} } {\text{ }}{\boldsymbol{V}}+{\bar {\boldsymbol{\eta}} _{\text{t}}} = {{\boldsymbol{l}}^{\text{T}}}{\boldsymbol{a}}+p\left(\sum\limits_\mu {\sum\limits_v {\sum\limits_w {{{\boldsymbol{\eta}} _{\mu ,v,w}}^{\text{T}}{{\boldsymbol{v}}_{\mu ,v,w}}} } } \right) . $

(9b)$ {\text{s}}{\text{.t}}{\text{. }}{\boldsymbol{B}}{\text{(}}{\boldsymbol{x,y,z}}{\text{)}}{{\boldsymbol{q}}_\alpha }{\text{ = }}{{\boldsymbol{f}}_\alpha } ; $

(9c)$ - {\sigma ^ - }{\boldsymbol{a}} \leqslant {{\boldsymbol{q}}_\alpha } \leqslant {\sigma ^+}{\boldsymbol{a}}; $

(9d)$ {{a_i}} \geqslant {\text{0}} ; $

(9e)$ \begin{split} {\overline {\boldsymbol{\eta}} _{\mu ,v,w}} =& | {{{\boldsymbol{\phi}} _{\mu ,v,w}}{\boldsymbol{e}} - {{\boldsymbol{\theta }}_{\mu ,v,w}} - {\phi _{{\text{max}}}}{\boldsymbol{e}}} |+| {{\boldsymbol{\phi}} _{\mu ,v,w}}{\boldsymbol{e}} - {{\boldsymbol{\theta}} _{\mu ,v,w}}+\\& {\phi _{{\text{max}}}}{\boldsymbol{e}} | - 2{\phi _{{\text{max}}}}, \;(\mu ,\nu ,w) \in \bf{Z} ;\end{split} $

(9f)$ \begin{split} {\phi _{x,\left( {i,v,w} \right)}} = &{\phi _{x,\left( {j,v,w} \right)}},{\phi _{y,\left( {\mu ,k,w} \right)}} = {\phi _{y,\left( {\mu ,l,w} \right)}}, \;(i,j) \in {\bf{Z}_\mu },\\& (k,l) \in {\bf{Z}_\nu }, (\mu ,\nu ,w) \in \bf{Z};\end{split} $

(9g)$\left( {{{\phi}} _{i,j}}-{{{\phi}} _{i,j-1}}\right) \leqslant {\phi _{{\text{t,max}}}} ; $

(9h)$ {{\boldsymbol{w}}_{{\text{upp}}}} - {{\boldsymbol{w}}_{{\text{low}}}} \geqslant {{\boldsymbol{w}}_{{\text{min}}}}; $

(9i)$ {{\boldsymbol{w}}_{{\text{upp}}}} \geqslant {\boldsymbol{X}} \geqslant {{\boldsymbol{w}}_{{\text{low}}}}. $

式中：$ {\bar {\boldsymbol{\eta}} _{\text{t}}} $为优化结构的总悬垂违反值；$ {{\boldsymbol{w}}_{{\text{low}}}} = [{\boldsymbol{w}}_{\mu - 1,v,w}^x, {\boldsymbol{w}}_{\mu ,v - 1,w}^y,{\boldsymbol{w}}_{\mu ,v,w - 1}^z] $和${{\boldsymbol{w}}_{{\text{upp}}}} = [{\boldsymbol{w}}_{\mu ,v,w}^x,{\boldsymbol{w}}_{\mu ,v,w}^y,{\boldsymbol{w}}_{\mu ,v,w}^z]$为边界坐标矩阵；${{\boldsymbol{w}}_{{\text{min}}}} = [{w_{x,{\text{min}}}}{\boldsymbol{e}},{w_{y,{\text{min}}}}{\boldsymbol{e}},{w_{z,{\text{min}}}}{\boldsymbol{e}}]$为立方体分区在x、y和z方向上的长度下限；$ p $为悬垂违反值的惩罚因子. 引入节点坐标变量使式(9)成为为非线性编程问题，本研究使用非线性编程求解器IPOPT处理非线性编程问题^[35].

式(9)是非线性非凸问题，初始解会对最终优化结果的质量造成显著影响^[36]. 本研究提出迭代的半整体方法，得到悬垂构件更少的式(9)初始解，提高最终优化结果的质量. 该方法从求解式(7)得到的弯曲打印面开始. 在已知打印面的情况下，由式(4)计算基结构和潜在构件的投影悬垂违反值$\beta $和$\gamma $，如果构件跨越多个分区，悬垂违反值选取各分区最大值. 将目标函数修改为约束条件5，重新求解布局优化问题式(2)，相应的悬垂值惩罚构件为

(10)$ \mathop {\min }\limits_{a,q} V = \sum\limits_{\boldsymbol{\mu}} {\sum\limits_{\boldsymbol{v}} {\sum\limits_w {{{\left( {{\boldsymbol{e}}+{p_0}{{\boldsymbol{\eta}} _{\mu ,v,w}}} \right)}^{\text{T}}}{{\boldsymbol{l}}_{u,v,w}}{{\boldsymbol{a}}_{\mu ,v,w}}} } } . $

式中：$ {p_0} $为悬垂违反值的惩罚因子；$ {{\boldsymbol{l}}_{u,v,w}} $为各构件长度的对角矩阵. 半整体优化方法得到具有更少悬垂构件的结构，通过式(7)获得进一步改进的打印面. 太多的迭代次数和过大的$ {p_0} $会降低结构性能，本研究设定最大迭代次数为5. 半整体优化问题是线性问题，可以用内点求解器^[34]求解，结合构件添加法^[21]进一步提高求解效率.

本研究提出的优化方案综合了结构设计、打印路径规划和迭代改进3个关键环节，不仅使打印出的结构满足力学性能要求，还能够在多轴3D打印过程中实现自支撑，减少支撑材料的使用. 如图9所示为所提方法的流程图. 1）利用布局优化方法对结构进行初步设计，确定构件的位置和尺寸，以满足设计域、载荷和边界条件的要求. 2）对于非自支撑部分，由式(7)进一步优化打印路径，提高结构的可打印性. 引入迭代的半整体优化方法，在已知打印面的基础上，重新计算基结构和潜在构件的投影悬垂违反值，据此调整结构布局. 4）再次运用式(7)细化打印面，实现更优的打印效果. 5）结合结构优化和打印路径规划的结果，进行整体优化，确保最终的结构设计在力学性能和打印可行性之间达到最佳平衡. 优化流程的设计旨在通过逐步迭代和多阶段优化，实现结构设计的最优化和打印路径的高效规划.

如图10所示，多轴3D打印机由机械臂以及安装在机械臂上的打印喷嘴和旋转打印平台组成. 打印平台可在打印过程中旋转、倾斜至不同角度，从而实现曲面打印. 多轴3D打印增加了局部打印位置上的最大悬垂角，有利于应对悬垂问题.

如图11所示，在三轴3D打印中，打印方向固定，各个切片层与水平面平行，打印喷头与已打印的结构不存在碰撞风险^[8]. 在多轴打印过程中，打印方向可以动态改变，在曲面凹度较大处容易发生碰撞(具有较大的打印旋转角$ {\phi _{\text{t}}} $)，通过限制打印喷头相对于某个参考方向的旋转角度来减少碰撞风险.

引入结构可打印比例（无悬垂体积与总体积的比值），量化结构在不依赖外部支撑条件下的成功打印程度. 总悬垂违反度$ {{\boldsymbol{\eta}} _{{\text{tot}}}} $代表结构的可打印性. 当$ \left| {{\boldsymbol{\eta}} _{{\text{tot}}}} \right|$低于阈值0.1时，认为结构可打印，表示只有当超过90%的结构体积可以无悬垂地打印时，才认为该结构是可打印的. 阈值设定基于对多轴3D打印技术的综合考量，包括打印材料的特性、打印速度、打印精度以及打印过程中支撑结构的可行性. 本研究中的优化算法使用 Python语言实现，算法开发和测试在 Anaconda 环境中进行，在实现过程中，使用多个第三方库，包括但不限于 NumPy、SciPy （用于数值计算），Mosek和IPOPT优化求解器.

惩罚因子$ {p_0} $决定了半整体优化对悬垂构件的惩罚程度. 如图12所示为旋转悬臂梁算例不同惩罚因子下的优化结果. 图中，设定$ {p_0} $∈[0.01, 100.00]，最大悬垂角为50°，竖直分区尺寸S_z=2.0，水平分区尺寸S_x、S_y=1.0. 可以看出，当$ {p_0} $=0.01和0.10时，优化结构中存在不可打印的构件；当$ {p_0} $=10.00或更大时，优化结构的体积显著增加. 在图12中，当$ {p_0} $∈[1.00, 7.00]时，优化过程能够产生结构可打印性与体积效率平衡较好的结果. 原因分析：过小的$ {p_0} $可能无法充分惩罚违反悬垂约束的构件，导致结构无法被打印出来；过大的$ {p_0} $会导致体积大幅增加，类似He等^[30]的硬性约束方法的优化结果. 相比名义优化结构，未经半整体优化获得的优化结构体积增加了20.20%($\Delta V $=20.20%)；相比$ {p_0} $∈[1.00, 7.00]的优化结果，未经半整体优化获得的优化结构的体积更大，这证明了半整体优化过程的作用.

如图13所示为当$ {p_0} $=5.0时，旋转悬臂梁打印切片路径. 优化结果得到打印面，再根据每层打印面对结构进行分层切片，得到打印路径.

如图14所示，选择悬臂梁算例分析最大悬垂角对优化结果的影响，设置$ {p_0} $=4.0，S_x=2.0，S_y、S_z=1.0. 综合考虑多轴3D打印机的技术参数，材料在打印过程中的支撑需求以及结构的力学性能要求，设置$ {\phi _{{\text{max}}}} $=45°、50°、55°和60°，以覆盖从较为保守到较为激进的打印角度，全面评估所提方法的有效性. 考虑悬垂约束的三轴优化方法的优化结果比名义优化结果的体积增加了117.19%($\Delta V $= 117.19%). 考虑悬垂约束后，设置不同最大悬垂角的多轴优化结果均为自支撑结构. 可以看到，随着最大悬垂角的减小，结构体积在增加，构件趋于竖直方向. 相比名义优化结构，设置不同最大悬垂角的多轴优化结构的体积增加均不超过4%，说明多轴优化方法带来了更高的设计灵活性.

如图15所示，分2组研究分区尺寸对优化结果的影响，设置$ {\phi _{{\text{max}}}} $=45°. 1）设置S_x、S_y=1.0， S_z由1.0变化到4.0. 由图可知，在竖直分区尺寸介于1.0与4.0之间，S_z=3.0的布局最简，体积较名义优化结构的增量最低($\Delta V $=0.35%). 在整体优化步骤中，由于结构节点被限制在相应分区内，分区层数的增加会限制整体优化的设计自由度，导致结构复杂且体积较大. 2）设置S_z=1.0，S_x、S_y由0.5变化到4.0. 由图可知，较大的分区水平尺寸会增加结构的复杂度，但水平分区尺寸的变化对结构体积的影响较小.

使用如图16所示的多轴3D打印系统进行优化结构的打印测试. 在多轴3D打印系统中，送丝和加热系统集成于KUKA机器人(KR 10 R1100-2)，熔融沉积成型(fused deposition modeling, FDM)打印喷头挤出线材，旋转打印平台提供2个额外的自由度(旋转和倾斜). 倾斜轴旋转范围为±120°，底座平台半径为200 mm，打印平台表面至倾斜中心轴的距离为170 mm. 在打印过程中，通过KUKA smartPAD控制多轴打印系统. 打印模型的设计域尺寸为60 mm×60 mm×120 mm. PLA打印线材直径为1.75 mm；喷嘴直径为0.6 mm；喷嘴温度设定为210 ℃；温床温度为50 ℃；层高为0.4 mm，打印速度为30 mm/s.

以图12的旋转悬臂梁为例，如图17所示，打印路径的生成包括3个步骤. 1）在Rhino软件^[37]中，基于优化结构^[38]中的节点坐标和构件横截面积生成STL格式的结构模型；2）使用获得的打印面对模型进行切片，生成打印路径轮廓层；3）使用锯齿形图案^[39]填充每层的内部打印路径，设置间距为0.4 mm，填充内缩长度为0.1 mm. 使用不同打印方向，打印过程及最终模型如图18所示. 优化打印面使打印过程无需支撑且不发生碰撞，但打印结构表面存在轻微的缺陷，如图18(e)中标记的C构件处. 缺陷产生的原因：1）同层打印方向的突然变化导致表出现刨削现象^[6]；2）多轴打印系统在多轴运动时存在恒定的机器误差，在90°倾斜角度下测量约为1 mm. 问题1）可以通过调整打印路径实现相邻局部打印方向的平滑过渡来解决^[40]；问题2）可以通过测量不同倾斜角度的运动误差，并在打印路径的生成过程中考虑误差校正解决. 虽然存在误差，但是所提方法成功克服了悬垂效应. 悬垂效应最明显的结构构件是图中标记的构件A和B，其悬垂角为71°. 为此，打印该构件时，旋转平台的倾斜角调整为59.3°. 2根构件在没有任何支撑结构的情况下被成功打印，体现了方法的有效性.

针对多轴3D打印，本研究提出考虑悬垂约束的3D桁架结构优化方法，采用两步优化方法得到优化结构及其弯曲打印面: 第一步识别得到优化结构的理想打印面，第二步采用整体优化同时优化结构及其打印面，提高结构性能并最大限度减少悬垂构件. 考虑整体问题的非线性和非凸性质，提出半整体方法提高了第二步初始解的质量. 选取数值算例进行优化分析，使用多轴3D打印系统对经典算例进行打印验证，证明所提方法的有效性. 1)考虑悬垂约束时，与基于三轴打印的方法相比，基于多轴打印的方法的灵活性更强. 2)分区尺寸的选择会直接影响优化结果. 当分区边界与节点紧密对齐时，可以通过迭代过程实现优化结果的改进，增强整体优化效果. 3)在打印验证过程中，局部构打印方向的变化导致轻微的表面缺陷，但悬垂效应被成功克服. 所提方法可以生成自支撑结构，带来的材料体积增加可忽略不计，证明了多轴3D打印桁架结构的设计和制造方面的优势及可行性.当前研究在突变曲率区域表面质量控制及多材料结构扩展方面仍存在局限. 后续将重点开展打印路径光顺优化算法开发，并探索工艺参数动态适配机制，以拓展该方法在复杂曲面构件中的工程应用.