混合B樣條實體模型的等幾何拓撲優化

楊佳明，趙 罡,2，王 偉,3，郭馬一，杜孝孝

混合B樣條實體模型的等幾何拓撲優化

楊佳明1，趙 罡1,2，王 偉1,3，郭馬一1，杜孝孝1

(1. 北京航空航天大學機械工程及自動化學院，北京 100191；2. 航空高端裝備智能制造技術工業和信息化部重點實驗室，北京 100191；3. 北京市高效綠色數控加工工藝及裝備工程技術研究中心，北京 100191)

等幾何拓撲優化方法將經典拓撲優化理論中的有限元分析過程更改為等幾何分析計算，從而提高了拓撲優化的效率與穩定性。針對現有的等幾何拓撲優化方法在處理復雜實體結構優化問題時具有一定的局限性，提出一種非結構化樣條實體等幾何拓撲優化方法。基于混合B樣條構造技術，在非結構化六面體網格上構造具有復雜結構的樣條實體，并將其作為拓撲優化問題的設計域。用于描述這一樣條實體的基函數被直接應用于材料密度分布的表達以及等幾何分析計算。在數值算例中，該方法表現出應用于復雜結構時的良好穩定性和魯棒性。研究成果對等幾何拓撲優化方法應用于實際工程問題具有一定的參考意義。

拓撲優化；等幾何分析；體參數化；B樣條；非結構化樣條

拓撲優化是一種重要的工程結構優化方法。目前常見的拓撲優化方法包括：均勻化方法[1]、固體各向同性懲罰微結構模型(solid isotropic microstructures with penalization，SIMP)法[2-3]、漸進結構優化方法[4]、水平集方法[5-7]、移動組件法[8-9]等。由于這些方法可以在設計域中按照給定約束條件高效地搜索最優的設計結構，因此被廣泛應用于各類工程問題。其中，AAGE等[10]將千兆體素級的SIMP法應用于全尺寸機翼的結構設計并將該成果發表于《Nature》期刊上，體現了拓撲優化技術的應用潛力。

這些經典的拓撲優化方法一般都會通過有限元方法來實現結構響應的計算。但傳統的有限元方法存在以下問題和局限性：①低階連續性的有限元有可能影響優化計算的準確性并導致數值計算的不穩定 性[11]；②繁瑣耗時的有限元前處理過程割裂了幾何設計與力學分析過程，從而使拓撲優化難以作為一種設計工具完全融入到現有的計算機輔助設計系統之中。

HUGHES等[12-13]推廣的等幾何分析方法作為有限元方法，被認為有可能突破上述局限性。其核心思想是直接利用CAD中描述幾何模型的樣條基函數來表達力學分析中的未知場。由于該方法統一了設計與分析過程中模型的基礎表達形式，避免了模型轉換引起的誤差，因此具有諸多優勢。國內外很多學者將其應用于拓撲優化，從而形成了等幾何拓撲優化方法。其中，SEO等[14]最先提出在拓撲優化中使用等幾何分析方法及裁剪樣條曲面；隨后，HASSANI等[15]提出了在點密度SIMP法基礎上引入NURBS基函數的等幾何拓撲優化方法；QIAN[16]將設計域嵌入到一個B樣條空間中，證明了基于樣條函數的材料密度分布表示方法不僅具有良好的效率和魯棒性，并且具有固有的過濾性質。國內WANG和BENSON[17]將等幾何分析引入到水平集方法中，并針對等幾何拓撲優化方法提出了一種高效的計算框架[18]；XIE等[19-20]探索了等幾何拓撲優化在移動組件法中的應用；GAO等[21-22]對基于NURBS的等幾何拓撲優化方法及相關應用開展了深入的研究。

可以看到，現有的等幾何拓撲優化研究大多基于經典的樣條構造方法，如B樣條、NURBS。而這些張量積樣條方法受制于矩形的拓撲結構。盡管可以使用裁剪、拼接等造型方法進行處理，但要保證復雜設計域模型的水密性依舊是一件極具挑戰的任務。現有的等幾何拓撲優化方法大多只能在簡單的矩形設計域上進行驗證，制約了在復雜問題上的探索。尤其對于三維實體問題，一些在復雜平面設計域問題中表現良好的方法難以直接推廣至實體層面，如ZHAO等[23]提出的基于非結構化T樣條曲面的等幾何拓撲優化方法。

針對上述問題，本文將一種可以在非結構化體網格上構造樣條實體的方法引入等幾何拓撲優化。由于允許奇異邊(點)的存在，因此可以對高虧格的復雜結構進行表達，突破傳統方法的矩形拓撲限制，更符合拓撲優化過程的需求。本文使用WEI等[24]提出的混合B樣條構造方法。此方法構造的樣條函數不僅具有優良的性質，如非負性、規范性、線性無關性，并且構造過程可靠，具有很好的適應性。本文基于此方法構造非結構化的樣條實體模型，并給出了在此類模型上實現等幾何拓撲優化的方法，通過算例驗證了其可行性，及對復雜實體結構的處理能力。該研究成果有望推動等幾何拓撲優化方法在實際工程問題中的應用與推廣。

1 混合B樣條實體構造方法

本節簡要回顧B樣條及非結構化六面體網格上的混合B樣條構造方法，相關細節可以參考文獻[24-25]。

1.1 B樣條

雙變量/三變量B樣條基函數可以通過單變量B樣條基函數的張量積得到。

1.2 混合B樣條構造方法

其中，為中的一部分控制頂點。

由節點插入算法[27]可知，，和0之間存在線性變換關系

其中，和可以直接通過節點插入算法得到。將其代入式(4)可以得到基函數之間的線性關系

圖2 三變量B樣條基函數的截斷((a)一個C2控制頂點；(b)基函數的等值面；(c) C2控制點周圍存在一個激活的C0控制頂點；(d)從中截取后得到的截斷基函數的等值面)

文獻[24]將上述情況稱為02構造。除此之外，根據激活的0，1控制頂點不同還存在012和02構造形式。相比于012和02，本文采用的02方法盡管在奇異點鄰域的連續性上存在不足，但具有構造形式簡單以及基函數局部線性無關等優勢。

1.3 C02混合B樣條實體構造方法

一個非結構化六面體網格中包含了頂點、邊、面和體單元4類元素。定義邊(頂點)的價為：包含這條邊(頂點)的體單元的個數。定義奇異邊為：價不為4的內部邊以及價不為1或2的邊界邊。定義奇異點為：奇異邊的端點。

如果一個體單元的頂點中包含了奇異點，則稱這個體單元為不規則單元。其余的體單元稱為規則單元。規則單元又可分為2類：1-鄰域內存在不規則單元的規則單元稱為規則過渡單元，否則稱為規則非過渡單元。圖3為一個包含三價奇異邊的簡單非結構化體網格，并通過顏色區分了不同的單元類型。

圖3 非結構化體網格及3種單元類型

1.3.1 不規則單元

其中，0包含了64個三-三次貝奇爾基函數。

1.3.2 規則非過渡單元

1.3.3 規則過渡單元

圖4 eRT與相鄰非規則單元之間不同的相對位置關系

1.3.4 統一構造形式

綜上，由式(13)、式(15)和式(17)，可將非結構化六面體網格上的02混合B樣條實體構造總結為

可以發現，3類單元均可寫為統一的形式，即

1.3.5 連續性

按照式(19)構造的樣條實體模型，在規則單元間為2連續，在非規則單元間為0連續，在非規則單元和規則單元間為0連續。02混合B樣條實體的連續性分析的簡要證明如下：

為了討論過渡單元與其相鄰規則單元間的連續性，引入以下結論：

引理[24]1. 對于任意一個固定的單元，其混合B樣條表達形式(式(12))等價于這個單元的B樣條表達形式(式(4))。

為了使模型插值于邊界，本文將邊界單元均按照非規則單元處理。因此，此時邊界曲面為0連續曲面。為了提高邊界曲面的連續性，可使用非結構化四邊形網格的012混合構造方法對其進行處理。當存在原始的B-rep模型時，也可通過調整邊界控制頂點使其恢復原始曲面的光順性。

2 等幾何拓撲優化

本節對混合B樣條實體的等幾何拓撲優化方法進行介紹。為了簡化問題的討論同時又不失一般性，本文著重考慮拓撲優化問題中的最小化柔順度問題。

2.1 基于混合B樣條基函數的材料密度分布

其中，0為材料本身的楊氏模量；min為設定的最小剛度值，用以避免剛度矩陣奇異的情況；為懲罰因子。

2.2 等幾何分析

在等幾何方法中，可以將靜態載荷下線彈性問題的剛度方程表示為

其中，為幾何映射(e)的雅可比矩陣。同時，可以通過式(23)對物理域坐標的偏導求得；則通過式(22)表示為

其中，0為楊氏模量等于1時的彈性系數矩陣。

使用高斯積分法來完成式(26)的數值積分運算，則可表示為

2.3 拓撲優化

2.4 靈敏度分析

則目標函數的靈敏度，即式(31)可寫為

2.5 優化流程

圖5 混合B樣條實體等幾何拓撲優化流程圖

3 數值算例

3.1 算例1

圖6(b)為用于構造此設計域樣條實體模型的非結構化六面體網格，通過彩色點標出了模型表面的部分奇異點。可以看到本文所使用的構造方法可以支持邊界奇異點，以及奇異點相鄰的情況。圖6(c)和6(d)分別給出了本文構造的02混合B樣條實體及其控制頂點分布。此模型的控制頂點數(同時也為設計變量數)為17 745，單元數為540。可以看到模型內部包含了一些奇異邊(點)。當直接使用傳統張量積樣條方法進行構造時，需要考慮多塊拼接、裁剪等問題。而本文使用的02混合構造方法在處理時則較為簡潔。

圖6 算例1的設計域及混合B樣條實體模型((a)設計域；(b)用于構造混合B樣條實體的非結構化六面體網格；(c)構造的混合B樣條實體模型及等幾何分析結果(x向位移)；(d)樣條實體模型的控制點)

等幾何拓撲優化結果如圖7(a)所示，最終目標函數值為60.283 8，迭代次數為78。作為對比，圖7(b)中展示了SIMP方法的優化結果。該模型的單元數為22 680，最終目標函數值為66.488 1。相比于SIMP方法，本文方法以較少的單元數得到了較好的優化結果(更低的柔順度)。主要是由于在SIMP方法中，材料密度分布通過離散單元的形式來表達，每一個單元對應一個單一的密度值。而本文方法以連續的樣條函數來表達材料密度分布，每個單元都具有內部的密度分布函數，且單元間密度分布連續，因此可以通過較少的單元表達復雜的材料密度分布結構。

圖7 算例1的拓撲優化結果((a)混合B樣條實體等幾何拓撲優化結果；(b) SIMP法優化結果)

圖8為本例中等幾何拓撲優化的收斂曲線以及部分中間結果。可以看到整個優化過程的收斂速度較快，且較為光滑，沒有出現明顯的振蕩。說明本文的方法具有較好的收斂性質。

圖8 算例1的優化收斂曲線

3.2 算例2

本例為一個復雜實體結構的拓撲優化問題。設計域及邊界條件如圖9(a)所示，圖中紅色圓柱面固定，3個角上的藍色平面施加垂直于平面向外的均布載荷，大小為1。本文將體積分數設為40%，設計變量的初始值為0.4。

圖9(b)為用于構造此設計域樣條實體模型的非結構化六面體網格。在圖9(c)中展示了在樣條實體構造過程中識別的不同單元類型。其中，紅色為非規則單元，白色為規則單元。規則單元間具有2連續性。圖9(d)為最終構造的02混合B樣條實體及等幾何分析結果。可以看到位移場在模型內部具有較高的連續性。該模型共有162 303個控制頂點，如圖10(a)所示。如果完全采用貝奇爾提取方法來處理此樣條實體模型，則控制頂點數為197 895，如圖10(b)所示。由于02構造方法在模型內部具有更高的連續性，因此使用了更少的控制頂點，其分布更為稀疏，同時減少了設計變量的規模。

圖9 算例2的設計域及混合B樣條實體模型((a)設計域；(b)用于構造混合B樣條實體的非結構化六面體網格；(c)不同的單元類型：紅色為非規則單元，白色為規則單元；(d)構造的混合B樣條實體模型及等幾何分析結果(x向位移))

圖11為本例中等幾何拓撲優化的結果，最終目標函數值為112.159，迭代次數為100。圖12為目標函數收斂曲線、體積分數收斂曲線，以及部分中間結果。可以看到，同算例1，在本例中優化收斂速度較快，無明顯振蕩，說明本文方法在處理復雜結構時也具有較好的穩定性和魯棒性。

圖10 混合B樣條實體與完全采用貝奇爾提取方法時的控制頂點分布對比，圖中每一個紅色點表示一個控制頂點，紅色點越密集的區域控制頂點數目越多((a)混合B樣條實體的控制頂點分布；(b)完全采用貝奇爾提取方法時的控制頂點分布)

圖11 算例2的等幾何拓撲優化結果

圖12 算例2的優化收斂曲線

3.3 算例3

本例為一個多工況復雜實體結構的拓撲優化問題。設計域及2個工況狀態如圖13(a)和圖13(b)所示。圖中4個紅色圓柱面固定。在工況1時兩側懸臂的藍色圓柱面區域分別受到軸正向和負向的載荷，大小為1。在工況2時兩側懸臂的藍色圓柱面區域分別受到軸正向和負向的載荷，大小為1。取體積分數約束為30%。

按照本文方法所構造的樣條實體如圖13(c)所示，圖中同時顯示了工況1下的等幾何分析結果。此模型的控制頂點分布如圖13(d)所示，總個數為110 283。可以看到由于此模型包含了多個孔洞結構，因此如果直接使用張量積NUBRS實體進行構造將會遇到很大困難。

圖13 算例3的設計域及混合B樣條實體模型((a)設計域及工況1；(b)設計域及工況2；(c)構造的混合B樣條實體模型及等幾何分析結果(z向位移)；(d)樣條實體模型的控制點)

在優化過程中本文設定目標函數為2種工況下柔順度的平均值。等幾何拓撲優化結果如圖14所示，最終目標函數值為141.173，其中工況1的柔順度為267.901，工況2為14.445，迭代次數為100次。此算例進一步證明了該方法在復雜實體結構拓撲優化問題中具有良好的適應性。

圖14 算例3的等幾何拓撲優化結果

4 結束語

本文針對現有等幾何拓撲優化方法難以處理復雜實體結構設計域這一問題，提出了一種基于混合B樣條構造的非結構化樣條實體等幾何拓撲優化方法。此方法可在非結構化六面體網格上構造具有一定內部連續性的樣條實體。本文將構造得到的樣條基函數直接用于描述拓撲優化問題中的材料密度分布，以及結構響應計算中的未知場。給出了等幾何拓撲優化表達式以及相應的靈敏度計算公式。最后，通過3個算例驗證了本文方法在處理復雜結構優化問題時的有效性和穩定性。在這一研究成果基礎之上，有望形成面向實際工程中復雜模型的等幾何拓撲優化設計工具。在后續研究中，如何提高奇異點鄰域的連續性將會是一個關鍵的研究方向。

[1] BENDS?E M P, Kikuchi N. Generating optimal topologies in structural design using a homogenization method[J]. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197-224.

[2] Zhou M, Rozvany G I N. The COC algorithm, part II: topological, geometry and generalized shape optimization[J]. Computer Methods in Applied Mechanics and Engineering, 1991, 89(1-3): 309-336.

[3] Bends?e M P, Sigmund O. Material interpolation schemes in topology optimization[J]. Archive of Applied Mechanics, 1999, 69(9-10): 635-654.

[4] Xie Y M, Steven G P. A simple evolutionary procedure for structural optimization[J]. Computers & Structures, 1993, 49(5): 885-896.

[5] Sethian J A, Wiegmann A. Structural boundary design via level set and immersed interface methods[J]. Journal of Computational Physics, 2000, 163(2): 489-528.

[6] Wang M Y, Wang X, Guo D. A level set method for structural topology optimization[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1-2): 227-246.

[7] Allaire G, Jouve F, Toader A M. Structural optimization using sensitivity analysis and a level-set method[J]. Journal of Computational Physics, 2004, 194(1): 363-393.

[8] Guo X, Zhang W S, Zhong W L. Doing topology optimization explicitly and geometrically - a new moving morphable components based framework[J]. Journal of Applied Mechanics, 2014, 81(8): 081009.

[9] Zhang W S, Song J F, Zhou J H, et al. Topology optimization with multiple materials via Moving Morphable Component (MMC) method[J]. International Journal for Numerical Methods in Engineering, 2017, 113(11): 1653-1675.

[10] Aage N, Andreassen E, Lazarov B S, et al. Giga-voxel computational morphogenesis for structural design[J]. Nature, 2017, 550: 84-86.

[11] Sigmund O, Petersson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima[J]. Structural and Multidisciplinary Optimization, 1998, 16(1): 68-75.

[12] 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-41): 4135-4195.

[13] Cottrell J A, Hughes T J R, Bazilevs Y. Isogeometric analysis: toward Integration of CAD and FEA[M]. Chichester: John Wiley & Sons, 2009: 1-18.

[14] Seo Y D, Kim H J, Youn S K. Isogeometric topology optimization using trimmed spline surfaces[J]. Computer Methods in Applied Mechanics and Engineering, 2010, 199(49-52): 3270-3296.

[15] Hassani B, Khanzadi M, Tavakkoli S M. An isogeometrical approach to structural topology optimization by optimality criteria[J]. Structural and Multidisciplinary Optimization, 2012, 45(2): 223-233.

[16] Qian X P. Topology optimization in B-spline space[J]. Computer Methods in Applied Mechanics and Engineering, 2013, 265(1): 15-35.

[17] Wang Y J, Benson D J. Isogeometric analysis for parameterized LSM-based structural topology optimization[J]. Computational Mechanics, 2016, 57(1): 19-35.

[18] Wang Y, Liao Z Y, Ye M, et al. An efficient isogeometric topology optimization using multilevel mesh, MGCG and local-update strategy[J]. Advances in Engineering Software, 2020, 139: 102733.

[19] Xie X D, WANG S T, Xu M M, et al. A new isogeometric topology optimization using moving morphable components based on R-functions and collocation schemes[J]. Computer Methods in Applied Mechanics and Engineering, 2018, 339(1): 61-90.

[20] Gai Y D, Zhu X F, Zhang Y J, et al. Explicit isogeometric topology optimization based on moving morphable voids with closed B-spline boundary curves[J]. Structural and Multidisciplinary Optimization, 2019, 61: 963-982.

[21] Gao J, Luo Z, Xiao M, et al. A NURBS-based Multi-Material Interpolation (N-MMI) for isogeometric topology optimization of structures[J]. Applied Mathematical Modelling, 2020, 81: 818-843.

[22] Gao J, Xue H P, Gao L, et al. Topology optimization for auxetic metamaterials based on isogeometric analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2019, 352(1): 211-236.

[23] Zhao G, Yang J M, Wang W, et al. T-splines based isogeometric topology optimization with arbitrarily shaped design domains[J]. Computer Modeling in Engineering and Sciences, 2020, 123(3): 1033-1059.

[24] Wei X D, Zhang Y J, Toshniwal D, et al. Blended B-spline construction on unstructured quadrilateral and hexahedral meshes with optimal convergence rates in isogeometric analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2018, 341(1): 609-639.

[25] Wei X D, Zhang Y J, Hughes T J R. Truncated hierarchical tricubic C0 spline construction on unstructured hexahedral meshes for isogeometric analysis applications[J]. Computers & Mathematics with Applications, 2017, 74(9): 2203-2220.

[26] Giannelli C, Jüttler B, Speleers H. THB-splines: the truncated basis for hierarchical splines[J]. Computer Aided Geometric Design, 2012, 29(7): 485-498.

[27] Boehm W. Inserting new knots into B-spline curves[J]. Computer Aided Design, 1980, 12(4): 199-201.

[28] 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.

[29] Bends?e M P, Sigmund O. Topology optimization: theory, methods, and applications[M]. Heidelberg: Springer, 2004: 9-28.

[30] Svanberg K. The method of moving asymptotes - a new method for structural optimization[J]. International Journal for Numerical Methods in Engineering, 1987, 24(2): 359-373.

Isogeometric topology optimization of blended B-spline solid model

YANG Jia-ming1, ZHAO Gang1,2, WANG Wei1,3, GUO Ma-yi1, DU Xiao-xiao1

(1. School of Mechanical Engineering & Automation, Beihang University, Beijing 100191, China;2. Key Laboratory of Aeronautics Smart Manufacturing, Ministry of Industry and Information Technology, Beijing 100191, China;3. Beijing Engineering Technological Research Center of High-Efficient & Green CNC Machining Process and Equipment, Beijing 100191, China)

For isogeometric topology optimization (ITO) methods, isogeometric analysis (IGA) is adopted for topology optimization to address the limitation of the finite element method, which can improve the efficiency and stability of the optimization. However, it is of great challenge for existing ITO methods to manage arbitrarily shaped design domains, especially in three-dimensional solid problems. Therefore, a new ITO method was proposed to handle unstructured solid models. A spline solid with complex structures was obtained from an unstructured hexahedral mesh based on the blended B-spline construction. The basis functions describing the unstructured spline solid were applied to the representation of density distribution and the calculation of IGA. Several examples proved the flexibility and robustness of the proposed method in dealing with complex structures. These results may shed light on the application of ITO in practical engineering problems.

topology optimization; isogeometric analysis; volume parameterization; B-spline; unstructured splines

TP 391

10.11996/JG.j.2095-302X.2021030501

A

2095-302X(2021)03-0501-10

2020-09-17；

2020-10-15

17 September，2020；

15 October，2020

國家自然科學基金項目(61972011，61572056)

National Natural Science Foundation of China (61972011, 61572056)

楊佳明(1995-)，男，北京人，博士研究生。主要研究方向為等幾何拓撲優化、T樣條。E-mail：yangjiaming@buaa.edu.cn

YANG Jia-ming (1995-), male, PhD candidate. His main research interests cover isogeometric topology optimization and T-splines. E-mail: yangjiaming@buaa.edu.cn

王 偉(1978-)，男，河北衡水人，副教授，博士。主要研究方向為CAD/CAE、智能計算等。E-mail：jrrt@buaa.edu.cn

WANG Wei (1978-), male, associate professor, Ph.D. His main research interests cover CAD/CAE, intelligence computation, etc. E-mail：jrrt@buaa.edu.cn