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

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

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

Sederberg等提出的T樣條曲面（T-splines）為解決NURBS的不足提供了契機[10-11]。相比較于NURBS曲面，T樣條曲面具有以下兩個優(yōu)點：（1）單張T樣條曲面可表示拓?fù)浣Y(jié)構(gòu)比較復(fù)雜的模型[12]。（2）局部細(xì)化[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提取算法有以下兩個關(guān)鍵步驟：（1）對參數(shù)域T網(wǎng)格進行Bézier分片，使得劃分后的參數(shù)域網(wǎng)格所對應(yīng)的空間曲面片（patch）的內(nèi)部C∞連續(xù)。（2）求取T樣條曲面每個混合函數(shù)所對應(yīng)的Bézier提取矩陣。文獻(xiàn)[16]以及文獻(xiàn)[18]對步驟（2）的具體實現(xiàn)過程進行了詳細(xì)說明，但都沒有給出步驟（1）的具體實現(xiàn)方法。為此，本文對參數(shù)域T網(wǎng)格的Bézier分片進行了研究，給出了一種參數(shù)域T網(wǎng)格的Bézier分片算法。該算法有利于后續(xù)進一步實現(xiàn)T樣條曲面的Bézier提取算法，最終達(dá)到基于T樣條曲面等幾何分析的目的。

T樣條曲面的基本概念

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

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

定義2：T樣條曲面上的控制頂點都對應(yīng)一個參數(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號控制頂點對應(yīng)于參數(shù)域上的a13點。

T樣條曲面的計算表達(dá)式為：

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

T樣條曲面的Bézier單元

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

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

因此，為了實現(xiàn)T樣條曲面的Bézier提取，必須首先將T樣條曲面進行分割，使得分割后的每個T樣條曲面片的內(nèi)部無限次可微，而邊界處則只能達(dá)到有限次可微，稱具備上述性質(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é)點線所對應(yīng)的空間曲線處連續(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ù)點，它們所對應(yīng)的控制頂點分別為p1～p5，s1～s7為參數(shù)域 T網(wǎng)格上的7個面片。其中控制頂點p2所對應(yīng)的混合函數(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）

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

為減少重復(fù)計算，并避免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：計算當(dāng)前參數(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型連接點，因此應(yīng)首先對這些T型參數(shù)點進行延伸。如圖6所示，參數(shù)點v2以及v3均向v正方向延伸兩個T網(wǎng)格距離，而參數(shù)點v4以及v5均向u負(fù)方向延伸兩個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)格中的每個面片所對應(yī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)格比較復(fù)雜，還會出現(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：遍歷與當(dāng)前面片索引號相同的隊列中的所有子面片，對于每一個子面片執(zhí)行STEP4。

STEP4：如果當(dāng)前子面片需要切分，則建立兩個新面片并將位于隊頭的當(dāng)前子面片執(zhí)行出隊操作，之后將新建兩個面片分別執(zhí)行入隊操作；如果當(dāng)前子面片不需要切分，則將位于隊頭的當(dāng)前子面片先執(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的延伸，導(dǎo)致面片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)到隊尾，這個過程相當(dāng)于建立了面片f1。之后，隊列中的第二個面片e2，顯然需要進行切分，因此應(yīng)新建兩個面片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ù)域面片均所對應(yīng)的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)格進行處理，之后采用文獻(xiàn)[16]提出的混合函數(shù)Bézier提取矩陣求取算法對T樣條模型進行Bézier提取，部分提取結(jié)果如圖11～圖13所示。

這里，以face模型為例，對T樣條模型以及Bézier提取模型的分割線進行說明。如圖11（a）所示，face模型中的每一小曲面片均對應(yīng)于參數(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提取算法。在此基礎(chǔ)上，未來可進一步研究基于T樣條曲面Bézier單元的等幾何分析算法。

參 考 文 獻(xiàn)

[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]. 機械工程學(xué)報, 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.