黏性阻尼系統(tǒng)時(shí)域響應(yīng)靈敏度及其一致性研究1)

張磊 張 嚴(yán) 丁 喆

(武漢科技大學(xué)冶金裝備及其控制教育部重點(diǎn)實(shí)驗(yàn)室,武漢 430081)

(武漢科技大學(xué)機(jī)械傳動(dòng)與制造工程湖北省重點(diǎn)實(shí)驗(yàn)室,武漢 430081)

引言

靈敏度分析用于定量預(yù)測(cè)結(jié)構(gòu)參數(shù)變化對(duì)系統(tǒng)響應(yīng)或模型的影響,在模型修正[1-2]、拓?fù)鋬?yōu)化[3-6]、損傷識(shí)別[7]、可靠性分析[8]等研究工作中具有重要意義.工程結(jié)構(gòu)在實(shí)際服役時(shí)不可避免的受到復(fù)雜外界環(huán)境的振動(dòng)或沖擊影響,往往要求進(jìn)行時(shí)域動(dòng)態(tài)響應(yīng)的分析,且經(jīng)常會(huì)伴隨著噪聲以及結(jié)構(gòu)系統(tǒng)自身的振動(dòng),這都對(duì)結(jié)構(gòu)動(dòng)力設(shè)計(jì)及優(yōu)化問(wèn)題造成了困難.因此開(kāi)發(fā)出高效準(zhǔn)確的時(shí)域響應(yīng)靈敏度分析方法非常必要.根據(jù)計(jì)算方式不同,靈敏度計(jì)算方法可主要分為[9]:有限差分法、直接微分法和伴隨變量法.有限差分法是一種直接而簡(jiǎn)單的方法,但是計(jì)算效率較低,通常作為衡量靈敏度求解方法精度的參考解.已有的研究表明[10],當(dāng)設(shè)計(jì)變量數(shù)目大于有效約束數(shù)目時(shí),伴隨變量法比直接微分法計(jì)算效率更高.

動(dòng)力響應(yīng)的靈敏度求解對(duì)結(jié)構(gòu)動(dòng)態(tài)性能的分析及優(yōu)化非常重要.目前為止,眾多學(xué)者已經(jīng)對(duì)該領(lǐng)域進(jìn)行了廣泛的研究.Choi等[11]在波導(dǎo)濾波器的優(yōu)化設(shè)計(jì)中,利用增廣拉格朗日法、材料導(dǎo)數(shù)概念和伴隨變量法,系統(tǒng)的推導(dǎo)了頻域解析靈敏度公式.周大為等[12]針對(duì)黏性阻尼系統(tǒng)的頻響振幅靈敏度求解問(wèn)題,并結(jié)合廣義模態(tài)截?cái)鄶U(kuò)增方法,提出了一種具有更高精度的伴隨變量法.Zhao等[13]提出了一種自適應(yīng)混合展開(kāi)法來(lái)確定需要計(jì)算的低階特征頻率和特征模數(shù),然后再采用伴隨變量法進(jìn)行靈敏度分析,以克服簡(jiǎn)諧激勵(lì)下的拓?fù)鋬?yōu)化問(wèn)題計(jì)算量大的困難.

上述研究討論了伴隨變量法求解頻率響應(yīng)靈敏度的問(wèn)題.然而,實(shí)際工程由于工作環(huán)境復(fù)雜,結(jié)構(gòu)系統(tǒng)時(shí)域響應(yīng)的靈敏度分析同樣非常重要.Elesin等[14]在非線性光學(xué)器件的拓?fù)鋬?yōu)化結(jié)構(gòu)設(shè)計(jì)中,采用了伴隨變量法求解時(shí)域響應(yīng)靈敏度.Mello等[15]在減少微機(jī)電系統(tǒng)響應(yīng)時(shí)間的拓?fù)鋬?yōu)化中,對(duì)于時(shí)域相關(guān)目標(biāo)函數(shù)的靈敏度求解同樣采用了伴隨變量法.Zhao 和Wang[16]研究了涉及時(shí)域動(dòng)力響應(yīng)的拓?fù)鋬?yōu)化問(wèn)題,其中對(duì)時(shí)域響應(yīng)相關(guān)的目標(biāo)函數(shù)靈敏度分析采用了伴隨變量法.此外,還有一些學(xué)者也采用伴隨變量法來(lái)進(jìn)行時(shí)域響應(yīng)問(wèn)題相關(guān)的靈敏度分析[17-19],用以優(yōu)化結(jié)構(gòu)在外載荷作用下或應(yīng)力約束下的機(jī)械性能,如最小化能量耗散、動(dòng)柔度和結(jié)構(gòu)總體積等.時(shí)域響應(yīng)靈敏度求解同時(shí)涉及設(shè)計(jì)變量的微分計(jì)算和時(shí)間域的離散化,因此微分和離散的先后順序可能對(duì)時(shí)域響應(yīng)靈敏度的計(jì)算結(jié)果產(chǎn)生影響.但上述研究均利用先微分后離散的伴隨變量法實(shí)施.然而,Le等[20]在進(jìn)行材料微觀結(jié)構(gòu)的拓?fù)鋬?yōu)化時(shí)指出先微分后離散的伴隨變量法會(huì)導(dǎo)致一致性問(wèn)題,而采用先離散后微分的伴隨變量法可大幅降低一致性誤差.

一致性誤差用來(lái)描述某靈敏度計(jì)算方法的結(jié)果與數(shù)值參考解間的偏差,它與計(jì)算量和實(shí)現(xiàn)難度被共同視為評(píng)價(jià)靈敏度算法性能的三項(xiàng)準(zhǔn)則[21].Jensen等[22]指出基于Newmark 方法的先微分后離散伴隨變量法會(huì)產(chǎn)生一致性誤差,而解決此問(wèn)題可用先離散后微分方法.胡智強(qiáng)和馬海濤[23]針對(duì)伴隨變量法的一致性問(wèn)題進(jìn)行了更為細(xì)致的研究,發(fā)現(xiàn)若能保證動(dòng)力響應(yīng)分析結(jié)果的精度,先微分后離散法的一致性誤差并不會(huì)影響靈敏度結(jié)果的可靠性.Yun 和Youn[24]針對(duì)黏彈性阻尼系統(tǒng),提出了一種先離散后微分的時(shí)域響應(yīng)靈敏度求解方法,并將其應(yīng)用于黏彈性系統(tǒng)微結(jié)構(gòu)拓?fù)鋬?yōu)化.文獻(xiàn)[25]通過(guò)基于密度的多材料拓?fù)鋬?yōu)化方法,并結(jié)合先離散后微分的伴隨變量法進(jìn)行靈敏度分析,提出了一種新穎的計(jì)算框架,用于在有限變形的時(shí)變載荷條件下設(shè)計(jì)最佳耗散(阻尼)超材料.文獻(xiàn)[26]基于有限變形理論,提出了一種針對(duì)動(dòng)態(tài)問(wèn)題的拓?fù)鋬?yōu)化方法,推導(dǎo)出了一種能在任意變形動(dòng)載荷下減小變形的結(jié)構(gòu).并利用先離散后微分的伴隨變量法建立了一個(gè)通用設(shè)計(jì)靈敏度方程,以準(zhǔn)確地反映結(jié)構(gòu)對(duì)動(dòng)載荷的響應(yīng).此外,先離散后微分方法也被擴(kuò)展至非黏性阻尼系統(tǒng)[27]和非線性系統(tǒng)[28].

上述方法的動(dòng)力響應(yīng)分析均基于Newmark 方法開(kāi)展,因此伴隨變量也是通過(guò)Newmark 積分方案進(jìn)行求解.當(dāng)時(shí)間步長(zhǎng)較大時(shí),Newmark 積分法可能不能精確反映某些高頻振動(dòng)的分量[29],進(jìn)而影響求解精度;減小時(shí)間步長(zhǎng)可提高計(jì)算精度,但會(huì)大幅增加計(jì)算時(shí)間,這也導(dǎo)致了基于Newmark 積分方案的伴隨變量法在一定程度上限制了靈敏度算法性能.改進(jìn)精細(xì)積分法(MPIM)[30]由于采用2N類算法,且載荷項(xiàng)用高斯-勒讓德積分近似,使其在求解結(jié)構(gòu)系統(tǒng)時(shí)域響應(yīng)時(shí)能夠同時(shí)保證高精度和高效率.Ding等[31]針對(duì)非黏性阻尼系統(tǒng)時(shí)域響應(yīng)靈敏度計(jì)算問(wèn)題,提出了基于MPIM 方法的先離散后微分和先微分后離散伴隨變量法,并且對(duì)兩種方法的計(jì)算性能進(jìn)行了詳細(xì)的比較.但非黏性阻尼系統(tǒng)時(shí)域響應(yīng)靈敏度誤差主要來(lái)源于卷積型阻尼模型的近似處理,與黏性阻尼系統(tǒng)存在本質(zhì)上的區(qū)別.因此,其結(jié)論并不能直接適用于黏性阻尼系統(tǒng).

綜上所述,針對(duì)黏性阻尼系統(tǒng)時(shí)域響應(yīng)靈敏度求解問(wèn)題,本文推導(dǎo)了基于MPIM 的先微分后離散和先離散后微分兩種伴隨變量方法的理論表達(dá).通過(guò)數(shù)值算例驗(yàn)證了兩種方法的有效性和準(zhǔn)確性,并與傳統(tǒng)基于Newmark 的方法進(jìn)行比較.討論了微分和離散的先后順序在不同積分方案條件下對(duì)時(shí)域響應(yīng)靈敏度計(jì)算性能的影響,為黏性阻尼系統(tǒng)時(shí)域梯度優(yōu)化算法奠定了基礎(chǔ).

1 問(wèn)題提出

1.1 靈敏度分析的伴隨變量法

對(duì)于N自由度黏性阻尼系統(tǒng),其運(yùn)動(dòng)方程可表示為

對(duì)于時(shí)域響應(yīng)優(yōu)化問(wèn)題,響應(yīng)函數(shù)(或目標(biāo)函數(shù))通常定義為由位移、速度和加速度組成的積分形式,即

在伴隨變量法中,通常引入伴隨變量 λ 構(gòu)造增廣多項(xiàng)式

式中,R(t) 為一恒等于零的多項(xiàng)式,通常與系統(tǒng)運(yùn)動(dòng)方程相關(guān).

1.2 靈敏度分析結(jié)果的準(zhǔn)確性與一致性

準(zhǔn)確性用來(lái)描述算法結(jié)果與解析解的相近程度,兩者的偏差用精度表示.而一致性用來(lái)描述算法結(jié)果與數(shù)值參考解的相近程度,兩者的偏差用一致性誤差表示.

考慮到文章的簡(jiǎn)潔性,本節(jié)以響應(yīng)函數(shù)r為例,分別以ra(t,p)和rb(t,p,Δt) 表示解析解和數(shù)值解.其中p為設(shè)計(jì)變量,t為時(shí)間,Δt為時(shí)間步長(zhǎng).響應(yīng)函數(shù)r對(duì)p靈敏度可定義為

顯然,響應(yīng)函數(shù)的解析解ra(t,p) 與數(shù)值解rb(t,p,Δt)之間存在偏差,這是由數(shù)值解存在的時(shí)間離散誤差而引起的偏差.基于此,由式(5) 和式(6)給出的靈敏度結(jié)果也會(huì)有不同.為區(qū)分靈敏度分析方法的準(zhǔn)確性與一致性,在相同的時(shí)間離散情況下,若某種靈敏度分析方法計(jì)算出的結(jié)果為h(t,p,Δp),則可定義兩種誤差為

按上述理論,ε 定義為計(jì)算精度,其數(shù)值反映靈敏度分析方法的準(zhǔn)確性;而 εc定義為一致性誤差,其數(shù)值反映靈敏度分析方法結(jié)果與數(shù)值參考解偏差.由于時(shí)域逐步積分方法所得動(dòng)力響應(yīng)解僅能保證在離散時(shí)間點(diǎn)處嚴(yán)格滿足動(dòng)力平衡方程,因此先微分后離散伴隨變量法在理論上會(huì)給出不一致的靈敏度求解結(jié)果,較大的一致性誤差,將無(wú)法準(zhǔn)確反映參數(shù)變化對(duì)響應(yīng)函數(shù)的影響,進(jìn)而造成后續(xù)計(jì)算和優(yōu)化結(jié)果的嚴(yán)重偏差.

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

本節(jié)將簡(jiǎn)要介紹后續(xù)內(nèi)容中涉及到的關(guān)于式(1)所表示的黏性阻尼系統(tǒng)動(dòng)力方程的求解方法,即Newmark 方法和MPIM 方法.

2.1 求解動(dòng)力響應(yīng)的Newmark 逐步積分法

利用Newmark 法進(jìn)行求解,其逐步積分列式為[32]

2.2 基于狀態(tài)空間的MPIM法

對(duì)于黏性阻尼系統(tǒng),首先將其運(yùn)動(dòng)方程轉(zhuǎn)化為狀態(tài)空間形式.引入狀態(tài)空間向量

于是,式(1)轉(zhuǎn)化為狀態(tài)空間表達(dá)式

在式(14)中,A為狀態(tài)空間矩陣,z(t) 和Q(t) 分別為狀態(tài)空間中的狀態(tài)向量和力向量.基于狀態(tài)空間表達(dá)式,求解時(shí)域響應(yīng)的MPIM 逐步積分遞推公式為[22]

其中,l為高斯求積節(jié)點(diǎn)數(shù),ωi和 ζi分別為高斯節(jié)點(diǎn)和高斯求積系數(shù),T(Δt) 為指數(shù)矩陣函數(shù),定義為

若已知運(yùn)動(dòng)的初始值z(mì)0,由式(17)迭代計(jì)算可得各時(shí)刻點(diǎn)位移、速度、加速度.當(dāng)時(shí)間步長(zhǎng)較大時(shí),MPIM 方法的計(jì)算精度高于Newmark 方法,而當(dāng)時(shí)間步長(zhǎng)較小時(shí),兩方法雖可得到相近的結(jié)果,但前者的計(jì)算效率更高.接下來(lái),將基于MPIM 方法,分別推導(dǎo)出針對(duì)黏性阻尼系統(tǒng)時(shí)域響應(yīng)靈敏度求解的先微分后離散和先離散后微分伴隨變量法.

3 基于MPIM 的先微分后離散法

式(2)表示的響應(yīng)函數(shù)在狀態(tài)空間內(nèi)可等效的表示為

3.1 響應(yīng)函數(shù)的微分

引入任意2N維伴隨變量向量λT(t),定義

由于式(14)在任意時(shí)刻均成立,因此響應(yīng)函數(shù)r等于恒成立.對(duì)于它們的靈敏度,同樣有

令式(20)對(duì)設(shè)計(jì)變量p求微分,得

由于式(22)對(duì)任意伴隨變量 λT均成立,同時(shí)為了消除上式中含時(shí)域響應(yīng)偏導(dǎo)數(shù)的各項(xiàng),可得如下關(guān)系

對(duì)式(23)進(jìn)行分部積分,并整理得

顯然,式(24)為一階常微分方程的終值問(wèn)題.引入變量代換s=T-t,式(24)即等價(jià)的變?yōu)槌Ｒ?jiàn)的一階常微分方程初值問(wèn)題,即

3.2 響應(yīng)函數(shù)的離散

結(jié)合式(21)、式(22)以及式(25),響應(yīng)函數(shù)r的靈敏度表達(dá)式即為

式中,動(dòng)力響應(yīng)z的在每個(gè)離散時(shí)間點(diǎn)已知,于是響應(yīng)函數(shù)的靈敏度可用梯形積分公式近似表示為

可以看到,為獲取響應(yīng)函數(shù)的時(shí)域靈敏度,本方法先對(duì)其進(jìn)行微分運(yùn)算以推導(dǎo)出伴隨方程,然后在時(shí)域內(nèi)對(duì)伴隨方程進(jìn)行離散求解.因此該方法稱為先微分后離散的伴隨變量法.

4 基于MPIM 的先離散后微分法

4.1 響應(yīng)函數(shù)的離散

對(duì)于先離散后微分的伴隨變量法,首先將式(19) 表示的狀態(tài)空間下的響應(yīng)函數(shù)在各時(shí)間點(diǎn)nΔt(0 ≤n≤T/Δt)進(jìn)行離散

同樣地,響應(yīng)函數(shù)r的靈敏度可以由梯形積分公式近似求解,即

利用式(17)可定義殘值多項(xiàng)式為

顯然,式(30)在任何時(shí)刻tn嚴(yán)格滿足Rn=0 .

4.2 響應(yīng)函數(shù)的微分

引入任意2N維伴隨變量向量,利用殘值多項(xiàng)式(30)構(gòu)造增廣方程

由于Ri恒等于零,因此=zn在任意時(shí)刻點(diǎn)也恒成立,對(duì)于靈敏度同樣有

進(jìn)一步有

將式(33)進(jìn)一步展開(kāi)并整理

由式(30)可知,殘值方程Ri對(duì)狀態(tài)空間向量zi的靈敏度可表示為

將式(36)代入式(35),可得伴隨變量表達(dá)式為

至此,按照逆向的順序可逐步求得滿足式(34)的所有伴隨變量.當(dāng)求取每一時(shí)刻點(diǎn)的所有伴隨變量后,代入式(34)可消去所有含 ?zi/?p的隱式導(dǎo)數(shù)項(xiàng),狀態(tài)空間向量zi對(duì)設(shè)計(jì)變量p的靈敏度即為

基于上述過(guò)程,式(29)中所有導(dǎo)數(shù)項(xiàng)均為已知,因此通過(guò)適當(dāng)?shù)那蠓e準(zhǔn)則可得到響應(yīng)函數(shù)r的近似靈敏度結(jié)果.

本方法先將目標(biāo)函數(shù)在時(shí)域內(nèi)進(jìn)行離散,再對(duì)由離散的狀態(tài)空間向量構(gòu)造的擴(kuò)展方程進(jìn)行微分求出伴隨變量,最終求出時(shí)域響應(yīng)靈敏度,因此將方法稱為先離散后微分的伴隨變量法.

5 數(shù)值算例

本章將通過(guò)兩個(gè)數(shù)值算例驗(yàn)證所提出的黏性阻尼系統(tǒng)時(shí)域響應(yīng)靈敏度求解方法的正確性和有效性,并將其與基于Newmark 方法[23]所得結(jié)果進(jìn)行比較,以討論各方法在計(jì)算精度、效率和一致性問(wèn)題等方面的差異.

5.1 單自由度彈簧與質(zhì)量系統(tǒng)

考慮的單自由度系統(tǒng)模型如圖1 所示.系統(tǒng)質(zhì)量m=1 kg,阻尼c=0.1 N·s/m,剛度k=1 N/s,初始條件 為:u0=1 m,u˙0=0 m/s.Newmark 積分方法的參數(shù):β=0.25,γ=0.5;MPIM 方法的參數(shù):L=4,Ne=15,l=2.

圖1 單自由度黏性阻尼彈簧-質(zhì)量系統(tǒng)Fig.1 Single-DOF spring-mass system with viscous damping model

5.1.1 自由振動(dòng)

在自由振動(dòng)情形下,圖1 所示的單自由度黏性阻尼系統(tǒng)的時(shí)域響應(yīng)及其靈敏度可以解析求解得到.因此,可以利用此情形下推導(dǎo)出的解析解與各方法得到的數(shù)值解進(jìn)行比較.本文中相對(duì)誤差定義為

式中,r′為靈敏度算法計(jì)算結(jié)果,為解析解或數(shù)值參考解.當(dāng)為解析解時(shí),下文簡(jiǎn)稱為相對(duì)誤差,當(dāng)為數(shù)值參考解時(shí),下文簡(jiǎn)稱為一致性相對(duì)誤差.

分別利用基于MPIM 和Newmark 積分方案的先離散后微分和先微分后離散伴隨變量法求解時(shí)域位移響應(yīng)關(guān)于剛度的靈敏度,圖2 為各方法在相同時(shí)間步長(zhǎng)(Δt=0.1 s)下前20 秒的時(shí)域靈敏度相對(duì)誤差.可以看到,四種方法均能得到令人滿意的結(jié)果,且在相同的時(shí)間步長(zhǎng)下,微分和離散的先后順序?qū)PIM 方法的影響顯著,而對(duì)Newmark 方法的影響相對(duì)較小.在所考慮的四種方法中,本文所提出的基于MPIM 的先微分后離散伴隨變量法具有最高的精度.

圖2 相同時(shí)間步長(zhǎng)下(Δ t=0.1 s)各方法的位移響應(yīng)關(guān)于剛度靈敏度的相對(duì)誤差Fig.2 The relative errors of the sensitivities of the displacement w.r.t.stiffness for different methods with the same time step

圖3 和圖4 分別展示了在不同時(shí)間步長(zhǎng)下,各方法求解時(shí)域響應(yīng)靈敏度的相對(duì)誤差.從圖3和圖4 可以看出,隨著時(shí)間步長(zhǎng)減小,各方法的相對(duì)誤差均逐漸降低.當(dāng)時(shí)間步長(zhǎng)相同時(shí),基于MPIM 的先微分后離散方法精度高于對(duì)應(yīng)的Newmark 方法,而基于MPIM 的先離散后微分方法精度則低于相應(yīng)的Newmark 方法.在各不同時(shí)間步長(zhǎng)下,本文所提出的基于MPIM 的先微分后離散伴隨變量法仍具有最高的精度.

圖3 不同時(shí)間步長(zhǎng)下先微分后離散法時(shí)域靈敏度的相對(duì)誤差Fig.3 The relative errors of the transient sensitivities for differentiatethen-discretize method with different time steps

圖4 不同時(shí)間步長(zhǎng)下先離散后微分法時(shí)域靈敏度的相對(duì)誤差Fig.4 The relative errors of the transient sensitivities for discretizethen-differentiate method with different time steps

接下來(lái),對(duì)比研究各方法的計(jì)算效率.表1 展示了在不同時(shí)間步長(zhǎng)下四種方法所花費(fèi)的計(jì)算時(shí)間.可以看到,隨著時(shí)間步長(zhǎng)減小,離散步數(shù)增多,各方法所需計(jì)算時(shí)間也逐漸增多.其中,兩種先微分后離散法效率更高,Newmark 先微分后離散法計(jì)算效率最高.

對(duì)于單自由度自由振動(dòng),可以推導(dǎo)出上述兩個(gè)響應(yīng)函數(shù)的靈敏度解析解.又由于上式中的各項(xiàng)均已求出,因此響應(yīng)函數(shù)的靈敏度也可由梯形法則近似求解.在不同的步數(shù)下,利用各方法求解兩個(gè)響應(yīng)函數(shù)靈敏度的相對(duì)誤差分別如圖5 和圖6 所示.首先可以看到,隨著時(shí)間步長(zhǎng)減小(步數(shù)增加),各方法計(jì)算得到的響應(yīng)函數(shù)靈敏度均降低,且具有二階收斂速度.其次,對(duì)于計(jì)算不同響應(yīng)函數(shù)的靈敏度,基于MPIM 的先微分后離散伴隨變量法相對(duì)誤差均小于相應(yīng)的先離散后微分方法,而基于Newmark 方法的伴隨變量法則對(duì)不同的響應(yīng)函數(shù)呈現(xiàn)出相反的結(jié)果.可以看出,時(shí)間步長(zhǎng)、積分方案和微分與離散的先后順序共同影響響應(yīng)函數(shù)靈敏度相對(duì)誤差.綜合考慮各方法對(duì)兩個(gè)響應(yīng)函數(shù)靈敏度的相對(duì)誤差,所提出的基于MPIM 的先微分后離散伴隨變量法更優(yōu).

圖5 響應(yīng)函數(shù) r1 關(guān)于剛度k 的時(shí)域靈敏度相對(duì)誤差隨離散步數(shù)變化Fig.5 The relative errors of sensitivity result for response functionr1 versus the number of time steps

圖6 響應(yīng)函數(shù) r2 關(guān)于剛度k 的時(shí)域靈敏度相對(duì)誤差隨離散步數(shù)變化Fig.6 The relative errors of sensitivity result for response functionr2 versus the number of time steps

5.1.2 受迫振動(dòng)

接著考慮受迫振動(dòng)情形.假設(shè)對(duì)質(zhì)量塊m施加外激勵(lì)f(t)=sint.對(duì)于單自由度受迫振動(dòng)和多自由度系統(tǒng),由于不易求得其時(shí)域響應(yīng)靈敏度的解析解,因此下文均利用有限差分法(FDM)計(jì)算的時(shí)域響應(yīng)靈敏度結(jié)果作為參考解,攝動(dòng)量取0.01%.

分別利用4 種方法計(jì)算位移響應(yīng)關(guān)于剛度和速度響應(yīng)關(guān)于質(zhì)量的靈敏度,其前20 s 的結(jié)果分別如圖7 和圖8 中小圖所示.從圖7 可以看出,在相同時(shí)間步長(zhǎng)(Δt=0.1 s)和同樣的離散與微分順序下,基于MPIM 伴隨變量法的一致性相對(duì)誤差總體上明顯小于對(duì)應(yīng)的基于Newmark 方法,這一現(xiàn)象隨著振動(dòng)的持續(xù)更加明顯.從圖8 中,也可以得到相似的結(jié)論.綜合考慮圖7 和圖8,從計(jì)算精度的角度看,基于MPIM 的先微分后離散方法的一致性相對(duì)誤差更低.

圖7 相同時(shí)間步長(zhǎng)下(Δ t=0.1 s)各方法位移響應(yīng)關(guān)于剛度du/dk的靈敏度一致性相對(duì)誤差Fig.7 The relative consistency errors of the sensitivities of the displacement w.r.t.stiffness d u/dk for different methods with the same time step (Δ t=0.1 s)

同時(shí),對(duì)比受迫振動(dòng)情形下各方法的計(jì)算效率.表2 展示了在不同時(shí)間步長(zhǎng)下四種方法所花費(fèi)的計(jì)算時(shí)間.結(jié)論同自由振動(dòng)情形一致,隨著時(shí)間步長(zhǎng)減小,離散步數(shù)增多,各方法所需計(jì)算時(shí)間也逐漸增多.其中,兩種先微分后離散法效率更高,Newmark 先微分后離散法計(jì)算效率最高.

表2 受迫振動(dòng)情形下各方法計(jì)算時(shí)間隨步長(zhǎng)變化對(duì)比Table 2 Comparisons of computational time for each method of forced vibration with different time steps

與自由振動(dòng)情形一樣,為提高計(jì)算精度,響應(yīng)函數(shù)靈敏度結(jié)果使用梯形近似積分求解.在不同步數(shù)下,利用各方法求解上述響應(yīng)函數(shù)靈敏度的一致性相對(duì)誤差結(jié)果分別如圖9 和圖10 所示.可以看到,當(dāng)步數(shù)相同時(shí),基于不同積分方案的伴隨變量法一致性誤差明顯不同,且改變MPIM 伴隨變量法微分與離散的先后次序后計(jì)算結(jié)果十分相近,因此對(duì)于本例而言理論上存在的一致性問(wèn)題并未對(duì)MPIM 先微分后離散法的計(jì)算精度產(chǎn)生明顯影響.相反,在求解靈敏度時(shí),Newmark 先微分后離散法相對(duì)誤差顯然低于Newmark 先離散后微分法,這說(shuō)明Newmark 先微分后離散法理論上存在的一致性問(wèn)題實(shí)際上提高了該方法的計(jì)算精度.由此可見(jiàn)積分方案同樣是影響靈敏度分析結(jié)果一致性誤差的重要因素.還可以觀察到,在計(jì)算兩類受迫振動(dòng)響應(yīng)函數(shù)靈敏度時(shí),四種方法計(jì)算結(jié)果的一致性誤差最終趨于相近,這說(shuō)明當(dāng)動(dòng)力響應(yīng)計(jì)算精度足夠高時(shí),一致性問(wèn)題不會(huì)影響方法的可靠性.但需要指出的是,當(dāng)時(shí)間步長(zhǎng)相同時(shí),不同積分方案所需要的計(jì)算時(shí)間相差巨大.將通過(guò)下面的軸向桿振動(dòng)問(wèn)題具體研究各方法的計(jì)算效率.

圖9 響應(yīng)函數(shù) r1 關(guān)于剛度k 的時(shí)域靈敏度一致性相對(duì)誤差隨離散步數(shù)變化Fig.9 The relative consistency errors of sensitivity result for response function r1 versus the number of time steps

圖10 響應(yīng)函數(shù) r3 關(guān)于剛度k 的時(shí)域靈敏度一致性相對(duì)誤差隨離散步數(shù)變化Fig.10 The relative consistency errors of sensitivity result for response function r3 versus the number of time steps

5.2 非比例黏性阻尼軸向桿問(wèn)題

考慮一個(gè)非比例黏性阻尼軸向桿振動(dòng)問(wèn)題.該軸向振動(dòng)桿簡(jiǎn)化模型如圖11 所示.桿的單元質(zhì)量矩陣和單元?jiǎng)偠染仃嚪謩e為

圖11 非比例黏性阻尼軸向桿振動(dòng)示意圖Fig.11 A fixed-free axially vibrating rod with non-proportionally viscous damping model

式中,ρ 表示密度,E表示材料楊氏模量,A為梁的橫截面積,le=L/N為單位桿單元的長(zhǎng)度,N為單元數(shù),各參數(shù)分別A=6.25×10-4m2,E=2.1×104N/m2,ρ=7.8×103kg/m3.桿的全長(zhǎng)為L(zhǎng)=4 m .據(jù)此,桿的整體質(zhì)量矩陣M和整體剛度矩陣K可組裝得到.

該軸向振動(dòng)桿由兩種不同的阻尼材料組成,假設(shè)其阻尼模型由兩種不同的比例阻尼模型描述:C1=αM+β1K和C2=αM+β2K,其中比例因子分別為 α=0.1,β1=0.001,β1=0.002 .因此,該軸向振動(dòng)桿整體可看作非比例黏性阻尼系統(tǒng).

首先,將上例中長(zhǎng)為L(zhǎng)的桿離散為80 個(gè)單元,自由振動(dòng)初始條件為u0=0,u˙0=[1,0,···,0]T.分別利用基于MPIM 和Newmark 的先離散后微分和先微分后離散方法計(jì)算自由端位移響應(yīng)對(duì)彈性模量E的靈敏度,其中離散時(shí)間步長(zhǎng)均為3 ms,持續(xù)時(shí)間為前0.36 s.本例中靈敏度的參考值采用基于模態(tài)疊加法的FDM 求出(設(shè)計(jì)變量的攝動(dòng)量為0.01%).圖12 為四種方法時(shí)域靈敏度結(jié)果相對(duì)于參考解的一致性相對(duì)誤差.可以看到MPIM 伴隨變量法計(jì)算精度明顯高于Newmark 伴隨變量法.同時(shí)可以看出,盡管MPIM 先微分后離散法在理論上存在一致性問(wèn)題,但在適當(dāng)小的時(shí)間步長(zhǎng)下,該方法計(jì)算結(jié)果一致性相對(duì)誤差都能控制在相對(duì)較小的范圍(E≈10-4),因此一致性問(wèn)題并不影響MPIM 先微分后離散法的準(zhǔn)確性和使用效果.

圖12 各方法的位移對(duì)彈性模量 d u/dE 的靈敏度一致性對(duì)誤差Fig.12 The relative consistency errors of the sensitivities of the displacement w.r.t.elastic modulus d u/dE for different methods

為了研究各方法的計(jì)算效率,可將軸向振動(dòng)桿按照要求劃分為不同單元數(shù)目.在時(shí)間步長(zhǎng)為3 ms時(shí),四種時(shí)域響應(yīng)靈敏度方法的計(jì)算時(shí)間隨自由度數(shù)目的變化如圖13 所示.可以看到,先微分后離散方法的計(jì)算時(shí)間要比先離散后微分少得多.這是由于先離散后微分方法需要計(jì)算的伴隨變量數(shù)目遠(yuǎn)多于相同步長(zhǎng)條件下的先微分后離散方法.對(duì)于本文所考慮黏性阻尼系統(tǒng),基于Newmark的先微分后離散方法具有最高的效率,而基于Newmark 的先離散后微分方法計(jì)算效率最低.而對(duì)于非黏性阻尼系統(tǒng),基于MPIM 積分方案的時(shí)域響應(yīng)靈敏度計(jì)算效率全面高于相同條件下基于Newmark 積分方案的方法.這是因?yàn)榉丘ば宰枘嵯到y(tǒng)時(shí)域響應(yīng)靈敏度誤差主要來(lái)源于卷積型阻尼模型的近似處理,而黏性阻尼系統(tǒng)不含卷積型阻尼模型.雖然基于MPIM 的先微分后離散方法的計(jì)算時(shí)間略高于基于Newmark 的先微分后離散方法,但綜合考慮計(jì)算精度,其仍為更推薦的時(shí)域靈敏度求解方法.

圖13 各方法計(jì)算時(shí)間隨自由度增加的變化Fig.13 Calculation time of different DOF for the compared methods of the axially vibrating rod

表3 各方法在不同時(shí)間步長(zhǎng)下響應(yīng)函數(shù)及其靈敏度計(jì)算精度及時(shí)間比較Table 3 Comparisons of the relative consistency errors for response function r1 and its sensitivity w.r.t.stiffness with different time steps using various methods

6 結(jié)論

(1)基于MPIM 的先離散后微分和先微分后離散方法均能得到準(zhǔn)確可靠的時(shí)域響應(yīng)靈敏度結(jié)果.因此本文推導(dǎo)出的兩種方法可有效的求解黏性阻尼系統(tǒng)時(shí)域響應(yīng)靈敏度.

(2)對(duì)于黏性阻尼系統(tǒng)時(shí)域響應(yīng)靈敏度求解問(wèn)題,在相同時(shí)間步長(zhǎng)和離散與微分順序下,基于MPIM 積分方案的計(jì)算精度高于基于Newmark 積分方案的結(jié)果.該差異由各積分方案自身求解精度所致.對(duì)于本文所給出的基于MPIM 的時(shí)域響應(yīng)靈敏度求解方法,先微分后離散方法的精度高于先離散后微分方法.

(3)對(duì)于黏性阻尼系統(tǒng)時(shí)域響應(yīng)靈敏度求解,基于MPIM 和Newmark 方法的先微分后離散方法均存在一致性誤差,改變微分和離散的先后次序?qū)煞椒ㄒ恢滦哉`差的作用不完全一致.通過(guò)減小時(shí)間步長(zhǎng)也可降低算法的一致性誤差.當(dāng)動(dòng)力響應(yīng)計(jì)算精度足夠高時(shí),一致性問(wèn)題并不會(huì)影響時(shí)域靈敏度求解方法使用的可靠性.

(4)雖然MPIM 先微分后離散法計(jì)算效率略低于Newmark 先微分后離散法,但在相同條件下,前者對(duì)時(shí)域響應(yīng)函數(shù)及靈敏度的計(jì)算精度遠(yuǎn)高于后者.綜合考慮精度、效率和一致性問(wèn)題,本文給出的基于MPIM 的先微分后離散伴隨變量法表現(xiàn)更優(yōu),最適合應(yīng)用于黏性阻尼系統(tǒng)時(shí)域梯度優(yōu)化算法.