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

西北工業(yè)大學(xué) 現(xiàn)代設(shè)計(jì)與集成制造技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室 西安 710072

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

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

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

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

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

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

1 IGA等幾何分析方法

對(duì)于平面線彈性問(wèn)題的IGA離散形式為：

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

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

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

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

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

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

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

式中：x為與邊界NURBS控制點(diǎn)相關(guān)的設(shè)計(jì)變量；ND為設(shè)計(jì)變量個(gè)數(shù)；f和gj分別為目標(biāo)函數(shù)和約束；NC為不等式約束個(gè)數(shù)。

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

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

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

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

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

故有：從而：

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

性質(zhì) 1：當(dāng) i=m 時(shí)，ΔPm=CmΔP0=0，這與 Pm為固定控制點(diǎn)相符。

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

性質(zhì) 3：當(dāng) i2＞i1時(shí)，Pi2＞Pi1， 考察擾動(dòng)后兩點(diǎn)Pi1+ΔPi1，Pi2+ΔPi2的位置關(guān)系，有：

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

3 數(shù)值算例

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

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

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

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

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

4 結(jié)論與展望

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

［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.