基于PDE靈敏度濾波器的算法研究*

孟換利，張岐良，王 杰

（1.中山大學(xué) 系統(tǒng)科學(xué)與工程學(xué)院，廣州 510006；2.天津航天長征火箭制造有限公司，天津 300462）

引言

為應(yīng)對越來越高的力學(xué)性能要求和越來越嚴(yán)酷的服役環(huán)境，在結(jié)構(gòu)設(shè)計中應(yīng)用拓?fù)鋬?yōu)化技術(shù)已成為一種必然選擇.在結(jié)構(gòu)的概念設(shè)計階段通過結(jié)構(gòu)拓?fù)鋬?yōu)化確定材料的最優(yōu)布局形式，可以幫助設(shè)計者發(fā)現(xiàn)結(jié)構(gòu)最佳的應(yīng)力傳導(dǎo)路徑，從而以最少的材料實現(xiàn)最優(yōu)的結(jié)構(gòu)性能.Deaton 和Grandhi[1]調(diào)查了2000年至2012年連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的研究現(xiàn)狀與應(yīng)用.Zhu 等[2]綜述了面向增材制造的拓?fù)鋬?yōu)化技術(shù)的典型應(yīng)用場景、研究熱點和挑戰(zhàn).衛(wèi)志軍等[3]對雙隔板進(jìn)行了拓?fù)鋬?yōu)化，可以有效地抑制流體晃蕩.

然而，連續(xù)體結(jié)構(gòu)的拓?fù)鋬?yōu)化結(jié)果往往存在棋盤格現(xiàn)象、數(shù)值不穩(wěn)定等問題，Sigmund[4]以及Bruns、Tortorelli[5]使用的靈敏度濾波器和密度濾波器可以有效地解決上述問題.研究者們在此基礎(chǔ)上發(fā)展了多種過濾方法，如灰度單元等效轉(zhuǎn)換方法[6]、動態(tài)參數(shù)調(diào)整方法[7]和灰度單元分層雙重懲罰方法[8]等.上述方法在一定程度上可以有效地抑制灰度單元、消除棋盤格現(xiàn)象，防止邊界灰度擴(kuò)散，但沒有關(guān)注過濾效率.與通過獲取鄰域的單元信息來修正中心單元的靈敏度濾波器不同，PDE 濾波器的實質(zhì)是將靈敏度過濾值表示成Neumann 邊界條件下的Helmholtz 方程的解，其求解只需考慮有限元離散的網(wǎng)格信息，避免了占用大量內(nèi)存，極大地提高了過濾效率.然而隨著網(wǎng)格的不斷細(xì)分，該濾波器的代數(shù)方程維數(shù)越來越高，采用以Gauss 消元法為基礎(chǔ)的直接法求解，會導(dǎo)致耗時過長和內(nèi)存需求量增加等問題，從而影響該濾波器的效率.此外，隨著計算機(jī)性能的提高和算法技術(shù)的發(fā)展，精細(xì)的數(shù)值模擬成為可能，使得待求解的線性方程組的維數(shù)越來越高，研發(fā)高效的線性求解器加速算法[9-11]已變得至關(guān)重要.

多重網(wǎng)格是一種高效的求解算法，近期在偏微分方程參數(shù)反演問題方面有許多應(yīng)用[12-13]，并廣泛地應(yīng)用到圖像處理[14]、流體模擬[15]和結(jié)構(gòu)拓?fù)鋬?yōu)化[16-17]等領(lǐng)域.為提高大規(guī)模結(jié)構(gòu)拓?fù)鋬?yōu)化靈敏度過濾的效率，本文基于體積約束下柔度值最小的固體各向同性懲罰微結(jié)構(gòu)(SIMP)插值模型，針對大規(guī)模PDE 靈敏度濾波器的求解問題，對其進(jìn)行有限元分析得到代數(shù)方程，再采用共軛梯度算法、多重網(wǎng)格算法以及多重網(wǎng)格預(yù)處理共軛梯度算法進(jìn)行求解，得到了過濾域的節(jié)點靈敏度值，通過節(jié)點到單元的映射得到了單元靈敏度過濾值.

1 優(yōu)化模型

在變密度法中，SIMP 模型是目前工程領(lǐng)域應(yīng)用最廣泛的材料插值模型，為了使模型更具有代表性，基于SIMP 插值模型，建立在體積約束下柔度值最小的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化模型：

其中，目標(biāo)函數(shù)C為結(jié)構(gòu)的總體柔度，F(xiàn)為力向量，U為位移，K為結(jié)構(gòu)的總剛度矩陣，xi和ui是第i個單元的相對密度和節(jié)點位移，ki和k0是單元剛度矩陣，V是優(yōu)化后的結(jié)構(gòu)體積，V?為結(jié)構(gòu)體積約束，V0為整個設(shè)計域的初始體積，f為優(yōu)化體積比，vi為第i個單元的體積，xmin和xmax為單元相對密度的最小和最大極限值，N為結(jié)構(gòu)離散單元總數(shù).

SIMP 法又稱冪律法，其材料插值模型以密度為設(shè)計變量，單元材料彈性模量和密度之間的關(guān)系為

其中E0和Emin是固體和空洞部分的材料彈性模量，Ei為第i個單元插值后的彈性模量.

2 拓?fù)鋬?yōu)化算法

2.1 總體框架

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化實現(xiàn)流程如圖1所示.首先，初始化優(yōu)化參數(shù)：定義問題的設(shè)計區(qū)域、約束和載荷.其次，建立有限元模型：劃分有限元網(wǎng)格，將結(jié)構(gòu)離散化，建立式(1)所示的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化數(shù)學(xué)模型.另外，有限元分析：計算剛度矩陣、節(jié)點位移、目標(biāo)柔度值和靈敏度.之后，PDE 靈敏度過濾：解決數(shù)值不穩(wěn)定問題和消除棋盤格現(xiàn)象.然后，采用二分法和OC 準(zhǔn)則更新Lagrange 乘子和設(shè)計變量，判斷是否滿足體積約束條件：滿足則輸出設(shè)計變量；若不滿足，則回到上一步再次更新Lagrange 乘子和設(shè)計變量.最后，檢查是否滿足收斂準(zhǔn)則：若滿足，則輸出拓?fù)鋬?yōu)化結(jié)果；若不滿足，則再次進(jìn)行有限元分析，進(jìn)行下一輪迭代.

圖1 基于SIMP 插值模型連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化實現(xiàn)流程圖Fig.1 The flow chart of continuum structure topology optimization based on the SIMP interpolation model

為解決傳統(tǒng)PDE 靈敏度濾波器過濾慢的問題，本文采用共軛梯度算法求解該濾波器的對稱正定代數(shù)方程，為進(jìn)一步解決網(wǎng)格數(shù)量增加導(dǎo)致系數(shù)矩陣維數(shù)增加和趨近最優(yōu)解時，該算法收斂速度減緩的問題，結(jié)合多重網(wǎng)格可以均勻衰減誤差，加快算法收斂速度的優(yōu)點，使用多重網(wǎng)格算法和預(yù)處理共軛梯度算法求解該方程，提升過濾效率.此外，在預(yù)處理共軛梯度算法中，采用多重網(wǎng)格預(yù)處理降低系數(shù)矩陣的條件數(shù).

2.2 PDE 靈敏度濾波器

本文采用PDE 靈敏度濾波器[18-19]，其核心思想是將過濾域的節(jié)點靈敏度值隱式表示為具有Neumann 邊界條件的Helmholtz 方程的隱式解，然后通過節(jié)點到單元的映射得到單元靈敏度過濾結(jié)果：

2.3 PDE 靈敏度濾波器的求解

PDE 靈敏度濾波器的實質(zhì)是以Neumann 邊界為條件的Helmholtz 偏微分方程，其求解步驟如下：

第1 步 有限元離散得到Helmholtz 方程(3)的代數(shù)方程：

其中KF是標(biāo)量問題的標(biāo)準(zhǔn)有限元剛度矩陣，TF是將單元相對密度映射到具有節(jié)點值的向量的矩陣.針對代數(shù)方程(5)的系數(shù)矩陣是對稱正定稀疏病態(tài)的特點，本文擬采用共軛梯度算法、多重網(wǎng)格算法和多重網(wǎng)格預(yù)處理共軛梯度算法求解該代數(shù)方程，從而得到過濾區(qū)域內(nèi)的節(jié)點靈敏度值.

第2 步 通過對過濾域的節(jié)點靈敏度值映射得到單元靈敏度過濾值：

2.3.1 多重網(wǎng)格算法

在固定網(wǎng)格中采用常規(guī)的迭代方法求解線性方程組，初始時收斂速度很快，趨近最優(yōu)解時，收斂速度慢，計算時間長.采用Fourier 級數(shù)表示誤差，得到兩種不同的誤差：波長較短、波形頻率高的高頻分量誤差和波長較長、波形頻率低的低頻分量誤差.由于波形頻率的高低是相對的，把細(xì)網(wǎng)格中的低頻波形放到粗網(wǎng)格中就有可能變成高頻波形.因此采用多重網(wǎng)格來消除誤差分量，使各種頻率的誤差得到比較均勻的衰減，從而加快迭代收斂速度.假設(shè)有L重網(wǎng)格，k=1 為最細(xì)的網(wǎng)格層，在k層上的多重網(wǎng)格迭代解，多重網(wǎng)格V 循環(huán)具體過程如下（偽代碼見算法1）：

步驟1 前光滑.采用常規(guī)的迭代方法消除誤差高頻分量，具體步驟為：在細(xì)網(wǎng)格上，采取2 ～ 3 次的阻尼Jacobi 迭代方法得到線性方程組的近似解.迭代公式的通式為

其中Sk為光滑子，Sk=ωDk，ω 為阻尼因子，Dk為Ak的主對角元素矩陣.

步驟2 粗網(wǎng)格校正.計算在粗網(wǎng)格上誤差的近似解，具體步驟為：首先計算細(xì)網(wǎng)格上的殘差然后使用限制算子把細(xì)網(wǎng)格上的殘差限制到粗網(wǎng)格上，得到粗網(wǎng)格上的殘差最后根據(jù)網(wǎng)格的層數(shù)求解線性方程組得到粗網(wǎng)格上誤差的近似解，若是最粗的網(wǎng)格層，則利用直接法求解否則利用多重網(wǎng)格V 循環(huán)求解其中粗網(wǎng)格算子Ak+1取Galerkin 算子.

步驟3 插值.計算細(xì)網(wǎng)格上的解，具體步驟為：先使用插值算子把粗網(wǎng)格上的近似誤差插值到細(xì)網(wǎng)格上然后使用公式得到在細(xì)網(wǎng)格的解.本文采用的插值算子為

步驟4 后光滑.為保證整個迭代對稱，在細(xì)網(wǎng)格上，采取具有相同參數(shù)的阻尼Jacobi 迭代法求解線性方程組.

算法1

2.3.2 多重網(wǎng)格預(yù)處理共軛梯度算法

考慮到代數(shù)方程的系數(shù)矩陣是大型、稀疏、對稱正定的，有時甚至為病態(tài)的，則重點研究了多重網(wǎng)格預(yù)處理共軛梯度算法（偽代碼見算法2），采用多重網(wǎng)格作為預(yù)處理方法降低系數(shù)矩陣的條件數(shù).

算法2

2.4 設(shè)計變量的更新

優(yōu)化準(zhǔn)則法（OC 準(zhǔn)則）具有收斂速度快、易于實現(xiàn)等特點.本文采取OC 準(zhǔn)則求解優(yōu)化模型(1)，基本思路為，先引入Lagrange 乘子把約束問題(1)轉(zhuǎn)化為無約束問題，再根據(jù)極值點存在的條件即K-T 條件以及不動點迭代得到設(shè)計變量的迭代公式為

2.5 收斂準(zhǔn)則

包括外循環(huán)和內(nèi)循環(huán)，外循環(huán)包括有限元分析和更新單元靈敏度，迭代收斂條件如下：

內(nèi)循環(huán)包括更新設(shè)計變量、Lagrange 乘子和過濾域的節(jié)點值.采用二分法和OC 準(zhǔn)則更新Lagrange 乘子和設(shè)計變量，以滿足模型(1)中的體積約束和設(shè)計變量的邊界約束，迭代收斂條件如下：

其中η為收斂控制因子，l2和l1為較大和較小的Lagrange 乘子.采用共軛梯度算法、多重網(wǎng)格算法以及多重網(wǎng)格預(yù)處理共軛梯度算法求解代數(shù)方程，構(gòu)造迭代收斂準(zhǔn)則為

其中τ為收斂控制因子，rk和r0為第k次和初始時代數(shù)方程的殘量.

3 數(shù)值算例

為了驗證本文拓?fù)鋬?yōu)化方法的可行性、可靠性和高效性，針對長960 mm、高320 mm 的懸臂梁和長960 mm、高480 mm Michell 結(jié)構(gòu)的棋盤格現(xiàn)象、數(shù)值不穩(wěn)定、靈敏度過濾慢等問題開展了數(shù)值實驗，如圖2 和3所示（算例優(yōu)化參數(shù)為體積約束比為0.5、彈性模量為1 GPa、Poisson 比為0.3、懲罰因子為3、移動極限為0.2、外載荷F為1 kN）.拓?fù)鋬?yōu)化程序基于文獻(xiàn)[16,20]的MATLAB 代碼進(jìn)行編寫，在PDE 濾波器的求解中引入了共軛梯度（CG）算法、多重網(wǎng)格（MG）算法以及多重網(wǎng)格預(yù)處理共軛梯度（MGCG）算法，并比較了三種算法對結(jié)構(gòu)拓?fù)鋬?yōu)化效率的影響，相關(guān)數(shù)值結(jié)果如下.

圖2 懸臂梁Fig.2 The cantilever beam

圖3 Michell 結(jié)構(gòu)Fig.3 The Michell structure

3.1 精度的影響

依次取4 種不同的誤差精度P(P=1E–2,1E–3,1E–4,1E–5)，采用CG 算法、網(wǎng)格層數(shù)為5 的MG 算法和MGCG 算法求解經(jīng)有限元分析Helmholtz 方程得到的代數(shù)方程，并對整個結(jié)構(gòu)分別劃分為480×160和480×240個單元的懸臂梁和Michell 結(jié)構(gòu)進(jìn)行半徑為6 的PDE 靈敏度過濾拓?fù)鋬?yōu)化實驗.其中，不同精度下懸臂梁和Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)、目標(biāo)柔度值C、整體迭代次數(shù)I以及整體優(yōu)化時間T的分布如圖4、圖5、表1、圖6、圖7所示.

圖6 懸臂梁結(jié)構(gòu)優(yōu)化時間分布圖Fig.6 The optimized time distribution of the cantilever beam

圖7 Michell 結(jié)構(gòu)優(yōu)化時間分布圖Fig.7 The optimized time distribution of the Michell structure

表1 兩種結(jié)構(gòu)的拓?fù)鋬?yōu)化迭代次數(shù)和柔度值Table 1 Total numbers of iterations and compliance for topology optimization of 2 structures

圖4 不同精度下CG、MG 和MGCG 算法的懸臂梁的優(yōu)化結(jié)構(gòu)：(a)P=1E?2 (CG)；(b)P=1E?2 (MG)；(c)P=1E?2 (MGCG)；(d)P=1E?3 (CG)；(e)P=1E?3(MG)；(f)P=1E?3 (MGCG)；(g)P=1E?4 (CG)；(h)P=1E?4 (MG)；(i)P=1E?4 (MGCG)；(j)P=1E?5 (CG)；(k)P=1E?5 (MG)；(l)P=1E?5 (MGCG)Fig.4 The optimized cantilever beams by CG,MG and MGCG algorithms with different degrees of precision:(a)P=1E?2 (CG);(b)P=1E?2 (MG);(c)P=1E?2 (MGCG);(d)P=1E?3 (CG);(e)P=1E?3 (MG);(f)P=1E?3 (MGCG);(g)P=1E?4 (CG);(h)P=1E?4 (MG);(i)P=1E?4 (MGCG);(j)P=1E?5 (CG);(k)P=1E?5 (MG);(l)P=1E?5 (MGCG)

圖5 不同精度下CG、MG 和MGCG 算法的Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)：(a)P=1E?2 (CG)；(b)P=1E?2 (MG)；(c)P=1E?2 (MGCG)；(d)P=1E?3 (CG)；(e)P=1E?3(MG)；(f)P=1E?3 (MGCG)；(g)P=1E?4 (CG)；(h)P=1E?4 (MG)；(i)P=1E?4 (MGCG)；(j)P=1E?5 (CG)；(k)P=1E?5 (MG)；(l)P=1E?5 (MGCG)Fig.5 The optimized Michell structures by CG,MG and MGCG algorithms with different degrees of precision:(a)P=1E?2 (CG);(b)P=1E?2 (MG);(c)P=1E?2 (MGCG);(d)P=1E?3 (CG);(e)P=1E?3 (MG);(f)P=1E?3 (MGCG);(g)P=1E?4 (CG);(h)P=1E?4 (MG);(i)P=1E?4 (MGCG);(j)P=1E?5 (CG);(k)P=1E?5 (MG);(l)P=1E?5 (MGCG)

從圖4、圖5 和表1 可以看出，在過濾半徑為6 的PDE 靈敏度濾波器中使用不同精度的CG 算法、MG 算法和MGCG 算法得到的懸臂梁和Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)圖基本相同，并且目標(biāo)函數(shù)值非常接近，至多相差1%，說明PDE 濾波器求解算法的精度對拓?fù)鋬?yōu)化的結(jié)果影響不大.此外，在PDE 濾波器中使用MGCG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的迭代次數(shù)整體低于其他兩種算法.

從圖6 和圖7 可以看出，同一精度下，在PDE 濾波器中使用MGCG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的時間整體低于CG 算法和MGCG 算法，且隨著精度的增加，在PDE 濾波器中使用CG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的時間逐漸降低并趨于穩(wěn)定，在PDE 濾波器中使用MG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的時間逐漸增加并趨于穩(wěn)定.此外，在PDE 濾波器中采用MGCG 算法對懸臂梁和Michell 結(jié)構(gòu)拓?fù)鋬?yōu)化的時間分別比CG 算法少了11% ～ 34%和65% ～ 80%，比MG 算法少了29% ～ 46%，3%和77%.因此，從結(jié)構(gòu)的整體優(yōu)化效率看，MGCG 算法高于CG 算法和MG 算法.

3.2 過濾半徑的影響

依次取4 種不同的過濾半徑R(R=3，6，9，12),采用CG 算法、網(wǎng)格層數(shù)為5 的MG 算法以及MGCG 算法求解有限元分析Helmholtz 方程得到的代數(shù)方程，并對整個結(jié)構(gòu)分別劃分為480×160和480×240個單元的懸臂梁和Michell 結(jié)構(gòu)進(jìn)行控制精度為1E–4 的PDE 靈敏度過濾的拓?fù)鋬?yōu)化實驗.其中，不同過濾半徑下懸臂梁和Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)、目標(biāo)柔度值C、整體迭代次數(shù)I以及整體優(yōu)化時間T的分布分別如圖8、圖9、表2、圖10、圖11所示.

圖8 不同過濾半徑下CG、MG 和MGCG 算法的懸臂梁的優(yōu)化結(jié)構(gòu)：(a)R=3 (CG)；(b)R=3 (MG)；(c)R=3 (MGCG)；(d)R=6 (CG)；(e)R=6 (MG)；(f)R=6 (MGCG)；(g)R=9 (CG)；(h)R=9 (MG)；(i)R=9 (MGCG)；(j)R=12 (CG)；(k)R=12 (MG)；(l)R=12 (MGCG)Fig.8 The optimized cantilever beams by CG,MG and MGCG algorithms with different filter radius:(a)R=3 (CG);(b)R=3 (MG);(c)R=3 (MGCG);(d)R=6 (CG);(e)R=6 (MG);(f)R=6 (MGCG);(g)R=9 (CG);(h)R=9 (MG);(i)R=9 (MGCG);(j)R=12 (CG);(k)R=12 (MG);(l)R=12 (MGCG)

圖9 不同過濾半徑下CG、MG 和MGCG 算法的Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)：(a)R=3 (CG)；(b)R=3 (MG)；(c)R=3 (MGCG)；(d)R=6 (CG)；(e)R=6 (MG)；(f)R=6 (MGCG)；(g)R=9 (CG)；(h)R=9 (MG)；(i)R=9 (MGCG)；(j)R=12 (CG)；(k)R=12 (MG)；(l)R=12 (MGCG)Fig.9 The optimized Michell structures by CG,MG and MGCG algorithms with different filter radius:(a)R=3 (CG);(b)R=3 (MG);(c)R=3 (MGCG);(d)R=6 (CG);(e)R=6 (MG);(f)R=6 (MGCG);(g)R=9 (CG);(h)R=9 (MG);(i)R=9 (MGCG);(j)R=12 (CG);(k)R=12 (MG);(l)R=12 (MGCG)

圖10 懸臂梁結(jié)構(gòu)優(yōu)化時間分布圖(不同過濾半徑)Fig.10 The optimized time distribution of the cantilever beam(different filter radius)

圖11 Michell 結(jié)構(gòu)優(yōu)化時間分布圖(不同過濾半徑)Fig.11 The optimized time distribution of the Michell structure(different filter radius)

從圖8、圖9 和表2 可以看出，在PDE 靈敏度濾波器中使用精度為1E–4 的CG 算法、MG 算法和MGCG 算法得到的懸臂梁和Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)圖、目標(biāo)函數(shù)值的特征隨著過濾半徑的變化呈現(xiàn)類似的變化規(guī)律.隨著過濾半徑的增大，拓?fù)鋬?yōu)化結(jié)構(gòu)的清晰度越來越低，細(xì)支結(jié)構(gòu)越來越少，說明該兩種算法對結(jié)構(gòu)承力路徑的識別精度與過濾半徑間體現(xiàn)出負(fù)相關(guān)性，且柔度值越來越大.且從表2 可以看出，隨著過濾半徑的改變，在PDE 濾波器中使用MG 算法和MGCG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的迭代次數(shù)整體低于CG 算法.

表2 兩種結(jié)構(gòu)的拓?fù)鋬?yōu)化迭代次數(shù)和柔度值(不同過濾半徑)Table 2 Total numbers of iterations and compliance for topology optimization of 2 structures (different filter radius)

從圖10 和圖11 可以看出，同一過濾半徑下，在PDE 濾波器中使用MG 算法和MGCG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的時間整體低于CG 算法，且隨著過濾半徑增大，MGCG 算法的優(yōu)勢逐漸降低.此外，在PDE 濾波器中采用MGCG 算法對懸臂梁和Michell 結(jié)構(gòu)拓?fù)鋬?yōu)化的時間分別比CG 算法少了11%～89%和50%～91%，比MG 算法少了7%～61%，7%和70%.因此，從結(jié)構(gòu)的整體優(yōu)化效率看，MGCG 算法高于CG 算法和MG 算法.

3.3 網(wǎng)格數(shù)量的影響

分別對網(wǎng)格劃分為960×320，過濾半徑為12，層數(shù)為6；網(wǎng)格劃分為480×160，過濾半徑為6，層數(shù)為5；網(wǎng)格劃分為240×80，過濾半徑為3，層數(shù)為4 的懸臂梁和網(wǎng)格劃分為960×480，過濾半徑為12，層數(shù)為6；網(wǎng)格劃分為480×240，過濾半徑為6，層數(shù)為5；網(wǎng)格劃分為240×120，過濾半徑為3，層數(shù)為4 的Michell 結(jié)構(gòu)進(jìn)行控制精度為1E–4 的PDE 靈敏度過濾的拓?fù)鋬?yōu)化實驗.其中，不同網(wǎng)格數(shù)量下懸臂梁和Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)、目標(biāo)柔度值C、整體迭代次數(shù)I以及整體優(yōu)化時間T的分布分別如圖12、圖13、表3、圖14、圖15所示.

從圖12、圖13 和表3 可以看出，在PDE 靈敏度濾波器中使用精度為1E–4 的CG 算法、MG 算法和MGCG 算法得到不同網(wǎng)格數(shù)量的懸臂梁和Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)基本相同，并且隨著網(wǎng)格數(shù)量增加，柔度值越來越大.可見PDE 濾波器求解算法對不同網(wǎng)格數(shù)量的優(yōu)化結(jié)構(gòu)影響不大，且可有效地消除網(wǎng)格依賴性和棋盤格現(xiàn)象.此外，在PDE 濾波器中使用MGCG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的迭代次數(shù)整體低于其他兩種算法.

表3 兩種結(jié)構(gòu)的拓?fù)鋬?yōu)化迭代次數(shù)和柔度值(不同網(wǎng)格數(shù)量)Table 3 Total numbers of iterations and compliance for topology optimization of 2 structures (different grids)

圖12 不同網(wǎng)格數(shù)量下CG、MG 和MGCG 算法的懸臂梁的優(yōu)化結(jié)構(gòu)：(a)網(wǎng)格為240 × 80 (CG)；(b)網(wǎng)格為240 × 80 (MG)；(c)網(wǎng)格為240 × 80 (MGCG)；(d)網(wǎng)格為480 × 160 (CG)；(e)網(wǎng)格為480 × 160 (MG)；(f)網(wǎng)格為480 × 160 (MGCG)；(g)網(wǎng)格為960 × 320 (CG)；(h)網(wǎng)格為960 × 320 (MG)；(i)網(wǎng)格為960 × 320 (MGCG)Fig.12 The optimized cantilever beams by CG,MG and MGCG algorithms with different grids:(a)grid 240 × 80 (CG);(b)grid 240 × 80 (MG);(c)grid 240 × 80 (MGCG);(d)grid 480 × 160 (CG);(e)grid 480 × 160 (MG);(f)grid 480 × 160 (MGCG);(g)grid 960 × 320 (CG);(h)grid 960 × 320 (MG);(i)grid 960 × 320 (MGCG)

圖13 不同網(wǎng)格數(shù)量下CG、MG 和MGCG 算法的Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)：(a)網(wǎng)格為240 × 120 (CG)；(b)網(wǎng)格為240 × 120 (MG)；(c)網(wǎng)格為240 × 120 (MGCG)；(d)網(wǎng)格為480 × 240 (CG)；(e)網(wǎng)格為480 × 240 (MG)；(f)網(wǎng)格為480 × 240 (MGCG)；(g)網(wǎng)格為960 × 480 (CG)；(h)網(wǎng)格為960 × 480 (MG)；(i)網(wǎng)格為960 × 480 (MGCG)Fig.13 The optimized Michell structures by CG,MG and MGCG algorithms with different grids:(a)grid 240 × 120 (CG);(b)grid 240 × 120 (MG);(c)grid 240 × 120 (MGCG);(d)grid 480 × 240 (CG);(e)grid 480 × 240 (MG);(f)grid 480 × 240 (MGCG);(g)grid 960 × 480 (CG);(h)grid 960 × 480 (MG);(i)grid 960 × 480 (MGCG)

從圖14 和圖15 可以看出，同一網(wǎng)格數(shù)量下，在PDE 濾波器中使用MGCG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的時間整體低于CG 算法和MGCG 算法.特別地，在PDE 濾波器中采用MGCG 算法對懸臂梁和Michell 結(jié)構(gòu)拓?fù)鋬?yōu)化的時間分別比CG 算法少了3% ～ 25%和64% ～ 73%，比MG 算法少了15% ～ 37%，25%和76%.因此，從結(jié)構(gòu)的整體優(yōu)化效率看，MGCG 算法高于CG 算法和MG 算法.

圖14 懸臂梁結(jié)構(gòu)優(yōu)化時間分布圖(不同網(wǎng)格數(shù)量)Fig.14 The optimized time distribution of the cantilever beam(different grids)

圖15 Michell 結(jié)構(gòu)優(yōu)化時間分布圖(不同網(wǎng)格數(shù)量)Fig.15 The optimized time distribution of the Michell structure(different grids)

4 結(jié)論

本文在PDE 靈敏度濾波器中引入了CG 算法、MG 算法和MGCG 算法，對懸臂梁和Michell 結(jié)構(gòu)進(jìn)行拓?fù)鋬?yōu)化，并比較了算法精度、過濾半徑和網(wǎng)格數(shù)量對優(yōu)化結(jié)果和優(yōu)化效率的影響.實驗結(jié)果表明：

1）不同精度的三種算法得到的優(yōu)化結(jié)構(gòu)、目標(biāo)函數(shù)值基本相同.隨著精度的提升，MG 算法的結(jié)構(gòu)優(yōu)化時間不斷增加，CG 算法總體趨于穩(wěn)定，MGCG 算法有微小的波動，且MGCG 算法結(jié)構(gòu)優(yōu)化時間整體少于其他兩種算法.

2）隨著過濾半徑的增大，CG 算法和MG 算法的結(jié)構(gòu)優(yōu)化時間大幅減少，整體呈下降趨勢，MGCG 算法有微小的波動，但優(yōu)化時間整體少于其他兩種算法，且過濾半徑越小，優(yōu)勢越明顯.

3）不同網(wǎng)格數(shù)量的三種算法得到的優(yōu)化結(jié)構(gòu)、目標(biāo)函數(shù)值基本相同.隨著網(wǎng)格數(shù)量的增加，三種算法的結(jié)構(gòu)優(yōu)化時間和迭代次數(shù)整體呈上升趨勢，且隨著網(wǎng)格數(shù)量的增加，MGCG 算法的優(yōu)勢越明顯.

因此，從結(jié)構(gòu)的整體優(yōu)化效率看，MGCG 算法優(yōu)于CG 算法和MG 算法.