□ 李 楊 □ 張衛(wèi)紅 □ 蔡守宇

西北工業(yè)大學 現(xiàn)代設計與集成制造技術教育部重點實驗室 西安 710072

結構形狀優(yōu)化旨在通過修改結構內(nèi)外邊界幾何形狀，改善結構特性（如降低應力集中、提高剛度、減輕結構質(zhì)量），主要包括三步：參數(shù)化幾何模型的建立；結構響應分析和靈敏度計算；優(yōu)化算法的實現(xiàn)。在優(yōu)化過程中，用于描述結構幾何形狀的參數(shù)作為設計變量，結構質(zhì)量與響應（如應力、柔順度、頻率等）作為目標函數(shù)和約束。

傳統(tǒng)的形狀優(yōu)化方法采用計算機輔助幾何設計（CAD）進行幾何建模，即利用Bezier、B樣條或非均勻有理B樣條（NURBS）參數(shù)化描述。由于CAD模型無法直接用于有限元分析（FEA），需要劃分有限元網(wǎng)格后 進 行 分 析 和 優(yōu) 化［1］， 幾 何 模 型 （CAD） 和 分 析 模 型（FEA）相互孤立，如圖 1（a）所示。

IGA方法將NURBS基函數(shù)不僅用于參數(shù)化幾何建模，而且用于力學場計算［2］，實現(xiàn)了幾何模型（CAD）、分析模型（IGA）和優(yōu)化模型的無縫連接，如圖1（b）所示。由于CAD模型直接用于等幾何分析，幾何形狀的精確表示使分析精度提高。目前，基于IGA的形狀優(yōu)化已被成功應用于平面線彈性問題［3,4,5］、三維線彈性問題［6］、梁結構［7］和振動腔問題［8］。 基于 T 樣條和Coons曲面的幾何參數(shù)化方法也被應用于等幾何形狀優(yōu)化中［9,10］。

▲圖1 優(yōu)化過程的比較示意圖

然而，優(yōu)化迭代過程中NURBS控制點的大幅變化，時常使相鄰控制點距離過近或過遠，甚至出現(xiàn)網(wǎng)格重疊畸形，需要在迭代過程中進行網(wǎng)格更新。文獻［3］采用內(nèi)部控制點的位置與邊界控制點的位置成固定比例的方法控制網(wǎng)格，這種方法雖然簡單但缺乏通用性。文獻［9］采用虛擬位移載荷法，把邊界控制點的變化量作為一個彈性輔助系統(tǒng)的位移邊界載荷處理，通過求解相應的線性方程組得到內(nèi)部控制點的變化量。方法較為通用，但是每次迭代之后相當于要額外進行一次有限元分析，效率較低。

本文借鑒文獻 ［11］空氣動力學中的網(wǎng)格調(diào)整策略，采用基于NURBS控制點間距的網(wǎng)格更新方法，進行IGA形狀優(yōu)化。首先簡要介紹IGA方法、相應的形狀優(yōu)化模型，然后建立網(wǎng)格更新方法，最后給出平面線彈性形狀優(yōu)化算例。

1 IGA等幾何分析方法

對于平面線彈性問題的IGA離散形式為：

式中：Ku=F為IGA離散的平衡方程；Ω為彈性體幾何域；Γ為面力邊界；K為總剛度矩陣；F為總載荷向量；u為節(jié)點位移向量；fb為體力載荷向量；fs為面力載荷向量；D 為彈性矩陣；為基函數(shù)矩陣，

上述基函數(shù)R可由雙變量NURBS基函數(shù)構造得到：

式中：ξ、η 為 NURBS 節(jié)點參數(shù)；p、q 為基函數(shù)階數(shù)；n、m 為基函數(shù)個數(shù)；ωi，j為權值。

類似于有限元分析，總剛度矩陣和總體載荷向量可以由單元剛度矩陣和單元載荷向量組裝而成。由此可見，等幾何分析得到的離散方程，形式上與有限元分析相同：NURBS控制點對應有限元分析中的節(jié)點，NURBS物理空間網(wǎng)格對應有限元分析中的網(wǎng)格，NURBS控制點上的位移對應有限元分析中的節(jié)點位移等。

2 基于等幾何分析的形狀優(yōu)化

2.1 形狀優(yōu)化模型

形狀優(yōu)化就是通過改變連續(xù)體結構內(nèi)外邊界形狀以改善結構特性（如減輕結構質(zhì)量、降低應力集中等）的一種優(yōu)化過程?；贗GA的形狀優(yōu)化問題一般可以描述為：

式中：x為與邊界NURBS控制點相關的設計變量；ND為設計變量個數(shù)；f和gj分別為目標函數(shù)和約束；NC為不等式約束個數(shù)。

2.2 網(wǎng)格更新方法

基于IGA的形狀優(yōu)化在迭代過程中，常會出現(xiàn)設計變量的大變形。如圖2（a），優(yōu)化迭代之前NURBS控制網(wǎng)格中控制點P在Y的正方向有大變形，那么迭代之后可能會出現(xiàn)網(wǎng)格重疊［如圖2（b）］。這樣變形之后的網(wǎng)格是畸形網(wǎng)格，無法繼續(xù)進行分析。為保證網(wǎng)格的協(xié)調(diào)性，采用一種基于控制點間距的網(wǎng)格更新方法。

▲圖2 NURBS控制網(wǎng)格大變形示意圖

考慮四分之一厚壁圓筒截面的NURBS控制點網(wǎng)格，如圖3所示。選擇實心圓代表的控制點作為設計變量，僅限制其在虛線（邊界法向）的方向上擾動。為方便說明，將圖中所示控制網(wǎng)格中的虛線上“一列”提取出來說明。 P0，P1，...，Pm為 m+1 個控制點的位置，其中，P0為邊界設計變量控制點（Variable Control Point），Pm為固定控制點（Fix Control Point），其余為連接控制點（Linked Control Point）。ΔP0，ΔP1，...，ΔPm為對應的控制點變形量，di=dist（Pi，Pi-1）（1≤i≤m）為相鄰兩控制點之間的距離。建立如下關系：

▲圖3 網(wǎng)格更新方法示意圖

故有：從而：

下面給出式（7）的3條性質(zhì)：

性質(zhì) 1：當 i=m 時，ΔPm=CmΔP0=0，這與 Pm為固定控制點相符。

性質(zhì) 2：若 ΔPi1=Ci1ΔP0，ΔPi2=Ci2ΔP0，則有當 i2＞i1 時，Ci2＜Ci1，所以 ΔPi1＞ΔPi2，說明:距離邊界可變控制點P0越近的點，變形量越大；距離P0越遠的點，變形量越小。

性質(zhì) 3：當 i2＞i1時，Pi2＞Pi1， 考察擾動后兩點Pi1+ΔPi1，Pi2+ΔPi2的位置關系，有：

又 ΔP0＜dsum，故（Pi2+ΔPi2）-（Pi1+ΔPi1）＞0 即 Pi2+ΔPi2＞Pi1+ΔPi1。也就是任意兩控制點在網(wǎng)格更新之后，都保持變形前的拓撲關系，不會出現(xiàn)網(wǎng)格重疊現(xiàn)象。

3 數(shù)值算例

基于上述的網(wǎng)格更新方法，下面給出等幾何形狀優(yōu)化的數(shù)值算例。優(yōu)化算法采用移動漸近線算法（MMA）［12］。 收斂準則為：

▲圖4 四分之一帶孔平板初始NURBS模型

▲圖5 優(yōu)化前后von-Mises應力云圖對比

式中：fi為第i次迭代之后的目標值；f0為目標初始值。

本文中控制點權值在優(yōu)化過程中保持不變。考慮平板應力集中問題，由于對稱性，僅考慮四分之一帶孔平板。問題初始模型定義如圖4所示，板兩邊受均布拉力，楊氏模量為210 GPa，泊松比為0.3，平板厚度為1 mm。在體分比99%的約束下，優(yōu)化目標為最大von-Mises應力最小。初始形狀下von-Mises應力云圖如圖5（a）所示，優(yōu)化后形狀和von-Mises應力云圖如圖5（b）所示。初始最大von-Mises應力為27.245 2 MPa，優(yōu)化后為15.744 2 MPa，應力集中明顯降低。

4 結論與展望

等幾何分析通過NURBS將CAD、FEA和結構形狀優(yōu)化結合起來，省去了繁瑣的模型轉(zhuǎn)換，擁有諸多優(yōu)點，有著極大的發(fā)展前景。此外，等幾何分析對于解決一些高階問題（如板殼問題）具有獨特的優(yōu)勢。本文采用一種基于NURBS控制點間距的網(wǎng)格更新方法，成功應用于等幾何形狀優(yōu)化。數(shù)值算例表明，此方法簡單有效，下一步將結合此方法進行三維的等幾何形狀優(yōu)化研究，并進一步嘗試等幾何形狀優(yōu)化在板殼問題中的應用。

［1］ V Braibant and C Fleury.Shape Optimal Design Using B-splines ［J］.Computer Methods in Applied Mechanics and Engineering,1984,44：247-267.

［2］ Hughes TJR,Cottrell JA,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.

［3］ Wall WA,Frenzel MA,Cyron C.Isogeometric Structural Shape Optimization ［J］.Computer Methods in Applied Mechanics and Engineering,2008,197（33-40）： 2976-2988.

［4］ Cho S,Ha S-H.Isogeometric Shape Design Optimization：Exact Geometry and Enhanced Sensitivity［J］.Structural and Multidisciplinary Optimization,2009,38（1）： 53-70.

［5］ Xiaoping Q.Full Analytical Sensitivities in NURBS Based Isogeometric Shape Optimization ［J］.Computer Methods in Applied Mechanics and Engineering,2010,199 （29-32）：2059-2071.

［6］ Hassani B,Tavakkoli SM,Moghadam NZ.Application of Isogeometric Analysis in Structural Shape Optimization ［J］.Scientia Iranica,2011,18（4）： 846-852.

［7］ Nagy AP,Abdalla MM,Gürdal Z.Isogeometric Sizing and Shape Optimisation of Beam Structures ［J］.Computer Methods in Applied Mechanics and Engineering,2010，199（17-20）： 1216-1230.

［8］ Manh ND,Evgrafov A,Gersborg AR,et al.Isogeometric Shape Optimization of Vibrating Membranes ［J］.Computer Methods in Applied Mechanics and Engineering,2011,200（13-16）： 1343-1353.

［9］ Ha S-H,Choi K,Cho S.Numerical Method for Sshape Optimization Using T-spline Based Isogeometric Method［J］.Structural and Multidisciplinary Optimization,2010,42 （3）：417-428.

［10］ Qian X,Sigmund O.Isogeometric Shape Optimization of Photonic Crystals via Coons Patches ［J］.Computer Methods in Applied Mechanics and Engineering， 2011，200（25-28）：2237-2255.

［11］ Pandya MJ,Baysal O.Gradient-based Aerodynamic Shape Optimization Using Alternating Direction Implicit Method［J］.Journal of Aircraft,1997,34（3）： 346-352.

［12］ K.Svanberg,The Method of Moving Asymptotes： A New Method for Structural Optimization［J］.International Journal of Numerical Methods in Engineering， 1987，24：359-373.