基于免組裝有限元的懸架控制臂拓?fù)鋬?yōu)化

卞翔+方宗德

摘 要：為了提高拓?fù)鋬?yōu)化在實際工程應(yīng)用中的速度，將一種基于拓?fù)潇`敏度的拓?fù)鋬?yōu)化算法與免組裝有限元法結(jié)合，對汽車懸架控制臂進(jìn)行多工況拓?fù)鋬?yōu)化.該方法使用一致的體素單元劃分網(wǎng)格，用壓縮共軛梯度法（DCG）加速有限元的求解，得到結(jié)構(gòu)的應(yīng)力和應(yīng)變場，從而獲得拓?fù)潇`敏度場，用來控制結(jié)構(gòu)拓?fù)渥兓?通過添加多工況優(yōu)化以及加工約束功能，使該方法更符合工程應(yīng)用的要求.通過對比，該方法和商用軟件能夠獲得性能接近、形狀相似的優(yōu)化結(jié)構(gòu)，而整個優(yōu)化過程的速度有顯著提高.

關(guān)鍵詞：懸架控制臂；多工況拓?fù)鋬?yōu)化；免組裝有限元；體素化；壓縮共軛梯度法；拓?fù)潇`敏度

中圖分類號：U463.33 文獻(xiàn)標(biāo)志碼：A

Topology Optimization of Suspension Control Arm Based on Assembly-free Finite Element Method BIAN Xiang， FANG Zongde

（Department of Mechanical Engineering， Northwestern Polytechnic University， Xian 710072， China）

Abstract： A suspension control arm was optimized by topology optimization under multi-working condition. In order to improve the speed of topology optimization for practical engineering problems， a topology optimization method based on topological sensitivity was combined with an assembly-free finite element method. The structure was discretized via uniform voxels. A deflated conjugate gradient （DCG） method was used to accelerate the finite element analysis. After the finite element analysis， the stress and strain field can be used to generate the topological sensitivity field， and the topological sensitivity field can be used to control the change of topology. Multi-working condition optimization and manufacturing constraints were used to make the method more suitable for practical problem. The comparison shows that the method can result in the close performance and similar shape of optimized structure to the commercial software. Meanwhile， the speed of the entire optimization is increased significantly.

Key words：suspension control arm； topology optimization under multi-working condition； assembly-free finite element method； voxelization； deflated conjugate gradient； topological sensitivity

輕量化設(shè)計是汽車設(shè)計中的一項重要內(nèi)容，對于降低重量、提高性能、降低油耗等有重要作用.拓?fù)鋬?yōu)化技術(shù)在汽車結(jié)構(gòu)的輕量化設(shè)計中運(yùn)用越來越多，拓?fù)鋬?yōu)化能在結(jié)構(gòu)設(shè)計的初始階段提供一個概念性設(shè)計，使結(jié)構(gòu)在布局上采用最優(yōu)方案，改變了以往的設(shè)計、校核、修改這樣一個不斷反復(fù)的開發(fā)流程，可以提高設(shè)計效率和質(zhì)量，縮短研發(fā)周期，降低成本.

近年來，拓?fù)鋬?yōu)化技術(shù)逐漸成熟，國內(nèi)外研究較多的有固體各向同性材料懲罰法（Solid Isotropic Material with Penalization，SIMP）[1]，進(jìn)化結(jié)構(gòu)優(yōu)化法（Evolutionary Structural Optimization，ESO）[2-3]，水平集法（Level-Set）[4-5]等.其中SIMP方法基于有限元法，給每個單元賦予偽密度，通過優(yōu)化各單元的偽密度來達(dá)到優(yōu)化目標(biāo)[6].由于其概念簡單，算法容易實現(xiàn)，被推廣應(yīng)用于多材料、多工況、多物理場等多種實際問題中，成為主流的拓?fù)鋬?yōu)化方法之一[7].目前，常用的商用拓?fù)鋬?yōu)化軟件如 Hyperworks 和 ANSYS 都采用SIMP方法進(jìn)行拓?fù)鋬?yōu)化.國內(nèi)有許多研究使用Hyperworks進(jìn)行汽車零部件的拓?fù)鋬?yōu)化設(shè)計，文獻(xiàn)[8-11]使用Hyperworks對轉(zhuǎn)向節(jié)、懸架控制臂、發(fā)動機(jī)懸置支架、車架等汽車結(jié)構(gòu)進(jìn)行了拓?fù)鋬?yōu)化設(shè)計.

如今，拓?fù)鋬?yōu)化技術(shù)發(fā)展的挑戰(zhàn)之一是計算速度[12]，尤其是對于大型有限元模型，比如對于上百萬自由度的優(yōu)化問題需要數(shù)小時甚至數(shù)天來完成.文獻(xiàn)[12-14]提出了一種基于拓?fù)潇`敏度的拓?fù)鋬?yōu)化算法，對多目標(biāo)優(yōu)化問題帕雷托最優(yōu)解的計算十分高效，該方法將拓?fù)潇`敏度作為水平集，用來控制孔洞的形成，得到的結(jié)果相比于SIMP方法更加清晰、明確，不需要網(wǎng)格過濾等技術(shù).但是這種方法的研究還處于初步階段，沒有商用軟件全面的功能，還不能滿足實際的工程應(yīng)用.

本文基于C++語言，將這種基于拓?fù)潇`敏度的拓?fù)鋬?yōu)化算法與一種免組裝（Assembly-free）有限元分析方法[15-16]相結(jié)合，并用壓縮共軛梯度法[17]進(jìn)一步加速拓?fù)鋬?yōu)化的速度，通過添加多工況優(yōu)化、加工約束等滿足實際工程應(yīng)用的功能，用于汽車零部件的結(jié)構(gòu)優(yōu)化設(shè)計.以懸架控制臂的優(yōu)化為例，通過與商用軟件Hyperworks對比，證明該方法的有效性.

1 免組裝有限元模型

1.1 懸架控制臂體素化模型

一些商用軟件需要花費(fèi)大量時間在網(wǎng)格劃分上，需要人工進(jìn)行簡化幾何、拆解幾何和改變網(wǎng)格尺寸等操作，才能獲得滿足軟件要求的網(wǎng)格模型.

體素化[18]是一種簡單的有限元網(wǎng)格離散方法，通過統(tǒng)一的六面體單元（體素）來劃分幾何模型.由于使用了形狀和尺寸一樣的單元，體素化網(wǎng)格具有穩(wěn)健性和內(nèi)存占用空間低的優(yōu)點[19]，能夠配合免組裝方法，提高有限元分析的計算速度.

體素化模型由于使用一致的網(wǎng)格，可以通過編程實現(xiàn)快速的自動網(wǎng)格劃分，圖 1 （a）所示是懸架控制臂的幾何模型，圖1（b）為相應(yīng)的體素網(wǎng)格模型，單元數(shù)量為79 000，自由度為280 000，自動網(wǎng)格劃分時間僅為4 s.

由于網(wǎng)格一致性的限制，體素化模型不能夠完全精確體現(xiàn)所有的幾何細(xì)節(jié).但是考慮到拓?fù)鋬?yōu)化是設(shè)計初期的概念設(shè)計階段，其主要目的是為后期設(shè)計提供參考，在多數(shù)情況下，尤其是對于諸如柔度、頻率等全局參數(shù)的優(yōu)化變量，一些幾何細(xì)節(jié)并不重要.對于本文考慮的懸架控制臂模型，3個連接端這幾個細(xì)節(jié)部分都屬非設(shè)計區(qū)域，細(xì)節(jié)的計算精度并不影響主體設(shè)計區(qū)域的拓?fù)鋬?yōu)化，所以體素化網(wǎng)格適用于本文對懸架控制臂的拓?fù)鋬?yōu)化.

1.2 免組裝的壓縮共軛梯度法

拓?fù)鋬?yōu)化中的迭代過程需要進(jìn)行多次有限元求解，本文將基于壓縮共軛梯度法的免組裝有限元算法引入到拓?fù)鋬?yōu)化過程中，可以提高每次有限元分析的速度，從而大幅縮短整個拓?fù)鋬?yōu)化過程的時間.

免組裝（Assembly-free）有限元法，即總剛矩陣免組裝的有限元法，最初是Hughes等[15]在1983年提出的，隨著并行算法的發(fā)展，這種方法也得到了改進(jìn)[17].免組裝有限元方法的基本概念是單元剛度矩陣不用組裝成總體剛度矩陣，而是在單元層面進(jìn)行稀疏矩陣向量乘法.換而言之是將“組裝然后相乘”的過程：

變?yōu)椤跋喑巳缓蠼M裝”：

式中：K為總體剛度矩陣；d為位移向量.

共軛梯度法Conjugate Gradient （CG）是求解大型稀疏線性方程的一種迭代算法，迭代過程中最占用計算時間的是矩陣向量乘法運(yùn)算Kd，其大部分的內(nèi)存占用來自存儲和提取剛度矩陣.如上節(jié)所述，由于結(jié)構(gòu)使用統(tǒng)一的單元進(jìn)行網(wǎng)格劃化，所有單元的形狀和尺寸是一樣的，所以單元剛度矩陣是一樣的，不需要組裝并存儲總體剛度矩陣K，而只需要儲存單個單元剛度矩陣Ke.這將大大減少內(nèi)存占用，提高矩陣向量乘法運(yùn)算的速度，從而提高迭代算法的計算速度.

壓縮共軛梯度法Deflated Conjugate Gradient （DCG）是一種加速的迭代算法[20].它將有限元網(wǎng)格的點劃分成少數(shù)的幾塊，把每塊網(wǎng)格當(dāng)作剛體處理，構(gòu)造出壓縮空間.在某塊里的一個點的位移表示為

式中：（u0，v0，w0， θx， θy， θz）T是此塊在6個自由度上的剛體運(yùn)動.（x，y，z）是該塊中一點相對此塊幾何中心的坐標(biāo).通過所有點的相對坐標(biāo)，構(gòu)造出壓縮矩陣W，于是有：

式中：d有3N個自由度（N是所有點個數(shù)）；λ是6G個與塊相關(guān)的自由度（G是分塊的個數(shù)）.可以使用壓縮空間矩陣W進(jìn)行共軛梯度運(yùn)算.（具體推導(dǎo)及理論分析可參考文獻(xiàn)[21]）.由于分塊的數(shù)量G遠(yuǎn)小于節(jié)點數(shù)量N，所以壓縮矩陣的尺度比剛度矩陣的小，在迭代過程中能夠提高矩陣向量乘法運(yùn)算的速度，從而對迭代過程起到加速作用.同時，這種壓縮共軛梯度算法也能使用Assembly-free方法，就像在式（2）中的Kd運(yùn)算不需要組裝總體剛度矩陣，總體壓縮矩陣W也不需要組裝，而是采用“相乘然后組裝”的方式：

這種免組裝的壓縮共軛梯度方法能夠高效地處理高達(dá)幾百萬自由度的大型有限元問題[17].

如圖 2所示，本文將懸架控制臂的網(wǎng)格模型分為200個組.使用免組裝的壓縮共軛梯度法求解靜力學(xué)問題，并得到位移和應(yīng)力場.

2 多工況懸架控制臂拓?fù)鋬?yōu)化

2.1 拓?fù)潇`敏度控制的水平集優(yōu)化算法

不同于利用偽密度的SIMP方法，本文使用的拓?fù)鋬?yōu)化方法基于拓?fù)潇`敏度.拓?fù)潇`敏度是當(dāng)拓?fù)渖习l(fā)生極小的改變時，目標(biāo)量的變化率.這是由Eschenauer[22]最早研究的，隨后諸多學(xué)者[23-25]對其進(jìn)行了擴(kuò)展研究.

這里通過圖 3中的二維例子對拓?fù)潇`敏度進(jìn)行闡述，研究的目標(biāo)量是Q，假設(shè)結(jié)構(gòu)域某處p去除一個極小的半徑為r的小孔，結(jié)構(gòu)的變形將會發(fā)生改變，所關(guān)心的目標(biāo)量Q也會改變，拓?fù)潇`敏度TQ定義為[22]：

本文的目標(biāo)量是柔度J，從公式（6）可以推導(dǎo)出關(guān)于柔度拓?fù)潇`敏度[26]的解析表達(dá)式：

式中：TJ為關(guān)于柔度J的拓?fù)潇`敏度；p表示網(wǎng)格中的某點；σ為應(yīng)力張量；ε為應(yīng)變張量；v為泊松比.

設(shè)計域內(nèi)所有點的拓?fù)潇`敏度構(gòu)成一個靈敏度場，由定義可知，這個靈敏度場中具有相對較高值的區(qū)域表示這個部分對于目標(biāo)量Q相對重要.使用此拓?fù)潇`敏度場作為水平集，可以用來引入孔洞，確定新的結(jié)構(gòu)域.例如圖 4中，引入了一個任意閾值為τ=0.02的“切割”平面，通過式（8）可以確定圖 5所示的拓?fù)浣Y(jié)構(gòu)Ωτ：

式中的閾值τ是通過當(dāng)前所需的體積比確定的.結(jié)構(gòu)域Ωτ是所有拓?fù)潇`敏度超過τ的點的集合.

然而直接得到的結(jié)構(gòu)域Ωτ并不一定是pareto最優(yōu)解[13]，需要反復(fù)進(jìn)行以下3步：1）對生成的結(jié)構(gòu)域進(jìn)行有限元分析；2）重新計算拓?fù)潇`敏度；3）根據(jù)當(dāng)前目標(biāo)體積分?jǐn)?shù)確定新的閾值τ，并生成新結(jié)構(gòu)域.通常這個過程需要3到4次迭代能夠收斂[12].收斂之后，可以繼續(xù)降低當(dāng)前目標(biāo)體積分?jǐn)?shù)，重復(fù)上述過程，直至得到最終設(shè)計要求的體積分?jǐn)?shù)f V0.總體的流程如圖 6所示，優(yōu)化算法的步驟概括為：

①從初始結(jié)構(gòu)Ω=Ω0開始，初始體積V=V0；

②對初始結(jié)構(gòu)Ω0進(jìn)行有限元分析，通過公式（7）計算拓?fù)潇`敏度；

③按給定的體積減少步長ΔV，確定當(dāng)前目標(biāo)體積V=V-ΔV；

④按當(dāng)前目標(biāo)體積和拓?fù)潇`敏度場，確定水平集參數(shù)τ，使新結(jié)構(gòu)Ωτ的體積等于當(dāng)前目標(biāo)體積V；

⑤進(jìn)行迭代處理：對新結(jié)構(gòu)Ωτ進(jìn)行有限元分析、計算拓?fù)潇`敏度、得到新結(jié)構(gòu)……直到目標(biāo)函數(shù)J收斂；

⑥判斷當(dāng)前的體積V是否達(dá)到最終設(shè)計要求的體積fV0，如果沒有達(dá)到，繼續(xù)按照步長ΔV降低目標(biāo)體積分?jǐn)?shù)，返回步驟②；

⑦當(dāng)前體積達(dá)到設(shè)計要求的體積分?jǐn)?shù)時，算法終止.

2.2 多工況優(yōu)化模型

多工況條件下懸架控制臂的拓?fù)鋬?yōu)化模型為

式中：J為綜合目標(biāo)函數(shù)；m為工況總數(shù)；Ji為第i個工況的柔度；ωi為第i個工況的權(quán)重；V為優(yōu)化后的體積；V0為結(jié)構(gòu)初始體積；f為體積約束的百分比.

根據(jù)各工況的權(quán)重系數(shù)，可以構(gòu)造一個新的多工況條件下的拓?fù)潇`敏度：

式中：Ti為第i個工況下的拓?fù)潇`敏度.

各工況的權(quán)重系數(shù)可以通過經(jīng)驗法、層次分析法等確定.由于如何選擇權(quán)重系數(shù)并不是本文研究重點，本文重在使用相同的工況和權(quán)重系數(shù)設(shè)置與基于SIMP方法的商用軟件作對比.所以本文考慮表 1所示的穩(wěn)態(tài)轉(zhuǎn)向和直線制動兩種最基本的工況[27]，假設(shè)兩種工況的權(quán)重系數(shù)都為0.5.

對于本文考慮的麥弗遜懸架的下控制臂，在前襯套A和后襯套B通過鉸接副與車架相連，在外球銷點C通過球形副與轉(zhuǎn)向節(jié)相連.主要在加速、制動時承受縱向力Fx，以及在轉(zhuǎn)向時承受側(cè)向力Fy.對于垂向力，下控制臂只是抵消前后橡膠襯套扭轉(zhuǎn)變形時的一些結(jié)構(gòu)反力，而垂向力主要由懸架彈簧來承受，控制臂承受的垂向力的數(shù)量級遠(yuǎn)小于縱向力及側(cè)向力，所以在分析下控制臂時通常不考慮垂向力Fz.分析麥弗遜懸架控制臂時通常固定前襯套X，Y，Z 3個方向平動自由度，后襯套Y，Z方向平動自由度，外球銷點Z方向平動自由度.制動、轉(zhuǎn)向時的縱向、側(cè)向力分為兩個工況施加到外球銷點上.

2.3 加工約束

通過拓?fù)鋬?yōu)化方法進(jìn)行性能最優(yōu)化設(shè)計的結(jié)果常常會出現(xiàn)中空的復(fù)雜結(jié)構(gòu)，不適合傳統(tǒng)加工方法.目前，大多數(shù)汽車零部件還是采用傳統(tǒng)的方法加工，所以此懸架控制臂的拓?fù)鋬?yōu)化中需要考慮加工制造約束.為了使控制臂結(jié)構(gòu)易于制造加工，在Z方向添加方向約束.具體方法就是對Z方向的各單元的拓?fù)潇`敏度進(jìn)行額外處理，使外側(cè)單元的拓?fù)潇`敏度不大于內(nèi)側(cè)單元的拓?fù)潇`敏度，這樣就能保證在拔模方向上不會出現(xiàn)中空的結(jié)構(gòu).值得注意的是，由于使用了體素單元劃分網(wǎng)格，網(wǎng)格的排列十分工整，有利于編程實現(xiàn)加工制造約束.

3 優(yōu)化結(jié)果

使用本文方法基于C++語言開發(fā)的程序和Hyperworks分別進(jìn)行拓?fù)鋬?yōu)化，設(shè)置相同的材料參數(shù)，彈性模量210 GPa，泊松比0.3，使用前文所述的工況設(shè)置，網(wǎng)格模型自由度為280 000，目標(biāo)體積分?jǐn)?shù)為50%.優(yōu)化后的結(jié)果分別如圖 7 （a） （b）所示.圖 8給出了優(yōu)化過程中兩種工況下的柔度變化曲線，以及優(yōu)化過程中的一些拓?fù)浣Y(jié)構(gòu)和應(yīng)力云圖.

兩種方法的最優(yōu)拓?fù)渎杂胁煌艘恍┘?xì)小結(jié)構(gòu)的區(qū)別，主要的材料分布趨勢是近似的.根據(jù)最優(yōu)拓?fù)涮峁┑膮⒖迹瑢Y(jié)果進(jìn)行幾何重構(gòu)后如圖9所示，性能對比如表2所示.

通過表 2中各工況的柔度、最大位移、最大應(yīng)力的對比可知，兩種結(jié)構(gòu)的性能接近，且本文方法優(yōu)化的結(jié)構(gòu)性能略優(yōu)于Hyperworks優(yōu)化的結(jié)果.而在計算時間上，兩種方法有很大不同，本文方法是Hyperworks的26 %，可見對于大自由度的模型，本文方法有明顯的速度優(yōu)勢.這主要是由于DCG算法對有限元求解的加速作用.圖10表明了使用不同分組數(shù)的DCG算法的收斂性，分組數(shù)量的不同導(dǎo)致DCG運(yùn)算收斂所需的迭代次數(shù)不同.如圖所示，分成200組比未分組的CG運(yùn)算所需的迭代次數(shù)大幅降低，這意味著單次有限元分析所需的時間大幅降低，從而縮短了整個拓?fù)鋬?yōu)化過程的時間.不同分組數(shù)下相應(yīng)的拓?fù)鋬?yōu)化時間如表3所示.

迭代次數(shù)

此外，對于復(fù)雜結(jié)構(gòu)，商用軟件需要人工進(jìn)行幾何分解、網(wǎng)格劃分、網(wǎng)格質(zhì)量檢查等步驟，直到滿足商用軟件的特定要求后才能開始拓?fù)鋬?yōu)化.而本文方法采用體素化網(wǎng)格實現(xiàn)自動網(wǎng)格劃分，在拓?fù)鋬?yōu)化的前處理階段就可以節(jié)省大量的時間.

4 結(jié) 論

本文將一種基于拓?fù)潇`敏度的拓?fù)鋬?yōu)化算法與免組裝壓縮共軛梯度法結(jié)合，提高了大型三維拓?fù)鋬?yōu)化的速度，通過添加多工況優(yōu)化和加工約束功能，使其更適合實際工程應(yīng)用.通過懸架控制臂的多工況優(yōu)化，與基于SIMP方法的商用軟件對比，結(jié)果表明，該算法具有明顯的速度優(yōu)勢，且能獲得性能相近的結(jié)構(gòu).由于拓?fù)鋬?yōu)化是初始的概念設(shè)計，體素化網(wǎng)格的精度局限性對拓?fù)鋬?yōu)化結(jié)果并不會產(chǎn)生很大影響，而速度優(yōu)勢可以在前處理階段以及拓?fù)鋬?yōu)化階段得到很好的發(fā)揮.在未來的研究中，可以拓展更全面的優(yōu)化功能，使該算法能更好地滿足實際工程應(yīng)用的要求.

參考文獻(xiàn)

[1] SIGMUND O. A 99 line topology optimization code written in Matlab[J]. Structural & Multidisciplinary Optimization， 2001， 21（2）：120-127.

[2] HUANG X， XIE Y M. A new look at ESO and BESO optimization methods[J]. Structural & Multidisciplinary Optimization， 2008， 35（1）：89-92.

[3] MUNK D J， VIO G A， STEVEN G P. Topology and shape optimization methods using evolutionary algorithms： a review[J]. Structural & Multidisciplinary Optimization， 2015， 52（3）：613-631.

[4] MEI Y， WANG X. A level set method for structural topology optimization and its applications[J]. Advances in Engineering Software， 2004， 35（7）：415-441.

[5] DIJK N P V， MAUTE K， LANGELAAR M， et al. Level-set methods for structural topology optimization： a review[J]. Structural & Multidisciplinary Optimization， 2013， 48（3）：437-472.

[6] BENDSE M P， SIGMUND O. Material interpolation schemes in topology optimization[J]. Archive of Applied Mechanics， 1999， 69（9/10）： 635-654.

[7] ROZVANY G I N. A critical review of established methods of structural topology optimization[J]. Structural & Multidisciplinary Optimization， 2009， 37（3）：217-237.

[8] 祝小元，方宗德，申閃閃，等. 汽車懸架控制臂的多目標(biāo)拓?fù)鋬?yōu)化[J]. 汽車工程， 2011， 33（2）： 138-141.

ZHU Xiaoyuan， FANG Zongde， SHEN Shanshan， et al. Multi-objective topology optimization for the control arm of vehicle suspension[J]. Automotive Engineering， 2011， 33（2）： 138-141. （In Chinese）

[9] 蘭鳳崇，張浩鍇，王家豪，等. 汽車轉(zhuǎn)向節(jié)拓?fù)鋬?yōu)化方法研究及應(yīng)用[J]. 汽車工程， 2014， 36（4）： 464-468.

LAN Fengchong， ZHANG Haokai，WANG Jiahao， et al. Study and application of topology optimization technique for vehicle steering knuckles[J]. Automotive Engineering， 2014， 36（4）： 464-468. （In Chinese）

[10]朱劍峰，林逸，施國標(biāo)，等. 考慮工程約束的發(fā)動機(jī)懸置支架拓?fù)鋬?yōu)化[J]. 汽車工程， 2014， 36（12）： 1508-1512.

ZHU Jianfeng， LIN Yi， SHI Guobiao， et al. Topology optimization of engine mount bracket with consideration of engineering constraints[J]. Automotive Engineering， 2014， 36（12）： 1508-1512. （In Chinese）

[11]張偉， 侯文彬， 胡平. 基于拓?fù)鋬?yōu)化的電動汽車白車身優(yōu)化設(shè)計[J]. 湖南大學(xué)學(xué)報：自然科學(xué)版， 2014，41（10）：42-48.

ZHANG Wei， HOU Wenbin， HU Ping. The body in white optimization of an electric vehicle using topology optimization[J]. Journal of Hunan University：Natural Sciences， 2014，41（10）：42-48. （In Chinese）

[12]SURESH K. Efficient generation of large-scale pareto-optimal topologies[J]. Structural & Multidisciplinary Optimization， 2013， 47（1）：49-61.

[13]SURESH K. A 199-line matlab code for pareto-optimal tracing in topology optimization[J]. Structural & Multidisciplinary Optimization， 2010， 42（5）：665-679.

[14]INNA T， KRISHNAN S. Efficient generation of pareto-optimal topologies for compliance optimization[J]. International Journal for Numerical Methods in Engineering， 2011， 87（12）：1207-1228.

[15]HUGHES T J R， LEVIT I， WINGET J. An element-by-element solution algorithm for problems of structural and solid mechanics[J]. Computer Methods in Applied Mechanics & Engineering， 1983， 36（2）：241-254.

[16]YADAV P， SURESH K. Assembly-free large-scale modal analysis on the graphics-programmable unit[J]. Journal of Computing & Information Science in Engineering， 2013， 13（1）：1-11.

[17]YADAV P， SURESH K. Large scale finite element analysis via assembly-free deflated conjugate gradient[J]. Journal of Computing & Information Science in Engineering， 2014， 14（4）： 1-9.

[18]KARABASSI E A， PAPAIOANNOU G， THEOHARIS T. A fast depth-buffer-based voxelization algorithm[J]. Journal of Graphics Tools， 1999， 4（4）： 5-10.

[19]DSTER A， PARVIZIAN J， YANG Z， et al. The finite cell method for three-dimensional problems of solid mechanics[J]. Computer Methods in Applied Mechanics & Engineering， 2008， 197（45/48）：3768-3782.

[20]AUBRY R， MUT F， DEY S. Deflated preconditioned conjugate gradient solvers for linear elasticity[J]. International Journal for Numerical Methods in Engineering， 2011， 88（11）： 1112-1127.

[21]SAAD Y， YEUNG M， ERHEL J， et al. A deflated version of the conjugate gradient algorithm[J]. Siam Journal on Scientific Computing， 1998， 21（5）： 1909-1926.

[22]ESCHENAUER H A， KOBELEV V V， SCHUMACHER A. Bubble method for topology and shape optimization of structures[J]. Structural & Multidisciplinary Optimization， 1994， 8（8）：42-51.

[23]PARS J， NAVARRINA F， COLOMINAS I， et al. Topology optimization of continuum structures with local and global stress constraints[J]. Structural & Multidisciplinary Optimization， 2010， 39（4）：419-437.

[24]STOLPE M， STIDSEN T. A hierarchical method for discrete structural topology design problems with local stress and displacement constraints[J]. International Journal for Numerical Methods in Engineering， 2007， 69（5）：1060-1084.

[25]RAMANI A. Multi-material topology optimization with strength constraints[J]. Structural & Multidisciplinary Optimization， 2011， 43（5）：597-615.

[26]GIUSTI S M， NOVOTNY A A， PADRA C. Topological sensitivity analysis of inclusion in two-dimensional linear elasticity[J]. Engineering Analysis with Boundary Elements， 2008， 32（11）：926-935.

[27]上官文斌，蔣翠翠，潘孝勇. 汽車懸架控制臂的拓?fù)鋬?yōu)化與性能計算[J]. 汽車工程， 2008， 30（8）： 709-712.

SHANGGUAN Wenbin， JIANG Cuicui， PAN Xiaoyong. Topology optimization and performance calculation for control arm of vehicle suspension[J]. Automotive Engineering， 2008， 30（8）： 709-712. （In Chinese）