基于等幾何分析的板殼結(jié)構(gòu)形狀-拓撲協(xié)同優(yōu)化

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

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

0 引言

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

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

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

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

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

1 板殼結(jié)構(gòu)的等幾何分析

1.1 NURBS的基本原理

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

(1)

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

(2)

式中，Pi，j是控制點。

1.2 基于NURBS的退化殼單元

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

(3)

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

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

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

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

(4)

ifi×v3i=0,

(5)

彈性矩陣的表達式如下

(6)

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

K=∑Ke

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

(7)

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

(8)

(9)

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

2 形狀-拓撲協(xié)同優(yōu)化

2.1 等幾何形狀優(yōu)化

體積約束下，最小應變能為目標函數(shù)的形狀優(yōu)化問題可描述為

(10)

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

2.1.1 解析靈敏度

目標函數(shù)對設(shè)計變量的偏導數(shù)為

(11)

(12)

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

2.1.2 多水平模型

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

圖2 多水平模型的網(wǎng)格與控制點

(13)

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

(14)

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

(15)

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

(16)

2.2 等幾何拓撲優(yōu)化

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

(17)

設(shè)計變量為控制點密度xc，構(gòu)造增廣拉格朗日函數(shù)，柔順度可表示為

C=FTU+λT(KU-F)

(18)

(19)

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

(20)

密度過濾的表達式如下

(21)

其中，

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

(22)

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

為了消除過濾帶來的灰度單元，之后通過投影實現(xiàn)清晰的優(yōu)化構(gòu)型。本文采用正切型Heavi-side函數(shù)，表達式如下

(23)

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

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

(24)

2.3 等幾何形狀-拓撲協(xié)同優(yōu)化

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

圖3 基于等幾何分析的形狀-拓撲協(xié)同優(yōu)化流程圖

3 算例

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

3.1 MBB梁

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

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

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

(a)等幾何形狀-拓撲協(xié)同優(yōu)化方法

(a)形狀優(yōu)化結(jié)果C=0.112 8

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

經(jīng)典等幾何拓撲優(yōu)化方法下的迭代曲線與中間設(shè)計如圖8所示。與圖6協(xié)同優(yōu)化方法下的優(yōu)化結(jié)果對比發(fā)現(xiàn)，在體積基本相同的條件下，應變能相對減小了45.89%，得到了更高性能的結(jié)構(gòu)，驗證了該方法的高效性。為了更清晰地展示結(jié)構(gòu)形狀變化，所提方法下優(yōu)化前后的控制點y向坐標見表1?？擅黠@地觀察到控制點變化的對稱性，下邊界大幅度下移，結(jié)構(gòu)整體加寬。

表1 MBB梁優(yōu)化前后控制點坐標

圖8 MBB梁的經(jīng)典等幾何拓撲優(yōu)化結(jié)果C=1.110 6

3.2 集中載荷下的自由曲面

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

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

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

(a)等幾何形狀-拓撲協(xié)同優(yōu)化方法

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

表2 自由曲面優(yōu)化前后控制點坐標

圖12 自由曲面的經(jīng)典等幾何拓撲優(yōu)化結(jié)果C=769.908 4

3.3 均布載荷下的半圓柱殼

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

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

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

(a)等幾何形狀-拓撲協(xié)同優(yōu)化方法

(a)形狀優(yōu)化結(jié)果C=1 037.260 2

圖16 半圓柱殼的經(jīng)典等幾何拓撲優(yōu)化結(jié)果C=4 456.149 7

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

表3 半圓柱殼優(yōu)化前后控制點坐標

4 結(jié)束語

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