异质材料模型是指在内部由不同的材料组分或拓扑结构组成的实体模型，典型例子包括多材料模型、梯度渐变材料模型和多孔结构等^[1]. 其中，多材料模型具有离散的材料分布类型，通常用有限的整数值来标识各均质材料组分；梯度渐变材料模型具有在空间上逐渐变化的材料分布特性，需要定义相关的材料分布函数来拟合材料组分在欧氏空间的分布信息.

计算机辅助设计的传统建模方法侧重于表示模型的几何形状和拓扑结构，忽略了密度变化、材料分布和复杂拓扑结构等模型内部信息，限制了异质材料模型在制造和工程中的应用，因此适用于设计、分析和制造异质材料模型的系统性建模方案对于充分发挥异质材料模型的卓越性能至关重要. 与传统的均质材料模型相比，表征异质材料模型的重要属性包括几何结构、材料分布、微观结构等.

增材制造技术通过将三维模型离散为一簇二维分层切片，利用不同的制造原理、制造设备及材料性能逐层实施加工制造. 对于任意复杂的几何结构和内部材料分布，选取适当的制造设备和制备工艺并完成相应的基本工艺规划步骤：选择制造方向、离散分层切片、逐层路径规划，可得满足目标需求的制造结果. 现有的切片方法研究主要集中于从不同的表面模型中分层得到精确的切片轮廓，不包含材料分布、密度变化等内部信息.

若提升相关制造工艺来实现多种材料组分在模型内部不同位置的连续变化，则增材制造技术能够通过定制材料性能制造出异质材料物件. 适用于制造异质材料的增材制造技术有很多种，这些技术的不同之处在于所使用的材料性能、机器结构及制造原理^[2]. 立体光刻技术通过紫外线光逐层固化均匀分布的光聚合物液体树脂，沿着制造方向改变工艺参数，能够固化得到材料组分呈现梯度变化的异质材料物件. 熔融沉积成型技术的制造原理是将热塑性塑料细丝通过加热喷嘴熔融并挤出，当加热喷嘴允许不同热塑性材料的任意混合时，装配多个喷嘴的FDM机器可以制造出具有梯度材料特性的物体. 基于激光的增材制造技术可以根据使用的机器类型和材料性能来调整工艺参数，制造得到具有相当结构强度和制造精度的异质金属物件. 与传统的制造技术相比，这些适用于制造异质材料物件的增材制造技术均提供了灵活的工艺原理，设计并制造复杂的异质材料分布类型.

本文着重于几何结构和材料分布这2种基本属性，提出基于三维T样条的异质材料实体重建框架. 基于三维T型控制网格的自适应结构，重建过程通过最小化相应的能量泛函来逐步拟合四面体网格的几何结构及材料属性分布，达到降低拟合误差的目的.

材料属性在实体模型的传统表达中添加了额外的维度，使得高效地表达实体内部的材料分布信息具有一定的难度和复杂性. 传统的CAD建模方法侧重于表示模型的几何形状和拓扑结构，没有直接定义材料信息.

本质上，异质材料建模过程是将材料属性在几何结构内的分布信息表示为从欧式空间到材料空间的材料分布函数. 大量的文献研究了传统几何表示法在异质材料实体上的扩展. Kumar等^[3]使用与几何相关的材料分布定义来建模异质材料模型，将经典的几何模型表达方法直接用于材料属性表达. 网格化模型^[4-5]，如四面体网格和六面体网格，已经被用来生成具有材料属性的高质量体网格模型. 基于网格化的模型表达缺乏高阶连续性. 为了更精确地拟合材料属性的非规则变化，需要生成大量的网格结构，增加了模型的复杂度，产生了一定的数据冗余.

随着等几何分析技术^[6]的提出和发展，具有指定材料属性的三变量样条模型实现了对材料分布的高阶连续性表达，可以通过等几何分析技术直接应用于分析过程. Hua等^[7-8]定义三变量单纯形样条模型来表示具有任意拓扑结构的复杂异质材料模型，但是计算这种三变量单纯形样条的混合函数及其导数不简单高效. 三维B样条模型和三维NURBS模型被广泛地应用于构建具有高阶连续性的异质材料模型. Zhang等^[9]提出基于三维NURBS模型的建模方法，用于构建患者的特定血管模型. Sasaki等^[10]介绍基于三维B样条的异质材料实体建模方法，在后续应用中对其直接切片，用于增材制造. B样条模型和NURBS模型都定义在一个规则的张量积域上，这使得它们在拓扑结构方面受到极大的限制. 当细化一个特定的区域时，需要在其邻近区域引入许多额外的自由度来维持张量积拓扑结构. 为了克服这一限制，T样条^[11-12]被提出并定义为NURBS模型的一种特殊结构，通过在参数域中引入T型网格结构实现了局部细化. 三维T样条作为T样条的衍生模型，它可以对具有复杂几何结构的实体实施基于局部细化的直接建模. 到目前为止，基于三维T样条的几何实体建模研究包含：基于体参数化方法来构造复杂三维T样条的方法^[13]；基于广义多立方体结构来参数化构造三维T样条的方法^[14-15]；针对亏格为零的几何实体提出有理三维T样条实体构造方法^[16]，将其推广为适用于任意亏格的通用表示方法^[17]. 尽管这些研究对异质材料建模过程的讨论较少，但为利用三维T样条更高效地表示实体内部材料属性分布提供了必要的理论基础.

在异质材料模型的建模过程中，当材料分布的定义与几何结构无关时，被归类为与几何无关的材料定义，可以依据材料分布函数直接计算材料组分信息. 对于这类异质材料模型，几何结构和材料分布之间的独立无关性可以允许定义高度复杂的几何结构或材料分布，几何结构的表示不会指导或限制材料分布的表示，反之亦然. Martin等^[18]通过使用不同的三变量样条，独立地表示实体模型的几何结构和材料分布. 在该框架下，几何模型和材料模型对应的2个控制网格共享相同的参数域，但是它们的控制点属性、基函数阶次和节点矢量等信息都是独立的. 利用该框架，通过使用共享相同控制网格的2类三维T样条，即几何三维T样条 $G({{x}})$和材料三维T样条 $M({{x}})$，分别表示异质材料实体模型的几何结构和材料分布，不仅可以保证重建模型具有精确的几何结构，而且能够保持材料分布的高阶连续性. 几何三维T样条 $G({{x}})$和材料三维T样条 $M({{x}})$分别定义如下：

(1) $\begin{split} G({{x}}) = & G(u,v,w) = \sum\limits_{i = 0}^{{{{N}}_G}} {{{{P}}_i}{\phi _i}(u,v,w)} = \\ & \sum\limits_{i = 0}^{{{{N}}_G}} {{{{P}}_i}{{N}}_i^u(u){{N}}_i^v(v){{N}}_i^w(w)} , \end{split} $

(2) $\begin{split} M({{x}}) = & M(u,v,w) = \sum\limits_{i = 0}^{{{{N}}_{{M}}}} {{{{M}}_i}{\phi _i}(u,v,w)} = \\ & \sum\limits_{i = 0}^{{{{N}}_{{M}}}} {{{{M}}_i}{{N}}_i^u(u){{N}}_i^v(v){{N}}_i^w(w)} . \end{split} $

式中： ${{P}}_{i}$为 $ G\left({{x}}\right) $中与基函数 $ {\varphi }_{i}(u,v,w) $对应的控制点，表示为欧氏空间的三维点坐标； ${{M}}_{i}$为 $ M\left({{x}}\right) $中与基函数 $ {\varphi }_{i}(u,v,w) $对应的控制点，表示为反映材料属性分布的连续单值函数； $ {N}_{G} $和 $ {N}_{M} $为相应的控制点数目. 在参数域内， ${{x}}=\left[u,v,w\right]$， ${{N}}_{i}^{u}$、 ${{N}}_{i}^{v}$和 ${{N}}_{i}^{w}$为基于控制点的局部节点矢量定义的单变量B样条基函数.

为了重建得到异质材料实体的几何结构 $ G\left({{x}}\right) $和材料分布 $ M\left({{x}}\right) $，在图1（a）的四面体网格与图1（d）的单元立方体参数域之间建立体参数化映射. 将四面体网格模型的边界三角形网格由图1（a）的12条灰色线段分割成6个边界网格区域，分别对应于参数立方体的6个外表面；如图1（d）所示，利用体参数化映射方法，将这12条线映射到参数立方体的对应边上，并将四面体网格的内部网格点映射到参数立方体内部，得到相应的参数坐标 $ ({u}_{i},{v}_{i},{w}_{i}) $.

在体参数化的过程中，在四面体网格的边界三角网格通过求解伞算子来定义表面参数化映射：

(3) $ \Delta {f}_{i}={\displaystyle \sum\limits _{j\in {N}_{i}}{\omega }_{ij}({f}_{j}-{f}_{i})}=0.$

式中：函数值 $ {f}_{i} $为顶点 $ {v}_{i} $处相应的参数值； $ {N}_{i} $为边界三角网格中顶点 $ {v}_{i} $的所有邻接顶点的集合； $ {\omega }_{ij} $为分配给边 $ {v}_{i}{v}_{j} $的权重，

(4) ${\omega _{ij}} = \left[\tan \;\left({\theta }/{2}\right) + \tan\; \left({\beta }/{2}\right)\right]/\left\| {{{{v}}_j} - \left. {{{{v}}_i}} \right\|} \right.,$

其中 $ \theta $和 $ \beta $为三角网格中边 $ {v}_{i}{v}_{j} $对应的二面角.

通过将如图1（a）所示的 $ 12 $条灰色线段处的参数值作为约束条件，求解伞算子可以得到边界三角网格中每一顶点的参数值. 将参数化的边界三角网格作为约束条件，求解线性方程组^[19]可以计算出四面体网格的每个内部顶点的参数值，完成体参数化映射并得到如图1（d）所示的立方体参数域. 线性方程组中的每个线性方程式对应于一个内部顶点，根据四面体网格的几何结构可以定义为

(5) $\left( {\sum\limits_{k \in {T_k}} {\frac{{S_{k,i}^2}}{{9{V_k}}}} } \right){f_i} + \sum\limits_{j \in {N_i}} {\left( {\sum\limits_{l \in {T_{ij}}} {\frac{{{S_{l,i}}{S_{l,j}}}}{{9{V_k}}}\cos\; {\theta _{l,ij}}} } \right)} {f_j} = 0.$

在四面体网格中， $ {T}_{k} $为顶点 $ {v}_{i} $的所有邻接四面体的集合， $ {N}_{i} $为顶点 $ {v}_{i} $的所有邻接顶点的集合， $ {T}_{ij} $为边 $ {v}_{i}{v}_{j} $的所有邻接四面体的集合， $ {V}_{k} $为顶点 $ {v}_{i} $邻接的第 $ k $个四面体的体积， $ {S}_{k,i} $为其中与顶点 $ {v}_{i} $相对的三角形的面积， $ {V}_{l} $为与边 $ {v}_{i}{v}_{j} $邻接的第 $ l $个四面体的体积， $ {S}_{l,i} $和 $ {S}_{l,j}$分别为与顶点 $ {v}_{i} $和 $ {v}_{j} $相对的三角形 $ {F}_{l,i} $和 $ {F}_{l,j} $的面积， $ {\theta }_{l,ij} $为 $ {F}_{l,i} $和 $ {F}_{l,j} $构成的二面角.

基于四面体网格的异质材料实体模型表示为一簇四面体网格单元的集合，每个四面体网格单元分别表示为一簇顶点的集合，每个顶点内可以同时存储相应的几何坐标和材料组分信息. 为了精确表示异质材料实体，四面体网格模型需要在几何特征和材料特征处生成大量的细分网格结构. 基于几何特征和材料特征得到的网格单元自适应分布可以指导三维T型控制网格的自适应初始化过程.

在三维T型控制网格的自适应初始化过程中，以参数域包围的单位立方体网格为初始控制网格实施均匀细分. 如图2所示，在均匀细分的控制网格内部实施基于八叉树结构的自适应细分过程，直至细分得到的每个子网格单元内部至多包含一定数目的四面体网格单元，得到可以精确表示几何特征和材料特征分布情况的三维T型控制网格.

在三维T型控制网格中计算初始几何三维T样条 $ {G}^{0}\left({{x}}\right) $和初始材料三维T样条 $ {M}^{0}\left({{x}}\right) $的节点矢量和基函数，初始化 $ {G}^{0}\left({{x}}\right) $和 $ {M}^{0}\left({{x}}\right) $的控制点分别为 $ {{P}}_{i}^{0}={I}_{G}\left({{{x}}}_{i}^{0}\right) $和 $ {{M}}_{i}^{0}={I}_{M}\left({{{x}}}_{i}^{0}\right) $. 其中， ${I}_{G}\left({{x}}\right):{\bf{R}}^{3}\to$ ${\bf{R}}^{3}$和 ${I}_{M}\left({{x}}\right):{\bf{R}}^{3}\to {\bf{R}}$分别表示从参数域 ${{x}}=\left[u,v,w\right]$到输入四面体网格及其材料分布之间的逆参数化映射. $ {{{x}}}_{i}^{0} $为在三维T型控制网格中使用Greville Abscissae定义得到的参数值，满足参数域内的恒等映射 ${\displaystyle\sum }_{i=0}^{N}{{{x}}}_{i}^{0}{\varphi }_{i}\left({{{x}}}_{i}^{0}\right)={{{x}}}_{i}^{0}$. 利用初始化控制点 ${{P}}_{i}^{0}= {I}_{G}\left({{{x}}}_{i}^{0}\right)$和 $ {{M}}_{i}^{0}={I}_{M}\left({{{x}}}_{i}^{0}\right) $得到的 $ {G}^{0}\left({{x}}\right) $和 $ {M}^{0}\left({{x}}\right) $具有较高的初始拟合精度，大大提升了逐步迭代拟合过程的收敛速度.

最小化能量泛函的方法被广泛应用在图像配准领域. 与异质材料的实体重建过程不同，图像配准通过空间变换来构建2幅图像之间的最佳对应映射. 其中，二维图像配准^[20-21]和三维图像配准^[22-23]等非刚性图像配准方法均定义了相应的配准能量泛函来用于捕获大变形或拓扑结构变化. 本文将能量泛函的定义扩展到异质材料实体的重建过程，利用能量泛函来度量输入的四面体网格和几何三维T样条 $ G\left({{x}}\right) $或者材料三维T样条 $ M\left({{x}}\right) $之间对应于几何结构或者材料分布信息的拟合误差. 通过最小化能量泛函来进行控制点的逐步迭代优化，以达到降低拟合误差的目标. 基于三维T样条的局部控制性质，求解过程中不需要构造大的矩阵，在保证拟合精度的前提下，提高了计算效率.

在几何建模过程中使用的能量泛函 $ E\left(G\left({{x}},t\right)\right) $和在材料建模过程中使用的能量泛函 $ E\left(M\left({{x}},t\right)\right) $分别定义如下。

(6) $ \begin{split} E\left( {G({{x}},t)} \right) =& {A_G}({{x}},t) + {\lambda _1}\displaystyle\int_\varOmega{B_G}({{x}},t){\rm{d}}\varOmega +\\ &{\lambda _2} \displaystyle\int_\varOmega {C_G}({{x}},t) {\rm{d}}\varOmega + {\lambda _3}{D_G}({{x}},t). \end{split} $

式中：

$ {\begin{array}{l} \!\!\!\!\!\!\!{A_G}({{x}},t) = \displaystyle\int_\varOmega {{g_G}({{x}},t)\left\| {{I_G}({{x}}) - G({{x}},t)} \right\|} _2^2{\rm{d}}\varOmega ,\\ \!\!\!\!\!\!\!{B_G}({{x}},t) = { {\left\| {{G_{,u}}({{x}},t)} \right\|_2^2 + \left\| {{G_{,v}}({{x}},t)} \right\|_2^2 + \left\| {{G_{,w}}({{x}},t)} \right\|_2^2} }, \\ \!\!\!\!\!\!\!{C_G}({{x}},t)\!\! =\! {\left\| \!{{G_{,u}}({{x}},\!t\!)} \right\|_2^2\left\| {{G_{,v}}({{x}},t)} \right\|_2^2 - {\left\langle\! {{G_{,u}}({{x}},t),{G_{,v}}({{x}},t)} \right\rangle }^2} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\; {\left\| {{G_{,v}}({{x}},t)} \right\|_2^2\left\| {{G_{,w}}({{x}},t)} \right\|_2^2 - {\left\langle {{G_{,v}}({{x}},t),{G_{,w}}({{x}},t)} \right\rangle }^2}\! + \!\!\\ \;\;\;\;\;\;\;\;\;\;\;\;\; {\left\| {{G_{,w}}({{x}},t)} \right\|_2^2\left\| {{G_{,u}}({{x}},t)} \right\|_2^2 - {\left\langle {{G_{,w}}({{x}},t),{G_{,u}}({{x}},t)} \right\rangle } ^2} ,\\ \!\!\!\!\!\!\!{D_G}({{x}},t) = \displaystyle\int_{\partial {\varOmega _u}} {\dfrac{{\left\| {{G_{,v}}({{x}},t)} \right\|_2^2\left\| {{G_{,w}}({{x}},t)} \right\|_2^2}}{{{{ {\left\langle {{G_{,v}}({{x}},t) \times {G_{,w}}({{x}},t),{N_{{\rm{tet}}}}({{x}})} \right\rangle } }^2}}}{\rm{d}}\partial {\varOmega _u}} + \\ \;\;\;\;\;\;\;\;\;\displaystyle\int_{\partial {\varOmega _v}} {\dfrac{{\left\| {{G_{,w}}({{x}},t)} \right\|_2^2\left\| {{G_{,u}}({{x}},t)} \right\|_2^2}}{{{{ {\left\langle {{G_{,w}}({{x}},t) \times {G_{,u}}({{x}},t),{N_{{\rm{tet}}}}({{x}})} \right\rangle } }^2}}}{\rm{d}}\partial {\varOmega _v}} + \\ \;\;\;\;\;\;\;\;\;\displaystyle\int_{\partial {\varOmega _w}} {\frac{{\left\| {{G_{,u}}({{x}},t)} \right\|_2^2\left\| {{G_{,v}}({{x}},t)} \right\|_2^2}}{{{{\left\langle {{G_{,u}}({{x}},t) \times {G_{,v}}({{x}},t),{N_{{\rm{tet}}}}({{x}})} \right\rangle }^2}}}{\rm{d}}\partial {\varOmega _w}} . \end{array} }$

(7) $\begin{split} E\left( {M({{x}},t)} \right) =& {A_M}({{x}},t) + {\lambda _1} \displaystyle\int_\varOmega {B_M}({{x}},t){\rm{d}}\varOmega +\\ &{\lambda _2} \displaystyle\int_\varOmega {C_M}({{x}},t){\rm{d}}\varOmega.\\ \end{split} $

式中：

$ {\begin{array}{l} \!\!\!\!{A_M}({{x}},t) = \displaystyle\int_\varOmega {{g_M}({{x}},t)\left\| {{I_M}({{x}}) - M({{x}},t)} \right\|} _2^2{\rm{d}}\varOmega, \\ \!\!\!\!{B_M}({{x}},t) = {\left\| {{M_{,u}}({{x}},t)} \right\|_2^2 + \left\| {{M_{,v}}({{x}},t)} \right\|_2^2 + \left\| {{M_{,w}}({{x}},t)} \right\|_2^2} ,\\ \!\!\!\!{C_M}({{x}},t) = { {\left\| {{M_{,u}}({{x}},\!t)} \right\|_2^2\left\| {{M_{,v}}({{x}},t)} \right\|_2^2 \!-\! {\left\langle {{M_{,\!u}}({{x}},t),{M_{,v}}({{x}},t)} \right\rangle }^2} \!+} \\ \;\;\;\;\;\; {\left\| {{M_{,v}}({{x}},t)} \right\|_2^2\left\| {{M_{,w}}({{x}},t)} \right\|_2^2 - {\left\langle {{M_{,v}}({{x}},t),{M_{,w}}({{x}},t)} \right\rangle }^2} +\\ \;\;\;\;\;\; { { {\left\| {{M_{,w}}({{x}},t)} \right\|_2^2\left\| {{M_{,u}}({{x}},t)} \right\|_2^2 - {\left\langle {{M_{,w}}({{x}},t),{M_{,u}}({{x}},t)} \right\rangle }^2}} }. \end{array} } $

式中： $ {\lambda }_{1} $、 $ {\lambda }_{2} $、 $ {\lambda }_{3} $为预先给定的约束系数， $ ，G\left({{x}},t\right) $和 $ M\left({{x}},t\right) $分别表示随时间 $ t $动态更新的几何三维T样条和材料三维T样条， $ ，{G}_{,u}\left({{x}},t\right) $、 $ {G}_{,v}\left({{x}},t\right) $、 $ {G}_{,w}\left({{x}},t\right) $及 $ {M}_{,u}\left({{x}},t\right) $、 $ {M}_{,v}\left({{x}},t\right) $、 $ {M}_{,w}\left({{x}},t\right) $分别为 $ G\left({{x}},t\right) $和 $M \left({{x}},t\right)$关于参数 $ u $、 $ v $、 $ w $的一阶偏导，<·, ·>为内积算子.

速度约束项 $ g\left({{x}},t\right) $可以保证逐步迭代拟合过程在 $ G\left({{x}},t\right) $或 $ M\left({{x}},t\right) $的梯度较大的区域内延缓最小化过程，同时在梯度较小的区域内加速最小化过程，具体定义为

(8) $ {g_G}\left( \!{{{x}},t} \right)\! =\!\left(\! { {{\sqrt {\gamma + \left\| {{G_{,\!u}}\left( {{{x}},\!t} \right)} \right\|_2^2 \!+\! \left\| {{G_{,v}}\left( {{{x}},\!t} \right)} \right\|_2^2\! +\! \left\| {{G_{,w}}\left( {{{x}},t} \right)} \right\|_2^2} }}} \right)^{-1}, $

(9) $ {g_M}\left(\! {{{x}},\!t} \right) \!=\left( {\! {{\sqrt {\gamma + \left\| {{M_{,\!u}}\left( \!{{{x}},t} \right)} \!\right\|_2^2\! + \!\left\| {{M_{,v}}\left( \!{{{x}},\!t} \right)} \right\|_2^2\! +\! \left\| {{M_{,w}}\left( {{{x}},t} \!\right)} \!\right\|_2^2} }}} \!\right)^{-1}. $

式中： $ \gamma $为预先给定的非零常值，用于避免除零的情况出现，通常设置为 $ 1{0}^{-2} $.

在以上2个能量泛函 $ E\left(G\left({{x}},t\right)\right) $和 $ E\left(M\left({{x}},t\right)\right) $中，第1项是逐步迭代拟合过程中动态更新的三维T样条 $ G\left({{x}},t\right) $或者 $ M\left({{x}},t\right) $与输入四面体网格模型之间的拟合误差的平方和. 与约束系数 $ {\lambda }_{1} $和 $ {\lambda }_{2} $分别对应的第2项和第3项是在每个能量泛函中引入的正则约束项，通过实施这些约束，以确保控制网格在逐步迭代优化过程中实现平滑变形，防止发生歧义变形. 其中，第2项为一阶正则约束项，保证了 $ G\left({{x}},t\right) $或 $ M\left({{x}},t\right) $在 $ u $、 $ v $、 $ w $方向上的平稳变化；第3项为2阶正则约束项，保证了优化过程中控制网格单元中每个面面积的平稳变化. 在几何建模过程的能量泛函 $ E\left(G\left({{x}},t\right)\right) $中，与约束系数 $ {\lambda }_{3} $对应的第4项是在控制网格边界引入的额外正则约束项，测量四面体网格模型与 $ G\left({{x}},t\right) $之间的边界法向量差，保证了优化过程中几何特征的拟合精度. 约束系数 $ {\lambda }_{1} $、 $ {\lambda }_{2} $和 $ {\lambda }_{3} $均为经验参数，数值对具体的优化过程不敏感，因此在优化过程开始实施时，通常将它们设置为常量，例如 $ 1{0}^{-3} $.

如图3（a）~（c）所示，几何模型的重建过程是将 $ G\left({{x}}\right) $拟合到输入四面体网格模型并保证拟合精度. 如图3（d）~（f）所示，材料模型的重建过程是将欧式空间中的每个点映射至材料空间，得到材料组分，以重建模型内部连续变化的材料组分. 以材料三维T样条 $ M\left({{x}}\right) $作为材料分布函数来表示异质材料分布，得到从参数域到材料空间的材料分布函数.

在最小化能量泛函的过程中，计算能量泛函对控制点的微分结果，可得

(10) $\frac{{{\rm{d}}{{{P}}_i}(t)}}{{{\rm{d}}t}} = - {\rm{\delta}} {E_i}(G({{x}},t)),$

(11) $\frac{{{\rm{d}}{{{M}}_i}(t)}}{{{\rm{d}}t}} = - {\rm{\delta}} {E_i}(M({{x}},t)).$

引入时间步长 $ \varepsilon $，通过欧拉法对方程的左边进行近似，将其转化为常微分方程，可得

(12) $\frac{{{{P}}_i^{s + 1} - {{P}}_i^s}}{\varepsilon } = - {\text{δ}} {E_i}({G^s}({{x}},t)),$

(13) $\frac{{{{M}}_i^{s + 1} - {{M}}_i^s}}{\varepsilon } = - {\text{δ}} {E_i}({M^s}({{x}},t)).$

式中： $ \varepsilon $为逐步迭代拟合过程中的时间步长.

(14) $ \begin{split} &{{{\delta E}}_{{{i}}}}\left( {{{{G}}^{{s}}}\left( {{{{ x}, t}}} \right)} \right) = - 2{R_G}\left( {{{{ x}, t}}} \right) + {\lambda _1}\int_{ {\varOmega }} {{S_G}} \left( {{{{ x}, t}}} \right){\rm{d}}\varOmega + \\ &\;\;\;\;\;\;\;\;\;\;\;\;\;{\lambda _2}\int_\varOmega {{T_G}} \left( {{{{ x}, t}}} \right){\rm{d}}\varOmega + {\lambda _3} \left[ { \int_{\partial {\varOmega _u}}{{U_G}} \left( {{{{ x}, t}}} \right){\rm{d}}\partial {\varOmega _u} + } \right.\\ &\;\;\;\;\;\;\;\;\;\;\;\;\left. {\int_{\partial {\varOmega _v}}{{V_G}} \left( {{{{ x}, t}}} \right){\rm{d}}\partial {\varOmega _v} + \int_{\partial {\varOmega _w}}{{W_G}} \left( {{{{ x}, t}}} \right){\rm{d}}\partial {\varOmega _w}} \right].\\ \end{split} $

式中：

$ {\begin{split} &\!\!\!R_G( x, t) = \displaystyle\int_\varOmega {{g_G}({{x}},t)\left( {{I_G}({{x}}) - {G^s}({{x}},t)} \right)} {\rm{d}}\varOmega , \\ &\!\!\!S_G( x, t)= 2G^s_{,u}({{x}},t){\phi _{i,u}}({{x}}) +\\ & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;2G^s_{,v}({{x}},t){\phi _{i,v}}({{x}}) + 2G^s_{,w}({{x}},t){\phi _{i,w}}({{x}}) , \\&\!\!\! T_G( x, t)= \\ & {2\left\|\! {G^s_{,u}({{x}},t) \!\times G^s_{,v}({{x}},t)} \!\right\|\left(\! {{\phi _{i,u}}({{x}})G^s_{,v}\!({{x}},t) \!+\! {\phi _{i,v}}({{x}})G^s_{,u}({{x}},t)} \right)} \! + \!\!\\ & {2\left\| {G^s_{,v}({{x}},t) \times G^s_{,w}({{x}},t)}\! \right\|\left(\! {{\phi _{i,v}}({{x}})G^s_{,w}({{x}},t)\! + \!{\phi _{i,w}}({{x}})G^s_{,v}({{x}},t)} \right)} \! +\! \\& {{2\left\| {G^s_{,w}({{x}},t)\! \times G^s_{,u}\!({{x}},t)} \right\|\left(\! {{\phi _{i,w}}\!({{x}})G^s_{,u}({{x}},t)\! + \!{\phi _{i,u}}({{x}})G^s_{,w}({{x}},t)}\! \right)} } , \\ &\!\!\!U_G({{x}}, t)= - 2 \left[ {{{\left\| {G^s_{,v}({{x}},t)} \right\|_2^2\left\| {G^s_{,w}({{x}},t)} \right\|_2^2}} } \right. \times \\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {{\phi _{i,v}}({{x}})G^s_{,w}({{x}},t) + {\phi _{i,w}}({{x}})G^s_{,v}({{x}},t)} \right) \Bigg / \\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{ {\left\langle {G^s_{,v}({{x}},t) \times G^s_{,w}({{x}},t),{N_{{\rm{tet}}}}({{x}})} \right\rangle } }^3 {}} \right] , \\ \end{split}} $

$ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\begin{split} \!V_G({{x}}, t)= & - 2 \left[ {{{\left\| {G^s_{,w}({{x}},t)} \right\|_2^2\left\| {G^s_{,u}({{x}},t)} \right\|_2^2}} } \times \right.\\ &\left( {{\phi _{i,w}}({{x}})G^s_{,u}({{x}},t) + {\phi _{i,u}}({{x}})G^s_{,w}({{x}},t)} \right) \Bigg / \\ & \left. {{{{{ {\left\langle {G^s_{,w}({{x}},t) \times G^s_{,u}({{x}},t),{N_{{\rm{tet}}}}({{x}})} \right\rangle } }^3}}}} \right], \\ W_G({{x}}, t)= & - 2 \left[ {{{\left\| {G^s_{,u}({{x}},t)} \right\|_2^2\left\| {G^s_{,v}({{x}},t)} \right\|_2^2}}} \times \right.\\ &\left( {{\phi _{i,u}}({{x}})G^s_{,v}({{x}},t) + {\phi _{i,v}}({{x}})G^s_{,u}({{x}},t)} \right)\Bigg / \\ &\left. {{{{{ {\left\langle {G^s_{,u}({{x}},t) \times G^s_{,v}({{x}},t),{N_{{\rm{tet}}}}({{x}})} \right\rangle }}^3}}} } \right] .\end{split} $

(15) $ \begin{split} \delta {E_i}\left( {{M^s}({{x}},t)} \right) =& - 2{R_M}({{x}},t) + {\lambda _1}\int_\varOmega {{S_M}({{x}},t){\rm{d}}\varOmega + } \\ &{\lambda _2}\int_\varOmega {{T_M}({{x}},t) } {\rm{d}}\varOmega .\\ \end{split} $

式中：

$ {\begin{array}{l} \!\!\!R_M( x, t)= \int_\varOmega {{g_M}({{x}},t)\left( {{I_M}({{x}}) - {M^s}({{x}},t)} \right)} {\rm{d}}\varOmega , \\ \!\!\!S_M( x, t)= 2M^s_{,u}({{x}},t){\phi _{i,u}}({{x}}) + 2M^s_{,v}({{x}},t){\phi _{i,v}}({{x}}) + \\ \;\;\;\;2M^s_{,w}({{x}},t){\phi _{i,w}}({{x}}) , \\ \!\!\!\!T_M( x,\! t)= \\ \!\!\!2\left\| {M^s_{,u}({{x}},t) \times M^s_{,v}({{x}},t)} \right\| \left( {{\phi _{i,u}}({{x}})M^s_{,v}({{x}},t)\! +\! {\phi _{i,v}}({{x}})M^s_{,u}({{x}},t)} \right) \! +\! \\ \!\!\! 2\left\|\! {M^s_{,v}({{x}}\!,t) \!\times \!M^s_{,w}({{x}},t)}\! \right\| \left( \!{{\phi _{i,v}}({{x}})M^s_{,w}\!({{x}},t)\! + \!{\phi _{i,w}}\!({{x}})M^s_{,v}({{x}},t)} \!\right) \! +\! \\ \!\!\!2\left\|\! {M^s_{,\!w}\!({{x}},t)\! \times \!M^s_{,\!u}\!({{x}},t)\!} \right\|\!\times\! \left(\! {{\phi _{i,\!w}}({{x}})M^s_{,\!u}(\!{{x}},t) \!+\! {\phi _{i,u}}\!({{x}})M^s_{,\!w}({{x}},\!t)}\! \right) . \\ \end{array} } $

通过计算

(16) ${{{P}}^{s + 1}} = {{{P}}^s} - \varepsilon {\rm{\delta}} {{E}}({G^s}({{x}},t)),$

(17) ${{{M}}^{s + 1}} = {{{M}}^s} - \varepsilon \delta {{E}}({M^s}({{x}},t)),$

可以动态更新控制点，生成拟合精度更高的三维T样条. 在进一步实施迭代拟合步骤之前，计算当前的拟合相似度，与上一步的拟合相似度进行比较. 若相似度的变化低于给定容差，则迭代循环终止；反之，则进一步实施迭代拟合，直至满足给定容差.

图1（b）中的 $ G\left({{x}}\right) $以及图1（f）中的 $ M\left({{x}}\right) $是由图1（e）中的自适应三维T型网格逐步迭代重建得到的几何三维T样条和材料三维T样条. 由于 $ G\left({{x}}\right) $与 $ M\left({{x}}\right) $具有相同的参数域，如图1（c）所示，可以通过参数域的中间运算，得到材料分布在欧式空间的映射结果.

按照三维模型的表达方式，对分层切片方法的研究主要分为2类：网格模型切片和直接切片方法. 基于网格的模型易于建模和操作，在选择制造方向、添加支持结构等过程具有更好的灵活性. 网格模型缺乏高阶连续性，在网格化过程中需要添加大量的网格单元来表示复杂的几何特征；随着引入复杂的异质材料属性，数据量过大的情况会更加严重. 相反，直接切片方法直接应用于精确的模型表示，在不涉及网格化过程的情况下，不仅可以减少预处理时间，而且能够保持模型从原始输入到分层切片结果乃至最终制造产品的一致性.

随着各种制造设备和工艺的发展，适用于异质材料模型的制造方法得到了广泛的关注和研究. 为了实现异质材料实体模型在增材制造领域的应用，制造方向选择、离散分层切片及逐层路径规划等相关工艺步骤亟需研究. 本文采用适用于异质材料实体模型的直接切片算法^[24]，可以沿着给定的制造方向实施直接分层切片，得到基于三角网格表示的二维切片结果.

以股骨下端模型为例，选择Z轴为制造方向，具体说明分层切片的实施过程. 对于任意给定的Z轴坐标值 $ \widehat{z} $，在参数域包围的立方体网格中实施基于八叉树结构的自适应细分过程，得到包围 $ \widehat{z} $等值面的分层T网格；基于分层T网格实施直接切片，得到参数域中的切片结果，利用几何三维T样条 $ G\left({{x}}\right) $映射得到欧氏空间中的分层切片结果.

基于八叉树结构自适应细分生成分层T网格包含以下步骤.

1）综合考虑目标模型的几何结构及实际制造过程中的工艺参数、目标制造精度等，给定细分精度 $ \tau $.

2）在参数域包围的立方体网格中进行初始化均匀细分，得到初始分层T网格. 初始化均匀细分的阶数是一个经验参数值，需要综合考虑目标模型的几何结构及实际制造过程中的工艺参数等进行实验性设置. 细分阶数越高，算法实施过程中的数据存储量越大，但后续步骤的计算过程会越快；反之，细分阶数越低，算法实施过程的数据存储量越小，但后续步骤的计算过程会越慢. 在本例的具体实施过程中，初始化均匀细分的阶数设置为1.

3）在当前分层T网格的每个叶网格内，利用 $ G\left({{x}}\right) $计算得到该网格单元中8个网格顶点处的欧式坐标. 在每个网格顶点处，根据Z轴坐标值是否大于（或小于） $\; \widehat{z} $进行分类标记，将当前分层T网格的所有叶网格分类为边界网格(见图4的灰色网格单元)或者内部网格(见图4的透明网格单元).

4）以标记后的分层T网格为基础，利用 $ G\left({{x}}\right) $计算得到每个边界网格在欧式空间的体积 $ V $. 基于八叉树结构的自适应细分过程，将体积 $ V $大于细分精度 $ \tau $的边界网格细分为8个子网格单元.

步骤3）的叶网格分类标记过程以及步骤4）的自适应细分过程交替迭代实施，直至分层T网格的每个边界网格均小于细分精度 $ \tau $. 以股骨下端模型为例，给定Z轴坐标值 $ \widehat{z}=65 $，图4（a）~（c）所示的灰色边界网格单元为满足3种不同的细分精度 $ \tau =10{\text{、}}1 $和 $ 0.1 $要求的分层T网格.

分层T网格的边界网格内，8个网格顶点包含不同的材料属性，因而可以找到边界网格内部的变化边. 应用迭代线性插值方法，计算变化边上Z轴坐标值等于 $ \widehat{z} $的交点. 如图5（a）~（c）所示，应用移动立方体方法^[25]连接每一个边界网格内的所有交点，得到参数域中与 $ \widehat{z} $等值面对应的参数三角网格面. 如图5（d）~（f）所示，利用几何三维T样条 $ G\left({{x}}\right) $将图5（a）~（c）的参数三角网格面映射到欧式空间，可得相应的三角网格分层切片结果.

以具有异质材料信息的四面体网格为输入模型，在C++中实现了模型重建与切片流程，以股骨下端和上端模型及橡皮鸭模型为例进行数值实验.

在股骨下端和上端模型的四面体网格中，每个网格顶点的材料属性依据几何结构定义得到：在边界三角网格中，利用热运动方程计算得到每个网格顶点到模型底面的测地距离；以基于边界三角网格的测地距离场为约束条件，通过式 $ \left(4\right) $计算得到四面体网格的每个内部网格点到模型底面的测地距离，完成如图1（a）、6（a）所示的材料属性定义. 如图7（a）所示，橡皮鸭模型中每个网格顶点的材料属性定义为该网格顶点到一个指定的内部网格顶点之间的欧氏距离.

在基于最小化能量泛函的重建过程中，在每一次迭代拟合步骤结束之后，需要对当前的拟合精度采用相似度（similarity ratio, RS）进行度量。较高的RS对应较好的拟合结果，RS=100%对应完美的拟合结果. 几何三维T样条 $ G\left({{x}}\right) $和材料三维T样条 $ M\left({{x}}\right) $的相似度分别定义为

(18) ${\rm{R}}{{\rm{S}}_G} = \left( {1 - \frac{{{{\left\| {{I_G}({{x}}) - G({{x}})} \right\|}_2}}}{{{{\left\| {{I_G}({{x}}) - {G^0}({{x}})} \right\|}_2}}}} \right) \times 100{\text{%}} ,$

(19) ${\rm{R}}{{\rm{S}}_M} = \left( {1 - \frac{{{{\left\| {{I_M}({{x}}) - M({{x}})} \right\|}_2}}}{{{{\left\| {{I_M}({{x}}) - {M^0}({{x}})} \right\|}_2}}}} \right) \times 100{\text{%}} .$

式中： $ G\left({{x}}\right) $和 $ M\left({{x}}\right) $分别为迭代拟合得到的几何三维T样条和材料三维T样条， $ {I}_{G}\left({{x}}\right) $和 $ {I}_{M}\left({{x}}\right) $分别为从参数域到输入四面体网格及其材料分布的逆参数化映射， $ {G}^{0}\left({{x}}\right) $和 $ {M}^{0}\left({{x}}\right) $分别为逐步迭代拟合过程的初始三维T样条.

在三维T型控制网格的自适应初始化过程中，初始化控制网格直至自适应细分得到的每个子网格单元内部至多包含5个四面体网格单元. 在渐进式重建过程中， $ \varepsilon $是一个经验值，可以通过观察优化过程的收敛性而设定，以避免在优化过程最终收敛之前产生异常中断. 在积分方程 ${\delta E}_{i} \left({G}^{s}\left({{x}},t\right)\right)$和 ${\delta E}_{i}\left({M}^{s}\left({{x}},t\right)\right)$的求解过程中，基于高斯求积规则来计算相应的加权平均函数值.

3个实例的重建结果如表1所示。表中，n_c为控制点数目，n_i为逐步迭代拟合次数. 基于最小化能量泛函的渐进式重建框架能够通过逐步迭代拟合来重建复杂的几何特征和材料特征. 如图2（c）、6（c）、7（c）所示为三维控制网格中引入的T型节点. 这些T型节点与目标四面体网格的几何特征和材料特征之间存在着精确匹配，仅在特征复杂的区域添加，其余控制网格保持粗糙.

三维NURBS模型定义在基于均匀细分的三维控制控制网格上. 如表2所示，在异质材料实体建模过程中需要大量的冗余控制点来满足张量积拓扑约束，控制点的数量明显高于三维T样条. 尤其是当几何结构或者材料分布中存在复杂特征时，三维NURBS模型需要在拟合精度和数据量之间进行权衡；三维T样条模型能够利用局部控制性质，在某些复杂特征分布的区域中自适应地插入节点，以改善模型的自由度，从而可以在达到相似甚至更优的拟合精度的前提下减少冗余控制点的数量，使得表达方式更加紧凑.

在实施直接切片流程之前，选择Z轴作为切片方向，使得分层切片簇沿着Z轴方向均匀分布. 任意给定一簇Z轴坐标值 $ \left\{{z}_{i}\right\} $，如图8（a）、9（a）、10（a）所示，基于八叉树结构自适应细分，得到与每个 $ {z}_{i} $等值面对应的分层T网格；如图8（b）、9（b）、10（b）所示，采用移动立方体方法遍历边界网格，生成一簇参数三角网格面；如图8（c）、9（c）、10（c）所示，利用几何三维T样条 $ G\left({{x}}\right) $，将参数三角网格面簇映射到欧式空间，得到相应的分层切片结果.

基于4.2节的直接切片方法，如图11~13所示，沿着材料三维T样条 $ M\left({{x}}\right) $的材料梯度方向对异质材料属性实施等材料面提取过程. 由于 $ G\left({{x}}\right) $及 $ M\left({{x}}\right) $共享相同的参数域，在具体实施等材料面提取过程中，对于任意给定的材料值 $ m\in \left[0,M\right] $，以 $ M\left({{x}}\right) $为基础，在参数域中自适应细分，得到包围 $ M\left({{x}}\right)=m $这一等材料面的分层T网格；基于分层T网格实施移动立方体算法，得到参数域中的等材料面，利用 $ G\left({{x}}\right) $映射得到欧氏空间中的等材料面结果，如图11（b ）、12（b）、13（b）所示.

本文提出用于异质材料实体的建模和切片工艺，将几何结构和材料分布由基于相同三维T型控制网格的2类三维T样条定义. 通过最小化能量泛函来逐步迭代拟合控制点，降低构造的三维T样条与输入的四面体网格之间的拟合误差. 针对异质材料实体在增材制造中的切片工艺，采用基于八叉树结构的自适应细分过程，在参数域生成分层T网格，应用移动立方体方法生成三角网格切片结果. 数值结果表明，三维T型控制网格能够局部自适应控制几何结构和材料分布，使得异质材料实体的重建过程和切片工艺可以得到较高的拟合和切片精度.

本文所采用的体参数化方式很难捕获尖锐突起之类的复杂几何特征，使得材料建模过程无法进一步提高拟合精度. 采用基于Polycube的复杂参数域结构，实施适用于亏格非零的复杂实体模型的体参数化过程. 在直接切片过程中，制造方向是预先给定的. 虽然这种方式对整个直接切片过程没有影响，但是为了实现设计-分析的无缝集成，如何自动获得最佳制造方向是亟待解决的问题.