結(jié)構(gòu)動力響應(yīng)靈敏度分析伴隨法一致性問題研究

第一作者胡智強男，博士生，1984年12月生

通信作者馬海濤男，教授，博士生導(dǎo)師，1962年生

結(jié)構(gòu)動力響應(yīng)靈敏度分析伴隨法一致性問題研究

胡智強1,馬海濤1,2

(1.華南理工大學(xué)土木與交通學(xué)院亞熱帶建筑科學(xué)國家重點實驗室，廣州510640；2.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連116023)

摘要：研究結(jié)構(gòu)瞬態(tài)動力響應(yīng)靈敏度分析伴隨法可能存在的一致性問題，具體考慮先微分后離散與先離散后微分兩種敏度分析方法的計算精度、收斂速度及一致性等?；趧恿憫?yīng)分析的時域顯式法基本思想，以更簡潔方式推導(dǎo)先離散后微分的伴隨法計算列式。結(jié)果表明，先微分后離散伴隨法一致性問題是由計算中所用動力響應(yīng)結(jié)果僅在離散時間點滿足運動方程產(chǎn)生的，存在的一致性問題不影響該方法的可靠性及應(yīng)用。

關(guān)鍵詞：靈敏度分析；伴隨法；一致性；瞬態(tài)動力響應(yīng)

基金項目：國家自然科學(xué)基金(11372004)；工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室(大連理工大學(xué))(GZ1305)

收稿日期：2014-07-16修改稿收到日期：2014-09-25

中圖分類號:TU311.3；O342

文獻標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2015.20.028

Abstract:The inconsistency issue of adjoint variable methods (AVMs) for sensitivity analysis of transient dynamic responses was investigated. The differentiate-then-discretize and discretize-then-differentiate approaches were considered, focusing on their computational accuracy, convergence rates and result consistency. Based on the basic idea of the explicit time-domain method for dynamic analysis, a concise discretize-then-differentiate AVM formulation was presented. It is found that the inconsistency of the differentiate-then-discretize approach is caused by the fact that numerical solutions for dynamic responses satisfy equations of motion only at integration points in the time domain. However, despite this consistency problem, this approach is still reliable for sensitivity analysis of dynamic responses.

On the consistency issue of adjoint methods for sensitivity analysis of dynamic responses

HUZhi-qiang1,MAHai-tao1,2(1. State Key Laboratory of Subtropical Building Science, School of Civil Engineering and Transportation, South China University of Technology, Guangzhou 510640, China;2. State Key Laboratory of Structural Analysis of Industrial Equipment, Dalian University of Technology, Dalian 116023, China)

Key words:sensitivity analysis; adjoint variable method; consistency; transient dynamic response

靈敏度分析用于定量預(yù)測結(jié)構(gòu)參數(shù)變化對響應(yīng)影響，在結(jié)構(gòu)優(yōu)化、系統(tǒng)辨識等研究中具有重要作用[1-4]。據(jù)計算方式不同，靈敏度分析方法主要分為直接法、伴隨法(伴隨變量方法)。對參數(shù)變量數(shù)目較多、所需考慮的響應(yīng)量較少的一類問題，伴隨法的計算效率高于直接法，因此應(yīng)用更廣[5-8]。Keulen等[9]依據(jù)精度與一致性、計算量及實現(xiàn)難度三項準(zhǔn)則，對線性問題的不同靈敏度分析算法進行總結(jié)與評述。

靈敏度方程推導(dǎo)會涉及對結(jié)構(gòu)參數(shù)的微分運算。對瞬態(tài)動力問題，響應(yīng)計算涉及時間域的離散化，因此靈敏度分析方法一般存在先微分后離散或先離散后微分兩種方式。對直接法，采用此兩種方式所得動力響應(yīng)靈敏度結(jié)果相同；但對伴隨法，微分、離散的先后次序會影響靈敏度分析結(jié)果。Jensen等[10]用單自由度模型研究先微分后離散、先離散后微分的兩種靈敏度分析伴隨法認(rèn)為，先微分后離散伴隨法會給出不一致的靈敏度結(jié)果，相應(yīng)誤差須通過適當(dāng)選擇時間步長及數(shù)值求積公式控制，而徹底解決此問題可用先離散后微分算法。

針對一致性問題，本文全面研究動力響應(yīng)靈敏度分析伴隨法的基本問題，包括一致性及收斂速度等。即總結(jié)一致性問題提法，簡要歸納瞬態(tài)動力響應(yīng)分析的Newmark-β方法及先微分后離散的靈敏度分析伴隨法，明確一致性問題根源。并基于動力響應(yīng)計算時域顯式法基本思想，給出先離散后微分的伴隨法靈敏度計算列式，針對單自由度、多自由度問題通過兩算例研究兩種伴隨法靈敏度計算精度及收斂性。

1問題提出

靈敏度分析算法選擇應(yīng)依據(jù)精度與一致性、計算量及實現(xiàn)難度三項準(zhǔn)則。文獻[9]將靈敏度分析結(jié)果誤差分為兩類，一類為離散系統(tǒng)的數(shù)值解與精確值偏差，用精度表示；一類為計算靈敏度與相對精確數(shù)值解的偏差，用不一致性表示。

為簡便，考慮離散結(jié)構(gòu)系統(tǒng)動力問題。以動力響應(yīng)Φ為例，分別用Φa(t,b)及Φd(t,b,Δt)表示精確解、數(shù)值解。其中b為某結(jié)構(gòu)參數(shù)；t為時間；Δt為時間步長。響應(yīng)靈敏度可定義為

(1)

(2)

由于數(shù)值解存在離散誤差，解析解Φa(t,b)與數(shù)值解Φd(t,b,Δt)之間存在差別，由式(1)、(2)給出的靈敏度結(jié)果會不同，誤差由時間離散誤差引起。為考察靈敏度分析方法的準(zhǔn)確性、一致性，按文獻[9-10]，在相同時間離散情況下，若某個靈敏度分析方法給出的結(jié)果為gs(t,b,Δt)，則可定義兩種誤差為

ε=gs(t,b,Δt)-ga(t,b)

(3)

εc=gs(t,b,Δt)-gd(t,b,Δt)

(4)

式中：ε為計算靈敏度與精確解差別；εc為計算靈敏度與目標(biāo)離散數(shù)值解差別。

顯然，誤差ε反映計算精度，其數(shù)值大小反映敏度分析方法的準(zhǔn)確性；而誤差εc為該靈敏度分析方法計算結(jié)果與所得數(shù)值解偏差，稱為一致性誤差。

若某靈敏度分析方法給出的結(jié)果存在較大一致性誤差，該結(jié)果無法準(zhǔn)確反映參數(shù)變化影響，會造成后續(xù)計算效率降低及最終結(jié)果偏差。雖常規(guī)基于先微分后離散的伴隨法計算結(jié)果有一致性誤差，但仍未見理論解釋。為避免不一致靈敏度分析方法產(chǎn)生的影響，文獻[7]在瞬態(tài)優(yōu)化中選擇用一致靈敏度分析算法，而文獻[10]則認(rèn)為可選先離散后微分的伴隨法。對此，本文研究一致性問題產(chǎn)生的根源及對結(jié)果影響，加深對該問題的理解、認(rèn)識，有助于靈敏度分析方法的理論、應(yīng)用研究。

2結(jié)構(gòu)動力響應(yīng)計算與靈敏度分析

線性離散結(jié)構(gòu)動力系統(tǒng)運動方程可表示為

(5)

2.1動力方程的數(shù)值解法

(6)

用Newmark-β法[11]，其逐步積分列式為

(7)

(8)

(9)

2.2動力響應(yīng)函數(shù)及靈敏度分析的伴隨法

對動力問題，將響應(yīng)函數(shù)取成一般積分形式，即

(10)

為計算響應(yīng)函數(shù)Φ的靈敏度，引入n維伴隨變量向量λ(t)，n為結(jié)構(gòu)自由度數(shù)，并定義函數(shù)Η為

(11)

由于式(5)在任意時刻均成立，響應(yīng)函數(shù)Φ恒等于Η。因此，兩函數(shù)偏導(dǎo)數(shù)也應(yīng)恒等。對式(11)微分，再進行分部積分，整理后得

(12)

為消除上式中含響應(yīng)偏導(dǎo)數(shù)的各項，引入關(guān)于伴隨向量λ(t)的微分方程為

(13)

設(shè)初始響應(yīng)與結(jié)構(gòu)參數(shù)無關(guān)，據(jù)式(12)后兩項引入λ(t)應(yīng)滿足的終值條件為

(14)

采用變量代換s=T-t，將式(13)、(14)的終值問題轉(zhuǎn)化成初值問題，即

(15)

(16)

求解式(15)可確定伴隨變量Λ(t)，再由式(12)得響應(yīng)函數(shù)Φ對結(jié)構(gòu)參數(shù)的靈敏度為

(17)

2.3靈敏度分析伴隨法應(yīng)用與一致性問題

(18)

由于式(18)用時刻點t=T-si處動力響應(yīng)，故采用Newmark-β法即可獲得不同時刻點伴隨向量(Λi)。此時，用梯形積分公式[12]可計算式(17)積分，即

(19)

式中：wi為數(shù)值積分權(quán)重，分別為w0=0.5Δt，w1=…=wN-1=Δt及wN=0.5Δt。

為獲得靈敏度計算公式，對連續(xù)響應(yīng)函數(shù)進行微分運算再推導(dǎo)獲得伴隨方程，方可在時域內(nèi)離散求解方程。因此將其稱為先微分后離散的伴隨法。

需要注意的是，按解析推導(dǎo)過程要求，伴隨向量計算采用的動力響應(yīng)結(jié)果應(yīng)在整個時域內(nèi)滿足平衡方程。但式(18)中所用動力響應(yīng)值僅在離散時間點處滿足式(6)，因此響應(yīng)函數(shù)Φ并不恒等于Η，以此為基礎(chǔ)所得伴隨向量及再由式(19)計算的靈敏度結(jié)果會存在偏差。將式(11)兩端同時對b求偏導(dǎo)，得

(20)

(21)

至此，以推導(dǎo)的靈敏度分析方法為基礎(chǔ)，采用數(shù)值計算方式可實現(xiàn)先微分后離散伴隨法；但該算法存在一致性問題，由于響應(yīng)的數(shù)值解無法滿足解析的假設(shè)條件，計算所得靈敏度結(jié)果會存在一致性誤差。本文將通過兩算例具體研究此誤差對靈敏度計算精度及收斂速度影響。

2.4先離散后微分的靈敏度分析伴隨法

文獻[10]認(rèn)為先微分后離散伴隨法會給出不一致的靈敏度結(jié)果，而先離散后微分的伴隨法則不存在一致性問題。本文基于動力響應(yīng)計算的時域顯式法[13]基本思想，以更簡潔方式給出先離散后微分的伴隨法計算列式。采用逐步積分法計算獲得各時刻點動力響應(yīng)后，可用梯形積分公式計算式(10)的響應(yīng)函數(shù)，即

(22)

(23)

(24)

式中：

(25)

(26)

(27)

(1≤i≤N-1)

(28)

(29)

(30)

由時域逐步積分法獲得各時刻點動力響應(yīng)后，對響應(yīng)函數(shù)計算各時刻點的伴隨向量，再用上式計算對不同結(jié)構(gòu)參數(shù)的靈敏度信息。

3數(shù)值算例

利用單、多自由度兩個算例，對比、討論兩種伴隨法的計算精度及收斂速度，說明一致性問題的存在不影響先微分后離散伴隨法的可靠性及應(yīng)用。

3.1單自由度系統(tǒng)

圖1　單自由度振子模型 Fig.1 The model of single dof oscillator

考慮三種不同形式響應(yīng)函數(shù)，即

對比響應(yīng)函數(shù)與靈敏度計算結(jié)果知，在相同時間離散情況下二者結(jié)果相近。隨時間步長Δt減小，兩種伴隨法靈敏度誤差均逐漸減小，且相對誤差與時間步長Δt的平方基本呈線性關(guān)系，即兩種伴隨法均具有二階收斂速度，并與計算響應(yīng)函數(shù)用Newmark-β法、梯形積分公式一致。

值得注意的是，雖先微分后離散伴隨法存在一致性問題，但其靈敏度計算精度優(yōu)于響應(yīng)函數(shù)結(jié)果(圖3、圖4)，該現(xiàn)象為本算例特有，不具一般性。對本算例而言，一致性問題并未產(chǎn)生不良影響，反使靈敏度結(jié)果精度有一定程度提高。

圖2 單自由度振子響應(yīng)函數(shù)Φ1及靈敏度相對誤差Fig.2RelativeerrorsofresponsefunctionΦ1anditssensitivity圖3 單自由度振子響應(yīng)函數(shù)Φ2及靈敏度相對誤差Fig.3RelativeerrorsofresponsefunctionΦ2anditssensitivity圖4 單自由度振子響應(yīng)函數(shù)Φ3及靈敏度相對誤差Fig.4RelativeerrorsofresponsefunctionΦ3anditssensitivity

3.2三層剪切型結(jié)構(gòu)

圖5　三層剪切型結(jié)構(gòu) Fig.5 Three-storey shear structure

取兩種形式的響應(yīng)函數(shù)，即

考慮持時T=25T3，T3=0.607 5為第三階自振周期。取不同離散時間步長Δt=T3/q，q為每個周期離散時間步數(shù)。采用兩種伴隨法分別計算響應(yīng)函數(shù)對參數(shù)k2，k3的靈敏度，并與參考解比較，相對誤差結(jié)果見圖6～圖9。動力響應(yīng)及伴隨向量計算仍用Newmark-β法，參數(shù)取γ=0.5，β=0.25。為獲得準(zhǔn)確可靠參考解，除采用較小時間步長外，用理查森(Richardson)外插法[14]計算參考解。

圖6 三層剪切型結(jié)構(gòu)Φ1及對k2靈敏度相對誤差Fig.6RelativeerrorsofresponsefunctionΦ1anditssensitivitywithrespecttok2圖7 三層剪切型結(jié)構(gòu)Φ1及對k3靈敏度相對誤差Fig.7RelativeerrorsofresponsefunctionΦ1anditssensitivitywithrespecttok3圖8 三層剪切型結(jié)構(gòu)Φ2及對k2靈敏度相對誤差Fig.8RelativeerrorsofresponsefunctionΦ2anditssensitivitywithrespecttok2

在相同時間步長Δt下，由圖6、圖7看出，對響應(yīng)函數(shù)Φ1，一致性誤差導(dǎo)致先微分后離散伴隨法靈敏度計算精度降低；而對響應(yīng)函數(shù)Φ2，一致性誤差可使靈敏度分析精度略有提高(圖8)或降低(圖9)。因此，雖一致性誤差對靈敏度計算精度有影響，但影響作用不確定。另外，本算例靈敏度計算結(jié)果精度低于響應(yīng)函數(shù)計算精度。隨時間步數(shù)增加，兩種伴隨法計算精度均有提高，且與響應(yīng)函數(shù)的計算具有相同二階收斂速度。由此認(rèn)為，時間離散誤差仍為影響靈敏度計算精度的重要因素，一致性并不影響其計算收斂速度，該方法仍可靠。

圖9　三層剪切型結(jié)構(gòu)Φ 2及對k 3靈敏度相對誤差 Fig.9 Relative errors of response function Φ 2 and its sensitivity with respect to k 3

4結(jié)論

通過研究結(jié)構(gòu)動力響應(yīng)靈敏度分析伴隨法基本問題，包括算法的一致性、收斂速度及精度等，并由單、多自由度算例對比、討論先微分后離散與先離散后微分兩種伴隨法的計算精度及收斂速度，結(jié)論如下：

(1)由于時域逐步積分方法所得動力響應(yīng)解僅能保證在離散時間點處嚴(yán)格滿足動力平衡方程，因此會導(dǎo)致先微分后離散伴隨法出現(xiàn)一致性問題。

(2)時間離散誤差是影響兩種伴隨法計算精度的重要因素，且兩種方法二階收斂速度相同。對先微分后離散伴隨法，若能保證動力分析結(jié)果精度，一致性問題不會影響其可靠性。

(3)先微分后離散伴隨法的理論推導(dǎo)及計算列式較簡單，計算效率較高，雖存在一致性問題，但仍為較好的結(jié)構(gòu)動力響應(yīng)靈敏度分析方法。

參考文獻

[1]陳太聰,韓大建,蘇成. 參數(shù)靈敏度分析的神經(jīng)網(wǎng)絡(luò)方法及其工程應(yīng)用[J]. 計算力學(xué)學(xué)報, 2004, 21(6): 752-756.

CHEN Tai-cong, HAN Da-jian, SU Cheng. Neural network method in parameter sensitivity analysis and its application in engineering [J]. Chinese Journal of Computational Mechanics, 2004, 21(6): 752-756.

[2]陳鋼,趙國忠,顧元憲. 聲場-結(jié)構(gòu)耦合系統(tǒng)靈敏度分析及優(yōu)化設(shè)計研究[J]. 振動與沖擊, 2007, 26(4): 86-89.

CHEN Gang, ZHAO Guo-zhong, GU Yuan-xian. Sensitivity analysis and design optimization method for acoustic-structural coupled systems[J]. Journal of Vibration and Shock, 2007, 26(4): 86-89.

[3]毛玉明,郭杏林,趙巖,等. 基于靈敏度分析的結(jié)構(gòu)動態(tài)載荷識別研究[J]. 振動與沖擊, 2010, 29(10): 1-3.

MAO Yu-ming, GUO Xing-lin, ZHAO Yan, et al. Force identification based on sensitivity analysis method [J]. Journal of Vibration and Shock, 2010, 29(10): 1-3.

[4]趙杰,李峰,劉錄. 基于有限元的超高壓管線系統(tǒng)振動特性靈敏度分析 [J]. 振動與沖擊, 2014, 33(10): 148-151.

ZHAO Jie, LI Feng, LIU Lu. Vibration characteristic sensitivity analysis based on finite element for ultra-high pressure pipeline systems[J].Journal of Vibration and Shock, 2014, 33(10): 148-151.

[5]Dahl J, Jensen J S, Sigmund O. Topology optimization for transient wave propagation problems in one dimension[J]. Structural and Multidisciplinary Optimization, 2008, 36:585-595.

[6]Kang Z, Zhang X P, Jiang S G,et al. On topology optimization of damping layer in shell structures under harmonic excitations[J]. Structural and Multidisciplinary Optimization, 2012, 46: 51-67.

[7]Le C, Bruns T E, Tortorelli D A. Material microstructure optimization for linear elastodynamic energy wave management[J]. Journal of the Mechanics and Physics of Solids,2012, 60: 351-378.

[8]Kang B S, Park G J, Arora J S. A review of optimization of structures subjected to transient loads[J]. Structural and Multidisciplinary Optimization, 2006, 31: 81-95.

[9]Keulen F V, Haftka R T, Kim N H. Review of options for structural design sensitivity analysis,part 1: linear systems [J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194: 3213-3243.

[10]Jensen J S, Nakshatrala P B,Tortorelli D A. On the consistency of adjoint sensitivity analysis for structural optimization of linear dynamic problems[J]. Structural and Multidisciplinary Optimization, 2014, 49(5): 831-837.

[11]Bathe K J. Finite element procedures[M]. Prentice-Hall, Englewood Cliffs, New Jersey, 1996.

[12]Kreyszig E. Advanced engineering mathematics (7th edn) [M]. New York:Wiley, 1993.

[13]蘇成,徐瑞.非平穩(wěn)激勵下結(jié)構(gòu)隨機振動時域分析法[J]. 工程力學(xué), 2010, 27(12): 77-83.

SU Cheng, XU Rui. Random vibration analysis of structures subjected to non-stationary excitations by time domain method [J]. Engineering Mechanics, 2010, 27(12): 77-83.

[14]Cook R D, Malkus D S, Plesha M E, et al. Concepts and applications of finite element analysis(4th edn)[M]. New York: Wiley,2000.

附錄 A

Vi=TVi-1+Q1ui-1+Q2ui

(A-1)

式中：

(A-2)

重復(fù)應(yīng)用式(A-1)，可得第i時刻位移、速度響應(yīng)計算式如下

(A-3)

式中： Ai,k為系數(shù)矩陣，可表示為

(A-4)

因此，各時刻點處的位移、速度響應(yīng)可用對應(yīng)時刻點及之前全部時刻點加速度顯式表達(dá)成(A-3)形式。