黏性阻尼系統時域響應靈敏度及其一致性研究1)

張磊 張 嚴 丁 喆

(武漢科技大學冶金裝備及其控制教育部重點實驗室,武漢 430081)

(武漢科技大學機械傳動與制造工程湖北省重點實驗室,武漢 430081)

引言

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

動力響應的靈敏度求解對結構動態性能的分析及優化非常重要.目前為止,眾多學者已經對該領域進行了廣泛的研究.Choi等[11]在波導濾波器的優化設計中,利用增廣拉格朗日法、材料導數概念和伴隨變量法,系統的推導了頻域解析靈敏度公式.周大為等[12]針對黏性阻尼系統的頻響振幅靈敏度求解問題,并結合廣義模態截斷擴增方法,提出了一種具有更高精度的伴隨變量法.Zhao等[13]提出了一種自適應混合展開法來確定需要計算的低階特征頻率和特征模數,然后再采用伴隨變量法進行靈敏度分析,以克服簡諧激勵下的拓撲優化問題計算量大的困難.

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

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

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

綜上所述,針對黏性阻尼系統時域響應靈敏度求解問題,本文推導了基于MPIM 的先微分后離散和先離散后微分兩種伴隨變量方法的理論表達.通過數值算例驗證了兩種方法的有效性和準確性,并與傳統基于Newmark 的方法進行比較.討論了微分和離散的先后順序在不同積分方案條件下對時域響應靈敏度計算性能的影響,為黏性阻尼系統時域梯度優化算法奠定了基礎.

1 問題提出

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

對于N自由度黏性阻尼系統,其運動方程可表示為

對于時域響應優化問題,響應函數(或目標函數)通常定義為由位移、速度和加速度組成的積分形式,即

在伴隨變量法中,通常引入伴隨變量 λ 構造增廣多項式

式中,R(t) 為一恒等于零的多項式,通常與系統運動方程相關.

1.2 靈敏度分析結果的準確性與一致性

準確性用來描述算法結果與解析解的相近程度,兩者的偏差用精度表示.而一致性用來描述算法結果與數值參考解的相近程度,兩者的偏差用一致性誤差表示.

考慮到文章的簡潔性,本節以響應函數r為例,分別以ra(t,p)和rb(t,p,Δt) 表示解析解和數值解.其中p為設計變量,t為時間,Δt為時間步長.響應函數r對p靈敏度可定義為

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

按上述理論,ε 定義為計算精度,其數值反映靈敏度分析方法的準確性;而 εc定義為一致性誤差,其數值反映靈敏度分析方法結果與數值參考解偏差.由于時域逐步積分方法所得動力響應解僅能保證在離散時間點處嚴格滿足動力平衡方程,因此先微分后離散伴隨變量法在理論上會給出不一致的靈敏度求解結果,較大的一致性誤差,將無法準確反映參數變化對響應函數的影響,進而造成后續計算和優化結果的嚴重偏差.

2 動力方程的數值解法

本節將簡要介紹后續內容中涉及到的關于式(1)所表示的黏性阻尼系統動力方程的求解方法,即Newmark 方法和MPIM 方法.

2.1 求解動力響應的Newmark 逐步積分法

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

2.2 基于狀態空間的MPIM法

對于黏性阻尼系統,首先將其運動方程轉化為狀態空間形式.引入狀態空間向量

于是,式(1)轉化為狀態空間表達式

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

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

若已知運動的初始值z0,由式(17)迭代計算可得各時刻點位移、速度、加速度.當時間步長較大時,MPIM 方法的計算精度高于Newmark 方法,而當時間步長較小時,兩方法雖可得到相近的結果,但前者的計算效率更高.接下來,將基于MPIM 方法,分別推導出針對黏性阻尼系統時域響應靈敏度求解的先微分后離散和先離散后微分伴隨變量法.

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

式(2)表示的響應函數在狀態空間內可等效的表示為

3.1 響應函數的微分

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

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

令式(20)對設計變量p求微分,得

由于式(22)對任意伴隨變量 λT均成立,同時為了消除上式中含時域響應偏導數的各項,可得如下關系

對式(23)進行分部積分,并整理得

顯然,式(24)為一階常微分方程的終值問題.引入變量代換s=T-t,式(24)即等價的變為常見的一階常微分方程初值問題,即

3.2 響應函數的離散

結合式(21)、式(22)以及式(25),響應函數r的靈敏度表達式即為

式中,動力響應z的在每個離散時間點已知,于是響應函數的靈敏度可用梯形積分公式近似表示為

可以看到,為獲取響應函數的時域靈敏度,本方法先對其進行微分運算以推導出伴隨方程,然后在時域內對伴隨方程進行離散求解.因此該方法稱為先微分后離散的伴隨變量法.

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

4.1 響應函數的離散

對于先離散后微分的伴隨變量法,首先將式(19) 表示的狀態空間下的響應函數在各時間點nΔt(0 ≤n≤T/Δt)進行離散

同樣地,響應函數r的靈敏度可以由梯形積分公式近似求解,即

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

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

4.2 響應函數的微分

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

由于Ri恒等于零,因此=zn在任意時刻點也恒成立,對于靈敏度同樣有

進一步有

將式(33)進一步展開并整理

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

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

至此,按照逆向的順序可逐步求得滿足式(34)的所有伴隨變量.當求取每一時刻點的所有伴隨變量后,代入式(34)可消去所有含 ?zi/?p的隱式導數項,狀態空間向量zi對設計變量p的靈敏度即為

基于上述過程,式(29)中所有導數項均為已知,因此通過適當的求積準則可得到響應函數r的近似靈敏度結果.

本方法先將目標函數在時域內進行離散,再對由離散的狀態空間向量構造的擴展方程進行微分求出伴隨變量,最終求出時域響應靈敏度,因此將方法稱為先離散后微分的伴隨變量法.

5 數值算例

本章將通過兩個數值算例驗證所提出的黏性阻尼系統時域響應靈敏度求解方法的正確性和有效性,并將其與基于Newmark 方法[23]所得結果進行比較,以討論各方法在計算精度、效率和一致性問題等方面的差異.

5.1 單自由度彈簧與質量系統

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

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

5.1.1 自由振動

在自由振動情形下,圖1 所示的單自由度黏性阻尼系統的時域響應及其靈敏度可以解析求解得到.因此,可以利用此情形下推導出的解析解與各方法得到的數值解進行比較.本文中相對誤差定義為

式中,r′為靈敏度算法計算結果,為解析解或數值參考解.當為解析解時,下文簡稱為相對誤差,當為數值參考解時,下文簡稱為一致性相對誤差.

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

圖2 相同時間步長下(Δ t=0.1 s)各方法的位移響應關于剛度靈敏度的相對誤差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 分別展示了在不同時間步長下,各方法求解時域響應靈敏度的相對誤差.從圖3和圖4 可以看出,隨著時間步長減小,各方法的相對誤差均逐漸降低.當時間步長相同時,基于MPIM 的先微分后離散方法精度高于對應的Newmark 方法,而基于MPIM 的先離散后微分方法精度則低于相應的Newmark 方法.在各不同時間步長下,本文所提出的基于MPIM 的先微分后離散伴隨變量法仍具有最高的精度.

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

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

接下來,對比研究各方法的計算效率.表1 展示了在不同時間步長下四種方法所花費的計算時間.可以看到,隨著時間步長減小,離散步數增多,各方法所需計算時間也逐漸增多.其中,兩種先微分后離散法效率更高,Newmark 先微分后離散法計算效率最高.

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

圖5 響應函數 r1 關于剛度k 的時域靈敏度相對誤差隨離散步數變化Fig.5 The relative errors of sensitivity result for response functionr1 versus the number of time steps

圖6 響應函數 r2 關于剛度k 的時域靈敏度相對誤差隨離散步數變化Fig.6 The relative errors of sensitivity result for response functionr2 versus the number of time steps

5.1.2 受迫振動

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

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

圖7 相同時間步長下(Δ t=0.1 s)各方法位移響應關于剛度du/dk的靈敏度一致性相對誤差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)

同時,對比受迫振動情形下各方法的計算效率.表2 展示了在不同時間步長下四種方法所花費的計算時間.結論同自由振動情形一致,隨著時間步長減小,離散步數增多,各方法所需計算時間也逐漸增多.其中,兩種先微分后離散法效率更高,Newmark 先微分后離散法計算效率最高.

表2 受迫振動情形下各方法計算時間隨步長變化對比Table 2 Comparisons of computational time for each method of forced vibration with different time steps

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

圖9 響應函數 r1 關于剛度k 的時域靈敏度一致性相對誤差隨離散步數變化Fig.9 The relative consistency errors of sensitivity result for response function r1 versus the number of time steps

圖10 響應函數 r3 關于剛度k 的時域靈敏度一致性相對誤差隨離散步數變化Fig.10 The relative consistency errors of sensitivity result for response function r3 versus the number of time steps

5.2 非比例黏性阻尼軸向桿問題

考慮一個非比例黏性阻尼軸向桿振動問題.該軸向振動桿簡化模型如圖11 所示.桿的單元質量矩陣和單元剛度矩陣分別為

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

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

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

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

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

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

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

表3 各方法在不同時間步長下響應函數及其靈敏度計算精度及時間比較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 結論

(1)基于MPIM 的先離散后微分和先微分后離散方法均能得到準確可靠的時域響應靈敏度結果.因此本文推導出的兩種方法可有效的求解黏性阻尼系統時域響應靈敏度.

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

(3)對于黏性阻尼系統時域響應靈敏度求解,基于MPIM 和Newmark 方法的先微分后離散方法均存在一致性誤差,改變微分和離散的先后次序對兩方法一致性誤差的作用不完全一致.通過減小時間步長也可降低算法的一致性誤差.當動力響應計算精度足夠高時,一致性問題并不會影響時域靈敏度求解方法使用的可靠性.

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