史俊磊， 張 磊， 丁 喆

(1. 武漢科技大學(xué) 冶金裝備及其控制教育部重點實驗室， 武漢 430081；2. 武漢科技大學(xué) 機械傳動與制造工程湖北省重點實驗室， 武漢 430081)

靈敏度分析可以定量地確定系統(tǒng)參數(shù)改變對系統(tǒng)動態(tài)特性的影響，在結(jié)構(gòu)優(yōu)化、模型修正和參數(shù)識別等領(lǐng)域有著廣泛應(yīng)用[1-2]。黏性阻尼模型由于其模型簡單和計算方便的優(yōu)點已被廣泛研究和應(yīng)用[3]，但隨著各類智能材料和高阻尼復(fù)合材料的發(fā)展，原有的黏性阻尼模型假設(shè)已不能準(zhǔn)確地描述該類結(jié)構(gòu)系統(tǒng)的阻尼特性。卷積型非黏滯阻尼模型由于假設(shè)阻尼力與質(zhì)點速度的時間歷程相關(guān)，其數(shù)學(xué)表達式為阻尼核函數(shù)與質(zhì)點速度的卷積,因而能夠準(zhǔn)確地描述黏彈性阻尼材料的遲滯效應(yīng)。目前，非黏滯阻尼系統(tǒng)的研究主要針對其特征值分析、時域響應(yīng)分析和參數(shù)識別等[4-7]方面，其靈敏度分析也主要集中在特征靈敏度分析和時域響應(yīng)靈敏度分析[8-9]，而非黏滯阻尼系統(tǒng)頻響函數(shù)(frequency response function, FRF)靈敏度的研究相對較為缺乏。

FRF靈敏度分析可以評估系統(tǒng)參數(shù)對FRF的影響程度。一般來說，結(jié)構(gòu)系統(tǒng)靈敏度的計算方法分為三類：有限差分法[10]、直接微分法和伴隨變量法[11]。其中，有限差分法理論簡單，但在某些情況下計算效率較低。當(dāng)設(shè)計變量數(shù)目大于約束數(shù)目時，伴隨變量法具有更高的計算效率，反之則直接微分法的計算效率較高。針對FRF靈敏度分析，Lin[12]提出了一種加權(quán)FRF靈敏度法，該方法即使在模型誤差較大的情況下，仍具有良好的收斂性能。Qu[13]基于冪級數(shù)展開和模態(tài)疊加法提出無阻尼系統(tǒng)FRF的模態(tài)加速法以及FRF靈敏度的雙模態(tài)加速法，減小該系統(tǒng)FRF及其靈敏度的模態(tài)截斷誤差。Wu等[14]基于Sturm定理得到參與FRF靈敏度計算的模態(tài)數(shù)目，并利用激勵頻率收斂冪級數(shù)的部分和近似被截斷高階模態(tài)影響，進而推導(dǎo)出無阻尼系統(tǒng)的FRF靈敏度表達式。Wu等[15]根據(jù)Sherman-Morrison公式，提出了一種僅利用初始FRF和參數(shù)攝動來計算靈敏度的方法，該方法避免了重復(fù)有限元分析，提高了計算效率。Erik等[16]針對黏滯阻尼系統(tǒng)，提出了一種高計算效率和高精度計算傳遞函數(shù)靈敏度的方法。Shadan等[17]利用靈敏度方法，將FRF變化與結(jié)構(gòu)參數(shù)變化關(guān)聯(lián)，提出了一個靈敏度方程來減小FRF數(shù)據(jù)不完全性的不利影響。Xiao等[18]針對非經(jīng)典阻尼系統(tǒng)，在僅使用感興趣頻率范圍內(nèi)模態(tài)的情況下推導(dǎo)出FRF靈敏度的表達式，該方法在不涉及狀態(tài)空間方程的基礎(chǔ)上修正模態(tài)截斷誤差。Zhu等[19]考慮多重擾動對系統(tǒng)動剛度矩陣的影響，采用Sherman-Morrison-Woodbury和有限差分法對系統(tǒng)FRF靈敏度進行了求解。上述FRF靈敏度求解方法均建立在無阻尼或黏性阻尼模型的基礎(chǔ)上，無法直接應(yīng)用于非黏滯阻尼模型系統(tǒng)。因此亟需提出一種適用于非黏滯阻尼系統(tǒng)FRF靈敏度的計算方法。

本文研究具有卷積型非黏滯阻尼模型結(jié)構(gòu)系統(tǒng)FRF靈敏度問題。非黏滯阻尼系統(tǒng)的特征值求解復(fù)雜，尤其對于大型結(jié)構(gòu)系統(tǒng)，特征值分析過程非常耗時。為提高計算效率，通常只利用有限數(shù)目的低階模態(tài)近似求解其動力響應(yīng)和靈敏度。但忽略高階模態(tài)不可避免地引入模態(tài)截斷誤差，對非黏滯阻尼系統(tǒng)動力響應(yīng)和靈敏度的計算精度影響尤為顯著。因此必須提出一種適用于非黏滯阻尼結(jié)構(gòu)系統(tǒng)FRF靈敏度求解的改進方法，使其能夠補償高階截斷模態(tài)的影響，進而提高系統(tǒng)FRF靈敏度計算的精度。解決該問題的難點和關(guān)鍵在于建立非黏滯阻尼系統(tǒng)模態(tài)和系統(tǒng)矩陣間的關(guān)系。為消除非黏滯阻尼系統(tǒng)FRF靈敏度的模態(tài)截斷誤差，本文利用Neumann級數(shù)展開定理建立結(jié)構(gòu)系統(tǒng)模態(tài)與系統(tǒng)矩陣間的關(guān)系，考慮其第一項和前兩項近似，將系統(tǒng)高階模態(tài)對FRF靈敏度的貢獻用系統(tǒng)矩陣和低階模態(tài)進行表達，推導(dǎo)非黏滯阻尼系統(tǒng)FRF靈敏度一階近似法和二階近似方法的表達式。通過數(shù)值算例驗證了所提出方法的有效性和準(zhǔn)確性。結(jié)果表明，與傳統(tǒng)方法相比，本文所提出的改進計算方法兼顧了計算精度和效率，更適用于非黏滯阻尼結(jié)構(gòu)系統(tǒng)FRF靈敏度的求解。

1 理論背景

卷積型非黏滯阻尼模型的阻尼力與質(zhì)點速度的時間歷程相關(guān)，數(shù)學(xué)表達式為質(zhì)點速度與某一阻尼核函數(shù)的卷積

(1)

g(t)=Cnvg(t)

(2)

式中：t為時間；τ為遲滯時間變量；Fd(t)∈N為阻尼力向量；N為速度向量；Cnv∈N×N為非黏滯阻尼系數(shù)矩陣；N為系統(tǒng)的自由度；g(t)為阻尼核函數(shù)。式(1)所表示的卷積型非黏滯阻尼模型適用于線性系統(tǒng)，它不僅考慮了當(dāng)前速度的影響，還考慮了非局部時間的速度影響，更加準(zhǔn)確地描述了黏彈性阻尼材料的遲滯效應(yīng)[20]。

理論上，任何使系統(tǒng)能量耗散函數(shù)非負的表達式都可作為卷積型非黏滯阻尼模型的核函數(shù)。因此阻尼核函數(shù)不唯一，通過選取不同阻尼核函數(shù)可適應(yīng)不同黏彈性材料的阻尼特性，常用的阻尼核函數(shù)模型有：Biot阻尼模型[21]、指數(shù)阻尼模型和Golla-Hughes-McTavish(GHM)阻尼模型等。特別地，當(dāng)選取的核函數(shù)滿足g(t-τ)=Cδ(t-τ)(C∈N×N是黏性阻尼矩陣，δ(t)為狄拉克函數(shù))時，卷積型非黏滯阻尼模型將退化為傳統(tǒng)的黏性阻尼模型。因此，黏性阻尼模型可看作該非黏滯阻尼模型的一種特殊形式。

含卷積型非黏滯阻尼的N自由度系統(tǒng)時域運動方程可表示為

(3)

其初始條件為

(4)

式中：M和K∈N×N分別是質(zhì)量矩陣和剛度矩陣；f(t)∈N是外激勵向量；x(t)和N分別是位移和加速度向量。

對式(3)進行拉普拉斯變換得

D(s)X(s)=F(s)

(5)

(6)

式中：λj和φj分別為第j階特征值和特征向量；m為非黏滯阻尼系統(tǒng)的特征值數(shù)目m由兩部分組成

m=2N+P，P>0

(7)

式中：2N為彈性模態(tài)數(shù)目，彈性模態(tài)是以復(fù)共軛對的形式出現(xiàn)；P為非黏性模態(tài)數(shù)目，非黏性模態(tài)的形式為實數(shù)。由式(7)可知，雖然系統(tǒng)有N自由度，但其特征解的個數(shù)大于2N。這是非黏滯阻尼系統(tǒng)和黏性阻尼系統(tǒng)的一個重要區(qū)別，也是導(dǎo)致非黏滯阻尼系統(tǒng)特征值求解耗時的原因之一。

2 非黏滯阻尼系統(tǒng)FRF靈敏度的求解方法

本章將根據(jù)上述非黏滯阻尼系統(tǒng)的運動方程，推導(dǎo)出傳統(tǒng)求解非黏滯阻尼系統(tǒng)FRF靈敏度的表達式。

2.1 直接頻響法

由式(5)可知

X(s)=D-1(s)F(s)=H(s)F(s)

(8)

其中，H(s)為系統(tǒng)FRF，它可表示為

H(s)=D-1(s)=(s2M+sG(s)+K)-1

(9)

該方法稱為直接頻響法(direct frequency response method, DFRM)，將式(9)關(guān)于設(shè)計變量p求導(dǎo)得到系統(tǒng)FRF靈敏度的表達式

(10)

2.2 模態(tài)疊加法

傳統(tǒng)求解系統(tǒng)FRF靈敏度的方法除了DFRM外還有模態(tài)疊加法(modal superposition method，MSM)。MSM是在系統(tǒng)特征方程求出系統(tǒng)模態(tài)的基礎(chǔ)上，利用留數(shù)定理推導(dǎo)出系統(tǒng)FRF矩陣H(s)。當(dāng)系統(tǒng)全部模態(tài)參與計算時，MSM可以獲得精確解，但其計算效率較低，并且在某些情況下，無法求出系統(tǒng)全部模態(tài)。因此，通常只考慮低階模態(tài)，忽略高階模態(tài)對FRF的影響

(11)

(12)

H(s)H(s)-1=I

(13)

式中，I為單位矩陣。同時對式(13)兩邊關(guān)于設(shè)計變量p求導(dǎo)，并考慮H(s)=D(s)-1得到

(14)

將式(11)代入式(14)得到系統(tǒng)FRF靈敏度表達式

(15)

由式(15)可知，MSM計算出FRF靈敏度的表達式中存在兩個模態(tài)截斷項。因此，模態(tài)截斷對系統(tǒng)FRF靈敏度分析的影響較大。由模態(tài)截斷引起的模態(tài)截斷誤差可以定量地表示為

(16)

當(dāng)L?m時，MSM計算出的FRF靈敏度表達式不僅不準(zhǔn)確，且在某些特殊情況下，計算結(jié)果是錯誤的。為保證系統(tǒng)FRF靈敏度的計算精度，MSM需要保留較高階模態(tài)，但相應(yīng)的計算效率降低。因此，需要提出一種兼顧計算精度和效率的非黏滯阻尼系統(tǒng)FRF靈敏度的計算方法。

3 非黏滯阻尼系統(tǒng)FRF靈敏度的改進計算方法

在本章將提出兩種考慮高階模態(tài)影響的FRF靈敏度的改進計算方法。

3.1 非黏滯阻尼系統(tǒng)模態(tài)與系統(tǒng)矩陣關(guān)系

本小節(jié)介紹一種利用Neumann級數(shù)展開定理建立的非黏滯阻尼系統(tǒng)模態(tài)與系統(tǒng)矩陣間的關(guān)系[22]。

Neumann級數(shù)展開定理的表達式為

(I+A)-1=I-A+A2-A3+A4-…

(17)

將式(11)表示為矩陣和向量相乘的形式，并考慮所有階模態(tài)得

H(s)=-UΘ-1Λ-1(I-sΛ-1)-1UT

(18)

其中I∈m×m，Λ=diag[λ1,λ2,…,λm]，U=[φ1,φ2,…,φm]和Θ=diag[θ1,θ2,…,θm]。利用Neumann級數(shù)展開定理將式(I-sΛ-1)-1展開得

(19)

將式(19)代入式(18)，可以得到

H(s)=-UΘ-1Λ-1UT-sUΘ-1Λ-2UT-

s2UΘ-1Λ-3UT-s3UΘ-1Λ-4UT-…=

(20)

由式(9)可知，F(xiàn)RF還可表示為

H(s)=D-1(s)=(s2M+sG(s)+K)-1=

[I+sK-1(G(s)+sM)]-1K-1

(21)

其中I∈N×N，根據(jù)Neumann級數(shù)展開定理將[I+sK-1(G(s)+sM)]-1展開并代入到式(21)可以得到

H(s)=K-1-sK-1(G(s)+sM)K-1+[sK-1(G(s)+sM)]2K-1-[sK-1(G(s)+sM)]3K-1+…

(22)

將式(22)展開，并按s的次數(shù)由低到高排列可以得到

H(s)=K-1-sK-1G(s)K-1+s2K-1G(s)K-1G(s)K-1-s2K-1MK-1+…=K-1-sK-1G(s)K-1+s2K-1[G(s)K-1G(s)-M]K-1+…=

(23)

其中：

(24)

因此，非黏滯阻尼系統(tǒng)矩陣與模態(tài)間關(guān)系

-UΘ-1Λ-k-1UT=Γk

(25)

由式(25)可知，由于阻尼核函數(shù)G(s)與頻率相關(guān)，當(dāng)k≥2時，所有展開項都會受前項的影響，故其不能再進一步近似等效。因此，考慮非黏滯阻尼系統(tǒng)矩陣與模態(tài)間的前兩項近似，并引出其關(guān)系

-UΘ-1Λ-1UT=K-1

(26)

-UΘ-1Λ-1UT-UΘ-1Λ-2UT=K-1-K-1G0K-1

(27)

3.2 FRF靈敏度改進計算方法

在上述關(guān)系的基礎(chǔ)上，本小節(jié)將推導(dǎo)考慮高階模態(tài)影響的系統(tǒng)FRF靈敏度的表達式。將式(18)分為可用的低階模態(tài)和被截斷的高階模態(tài)兩部分

(28)

式中，H表示被截斷的高階模態(tài)數(shù)目。利用Neumann級數(shù)展開將式(28)右側(cè)的第一項展開并取前兩階近似得

(29)

聯(lián)立式(26)～(29)將被截斷的高階模態(tài)用系統(tǒng)低階模態(tài)和系統(tǒng)矩陣表示為

(30)

由于高階模態(tài)對系統(tǒng)FRF的影響是利用非黏滯阻尼系統(tǒng)模態(tài)與系統(tǒng)矩陣關(guān)系的前兩項近似表達，因此稱該方法為二階近似法(second-order approximation method, SAM)。聯(lián)立式(11)和(30)得到考慮高階模態(tài)影響的系統(tǒng)FRF表達式

(31)

通過合并同類項，式(31)可進一步化簡為

(32)

將式(32)代入式(14)得到利用SAM方法計算的非黏滯阻尼系統(tǒng)FRF靈敏度表達式

(33)

SAM計算得到頻響函數(shù)對應(yīng)的模態(tài)截斷誤差為

(34)

其中，下標(biāo)H-SAM表示利用SAM計算出系統(tǒng)頻響函數(shù)，其頻響函數(shù)靈敏度的模態(tài)截斷誤差為

(35)

其中，下標(biāo)Hp-SAM表示利用SAM計算出系統(tǒng)頻響函數(shù)關(guān)于p的靈敏度，如果用非黏滯阻尼系統(tǒng)模態(tài)與系統(tǒng)矩陣關(guān)系第一項近似FRF高階模態(tài)的貢獻量可以得到

(36)

聯(lián)立式(11)和(34)并化簡可以得到一階近似法(first-order approximation method，F(xiàn)AM)的表達式

(37)

將式(35)代入式(14)得到由FAM計算出的非黏滯阻尼系統(tǒng)FRF靈敏度表達式

(38)

FAM計算得到頻響函數(shù)對應(yīng)的模態(tài)截斷誤差為

(39)

其中，下標(biāo)H-FAM表示利用FAM計算出系統(tǒng)頻響函數(shù)，其頻響函數(shù)靈敏度的模態(tài)截斷誤差為

(40)

其中，下標(biāo)Hp-FAM表示利用FAM計算出系統(tǒng)頻響函數(shù)關(guān)于p的靈敏度，與MSM相比，F(xiàn)AM和SAM方法考慮高階模態(tài)對非黏滯阻尼系統(tǒng)FRF靈敏度的影響，但計算更加復(fù)雜。為了更清晰地表達FAM和SAM的計算流程并便于編程，其算法流程圖如圖1所示。

圖1 FAM和SAM算法流程圖Fig.1 Algorithm flow chart of FAM and SAM

4 例子與結(jié)論

4.1 質(zhì)量-彈簧系統(tǒng)

在本節(jié)，通過如圖2所示的含非黏滯阻尼模型的N自由度系統(tǒng)(N可以根據(jù)需要進行調(diào)整)驗證本文提出的非黏滯阻尼系統(tǒng)FRF靈敏度的FAM和SAM的有效性，并將其與DFRM和MSM進行對比，討論四種方法的計算精度和效率。

本例選取雙指數(shù)阻尼核函數(shù)，含該阻尼模型的N自由度運動方程表達式為

(41)

式中：c為非黏滯阻尼參數(shù)；μ1和μ2均為松弛參數(shù)。

圖2 含非黏滯阻尼器的N自由度系統(tǒng)Fig.2 N DOF non-viscously damped system

4.1.1 計算精度

本小節(jié)選取自由度數(shù)目N=4，系統(tǒng)參數(shù)：m1=100 kg，m2=200 kg，m3=300 kg,m4=400 kg和ke=50 000 N/m。非黏滯阻尼模型的參數(shù)為：c=50 Ns/m，μ1=50 s-1和μ2=70 s-1。

通過求解狀態(tài)空間方程的特征方程可以得到系統(tǒng)的特征值，該系統(tǒng)的特征值如表1所示。可以注意到，本例中的特征值個數(shù)為16，其中彈性模態(tài)的個數(shù)為8，非黏性模態(tài)的個數(shù)P=r×q=4×2=8(其中r為阻尼系數(shù)矩陣的秩，q為頻域內(nèi)阻尼核函數(shù)分母的最高次冪)。由表1可知，彈性模態(tài)的形式為復(fù)數(shù)，其實部表示系統(tǒng)的衰減系數(shù)，虛部表示系統(tǒng)的阻尼振動頻率。由于彈性模態(tài)的實部與非黏性模態(tài)均是負數(shù)，所以該系統(tǒng)是穩(wěn)定的。

表1 四自由度非黏滯阻尼系統(tǒng)的特征解

為了驗證本文提出方法的有效性和準(zhǔn)確性，首先選取非黏滯阻尼參數(shù)c為設(shè)計變量。分別利用DFRM，MSM，F(xiàn)AM和SAM計算該四自由度非黏滯阻尼系統(tǒng)在0～40 rad/s范圍內(nèi)FRF靈敏度。其中，將DFRM計算得到的FRF靈敏度的結(jié)果視為參考值。為了方便起見，圖中只展示系統(tǒng)FRF靈敏度矩陣中第一對角元素H11。H11表示輸出和輸入均為1號自由度時得到的系統(tǒng)FRF。本節(jié)中各方法所使用的頻率步長均為0.2 rad/s。

圖3為三種方法選取全部彈性模態(tài)計算的FRF靈敏度相對誤差的結(jié)果。其中，DFRM計算的實際頻響函數(shù)靈敏度也展示在圖3中。本文中，頻響函數(shù)靈敏度相對誤差ER的計算公式為

(42)

從圖3可以看出，MSM，F(xiàn)AM和SAM計算出的FRF靈敏度在選取頻率范圍內(nèi)誤差均在允許范圍內(nèi)。FAM和SAM的振幅相對誤差明顯小于MSM的振幅相對誤差。SAM的相位相對誤差明顯小于FAM的相位相對誤差。因此，在所比較的三種FRF靈敏度計算方法中，SAM計算精度最高。

圖3 全部彈性模態(tài)近似頻響函數(shù)靈敏度的相對誤差Fig.3 The relative error of FRF sensitivity matrixapproximated in terms of all elastic modes.

為分析模態(tài)截斷誤差對三種方法計算精度的影響，選取第一對(前兩階)彈性模態(tài)參與計算FRF的靈敏度，這就意味著系統(tǒng)的模態(tài)從第三階開始被截斷，其結(jié)果如圖4所示?？梢钥闯觯?dāng)頻率小于10 rad/s時，MSM的計算結(jié)果與參考值仍存在較明顯的誤差，而FAM和SAM的計算結(jié)果相較更接近參考值。當(dāng)頻率大于10 rad/s時，三種方法的FRF靈敏度結(jié)果均與參考值存在顯著差異。

圖4 第一對彈性模態(tài)近似頻響函數(shù)靈敏度Fig.4 The sensitivity of FRF matrix approximated in terms ofthe first pair elastic modes.

圖5是三種方法選取前兩對彈性模態(tài)計算出的FRF靈敏度的結(jié)果。與圖4結(jié)果類似，F(xiàn)AM和SAM的計算精度仍高于相同情形下MSM得到的結(jié)果。由于選取的最高階模態(tài)頻率為16.71 rad/s，因此FAM和SAM的計算結(jié)果在該頻率以內(nèi)的計算結(jié)果與參考值一致，而MSM計算的靈敏度存在較為明顯的誤差。

圖5 前兩對彈性模態(tài)近似頻響函數(shù)靈敏度Fig.5 The sensitivity of FRF matrix approximated in terms ofthe first two pair elastic modes.

進一步增加參與計算FRF靈敏度的模態(tài)數(shù)目，考慮前三對彈性模態(tài)的各方法計算結(jié)果如圖6所示?？梢钥闯?，相較于僅考慮第一對和前兩對的靈敏度計算結(jié)果，三種方法的FRF靈敏度計算結(jié)果均有所提高。通過對比圖4～圖6可以看出，當(dāng)利用本文提出的FAM和SAM計算非黏滯阻尼系統(tǒng)FRF靈敏度時，若希望得到一定頻率范圍內(nèi)的準(zhǔn)確結(jié)果，則需要該頻率范圍內(nèi)的所有彈性模態(tài)參與計算。

圖6 前三對彈性模態(tài)近似頻響函數(shù)靈敏度Fig.6 The sensitivity of FRF matrix approximated in terms ofthe first three pair elastic modes

為了更加直觀地展示三種方法的準(zhǔn)確性，圖7為三種方法選取前三對彈性模態(tài)計算出FRF靈敏度的相對誤差?？梢钥闯?，各方法的幅值相對誤差高于相應(yīng)的相位誤差。當(dāng)選取相同階模態(tài)時，F(xiàn)AM和SAM的FRF靈敏度計算結(jié)果幾乎相同，但均明顯高于MSM的結(jié)果。圖8給出SAM選取不同階模態(tài)計算出FRF靈敏度相對誤差?？梢钥闯觯黾訁⑴c計算的模態(tài)數(shù)目可大幅降低其FRF靈敏度的相對誤差。結(jié)果表明，相較于傳統(tǒng)的MSM方法，由于本文所提出的FAM和SAM方法考慮了非黏滯阻尼系統(tǒng)高階模態(tài)截斷誤差的影響，因此能夠大幅提高非黏滯阻尼系統(tǒng)FRF靈敏度的求解精度。

圖7 選取前三對彈性模態(tài)計算出頻響函數(shù)靈敏度的相對誤差

圖8 SAM使用不同階彈性模態(tài)計算FRF靈敏度相對誤差

4.1.2 計算效率

為了比較DFRM，MSM，F(xiàn)AM和SAM的計算效率， 本小節(jié)通過圖2所示的N自由度非黏滯阻尼系統(tǒng)，比較四種方法求解該系統(tǒng)頻響函數(shù)H11靈敏度所消耗的時間。假定頻率范圍為0～5 rad/s，頻率步長為0.1 rad/s，選取系統(tǒng)自由度N的變化范圍0～1 500。

圖9表示使用不同計算方法的FRF靈敏度計算時間隨自由度增加的變化曲線?？梢钥闯觯?dāng)自由度小于600時，四種方法的計算時間相差不大。但隨著自由度數(shù)目繼續(xù)增加，DFRM的計算時間增加迅速，而FAM和SAM的計算時間增長相對平穩(wěn)。這是因為在每個外激勵頻率下DFRM要對系統(tǒng)動剛度矩陣求逆，而系統(tǒng)自由度增加必然會引起系統(tǒng)動剛度矩陣維數(shù)的增加，系統(tǒng)動剛度矩陣求逆消耗的時間就越長。而FAM和SAM避免了矩陣求逆，僅在MSM的基礎(chǔ)上考慮高階模態(tài)對系統(tǒng)頻響函數(shù)靈敏度的影響。因此，對于大自由度系統(tǒng)，F(xiàn)AM和SAM相較于DFRM在計算效率上具有明顯的優(yōu)勢，且FAM的效率稍高于SAM。這是因為FAM和SAM分別考慮利用系統(tǒng)矩陣與模態(tài)關(guān)系的第一項和前兩項近似高階模態(tài)對系統(tǒng)頻響函數(shù)靈敏度的影響。

圖9 四種方法計算時間隨自由度數(shù)目變化圖Fig.9 Computational time of different DOF with four methods

4.2 軸向振動桿

本節(jié)考慮一個具有雙指數(shù)阻尼的固定-自由軸向振動桿。如圖10所示，該軸向振動桿可根據(jù)需求劃分為N+1個單元。軸向振動桿單元的單元剛度矩陣和單元質(zhì)量矩陣如下

(43)

式中：ρ表示密度；A為桿的橫截面積；E表示材料的楊氏模量；le=L/N為單元桿的長度。各參數(shù)分別為A=6.25×10-4m2，E=2.1×1011N/m2,ρ=7.8×103kg/m3，桿的總長L=4 m。整體質(zhì)量矩陣M和整體剛度矩陣K可以根據(jù)經(jīng)典的有限元方法組裝得到。

假設(shè)軸向振動桿含雙指數(shù)阻尼模型的表達式為

Cnv=μ1e-μ1(t)Cnv1+μ2e-μ2(t)Cnv2

(44)

其中，松弛因子μ1和μ2為

(45)

Cnv1=αM,Cnv2=βK

(46)

其中

(47)

圖10 含雙指數(shù)阻尼模型的軸向振動桿示意圖Fig.10 Axially vibrating rod with double-exponentialdamping model

4.2.1 計算精度

為了討論不同方法的計算精度，本節(jié)將該軸向振動桿劃分為80個單元，即N=80。分別使用DFRM、MSM、FAM和SAM計算該軸向振動桿自由端的頻響函數(shù)H11關(guān)于ρ的靈敏度分析，選取的頻率范圍為：0～10 000 rad/s，選取的頻率步長為10 rad/s。本算例仍然以DFRM計算出的頻響函數(shù)靈敏度為參考值。

圖11給出使用前兩對彈性模態(tài)計算出系統(tǒng)頻響函數(shù)靈敏度，各方法相對于DFRM的相對誤差如圖12所示。在所選頻率范圍內(nèi)，F(xiàn)AM和SAM計算出頻響函數(shù)靈敏度的精度優(yōu)于相同情況下MSM的計算精度。對于FAM和SAM，二者幅值的相對誤差相差無幾，而SAM的相位相對誤差明顯比FAM的相位相對誤差小。因此，就計算精度而言，SAM的計算精度高于FAM計算精度。這與質(zhì)量-彈簧系統(tǒng)算例得到的結(jié)論相同。

圖11 前兩對彈性模態(tài)近似頻響函數(shù)靈敏度

4.2.2 計算效率

本小節(jié)通過改變軸向振動桿的自由度討論DFRM、MSM、FAM和SAM的計算效率。將軸向振動桿自由度劃分為N=0～800，自由度間隔為80。假定頻率范圍為0～10 000 rad/s，頻率步長為100 rad/s。圖13是四種方法計算效率對比圖?？梢钥闯?，F(xiàn)AM和SAM的計算效率明顯優(yōu)于DFRM，而FAM的計算效率高于SAM。各方法的計算效率排序與質(zhì)量-彈簧系統(tǒng)算例得到的結(jié)論相同。

圖12 前兩對彈性模態(tài)近似頻響函數(shù)靈敏度的相對誤差Fig.12 The relative error of FRF sensitivity approximated in terms

圖13 不同方法計算時間隨自由度數(shù)目變化圖Fig.13 Computational time of different DOF with different methods

5 結(jié) 論

本文提出了具有卷積型非黏滯阻尼模型結(jié)構(gòu)系統(tǒng)頻響函數(shù)靈敏度的改進計算方法。利用復(fù)模態(tài)疊加法和直接微分法推導(dǎo)出了非黏滯阻尼系統(tǒng)頻響函數(shù)靈敏度表達式。但若只用有限數(shù)目的低階模態(tài)求解其頻響函數(shù)靈敏度會引入模態(tài)截斷誤差，在某些情形下會嚴重影響非黏滯阻尼系統(tǒng)頻響函數(shù)靈敏度的求解精度。解決該問題的關(guān)鍵在于建立的系統(tǒng)模態(tài)和系統(tǒng)矩陣間的關(guān)系。本文利用Neumann級數(shù)展開定理，將系統(tǒng)高階模態(tài)對系統(tǒng)頻響函數(shù)靈敏度的貢獻用系統(tǒng)矩陣和低階模態(tài)進行表達，進而推導(dǎo)出了兩種考慮高階模態(tài)影響的系統(tǒng)頻響函數(shù)靈敏度的改進計算方法：一階近似法和二階近似法。通過數(shù)值算例驗證兩種方法的有效性和準(zhǔn)確性，結(jié)果表明：

(1) 一階近似法和二階近似法能得到滿意的系統(tǒng)頻響函數(shù)靈敏度計算結(jié)果，因此本文提出的兩種改進計算方法可有效地求解非黏滯阻尼系統(tǒng)頻響函數(shù)靈敏度。

(2) 一階近似法和二階近似法由于考慮了高階模態(tài)對系統(tǒng)頻響函數(shù)靈敏度的影響，因此其計算精度優(yōu)于模態(tài)疊加法；從計算效率角度看，改進的一階近似法和二階近似法避免了矩陣求逆，因此其計算效率高于直接頻響法。綜上，本文提出的兩種改進計算方法兼顧了計算精度和效率，是求解多自由度非黏滯阻尼系統(tǒng)頻響函數(shù)靈敏度問題的理想選擇。

(3) 由于一階近似法和二階近似法分別考慮利用系統(tǒng)矩陣與模態(tài)關(guān)系的第一項和前兩項近似高階模態(tài)對系統(tǒng)頻響函數(shù)靈敏度的影響，因此二階近似法比一階近似法的計算精度更高，而計算效率卻低于一階近似法。綜合考慮計算精度和效率，建議根據(jù)需要選取不同計算方法：當(dāng)需要較高計算精度時，選取二階近似法；當(dāng)需要較高計算效率時，選取一階近似法。