徐家琪，馬永其,2

(1.上海大學(xué) 上海市應(yīng)用數(shù)學(xué)與力學(xué)研究所，上海 200072；2.上海大學(xué) 力學(xué)系，上海 200444)

結(jié)構(gòu)拓?fù)鋬?yōu)化的源頭可追溯至20世紀(jì)初Michell等[1]提出的應(yīng)力約束下具有最小重量桁架的結(jié)構(gòu)優(yōu)化，它的出現(xiàn)使人們擺脫了原有依靠經(jīng)驗(yàn)獲得最優(yōu)結(jié)構(gòu)的思維模式，開(kāi)始依靠理論計(jì)算獲得經(jīng)濟(jì)且實(shí)用的結(jié)構(gòu)。連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化是在不知道結(jié)構(gòu)拓?fù)湫螤畹那疤嵯?，根?jù)已知邊界條件和荷載條件確定比較合理的結(jié)構(gòu)拓?fù)湫问剑c離散體結(jié)構(gòu)拓?fù)鋬?yōu)化相比，設(shè)計(jì)變量多，模型復(fù)雜。自1988年Bends?e等[2]提出了均勻化方法(Homogenization Method)開(kāi)啟連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化研究以來(lái)，針對(duì)連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的研究逐漸廣泛和深入。Bends?e等[3]在先前研究的基礎(chǔ)上提出了變密度法，Sui等[4]提出了獨(dú)立連續(xù)映射法(Independent Continuous Mapping，ICM)，Xie等[5]提出了漸進(jìn)結(jié)構(gòu)優(yōu)化法(Evolutionary Structure Optimization，ESO)等用來(lái)進(jìn)行實(shí)現(xiàn)連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的計(jì)算。

研究連續(xù)體的結(jié)構(gòu)拓?fù)鋬?yōu)化除了集中在靜力學(xué)領(lǐng)域以外，不少學(xué)者也進(jìn)行了動(dòng)力學(xué)領(lǐng)域的研究，主要包括兩個(gè)方面，一是以特征頻率作為目標(biāo)函數(shù)或是約束條件的結(jié)構(gòu)拓?fù)鋬?yōu)化。如Xie等[6]采用ESO方法以特征頻率作為約束條件進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化；邱海等[7]對(duì)具有頻率約束的板結(jié)構(gòu)進(jìn)行了拓?fù)鋬?yōu)化研究；Díaaz等[8]采用均勻化方法以結(jié)構(gòu)特征頻率最大化為目標(biāo)函數(shù)進(jìn)行了結(jié)構(gòu)拓?fù)鋬?yōu)化；Pedersen[9]采用變密度法以結(jié)構(gòu)特征頻率最大化為目標(biāo)函數(shù)進(jìn)行了結(jié)構(gòu)拓?fù)鋬?yōu)化。魏星等[10]基于自然單元法，研究了功能梯度板頻率優(yōu)化。問(wèn)題二是在動(dòng)載荷作用下，以結(jié)構(gòu)柔度為目標(biāo)函數(shù)的結(jié)構(gòu)拓?fù)鋬?yōu)化。如Ma等[11]基于均勻化方法研究了結(jié)構(gòu)受激勵(lì)載荷的動(dòng)態(tài)響應(yīng)問(wèn)題，Jog等[12]定義了結(jié)構(gòu)的動(dòng)力學(xué)柔度研究了帶有周期載荷的結(jié)構(gòu)拓?fù)鋬?yōu)化，徐斌等[13]采用ESO方法研究了頻率激勵(lì)載荷下薄板的結(jié)構(gòu)拓?fù)鋬?yōu)化。劉虎等[14]將結(jié)構(gòu)指定位置穩(wěn)態(tài)階段位移響應(yīng)幅值作為目標(biāo)函數(shù)，研究了簡(jiǎn)諧力激勵(lì)下結(jié)構(gòu)拓?fù)鋬?yōu)化問(wèn)題。

大多數(shù)的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化計(jì)算都是基于有限元方法進(jìn)行的，由于有限元方法本身對(duì)網(wǎng)格的依賴性，網(wǎng)格劃分、網(wǎng)格數(shù)量和形式都會(huì)對(duì)連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果產(chǎn)生極大影響，且會(huì)使優(yōu)化結(jié)果產(chǎn)生棋盤(pán)格現(xiàn)象[15]。因此避免棋盤(pán)格現(xiàn)象的產(chǎn)生，同時(shí)提高計(jì)算精度與計(jì)算效率是目前主要的研究方向[16-18]。

近年來(lái)，隨著無(wú)網(wǎng)格方法的興起，人們開(kāi)始采用無(wú)網(wǎng)格方法取代有限元方法以避免優(yōu)化過(guò)程中棋盤(pán)格現(xiàn)象的產(chǎn)生。目前，已經(jīng)發(fā)展的無(wú)網(wǎng)格方法有十幾種，主要包括光滑粒子法[19]，無(wú)單元Galerkin方法[20](EFG)，重構(gòu)核粒子法[21]，無(wú)網(wǎng)格點(diǎn)插值法[22](PIM)，自然單元法[23](NEM)等。Zhou等[24]采用了重構(gòu)核粒子法對(duì)連續(xù)體結(jié)構(gòu)進(jìn)行了靜力學(xué)拓?fù)鋬?yōu)化研究。Zheng等[25-26]采用無(wú)單元Galerkin方法以及徑向點(diǎn)插值法分別進(jìn)行了靜力學(xué)連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化研究。Wu等[27]將改進(jìn)的無(wú)單元Galerkin方法應(yīng)用于連續(xù)體拓?fù)鋬?yōu)化。在動(dòng)力學(xué)連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化研究方面，鄭娟等[28-29]采用無(wú)網(wǎng)格徑向插值點(diǎn)法、無(wú)單元Garlerkin方法，進(jìn)行了頻率激勵(lì)載荷下連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的研究，使優(yōu)化結(jié)果避免棋盤(pán)格現(xiàn)象的產(chǎn)生。但是提高計(jì)算精度與計(jì)算效率依然是人們努力研究的方向。大部分無(wú)網(wǎng)格方法相關(guān)的形函數(shù)不滿足Delta性質(zhì)，不能直接施加本證邊界條件，必須通過(guò)增加Lagrange乘子法、罰函數(shù)方法等附加計(jì)算，才能施加本證邊界條件，不可避免地增加了計(jì)算成本，影響計(jì)算效率。自然單元法是基于Voronoi圖和Delaunay三角形概念進(jìn)行插值的無(wú)網(wǎng)格方法，其形函數(shù)簡(jiǎn)單，計(jì)算效率高，邊界上滿足Delta插值性質(zhì)，可以直接施加本征邊界條件，避免增加附加計(jì)算成本。

本文利用自然單元法求解效率高的特點(diǎn)，將動(dòng)力學(xué)自然單元法應(yīng)用于頻率激勵(lì)載荷作用下的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化計(jì)算，以節(jié)點(diǎn)密度作為設(shè)計(jì)變量，結(jié)構(gòu)動(dòng)柔度最小作為目標(biāo)函數(shù)，施加體積約束，采用基于變密度法的各向同性固體微結(jié)構(gòu)懲罰模型(Solid Isotropic Microstructure with Penalization，SIMP)，建立拓?fù)鋬?yōu)化模型并采用優(yōu)化準(zhǔn)則法(OC)求解優(yōu)化模型，在避免棋盤(pán)格現(xiàn)象產(chǎn)生的同時(shí)提高計(jì)算的效率和精度，并通過(guò)算例說(shuō)明該方法的有效性和優(yōu)越性。

1 動(dòng)力學(xué)問(wèn)題自然單元法

1.1 Voronoi圖與Delaunay結(jié)構(gòu)

自然單元法(Natural Element method，NEM)是基于Voronoi圖及其對(duì)偶圖形Delaunay三角剖分來(lái)建立插值函數(shù)的。

在二維空間R2中分布著一系列點(diǎn)所組成的點(diǎn)集S={x1,x2,…,xn}。點(diǎn)xi的Voronoi晶胞定義為Pi。當(dāng)d(x,xi)代表空間中任意一點(diǎn)到xi的距離時(shí)，Pi可定義為

Pi={x∈R2|d(x,xi)

(1)

如圖1所示為包含7個(gè)點(diǎn)的Voronoi圖。將Voronoi圖中擁有公共邊的兩個(gè)晶胞所對(duì)應(yīng)的節(jié)點(diǎn)相連結(jié)，就可以構(gòu)成Delaunay三角剖分圖。與一般的三角剖分不同的是，Delaunay三角形具有空外接圓以及最大化最小角的特性。該特性為尋找插值點(diǎn)自然鄰點(diǎn)，建立插值函數(shù)提供了依據(jù)。繪制Delaunay三角形的外接圓，并找出外接圓的圓心，如圖2所示，此時(shí)Delaunay三角形的外接圓圓心即為Voronoi晶胞的頂點(diǎn)。

圖1 帶有7節(jié)點(diǎn)Vornonoi圖Fig.1 Voronoi diagram with seven nodes

圖2 Delaunay三角形外接圓Fig.2 Natural neighbor circumcircles

1.2 自然鄰點(diǎn)插值

如圖3所示，若點(diǎn)x為插值點(diǎn)，該點(diǎn)位于點(diǎn)1,6,7和點(diǎn)1,5,6組成的兩個(gè)三角形的外接圓中，則稱點(diǎn)1,5,6,7為插值點(diǎn)x的自然鄰點(diǎn)。插值點(diǎn)x的二階Voronoi晶胞如圖4所示。確定插值點(diǎn)的自然鄰點(diǎn)之后，即可構(gòu)造插值函數(shù)式(2)

(2)

式中:h(x)為插值點(diǎn)的物理量N為插值點(diǎn)的自然鄰點(diǎn)的個(gè)數(shù);hI為插值點(diǎn)的第I個(gè)自然鄰點(diǎn)的物理量;φI為自然鄰點(diǎn)I所對(duì)應(yīng)的形函數(shù)。形函數(shù)的構(gòu)成方法可以分為Sibson插值和non-Sibson插值。其中Sibson插值方法為

(3)

(4)

插值點(diǎn)x關(guān)于第1個(gè)自然鄰點(diǎn)的形函數(shù)則可表示為

(5)

圖3 插值點(diǎn)x的自然鄰點(diǎn)Fig.3 Natural neighbor nodes of x

圖4 插值點(diǎn)x的二階Voronoi晶胞Fig.4 Second order Voronoi diagram of x

1.3 動(dòng)力學(xué)問(wèn)題的平衡方程

二維彈性動(dòng)力學(xué)問(wèn)題的平衡方程為

(6)

(7)

(8)

(9)

(10)

(11)

式中:B為應(yīng)變矩陣;D平面應(yīng)力情況下的為彈性矩陣;ρm為質(zhì)量密度;Φ為插值點(diǎn)所對(duì)應(yīng)的形函數(shù)。令U(t)=UAeiwt,F(t)=FAeiwt，式(9)可表示為

KUAeiwt-ω2MUAeiwt=FAeiwt

(12)

同時(shí)約去eiwt，則

(K-ω2M)UA=FA

(13)

式中:ω為激勵(lì)頻率;UA為位移幅值;FA為頻率激勵(lì)載荷幅值。

2 頻率激勵(lì)載荷拓?fù)鋬?yōu)化

2.1 SIMP模型

SIMP模型的基本思想是將連續(xù)體假想為一種密度可變的材料，每個(gè)節(jié)點(diǎn)的相對(duì)密度可看作在0～1之間變化。對(duì)彈性矩陣及質(zhì)量矩陣進(jìn)行插值的SIMP模型為

D=ρPD0

(14)

M=ρQM0

(15)

式中:P，Q分別為針對(duì)剛度矩陣和質(zhì)量矩陣的懲罰因子。當(dāng)相對(duì)密度ρ=1時(shí)，表示插值點(diǎn)具有實(shí)體材料，當(dāng)相對(duì)密度ρ=0時(shí)，意味著該處沒(méi)有材料存在。

對(duì)于懲罰因子P來(lái)說(shuō)，上式矩陣D中包含材料楊氏模量E和泊松比μ。根據(jù)Hashin-Shtrikman上下限定理，材料的體積模量K和剪切模量G需滿足下列不等式

(16)

式中:K0、G0分別為初始的體積模量和剪切模量。推導(dǎo)該不等式，可得出懲罰因子P的下限，與材料的泊松比μ有關(guān)。

(17)

對(duì)于懲罰因子Q，若其值小于P，在相對(duì)密度低的區(qū)域，剛度矩陣與質(zhì)量矩陣的比值較小。為保證衰減的同步性，這里采用文獻(xiàn)[28]的做法，采用相等的P和Q。

2.2 拓?fù)鋬?yōu)化模型

拓?fù)鋬?yōu)化的目標(biāo)函數(shù)設(shè)置為結(jié)構(gòu)的動(dòng)柔度最小，設(shè)計(jì)變量為節(jié)點(diǎn)相對(duì)密度，施加體積約束。拓?fù)鋬?yōu)化的數(shù)值模型可表示為

(18)

(19)

(20)

3 靈敏度分析

求解優(yōu)化模型涉及的目標(biāo)函數(shù)以及約束的靈敏度，通過(guò)伴隨分析法進(jìn)行求解。

(21)

式中,A為任意實(shí)向量，此時(shí)，右邊第二項(xiàng)為0。對(duì)上式求導(dǎo)為

(22)

(23)

(24)

結(jié)構(gòu)的體積表示為式(25)，結(jié)構(gòu)體積的靈敏度表示為式(26)。

(25)

(26)

4 算 例

4.1 懸臂板

考慮如圖5所示的懸臂版，該板長(zhǎng)為10 m，寬為6 m。右側(cè)中點(diǎn)處加載一個(gè)豎直向下的頻率激勵(lì)載荷P=5eiwtkN，激勵(lì)載荷幅值為5 kN。結(jié)構(gòu)的彈性模量E為3×107Pa，泊松比μ為0.3，結(jié)構(gòu)的質(zhì)量密度為1 kg/m3，剛度矩陣以及質(zhì)量矩陣的懲罰因子同時(shí)設(shè)置為4。將結(jié)構(gòu)離散為25×15共375個(gè)節(jié)點(diǎn)，剖分出672個(gè)Delaunay三角形，如圖6所示。在每個(gè)三角形內(nèi)插入3個(gè)高斯積分點(diǎn)。設(shè)置體積約束為50%，收斂條件設(shè)置為設(shè)計(jì)變量的最大變化小于0.01。

圖5 懸臂板結(jié)構(gòu)Fig.5 The structure of cantilever plate

圖6 懸臂板Delaunay三角剖分Fig.6 Delaunay triangles of cantilever plate

當(dāng)頻率為從0 Hz逐漸增加到40 Hz時(shí)的拓?fù)鋬?yōu)化圖形如圖7所示，頻率為50 Hz的拓?fù)鋬?yōu)化圖形如圖9所示。當(dāng)ω=0 Hz時(shí)，相當(dāng)于靜力學(xué)拓?fù)鋬?yōu)化。激勵(lì)頻率不同，所獲得的拓?fù)鋬?yōu)化圖形相似。將優(yōu)化百分比定義為柔度減少量占原有柔度的比值，當(dāng)激勵(lì)頻率增大時(shí)，優(yōu)化前動(dòng)柔度增大，優(yōu)化后的動(dòng)柔度也隨之增大，目標(biāo)函數(shù)變化如表1所示。當(dāng)激勵(lì)頻率越小時(shí)，優(yōu)化百分比越大，說(shuō)明獲得的優(yōu)化結(jié)果越好。

(a)ω=0 Hz

(b)ω=10 Hz

(c)ω=20 Hz

(d)ω=30 Hz

(e)ω=40 Hz圖7 不同頻率激勵(lì)載荷下懸臂板的拓?fù)鋬?yōu)化結(jié)果Fig.7 Topology optimization results of cantilever plate under different frequencies

表1 不同頻率下懸臂板目標(biāo)函數(shù)值的優(yōu)化結(jié)果Tab.1 Target function values of the cantilever plate at different frequencies

在結(jié)構(gòu)參數(shù)，約束條件，收斂條件，節(jié)點(diǎn)數(shù)設(shè)置相同的情況下，分別采用無(wú)單元Galerkin方法(EFG)和自然單元法(NEM)進(jìn)行激勵(lì)頻率ω=50 Hz時(shí)的結(jié)構(gòu)拓?fù)鋬?yōu)化。采用EFG進(jìn)行計(jì)算時(shí)，在結(jié)構(gòu)中分布40個(gè)背景積分網(wǎng)格，在每個(gè)積分網(wǎng)格中插入16個(gè)高斯積分點(diǎn)，影響域設(shè)置為1.5，以獲相似的初始目標(biāo)函數(shù)。所得的拓?fù)鋱D形如圖8、圖9所示。兩種計(jì)算方法的收斂曲線如圖10所示。NEM優(yōu)化單步耗時(shí)為24.26 s，結(jié)構(gòu)動(dòng)柔度由184 599減小到2 673。為EFG優(yōu)化單步耗時(shí)為30.13 s，結(jié)構(gòu)動(dòng)柔度由185 303減小到8 250。從計(jì)算結(jié)果看，EFG和NEM方法獲得的優(yōu)化結(jié)果均無(wú)棋盤(pán)格現(xiàn)象出現(xiàn)，就優(yōu)化效果來(lái)說(shuō)，NEM方法圖形較清晰，效果較好。從計(jì)算效率看，在初始動(dòng)柔度相似的情況下，NEM單步耗時(shí)和總耗時(shí)均較少，說(shuō)明自然單元法求解速度更快。

圖8 EFG激勵(lì)頻率ω=50 Hz懸臂板拓?fù)鋬?yōu)化結(jié)果Fig.8 Topology optimization result by using EFG(ω=50 Hz)

圖9 NEM激勵(lì)頻率ω=50 Hz懸臂板拓?fù)鋬?yōu)化結(jié)果Fig.9 Topology optimization result by using NEM(ω=50 Hz)

圖10 NEM與EFG拓?fù)鋬?yōu)化收斂圖Fig.10 Optimization convergence curves of EFG and NEM

造成以上現(xiàn)象的原因主要是由于EFG的形函數(shù)較為復(fù)雜，背景積分網(wǎng)格獨(dú)立于節(jié)點(diǎn)，積分網(wǎng)格的數(shù)量以及和排布方式對(duì)于計(jì)算結(jié)果影響較大，影響域的大小對(duì)計(jì)算結(jié)果以及計(jì)算速度也都會(huì)產(chǎn)生影響。而NEM的形函數(shù)計(jì)算較為簡(jiǎn)單，可以直接施加本征邊界條件，計(jì)算時(shí)將結(jié)構(gòu)自動(dòng)劃分為Delaunay三角形作為積分網(wǎng)格，三角形中分布高斯積分點(diǎn)后，高斯積分點(diǎn)依照Delaunay規(guī)則自動(dòng)尋找自然鄰點(diǎn)，排除了人為設(shè)置的影響域?qū)τ?jì)算結(jié)果的影響。

4.2 簡(jiǎn)支板

如圖11所示簡(jiǎn)支板，長(zhǎng)為10 m，寬為4 m，下邊界中點(diǎn)處加載頻率激勵(lì)載荷P=5eiwtkN，激勵(lì)幅值為5 kN。結(jié)構(gòu)的彈性模量為3×107Pa，泊松比為0.3，結(jié)構(gòu)的質(zhì)量密度為1 kg/m3，剛度矩陣以及質(zhì)量矩陣的懲罰因子同時(shí)設(shè)置為4。將結(jié)構(gòu)離散為29×11共319個(gè)節(jié)點(diǎn)，并剖分成560個(gè)Delaunay三角形，如圖12所示。在每個(gè)三角形內(nèi)插入3個(gè)高斯積分點(diǎn)。設(shè)置體積約束為50%，收斂條件設(shè)置為設(shè)計(jì)變量的最大變化小于0.01。

圖11 簡(jiǎn)支板結(jié)構(gòu)Fig.11 Structure of simply supported plate

圖12 簡(jiǎn)支板Delaunay三角形剖分Fig.12 Delaunay triangles of cantilever plate simply supported plate

如圖13所示為ω為50 Hz時(shí)，簡(jiǎn)支梁的優(yōu)化過(guò)程，在不斷迭代的過(guò)程中，結(jié)構(gòu)外形趨于光滑。在激勵(lì)頻率ω分別為50 Hz和100 Hz時(shí)，優(yōu)化結(jié)果如圖14、圖15所示。

(a)第5步

(b)第10步

(c)第30步

(d)第50步

(e)第55步

(e)第60步圖13 ω為50 Hz的拓?fù)鋬?yōu)化過(guò)程Fig.13 Topology optimization process of simply supported plate(ω=50 Hz)

圖14 激勵(lì)頻率ω為50 Hz拓?fù)鋬?yōu)化結(jié)果Fig.14 Topology optimization result under ω=50 Hz

圖15 激勵(lì)頻率ω為100 Hz拓?fù)鋬?yōu)化結(jié)果Fig.15 Topology optimization result under ω=100 Hz

5 結(jié) 論

本文進(jìn)行了頻率激勵(lì)作用下，基于動(dòng)力學(xué)自然單元法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化計(jì)算。應(yīng)用自然單元法，將結(jié)構(gòu)離散成節(jié)點(diǎn)，并將所離散出的節(jié)點(diǎn)的相對(duì)密度作為設(shè)計(jì)變量，利用SIMP模型，建立動(dòng)柔度最小化作為目標(biāo)函數(shù)并施加體積約束的拓?fù)鋬?yōu)化模型。采用伴隨分析法計(jì)算結(jié)構(gòu)的靈敏度，再利用優(yōu)化準(zhǔn)則法對(duì)模型進(jìn)行求解。

通過(guò)對(duì)懸臂板頻率激勵(lì)下結(jié)構(gòu)拓?fù)鋬?yōu)化的計(jì)算表明，當(dāng)激勵(lì)頻率不同時(shí)，得到了相似的拓?fù)鋬?yōu)化圖形。當(dāng)激勵(lì)頻率增大，初始動(dòng)柔度增大，優(yōu)化結(jié)果的柔度也隨之增大，優(yōu)化百分比隨之減小，說(shuō)明頻率激勵(lì)較低時(shí)，可獲得更好的優(yōu)化效果。在初始動(dòng)柔度相似的情況下，采用EFG和NEM方法分別進(jìn)行拓?fù)鋬?yōu)化，由于形函數(shù)構(gòu)造方式簡(jiǎn)單，且可直接施加本征邊界條件，相比于EFG方法，NEM顯示出具有更高的計(jì)算效率。同時(shí)，采用NEM方法可獲得更小的目標(biāo)函數(shù)，說(shuō)明該方法的優(yōu)化效果更好。通過(guò)對(duì)受頻率激勵(lì)的簡(jiǎn)支板進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化計(jì)算，表明隨著迭代步數(shù)的增加，結(jié)構(gòu)外形趨于光滑，最終獲得了清晰的結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果，進(jìn)一步說(shuō)明了該方法的可行性和有效性。

可以看出，動(dòng)力學(xué)自然單元法進(jìn)行頻率激勵(lì)下的結(jié)構(gòu)拓?fù)鋬?yōu)化計(jì)算，不僅可以克服有限元方法進(jìn)行拓?fù)鋬?yōu)化計(jì)算的網(wǎng)格依賴性以及拓?fù)鋬?yōu)化中經(jīng)常出現(xiàn)的棋盤(pán)格等數(shù)值不穩(wěn)定現(xiàn)象，而且相比其它無(wú)網(wǎng)格方法，具有更高的計(jì)算效率，說(shuō)明自然單元法應(yīng)用于頻率激勵(lì)拓?fù)鋬?yōu)化計(jì)算具有的可行性和優(yōu)越性。