孟換利，張岐良，王 杰

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

引言

為應(yīng)對(duì)越來(lái)越高的力學(xué)性能要求和越來(lái)越嚴(yán)酷的服役環(huán)境，在結(jié)構(gòu)設(shè)計(jì)中應(yīng)用拓?fù)鋬?yōu)化技術(shù)已成為一種必然選擇.在結(jié)構(gòu)的概念設(shè)計(jì)階段通過(guò)結(jié)構(gòu)拓?fù)鋬?yōu)化確定材料的最優(yōu)布局形式，可以幫助設(shè)計(jì)者發(fā)現(xiàn)結(jié)構(gòu)最佳的應(yīng)力傳導(dǎo)路徑，從而以最少的材料實(shí)現(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)用場(chǎng)景、研究熱點(diǎn)和挑戰(zhàn).衛(wèi)志軍等[3]對(duì)雙隔板進(jìn)行了拓?fù)鋬?yōu)化，可以有效地抑制流體晃蕩.

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

多重網(wǎng)格是一種高效的求解算法，近期在偏微分方程參數(shù)反演問(wèn)題方面有許多應(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)化靈敏度過(guò)濾的效率，本文基于體積約束下柔度值最小的固體各向同性懲罰微結(jié)構(gòu)(SIMP)插值模型，針對(duì)大規(guī)模PDE 靈敏度濾波器的求解問(wèn)題，對(duì)其進(jìn)行有限元分析得到代數(shù)方程，再采用共軛梯度算法、多重網(wǎng)格算法以及多重網(wǎng)格預(yù)處理共軛梯度算法進(jìn)行求解，得到了過(guò)濾域的節(jié)點(diǎn)靈敏度值，通過(guò)節(jié)點(diǎn)到單元的映射得到了單元靈敏度過(guò)濾值.

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個(gè)單元的相對(duì)密度和節(jié)點(diǎn)位移，ki和k0是單元?jiǎng)偠染仃?，V是優(yōu)化后的結(jié)構(gòu)體積，V?為結(jié)構(gòu)體積約束，V0為整個(gè)設(shè)計(jì)域的初始體積，f為優(yōu)化體積比，vi為第i個(gè)單元的體積，xmin和xmax為單元相對(duì)密度的最小和最大極限值，N為結(jié)構(gòu)離散單元總數(shù).

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

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

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

2.1 總體框架

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

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

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

2.2 PDE 靈敏度濾波器

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

2.3 PDE 靈敏度濾波器的求解

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

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

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

第2 步 通過(guò)對(duì)過(guò)濾域的節(jié)點(diǎn)靈敏度值映射得到單元靈敏度過(guò)濾值：

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

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

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

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

步驟2 粗網(wǎng)格校正.計(jì)算在粗網(wǎng)格上誤差的近似解，具體步驟為：首先計(jì)算細(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 插值.計(jì)算細(xì)網(wǎng)格上的解，具體步驟為：先使用插值算子把粗網(wǎng)格上的近似誤差插值到細(xì)網(wǎng)格上然后使用公式得到在細(xì)網(wǎng)格的解.本文采用的插值算子為

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

算法1

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

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

算法2

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

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

2.5 收斂準(zhǔn)則

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

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

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

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

3 數(shù)值算例

為了驗(yàn)證本文拓?fù)鋬?yōu)化方法的可行性、可靠性和高效性，針對(duì)長(zhǎng)960 mm、高320 mm 的懸臂梁和長(zhǎng)960 mm、高480 mm Michell 結(jié)構(gòu)的棋盤(pán)格現(xiàn)象、數(shù)值不穩(wěn)定、靈敏度過(guò)濾慢等問(wèn)題開(kāi)展了數(shù)值實(shí)驗(yàn)，如圖2 和3所示（算例優(yōu)化參數(shù)為體積約束比為0.5、彈性模量為1 GPa、Poisson 比為0.3、懲罰因子為3、移動(dòng)極限為0.2、外載荷F為1 kN）.拓?fù)鋬?yōu)化程序基于文獻(xiàn)[16,20]的MATLAB 代碼進(jìn)行編寫(xiě)，在PDE 濾波器的求解中引入了共軛梯度（CG）算法、多重網(wǎng)格（MG）算法以及多重網(wǎng)格預(yù)處理共軛梯度（MGCG）算法，并比較了三種算法對(duì)結(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ù)方程，并對(duì)整個(gè)結(jié)構(gòu)分別劃分為480×160和480×240個(gè)單元的懸臂梁和Michell 結(jié)構(gòu)進(jìn)行半徑為6 的PDE 靈敏度過(guò)濾拓?fù)鋬?yōu)化實(shí)驗(yàn).其中，不同精度下懸臂梁和Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)、目標(biāo)柔度值C、整體迭代次數(shù)I以及整體優(yōu)化時(shí)間T的分布如圖4、圖5、表1、圖6、圖7所示.

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

圖7 Michell 結(jié)構(gòu)優(yōu)化時(shí)間分布圖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 可以看出，在過(guò)濾半徑為6 的PDE 靈敏度濾波器中使用不同精度的CG 算法、MG 算法和MGCG 算法得到的懸臂梁和Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)圖基本相同，并且目標(biāo)函數(shù)值非常接近，至多相差1%，說(shuō)明PDE 濾波器求解算法的精度對(duì)拓?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)化的時(shí)間整體低于CG 算法和MGCG 算法，且隨著精度的增加，在PDE 濾波器中使用CG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的時(shí)間逐漸降低并趨于穩(wěn)定，在PDE 濾波器中使用MG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的時(shí)間逐漸增加并趨于穩(wěn)定.此外，在PDE 濾波器中采用MGCG 算法對(duì)懸臂梁和Michell 結(jié)構(gòu)拓?fù)鋬?yōu)化的時(shí)間分別比CG 算法少了11% ～ 34%和65% ～ 80%，比MG 算法少了29% ～ 46%，3%和77%.因此，從結(jié)構(gòu)的整體優(yōu)化效率看，MGCG 算法高于CG 算法和MG 算法.

3.2 過(guò)濾半徑的影響

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

圖8 不同過(guò)濾半徑下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 不同過(guò)濾半徑下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)化時(shí)間分布圖(不同過(guò)濾半徑)Fig.10 The optimized time distribution of the cantilever beam(different filter radius)

圖11 Michell 結(jié)構(gòu)優(yōu)化時(shí)間分布圖(不同過(guò)濾半徑)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ù)值的特征隨著過(guò)濾半徑的變化呈現(xiàn)類(lèi)似的變化規(guī)律.隨著過(guò)濾半徑的增大，拓?fù)鋬?yōu)化結(jié)構(gòu)的清晰度越來(lái)越低，細(xì)支結(jié)構(gòu)越來(lái)越少，說(shuō)明該兩種算法對(duì)結(jié)構(gòu)承力路徑的識(shí)別精度與過(guò)濾半徑間體現(xiàn)出負(fù)相關(guān)性，且柔度值越來(lái)越大.且從表2 可以看出，隨著過(guò)濾半徑的改變，在PDE 濾波器中使用MG 算法和MGCG 算法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化的迭代次數(shù)整體低于CG 算法.

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

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

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

分別對(duì)網(wǎng)格劃分為960×320，過(guò)濾半徑為12，層數(shù)為6；網(wǎng)格劃分為480×160，過(guò)濾半徑為6，層數(shù)為5；網(wǎng)格劃分為240×80，過(guò)濾半徑為3，層數(shù)為4 的懸臂梁和網(wǎng)格劃分為960×480，過(guò)濾半徑為12，層數(shù)為6；網(wǎng)格劃分為480×240，過(guò)濾半徑為6，層數(shù)為5；網(wǎng)格劃分為240×120，過(guò)濾半徑為3，層數(shù)為4 的Michell 結(jié)構(gòu)進(jìn)行控制精度為1E–4 的PDE 靈敏度過(guò)濾的拓?fù)鋬?yōu)化實(shí)驗(yàn).其中，不同網(wǎng)格數(shù)量下懸臂梁和Michell 結(jié)構(gòu)的優(yōu)化結(jié)構(gòu)、目標(biāo)柔度值C、整體迭代次數(shù)I以及整體優(yōu)化時(shí)間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ù)量增加，柔度值越來(lái)越大.可見(jiàn)PDE 濾波器求解算法對(duì)不同網(wǎng)格數(shù)量的優(yōu)化結(jié)構(gòu)影響不大，且可有效地消除網(wǎng)格依賴(lài)性和棋盤(pán)格現(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)化的時(shí)間整體低于CG 算法和MGCG 算法.特別地，在PDE 濾波器中采用MGCG 算法對(duì)懸臂梁和Michell 結(jié)構(gòu)拓?fù)鋬?yōu)化的時(shí)間分別比CG 算法少了3% ～ 25%和64% ～ 73%，比MG 算法少了15% ～ 37%，25%和76%.因此，從結(jié)構(gòu)的整體優(yōu)化效率看，MGCG 算法高于CG 算法和MG 算法.

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

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

4 結(jié)論

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

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

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

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

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