周期性多孔結構特征值拓撲優化

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

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

周期性多孔結構是功構一體化的優良載體[1]，在航空航天、汽車、醫學植入等領域具有廣泛的應用前景。周期性多孔結構可設計性強，能根據功能特性需求，設計出不同形狀、尺寸和孔隙率的多孔構型。對于含有多孔結構的復雜機電系統，多孔結構的固有頻率要求遠離驅動系統的工作頻率。

在給定的材料屬性和邊界條件下，周期性多孔結構的固有頻率與拓撲構型相關，通過改變多孔結構的拓撲構型，可以被動地改變結構固有頻率，使其脫離工作環境下的激勵頻率，從而避免劇烈共振導致的結構失穩和破壞。特征值拓撲優化是改變多孔結構固有頻率的有效設計方法，與宏觀結構特征值拓撲優化[2]的不同之處在于，周期性多孔結構通過改變局部單胞的幾何特征，實現整體結構特征值的優化。

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

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

本文基于子結構動態凝聚方法，探討其動態縮減精度，建立了周期性多孔結構特征值最大化拓撲優化模型，實現了尺度關聯的二維、三維周期性多孔結構的拓撲優化，并揭示了細觀尺度多孔結構特征值拓撲優化的尺度效應。

1 周期性多孔結構幾何與物理描述

1.1 周期性多孔結構幾何描述

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

(a) 零水平面

(1)

式中：下標i為設計域D內的子結構編號；x為空間坐標向量；Ω為結構域；?Ω為結構邊界。

所有細觀結構水平集函數Φi(xi)的集合組成了全局水平集函數Φ(x)

(2)

局部水平集函數Φi在子結構內部動態演化的Hamilton-Jacobi方程如下

(3)

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

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

(4)

式中：φ為緊支徑向基函數插值矩陣；α為擴展系數向量，其中，k=1,…,N，N為細觀結構設計域緊支徑向基插值控制點的數量。

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

(5)

1.2 周期性多孔結構物理描述

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

圖2 子結構凝聚

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

(6)

(7)

由式(7)第二行可得主節點位移和從節點位移之間的關系

(8)

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

KsmUm+KssUs=Fs

假設Kss為正定矩陣，且結構內部的從節點不受外載荷作用，即Fs=0

(10)

由式(10)得到變化關系

(11)

式中，Tc為靜態凝聚的變換矩陣

(12)

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

(13)

其中

(14)

(15)

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

(16)

對式(10)求二階導數

(17)

將式(16)代入式(17)

(18)

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

(19)

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

(20)

由此得到子結構動態凝聚的變換矩陣

(21)

其中

經過子結構動態凝聚后的無阻尼自由振動平衡方程為

(23)

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

(24)

(25)

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

(26)

(27)

(28)

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

(a) 全自由度模型

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

表1 特征頻率對比

2 周期性多孔結構特征值拓撲優化模型

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

(29)

(30)

(31)

(32)

式中：ε為應變場；D為實體材料的彈性矩陣；ρ為材料質量密度。

由于特征值λ是針對多孔結構整體的特征值，優化模型式(30)中的目標函數并未寫成對每個多孔結構積分或累加的形式，但靈敏度表達式可寫成對子結構積分或累加的形式，其具體形式見靈敏度分析。

3 靈敏度分析與數值實施

3.1 靈敏度分析

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

(33)

(34)

(35)

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

(36)

(37)

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

(38)

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

(39)

根據質量矩陣的正交歸一化條件[23]

b(u,u)=1

(40)

式(38)可進一步簡化為

(41)

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

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

(42)

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

|?Φ|vndΩ

(43)

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

(44)

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

(45)

式中：下標i為子結構的編號；下標j為特征值的階次。

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

(46)

將特征值關于時間t的導數代入式(33)可得

(47)

通過鏈式法則對目標函數J(Φ)直接求其關于時間變量t的偏導數可以得到

(48)

對比式(47)和式(48)可得目標函數關于擴展系數α的導數

φi,k(xi)dΩidΩ

(49)

同理，體積約束對擴展系數α的導數為

(50)

3.2 數值實施

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

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

(a) 子結構內部

4 數值案例與結果分析

4.1 懸臂梁結構

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

圖5 懸臂梁設計域

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

(a) 最優子結構單胞

(a)

在迭代初期前6階頻率一直處于下降趨勢，這是因為初始化的結構體積分數大于指定的體積分數，迭代初期的約束條件不滿足，導致前6階頻率較大，如圖8所示。當體積分數約束得到滿足后，隨著優化過程的進行，前6階頻率逐漸收斂。

圖8 前6階頻率迭代圖

4.2 兩端固定梁結構

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

圖9 兩端固定二維梁設計域

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

圖10 多孔結構特征頻率變化

表2 特征值拓撲優化結果

針對結構具體的力學性能需求來設置多孔結構周期性排布的數量，從而避開結構的激振頻率。

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

圖11 計算效率對比

4.3 三維梁結構

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

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

圖13 三維周期性多孔結構優化結果

(a)

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

5 結 論

子結構法是高效、高精度的模型縮減方法，本文將參數化水平集方法和子結構動態凝聚相結合，研究了周期性多孔結構的特征值拓撲優化方法。研究表明，本方法能有效實現尺度關聯的二維和三維周期性多孔結構的特征值拓撲優化，且所設計的周期性多孔結構具有一定的尺度效應。