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

通信作者馬海濤男，教授，博士生導師，1962年生

結構動力響應靈敏度分析伴隨法一致性問題研究

胡智強1,馬海濤1,2

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

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

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

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

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

中圖分類號:TU311.3；O342

文獻標志碼：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ōu)化、系統(tǒng)辨識等研究中具有重要作用[1-4]。據計算方式不同，靈敏度分析方法主要分為直接法、伴隨法(伴隨變量方法)。對參數變量數目較多、所需考慮的響應量較少的一類問題，伴隨法的計算效率高于直接法，因此應用更廣[5-8]。Keulen等[9]依據精度與一致性、計算量及實現(xiàn)難度三項準則，對線性問題的不同靈敏度分析算法進行總結與評述。

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

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

1問題提出

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

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

(1)

(2)

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

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

(3)

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

(4)

式中：ε為計算靈敏度與精確解差別；εc為計算靈敏度與目標離散數值解差別。

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

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

2結構動力響應計算與靈敏度分析

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

(5)

2.1動力方程的數值解法

(6)

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

(7)

(8)

(9)

2.2動力響應函數及靈敏度分析的伴隨法

對動力問題，將響應函數取成一般積分形式，即

(10)

為計算響應函數Φ的靈敏度，引入n維伴隨變量向量λ(t)，n為結構自由度數，并定義函數Η為

(11)

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

(12)

為消除上式中含響應偏導數的各項，引入關于伴隨向量λ(t)的微分方程為

(13)

設初始響應與結構參數無關，據式(12)后兩項引入λ(t)應滿足的終值條件為

(14)

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

(15)

(16)

求解式(15)可確定伴隨變量Λ(t)，再由式(12)得響應函數Φ對結構參數的靈敏度為

(17)

2.3靈敏度分析伴隨法應用與一致性問題

(18)

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

(19)

式中：wi為數值積分權重，分別為w0=0.5Δt，w1=…=wN-1=Δt及wN=0.5Δt。

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

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

(20)

(21)

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

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

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

(22)

(23)

(24)

式中：

(25)

(26)

(27)

(1≤i≤N-1)

(28)

(29)

(30)

由時域逐步積分法獲得各時刻點動力響應后，對響應函數計算各時刻點的伴隨向量，再用上式計算對不同結構參數的靈敏度信息。

3數值算例

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

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

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

考慮三種不同形式響應函數，即

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

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

圖2 單自由度振子響應函數Φ1及靈敏度相對誤差Fig.2RelativeerrorsofresponsefunctionΦ1anditssensitivity圖3 單自由度振子響應函數Φ2及靈敏度相對誤差Fig.3RelativeerrorsofresponsefunctionΦ2anditssensitivity圖4 單自由度振子響應函數Φ3及靈敏度相對誤差Fig.4RelativeerrorsofresponsefunctionΦ3anditssensitivity

3.2三層剪切型結構

圖5　三層剪切型結構 Fig.5 Three-storey shear structure

取兩種形式的響應函數，即

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

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

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

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

4結論

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

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

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

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

參考文獻

[1]陳太聰,韓大建,蘇成. 參數靈敏度分析的神經網絡方法及其工程應用[J]. 計算力學學報, 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]陳鋼,趙國忠,顧元憲. 聲場-結構耦合系統(tǒng)靈敏度分析及優(yōu)化設計研究[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]毛玉明,郭杏林,趙巖,等. 基于靈敏度分析的結構動態(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)激勵下結構隨機振動時域分析法[J]. 工程力學, 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)

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

(A-3)

式中： Ai,k為系數矩陣，可表示為

(A-4)

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