周期性多孔結(jié)構(gòu)特征值拓?fù)鋬?yōu)化

付君健， 張 躍， 杜義賢， 高 亮

(1. 水電機(jī)械設(shè)備設(shè)計與維護(hù)湖北省重點(diǎn)實驗室， 湖北 宜昌 443002；2. 三峽大學(xué) 機(jī)械與動力學(xué)院， 湖北 宜昌 443002； 3. 華中科技大學(xué) 數(shù)字制造裝備與技術(shù)國家重點(diǎn)實驗室， 武漢 430074)

周期性多孔結(jié)構(gòu)是功構(gòu)一體化的優(yōu)良載體[1]，在航空航天、汽車、醫(yī)學(xué)植入等領(lǐng)域具有廣泛的應(yīng)用前景。周期性多孔結(jié)構(gòu)可設(shè)計性強(qiáng)，能根據(jù)功能特性需求，設(shè)計出不同形狀、尺寸和孔隙率的多孔構(gòu)型。對于含有多孔結(jié)構(gòu)的復(fù)雜機(jī)電系統(tǒng)，多孔結(jié)構(gòu)的固有頻率要求遠(yuǎn)離驅(qū)動系統(tǒng)的工作頻率。

在給定的材料屬性和邊界條件下，周期性多孔結(jié)構(gòu)的固有頻率與拓?fù)錁?gòu)型相關(guān)，通過改變多孔結(jié)構(gòu)的拓?fù)錁?gòu)型，可以被動地改變結(jié)構(gòu)固有頻率，使其脫離工作環(huán)境下的激勵頻率，從而避免劇烈共振導(dǎo)致的結(jié)構(gòu)失穩(wěn)和破壞。特征值拓?fù)鋬?yōu)化是改變多孔結(jié)構(gòu)固有頻率的有效設(shè)計方法，與宏觀結(jié)構(gòu)特征值拓?fù)鋬?yōu)化[2]的不同之處在于，周期性多孔結(jié)構(gòu)通過改變局部單胞的幾何特征，實現(xiàn)整體結(jié)構(gòu)特征值的優(yōu)化。

由于存在尺度跨越，特征值有限元求解方法對周期性多孔結(jié)構(gòu)拓?fù)鋬?yōu)化有重要影響。對于具有相同自由度數(shù)量的結(jié)構(gòu)系統(tǒng)，求解頻率和振型的特征值問題，其計算時間將比靜力分析高出一個數(shù)量級[3]。如果是局部幾何特征較多的多孔結(jié)構(gòu)，計算量將更大。由于尺度小、特性多、計算量大，多孔結(jié)構(gòu)宏觀等效屬性計算一般需采用模型縮減方法。均勻化法是常用的周期性多孔結(jié)構(gòu)宏觀等效屬性計算方法[4]，主要用于計算宏觀等效彈性矩陣，等效彈性矩陣則用于計算等效剛度矩陣和等效質(zhì)量矩陣[5]。均勻化法的假設(shè)之一為尺度分離，即多孔微結(jié)構(gòu)尺寸遠(yuǎn)遠(yuǎn)小于宏觀結(jié)構(gòu)尺寸。在實際應(yīng)用中，均勻化法更適合小尺度、大規(guī)模的周期性多孔結(jié)構(gòu)分析與優(yōu)化。對于有尺度要求的多孔結(jié)構(gòu)設(shè)計，需采用尺度關(guān)聯(lián)的模型縮減方法進(jìn)行宏觀等效屬性計算，如子結(jié)構(gòu)法[6]、尺度關(guān)聯(lián)均勻化法[7]、多尺度有限元法[8]，或者采用迭代法進(jìn)行求解[9]。

子結(jié)構(gòu)法作為一種較為成熟的模型縮減方法，近年來逐漸被用于細(xì)觀尺度周期性多孔結(jié)構(gòu)的拓?fù)鋬?yōu)化[10]。子結(jié)構(gòu)法具備尺度關(guān)聯(lián)的特性，適用于尺度關(guān)聯(lián)、中等規(guī)模多孔結(jié)構(gòu)分析與優(yōu)化。目前，基于子結(jié)構(gòu)法的多孔結(jié)構(gòu)拓?fù)鋬?yōu)化研究主要集中于靜力學(xué)問題[11-12]，實現(xiàn)周期性和多層級多孔結(jié)構(gòu)剛度拓?fù)鋬?yōu)化。雖然已有研究將子結(jié)構(gòu)法用于多層級結(jié)構(gòu)一階特征值拓?fù)鋬?yōu)化[13]，并取得了一定效果，但子結(jié)構(gòu)靜態(tài)凝聚求解高階特征值會產(chǎn)生較大誤差。因為子結(jié)構(gòu)靜態(tài)凝聚忽略了動力學(xué)方程中的慣性項，導(dǎo)致縮減的質(zhì)量矩陣為近似值，只有在頻率較低時才有較高的縮減精度，隨著頻率的增加，數(shù)值誤差將逐漸增大。O’Callahan[14]提出了針對動力學(xué)問題的改進(jìn)縮減系統(tǒng)(improved reduced system, IRS)，保證了子結(jié)構(gòu)法的動態(tài)縮減精度，但采用子結(jié)構(gòu)法進(jìn)行多孔結(jié)構(gòu)特征值分析與優(yōu)化仍存在諸多問題有待探索，如集中質(zhì)量施加、特征值尺度問題、特征值拓?fù)鋬?yōu)化敏度分析等。

本文基于子結(jié)構(gòu)動態(tài)凝聚方法，探討其動態(tài)縮減精度，建立了周期性多孔結(jié)構(gòu)特征值最大化拓?fù)鋬?yōu)化模型，實現(xiàn)了尺度關(guān)聯(lián)的二維、三維周期性多孔結(jié)構(gòu)的拓?fù)鋬?yōu)化，并揭示了細(xì)觀尺度多孔結(jié)構(gòu)特征值拓?fù)鋬?yōu)化的尺度效應(yīng)。

1 周期性多孔結(jié)構(gòu)幾何與物理描述

1.1 周期性多孔結(jié)構(gòu)幾何描述

針對周期性多孔結(jié)構(gòu)幾何描述，提出局部水平集函數(shù)(local level set function, LLSF)的概念，即在全局設(shè)計域中定義多個局部水平集函數(shù)描述細(xì)觀尺度下的周期性多孔結(jié)構(gòu)，將多孔結(jié)構(gòu)幾何邊界隱式地嵌入高一維的函數(shù)。局部水平集函數(shù)使得多孔結(jié)構(gòu)的幾何描述更加簡單，也使得多孔結(jié)構(gòu)拓?fù)鋬?yōu)化具有更多的設(shè)計自由[15]。如圖1所示，將二維周期性結(jié)構(gòu)設(shè)計域離散為M個細(xì)觀結(jié)構(gòu)，在每個細(xì)觀結(jié)構(gòu)內(nèi)部定義一個局部水平集函數(shù)，Φi表示其中一個細(xì)觀結(jié)構(gòu)的水平集函數(shù)

(a) 零水平面

(1)

式中：下標(biāo)i為設(shè)計域D內(nèi)的子結(jié)構(gòu)編號；x為空間坐標(biāo)向量；Ω為結(jié)構(gòu)域；?Ω為結(jié)構(gòu)邊界。

所有細(xì)觀結(jié)構(gòu)水平集函數(shù)Φi(xi)的集合組成了全局水平集函數(shù)Φ(x)

(2)

局部水平集函數(shù)Φi在子結(jié)構(gòu)內(nèi)部動態(tài)演化的Hamilton-Jacobi方程如下

(3)

式中：vn, i為局部水平集函數(shù)的法向速度場；t為時間變量。

為克服傳統(tǒng)水平集方法的數(shù)值問題，便于與梯度優(yōu)化算法結(jié)合，本文采用緊支徑向基函數(shù)插值替代離散的水平集函數(shù)，形成參數(shù)化水平集法[16]。在參數(shù)化水平集法中引入局部水平集函數(shù)，可有效避免初始擴(kuò)展系數(shù)求解困難的問題[17-18]。

(4)

式中：φ為緊支徑向基函數(shù)插值矩陣；α為擴(kuò)展系數(shù)向量，其中，k=1,…,N，N為細(xì)觀結(jié)構(gòu)設(shè)計域緊支徑向基插值控制點(diǎn)的數(shù)量。

根據(jù)式(4)，參數(shù)化之后的局部水平集函數(shù)法向速度場為

(5)

1.2 周期性多孔結(jié)構(gòu)物理描述

本文采用子結(jié)構(gòu)法[19]描述多孔結(jié)構(gòu)物理模型，如圖2所示。將多孔結(jié)構(gòu)單胞凝聚為具有多個節(jié)點(diǎn)的超單元。與周期性多孔結(jié)構(gòu)的幾何描述定義類似，子結(jié)構(gòu)法的基本思想是將計算域Ω分解為若干子區(qū)域Ωi，先求解邊界上的數(shù)值信息，然后對子區(qū)域內(nèi)部“各個擊破”。在有限元分析中采用子結(jié)構(gòu)法可以實現(xiàn)多孔結(jié)構(gòu)物理模型和幾何模型的匹配，縮小有限元分析的計算規(guī)模，降低計算內(nèi)存的消耗，提高周期性多孔結(jié)構(gòu)特征值分析與拓?fù)鋬?yōu)化的效率。

圖2 子結(jié)構(gòu)凝聚

特征值問題的有限元平衡方程為一個無阻尼自由振動系統(tǒng)，為了展示子結(jié)構(gòu)法動態(tài)凝聚也適用于考慮載荷的動力學(xué)問題，在特征值問題中引入載荷向量F=0，則特征值問題平衡方程的矩陣形式為

(6)

(7)

由式(7)第二行可得主節(jié)點(diǎn)位移和從節(jié)點(diǎn)位移之間的關(guān)系

(8)

對于靜態(tài)凝聚，忽略式(8)中所有的慣性項

KsmUm+KssUs=Fs

假設(shè)Kss為正定矩陣，且結(jié)構(gòu)內(nèi)部的從節(jié)點(diǎn)不受外載荷作用，即Fs=0

(10)

由式(10)得到變化關(guān)系

(11)

式中，Tc為靜態(tài)凝聚的變換矩陣

(12)

根據(jù)式(11)和式(12)，縮減后的平衡方程為

(13)

其中

(14)

(15)

由于式(9)忽略了全部慣性項，對于特征值問題，式(13)的縮減精度存在較大誤差，需要對此進(jìn)行改進(jìn)。由式(13)可得

(16)

對式(10)求二階導(dǎo)數(shù)

(17)

將式(16)代入式(17)

(18)

將式(16)和式(18)代入式(8)

(19)

對于無阻尼自由振動問題，主、從節(jié)點(diǎn)上的載荷Fs和Fm均為零，式(19)可進(jìn)一步簡化為

(20)

由此得到子結(jié)構(gòu)動態(tài)凝聚的變換矩陣

(21)

其中

經(jīng)過子結(jié)構(gòu)動態(tài)凝聚后的無阻尼自由振動平衡方程為

(23)

其中，縮減后的質(zhì)量矩陣和剛度矩陣分別為

(24)

(25)

單一結(jié)構(gòu)的凝聚可通過式(23)實現(xiàn)，周期性多孔結(jié)構(gòu)則涉及多個子結(jié)構(gòu)的凝聚。由于周期性排布的假設(shè)條件，含有重復(fù)結(jié)構(gòu)特征的子結(jié)構(gòu)只需進(jìn)行一次凝聚，多個重復(fù)子結(jié)構(gòu)的凝聚過程可通過對單個凝聚子結(jié)構(gòu)的組裝完成。由于載荷向量一般取零，多個重復(fù)子結(jié)構(gòu)組裝而成的縮減平衡方程為

(26)

(27)

(28)

為驗證子結(jié)構(gòu)動態(tài)凝聚求解特征值問題的精度，與子結(jié)構(gòu)靜態(tài)凝聚、全自由度有限元模型進(jìn)行對比分析。如圖3所示為對比分析的模型，尺寸為1 m×1 m二維正方形板，厚度為0.1 m，左下角和右下角固定約束，離散為40×40個雙線性四邊形單元，單元尺寸為0.025 m×0.025 m。材料彈性模量為200 GPa，泊松比為0.3，質(zhì)量密度為7 800 kg/m3。

(a) 全自由度模型

表1為有限元分析的特征頻率值，與全自由度有限元分析對比可知，子結(jié)構(gòu)靜態(tài)凝聚求解的特征值誤差較大，隨著階數(shù)的增大，特征頻率誤差越來越大。子結(jié)構(gòu)動態(tài)凝聚求解的特征值誤差較小，且誤差不會隨著階數(shù)的增加而放大。因此，子結(jié)構(gòu)動態(tài)凝聚能保證特征值求解的精度，同時能縮減結(jié)構(gòu)自由度數(shù)量，提高特征值問題的求解效率。

表1 特征頻率對比

2 周期性多孔結(jié)構(gòu)特征值拓?fù)鋬?yōu)化模型

多孔結(jié)構(gòu)的共振很大程度取決于前幾階頻率，劇烈的共振一般發(fā)生在外界激勵頻率和某一階固有頻率接近的時候，所以在多孔結(jié)構(gòu)設(shè)計時，必須對結(jié)構(gòu)的前幾階固有頻率進(jìn)行限制或改變。本文以前幾階加權(quán)特征值最大化為目標(biāo)函數(shù)，構(gòu)造如式所示的p-norm形式的函數(shù)，并以此緩解模態(tài)轉(zhuǎn)換引起的目標(biāo)函數(shù)振蕩問題[20]。該目標(biāo)函數(shù)可通過改變特征值權(quán)重ωj和p-norm函數(shù)的p值，靈活調(diào)整對全部或特定階次特征值的影響。

(29)

(30)

(31)

(32)

式中：ε為應(yīng)變場；D為實體材料的彈性矩陣；ρ為材料質(zhì)量密度。

由于特征值λ是針對多孔結(jié)構(gòu)整體的特征值，優(yōu)化模型式(30)中的目標(biāo)函數(shù)并未寫成對每個多孔結(jié)構(gòu)積分或累加的形式，但靈敏度表達(dá)式可寫成對子結(jié)構(gòu)積分或累加的形式，其具體形式見靈敏度分析。

3 靈敏度分析與數(shù)值實施

3.1 靈敏度分析

采用形狀導(dǎo)數(shù)來推導(dǎo)目標(biāo)函數(shù)和約束條件關(guān)于設(shè)計變量(擴(kuò)展系數(shù))的靈敏度。根據(jù)形狀導(dǎo)數(shù)的定義及其引理[21]，目標(biāo)函數(shù)、能量雙線性形式和載荷雙線性形式關(guān)于時間變量t的導(dǎo)數(shù)分別為：

(33)

(34)

(35)

將式(30)中的平衡方程兩邊對時間t求偏導(dǎo)得到

(36)

(37)

將式(34)、式(35)、式(37)代入式(36)可得

(38)

通過構(gòu)造伴隨方程可得，該優(yōu)化問題為自伴隨問題[22]，即v=u。

(39)

根據(jù)質(zhì)量矩陣的正交歸一化條件[23]

b(u,u)=1

(40)

式(38)可進(jìn)一步簡化為

(41)

式(41)仍存在邊界積分形式，引入式(42)

dΓ=δ(Φ)|?Φ|dΩ

(42)

式(41)的邊界積分轉(zhuǎn)換為體積積分

|?Φ|vndΩ

(43)

由于本文采用子結(jié)構(gòu)法，式(43)可寫成對每個子結(jié)構(gòu)積分的形式

(44)

γi,j(ui,λj)=εT(ui,j)Dε(ui,j)-λjρiui,jui,j

(45)

式中：下標(biāo)i為子結(jié)構(gòu)的編號；下標(biāo)j為特征值的階次。

將式(5)中的局部水平集函數(shù)的速度場vn,i代入式(44)可得

(46)

將特征值關(guān)于時間t的導(dǎo)數(shù)代入式(33)可得

(47)

通過鏈?zhǔn)椒▌t對目標(biāo)函數(shù)J(Φ)直接求其關(guān)于時間變量t的偏導(dǎo)數(shù)可以得到

(48)

對比式(47)和式(48)可得目標(biāo)函數(shù)關(guān)于擴(kuò)展系數(shù)α的導(dǎo)數(shù)

φi,k(xi)dΩidΩ

(49)

同理，體積約束對擴(kuò)展系數(shù)α的導(dǎo)數(shù)為

(50)

3.2 數(shù)值實施

結(jié)構(gòu)特征值的求解是在縮減的全局有限元模型上實現(xiàn)，該過程暫不涉及子結(jié)構(gòu)內(nèi)部的擴(kuò)展。在求解敏度信息時，需要求解子結(jié)構(gòu)內(nèi)部的特征向量。子結(jié)構(gòu)內(nèi)部特征向量求解可順序求解或者CPU并行求解，CPU并行求解在子結(jié)構(gòu)劃分?jǐn)?shù)量較多時能節(jié)約大量時間。在求解時為避免重復(fù)計算，在子結(jié)構(gòu)動態(tài)縮減時，可將式(20)中相關(guān)的數(shù)據(jù)重復(fù)使用。此外，在敏度計算時，子結(jié)構(gòu)內(nèi)部特征向量的求解可重復(fù)使用子結(jié)構(gòu)動態(tài)縮減中式(21)的值，這樣TIRS只需要計算一次。

為防止優(yōu)化過程中隨著結(jié)構(gòu)材料的不斷刪除，出現(xiàn)無限多個特征值，需要在結(jié)構(gòu)的特定位置加入所謂的非結(jié)構(gòu)集中質(zhì)量。如果不加入非結(jié)構(gòu)集中質(zhì)量，在結(jié)構(gòu)的特定位置就不可能獲得期望的最優(yōu)拓?fù)浣Y(jié)構(gòu)，從而影響到整個結(jié)構(gòu)頻率的最大化[24]。在周期性多孔結(jié)構(gòu)設(shè)計中，單方向的周期數(shù)量可能為奇數(shù)或者偶數(shù)，非結(jié)構(gòu)集中質(zhì)量的施加位置視具體情況而定，非結(jié)構(gòu)集中質(zhì)量可施加在如圖4(a)所示的子結(jié)構(gòu)內(nèi)部自由度上，也可施加在如圖4(b)所示的縮減結(jié)構(gòu)的自由度上。加在子結(jié)構(gòu)內(nèi)部時，需在子結(jié)構(gòu)內(nèi)部質(zhì)量矩陣自由度所在位置的對角線上添加一個很大的數(shù)值，大小約為允許設(shè)計質(zhì)量的10%左右，非結(jié)構(gòu)集中質(zhì)量會參與子結(jié)構(gòu)的動態(tài)凝聚。加在縮減的結(jié)構(gòu)自由度上時，需在子結(jié)構(gòu)凝聚和組裝完成之后，在縮減結(jié)構(gòu)的質(zhì)量矩陣自由度所在位置的對角線上添加設(shè)計質(zhì)量10%左右的數(shù)值。

(a) 子結(jié)構(gòu)內(nèi)部

4 數(shù)值案例與結(jié)果分析

4.1 懸臂梁結(jié)構(gòu)

如圖5所示的二維懸臂梁結(jié)構(gòu)，梁結(jié)構(gòu)設(shè)計域長寬比L∶H= 2 m∶1 m，厚度為0.1 m，將宏觀結(jié)構(gòu)設(shè)計域劃分為20×10個子結(jié)構(gòu)，子結(jié)構(gòu)內(nèi)部劃分為40×40個雙線性四邊形單元，單元的尺寸為0.025 m×0.025 m，材料彈性模量E=200 GPa，泊松比v= 0.3，質(zhì)量密度7 800 kg/m3，在宏觀結(jié)構(gòu)右邊界中點(diǎn)處施加大小為設(shè)計質(zhì)量10%的集中質(zhì)量，以多孔結(jié)構(gòu)體積分?jǐn)?shù)f= 0.5為約束，進(jìn)行二維周期性多孔結(jié)構(gòu)設(shè)計。

圖5 懸臂梁設(shè)計域

圖6所示為經(jīng)過優(yōu)化的子結(jié)構(gòu)單胞、水平集函數(shù)及其組裝而成的宏觀結(jié)構(gòu)，經(jīng)過162次迭代后，其最優(yōu)的目標(biāo)函數(shù)值收斂到16.103，體積分?jǐn)?shù)收斂到0.5，前6階頻率分別為0.398 Hz、31.621 Hz、55.398 Hz、60.04 Hz、94.633 Hz、110.121 Hz。圖7所示為前6階特征頻率對應(yīng)的特征向量圖，即結(jié)構(gòu)的模態(tài)振型。值得注意的是，數(shù)值算例中周期性多孔結(jié)構(gòu)一階固有頻率較低，通過數(shù)值試驗分析可知，引入非結(jié)構(gòu)集中質(zhì)量是導(dǎo)致該現(xiàn)象的主要原因。因為在整個宏觀設(shè)計域中僅定義了子結(jié)構(gòu)的局部水平集函數(shù)，因此只需要對一個子結(jié)構(gòu)單胞進(jìn)行求解，即可得到如圖6(c)所示的20×10個周期性分布的子結(jié)構(gòu)。從圖6(c)可以看出，最優(yōu)的周期性結(jié)構(gòu)在集中質(zhì)量施加點(diǎn)處存在材料，這是由于子結(jié)構(gòu)動態(tài)縮減的尺度關(guān)聯(lián)性所帶來的益處。

(a) 最優(yōu)子結(jié)構(gòu)單胞

(a)

在迭代初期前6階頻率一直處于下降趨勢，這是因為初始化的結(jié)構(gòu)體積分?jǐn)?shù)大于指定的體積分?jǐn)?shù)，迭代初期的約束條件不滿足，導(dǎo)致前6階頻率較大，如圖8所示。當(dāng)體積分?jǐn)?shù)約束得到滿足后，隨著優(yōu)化過程的進(jìn)行，前6階頻率逐漸收斂。

圖8 前6階頻率迭代圖

4.2 兩端固定梁結(jié)構(gòu)

為驗證本方法的高效性和特征值問題的尺度效應(yīng)，考慮如圖9所示的兩端固定的二維梁結(jié)構(gòu)，將本方法與沒有模型縮減的全自由度拓?fù)鋬?yōu)化[25]的計算時間進(jìn)行了對比。梁結(jié)構(gòu)設(shè)計域長寬比L∶H= 7 m∶1 m，厚度為0.1 m，將宏觀結(jié)構(gòu)設(shè)計域分別劃分為14×2、28×4、56×8和84×12個子結(jié)構(gòu)，子結(jié)構(gòu)內(nèi)部劃分為40×40個雙線性四邊形單元，單元的尺寸為0.025 m×0.025 m，材料彈性模量E=200 GPa，泊松比v= 0.3，質(zhì)量密度7 800 kg/m3。在結(jié)構(gòu)中心施加大小為設(shè)計質(zhì)量10%的集中質(zhì)量，以子結(jié)構(gòu)體積分?jǐn)?shù)f= 0.5為約束，進(jìn)行二維周期性多孔結(jié)構(gòu)設(shè)計。

圖9 兩端固定二維梁設(shè)計域

當(dāng)宏觀結(jié)構(gòu)離散為14×2、28×4、56×8和84×12個子結(jié)構(gòu)時，其結(jié)構(gòu)自由度數(shù)量分別約為9萬、36萬、144萬和323萬，采用子結(jié)構(gòu)法分別將其縮減至約為0.57萬、2萬、7.6萬、17萬。如表2所示為對應(yīng)4種子結(jié)構(gòu)數(shù)量下的拓?fù)鋬?yōu)化結(jié)果。其中最優(yōu)子結(jié)構(gòu)拓?fù)錁?gòu)型的特征在于上下兩端有橫向分布的主承載結(jié)構(gòu)，主承載結(jié)構(gòu)之間有交叉的結(jié)構(gòu)。隨著子結(jié)構(gòu)數(shù)量的增多，最優(yōu)子結(jié)構(gòu)的拓?fù)錁?gòu)型在局部展現(xiàn)出一定的變化，其目標(biāo)函數(shù)也隨之下降。圖10所示對應(yīng)4種子結(jié)構(gòu)數(shù)量下周期多孔結(jié)構(gòu)特征頻率的變化趨勢，在固定大小的設(shè)計域內(nèi)，當(dāng)離散的子結(jié)構(gòu)數(shù)量越多時，優(yōu)化結(jié)果表現(xiàn)出一定的尺度效應(yīng)，其最優(yōu)多孔結(jié)構(gòu)的前6階特征頻率整體呈下降趨勢。因此，在實際工程應(yīng)用中，可

圖10 多孔結(jié)構(gòu)特征頻率變化

表2 特征值拓?fù)鋬?yōu)化結(jié)果

針對結(jié)構(gòu)具體的力學(xué)性能需求來設(shè)置多孔結(jié)構(gòu)周期性排布的數(shù)量，從而避開結(jié)構(gòu)的激振頻率。

在拓?fù)鋬?yōu)化的計算效率方面，如圖10所示，對比本文子結(jié)構(gòu)法與全自由度模型，當(dāng)子結(jié)構(gòu)數(shù)量為14×2時，在每一步迭代中，兩種方法在有限元分析和拓?fù)鋬?yōu)化的平均計算時間相當(dāng)。當(dāng)子結(jié)構(gòu)數(shù)量大于14×2后，在每一步迭代中，本文方法在有限元分析和拓?fù)鋬?yōu)化方面的平均計算時間大幅下降。當(dāng)子結(jié)構(gòu)數(shù)量為84×12時，本文方法優(yōu)化更新的計算時間僅為全自由度模型的19.4%，有限元分析的計算時間僅為全自由模型的31.9%，平均每步的迭代時間約為全自由度模型的30.5%。其原因在于，本文采用了子結(jié)構(gòu)動態(tài)凝聚縮減了自由度數(shù)量，采用了局部水平集函數(shù)描述多孔結(jié)構(gòu)降低了幾何描述上的復(fù)雜度，從幾何模型和物理模型角度共同提升了周期性多孔結(jié)構(gòu)特征值拓?fù)鋬?yōu)化的效率。此外，從圖11可知，有限元分析占據(jù)了拓?fù)鋬?yōu)化大部分的計算時間，研究高效的有限元分析方法將有助于提升拓?fù)鋬?yōu)化的計算效率。

圖11 計算效率對比

4.3 三維梁結(jié)構(gòu)

如圖12所示的兩端固定的三維梁結(jié)構(gòu)，結(jié)構(gòu)設(shè)計域長寬高比例為L∶W∶H= 7 m∶1 m∶1 m，將結(jié)構(gòu)設(shè)計域劃分為14×2×2個三維子結(jié)構(gòu)，每個子結(jié)構(gòu)離散為10×10×10個八節(jié)點(diǎn)六面體單元，單元的尺寸為0.1m× 0.1 m × 0.1 m，材料彈性模量E= 200 GPa，泊松比v= 0.3，質(zhì)量密度7 800 kg/m3。在結(jié)構(gòu)中心施加大小為設(shè)計質(zhì)量10%的集中質(zhì)量，以多孔結(jié)構(gòu)體積分?jǐn)?shù)f= 0.4為約束，進(jìn)行三維周期性多孔結(jié)構(gòu)設(shè)計。

三維結(jié)構(gòu)的有限元分析和子結(jié)構(gòu)的凝聚相對于二維問題具有更大的計算量，本算例采用三維子結(jié)構(gòu)法將結(jié)構(gòu)自由度從18.6萬左右縮減至約4.2萬。如圖13所示為經(jīng)過優(yōu)化的三維周期性結(jié)構(gòu)及其子結(jié)構(gòu)單胞，從子結(jié)構(gòu)單胞拓?fù)錁?gòu)型可知，組裝后的周期性多孔結(jié)構(gòu)在集中質(zhì)量點(diǎn)處必然存在材料，這是子結(jié)構(gòu)法尺度關(guān)聯(lián)所帶來的益處。從圖14可知，拓?fù)鋬?yōu)化模型經(jīng)過93次迭代后，其最優(yōu)的目標(biāo)函數(shù)值收斂至36.564 7，體積分?jǐn)?shù)收斂至0.4。圖15所示為前6階頻率迭代圖，分別收斂至1.558 Hz、36.049 Hz、68.633 Hz、82.604 Hz、82.814 Hz、106.298 Hz。

圖13 三維周期性多孔結(jié)構(gòu)優(yōu)化結(jié)果

(a)

圖15 前6階特征頻率迭代圖

5 結(jié) 論

子結(jié)構(gòu)法是高效、高精度的模型縮減方法，本文將參數(shù)化水平集方法和子結(jié)構(gòu)動態(tài)凝聚相結(jié)合，研究了周期性多孔結(jié)構(gòu)的特征值拓?fù)鋬?yōu)化方法。研究表明，本方法能有效實現(xiàn)尺度關(guān)聯(lián)的二維和三維周期性多孔結(jié)構(gòu)的特征值拓?fù)鋬?yōu)化，且所設(shè)計的周期性多孔結(jié)構(gòu)具有一定的尺度效應(yīng)。