杜義賢，羅明亮，付君健,2，田啟華，周祥曼，孫鵬飛

(1.三峽大學(xué)機(jī)械與動(dòng)力學(xué)院，湖北 宜昌 443002；2.水電機(jī)械設(shè)備設(shè)計(jì)與維護(hù)湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 宜昌 443002)

0 引言

周期性多孔結(jié)構(gòu)是功構(gòu)一體化的優(yōu)良載體，具有高比表面積、高比剛度、高比強(qiáng)度、隔熱、隔音等特性[1-2]。在電化學(xué)和生物工程領(lǐng)域，周期性多孔結(jié)構(gòu)比表面積大小對(duì)結(jié)構(gòu)性能具有重要影響。例如，在鋰電池領(lǐng)域，電極中含有多孔式的集流體結(jié)構(gòu)，其表面包覆有活性材料，集流體結(jié)構(gòu)比表面積大小決定活性材料的分布[3]，對(duì)鋰離子擴(kuò)散系數(shù)、電子導(dǎo)電率、鋰離子存儲(chǔ)空間有重要影響。在化學(xué)工程領(lǐng)域，為了增大催化劑與反應(yīng)物的接觸面積，多孔催化劑載體結(jié)構(gòu)需要有較高的比表面積[4]。在生物工程領(lǐng)域，較高比表面積和合適孔徑大小的骨支架有利于細(xì)胞在其中生長[5]。因此，比表面積大小是周期性多孔結(jié)構(gòu)重要的工程設(shè)計(jì)參數(shù)。拓?fù)鋬?yōu)化通過調(diào)控結(jié)構(gòu)的參數(shù)，將結(jié)構(gòu)參數(shù)與性能聯(lián)系起來[6]，從而獲得高比表面積的拓?fù)錁?gòu)型。

拓?fù)鋬?yōu)化是周期性多孔結(jié)構(gòu)的重要設(shè)計(jì)方法[7]。通過施加周期性邊界條件約束，拓?fù)鋬?yōu)化可在設(shè)計(jì)域內(nèi)尋找滿足目標(biāo)函數(shù)和約束條件的最優(yōu)材料分布[8]。提高周期性多孔結(jié)構(gòu)比表面積的拓?fù)鋬?yōu)化設(shè)計(jì)方法可分為兩類：一類是基于均勻化[9-10]的設(shè)計(jì)方法；另一類是基于宏觀周期性約束的設(shè)計(jì)方法[11]。例如，采用均勻化法計(jì)算復(fù)合材料的性能參數(shù)，實(shí)現(xiàn)單胞內(nèi)材料的重新分配[12]，有效提高了周期性多孔結(jié)構(gòu)的比表面積；在單一體積約束下，基于能量均勻化法，對(duì)結(jié)構(gòu)等效屬性評(píng)價(jià)[13]，可得到具有高比表面積的周期性多孔結(jié)構(gòu)。但是，采用均勻化方法進(jìn)行高比表面積周期性多孔結(jié)構(gòu)設(shè)計(jì)，其尺度分離假設(shè)會(huì)帶來微結(jié)構(gòu)單胞間材料不連通的問題，不具備制造性。此外，將宏觀設(shè)計(jì)域分解為若干子區(qū)域，利用子結(jié)構(gòu)凝聚構(gòu)建超單元計(jì)算模型減少有限元計(jì)算量[14-15]，組裝子結(jié)構(gòu)可得到高比面積周期性多孔結(jié)構(gòu)。對(duì)宏觀周期性多孔結(jié)構(gòu)的最大尺寸進(jìn)行限制，可進(jìn)一步提高其比表面積。然而，基于宏觀周期性約束的高比表面積多孔結(jié)構(gòu)設(shè)計(jì)方法存在大量約束問題，不便于計(jì)算。

本文提出了一種高比表面積周期性多孔結(jié)構(gòu)拓?fù)鋬?yōu)化方法，引入局部體積約束，使設(shè)計(jì)域內(nèi)材料進(jìn)一步分散，顯著提高了結(jié)構(gòu)的比表面積。通過p-norm函數(shù)將多體積約束凝聚為單一體積約束，解決了宏觀周期性約束產(chǎn)生的大量約束問題，提高了拓?fù)鋬?yōu)化的求解效率。數(shù)值算例驗(yàn)證了本文方法的有效性。

1 周期性多孔結(jié)構(gòu)的定義

本文采用宏觀周期性約束方法，將整體設(shè)計(jì)域劃分為若干個(gè)相同大小的子區(qū)域，如圖1所示。所有子區(qū)域中相同位置單元的相對(duì)密度保持一致，從而使各子區(qū)域具有相同的拓?fù)湫问剑员ＷC結(jié)構(gòu)的周期性[11]。各單元密度關(guān)系的數(shù)學(xué)表達(dá)式為：

x1,j=x2,j=…=xm,j
xi,j∈[0,1](i=1,2,…,m;j=1,2,…,n)

(1)

式中，xi,j為設(shè)計(jì)變量，表示第i個(gè)子區(qū)域內(nèi)第j個(gè)單元的密度；m為總設(shè)計(jì)域劃分的子區(qū)域個(gè)數(shù)，n為子區(qū)域內(nèi)單元數(shù)量。

圖1 二維周期性結(jié)構(gòu)設(shè)計(jì)域

2 局部體積約束

假定拓?fù)鋬?yōu)化設(shè)計(jì)域內(nèi)材料分布由邏輯值ρe表示，ρe=1或0代表實(shí)體單元或孔洞單元。為了限制設(shè)計(jì)域內(nèi)材料積累形成大的實(shí)體區(qū)域，引入局部體積約束使設(shè)計(jì)域內(nèi)材料進(jìn)一步分散。

(2)

(3)

圖2 局部體積約束

(4)

(5)

式(5)可重新寫為：

(6)

式中，N是整體設(shè)計(jì)域內(nèi)單元總數(shù)量。p越大，每個(gè)單元的約束就越強(qiáng)，同時(shí)也增加了問題的非線性。

引入一個(gè)數(shù)值從0～1連續(xù)變化的單元密度xe作為拓?fù)鋬?yōu)化設(shè)計(jì)變量，為去除中間密度單元產(chǎn)生的棋盤格現(xiàn)象[17]，通過局部濾波器計(jì)算相鄰單元的加權(quán)平均值對(duì)xe過濾：

(7)

(8)

式中，re為敏度過濾半徑，其數(shù)值小于圖2中集合Ne的半徑RΩ；oi和oe為單元中心；式(7)中權(quán)重因子Wi,e大小和oi、oe兩單元中心距離有關(guān)：

(9)

(10)

式中，參數(shù)β控制閾值函數(shù)斜率，如圖3所示。當(dāng)β越大時(shí)，函數(shù)值越接近0或1。如果直接應(yīng)用一個(gè)較大β值，會(huì)導(dǎo)致高度非線性解；因此，本文從β=1開始迭代，經(jīng)過一定迭代次數(shù)后，將其數(shù)值翻倍，這個(gè)過程稱為參數(shù)擴(kuò)展，可提高優(yōu)化求解的收斂性[18]。

圖3 不同參數(shù)β控制的Heaviside函數(shù)

3 拓?fù)鋬?yōu)化數(shù)學(xué)模型

基于改進(jìn)固體各向同性(Modified Solid Isotropic Material with Penalization, modified SIMP)的變密度法材料插值模型[19]，建立單元楊氏模量Ee數(shù)學(xué)表達(dá)式：

(11)

式中，E0為實(shí)體單元?jiǎng)偠戎?；Emin為一非常小的數(shù)值，代表孔洞單元?jiǎng)偠戎?，以防止整體剛度矩陣產(chǎn)生奇異值；γ為單元密度的懲罰因子，對(duì)拓?fù)鋬?yōu)化中具有中間密度的單元進(jìn)行懲罰，使其收斂于指定的密度上下界，從而抑制灰度單元。任意單元?jiǎng)偠染仃嘖e為：

Ke=Ee(ρe)k0

(12)

式中，k0為實(shí)體單元?jiǎng)偠染仃嚒?/p>

引入局部體積約束對(duì)設(shè)計(jì)域內(nèi)的材料分布進(jìn)行限制，基于改進(jìn)SIMP插值模型，以單元相對(duì)密度xe為設(shè)計(jì)變量，構(gòu)造多孔結(jié)構(gòu)周期性約束條件，以結(jié)構(gòu)剛度最大化為優(yōu)化目標(biāo)，建立拓?fù)鋬?yōu)化數(shù)學(xué)模型：

find:xe= {x1,j,x2,j,x3,j,…,xi,j}

(i= 1,2,…,m;j= 1,2,…,n)

(13)

式中，c、K、U、F分別為整體柔度、剛度矩陣、位移向量、載荷向量；g(x)為局部體積約束方程，控制結(jié)構(gòu)比表面積大小。

4 敏度分析

基于式(13)的拓?fù)鋬?yōu)化數(shù)學(xué)模型，以單元相對(duì)密度xe為設(shè)計(jì)變量，采用有限元移動(dòng)漸近線方法[20](Method of Moving Asymptotes, MMA)來更新設(shè)計(jì)變量，需要分別計(jì)算目標(biāo)函數(shù)柔度c和局部體積約束方程g(x)對(duì)設(shè)計(jì)變量xe的一階導(dǎo)數(shù)，用鏈?zhǔn)椒▌t計(jì)算如下：

(14)

(15)

由式(11)、式(13)得：

(16)

式(14)、式(15)中其他部分導(dǎo)數(shù)為：

(17)

(18)

(19)

(20)

為保證式(13)中每個(gè)子區(qū)域第j個(gè)單元的密度相等，需要對(duì)敏度平均：

(21)

5 數(shù)值算例

本文采用懸臂梁和Michell梁拓?fù)鋬?yōu)化算例，分別計(jì)算和對(duì)比有、無局部體積約束的周期性多孔結(jié)構(gòu)表面積，驗(yàn)證本文提高結(jié)構(gòu)比表面積方法的有效性。在三維空間中，比表面積是指多孔結(jié)構(gòu)單位質(zhì)量所具有的孔洞總面積[21]，可定義為：面積/體積；在二維平面中，比表面積是指多孔結(jié)構(gòu)單位面積所具有的孔洞總周長，可定義為：長度/面積。材料彈性模量E=1，泊松比u=0.3，單元大小默認(rèn)為1(所有單位均為無量綱值)，以結(jié)構(gòu)剛度最大為目標(biāo)建立優(yōu)化模型，采用MMA算法進(jìn)行求解。

5.1 懸臂梁

圖4所示為懸臂梁結(jié)構(gòu)，其設(shè)計(jì)域豎直方向和水平方向的長度分別為：H=100，L=200。左側(cè)邊界固定，右側(cè)邊界中間節(jié)點(diǎn)加載一豎直向下的集中載荷F。有限元單元為四節(jié)點(diǎn)矩形單元，采用200×100網(wǎng)格離散設(shè)計(jì)域，劃分為4×2大小的子區(qū)域，總體體積分?jǐn)?shù)為0.45。

圖4 懸臂梁結(jié)構(gòu)的設(shè)計(jì)域示意圖

取局部體積分?jǐn)?shù)上限α=0.6，有、無局部體積約束的周期性多孔結(jié)構(gòu)對(duì)比如圖5所示，其迭代過程如圖6所示。

(a) 無局部體積約束 (b) 有局部體積約束圖5 懸臂梁拓?fù)鋬?yōu)化結(jié)構(gòu)對(duì)比圖

(a) 無局部體積約束 (b) 有局部體積約束圖6 懸臂梁拓?fù)鋬?yōu)化迭代曲線

由圖5可以得出，無局部體積約束的結(jié)構(gòu)桿件較為粗壯；有局部體積約束的結(jié)構(gòu)最大尺寸減小，拓?fù)鋬?yōu)化構(gòu)型細(xì)節(jié)特征增加，孔徑變小，孔洞數(shù)量增多。圖6為懸臂梁拓?fù)鋬?yōu)化迭代曲線，圖6a為無局部體積約束的迭代曲線，整個(gè)優(yōu)化過程收斂平穩(wěn)，目標(biāo)函數(shù)柔度c最終收斂于129；圖6b給出了局部體積分?jǐn)?shù)上限α=0.6、p-norm函數(shù)控制參數(shù)p為16時(shí)的迭代曲線，Heaviside函數(shù)控制參數(shù)β每迭代40次數(shù)值翻倍，因而目標(biāo)函數(shù)每隔40步會(huì)出現(xiàn)短暫的波動(dòng)，目標(biāo)函數(shù)柔度c最終收斂于173。因此，對(duì)比有、無局部體積約束的周期性結(jié)構(gòu)拓?fù)錁?gòu)型可以得出，有局部體積約束的結(jié)構(gòu)細(xì)桿特征增加，結(jié)構(gòu)柔度變大，實(shí)現(xiàn)了最大尺寸約束，具有更高的比表面積。

因周期性多孔結(jié)構(gòu)每個(gè)單胞完全一樣，為減少計(jì)算量，測(cè)量、對(duì)比其單胞的表面積即可。拓?fù)鋬?yōu)化原始構(gòu)型存在孔洞邊緣“模糊”、表面不光滑的問題，往往不能直接精確繪制其孔洞邊緣，從而無法準(zhǔn)確測(cè)量比表面積。為解決該問題，采用水平集方法(Level Set Method)[22-23]對(duì)周期性多孔結(jié)構(gòu)的原始單胞進(jìn)行后處理。以單胞拓?fù)錁?gòu)型圖的像素為單位測(cè)量孔洞總周長，即為表面積，結(jié)果見表1。在懸臂梁算例中，具有局部體積約束的拓?fù)錁?gòu)型比表面積(比表面積與表面積成正比)提高了約300%，驗(yàn)證了本文方法的有效性。

表1 懸臂梁單胞拓?fù)錁?gòu)型對(duì)比

5.2 Michell梁

圖7所示為Michell梁結(jié)構(gòu)，其設(shè)計(jì)域豎直方向和水平方向的長度分別為：H=100，L=200。左側(cè)邊界底部為簡支支撐，右側(cè)邊界底部為固定支撐，底部中間加載一豎直向下的集中載荷F。有限元單元為四節(jié)點(diǎn)矩形單元，采用200×100網(wǎng)格離散設(shè)計(jì)域，劃分為4×2大小的子區(qū)域，總體體積分?jǐn)?shù)為0.45。

圖7 Michell梁結(jié)構(gòu)的設(shè)計(jì)域示意圖

取局部材料體積分?jǐn)?shù)上限α=0.6，有、無局部體積約束的周期性多孔結(jié)構(gòu)對(duì)比如圖8所示，其迭代過程如圖9所示。

(a) 無局部體積約束 (b) 有局部體積約束圖8 Michell梁拓?fù)鋬?yōu)化結(jié)構(gòu)對(duì)比圖

(a) 無局部體積約束 (b) 有局部體積約束圖9 Michell梁拓?fù)鋬?yōu)化迭代曲線

由圖8可以得出，在Michell梁算例中，有局部體積約束的周期性拓?fù)錁?gòu)型具有更多的細(xì)節(jié)特征，孔洞數(shù)量增多。圖9 Michell梁拓?fù)鋬?yōu)化迭代曲線中，無局部體積約束的周期性結(jié)構(gòu)柔度c最終收斂于28；有局部體積約束的周期性結(jié)構(gòu)柔度c最終收斂于32。因此，有局部體積約束的結(jié)構(gòu)細(xì)桿特征增加，結(jié)構(gòu)柔度變大，實(shí)現(xiàn)了最大尺寸約束，具有更高的比表面積。單胞拓?fù)錁?gòu)型及表面積大小具體見表2。在Michell梁算例中，本文提出的局部體積約束方法將拓?fù)錁?gòu)型的比表面積提高了約200%，驗(yàn)證了該方法的有效性。

表2 Michell梁單胞拓?fù)錁?gòu)型對(duì)比

6 結(jié)論

本文提出了一種高比表面積周期性多孔結(jié)構(gòu)拓?fù)鋬?yōu)化方法，通過引入局部體積約束，使設(shè)計(jì)域內(nèi)材料進(jìn)一步分散，對(duì)周期性多孔結(jié)構(gòu)的最大尺寸實(shí)現(xiàn)了有效控制，顯著提高了結(jié)構(gòu)的比表面積。通過p-norm函數(shù)將多體積約束凝聚為單一體積約束，解決了局部體積約束產(chǎn)生的大量約束問題，提高了拓?fù)鋬?yōu)化的求解效率。使用水平集方法對(duì)周期性結(jié)構(gòu)單胞進(jìn)行后處理，得到光滑邊界的拓?fù)錁?gòu)型，從而在數(shù)據(jù)分析軟件中精確測(cè)量了結(jié)構(gòu)的表面積。拓?fù)鋬?yōu)化經(jīng)典二維數(shù)值算例驗(yàn)證了本文方法的有效性。