(北京航空航天大学机械工程及自动化学院，北京 100191）

等几何分析（isogeometry analysis）是一种基于样条理论的数值计算方法[1]。它的核心在于将CAD建模工具和有限元求解方式结合起来，使几何模型与分析网格模型表达统一，弥补了几何模型与分析模型巨大差异的缺陷[2]。基于此，等几何分析具有较高的效率和更精确的分析结果，目前已成功地应用于振动分析[3-4]、流体力学[5-6]、板壳[7-8]以及电磁[9]等领域。

非均匀有理B样条曲面（Non-Uniform Rational B-Splines，NURBS）是应用较为广泛的一种CAD设计模型。一些学者对基于NURBS模型的等几何分析进行了研究，取得了一些成果。但NURBS模型严格的矩形拓扑限制使得其在进行等几何分析时依旧存在以下不足：首先，使用NURBS曲面表达拓扑结构比较复杂的模型时，不可避免地会出现多张NURBS曲面的裁剪和拼接，增大了等几何分析中网格划分的难度；其次，NURBS曲面无法进行局部细化，因此分析者对分析模型的某个局部进行加密时，会不可避免地引入冗余的控制顶点，从而降低计算效率。

Sederberg等提出的T样条曲面（T-splines）为解决NURBS的不足提供了契机[10-11]。相比较于NURBS曲面，T样条曲面具有以下两个优点：（1）单张T样条曲面可表示拓扑结构比较复杂的模型[12]。（2）局部细化[13-15]。因此，T样条曲面是一种比较适合于进行等几何分析的模型。

为了推动基于T样条曲面的等几何分析发展，Scott等提出了T样条曲面的Bézier提取算法[16]。Nguyen将基于Bézier单元的等几何分析结果与传统的有限元分析进行了比较[17]。之后，Scott等又提出了带奇异点的T样条曲面的Bézier提取方法[18]。

T样条曲面的Bézier提取算法有以下两个关键步骤：（1）对参数域T网格进行Bézier分片，使得划分后的参数域网格所对应的空间曲面片（patch）的内部C∞连续。（2）求取T样条曲面每个混合函数所对应的Bézier提取矩阵。文献[16]以及文献[18]对步骤（2）的具体实现过程进行了详细说明，但都没有给出步骤（1）的具体实现方法。为此，本文对参数域T网格的Bézier分片进行了研究，给出了一种参数域T网格的Bézier分片算法。该算法有利于后续进一步实现T样条曲面的Bézier提取算法，最终达到基于T样条曲面等几何分析的目的。

T样条曲面的基本概念

T样条曲面可看作NURBS曲面的一般化，它允许控制顶点呈T型连接。本节将会对T样条曲面的基本概念进行简要介绍。

定义1：T样条曲面的控制网格（T-mesh）称为空间T网格，它包含了T样条曲面的所有控制顶点空间位置信息以及拓扑连接信息。

定义2：T样条曲面上的控制顶点都对应一个参数点，由于这些参数点连接而成的参数域网格与空间T网格同胚，因此称该参数域网格为参数域T网格。

空间T网格与参数域T网格分别如图1（a）和（b）所示，第13号控制顶点对应于参数域上的a13点。

T样条曲面的计算表达式为：

其中Bi(u,v)=Ni,k(u)Ni,p(v)为控制顶点pi所对应的混合函数；ωi为控制顶点pi的权值；Ni，k（u），Ni，p（v）分别为k次及p次B样条基函数。T样条曲面的具体计算过程见文献[10]。

T样条曲面的Bézier单元

T样条曲面在计算过程中以参数域T网格上的面片作为最小单元，但参数域T网格上的每个面片所对应的T样条曲面片（patch）不一定能用一张有理Bézier曲面表示，这里以一个具体实例对该现象进行说明。

图2为一个双三次T样条曲面的参数域T网格。v1～v6为参数域T网格中的6个参数点，它们所对应的控制顶点分别为p1～p6，s1为参数域T网格中的一个面片，它包含v2～v65个参数点，阴影区域为控制顶点p1所对应的混合函数N1，3（u）N1，3（v）的非零区域。由B样条基函数的可微性可知，面片s1所对应的T样条曲面片内部不能达到C∞连续，而有理Bézier曲面的内部可以达到C∞连续。因此，面片s1所对应的T样条曲面片不能用一张有理Bézier曲面表示。

因此，为了实现T样条曲面的Bézier提取，必须首先将T样条曲面进行分割，使得分割后的每个T样条曲面片的内部无限次可微，而边界处则只能达到有限次可微，称具备上述性质的T样条曲面片为T样条曲面的一个 Bézier单元。

定义3：T样条曲面划分Bézier单元的过程称为T样条曲面的Bézier分割。

T样条曲面的Bézier分割并不是直接在T样条曲面上进行切分得到，而是通过对参数域T网格进行分割而间接得到的。由B样条基函数的连续性可知，T样条曲面会在混合函数的节点线所对应的空间曲线处连续性降低。因此，只需要对参数域T网格按照T样条混合函数的节点线进行分割，则可以间接实现T样条曲面的 Bézier分割。

图1 T样条曲面的空间T网格和参数域T网格Fig.1 T-mesh in Cartesian space and parametric space

参数域T网格Bézier分割算法

定义4：T样条曲面在参数域T网格上间接实现T样条曲面Bézier分割的过程称为参数域T网格的Bézier分割。

图3为一个双三次T样条曲面参数域网格，v1～v5为参数域T网格中的5个参数点，它们所对应的控制顶点分别为p1～p5，s1～s7为参数域 T网格上的7个面片。其中控制顶点p2所对应的混合函数B2（u,v）的节点线分布如图4所示。因此，参数域上的面片至少需要切分成如图5所示的参数域网格。

图2 参数域T网格1Fig.2 T-mesh one in parametric space

图3 参数域T网格2Fig.3 T-mesh two in parametric space

图4 混合函数B2（u,v）的节点线Fig.4 Knot line of blending function B2（u,v）

图5 混合函数B2（u,v）的节点线分割效果图Fig.5 T-mesh segmented by the knot line of blending function B2（u,v）

显然，如果采用蛮力算法求取每一个控制顶点所对应的混合函数的节点线分布，最终虽能够实现参数域T网格的Bézier分割，但会出现很多重复的计算，并且在切分过程中可能出现如图5右上角所示的L型点。

为减少重复计算，并避免L型分割点的产生。本文给出了包含如下4个步骤的参数域T网格Bézier分割算法。

算法1：参数域T网格Bézier分割算法。

STEP1：将参数域T网格中所有的T型参数点进行延伸，通过T型参数点的延伸将参数域T网格中面片进行切分。

STEP2：遍历参数域上的所有参数点，对于每一个参数点，执行STEP3 ～ STEP4。

STEP3：计算当前参数点上所耦合的混合函数的非零区域。

STEP4：如果STEP3中得到的非零区域不能完全覆盖参数域T网格中的某个面片，则切分该面片，否则不切分该面片。

下面以如图3所示的参数域T网格为例对参数域T网格Bézier分割算法进行解释。

首先处理T型参数点，由于图3所示的参数域T网格中参数点v2～v5均为T型连接点，因此应首先对这些T型参数点进行延伸。如图6所示，参数点v2以及v3均向v正方向延伸两个T网格距离，而参数点v4以及v5均向u负方向延伸两个T网格。经过这个延伸过程，面片s3切分成子面片f1和f2；面片s4切分成子面片f3和f4；而面片s6切分成子面片f5、f6、f7以及f8。

其次计算参数点所耦合的混合函数非零区域，如图7所示，v1为参数域T网格中的十字型参数点，参数点v1所耦合的混合函数非零区域为见图7的阴影部分。

由于参数点v1耦合的混合函数的非零区域完全覆盖了面片f1和f3，因此这两个面片不需要切分。而面片s5不能被参数点v1所耦合的混合函数非零区域完全覆盖，因此需要将图3中的面片s5切分为如图8所示的两个子面片f9和f10。

图6 T型参数点延伸后的参数域T网格Fig.6 T-mesh generated by extending T-vertex

图7 参数点v1所耦合的混合函数的非零区域Fig.7 Nonzero region of the blending function corresponding to v1

经过上述4个步骤，参数域T网格中的每个面片所对应的T样条曲面均可以由一张有理Bézier曲面表示，也就实现了参数域T网格Bézier分割过程。

从上述实例中可以看到，参数域T网格Bézier分割算法可能将一个面片分成两个面片，如图3中的s3；也可能分成四个面片，如图3中的s6。事实上，如果参数域T网格比较复杂，还会出现一分为三、一分为六等多种情形。为了实现面片切分算法，本文使用如下所示的BézierFace结构体来处理面片切分问题。这个结构体中主要包含了每一个BézierFace范围信息，以及该BézierFace所属的面片的索引号信息。

};

基于上述结构，可以使用算法2实现参数域T网格中面片的具体切分过程。

算法2：参数域T网格Bézier分割中的面片切分算法。

图8 参数域Bézier分割网格Fig.8 Bézier segmentation of the T-mesh

STEP1：建立一个队列数组，该数组的大小与T样条曲面的参数域T网格中的面片数一致，并使用与队列索引号相同的面片信息初始化该队列。

图9 参数域T网格及其延伸网格Fig.9 T-mesh and its extending T-mesh

图10 参数域Bé zier分割网格Fig.10 Bé zier segmentation of the T-mesh

STEP2：在执行算法1的过程中，如果出现面片需要切分的情形，则执行 STEP3～STEP4。

STEP3：遍历与当前面片索引号相同的队列中的所有子面片，对于每一个子面片执行STEP4。

STEP4：如果当前子面片需要切分，则建立两个新面片并将位于队头的当前子面片执行出队操作，之后将新建两个面片分别执行入队操作；如果当前子面片不需要切分，则将位于队头的当前子面片先执行出队操作，随后执行入队操作。

下面以一个具体实例对算法2进行解释。图8为一个双三次T样条曲面参数域T网格。v1、v3以及v4为参数域T网格中的3个T型参数点，s1为参数域T网格中的一个面片，其umin、umax、vmin以及vmax分别为 1、4、3和4。参数域T网格将T型参数点进行延伸后得到如图9所示的参数域延伸T网格。结合图9（a）和（b）可知，由于T型参数点v4的延伸，导致面片s1被切分成面片e1以及e2。在具体程序实现时，首先建立两个新面片e1和e2，其中新面片 e1的umin、umax、vmin以及vmax分别为3、4、3和4；新面片e2的umin、umax、vmin以及vmax分别为 1、3、3和4。之后将队列中仅有的面片s1进行出队操作，最后对新面片e1和e2依次执行入队操作，此时队列中存在e1和e2两个面片。

图10中的阴影部分为十字型参数点v2所耦合的混合函数的非零区域，这个区域完全覆盖面片e1，但并未完全覆盖图9（b）面片e2。因此，首先对队列中的e1执行出队操作，接着对e1执行入队操作，此时面片e1由队头转到队尾，这个过程相当于建立了面片f1。之后，队列中的第二个面片e2，显然需要进行切分，因此应新建两个面片f2和f3，其中新面片f2的umin、umax、vmin以及vmax分别为 2、3、3和 4；新面片f3的umin、umax、vmin以及vmax分别为 1、2、3 和 4，之后对面片e2进行出队操作，最后对面片f2和f3依次执行入队操作，此时队列中存在f1、f2及f33个面片。

参数域T网格经过上述算法的分割后，每一个参数域面片均所对应的T样条曲面片均为一个Bézier单元，这就间接实现了T样条曲面的Bézier分割。

实例验证

为了验证本文所提出的参数域T网格的Bézier分割算法的有效性，在Visual Studio 2012软件开发平台上实现了T样条曲面的Bézier提取算法。该算法首先采用本文所提出的参数域T网格的Bézier分割算法对T样条模型的参数域T网格进行处理，之后采用文献[16]提出的混合函数Bézier提取矩阵求取算法对T样条模型进行Bézier提取，部分提取结果如图11～图13所示。

这里，以face模型为例，对T样条模型以及Bézier提取模型的分割线进行说明。如图11（a）所示，face模型中的每一小曲面片均对应于参数域T网格中的一个面片，该模型经过Bézier提取算法之后，得到了如图11（b）所示的模型，该模型上的每一个面片均为T样条曲面的一个Bézier单元，cat模型与 venus模型与之类似。上述实例也反映了本文算法的可行性。

图11 face模型Fig.11 Face model

图12 cat模型Fig.12 Cat model

图13 venus模型Fig.13 Venus model

结论

本文提出了一种参数域T网格的Bézier分片算法，该算法只需要借助队列这一常见的数据结构就能间接实现T样条曲面的Bézier分割，有利于后续完整的实现T样条曲面的Bézier提取算法。在此基础上，未来可进一步研究基于T样条曲面Bézier单元的等几何分析算法。

参 考 文 献

[1] HUGHES T J R, COTTRELL J A,BAZILEVS Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[J]. Computer Methods in Applied Mechanics and Engineering, 2005,194(39):4135-4195.

[2] 吴紫俊, 黄正东, 左兵权, 等. 等几何分析研究概述[J]. 机械工程学报, 2015,51(5):114-129.WU Zijun, HUANG Zhengdong, ZUO Bingquan, et al. Perspectives on isogeometric analysis[J]. Journal of Mechanical Engineering,2015, 51(5):114-129.

[3] BAZILEVS Y, CALO V M, COTTRELL J A, et al. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows[J]. Computer Methods in Applied Mechanics and Engineering, 2007,197(1):173-201.

[4] AKKERMAN I, BAZILEVS Y, CALO V M, et al. The role of continuity in residualbased variational multiscale modeling of turbulence[J]. Computational Mechanics, 2008,41(3): 371-378.

[5] BAZILEVS Y, CALO V M, ZHANG Y, et al. Isogeometric fluid-structure interaction analysis with applications to arterial blood flow[J].Computational Mechanics, 2006,38(4/5):310-322.

[6] ZHANG Y, BAZILEVS Y, GOSWAMI S, et al. Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow[J]. Computer Methods in Applied Mechanics and Engineering, 2007,196(29):2943-2959.

[7] BENSON D J, BAZILEVS Y, HSU M C,et al. Isogeometric shell analysis: the Reissner-Mindlin shell[J]. Computer Methods in Applied Mechanics and Engineering, 2010,199(5):276-289.

[8] DORNISCH W, KLINKEL S, SIMEON B. Isogeometric Reissner-Mindlin shell analysis with exactly calculated director vectors[J].Computer Methods in Applied Mechanics and Engineering, 2013,253:491-504.

[9] BUFFA A, SANGALLI G, VÁZQUEZ R. Isogeometric analysis in electromagnetics:B-splines approximation[J]. Computer Methods in Applied Mechanics and Engineering,2010,199(17):1143-1152.

[10] SEDERBERG T W, ZHENG J,BAKENOV A, et al. T-splines and T-NURCCs[J].ACM Transactions on Graphics, 2003,22(3):477-484.

[11] SEDERBERG T W, CARDON D L,FINNIGAN G T, et al. T-spline simplification and local refinement[J]. ACM Transactions on Graphics, 2004, 23(3): 276-283.

[12] GAN W F, FU J Z, SHEN H Y, et al.Five-axis tool path generation in CNC machining of T-spline surfaces[J]. Computer-Aided Design,2014,52:51-63.

[13] DÖRFEL M R, JÜTTLER B,SIMEON B. Adaptive isogeometric analysis by local h-refinement with T-splines[J]. Computer Methods in Applied Mechanics and Engineering,2010,199(5): 264-275.

[14] SCOTT M A, LI X, SEDERBERG T W, et al. Local refinement of analysis-suitable T-splines[J]. Computer Methods in Applied Mechanics and Engineering, 2012,213:206-222.

[15] WANG A, ZHAO G, LI Y D. An influence-knot set based new local refinement algorithm for T-spline surfaces[J]. Expert Systems with Applications, 2014,41(8):3915-3921.

[16] SCOTT M A, BORDEN M J,VERHOOSEL C V, et al. Isogeometric finite element data structures based on Bézier extraction of T-splines[J]. International Journal for Numerical Methods in Engineering,2011,88(2):126-156.

[17] NGUYEN T N. Isogeometric finite element analysis based on Bézier extraction of NURBS and T-Splines[D].Trondheim:Norwegian University of Science and Technology, 2011.

[18] SCOTT M A, SIMPSON R N, EVANS J A, et al. Isogeometric boundary element analysis using unstructured T-splines[J]. Computer Methods in Applied Mechanics and Engineering,2013,254:197-221.