結構優(yōu)化半解析靈敏度分析的改進算法

陳飆松， 張力丹， 魯一南， 程耿東

（1．大連理工大學 工業(yè)裝備結構分析國家重點實驗室，大連116024；2．大連理工大學 江蘇研究院，常州213164）

1 引 言

高效的靈敏度分析是提高優(yōu)化效率的關鍵［1］，已有的靈敏度求解方法包括全局差分法、解析法與半解析法。全局差分法［2］只要求有限元程序提供相鄰設計的結構性能，不牽涉這些設計的有限元模型和性能的計算方法和程序細節(jié)，可把有限元軟件作為黑箱，十分簡單，但每做一次敏度分析，針對n個設計變量要做n＋1次完整結構分析，效率很低［3］。解析法［4，5］通過推導得到靈敏度公式，計算效率和精度都較高，但推導這些公式并編程通常復雜且困難，尤其在形狀優(yōu)化中，形狀設計變量與有限元單元列式的關系不僅復雜而且依賴于單元類型，推導和編程均不具通用性［6］，不僅給程序開發(fā)造成困難，同時對于單元豐富的商用有限元程序，用戶也無法進行二次開發(fā)［7］。為此，程耿東等［8］提出了半解析靈敏度分析方法，為將有限元和結構優(yōu)化高效結合提供了途徑。該方法通用性強，效率高，廣泛地應用于多種結構優(yōu)化問題中。如Kiendl等［9］將半解析靈敏度分析方法及敏度加權應用于殼結構的等幾何形狀優(yōu)化分析中；陳敏等［10］在山地車車架結構的輕量化設計上使用了半解析靈敏度分析方法。傳統(tǒng)半解析靈敏度分析方法基于單元級別，需提取攝動前后兩組單元矩陣，影響計算精度和效率。

因此，本文對算法提出如下改進。（1）從總體角度考慮半解析靈敏度分析公式的計算方法，避免了攝動前單元剛度陣和質量陣等的提取，減少編程工作量和節(jié)省矩陣讀取時間的同時，減少了舍入誤差。（2）將算法推廣到應力、自振頻率、屈曲臨界荷載以及瞬態(tài)響應等力學量的靈敏度計算。以梁單元與殼單元為例，構建測試算例。結果表明，該方法可適用于多種分析類型，設計變量步長有更大的數值穩(wěn)定區(qū)域，可為采用非結構化網格建模結構的形狀優(yōu)化提供新途徑。（3）從總體角度推導的半解析靈敏度公式，在結合自主有限元程序與商用有限元程序進行實現時均具有較大優(yōu)勢。本文的公式推導以及所組織的算例測試主要針對基于自主有限元程序的半解析靈敏度分析方法，同時給出此方法結合商用軟件的應用，值得說明的是，對于靜力靈敏度問題，結合商用軟件使用時，無需提取任何矩陣。

2 傳統(tǒng)半解析靈敏度分析方法

半解析靈敏度方法將荷載和剛度陣等對于設計變量的導數用差分替代。靜力問題位移靈敏度分析公式為

式中u為結構整體位移向量，d為設計變量，F為施加在結構上的外荷載，K為結構總剛，Q為擬荷載向量。對于荷載與設計變量無關的情況，擬荷載可表示為

其中Ke為單元剛度矩陣，ue為單元位移向量。

對于靜力問題，式（2）需要提取攝動前后兩組單元剛度矩陣，再進行差分求解；而對于特征值問題，還需提取兩組單元質量陣或幾何剛度陣，從而增加編程工作量和數據輸入／輸出的工作量。在結合現有的商用軟件進行半解析靈敏度分析時，從軟件讀出的單元剛度陣往往只有固定位數的有效數字，按上述方法計算時有效位數有時會損失很多，嚴重影響計算精度。

3 改進半解析靈敏度分析方法

3．1 靜力問題

將傳統(tǒng)方法進行改進設計，使其程序實現過程更為簡便，首先對靜力問題進行推導。從總體角度考慮靜力問題位移靈敏度求解。對于結構的靜力分析問題：

結構位移的半解析法靈敏度分析公式如下，

由式（5）可知，F為外荷載。式（4～6）直接基于總剛進行求解，更適用于單剛提取困難、提取精度較低且提取過程耗時的商用有限元程序。在基于自主有限元程序進行半解析靈敏度分析時，因單剛較易獲得，可將式（6）進行如下改進，即F′分解為單元一級進行計算并累加于總體，從而避免總剛的集成。

該方法從總體角度入手，在基于自主有限元程序進行實現時，使用改進公式（4～7），避免初始剛度矩陣的提取，以及總體剛度陣的集成，并可減小剛度矩陣相近數值相減造成的有效數字損失；在結合商用軟件實現時，使用總體公式（4～6），F′的計算可將初始所有節(jié)點的位移u作為結構的位移全約束，通過商用軟件提取外荷載得到，從而無需讀出任何剛度陣，同時也避免了單元剛度矩陣讀取的有效數字位數對計算精度的影響。

在得到靜力位移靈敏度的條件下，還可將其應用于靜應力靈敏度的求解，結構的應力為

其基于總體的半解析靈敏度方法的計算可以近似為

式中 D為整體彈性矩陣，B為整體應變矩陣，u為整體結構位移，d為設計變量，將位移乘入即可得到

其中，De為單元彈性矩陣，Be為單元應變矩陣，ue為單元位移，由式（12）可知，σ為初始結構的應力。根據式（10），改進方法計算靜應力靈敏度時，無需提取原設計的單元彈性矩陣與單元應變矩陣，而用初始應力代替?；谧灾饔邢拊绦蜻M行的半解析靈敏度分析，σ′的計算可先形成攝動后的De·Be（d＋Δd），將攝動前的結構位移分解到單元上，代入攝動后的應力求解流程中，對于σ－的計算，可先求解出位移相對于設計變量的靈敏度，再將其作為結構位移并分解到單元上回代于攝動前的應力計算流程，得到的應力即為σ－。在結合商用軟件使用時，可將初始結構位移及求得的位移靈敏度分別作為結構的位移約束，通過商用有限元軟件直接計算應力得到σ′與σ－，無需提取任何矩陣。

3．2 特征值問題

對于自振分析非重特征值靈敏度問題［11］，基于總體的半解析式為

將總體特征向量乘入可得

將式（28）改寫為總體差分格式，并將加速度、速度及位移乘入差分近似的公式中，可得到f－為

考慮自振分析的特性，對質量陣歸一化：

則有

式（20）更適用于總剛讀取效率高于單剛的商用有限元軟件。基于自主有限元程序進行半解析靈敏度分析時，為避免總體矩陣的集成，可將式（16，17）進行如下改進，

則式（20）可寫為

同理，對于屈曲分析非重特征值靈敏度問題［12］可得

以及

3．3 瞬態(tài)問題

瞬態(tài)響應方程［13］的表達式為Mu¨＋Cu·＋Ku＝f（t） （26）式中 M為總體質量陣，C為總體阻尼陣，K為總體剛度矩陣，u¨為總體加速度，u·為總體速度，u為總體位移，f（t）為施加在結構上的時變荷載，對于荷載與設計變量無關的情況，瞬態(tài)問題靈敏度求解公式［14］為

其中，

由式（31），f為施加在結構上的時變荷載，即瞬態(tài)問題半解析法公式的右端項可寫為

在計算過程中，先計算式（27）的右端項f－，并將其作為新的時變載荷，代入瞬態(tài)問題求解方法［15］中，如 NewMark法［16］和 Willson－θ 法［17］等，求解初值為0的瞬態(tài)問題，得到的位移、速度和加速度值即為結構的位移、速度和加速度靈敏度值。

4 算例測試

4．1 測試模型

現構建兩種有限元模型，對靜力分析的位移、應力靈敏度、自振分析特征值靈敏度以及瞬態(tài)分析位移靈敏度進行算例測試。如圖1所示，框架尺寸h＝6m，b＝4m，梁矩形截面為0．15m×0．06m，彈性模量E＝2×1011Pa，泊松比ν＝0．3，密度ρ＝7．8×103kg／m3。采用歐拉伯努利梁單元進行建模。每段梁結構劃分20個梁單元，共300個梁單元。該模型用于以下分析類型。

（1）靜位移。底部固定，在其頂部左端施加沿x軸正方向大小為0．1N的集中荷載，取頂部右端點x方向位移對于b的靈敏度，

（2）靜應力。邊界條件同靜位移，取左側柱根部單元高斯點最大應力對于b的靈敏度。

圖1 框架模型Fig．1 Model of frame

（3）自振。底部固定，取基頻對于b的靈敏度。

（4）瞬態(tài)。底部固定，荷載位置與方向同靜位移，取頂部右端點x方向對于b的靈敏度。

如圖2所示，橢球殼結構，長軸a＝4m，短軸b＝3m，彈性模量E ＝2×1011Pa，泊松比ν＝0．3，密度ρ＝7．8×103kg／m3。采用柯西霍夫理論四邊形殼單元進行建模，網格剖分如圖2所示，共300個單元。該模型將用于以下分析類型。

（1）靜位移。底部固定，在其球殼頂部中央節(jié)點，施加沿z軸負方向大小為0．1N的集中荷載，取頂部z方向位移對于b的靈敏度。

（2）靜應力。邊界條件同靜位移，取頂部單元高斯點最大應力對于b的靈敏度。

（3）自振。底部固定，取基頻對于b的靈敏度。

（4）瞬態(tài)。施加于結構上的荷載位置與方向同靜位移，取頂部z方向位移對于b的靈敏度。

圖2 橢球殼模型Fig．2 Model of elliptical shell

4．2 測試結果

基于文獻［3］對于相對靈敏度的定義，設全局差分法穩(wěn)定值為ΔF0GFD，分別定義GFD方法、SA方法與MSA方法的相對靈敏度為

其符號說明如下，

GFD（global finite difference）為全局差分法；

SA（semi－analysis）為將單元剛度陣直接進行差分的傳統(tǒng)半解析方法；

MSA（modified semi－analysis）為基于總體的改進半解析方法，此處采用自主有限元程序進行實現（攝動后部分采用單元一級累加于總體，避免總體矩陣的集成）。

圖3～圖6分別為靜位移、靜應力、自振特征值和瞬態(tài)位移的相對靈敏度隨攝動步長與設計變量之比的負對數δ＝－log（ΔL／L）的變化規(guī)律。梁單元建模的框架結構（簡稱梁結構）測試結果繪于圖（a），用表示。殼單元建模的板殼結構（簡稱板殼結構）測試結果繪于圖（b），用表示。

從四種靈敏度分析結果可以看出，半解析方法（MSA與SA）在梁結構靜應力靈敏度的求解上，數值穩(wěn)定性較為突出，如圖4（a）所示，GFD穩(wěn)定范圍為1e－4～1e－7，SA方法為1e－5～1e－10，MSA方法為1e－5～1e－11。與SA相比，MSA方法具有較長的數值穩(wěn)定區(qū)間，但無法縮小有限單元剛體轉動造成的誤差［17］，因此不適用于大攝動步長。

從梁結構與板殼結構（非類梁結構）的結果對比可以看出，對于板殼結構，三種方法的數值穩(wěn)定范圍均較長且差異不大，而對于梁結構，三者具有較大差別。從整體趨勢來看，GFD法的穩(wěn)定值出現于大攝動步長，而對于SA或MSA方法，則在較小攝動步長達到穩(wěn)定，尤其對于MSA方法，緩解了相近數值相減造成的誤差，在步長更小的位置仍能保證數值精度。對于形狀設計變量問題，該特性能保證攝動前后網格的一致性，對于優(yōu)化迭代計算較為有利。

圖3 靜力位移相對靈敏度隨攝動步長的變化Fig．3 Static displacement scaled sensitivity fluctuated with the perturbation step length

圖4 靜應力相對靈敏度隨攝動步長的變化Fig．4 Static stress scaled sensitivity fluctuated with the perturbation step length

圖5 自振特征值相對靈敏度隨攝動步長的變化Fig．5 Eigenvalue scaled sensitivity of natural frequency analysis fluctuated with the perturbation step length

圖6 瞬態(tài)位移相對靈敏度隨攝動步長的變化Fig．6 Displacement scaled sensitivity of transient analysis fluctuated with the perturbation step length

5 結 論

（1）本文將傳統(tǒng)基于單元級別的靈敏度求解公式轉為結構層次進行考慮，得到無需提取初矩陣的程序實現方法。

（2）在靜力問題測試的基礎上，將半解析方法應用于自振和瞬態(tài)分析等多種分析類型，測試算例表明，該方法在多種分析類型上都可以得到針對攝動步長較大的穩(wěn)定區(qū)域。

（3）測試算例為改進方法在自主有限元程序上的應用。對于結合商用有限元軟件進行靈敏度分析，應將攝動后部分直接在結構層次上進行求解。相較傳統(tǒng)讀出攝動前后兩組單剛的方法，在很大程度上降低了矩陣讀取時間，提高了總體計算效率，且能避免商用軟件提供的單剛具有較少有效數字位數對計算精度造成的嚴重影響。