(北京航空航天大學機械工程及自動化學院，北京 100191）

等幾何分析（isogeometry analysis）是一種基于樣條理論的數(shù)值計算方法[1]。它的核心在于將CAD建模工具和有限元求解方式結(jié)合起來，使幾何模型與分析網(wǎng)格模型表達統(tǒng)一，彌補了幾何模型與分析模型巨大差異的缺陷[2]?；诖耍葞缀畏治鼍哂休^高的效率和更精確的分析結(jié)果，目前已成功地應用于振動分析[3-4]、流體力學[5-6]、板殼[7-8]以及電磁[9]等領域。

非均勻有理B樣條曲面（Non-Uniform Rational B-Splines，NURBS）是應用較為廣泛的一種CAD設計模型。一些學者對基于NURBS模型的等幾何分析進行了研究，取得了一些成果。但NURBS模型嚴格的矩形拓撲限制使得其在進行等幾何分析時依舊存在以下不足：首先，使用NURBS曲面表達拓撲結(jié)構(gòu)比較復雜的模型時，不可避免地會出現(xiàn)多張NURBS曲面的裁剪和拼接，增大了等幾何分析中網(wǎng)格劃分的難度；其次，NURBS曲面無法進行局部細化，因此分析者對分析模型的某個局部進行加密時，會不可避免地引入冗余的控制頂點，從而降低計算效率。

Sederberg等提出的T樣條曲面（T-splines）為解決NURBS的不足提供了契機[10-11]。相比較于NURBS曲面，T樣條曲面具有以下兩個優(yōu)點：（1）單張T樣條曲面可表示拓撲結(jié)構(gòu)比較復雜的模型[12]。（2）局部細化[13-15]。因此，T樣條曲面是一種比較適合于進行等幾何分析的模型。

為了推動基于T樣條曲面的等幾何分析發(fā)展，Scott等提出了T樣條曲面的Bézier提取算法[16]。Nguyen將基于Bézier單元的等幾何分析結(jié)果與傳統(tǒng)的有限元分析進行了比較[17]。之后，Scott等又提出了帶奇異點的T樣條曲面的Bézier提取方法[18]。

T樣條曲面的Bézier提取算法有以下兩個關鍵步驟：（1）對參數(shù)域T網(wǎng)格進行Bézier分片，使得劃分后的參數(shù)域網(wǎng)格所對應的空間曲面片（patch）的內(nèi)部C∞連續(xù)。（2）求取T樣條曲面每個混合函數(shù)所對應的Bézier提取矩陣。文獻[16]以及文獻[18]對步驟（2）的具體實現(xiàn)過程進行了詳細說明，但都沒有給出步驟（1）的具體實現(xiàn)方法。為此，本文對參數(shù)域T網(wǎng)格的Bézier分片進行了研究，給出了一種參數(shù)域T網(wǎng)格的Bézier分片算法。該算法有利于后續(xù)進一步實現(xiàn)T樣條曲面的Bézier提取算法，最終達到基于T樣條曲面等幾何分析的目的。

T樣條曲面的基本概念

T樣條曲面可看作NURBS曲面的一般化，它允許控制頂點呈T型連接。本節(jié)將會對T樣條曲面的基本概念進行簡要介紹。

定義1：T樣條曲面的控制網(wǎng)格（T-mesh）稱為空間T網(wǎng)格，它包含了T樣條曲面的所有控制頂點空間位置信息以及拓撲連接信息。

定義2：T樣條曲面上的控制頂點都對應一個參數(shù)點，由于這些參數(shù)點連接而成的參數(shù)域網(wǎng)格與空間T網(wǎng)格同胚，因此稱該參數(shù)域網(wǎng)格為參數(shù)域T網(wǎng)格。

空間T網(wǎng)格與參數(shù)域T網(wǎng)格分別如圖1（a）和（b）所示，第13號控制頂點對應于參數(shù)域上的a13點。

T樣條曲面的計算表達式為：

其中Bi(u,v)=Ni,k(u)Ni,p(v)為控制頂點pi所對應的混合函數(shù)；ωi為控制頂點pi的權(quán)值；Ni，k（u），Ni，p（v）分別為k次及p次B樣條基函數(shù)。T樣條曲面的具體計算過程見文獻[10]。

T樣條曲面的Bézier單元

T樣條曲面在計算過程中以參數(shù)域T網(wǎng)格上的面片作為最小單元，但參數(shù)域T網(wǎng)格上的每個面片所對應的T樣條曲面片（patch）不一定能用一張有理Bézier曲面表示，這里以一個具體實例對該現(xiàn)象進行說明。

圖2為一個雙三次T樣條曲面的參數(shù)域T網(wǎng)格。v1～v6為參數(shù)域T網(wǎng)格中的6個參數(shù)點，它們所對應的控制頂點分別為p1～p6，s1為參數(shù)域T網(wǎng)格中的一個面片，它包含v2～v65個參數(shù)點，陰影區(qū)域為控制頂點p1所對應的混合函數(shù)N1，3（u）N1，3（v）的非零區(qū)域。由B樣條基函數(shù)的可微性可知，面片s1所對應的T樣條曲面片內(nèi)部不能達到C∞連續(xù)，而有理Bézier曲面的內(nèi)部可以達到C∞連續(xù)。因此，面片s1所對應的T樣條曲面片不能用一張有理Bézier曲面表示。

因此，為了實現(xiàn)T樣條曲面的Bézier提取，必須首先將T樣條曲面進行分割，使得分割后的每個T樣條曲面片的內(nèi)部無限次可微，而邊界處則只能達到有限次可微，稱具備上述性質(zhì)的T樣條曲面片為T樣條曲面的一個 Bézier單元。

定義3：T樣條曲面劃分Bézier單元的過程稱為T樣條曲面的Bézier分割。

T樣條曲面的Bézier分割并不是直接在T樣條曲面上進行切分得到，而是通過對參數(shù)域T網(wǎng)格進行分割而間接得到的。由B樣條基函數(shù)的連續(xù)性可知，T樣條曲面會在混合函數(shù)的節(jié)點線所對應的空間曲線處連續(xù)性降低。因此，只需要對參數(shù)域T網(wǎng)格按照T樣條混合函數(shù)的節(jié)點線進行分割，則可以間接實現(xiàn)T樣條曲面的 Bézier分割。

圖1 T樣條曲面的空間T網(wǎng)格和參數(shù)域T網(wǎng)格Fig.1 T-mesh in Cartesian space and parametric space

參數(shù)域T網(wǎng)格Bézier分割算法

定義4：T樣條曲面在參數(shù)域T網(wǎng)格上間接實現(xiàn)T樣條曲面Bézier分割的過程稱為參數(shù)域T網(wǎng)格的Bézier分割。

圖3為一個雙三次T樣條曲面參數(shù)域網(wǎng)格，v1～v5為參數(shù)域T網(wǎng)格中的5個參數(shù)點，它們所對應的控制頂點分別為p1～p5，s1～s7為參數(shù)域 T網(wǎng)格上的7個面片。其中控制頂點p2所對應的混合函數(shù)B2（u,v）的節(jié)點線分布如圖4所示。因此，參數(shù)域上的面片至少需要切分成如圖5所示的參數(shù)域網(wǎng)格。

圖2 參數(shù)域T網(wǎng)格1Fig.2 T-mesh one in parametric space

圖3 參數(shù)域T網(wǎng)格2Fig.3 T-mesh two in parametric space

圖4 混合函數(shù)B2（u,v）的節(jié)點線Fig.4 Knot line of blending function B2（u,v）

圖5 混合函數(shù)B2（u,v）的節(jié)點線分割效果圖Fig.5 T-mesh segmented by the knot line of blending function B2（u,v）

顯然，如果采用蠻力算法求取每一個控制頂點所對應的混合函數(shù)的節(jié)點線分布，最終雖能夠?qū)崿F(xiàn)參數(shù)域T網(wǎng)格的Bézier分割，但會出現(xiàn)很多重復的計算，并且在切分過程中可能出現(xiàn)如圖5右上角所示的L型點。

為減少重復計算，并避免L型分割點的產(chǎn)生。本文給出了包含如下4個步驟的參數(shù)域T網(wǎng)格Bézier分割算法。

算法1：參數(shù)域T網(wǎng)格Bézier分割算法。

STEP1：將參數(shù)域T網(wǎng)格中所有的T型參數(shù)點進行延伸，通過T型參數(shù)點的延伸將參數(shù)域T網(wǎng)格中面片進行切分。

STEP2：遍歷參數(shù)域上的所有參數(shù)點，對于每一個參數(shù)點，執(zhí)行STEP3 ～ STEP4。

STEP3：計算當前參數(shù)點上所耦合的混合函數(shù)的非零區(qū)域。

STEP4：如果STEP3中得到的非零區(qū)域不能完全覆蓋參數(shù)域T網(wǎng)格中的某個面片，則切分該面片，否則不切分該面片。

下面以如圖3所示的參數(shù)域T網(wǎng)格為例對參數(shù)域T網(wǎng)格Bézier分割算法進行解釋。

首先處理T型參數(shù)點，由于圖3所示的參數(shù)域T網(wǎng)格中參數(shù)點v2～v5均為T型連接點，因此應首先對這些T型參數(shù)點進行延伸。如圖6所示，參數(shù)點v2以及v3均向v正方向延伸兩個T網(wǎng)格距離，而參數(shù)點v4以及v5均向u負方向延伸兩個T網(wǎng)格。經(jīng)過這個延伸過程，面片s3切分成子面片f1和f2；面片s4切分成子面片f3和f4；而面片s6切分成子面片f5、f6、f7以及f8。

其次計算參數(shù)點所耦合的混合函數(shù)非零區(qū)域，如圖7所示，v1為參數(shù)域T網(wǎng)格中的十字型參數(shù)點，參數(shù)點v1所耦合的混合函數(shù)非零區(qū)域為見圖7的陰影部分。

由于參數(shù)點v1耦合的混合函數(shù)的非零區(qū)域完全覆蓋了面片f1和f3，因此這兩個面片不需要切分。而面片s5不能被參數(shù)點v1所耦合的混合函數(shù)非零區(qū)域完全覆蓋，因此需要將圖3中的面片s5切分為如圖8所示的兩個子面片f9和f10。

圖6 T型參數(shù)點延伸后的參數(shù)域T網(wǎng)格Fig.6 T-mesh generated by extending T-vertex

圖7 參數(shù)點v1所耦合的混合函數(shù)的非零區(qū)域Fig.7 Nonzero region of the blending function corresponding to v1

經(jīng)過上述4個步驟，參數(shù)域T網(wǎng)格中的每個面片所對應的T樣條曲面均可以由一張有理Bézier曲面表示，也就實現(xiàn)了參數(shù)域T網(wǎng)格Bézier分割過程。

從上述實例中可以看到，參數(shù)域T網(wǎng)格Bézier分割算法可能將一個面片分成兩個面片，如圖3中的s3；也可能分成四個面片，如圖3中的s6。事實上，如果參數(shù)域T網(wǎng)格比較復雜，還會出現(xiàn)一分為三、一分為六等多種情形。為了實現(xiàn)面片切分算法，本文使用如下所示的BézierFace結(jié)構(gòu)體來處理面片切分問題。這個結(jié)構(gòu)體中主要包含了每一個BézierFace范圍信息，以及該BézierFace所屬的面片的索引號信息。

};

基于上述結(jié)構(gòu)，可以使用算法2實現(xiàn)參數(shù)域T網(wǎng)格中面片的具體切分過程。

算法2：參數(shù)域T網(wǎng)格Bézier分割中的面片切分算法。

圖8 參數(shù)域Bézier分割網(wǎng)格Fig.8 Bézier segmentation of the T-mesh

STEP1：建立一個隊列數(shù)組，該數(shù)組的大小與T樣條曲面的參數(shù)域T網(wǎng)格中的面片數(shù)一致，并使用與隊列索引號相同的面片信息初始化該隊列。

圖9 參數(shù)域T網(wǎng)格及其延伸網(wǎng)格Fig.9 T-mesh and its extending T-mesh

圖10 參數(shù)域Bé zier分割網(wǎng)格Fig.10 Bé zier segmentation of the T-mesh

STEP2：在執(zhí)行算法1的過程中，如果出現(xiàn)面片需要切分的情形，則執(zhí)行 STEP3～STEP4。

STEP3：遍歷與當前面片索引號相同的隊列中的所有子面片，對于每一個子面片執(zhí)行STEP4。

STEP4：如果當前子面片需要切分，則建立兩個新面片并將位于隊頭的當前子面片執(zhí)行出隊操作，之后將新建兩個面片分別執(zhí)行入隊操作；如果當前子面片不需要切分，則將位于隊頭的當前子面片先執(zhí)行出隊操作，隨后執(zhí)行入隊操作。

下面以一個具體實例對算法2進行解釋。圖8為一個雙三次T樣條曲面參數(shù)域T網(wǎng)格。v1、v3以及v4為參數(shù)域T網(wǎng)格中的3個T型參數(shù)點，s1為參數(shù)域T網(wǎng)格中的一個面片，其umin、umax、vmin以及vmax分別為 1、4、3和4。參數(shù)域T網(wǎng)格將T型參數(shù)點進行延伸后得到如圖9所示的參數(shù)域延伸T網(wǎng)格。結(jié)合圖9（a）和（b）可知，由于T型參數(shù)點v4的延伸，導致面片s1被切分成面片e1以及e2。在具體程序?qū)崿F(xiàn)時，首先建立兩個新面片e1和e2，其中新面片 e1的umin、umax、vmin以及vmax分別為3、4、3和4；新面片e2的umin、umax、vmin以及vmax分別為 1、3、3和4。之后將隊列中僅有的面片s1進行出隊操作，最后對新面片e1和e2依次執(zhí)行入隊操作，此時隊列中存在e1和e2兩個面片。

圖10中的陰影部分為十字型參數(shù)點v2所耦合的混合函數(shù)的非零區(qū)域，這個區(qū)域完全覆蓋面片e1，但并未完全覆蓋圖9（b）面片e2。因此，首先對隊列中的e1執(zhí)行出隊操作，接著對e1執(zhí)行入隊操作，此時面片e1由隊頭轉(zhuǎn)到隊尾，這個過程相當于建立了面片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依次執(zhí)行入隊操作，此時隊列中存在f1、f2及f33個面片。

參數(shù)域T網(wǎng)格經(jīng)過上述算法的分割后，每一個參數(shù)域面片均所對應的T樣條曲面片均為一個Bézier單元，這就間接實現(xiàn)了T樣條曲面的Bézier分割。

實例驗證

為了驗證本文所提出的參數(shù)域T網(wǎng)格的Bézier分割算法的有效性，在Visual Studio 2012軟件開發(fā)平臺上實現(xiàn)了T樣條曲面的Bézier提取算法。該算法首先采用本文所提出的參數(shù)域T網(wǎng)格的Bézier分割算法對T樣條模型的參數(shù)域T網(wǎng)格進行處理，之后采用文獻[16]提出的混合函數(shù)Bézier提取矩陣求取算法對T樣條模型進行Bézier提取，部分提取結(jié)果如圖11～圖13所示。

這里，以face模型為例，對T樣條模型以及Bézier提取模型的分割線進行說明。如圖11（a）所示，face模型中的每一小曲面片均對應于參數(shù)域T網(wǎng)格中的一個面片，該模型經(jīng)過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

結(jié)論

本文提出了一種參數(shù)域T網(wǎng)格的Bézier分片算法，該算法只需要借助隊列這一常見的數(shù)據(jù)結(jié)構(gòu)就能間接實現(xiàn)T樣條曲面的Bézier分割，有利于后續(xù)完整的實現(xiàn)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] 吳紫俊, 黃正東, 左兵權(quán), 等. 等幾何分析研究概述[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.