吳紫俊 肖人彬

摘要 針對傳統(tǒng)基于均勻化方法的周期性結(jié)構(gòu)設(shè)計宏微尺度分離問題，提出了基于有限元子結(jié)構(gòu)的周期性結(jié)構(gòu)設(shè)計方法。在該方法中，假設(shè)宏觀結(jié)構(gòu)由構(gòu)型相同的多個子結(jié)構(gòu)構(gòu)成。在宏觀層面，利用子結(jié)構(gòu)凝聚構(gòu)建了超單元計算模型;在微觀層面，采用SIMP優(yōu)化方法建立了子結(jié)構(gòu)內(nèi)部單元的拓?fù)錁?gòu)型設(shè)計方法。通過反求超單元模型計算結(jié)果，實現(xiàn)了宏微觀結(jié)構(gòu)的耦合計算。在耦合計算過程中，由于超單元之間具有公共的部分邊界，使得子結(jié)構(gòu)構(gòu)型間具有良好的材料連通性。同時，子結(jié)構(gòu)劃分源于宏觀設(shè)計域的有限元網(wǎng)格，確保了宏微結(jié)構(gòu)間具有明確的尺寸比例關(guān)系。最后，通過兩類梁的計算實例，驗證了所提方法的有效性。

關(guān) 鍵 詞 子結(jié)構(gòu);周期性;微結(jié)構(gòu);拓?fù)鋬?yōu)化

中圖分類號 TB383.3 ? ? 文獻(xiàn)標(biāo)志碼 A

0 引言

周期性結(jié)構(gòu)設(shè)計是現(xiàn)代結(jié)構(gòu)設(shè)計的重要分支[1]，已被廣泛應(yīng)用于多個領(lǐng)域，如航空[2]、汽車[3]等的結(jié)構(gòu)設(shè)計。周期性結(jié)構(gòu)不僅承載能力強(qiáng)，重量輕[4-5]，還具有隔熱[6]、隔音[7]等多功能物理特性，因其獨(dú)特的結(jié)構(gòu)組成形式和良好的設(shè)計加工性能，展現(xiàn)出了廣闊的研究和應(yīng)用前景[8-9]。周期性結(jié)構(gòu)設(shè)計是一種多尺度的拓?fù)鋬?yōu)化設(shè)計，方法上可分為2類：一類是基于均勻化的設(shè)計方法[10];另一類是基于宏觀周期性約束的設(shè)計方法[11]。

均勻化方法設(shè)定材料的宏觀結(jié)構(gòu)由微觀尺度的單胞周期性拓展而形成。周期性結(jié)構(gòu)的單胞具有周期性應(yīng)力和應(yīng)變邊界條件，從而獲取單胞內(nèi)材料的最優(yōu)分布，實現(xiàn)不同結(jié)構(gòu)性能的周期性單胞設(shè)計[12]。Vogiatzis等[13]基于水平集法提出了單相和多相負(fù)泊松比材料的設(shè)計方法，通過施加周期性的Dirichlet位移邊界條件得到周期性微結(jié)構(gòu)，并采用最速下降法構(gòu)建微結(jié)構(gòu)內(nèi)部構(gòu)型。Fantoni等[14]提出了一種多場漸近均勻化方法，用以分析具有周期性微結(jié)構(gòu)的壓電材料，利用周期攝動函數(shù)考慮具有代表性的微結(jié)構(gòu)非均勻性的影響。Hunag等[15]將雙向進(jìn)化結(jié)構(gòu)優(yōu)化方法推廣到具有最佳粘彈性特性的復(fù)合材料微觀結(jié)構(gòu)設(shè)計中，基于均勻化方法計算復(fù)合材料的性能參數(shù)，通過BESO方法實現(xiàn)單胞內(nèi)材料的重新分配。賈嬌等[16]研究了宏觀傳導(dǎo)條件對材料微結(jié)構(gòu)的影響，并提出了基于宏觀傳導(dǎo)條件的周期性結(jié)構(gòu)傳熱材料的研究，構(gòu)建了基于變密度法周期性傳熱結(jié)構(gòu)模型。杜義賢等[17]基于能量均勻化方法研究了抗剪性能更優(yōu)的周期性點(diǎn)陣結(jié)構(gòu)，并提出了基于經(jīng)典胞元理論的周期性點(diǎn)陣結(jié)構(gòu)的面內(nèi)等效切變模量的解析表達(dá)式。

在均勻化方法中，微觀單胞與宏觀結(jié)構(gòu)設(shè)計是相互耦合的設(shè)計過程，對微觀結(jié)構(gòu)的宏觀有效性材料進(jìn)行評價的同時，宏觀結(jié)構(gòu)的計算結(jié)果也影響微觀結(jié)構(gòu)的構(gòu)型優(yōu)化[18]。采用均勻化方法進(jìn)行周期性結(jié)構(gòu)設(shè)計，其尺度分離假設(shè)會帶來2個問題：一個是微結(jié)構(gòu)單胞間材料的連通性[19]，其微結(jié)構(gòu)單胞優(yōu)化是一個獨(dú)立的優(yōu)化過程，宏觀位移或應(yīng)力的變化往往引起微結(jié)構(gòu)間材料不連通，使得所優(yōu)化的結(jié)果只有理論構(gòu)型，不具備制造性[20];另一個是均勻化方法假設(shè)微結(jié)構(gòu)與宏觀結(jié)構(gòu)的尺寸比例約為10-9，當(dāng)該比例接近于1時往往得不到合適的優(yōu)化構(gòu)型[21-22]。

宏觀周期性約束的設(shè)計方法中，宏觀結(jié)構(gòu)被均勻地劃分成若干個具有特定長度比例大小的單元胞體。每個單元胞體同一位置上的單元具有相同的材料密度值和靈敏度值[23]。優(yōu)化變量為單元胞體中每個單元的材料密度。與均勻化方法相比，隨著單元胞體的細(xì)化，其優(yōu)化結(jié)果逐漸近似于均勻化方法設(shè)計結(jié)果。在基于BESO方法的周期性結(jié)構(gòu)設(shè)計中[24]，需要在整個設(shè)計域中求解具有周期性邊界約束的方程，存在計算量大的問題[25]。

為解決均勻化方法的尺度分離問題和宏觀周期性約束設(shè)計方法中計算量過大的問題[26]，本文在有限元子結(jié)構(gòu)方法的基礎(chǔ)上，提出了利用子結(jié)構(gòu)凝聚構(gòu)建宏觀結(jié)構(gòu)的超單元計算模型，采用SIMP優(yōu)化方法[27-28]建立子結(jié)構(gòu)內(nèi)部單元拓?fù)錁?gòu)型設(shè)計的周期性結(jié)構(gòu)設(shè)計方法。

1 子結(jié)構(gòu)凝聚及周期性結(jié)構(gòu)優(yōu)化問題定義

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

根據(jù)有限元子結(jié)構(gòu)的定義，有限元網(wǎng)格模型可以重新劃分為多個子結(jié)構(gòu)[29]，如圖1所示。

優(yōu)化目標(biāo)是宏觀周期性結(jié)構(gòu)的最小化柔度;設(shè)計變量為周期性子結(jié)構(gòu)中包含的n個單元的相對密度[ρj];[K]、F分別為子結(jié)構(gòu)對應(yīng)超單元組裝的整體剛度矩陣與外載荷向量，[Ubb，i]是第i個超單元的位移向量，[V]是當(dāng)前周期性子結(jié)構(gòu)所含的材料體積，N為宏觀結(jié)構(gòu)的子結(jié)構(gòu)總數(shù)。

2 基于子結(jié)構(gòu)的周期性結(jié)構(gòu)優(yōu)化模型

2.1 宏觀尺度模型

在周期性結(jié)構(gòu)的宏微結(jié)構(gòu)優(yōu)化中，宏觀結(jié)構(gòu)的剛度矩陣由式（5）建立。由于省略了內(nèi)部節(jié)點(diǎn)，可提高其計算效率。而子結(jié)構(gòu)的內(nèi)部節(jié)點(diǎn)的位移可根據(jù)式（2）獲得。此時，圖1中的子結(jié)構(gòu)模型的剛度矩陣可表示為

此時，所優(yōu)化的結(jié)構(gòu)將分布在整個宏觀設(shè)計域中。在宏觀結(jié)構(gòu)計算過程中，其構(gòu)建剛度矩陣的周期性子結(jié)構(gòu)均有相同構(gòu)型，對應(yīng)的每個超單元凝聚矩陣均相同，即[K*bb，t=K*bb，s，t≠s] 。

2.2 子結(jié)構(gòu)優(yōu)化模型

在子結(jié)構(gòu)構(gòu)型的優(yōu)化中，需要增加懲罰因子p對拓?fù)鋬?yōu)化設(shè)計變量的密度進(jìn)行懲罰，使其收斂于設(shè)定的密度上下界，從而獲得子結(jié)構(gòu)的“黑白”設(shè)計[30]：

式中，懲罰因子[p≥1]表明具有中間密度的子結(jié)構(gòu)[0<ρj<1]的剛度矩陣會被懲罰。該懲罰因子值越大，其優(yōu)化結(jié)構(gòu)桿系特征越明顯。式中，n表示該子結(jié)構(gòu)內(nèi)部的有限元單元總數(shù)，[Κ0]為實體材料下的有限元單元剛度矩陣，如圖2所示。

與傳統(tǒng)的SIMP方法不同，此時的懲罰因子p對子結(jié)構(gòu)剛度矩陣的懲罰，不是凝聚后的超單元對應(yīng)的剛度矩陣，而是子結(jié)構(gòu)內(nèi)所有單元對應(yīng)的有限元單元剛度矩陣。當(dāng)子結(jié)構(gòu)構(gòu)型優(yōu)化結(jié)束后，需對其進(jìn)行二次凝聚，從而在宏觀結(jié)構(gòu)中計算每個子結(jié)構(gòu)邊界的位移，實現(xiàn)宏微結(jié)構(gòu)的下一次迭代計算。

2.3 靈敏度計算及變量迭代

根據(jù)所建立的子結(jié)構(gòu)代理優(yōu)化模型，上述基于子結(jié)構(gòu)的周期性結(jié)構(gòu)設(shè)計問題可重新定義為

式中：[κ]為拉格朗日因子，該因子可通過分半算法強(qiáng)制使設(shè)計域中的材料滿足體積約束;體積導(dǎo)數(shù)[?V?ρj=vj]為每個子結(jié)構(gòu)的體積約束，在此與子結(jié)構(gòu)當(dāng)前的材料體積相等。

在優(yōu)化迭代計算中，通過OC方法進(jìn)行變量更新，每個子結(jié)構(gòu)的初始體積約束均相等，整個設(shè)計域在整個優(yōu)化過程中的材料體積與初始設(shè)置的體積約束相等。即材料體積始終不變，優(yōu)化過程只在周期性子結(jié)構(gòu)中尋找材料的最優(yōu)分配比例，達(dá)到最優(yōu)設(shè)計。對應(yīng)優(yōu)化框架如圖3所示。

3 數(shù)值計算實例

3.1 雙端固支梁結(jié)構(gòu)設(shè)計

以雙端固支梁的拓?fù)鋬?yōu)化設(shè)計為例[31]，設(shè)計域大小為[L1×L2=2×1]，其左右兩側(cè)固定，設(shè)計域中央承受[F=1]的垂直向壓力如圖4所示。所選用材料的彈性模量為[E0=1]，柏松比[υ0=0.3]。

把宏觀設(shè)計域劃分為3種不同數(shù)量的子結(jié)構(gòu)：[Nx×Ny=2×1]，[Nx×Ny=4×2]和[Nx×Ny=8×4];子結(jié)構(gòu)分別設(shè)置為[30×30]，[50×50]兩類不同單元數(shù)的單元集合。根據(jù)宏觀模型的定義，在優(yōu)化過程中，設(shè)體積約束設(shè)置為0.3，由于宏觀結(jié)構(gòu)中沒有設(shè)計變量，因此每個周期性子結(jié)構(gòu)的體積約束也為0.3。在周期性子結(jié)構(gòu)的優(yōu)化中，當(dāng)懲罰因子設(shè)置為1，所優(yōu)化的周期性子結(jié)構(gòu)拓?fù)浣Y(jié)構(gòu)中存在大量的中間密度單元，沒有明顯的桿系結(jié)構(gòu)特征。為了得到更為精細(xì)的局部結(jié)構(gòu)特征，在此把懲罰因子設(shè)置[p=3]。優(yōu)化迭代步長設(shè)置為0.02，為避免周期性子結(jié)構(gòu)中出現(xiàn)棋盤格現(xiàn)象，靈敏度的過濾半徑設(shè)置為子結(jié)構(gòu)內(nèi)部單元長度的1.5倍，當(dāng)相鄰兩次目標(biāo)函數(shù)值相對誤差小于0.000 1時所得拓?fù)錁?gòu)型即為最優(yōu)結(jié)構(gòu)。其優(yōu)化結(jié)果如表1所示。

從優(yōu)化結(jié)果中，可推斷出以下結(jié)論：1）相同的宏觀結(jié)構(gòu)子結(jié)構(gòu)劃分，優(yōu)化時所選取的子結(jié)構(gòu)大小不同，其最終的優(yōu)化拓?fù)錁?gòu)形也不同，其宏觀結(jié)構(gòu)的柔度值隨著子結(jié)構(gòu)單元數(shù)增多而減小;2）在同一大小的子結(jié)構(gòu)下，宏觀結(jié)構(gòu)子結(jié)構(gòu)劃分不同，最終優(yōu)化的周期性子結(jié)構(gòu)構(gòu)形基本相似，其宏觀結(jié)構(gòu)的柔度值隨著子結(jié)構(gòu)數(shù)量增多而逐漸增大。當(dāng)子結(jié)構(gòu)所含單元數(shù)增多時，意味著設(shè)計域的有限元網(wǎng)格更為細(xì)化，所得到的宏觀優(yōu)化結(jié)構(gòu)特征越精細(xì)，其優(yōu)化結(jié)構(gòu)的柔度值越小。

從優(yōu)化結(jié)果看，微結(jié)構(gòu)邊界間具有良好的材料連通性，改變了文獻(xiàn)[18]中由于應(yīng)變場突變引起的微結(jié)構(gòu)間材料不連通的情況，為其結(jié)構(gòu)加工制造提供了模型基礎(chǔ)。同時，設(shè)計過程中宏觀結(jié)構(gòu)與微觀構(gòu)型間具有明確的尺寸比例關(guān)系，保證了設(shè)計出的結(jié)構(gòu)性能與實際制造出的結(jié)構(gòu)性能的一致性，避免了均勻化方法中因尺度分離引起的結(jié)構(gòu)設(shè)計與制造的性能誤差。

3.2 懸臂梁結(jié)構(gòu)設(shè)計

在懸臂梁的拓?fù)鋬?yōu)化設(shè)計中，設(shè)置設(shè)計域大小為[L1×L2=2×1]，其左側(cè)固定，設(shè)計域右上方角點(diǎn)承受[F=1]的垂直向壓力如圖5所示。所選用材料的彈性模量為[E0=1]，柏松比[υ0=0.3]。

把宏觀設(shè)計域劃分為3種不同數(shù)量的子結(jié)構(gòu)：[Nx×Ny=2×1]，[Nx×Ny=4×2]和[Nx×Ny=8×4];子結(jié)構(gòu)設(shè)置為[40×40]的單元集合。為驗證懲罰因子對優(yōu)化構(gòu)型的影響，分別設(shè)置懲罰因子[p=1]和[p=2]。靈敏度的過濾半徑設(shè)置為子結(jié)構(gòu)內(nèi)部單元長度的1.2倍，優(yōu)化迭代步長設(shè)置為0.02，當(dāng)相鄰兩次目標(biāo)函數(shù)值相對誤差小于10-4時所得拓?fù)錁?gòu)型即為最優(yōu)結(jié)構(gòu)。其優(yōu)化結(jié)構(gòu)如表2所示。

從優(yōu)化結(jié)果中可以看出：當(dāng)懲罰因子設(shè)置為1，所優(yōu)化的周期性子結(jié)構(gòu)拓?fù)浣Y(jié)構(gòu)中存在大量的中間密度單元，沒有明顯的桿系結(jié)構(gòu)特征，隨著懲罰因子的增大，其桿系特征逐漸明顯。同時，最優(yōu)宏觀構(gòu)型的柔度值隨著懲罰因子和宏觀子結(jié)構(gòu)劃分?jǐn)?shù)的增大而增大。從優(yōu)化結(jié)構(gòu)來看，所得到的宏觀結(jié)構(gòu)及其對應(yīng)的微觀構(gòu)型，其材料具有很好的連通性，因此所得構(gòu)型具有良好的可制造性。

需要注意的是由于SIMP方法中懲罰因子的影響，使得子結(jié)構(gòu)內(nèi)部有限元單元的相對密度不是離散的0或1。盡管通過子結(jié)構(gòu)方法可獲得具有良好材料連通性的優(yōu)化結(jié)構(gòu)，但是所產(chǎn)生的中間密度單元給制造帶來了困難。

4 結(jié)論

本文提出了一種基于子結(jié)構(gòu)的周期性微結(jié)構(gòu)拓?fù)鋬?yōu)化方法，解決了傳統(tǒng)均勻化方法中尺度分離引起的微結(jié)構(gòu)間材料不連通和結(jié)構(gòu)設(shè)計與制造間性能誤差的問題，為設(shè)計制造一體化提供了一種參考方法。本文研究結(jié)論可歸納如下。

1）通過子結(jié)構(gòu)凝聚構(gòu)建宏觀結(jié)構(gòu)的超單元計算模型，把宏觀結(jié)構(gòu)尺寸信息帶入到了微結(jié)構(gòu)中，同時，通過引入超單元間的共同部分邊界確保了微結(jié)構(gòu)設(shè)計過程中材料的連通性。

2）優(yōu)化計算結(jié)果表明，宏觀子結(jié)構(gòu)劃分與子結(jié)構(gòu)內(nèi)部有限元網(wǎng)格越多，其優(yōu)化構(gòu)型的桿系特征越明顯。由于采用的變密度法SIMP，當(dāng)懲罰因子取值較小時，得到大量的中間密度單元。

3）采用SIMP方法，不可避免存在中間密度單元，給實際加工制造帶來困難。同時，所設(shè)計的宏微結(jié)構(gòu)模型邊界存在臺階現(xiàn)象，因此需要在可加工性方面做進(jìn)一步的改進(jìn)。

參考文獻(xiàn)：

[1] ? ?ZHU J，ZHANG W，BECKERS P. Integrated layout design of multi-component system[J]. International journal for numerical methods in engineering，2009，78（6）：631-651.

[2] ? ?LIU S，HU R，LI Q，et al. Topology optimization-based lightweight primary mirror design of a large-aperture space telescope[J]. Applied optics，2014，53（35）：8318-8325.

[18] ?XIA L，BREITKOPF P. Multiscale structural topology optimization with an approximate constitutive model for local material microstructure[J]. Computer Methods in Applied Mechanics and Engineering，2015，286：147-167.

[19] ?SIVAPURAM R，DUNNING P D，Kim H A. Simultaneous material and structural optimization by multiscale topology optimization[J]. Structural and multidisciplinary optimization，2016，54（5）：1267-1281.

[20] ?LI H，LUO Z，GAO L，et al. Topology optimization for functionally graded cellular composites with metamaterials by level sets[J]. Computer Methods in Applied Mechanics and Engineering，2018，328：340-364.

[21] ?ZHU Y，LI S，DU Z，et al. A novel asymptotic-analysis-based homogenisation approach towards fast design of infill graded microstructures[J]. Journal of the Mechanics and Physics of Solids，2019，124：612-633.

[22] ?LI Q，XU R，LIU J，et al. Topology optimization design of multi-scale structures with alterable microstructural length-width ratios[J]. Composite Structures，2019，230：111454.

[23] ?XIE Y M，ZUO Z H，Huang X，et al. Convergence of topological patterns of optimal periodic structures under multiple scales[J]. Structural and Multidisciplinary Optimization，2012，46（1）：41-50.

[24] ?XIA L，XIA Q，HUANG X，et al. Bi-directional evolutionary structural optimization on advanced structures and materials：a comprehensive review[J]. Archives of Computational Methods in Engineering，2018，25（2）：437-478.

[25] ?FU J，XIA L，GAO L，et al. Topology optimization of periodic structures with substructuring[J]. Journal of Mechanical Design，2019，141（7）：071403.

[26] ?HOU Y，ZHAO Q，Sapanathan T，et al. Parameter identifiability of ductile fracture criterion for DP steels using bi-level reduced surrogate model[J]. Engineering Failure Analysis，2019，100：300-311.

[27] ?BENDS?E M P，SIGMUND O. Material interpolation schemes in topology optimization[J]. Archive of applied mechanics，1999，69（9-10）：635-654.

[28] ?SIGMUND O. A 99 line topology optimization code written in Matlab[J]. Structural and multidisciplinary optimization，2001，21（2）：120-127.

[29] ?WU Z，XIA L，Wang S，et al. Topology optimization of hierarchical lattice structures with substructuring[J]. Computer Methods in Applied Mechanics and Engineering，2019，345：602-617.

[30] ?BENDS?E M P，SIGMUND O. Topology optimization：theory，methods and applications[M]. Berlin：Springer-Verlag，2003.

[31] ?ALEXANDERSEN J，LAZAROV B S. Topology optimization of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner[J]. Computer Methods in Applied Mechanics and Engineering，2015，290：156-182.

[責(zé)任編輯 ? ?楊 ? ?屹]