二維結(jié)構(gòu)拓?fù)鋬?yōu)化的GPU并行計(jì)算方法研究

劉家彤 ，王春潔 ，2，吳 健 ，付志方

（1.北京航空航天大學(xué)機(jī)械工程及自動(dòng)化學(xué)院，北京 100191；2.北京航空航天大學(xué)虛擬現(xiàn)實(shí)技術(shù)與系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100191）

1 引言

拓?fù)鋬?yōu)化致力于在設(shè)計(jì)域內(nèi)找到最優(yōu)結(jié)構(gòu)，使得結(jié)構(gòu)在滿足給定約束條件的情況下具有最佳的工作性能[1]。自從文獻(xiàn)[2]對(duì)這一問題進(jìn)行研究開始，拓?fù)鋬?yōu)化已經(jīng)在多個(gè)領(lǐng)域得到研究應(yīng)用，包括結(jié)構(gòu)剛度最大化設(shè)計(jì)[3-4]、柔性機(jī)構(gòu)設(shè)計(jì)[5-6]等。盡管拓?fù)鋬?yōu)化在過去幾十年得到了廣泛的應(yīng)用，但由于進(jìn)行拓?fù)鋬?yōu)化時(shí)涉及結(jié)構(gòu)有限元分析、大規(guī)模線性方程組的求解和結(jié)構(gòu)單元靈敏度求解等復(fù)雜計(jì)算過程非常費(fèi)時(shí)，提高計(jì)算效率是拓?fù)鋬?yōu)化設(shè)計(jì)中急需解決的重要問題之一。

由于其優(yōu)化過程穩(wěn)定，固體各項(xiàng)同性材料懲罰法（Solid Isotropic Material with Penalization，SIMP）是結(jié)構(gòu)拓?fù)鋬?yōu)化中常用的方法之一，經(jīng)過學(xué)者的努力，已經(jīng)應(yīng)用于解決各種優(yōu)化問題[7-8]。為了方便其應(yīng)用研究，文獻(xiàn)[9]給出了SIMP方法的MATLAB實(shí)現(xiàn)代碼，俗稱“99”行代碼，然而其迭代過程需要重復(fù)計(jì)算剛度矩陣導(dǎo)致其計(jì)算效率較低。因此，文獻(xiàn)[10]對(duì)其進(jìn)行改進(jìn)，提出了“88”行代碼，從而提高了計(jì)算效率。但是，對(duì)于大型結(jié)構(gòu)優(yōu)化問題，“88”行代碼計(jì)算能力尤顯不足，仍需要開發(fā)更加高效的算法。

GPU計(jì)算是由CPU和GPU協(xié)同工作實(shí)現(xiàn)的，CPU用于控制計(jì)算過程和實(shí)施存儲(chǔ)策略，而具體的計(jì)算過程則由GPU進(jìn)行[11]。CUDA（Compute Unified Device Architecture）是NVIDIA 推出的通用并行計(jì)算架構(gòu)。利用CUDA可以直接操作GPU內(nèi)存，從而使得GPU能夠進(jìn)行除圖形處理之外的其他一般化應(yīng)用[12-13]。由于GPU計(jì)算計(jì)算過程高度并行化，使得大規(guī)模的計(jì)算能夠在較短時(shí)間內(nèi)得到結(jié)果，近年來得到了廣泛的應(yīng)用和發(fā)展[14]。

基于GPU計(jì)算，以SIMP法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化為例，利用CUDA對(duì)其進(jìn)行并行化改進(jìn)，給出了一種GPU并行計(jì)算拓?fù)鋬?yōu)化方法，并與現(xiàn)有計(jì)算方法進(jìn)行了比較。

2 基于SIMP法的拓?fù)鋬?yōu)化方法

2.1 問題描述

基于SIMP法的拓?fù)鋬?yōu)化方法以單元密度ρe作為設(shè)計(jì)變量，在優(yōu)化過程中通過改變各單元密度大小來獲得最小的結(jié)構(gòu)柔度，其數(shù)學(xué)模型如下：

式中：c（ρe）—結(jié)構(gòu)柔度；U、F—結(jié)構(gòu)位移向量和加載的力向量；K—結(jié)構(gòu)總體剛度矩陣；ue—單元位移向量；k0—單位楊氏模量下的單元?jiǎng)偠染仃嚕籒—設(shè)計(jì)域的劃分單元數(shù)；ρ—由單元密度組成的設(shè)計(jì)變量向量；V（ρ）、V0—材料體積和設(shè)計(jì)域體積；f—給定的結(jié)構(gòu)體積分?jǐn)?shù)。

為了防止單元密度變?yōu)榱悖瑢?dǎo)致剛度矩陣奇異的情況，采用約束最小剛度的插值模型[15]計(jì)算材料屬性Ee，即：

式中：E0—材料的彈性模量；Emin—防止剛度矩陣奇異而賦予無材料單元很小的剛度值，p≥1是懲罰因子，取p=3。

2.2 優(yōu)化算法

采用了優(yōu)化準(zhǔn)則法對(duì)優(yōu)化問題進(jìn)行求解，準(zhǔn)則可表述為[16]：

式中：m—一個(gè)正的偏移量；η—數(shù)值衰減系數(shù)，中間變量Be表示如下：

其中拉格朗日乘子λ可由二分查找法確定，以滿足體積約束要求。目標(biāo)函數(shù)和材料體積對(duì)于單元密度的靈敏度可表示如下[10]：

2.3 過濾策略

為了保證拓?fù)鋬?yōu)化過程的穩(wěn)定性，需要對(duì)單元密度進(jìn)行過濾，過濾方法如下[17]：

式中：Ne—給定過濾半徑范圍內(nèi)的有限單元集合；駐（e，i）—當(dāng)前單元e和目標(biāo)單元i之間的中心點(diǎn)距離；rmin—過濾半徑；He，i—由 駐（e，i）和 rmin所決定的權(quán)值。

3 算法實(shí)現(xiàn)

主要在結(jié)構(gòu)有限元分析、靈敏度計(jì)算、單元密度過濾和優(yōu)化更新階段引入GPU計(jì)算，來實(shí)現(xiàn)對(duì)SIMP拓?fù)鋬?yōu)化算法的加速。

3.1 有限元分析

針對(duì)2維平面拓?fù)鋬?yōu)化問題，設(shè)計(jì)域?yàn)榈芽栕鴺?biāo)系內(nèi)的一個(gè)矩形區(qū)域，如圖1所示。對(duì)其進(jìn)行有限元?jiǎng)澐?，得到nelx×nely個(gè)矩形單元，每個(gè)單元包含四個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)與四個(gè)矩形單元鄰接（邊界和頂點(diǎn)處，分別為1個(gè)和2個(gè)）。

圖1 有限元網(wǎng)格劃分Fig.1 Finite Element Mesh Generation

傳統(tǒng)方法在進(jìn)行有限元分析時(shí)，以單個(gè)單元為基礎(chǔ)，先構(gòu)建出一個(gè)單元的剛度矩陣Ke，通過裝配得到整體剛度矩陣。為了通過GPU并行化計(jì)算對(duì)這一過程進(jìn)行加速，有兩種加速方案：針對(duì)單元的并行化處理及針對(duì)節(jié)點(diǎn)的并行化處理。針對(duì)單元的并行處理，即以單元為劃分，針對(duì)每個(gè)單元進(jìn)行剛度矩陣的組裝，由于每個(gè)單元的節(jié)點(diǎn)可能分屬于其他三個(gè)單元，對(duì)于其中任意一個(gè)節(jié)點(diǎn)需要將幾個(gè)單元計(jì)算得到的值進(jìn)行合成，這樣就出現(xiàn)了在并行過程中的讀寫競(jìng)爭(zhēng)問題，解決這一問題會(huì)引入其他計(jì)算而增加額外的計(jì)算時(shí)間。而針對(duì)節(jié)點(diǎn)的并行化處理只需一次計(jì)算得到節(jié)點(diǎn)處的剛度值，不存在讀寫競(jìng)爭(zhēng)問題，因此在裝配剛度矩陣時(shí)選擇了該處理方式。

為了進(jìn)一步說明以上兩種處理方式的不同，其剛度矩陣裝配示意圖，如圖2所示。圖中左側(cè)為任意一個(gè)單元，以順時(shí)針方向?qū)ζ涔?jié)點(diǎn)進(jìn)行編號(hào)，通過計(jì)算得到單元的剛度矩陣Ke，右圖是以節(jié)點(diǎn)e為中心的四個(gè)單元，以左上的第一個(gè)單元為例，節(jié)點(diǎn)e位于相對(duì)左圖中的3節(jié)點(diǎn)，對(duì)應(yīng)了Ke中第5，6行的剛度值，對(duì)另外三個(gè)單元進(jìn)行同樣的分析，就針對(duì)節(jié)點(diǎn)e對(duì)剛度矩陣Ke進(jìn)行了重組，即可將原本以單元為基礎(chǔ)的剛度矩陣轉(zhuǎn)變成為以節(jié)點(diǎn)為基礎(chǔ)的剛度數(shù)組

圖2 從單元到節(jié)點(diǎn)的單元?jiǎng)偠染仃囍亟M方式Fig.2 Reconstruction of Element Stiffness Matrix from Element Based to Node Based

3.2 GPU實(shí)現(xiàn)

圖3 線程塊網(wǎng)格劃分Fig.3 Grid of Thread Blocks Partition

將文獻(xiàn)[10]中剛度矩陣裝配的計(jì)算時(shí)間和文本GPU計(jì)算方法剛度矩陣裝配時(shí)間進(jìn)行比較，測(cè)試結(jié)果，如表1所示?？梢钥闯?，采用節(jié)點(diǎn)剛度方法的GPU計(jì)算方法進(jìn)行剛度矩陣計(jì)算速度得到了很大的提升，是文獻(xiàn)[10]中算法的速度的（3～5）倍。為了更直觀的表達(dá)，給出了兩者測(cè)試時(shí)間的柱狀圖，如圖4所示。

表1 剛度矩陣組裝測(cè)試結(jié)果Tab.1 Test Results of Stiffness Matrix Assembly

圖4 剛度矩陣裝配效率比較柱狀圖Fig.4 Histogram of Algorithm Efficiency Comparison of Stiffness Matrix Assemble

4 數(shù)值算例

為了驗(yàn)證所述算法在結(jié)構(gòu)優(yōu)化設(shè)計(jì)中的效率，以SIMP法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化為例，對(duì)兩端簡(jiǎn)支梁、懸臂梁兩個(gè)算例分別進(jìn)行優(yōu)化。實(shí)現(xiàn)平臺(tái)及硬件參數(shù)，如表2所示。

表2 算法實(shí)現(xiàn)平臺(tái)及參數(shù)Tab.2 Algorithm Implementation Platform and Parameters

4.1 兩端簡(jiǎn)支梁

對(duì)于兩端簡(jiǎn)支的二維梁結(jié)構(gòu)拓?fù)鋬?yōu)化問題，如圖5所示。根據(jù)對(duì)稱性，可取其中一半進(jìn)行優(yōu)化處理。問題的設(shè)計(jì)域?yàn)?00×100的矩形區(qū)域，左端邊界x方向固定，右下角y方向固定，左上角加載方向向下的載荷F，材料需用體積為40%。分別采用文獻(xiàn)[10]和提出的GPU計(jì)算方法以柔度為優(yōu)化目標(biāo)，對(duì)結(jié)構(gòu)進(jìn)行拓?fù)鋬?yōu)化計(jì)算，得到優(yōu)化結(jié)果，如圖6所示。由圖可知，所述的GPU計(jì)算方法可以得到與文獻(xiàn)[10]算法相同的優(yōu)化結(jié)果，證明了算法的有效性。

圖5 兩端簡(jiǎn)支的二維梁和問題簡(jiǎn)化圖Fig.5 Two Dimentional Simplified Beam with Two Ends and Simplified Graph

圖6 簡(jiǎn)支梁優(yōu)化結(jié)果對(duì)比Fig.6 Comparison of Optimization Results of Simply Supported Beam

將算法和文獻(xiàn)[10]的優(yōu)化時(shí)間進(jìn)行比較，結(jié)果，如表3所示。可以看出，算法計(jì)算速度得到了提升，是文獻(xiàn)[10]中算法的速度的（1.5～2）倍，說明了算法的高效性。為了更直觀的表達(dá)，給出了兩算法計(jì)算時(shí)間的柱狀圖，如圖7所示。

表3 簡(jiǎn)支梁算法效率比較Tab.3 Algorithm Efficiency Comparison of Simply Supported Beam

圖7 簡(jiǎn)支梁算法效率比較柱狀圖Fig.7 Histogram of Algorithm Efficiency Comparison of Simply Supported Beam

4.2 懸臂梁

對(duì)懸臂梁拓?fù)鋬?yōu)化問題，如圖8所示。問題的設(shè)計(jì)域?yàn)?60×100的矩形區(qū)域，左端邊界固定，右下角加載方向向下的載荷F，材料需用體積為40%。

圖8 懸臂梁?jiǎn)栴}簡(jiǎn)圖Fig.8 Simplified Graph of Cantilever Beam

分別采用文獻(xiàn)[10]和提出的GPU計(jì)算方法，柔度為優(yōu)化目標(biāo)，對(duì)結(jié)構(gòu)進(jìn)行拓?fù)鋬?yōu)化計(jì)算，得到優(yōu)化結(jié)果，如圖9所示?？梢钥闯?，所述的GPU計(jì)算方法可以得到與文獻(xiàn)[10]算法一致的優(yōu)化結(jié)果，證明算法是有效的。

圖9 懸臂梁優(yōu)化結(jié)果對(duì)比Fig.9 Comparison of Optimization Results of Cantilever Beam

將算法與文獻(xiàn)[10]的優(yōu)化時(shí)間進(jìn)行比較結(jié)果，如表4所示?？梢钥闯觯惴ㄓ?jì)算速度得到了提升，約是文獻(xiàn)[10]計(jì)算速度的1.5倍，說明了算法的高效性。為了更直觀的表達(dá)，給出了兩算法計(jì)算時(shí)間的柱狀圖，如圖10所示。

表4 懸臂梁算法效率比較Tab.4 Algorithm Efficiency Comparison of Cantilever Beam

圖10 懸臂梁算法效率比較柱狀圖Fig.10 Histogram of Algorithm Efficiency Comparison of Cantilever Beam

5 結(jié)論

研究分析了基于SIMP法的二維結(jié)構(gòu)拓?fù)鋬?yōu)化方法，利用GPU并行計(jì)算對(duì)算法進(jìn)行了改進(jìn)，通過算例分析，得到了以下結(jié)論：（1）通過對(duì)二維結(jié)構(gòu)化網(wǎng)格有限元分析過程的研究，提出了并行化有限元分析的計(jì)算方法，通過等效節(jié)點(diǎn)剛度的有限元分析法，加速了結(jié)構(gòu)剛度矩陣的計(jì)算過程，相比現(xiàn)有的計(jì)算方法，計(jì)算速度得到了很大提升；（2）對(duì)二維結(jié)構(gòu)拓?fù)鋬?yōu)化計(jì)算過程進(jìn)行了GPU層面的實(shí)現(xiàn)，并通過算例證實(shí)了方法的有效性，得到了比較理想的結(jié)構(gòu)優(yōu)化結(jié)果；（3）將所述的計(jì)算方法與文獻(xiàn)[10]的計(jì)算方法進(jìn)行了效率的比較，證明了GPU計(jì)算方法在解決拓?fù)鋬?yōu)化問題的高效性。