第一作者胡智強(qiáng)男，博士生，1984年12月生

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

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

胡智強(qiáng)1,馬海濤1,2

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

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

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

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

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

中圖分類(lèi)號(hào):TU311.3；O342

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

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

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

1問(wèn)題提出

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

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

(1)

(2)

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

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

(3)

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

(4)

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

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

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

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

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

(5)

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

(6)

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

(7)

(8)

(9)

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

對(duì)動(dòng)力問(wèn)題，將響應(yīng)函數(shù)取成一般積分形式，即

(10)

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

(11)

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

(12)

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

(13)

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

(14)

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

(15)

(16)

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

(17)

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

(18)

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

(19)

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

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

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

(20)

(21)

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

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

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

(22)

(23)

(24)

式中：

(25)

(26)

(27)

(1≤i≤N-1)

(28)

(29)

(30)

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

3數(shù)值算例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4結(jié)論

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

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

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

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

參考文獻(xiàn)

[1]陳太聰,韓大建,蘇成. 參數(shù)靈敏度分析的神經(jīng)網(wǎng)絡(luò)方法及其工程應(yīng)用[J]. 計(jì)算力學(xué)學(xué)報(bào), 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]陳鋼,趙國(guó)忠,顧元憲. 聲場(chǎng)-結(jié)構(gòu)耦合系統(tǒng)靈敏度分析及優(yōu)化設(shè)計(jì)研究[J]. 振動(dòng)與沖擊, 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)動(dòng)態(tài)載荷識(shí)別研究[J]. 振動(dòng)與沖擊, 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)振動(dòng)特性靈敏度分析 [J]. 振動(dòng)與沖擊, 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)激勵(lì)下結(jié)構(gòu)隨機(jī)振動(dòng)時(shí)域分析法[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時(shí)刻位移、速度響應(yīng)計(jì)算式如下

(A-3)

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

(A-4)

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