丁延冬 羅年猛 楊奧迪 王書亭 朱浩然 謝賢達

華中科技大學(xué)機械科學(xué)與工程學(xué)院，武漢，430074

0 引言

拓撲優(yōu)化是一種在給定條件約束的設(shè)計域進行結(jié)構(gòu)設(shè)計的有效工具，被廣泛應(yīng)用于散熱結(jié)構(gòu)[1]、柔性機構(gòu)[2]與彈性結(jié)構(gòu)[3]。經(jīng)典拓撲優(yōu)化方法主要有各向同性懲罰微結(jié)構(gòu)法(simple isotropic material with penalization，SIMP)[4]、進化法[5]、水平集法[6-7]、移動可變形組件/孔洞法[8-9]等。有限元法是傳統(tǒng)拓撲優(yōu)化的求解基礎(chǔ)，其求解精度對拓撲優(yōu)化過程的魯棒性至關(guān)重要，而有限元物理場離散形函數(shù)的C0連續(xù)無法滿足應(yīng)力約束等高階拓撲優(yōu)化問題的求解要求。

無論是傳統(tǒng)有限元拓撲優(yōu)化程序還是等幾何拓撲優(yōu)化程序，更精確的優(yōu)化結(jié)果均需要規(guī)模更大的有限元求解，因此優(yōu)化過程需消耗大量的計算成本[17]。通過改進有限元求解算法來提高拓撲優(yōu)化計算效率一直是拓撲優(yōu)化的研究熱點，主要方法有自適應(yīng)網(wǎng)格法[18-19]和多重網(wǎng)格法[20-22]。WANG等[23]將等幾何分析求解精度優(yōu)勢與多重網(wǎng)格求解效率的優(yōu)勢相結(jié)合，提出了多重網(wǎng)格等幾何優(yōu)化模型。然而，由于B樣條基函數(shù)空間的非單元一致性，上述多重網(wǎng)格等幾何拓撲優(yōu)化方法在每一層的拓撲優(yōu)化過程中需要進行所有B樣條單元剛度的計算和存儲，進而實現(xiàn)優(yōu)化過程中等幾何單元剛度矩陣的快速更新，導(dǎo)致多重等幾何拓撲優(yōu)化存在數(shù)據(jù)存儲負擔(dān)過重和預(yù)處理時間過長等問題。針對上述問題，可采用Bézier提取算子對B樣條單元數(shù)據(jù)進行預(yù)處理，實現(xiàn)任一層級網(wǎng)格B樣條基函數(shù)通過C0連續(xù)且單元空間保持的Bernstein基函數(shù)來等效表達。其中，Bézier提取算子作為一種基本的轉(zhuǎn)換方式，被用于保證B樣條基函數(shù)的高階連續(xù)性，已被廣泛應(yīng)用于等幾何分析框架中的T樣條[24]、截斷層次B樣條[25-27]與極樣條[28]。

本文提出基于Bézier單元剛度映射的多重網(wǎng)格等幾何拓撲優(yōu)化模型，以實現(xiàn)等幾何拓撲優(yōu)化數(shù)據(jù)的存儲結(jié)構(gòu)和預(yù)處理過程效率的優(yōu)化。通過Bézier提取算子將各層級的所有B樣條單元剛度矩陣的大規(guī)模存儲轉(zhuǎn)換為標準剛度矩陣與各層級Bézier提取矩陣的存儲，并將任一層級所有B樣條單元剛度矩陣的計算過程轉(zhuǎn)換為任一層級若干個(2或3)參數(shù)方向的Bézier提取矩陣的計算過程，進而建立任一層級任一B樣條單元剛度矩陣與標準Bézier單元剛度矩陣的映射關(guān)系。利用B樣條單元剛度矩陣與標準Bézier單元剛度矩陣之間的映射法則，顯著降低多重網(wǎng)格等幾何拓撲優(yōu)化的剛度矩陣內(nèi)存占用并減少剛度矩陣預(yù)處理流程的時間消耗。最后，通過二維與三維數(shù)值算例的優(yōu)化結(jié)果展示本文模型在內(nèi)存空間和預(yù)處理時間兩個方面的優(yōu)勢。

1 理論基礎(chǔ)

首先介紹Bernstein多項式，接著對開節(jié)點向量生成B樣條的遞歸定義進行回顧，最后基于Bézier分解算法引出Bézier提取算子。其中，Bézier提取算子是實現(xiàn)Bernstein多項式與B樣條之間映射的理論基礎(chǔ)。

1.1 Bernstein多項式

對于任意位于[0，1]內(nèi)的參數(shù)坐標ζ，相應(yīng)的p階Bernstein多項式由以下遞推式獲得：

Bi，p(ζ)=

(1)

其中，Bi，0(ζ)=1。階次p=1，2，3的Bernstein基函數(shù)如圖1所示。

(a)p=1 (b)p=2 (c)p=3

Bézier曲線可以由Bernstein多項式作為基函數(shù)并與控制點Pi的線性組合來表達，表達式為

(2)

1.2 B樣條

對于給定的開節(jié)點向量Ξ=(ζ1，ζ2，…，ζn+p+1)，B樣條的基函數(shù)由Cox-de Boor遞推公式定義：

Nj，p(ζ)=

(3)

其中，p、n分別表示B樣條基函數(shù)的階次和相應(yīng)控制點個數(shù)。通過組合B樣條基函數(shù)與控制點可生成B樣條曲線，B樣條曲線的表達式如下：

(4)

多維參數(shù)空間下，B樣條基函數(shù)由張量積結(jié)構(gòu)生成，因此，二維B樣條基函數(shù)由以下公式定義：

N(j1，j2)，(p1，p2)(ζ，η)=Nj1，p1(ζ)Nj2，p2(η)

(5)

其中，(ζ，η)表示二維參數(shù)空間下的參數(shù)坐標，jd(d=1，2)、pd(d=1，2)分別表示d維參數(shù)方向基函數(shù)的索引和階次。開節(jié)點向量(0,0,0,1,2,3,3,3)的B樣條基函數(shù)空間和各節(jié)點向量單元的局部B樣條基函數(shù)空間如圖2所示，圖3為圖2所示B樣條基函數(shù)生成的B樣條曲線。

圖2 開節(jié)點向量(0，0，0，1，2，3，3，3)的2階B樣條基函數(shù)

圖3 圖2中B樣條基函數(shù)所生成的B樣條曲線

1.3 Bézier提取算子

(6)

其中，αj的推導(dǎo)公式為

(7)

C(k)=

(8)

C(k)∈R(n+k-1)×(n+k)

整體的Bézier提取矩陣由C(k)(k=1，2，…，m)組成，其表達式為

Cext=C(1)C(2)…C(m)

(9)

N(ζ)=CextB(ζ)

(10)

Bézier曲線的控制點PB也可以通過Bézier提取算子從相應(yīng)B樣條控制點P獲取，如下式：

(11)

其中，PB的維度為(n+m)×d；P的維度為n×d。

借助于B樣條的張量積結(jié)構(gòu)，多變量的Bézier提取算子可方便地從單變量Bézier提取算子生成。雙變量情況下，Bézier提取算子的表達式如下：

C(ζ，η)=Cext(η)?Cext(ζ)

(12)

階次p=2時，圖2中任意單元的局部B樣條基函數(shù)均可映射為圖1b中具有空間保持特點的Bernstein基函數(shù)。圖4所示為圖2中單元3的B樣條基函數(shù)空間與Bernstein基函數(shù)空間之間的映射過程。圖5所示為B樣條基函數(shù)空間與相應(yīng)的Bernstein基函數(shù)空間，具體的映射關(guān)系表達如下：

圖4 圖2中單元3的基函數(shù)與圖1b的Bernstein基函數(shù)的映射關(guān)系

圖5 圖2中B樣條基函數(shù)與Bernstein基函數(shù)映射關(guān)系

(13)

2 多重網(wǎng)格等幾何拓撲優(yōu)化模型

首先介紹B樣條空間下等幾何拓撲優(yōu)化模型，在此基礎(chǔ)上對傳統(tǒng)多重網(wǎng)格等幾何拓撲優(yōu)化模型[23]進行回顧；接著，基于Bézier提取算子對傳統(tǒng)等幾何拓撲優(yōu)化模型進行重構(gòu)，獲得改進多重網(wǎng)格模型；最后，通過多重網(wǎng)格各層設(shè)計變量之間的映射關(guān)系，完成改進多重網(wǎng)格等幾何拓撲優(yōu)化模型的靈敏度分析。

2.1 傳統(tǒng)多重網(wǎng)格等幾何拓撲優(yōu)化模型

給定圖6所示的二維矩形設(shè)計域，等幾何拓撲優(yōu)化的設(shè)計目標是在約束條件下尋找最優(yōu)的材料分布，即尋找具有最小應(yīng)變能的最優(yōu)結(jié)構(gòu)，目標函數(shù)表示為

圖6 B樣條空間下的6×2網(wǎng)格的二維懸臂梁設(shè)計域

c(x)=fTu(x)

(14)

式中，x為所有B樣條單元的相對材料密度分布向量；c(x)為在材料分布x下的模型結(jié)構(gòu)柔度；f為外部力向量；u(x)為各個控制點上的位移向量。

綜合SIMP方法，B樣條空間下等幾何拓撲優(yōu)化的數(shù)學(xué)模型可以表示為

(15)

式中，ne為設(shè)計域離散成B樣條網(wǎng)格的單元個數(shù)；xe、ue、Ke(xe)分別為第e個單元的相對密度、位移向量和剛度矩陣；K(x)為全局剛度矩陣；V(x)為固體材料的體積；V0為設(shè)計域的總體積；vf為模型限制的材料體積分數(shù)；為x的可容許空間。

式(15)中的Ke(xe)可由下式獲得：

(16)

Ee(xe)=Emin+(xe)p(E0-Emin)xe∈[0，1]

(17)

式中，E0為實體材料的彈性模量；Emin為空白單元的彈性模量，通常取0

基于多重網(wǎng)格的等幾何拓撲優(yōu)化模型是一種平衡優(yōu)化效率與優(yōu)化精度的有效手段，具體過程為：首先，等幾何拓撲優(yōu)化在較粗的網(wǎng)格下進行優(yōu)化，此時迭代的速度較快、計算成本較低；然后，當設(shè)計變量變化較小時，將粗網(wǎng)格優(yōu)化結(jié)果映射為更細網(wǎng)格的初始設(shè)計，并在細網(wǎng)格上進行等幾何分析和優(yōu)化；最后，重復(fù)以上操作若干次，由最精細的網(wǎng)格生成等幾何拓撲優(yōu)化結(jié)果。

多重網(wǎng)格優(yōu)化框架如圖7所示，其中，等幾何拓撲優(yōu)化部分是基礎(chǔ)部分，是迭代過程中的迭代體；設(shè)計變量映射繼承部分是算法實現(xiàn)的核心，是迭代過程中完成代際繼承的關(guān)鍵。本文采用的設(shè)計變量映射方法為多維度擴展方法，具體實現(xiàn)過程為：通過Kronecker乘積方法將粗網(wǎng)格模型的優(yōu)化設(shè)計變量以1對4(二維)或8(三維)的方式等值映射至細網(wǎng)格模型，進而實現(xiàn)細網(wǎng)格模型相應(yīng)設(shè)計變量的初始化，如圖8所示。通過這樣的分層迭代過程與代際映射繼承的方式，在保證優(yōu)化精度的前提下顯著提高優(yōu)化求解的效率。

圖7 多重網(wǎng)格等幾何拓撲優(yōu)化的算法流程圖

(a)二維多重網(wǎng)格模型映射示意圖

將多重網(wǎng)格的框架與式(15)所示的數(shù)學(xué)模型相結(jié)合，得到多重網(wǎng)格等幾何拓撲優(yōu)化模型：

(18)

式中，l為多重網(wǎng)格模型的網(wǎng)格層級，l=1，2，…，L；xl的元素取值區(qū)間為[0，1]。

2.2 基于Bézier提取算子的多重網(wǎng)格等幾何拓撲優(yōu)化模型

首先，基于Bézier提取算子對傳統(tǒng)等幾何拓撲優(yōu)化模型進行重構(gòu)；接著，將改進等幾何拓撲優(yōu)化模型與圖7所示的多重網(wǎng)格的框架結(jié)合，提出B樣條單元剛度矩陣統(tǒng)一表達的多重網(wǎng)格等幾何拓撲優(yōu)化模型。

(19)

(20)

式(19)中，Be，Ber是關(guān)于Bernstein基函數(shù)的偏導(dǎo)數(shù)矩陣。根據(jù)圖5可知，任一B樣條單元的基函數(shù)均可使用相同的Bernstein基函數(shù)空間進行等效表達，即各個單元的Be，Ber均相等。

根據(jù)上述特性及式(20)，可推導(dǎo)得到圖9所示的不同Bézier單元剛度矩陣之間的關(guān)系：

圖9 Bernstein空間下6×2網(wǎng)格的二維懸臂梁設(shè)計域

(21)

通過Bézier提取矩陣，圖9所示的B樣條單元剛度矩陣表達式如下：

(22)

根據(jù)等幾何分析的等參特性，除了幾何空間可使用Bézier提取算子進行映射之外，B樣條和Bézier控制點的位移向量也可使用Bézier提取算子實現(xiàn)二者之間的轉(zhuǎn)換。B樣條控制點位移u到Bézier控制點的位移ub的轉(zhuǎn)換方程如下：

(23)

基于以上的分析推導(dǎo)，對式(15)多重網(wǎng)格等幾何拓撲優(yōu)化模型進行重構(gòu)，可得

(24)

相對于傳統(tǒng)多重網(wǎng)格模型，本文提出的多重網(wǎng)格等幾何模型的單元數(shù)據(jù)存儲結(jié)構(gòu)優(yōu)化如圖10所示，其中，L、n分別表示多重網(wǎng)格模型的層數(shù)與等幾何分析域參數(shù)空間的維度，KE0為各層統(tǒng)一的標準單元剛度矩陣。傳統(tǒng)的多重網(wǎng)格模型每一層都需要進行所有B樣條單元剛度矩陣的計算和存儲，而本文多重網(wǎng)格模型采用標準Bézier單元剛度矩陣與相應(yīng)的Bézier提取算子實現(xiàn)任意層級的任一B樣條單元剛度矩陣的統(tǒng)一表達，大幅減少了多重網(wǎng)格等幾何拓撲優(yōu)化的預(yù)處理時耗與單元數(shù)據(jù)內(nèi)存的消耗。

圖10 傳統(tǒng)與本文多重網(wǎng)格等幾何拓撲優(yōu)化模型的單元剛度矩陣內(nèi)存結(jié)構(gòu)對比

2.3 靈敏度分析

根據(jù)式(24)的數(shù)學(xué)模型以及圖7所示的算法流程，目標函數(shù)相對于第l層的任意一個設(shè)計變量xe,l的靈敏度表達如下：

(25)

為避免數(shù)值不穩(wěn)定性現(xiàn)象，本文采用靈敏度過濾器修正目標函數(shù)的靈敏度公式：

(26)

(27)

式中，wa,l為第l層的第e個單元與第a個單元之間的權(quán)重因子；dist(e，a)表示第a個單元與第e個單元的質(zhì)心之間的歐氏距離；rmin為由用戶定義的過濾器半徑。

3 數(shù)值算例

本文算例中的程序均運行于Windows 10 系統(tǒng)的MATLAB R2021a，所用計算機的配置為Intel(R)Xeon(R)Gold 5102@2.20 GHz的CPU與64 GB內(nèi)存。圖11所示為激光雷達測試平臺的支撐結(jié)構(gòu)。對于激光雷達的主支撐結(jié)構(gòu)，可使用拓撲優(yōu)化方法進行優(yōu)化設(shè)計，主支撐厚度較小時可用二維拓撲優(yōu)化以簡化設(shè)計與制造流程，主支撐厚度較大時可使用三維拓撲優(yōu)化。二維和三維的激光雷達主支撐結(jié)構(gòu)設(shè)計模型如圖12所示。

圖11 激光雷達測試平臺

圖12a中二維主支撐結(jié)構(gòu)受到右側(cè)上端向下的外力負載，分別在該二維主支撐結(jié)構(gòu)的計算條件下應(yīng)用傳統(tǒng)多重網(wǎng)格模型[23]與基于Bézier提取的改進多重網(wǎng)格模型。其中，等幾何網(wǎng)格的初始尺寸為40×20，多重網(wǎng)格的最大層數(shù)為4。二維主支撐結(jié)構(gòu)的多重網(wǎng)格等幾何拓撲優(yōu)化結(jié)果和優(yōu)化收斂曲線分別如圖13和圖14所示，其中，PDOF為自由度，k為迭代次數(shù)。從圖中可知，傳統(tǒng)多重網(wǎng)格等幾何模型[23]和所提的改進多重網(wǎng)格等幾何模型具有相同的優(yōu)化收斂過程和優(yōu)化結(jié)果，驗證了提出的多重網(wǎng)格B樣條單元剛度矩陣統(tǒng)一映射法則在二維優(yōu)化問題中的正確性。

(a)二維算例模型

圖13 二維主支撐結(jié)構(gòu)的傳統(tǒng)[23]與改進多重網(wǎng)格模型優(yōu)化結(jié)果對比

圖14 二維主支撐結(jié)構(gòu)的傳統(tǒng)[23]與改進多重網(wǎng)格模型收斂曲線對比

此外，上述兩種多重網(wǎng)格等幾何拓撲優(yōu)化模型的各層剛度矩陣計算內(nèi)存消耗與預(yù)處理時間對比結(jié)果如圖15所示。由圖15可知：在各層剛度矩陣計算的內(nèi)存消耗上，改進模型較傳統(tǒng)模型減少99.76%～99.96%；在各層的剛度矩陣計算的預(yù)處理時間上，改進模型較傳統(tǒng)模型縮短95.89%～99.62%。因此，對于二維拓撲優(yōu)化問題，所提出的多重網(wǎng)格剛度映射方法可大幅減少多重網(wǎng)格等幾何拓撲優(yōu)化模型各層等幾何單元剛度矩陣的內(nèi)存消耗與預(yù)處理時耗。

圖15 二維主支撐結(jié)構(gòu)各層等幾何剛度矩陣所需的內(nèi)存空間與預(yù)處理時間對比

當多重網(wǎng)格的最大層數(shù)均為4且最大網(wǎng)格尺寸變化時，傳統(tǒng)模型與改進模型進行二維主支撐結(jié)構(gòu)優(yōu)化所需的單元剛度矩陣計算總內(nèi)存消耗與總預(yù)處理時間對比如圖16所示。在不同網(wǎng)格尺寸下，改進模型剛度矩陣數(shù)據(jù)的總體內(nèi)存消耗較傳統(tǒng)模型可減少98.95%～99.48%，總體預(yù)處理時間可縮短99.93%～99.96%。相比于傳統(tǒng)模型，改進模型在多種多重網(wǎng)格情況下進行二維優(yōu)化問題求解所需的等幾何單元剛度矩陣的總體內(nèi)存消耗與總預(yù)處理時耗均可大幅減少。

圖16 二維主支撐結(jié)構(gòu)在不同最大網(wǎng)格尺寸下剛度矩陣計算的總體內(nèi)存消耗與總預(yù)處理時間變化對比

三維主支撐結(jié)構(gòu)算例模型如圖12b所示，結(jié)構(gòu)的右上部分受到負載的作用。其中，等幾何網(wǎng)格的初始尺寸為40×20×10，多重網(wǎng)格的最大層數(shù)為2。三維主支撐結(jié)構(gòu)等幾何拓撲優(yōu)化的優(yōu)化結(jié)果和優(yōu)化收斂曲線分別如圖17和圖18所示。從圖中可知，傳統(tǒng)模型和改進模型具有相同的優(yōu)化收斂過程和優(yōu)化結(jié)果，驗證了提出的多重網(wǎng)格B樣條單元剛度矩陣統(tǒng)一映射法則在三維優(yōu)化問題中的正確性。

圖17 三維主支撐結(jié)構(gòu)的傳統(tǒng)[23]與改進多重網(wǎng)格等幾何拓撲優(yōu)化模型優(yōu)化結(jié)果對比

圖18 三維主支撐結(jié)構(gòu)的傳統(tǒng)[23]與改進多重網(wǎng)格模型收斂曲線對比圖

上述兩種多重網(wǎng)格等幾何拓撲優(yōu)化模型的各層剛度矩陣計算內(nèi)存消耗與預(yù)處理時間對比結(jié)果如圖19所示。從圖19中可知：在各層的剛度矩陣計算的內(nèi)存消耗上，改進模型較傳統(tǒng)模型減少99.97%～99.99%；在各層的剛度矩陣計算的預(yù)處理時間上，改進模型較傳統(tǒng)模型縮短94.96%～94.99%。因此，對于三維拓撲優(yōu)化問題，本文提出的多重網(wǎng)格剛度映射方法可大幅減少多重網(wǎng)格等幾何拓撲優(yōu)化模型的各層等幾何單元剛度矩陣的內(nèi)存消耗與預(yù)處理時耗。

圖19 三維主支撐結(jié)構(gòu)各層等幾何剛度矩陣所需的內(nèi)存空間與預(yù)處理時間對比

當多重網(wǎng)格模型的最大層數(shù)為2且最大網(wǎng)格尺寸變化時，傳統(tǒng)模型與改進模型進行三維主支撐結(jié)構(gòu)優(yōu)化所需的單元剛度矩陣總內(nèi)存消耗與總預(yù)處理時間對比結(jié)果如圖20所示。在不同網(wǎng)格尺寸下，改進模型單元剛度所需的總內(nèi)存消耗較傳統(tǒng)模型減少99.98%～99.99%，所需的總體預(yù)處理時間縮短94.97%～95.52%。因此，本文改進模型在不同的多重網(wǎng)格情況下進行三維優(yōu)化問題求解所需的等幾何單元剛度矩陣的總體內(nèi)存消耗與總預(yù)處理時間上都較傳統(tǒng)模型有明顯優(yōu)勢。

圖20 三維主支撐結(jié)構(gòu)在不同網(wǎng)格尺寸下剛度矩陣計算的總體內(nèi)存消耗與總預(yù)處理時間變化對比

綜合優(yōu)化模型設(shè)計分析與實驗結(jié)果分析，本文提出的改進多重網(wǎng)格模型具有正確性以及內(nèi)存和預(yù)處理時間方面的高效性，具體表現(xiàn)為：本文方案在優(yōu)化過程與優(yōu)化結(jié)果上與傳統(tǒng)方案保持一致，證明了本文方案是傳統(tǒng)方案正確有效的替代方案；二維與三維模型的拓撲優(yōu)化剛度計算中，改進模型的各層內(nèi)存消耗與各層預(yù)處理時間都較傳統(tǒng)模型具有優(yōu)勢；二維與三維模型的拓撲優(yōu)化剛度計算中，改進模型的總內(nèi)存消耗與總預(yù)處理時間都較傳統(tǒng)模型具有優(yōu)勢。

4 結(jié)語

本文提出一種基于Bézier單元剛度映射的改進多重網(wǎng)格等幾何拓撲優(yōu)化模型，解決了傳統(tǒng)的多重網(wǎng)格等幾何拓撲優(yōu)化模型剛度矩陣計算過程存在的單元數(shù)據(jù)存儲內(nèi)存消耗過大和預(yù)處理時間過長的問題。通過Bézier提取算子將B樣條基函數(shù)等效表達為具有空間保持特點的Bernstein基函數(shù)，建立了標準Bézier單元剛度矩陣與各層Bézier提取算子進行各層級B樣條單元剛度矩陣的統(tǒng)一表達，實現(xiàn)了多重網(wǎng)格等幾何拓撲優(yōu)化內(nèi)存結(jié)構(gòu)的優(yōu)化與B樣條單元數(shù)據(jù)預(yù)處理過程的簡化。在二維和三維拓撲優(yōu)化問題中，本文提出的模型在數(shù)據(jù)存儲負擔(dān)與預(yù)處理效率上均得到大幅度的改善，進而驗證了本文提出的改進多重網(wǎng)格等幾何拓撲優(yōu)化模型的有效性。