基于等幾何分析的板殼結構形狀-拓撲協同優化

郝 鵬，張越一，金靈智，馮少軍，王 博

(大連理工大學工程力學系工業裝備結構分析優化與CAE軟件全國重點實驗室，大連 116023)

0 引言

板殼結構在航空航天領域已經獲得了廣泛關注，其輕質化是提升運載火箭結構系數的重要途徑[1]。在以往研究中，結構輕量化設計一直是一個熱門主題，國內外學者進行了大量探索，主要采用有限元方法作為分析手段。然而，對于復雜板殼結構，采用傳統有限元方法進行優化是低效的。Hughes等[2]提出了一種基于樣條函數的等幾何分析(IGA)方法，采用非均勻有理B樣條(NURBS)精確靈活描述幾何模型，實現與分析模型之間的無縫集成，分析效率更高。

結構優化是輕量化設計的主要技術手段，包括尺寸優化、形狀優化和拓撲優化。本文主要探究形狀優化與拓撲優化。在形狀優化早期，選取有限元離散化邊界節點的坐標作為設計變量[3]。但這種方法面臨著一個難題，即不可避免地會出現邊界不規則與非光滑現象。等幾何方法具有高精度、精確幾何控制、計算機輔助設計(CAD)與計算機輔助工程(CAE)無縫結合等優勢。因此，相關學者將其應用于結構形狀優化。關鍵步驟是利用控制點坐標作為設計變量直接描述形狀，平滑邊界的同時也使分析模型準確代表了結構的幾何模型。Qian[4]提出了一種計算NURBS控制點位置和權重的全解析靈敏度分析方法，將控制點和權系數同時作為設計變量，大幅提升了設計空間。面向加筋薄壁結構的形狀優化設計，Hao等[5]提出了基于等幾何的形狀優化與筋條布局同步優化框架，該方法具有一定的準確性、靈活性與穩健性。

拓撲優化方法能夠在規定的設計區域內使定量材料達到最佳分布形式。但目前大多數拓撲優化方法中的數值分析主要采用有限元法，具有以下限制[6]：1)有限元網格無法準確描述結構幾何，降低了數值精度。2)低階(C0)連續性影響了優化結果的準確性。3)高質量有限元網格實現起來較為困難。為了突破上述局限性，Feng等[7]提出了一種有效的B樣條參數化方法，用于殼體結構的筋條布局優化，采用有限元法和固殼耦合法進行分析，得到了清晰布局。隨著等幾何的發展，有學者提出了基于IGA的拓撲優化方法，Kang等[8]通過使用NURBS修剪曲線代表孔的方式進行殼結構的等幾何拓撲優化工作，得到了光滑連續的材料布局。Gupta等[9]提出了一種新的使用PHT(Polynomial splines over Hierarchical T-meshes)樣條的自適應等幾何拓撲優化方法。如上所述，等幾何思想已經在很大程度上用于結構優化。

由以往研究可知，形狀優化無法改變拓撲構型，而拓撲優化的初始設計域固定，因此兩種方法均存在一定的限制。為此，相關學者相繼提出了不同的形狀與拓撲組合優化形式。Ansola等[10]提出了一種形狀與拓撲的集成優化方法，執行了連續的兩步程序，形狀與拓撲交替進行。Hassani等[11]基于有限元分析對殼體結構同時進行形狀和拓撲優化。Zhu等[12]提出了一種耦合形狀與拓撲優化(CSTO)技術用于支撐結構的設計。Jiang等[13]提出了一種基于移動變形構件(MMC)方法的顯示形狀與拓撲同時優化方法，最終優化結果可實現與CAD系統的轉化。上述研究結果表明，采用一些組合優化形式更能使結構獲得良好性能。

基于以上工作，本文提出了一種基于等幾何分析的形狀與拓撲協同優化方法。采用NURBS精確描述幾何模型，提供高平滑度并實現對形狀的靈活控制。將控制點坐標定義為形狀變量，控制點密度定義為拓撲變量。由于結構形狀的改變，需要更加精細的網格進行描述，因此通過等幾何分析可避免有限元中重新劃分網格的復雜過程。同時，采用3場SIMP方法(固體各向同性材料密度懲罰模型)獲取清晰邊界并消除中間密度。與經典拓撲優化結果對比，新方法不僅保證了計算效率與計算精度，而且得到了更高性能的優化結構。

1 板殼結構的等幾何分析

1.1 NURBS的基本原理

本節主要介紹以B樣條為基礎的NURBS曲面[14]。基于B樣條基函數可以得到NURBS基函數，NURBS基函數的權重和非均勻節點矢量使其在描述幾何形狀時更加靈活。對于NURBS曲面，雙變量的分段有理基函數可定義為

(1)

其中，i=1，2，…，n+p+1；j=1，2，…，m+q+1，n和m分別是ξ方向與η方向上基函數的個數，p和q分別為ξ方向與η方向上基函數的階次，Ni，p(ξ)和Nj，q(η)分別是p階和q階的單變量B樣條基函數，wi，j表示Ni，p(ξ)Nj，q(η)的相關權重。則NURBS曲面可定義為

(2)

式中，Pi，j是控制點。

1.2 基于NURBS的退化殼單元

在以往的工作中，已經將基于NURBS的退化殼單元用于復雜殼體的等幾何靜力與屈曲分析[15]。在退化殼單元中，任一點的總體坐標可近似地表示為結點坐標的插值形式，即

(3)

其中，t為殼體厚度，v3i為殼體局部法向矢量，ζ為高度方向上參數坐標。

根據殼體理論的基本假設，變形前中面的法線變形后仍保持為直線，且忽略其長度的變化，因此殼體內任一點的位移可由中面對應點沿總體坐標x，y，z方向的3個位移分量(ui，vi，wi)及該對應點的法線繞與它相垂直的兩個正交矢量的轉動(θxi，θyi)所確定。基于二階基函數的退化殼單元局部坐標和位移如圖1所示。

圖1 基于二階基函數退化殼單元局部坐標系及控制變量

通過精確的方向矢量，位移矢量可表示為

(4)

ifi×v3i=0,

(5)

彈性矩陣的表達式如下

(6)

該矩陣取決于材料的楊氏模量E和泊松比ν，k是考慮剪應力沿厚度方向不均勻分布的影響而引入的修正系數。整理后可得到退化殼的總體剛度陣

K=∑Ke

=∑(∑BTDB|J1||J2|w1iw2iw3i)

(7)

其中，Β為局部坐標系下的應力應變矩陣，w1i，w2i，w3i分別表示3個方向上的權系數，J1，J2為雅可比矩陣，J1用于物理空間與參數空間的轉換，J2用于參數空間與母空間之間的轉換。J1和J2分別表示如下

(8)

(9)

有關退化殼單元剛度陣與雅可比矩陣的詳細推導過程可參考文獻[16]。

2 形狀-拓撲協同優化

2.1 等幾何形狀優化

體積約束下，最小應變能為目標函數的形狀優化問題可描述為

(10)

C為目標函數，U為結構總位移，F為載荷，V0為結構初始體積，f為規定體積約束。設計變量為控制點坐標xs，xsmin與xsmax分別為自定義的設計變量坐標變化上下限，用于規定設計域的變化邊界，因此，采用移動漸近線(MMA)算法更新形狀變量。

2.1.1 解析靈敏度

目標函數對設計變量的偏導數為

(11)

(12)

其中，由于矩陣B在不同條件下是不斷變化的，因此，?B/?xs的計算是一個繁瑣的過程。另外，可采用半解析法或有限差分法近似計算形狀變量靈敏度，但解析法能得到更準確的靈敏度信息，關于解析解的推導過程可參考文獻[16]。

2.1.2 多水平模型

采用多水平模型說明等幾何形狀優化下的靈敏度傳遞問題，如圖2所示，包括初始模型、分析模型與設計模型。初始模型具有較少控制點，通過細化增加控制點得到設計模型，進一步細化得到分析模型滿足數值分析精度。設計變量定義在設計模型中，但目標函數由分析模型求解。在IGA的h細化過程中，原控制點被消除并被新控制點取代，保持結構形狀不變的同時確保分析準確性。因此，可以實現靈敏度從設計模型到分析模型之間的傳遞。NURBS曲面的h細化過程可參考文獻[16]。

圖2 多水平模型的網格與控制點

(13)

wQi=αiwi+(1-αi)wi-1

(14)

新控制點的物理坐標可表示為

(15)

在h細化的過程中，只有控制點發生變化。因此，相對于目標函數的靈敏度傳遞被轉化為控制點坐標的靈敏度傳遞。相對于設計變量的新控制點靈敏度表示為

(16)

2.2 等幾何拓撲優化

在等幾何拓撲優化過程中，為獲得更加清晰且光滑的邊界，消除中間密度，采用3場SIMP方法，即在過濾密度的基礎上，通過投影方式實現過濾密度的清晰化。體積約束下，最小應變能為目標函數的拓撲優化問題可描述為

(17)

設計變量為控制點密度xc，構造增廣拉格朗日函數，柔順度可表示為

C=FTU+λT(KU-F)

(18)

(19)

定義伴隨方程FT+λTK=0，載荷與結構無關時?F/?xc=0，同時考慮剛度陣的對稱性Kλ=-F。由KU=F可知λ=-U，故靈敏度可表示為

(20)

密度過濾的表達式如下

(21)

其中，

Hei=max{rmin-Δ(e,i),0}

(22)

式中，rmin為過濾半徑，Δ(e，i)為兩控制點之間的距離。

為了消除過濾帶來的灰度單元，之后通過投影實現清晰的優化構型。本文采用正切型Heavi-side函數，表達式如下

(23)

其中，η為投影閾值，β控制著函數在η附近的陡峭程度。

采用優化準則法(OC)進行迭代求解，設計變量的迭代更新表達式為

(24)

2.3 等幾何形狀-拓撲協同優化

基于等幾何分析的形狀-拓撲協同優化流程圖如圖3所示。采用一種先等幾何形狀優化，后等幾何拓撲優化的協同優化方法，對求得的最佳結構形狀進行拓撲優化，求解新結構在體積約束下的最小應變能。

圖3 基于等幾何分析的形狀-拓撲協同優化流程圖

3 算例

在本章中，通過3個數值實例驗證所提出的方法。第1個是拓撲優化中的常用測試算例，即平面MBB梁問題。第2和第3個例子是具有不同類型載荷條件的曲面問題，以此來驗證所提出框架對復雜曲面的適用性。所有NURBS模型均采用3階，模型參數均是量綱為1，材料屬性如下：彈性模量E=1，泊松比ν=0.3。殼體結構的厚度為1，過濾半徑為2，Heaviside函數中β=2，η=0.5。目標函數為最小應變能，約束函數為使目標體積小于或等于初始體積與體分比的乘積。為了驗證形狀-拓撲協同優化框架在結構輕量化設計方面的高效性，3個算例均與經典等幾何拓撲優化結果進行對比。

3.1 MBB梁

首先考慮MBB梁，邊界與載荷條件如圖4所示，邊長分別為600與100，模型的初始體積為60 000。形狀優化方面：考慮對稱性，選取初始模型的1/2為設計域，取其左側4個控制點的y向坐標作為形狀變量，其中上邊界的1，2號控制點的y向坐標改變量上下限是[-50，50](初始y向坐標為50，則坐標范圍是[0，100])；下邊界3，4號控制點的y向坐標改變量的范圍是[-230，-50](初始y向坐標為-50，則坐標范圍是[-100，-280])；形狀優化體積約束為整體體積的2.0倍。拓撲優化方面：分析模型的控制點為100×25，每個控制點的密度設置為設計變量，拓撲優化體積約束為整體體積的0.25倍。因此，形狀-拓撲協同優化的總體約束為初始體積的0.5倍。作為對比算例，在經典等幾何拓撲優化方法中，保持與形狀-拓撲協同優化相同的體積約束，總體積的約束值是初始體積的一半。

圖4 集中載荷作用下的MBB梁

所提方法與經典等幾何拓撲優化方法下參數域內的密度分布函數(Density Distribution Function，DDF)圖在圖5中給出，可初步觀察到結構構型的變化，所提方法下的結構桿件個數明顯減少。迭代曲線及優化過程中的一些中間設計如圖6所示。從形狀優化的迭代進程中可知，結構設計域向下演化，整體框架逐漸加寬，收斂速度快，穩定性強。結構的體積與應變能均在前4步內急劇變化，在第19步達到最佳形狀。在拓撲優化結果中，應變能在前10步內明顯減少，之后細微調整直至收斂，最終為桁架狀結構。對于類似結構，Rong等[17]使用不同方法獲得的優化結果如圖7所示，該方法通過增減子域的方式來重塑設計域，之后在自適應設計域中執行拓撲優化。由圖7可知，所提方法得到了相似的構型。兩種方法下使用的單元類型不一致，因此無法直接比較目標函數值。但值得注意的是，Rong等方法下的設計變量個數為98，優化周期較長，而提出方法的設計變量個數為4，優化迭代次數較少，只采用了極少的形狀變量，且收斂速度快，穩定性強。

(a)等幾何形狀-拓撲協同優化方法

(a)形狀優化結果C=0.112 8

(a)MBB梁載荷與邊界條件

經典等幾何拓撲優化方法下的迭代曲線與中間設計如圖8所示。與圖6協同優化方法下的優化結果對比發現，在體積基本相同的條件下，應變能相對減小了45.89%，得到了更高性能的結構，驗證了該方法的高效性。為了更清晰地展示結構形狀變化，所提方法下優化前后的控制點y向坐標見表1。可明顯地觀察到控制點變化的對稱性，下邊界大幅度下移，結構整體加寬。

表1 MBB梁優化前后控制點坐標

圖8 MBB梁的經典等幾何拓撲優化結果C=1.110 6

3.2 集中載荷下的自由曲面

中心點受集中載荷的自由曲面問題如圖9所示，其投影邊長為200，初始體積為40 720，4個端點固支。形狀優化方面：根據對稱性，選取初始模型的1/8作為設計域，取3個坐標方向為形狀變量，1號控制點的z向坐標，改變量范圍是[-50，50]，2號控制點的z向坐標與x向坐標(根據對稱性，正交方向上的形狀變量則為y向坐標)，改變量范圍均是[-30，30]，形狀優化體積約束為整體體積的1倍。拓撲優化方面：分析模型的控制點為60×60，每個控制點的密度設置為設計變量，拓撲優化體積約束為整體體積的3/10。因此，形狀-拓撲協同優化的總體約束為初始體積的3/10。作為對比算例，在經典等幾何拓撲優化方法中，保持與形狀-拓撲協同優化相同的體積約束，總體積的約束值是初始體積的3/10。

圖9 集中載荷作用下的自由曲面

如圖10所示，分別得到兩種優化方法下參數域內的DDF圖，可見二者結構形狀差別不大。形狀-拓撲協同優化結果如圖11所示，形狀優化方面：中心控制點z向坐標增大，結構具有更大翹曲度，體分比存在振蕩，但并未達到目標約束，這說明獲得最佳結構形狀時所需材料并不多。拓撲優化方面：應變能與體分比均迅速達到穩定狀態，結構中心開孔數量逐漸增多且清晰。

(a)等幾何形狀-拓撲協同優化方法

經典等幾何拓撲優化結果如圖12所示，對比兩組優化結果發現，形狀-拓撲協同優化方法的最終優化結構具有更大的彎曲程度。雖然形狀優化中結構體積并未達到目標約束，但應變能相對減少了71.80%，節省材料的同時結構性能進一步提高。形狀-拓撲協同優化方法下優化前后的控制點坐標見表2。由表2可知，中心控制點z向坐標已達到最大值，四邊中點坐標向內向上移動，結構翹曲度增大。

表2 自由曲面優化前后控制點坐標

圖12 自由曲面的經典等幾何拓撲優化結果C=769.908 4

3.3 均布載荷下的半圓柱殼

最后考慮半圓柱殼問題，模型尺寸與邊界條件在圖13中給出，其直徑為6，高度為12，初始體積為113，方向向下的均布載荷作用在高度方向邊界上，4個端點固支。形狀優化方面：考慮對稱性，選取初始模型的1/4作為設計域，取其6個控制點的法向坐標作為形狀變量，其坐標改變量的上下限均為[-2，2]；形狀優化體積約束為整體體積的1.3倍。拓撲優化方面：分析模型的控制點為60×84，每個控制點的密度設置為設計變量，拓撲優化體積約束為整體體積的3/10。因此，形狀-拓撲協同優化的總體約束為初始體積的39/100。作為對比算例，在經典等幾何拓撲優化方法中，保持與形狀-拓撲協同優化相同的體積約束，總體積的約束值是初始體積的39/100。

圖13 均布載荷作用下的半圓柱殼

所提出方法與經典等幾何拓撲優化方法下參數域內的DDF圖在圖14中給出，各自的迭代曲線與一些中間結構分別如圖15與圖16所示。從形狀-拓撲協同優化結果中可知，形狀優化方面：設計域改變的初始階段存在振蕩，之后迅速收斂，結構整體沿控制點法向方向外擴，由圓弧形變為“門”字形。拓撲優化方面：隨著迭代進程的推進，兩側弧形結構開孔個數增加。與經典等幾何拓撲優化結果對比可發現，結構沒有中間連接部分，且兩側弧形開孔結構的傾斜角度減小。同時，對比兩種方法下的最終目標函數值可知，在體積基本相同的條件下，等幾何形狀-拓撲協同優化方法下的應變能相對減小31.87%，提高了優化效率。

(a)等幾何形狀-拓撲協同優化方法

(a)形狀優化結果C=1 037.260 2

圖16 半圓柱殼的經典等幾何拓撲優化結果C=4 456.149 7

為了更清晰地觀察結構形狀變化，形狀-拓撲協同優化方法下優化前后的控制點坐標在表3中給出。根據對稱性可知，上下邊界中控制點的x向、y向坐標變化一致，中線控制點的變動較大，結構整體弧度減小，最終演化為“門”形結構。

表3 半圓柱殼優化前后控制點坐標

4 結束語

本文提出了基于等幾何分析的殼體結構形狀-拓撲協同優化框架，可以找到板殼結構的最佳結構形狀與最佳材料布局。在該過程中，等幾何形狀優化方面，推導了解析靈敏度公式，采用多水平模型說明了形狀變化過程中的靈敏度傳遞問題；等幾何拓撲優化方面，推導了靈敏度公式，考慮了過濾與投影以防止出現棋盤格與網格依賴性問題。最后，在MBB梁的經典算例中，與前人不同組合優化方法以及與經典等幾何拓撲優化方法相比，所提出方法可以實現具有改進性能的優化解決方案；之后在受集中載荷自由曲面與受均布載荷半圓柱殼的數值算例中，通過與經典等幾何拓撲優化方法下的結果進行對比，驗證了所提出方法的有效性與高效性。