陳 龍 張高朋

上海理工大學(xué)機械工程學(xué)院，上海，200093

?

基于靈敏度分析的三維結(jié)構(gòu)等幾何形狀優(yōu)化方法

陳 龍 張高朋

上海理工大學(xué)機械工程學(xué)院，上海，200093

以等幾何分析為數(shù)值分析方法，針對三維結(jié)構(gòu)的線彈性問題，以優(yōu)化邊界控制點坐標(biāo)為設(shè)計變量，通過Coons體插值方法得到插值控制點對優(yōu)化邊界控制點的坐標(biāo)導(dǎo)數(shù)，推導(dǎo)了單元剛度矩陣對控制點坐標(biāo)變量的形狀優(yōu)化靈敏度方程，并對應(yīng)力靈敏度和位移靈敏度方程進行了推導(dǎo)。數(shù)值算例證明該方法具有較高的優(yōu)化效率，優(yōu)化結(jié)果良好。

等幾何分析；三維結(jié)構(gòu)；形狀優(yōu)化；靈敏度

0 引言

在產(chǎn)品優(yōu)化設(shè)計階段，根據(jù)模擬分析階段得到的仿真結(jié)果，對初始的CAD 模型進行優(yōu)化設(shè)計，得到滿足需求的最優(yōu)CAD 模型，可進一步提高CAD 產(chǎn)品的性能[1-3]。根據(jù)現(xiàn)有方法優(yōu)化的模型對象往往和最初的CAD模型有較大出入。

優(yōu)化分為尺寸優(yōu)化、形狀優(yōu)化和拓?fù)鋬?yōu)化三個層次[4]。形狀優(yōu)化以結(jié)構(gòu)的邊界形狀為優(yōu)化對象，單元節(jié)點坐標(biāo)為設(shè)計變量，通過不斷改變邊界形狀達(dá)到提高產(chǎn)品性能的目的[5]。在產(chǎn)品拓?fù)浣Y(jié)構(gòu)的基礎(chǔ)上，進行部分邊界形狀的定量化描述，屬于產(chǎn)品的詳細(xì)設(shè)計階段。在基于傳統(tǒng)有限元方法的形狀優(yōu)化中，通常以有限單元的節(jié)點坐標(biāo)為設(shè)計變量，該方法存在以下弊端：①往往需要大量的設(shè)計變量描述邊界形狀；②優(yōu)化迭代中，各節(jié)點獨立變動，容易造成網(wǎng)格畸變與扭曲；③幾何模型與分析模型不統(tǒng)一，使得數(shù)據(jù)傳遞只能單向進行。靈敏度分析是通過求解優(yōu)化方程中的目標(biāo)函數(shù)和約束函數(shù)對設(shè)計變量的導(dǎo)數(shù)來進行優(yōu)化迭代更新的[6-7]。

等幾何分析方法基于有限元離散思想，將CAD系統(tǒng)采用的非均勻有理B樣條用于構(gòu)建NURBS單元，統(tǒng)一了幾何建模和有限元分析的整個過程[8]。將等幾何分析用于結(jié)構(gòu)形狀優(yōu)化，即以幾何模型控制點作為設(shè)計變量，解決了傳統(tǒng)形狀優(yōu)化難以處理光滑邊界、網(wǎng)格畸變等問題。

在優(yōu)化設(shè)計變量的選擇方面，QIAN[9]將控制點坐標(biāo)和權(quán)重同時作為設(shè)計變量進行形狀優(yōu)化，推導(dǎo)了敏度方程；SONG等[10]將形狀優(yōu)化過程分為控制點坐標(biāo)優(yōu)化和權(quán)重優(yōu)化兩步，設(shè)計初期將控制點坐標(biāo)作為設(shè)計變量，迭代一定程度后，將權(quán)重作為設(shè)計變量。網(wǎng)格更新規(guī)則方面的研究較少，主要依靠控制點坐標(biāo)相互之間距離的比例關(guān)系進行網(wǎng)格更新，容易引起單元形狀的不規(guī)則。

為了解決網(wǎng)格迭代更新的問題，本文以三維模型為例，通過Coons插值得到模型內(nèi)部控制點坐標(biāo)，籍此由邊界控制點坐標(biāo)對設(shè)計變量的導(dǎo)數(shù)得到內(nèi)部控制點坐標(biāo)對設(shè)計變量的導(dǎo)數(shù)，維持了整體網(wǎng)格的規(guī)整性，最后通過數(shù)值算例驗證了算法的效率與正確性。

1 B樣條表達(dá)與等幾何分析方法簡介

1.1 B樣條及NURBS表達(dá)

定義節(jié)點向量ξ=(ξ1，ξ2，…，ξn+p+1)，則B樣條曲線可以定義

(1)

其中，Pi為第i個控制點的坐標(biāo)；Ni，p為第i個p階基函數(shù)?；瘮?shù)可由Cox-de Boor遞推公式得到，p=0時，

(2)

p>0時，

(3)

類似地，B樣條曲面可以定義為

(4)

B樣條實體定義為

(5)

NURBS樣條為非均勻有理B樣條，由B樣條發(fā)展而來。本文優(yōu)化對象為三維實體，在此僅給出NURBS實體表達(dá)式：

(6)

(7)

式中，p、q、r為NURBS實體三個方向上的節(jié)點向量的階次。

可以看出，若所有控制點處的權(quán)重為1，則NURBS退化為B樣條。

以單個單元為研究對象，設(shè)一單元的基函數(shù)向量序列為

R=(R1，R2，…，RL)

(8)

其中，第s(s=1，2，…，L；L為單個單元所含控制點數(shù)目)個控制點處基函數(shù)表達(dá)為

Rs=MiNjQk

s=(k-1)mn+(j-1)n+i

i=1，2，…，nj=1，2，…，mk=1，2，…，l

則有該單元幾何場的插值形式：

[x(ξ,η,ζ)y(ξ,η,ζ)z(ξ,η,ζ)]T

(9)

(10)

1.2 三維線彈性問題的等幾何分析

本文只針對線彈性問題展開研究。應(yīng)用最小勢能原理，得到

Kq=R

(11)

單元e的單元剛度矩陣為

Ke=∫VeBTDBdV

(12)

三維問題中，單元剛度矩陣Ke可以寫為

Ke=∫VeBTDBdV

(13)

式中，B為mnl×6階矩陣；D為6階方陣。

等幾何分析與傳統(tǒng)有限元有幾點區(qū)別。首先，基函數(shù)表達(dá)式不同，傳統(tǒng)有限元分析多用拉格朗日基函數(shù)來離散表達(dá)幾何模型和位移場，而等幾何分析采用B樣條作為基函數(shù)。其次，等參變換的不同，等幾何分析采用等參單元思想，但需要進行物理坐標(biāo)、參數(shù)坐標(biāo)和母單元坐標(biāo)之間的映射，而傳統(tǒng)有限元的等參變換僅涉及到物理坐標(biāo)和參數(shù)坐標(biāo)的變換。等幾何分析中，每對參數(shù)變換對應(yīng)一個雅可比變換，其中，參數(shù)坐標(biāo)到母單元坐標(biāo)變換所對應(yīng)的雅可比變換行列式為

|J2|=(ξi+1-ξi)(ηj+1-ηj)(ζk+1-ζk)/8

(14)

物理單元坐標(biāo)到參數(shù)空間坐標(biāo)變換對應(yīng)的雅可比變換矩陣為

(15)

那么母單元到物理單元映射的雅可比行列式為

|J|=|J1||J2|

(16)

為了描述方便，定義兩個矩陣

應(yīng)變矩陣為(矩陣中的零元素均省略)

(17)

則應(yīng)變矩陣B可由Г矩陣變換元素位置得到。

將式(9)兩邊同時對 (ξ，η，ζ)求導(dǎo)，等號左邊為

(18)

等號右邊為

(19)

由式(18)、式(19)可得

(20)

Г、M與J1之間有下式關(guān)系：

(21)

通過式(21)，可求得應(yīng)變矩陣B。

采用高斯積分法對式(13)進行求解，則單元剛度矩陣表達(dá)為

(22)

式中，nel為單元總個數(shù)；mi、mj、mk分別為i、j、k方向上的高斯積分點的個數(shù)；wijk為高斯點處權(quán)重值。

通過式(22)求得單元剛度矩陣后，通過組裝即可得到總體剛度矩陣。

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

2.1 優(yōu)化方程描述

結(jié)構(gòu)優(yōu)化的目標(biāo)一般為結(jié)構(gòu)的體積或面積最小、結(jié)構(gòu)的剛度最大或柔度最小，約束一般是結(jié)構(gòu)等效應(yīng)力小于給定值或最大位移小于給定值。本文優(yōu)化問題定義為，在一定的體積約束下使得結(jié)構(gòu)柔度達(dá)到最小。優(yōu)化方程如下：

(23)

其中，C為優(yōu)化目標(biāo)柔度；[αi]為優(yōu)化設(shè)計變量序列，在文中為優(yōu)化形狀邊界控制點坐標(biāo)；V為優(yōu)化迭代中模型體積；V*為模型體積上限值；[αi]min、[αi]max分別為設(shè)計變量的下限值和上限值，Ku=f為離散控制方程；K為剛度矩陣；u為位移列陣；f為載荷列陣。

2.2 優(yōu)化流程

對于形狀優(yōu)化，首先確定需要進行優(yōu)化的模型邊界曲線或邊界曲面，將位于其上的控制點坐標(biāo)以參數(shù)表示，該參數(shù)即為優(yōu)化設(shè)計變量，通過Coons插值方法生成模型控制點參數(shù)化坐標(biāo)矩陣。針對該矩陣，逐次求得其對各設(shè)計變量的導(dǎo)數(shù)矩陣。通過等幾何分析程序求得目標(biāo)函數(shù)值及約束函數(shù)值和目標(biāo)函數(shù)及約束函數(shù)對設(shè)計變量的靈敏度。將上述數(shù)值代入優(yōu)化算法進行設(shè)計變量的更新，從而得到新的設(shè)計變量值，然后進行優(yōu)化終止條件判斷，滿足條件則終止優(yōu)化，否則以新的設(shè)計變量進行等幾何分析，再次得到更新的設(shè)計變量，直至得到最優(yōu)的設(shè)計變量，整個流程如圖1所示。

圖1 基于等幾何分析的形狀優(yōu)化流程圖Fig.1 Flow chart of shape optimization based on isogeometric analysis

優(yōu)化算法是影響形狀優(yōu)化質(zhì)量和效率的關(guān)鍵，本文選用移動漸近線法(method of moving asymptotes，MMA )。MMA需要輸入優(yōu)化問題的目標(biāo)函數(shù)和約束方程對設(shè)計變量的導(dǎo)數(shù)，因此需要通過等幾何分析程序求得該導(dǎo)數(shù)，即敏度信息。每次優(yōu)化迭代中，MMA通過等幾何分析計算得到的敏度矩陣更新設(shè)計變量值。然后根據(jù)迭代終止條件判斷優(yōu)化結(jié)果是否滿足收斂，若不收斂，則繼續(xù)優(yōu)化，直至得到最優(yōu)設(shè)計變量值。

2.3 Coons插值法獲得控制點坐標(biāo)

形狀優(yōu)化的設(shè)計變量一般為邊界曲面上的控制點坐標(biāo)，為了得到等幾何分析所需的控制點網(wǎng)格，需要在邊界曲面內(nèi)部進行插值。如圖2所示，已知實體邊界面上的控制點，通過Coons體插值，得到實體內(nèi)部所有控制點坐標(biāo)。設(shè)某邊界面控制點為Pi，j，k(i=1，2，…，l；j=1，2，…，m；k=1，2，…，n)，設(shè)l為u方向控制點數(shù)目，m為v方向控制點數(shù)目，n為w方向上的控制點數(shù)目。令

u0=(i-1)/(m-1)u1=1-u0

v0=(j-1)/(n-1)v1=1-v0

w0=(k-1)/(n-1)w1=1-w0

圖2 孔斯體插值方法Fig.2 Coons body interpolation method

那么內(nèi)部控制點可由邊界控制點通過下式求得：

Pi,j,k=u1P0,j,k+u0Pl,j,k+v1Pi,0,k+v0Pi,m,k+

(24)

求得邊界控制點坐標(biāo)對設(shè)計變量的導(dǎo)數(shù)后，基于式(24)可以得到模型所有控制的點坐標(biāo)對設(shè)計變量的靈敏度：

(25)

2.4 單元剛度矩陣對設(shè)計變量的靈敏度方程推導(dǎo)

靈敏度分析首先要求得離散控制方程(式(22))對設(shè)計變量的導(dǎo)數(shù)。整體剛度矩陣K是由單元剛度矩陣組裝而來的，因此轉(zhuǎn)化為求單元剛度矩陣Ke對設(shè)計變量導(dǎo)數(shù)問題。對設(shè)計變量求導(dǎo)得

(26)

根據(jù)式(26)，需要計算B及|J|對設(shè)計變量的導(dǎo)數(shù)，首先計算B對αi的導(dǎo)數(shù)。

式(20)、式(21)兩邊同時對設(shè)計變量求導(dǎo)有

(27)

(28)

將式(28)代入式(27)得

(29)

(30)

通過式(30)求得?Γ/?αi后，將?Γ/?αi矩陣中元素進行互換位置即可得到?B/?αi。接著需要計算

(31)

本文模型所加載荷不受設(shè)計變量變化的影響，故單元載荷列陣fe對設(shè)計變量導(dǎo)數(shù)為0，即?fe/?αi=0。將式(29)、式(31)代入式(26)即可求得單元剛度矩陣對設(shè)計變量的導(dǎo)數(shù)。

求得單元剛度矩陣對設(shè)計變量靈敏度后，通過單元剛度矩陣組裝成總體剛度矩陣的方法，得到總體剛度矩陣對各個設(shè)計變量的靈敏度矩陣，進而求得優(yōu)化方程中目標(biāo)和約束方程對設(shè)計變量的靈敏度。

2.5 目標(biāo)函數(shù)和約束函數(shù)對設(shè)計變量的靈敏度方程推導(dǎo)

在已知優(yōu)化方程(式(23))的情況下，可以對下列靈敏度方程進行求解。

(1)結(jié)構(gòu)柔度對設(shè)計變量的靈敏度：

(32)

(2)體積對設(shè)計變量的靈敏度：

(33)

(3)指定點位移對設(shè)計變量的靈敏度：

(34)

其中，行向量e除第i個元素為1外，其余元素均為0。那么通過式(34)即可得到位移向量中第i個位移值對設(shè)計變量的導(dǎo)數(shù)。

(4)指定點應(yīng)力分量對設(shè)計變量的靈敏度：

(35)

(36)

其中，λ為Kλ=(DBe)T的解。

(5)等效應(yīng)力對設(shè)計變量敏度：

(37)

令

則有

(38)

等效應(yīng)力對設(shè)計變量的導(dǎo)數(shù)為

(39)

將式(36)中應(yīng)力分量對設(shè)計變量的敏度代入式(39)即可求得等效應(yīng)力對設(shè)計變量的敏度。

3 數(shù)值算例分析

3.1 變截面懸臂梁形狀優(yōu)化

首先為了表述方便，將模型控制點通過集合的形式進行描述，將坐標(biāo)分量符合條件L0≤L≤L1且W0≤W≤W1且H0≤H≤H1的控制點P(L，W，H)定義為{P(L，W，H)|L∈[L0，L1]，W∈[W0，W1] ，H∈[H0，H1]}，其中，L0、L1、W0、W1、H0、H1均為常數(shù)值。

圖3 模型初始形狀及邊界條件(懸臂梁)Fig.3 Initial shape of model and boundary conditions (cantilever beam)

如圖3所示，模型長度方向的節(jié)點向量u=(0，0，0，0，0.1，0.2，0.4，0.5，0.6，0.7，0.8，0.9，1，1，1，1) ，寬度方向節(jié)點向量v、高度方向節(jié)點向量w均為(0，0，0，1，1，1)?？刂泣c數(shù)目u×v×w=11×3×3。材料的彈性模量E=1.5 MPa，泊松比μ=0.3。載荷與邊界條件：模型左端面上控制點{P(L，W，H)|L=0，W∈[0，2]，H∈[0，2]}完全固定，右端面上控制點{P(L，W，H)|L=10，W∈[0，2]，H∈[0，2]}均受到沿H軸負(fù)方向的載荷F=-2 mN。優(yōu)化設(shè)計變量：取下底面上控制點{P(L，W，H)|L∈[0，10]，W∈[0，2]，H=0}的H坐標(biāo)為設(shè)計變量。設(shè)置對稱約束條件：控制點{P(L，W，H)|L∈[0，10]，W∈[0，2]，H=0}與控制點{P(L，W，H)|L∈[0，10]，W∈[0，2]，H=2}關(guān)于面H=1對稱。優(yōu)化目標(biāo)為柔度最小化，約束為體積V≤V0=40，優(yōu)化的最終形狀如圖4所示。圖5顯示了優(yōu)化迭代過程中，柔度和模型體積的變化。通過表1可知，在體積不增加的情況下，柔度降低了40%。

圖4 最終優(yōu)化結(jié)果(懸臂梁)Fig.4 Final optimization results(cantilever beam)

(a)柔度

(b)體積圖5 迭代曲線(懸臂梁)Fig.5 Iterations(cantilever beam)

柔度(J)體積(mm3)優(yōu)化前0．052840優(yōu)化后0．032240

3.2 兩端固支梁形狀優(yōu)化

該例沿用例3.1中關(guān)于控制點的表述方法。

圖6 模型初始形狀及邊界條件(固支梁)Fig.6 Initial shape of model and boundary conditions (built-in beam)

如圖6所示，模型節(jié)點向量u=(0，0，0，0，0.1，0.2，0.4，0.5，0.6，0.7，0.8，0.9，1，1，1，1) ，v=(0，0，0，1，1，1)，w=(0，0，0，0.33，0.66，1，1，1)?？刂泣c數(shù)目u×v×w=11×3×5。材料的彈性模量E=1.5 MPa，泊松比μ=0.3。載荷與邊界條件：左端面上控制點{P(L，W，H)|L=0，W∈[0，2]，H∈[0，4]}和右端面上控制點{P(L，W，H)|L=15，W∈[0，2]，H∈[0，2]}完全固定，上底面上控制點{P(L，W，H)|L∈[0，15]，W=1，H=4}受到沿H軸負(fù)方向的載荷F=-100 mN。優(yōu)化設(shè)計變量：控制點{P(L，W，H)|L∈[0，15]，W=0，H=[0，4]}的W坐標(biāo)為設(shè)計變量。設(shè)置優(yōu)化對稱約束，控制點{P(L，W，H)|L∈[0，15]，W=0，H=[0，4]}與控制點{P(L，W，H)|L∈[0，15]，W=2，H=[0，4]}關(guān)于面W=1對稱。優(yōu)化目標(biāo)為柔度最小化，約束為體積V≤V0=120，優(yōu)化的最終形狀如圖7所示。圖8顯示了優(yōu)化迭代過程中，柔度和模型體積的變化。通過表2可以得到，體積不增加的情況下，柔度降低了23%。

圖7 最終優(yōu)化結(jié)果(固支梁)Fig.7 Final optimization results(built-in beam)

4 結(jié)語

本文的形狀優(yōu)化分別以柔度和體積為目標(biāo)和約束，以單片NURBS三維模型為優(yōu)化對象，基于等幾何分析推導(dǎo)了形狀優(yōu)化的靈敏度方程，并通過具體實例證明了方法的可行性，結(jié)果表明優(yōu)化收斂速度較快，形狀邊界光滑。

將等幾何分析應(yīng)用于形狀優(yōu)化，統(tǒng)一了產(chǎn)品的設(shè)計、分析和優(yōu)化過程。NURBS基函數(shù)的應(yīng)用使得和CAD系統(tǒng)的交互更加便捷，可以得到無限光滑的模型邊界，省去了優(yōu)化迭代中復(fù)雜的網(wǎng)格重劃分操作。今后將根據(jù)結(jié)構(gòu)具體失效形式，選取局部應(yīng)力、最大等效應(yīng)力或局部位移為約束條件。在優(yōu)化模型拓?fù)鋸?fù)雜度方面，將以多片NURBS構(gòu)成的復(fù)雜拓?fù)浣Y(jié)構(gòu)為優(yōu)化對象。

(a)柔度

(b)體積圖8 迭代曲線(固支梁)Fig.8 Iterations(built-in beam)

柔度(J)體積(mm3)優(yōu)化前0．5886120優(yōu)化后0．4531120

當(dāng)然，Coons插值方法并不能針對所有模型生成適合等幾何分析的參數(shù)化模型，尤其是在模型比較復(fù)雜的時候，需要引入控制點優(yōu)化算法，而針對這些優(yōu)化的控制點，如何求解其靈敏度矩陣將是下一步研究的方向。

[1] QIAN X. Topology Optimization in B-spline Space[J]. Computer Methods in Applied Mechanics and Engineering, 2013, 265: 15-35.

[2] CHRISTENSEN P W, KLARBRING A. An Introduction to Structural Optimization[J]. Solid Mechanics & Its Applications, 2008, 153(1):121-146.

[3] WALL W A, FRENZEL M A, CYRON C. Isogeometric Structural Shape Optimization[J]. Computer Methods in Applied Mechanics & Engineering, 2008, 197(33/40): 2976-2988.

[4] HASSANI B, TAVAKKOLI S M, MOGHADAM N Z. Application of Isogeometric Analysis in Structural Shape Optimization[J]. Scientia Iranica, 2011, 18(4):846-852.

[5] NAGY A P, ABDALLA MM, GüRDAL Z. Isogeometric Sizing and Shape Optimization of Beam Structures[J]. Computer Methods in Applied Mechanics & Engineering, 2006, 199(17/20):1216-1230.

[6] 李初曄， 王衛(wèi)朝， 馬巖. 基于參數(shù)靈敏度的結(jié)構(gòu)性能優(yōu)化[J]. 中國機械工程， 2011, 22(4):397-402. LI Chuye, WANG Weichao, MA Yan. Structural Capability Optimization Based on Sensitivity of Parameters[J]. China Mechanical Engineering, 2011, 22(4): 397-402.

[7] 李蕾，張葆，李全超，等. 基于靈敏度數(shù)的薄板結(jié)構(gòu)加強筋布局優(yōu)化設(shè)計[J]. 中國機械工程，2016，27(9): 1143-1149. LI Lei, ZHANG Bao, LI Quanchao, et al. Stiffener Lay Out Optimization of Thin Plate Structures Based on Sensitivity Number[J]. China Mechanical Engineering, 2016, 27(9):1143-1149.

[8] NGUYEN V P, ANITESCU C, BORDAS S P A, et al. Isogeometric Analysis: an Overview and Computer Implementation Aspects[J]. Mathematics & Computers in Simulation, 2012, 117:89-116.

[9] QIAN X. Full Analytical Sensitivities in NURBS Based Isogeometric Shape Optimization[J]. Computer Methods in Applied Mechanics and Engineering, 2010，199(29/32): 2059-2071.

[10] SONG Y U，HUR J Y，YOUN S K. 2D Isogeometric Shape Optimization Considering both Control Point Positions and Weights as Design Variables[C]//10th World Congress on Structural and Multidisciplinary Optimization. Orlando, 2013:1-8.

(編輯 張 洋)

Three Dimensional Structural Shape Optimization Based on Sensitivity Analysis UsingIsogeometric Analysis

CHEN Long ZHANG Gaopeng

School of Mechanical Engineering,University of Shanghai for Science and Technology,Shanghai,200093

An isogeometric analysis method was used as numerical analysis method. As for three-dimensional structures of linear elastic problem, boundary control points that needed to be optimized were regarded as design variables, derivatives of internal control points as to design variables might be obtained by using Coons interpolation method. In shape optimization, sensitivity equation of element stiffness matrix as to control points coordinates was derived, and sensitivity equations of stresses and displacements were deduced. Numerical examples show that the method has high efficiency and good results.

isogeometric analysis; three-dimensional structure; shape optimization; sensitivity

2016-08-24

國家自然科學(xué)基金資助項目(51475309,61472111 )

TB121

10.3969/j.issn.1004-132X.2017.14.015

陳 龍，男，1978年生。上海理工大學(xué)機械工程學(xué)院副教授。主要研究方向為產(chǎn)品計算設(shè)計、圖像處理。發(fā)表論文30余篇。E-mail: wwwchenl@163.com。張高朋，男，1991年生。上海理工大學(xué)機械工程學(xué)院碩士研究生。