矩陣攝動和Newmark-β逐步積分相結(jié)合的不確定性結(jié)構(gòu)載荷識別

李曉旺， 黃 科， 魏廣威， 楊翔飛， 黃 磊

(北京機械設(shè)備研究所，北京 100854)

1 引 言

動載荷識別[1]是結(jié)構(gòu)動力學(xué)研究中的第二類逆問題，即已知系統(tǒng)特性和輸出響應(yīng)反求輸入激勵。快速準(zhǔn)確獲取結(jié)構(gòu)的受力情況可以為結(jié)構(gòu)的健康監(jiān)測提供重要參考。對于確定性結(jié)構(gòu)的載荷識別已經(jīng)涌現(xiàn)出大量算法，歸納起來主要有頻域法[2-4]、時域法[5-7]以及智能算法[8,9]。然而，在實際工程中也存在較多不確定性現(xiàn)象。這些不確定性與材料的理化性能、幾何特性、邊界條件、測量誤差和分析模型等因素有關(guān)。如在工程實際中某個梁或桁架制造完成后發(fā)現(xiàn)原材料內(nèi)部存在細(xì)微缺陷，或者工作人員在制作過程中對長度和外徑等尺寸測量不準(zhǔn)確。在這些情況下，就需要將結(jié)構(gòu)對應(yīng)的物理參數(shù)或尺寸參數(shù)作為不確定性參數(shù)來處理[10]。

目前專門針對不確定性結(jié)構(gòu)進行的載荷識別工作還比較少，現(xiàn)有的識別算法主要分為兩種思路，即蒙特卡洛模擬[11]和攝動理論[12-13]。蒙特卡洛模擬是最早出現(xiàn)的處理不確定性的方法，其原理簡單，只需要在區(qū)間范圍內(nèi)對不確定參數(shù)大量重復(fù)采樣，對每一個采樣點計算，最終獲得載荷的最大值和最小值。蒙特卡洛模擬的缺點是計算量過大，運算效率極低。孫興盛等[13]將所求載荷在不確定性參數(shù)的鄰域內(nèi)展開成一階級數(shù)，分別計算載荷中值和一階攝動量，從而建立了針對不確定性結(jié)構(gòu)進行激勵反演的矩陣攝動方法。然而，該攝動方法需要對載荷時間歷程的每一個時刻都計算一次中值和攝動值。由于在時域內(nèi)動態(tài)激勵的時間歷程一般都含有大量時間步，因此會導(dǎo)致迭代步驟過多，運算量較大。

為了彌補現(xiàn)有的不確定性結(jié)構(gòu)載荷識別方法存在的不足，本文在時域內(nèi)建立了一種矩陣攝動和Newmark-β逐步積分相結(jié)合的算法。首先借助矩陣攝動理論將動態(tài)載荷表示成中值和一階攝動量疊加的形式，然后引入Newmark-β逐步積分法分別反求載荷中值和一階攝動值，最后將載荷中值與一階攝動值進行加減運算，從而確定動態(tài)激勵的上下界。仿真算例驗證了該方法能夠高效準(zhǔn)確地反演出載荷邊界，并有力抵御噪聲干擾。

2 載荷識別算法

2.1 Newmark-β逐步積分算法

對于一個多自由度線彈性結(jié)構(gòu)，振動微分方程為

(1)

假設(shè)阻尼形式為瑞利阻尼[14]，則阻尼矩陣可以表示為

C=c1M+c2K

(2)

(0≤γ≤1)

(3)

(0≤β≤0.5)

(4)

ti + 1時刻的速度和位移可以表示為

(5)

(6)

將式(3)代入式(5)可得

(7)

將式(4)代入式(6)可得

(8)

對式(8)變形可得

(9)

將式(9)代入式(7)可得

(10)

將式(9，10)代入式(11)得

(11)

(12)

(13)

Cv=-C0(C+ΔtK)

(14)

將式(12～14)結(jié)合起來可得

(15)

(16)

令

(17)

式(16)可以簡化成振動離散方程為

Y=HF

(18)

式中

由于矩陣不適定的存在，在獲得振動離散方程之后，采用Tikhonov正則化技術(shù)[15]重構(gòu)式(18)的輸入載荷F。設(shè)系統(tǒng)誤差e為

e=Y-HF

(19)

引入罰函數(shù)的概念J為

J=(eHe)+λ(FHF)

(20)

當(dāng)J對F的一階導(dǎo)數(shù)為0，誤差e達(dá)到最小值。此時激勵F可以表示為

F=(HHH+λI)-1HHY

(21)

式中I為單位矩陣，λ為正則化參數(shù)，正則化參數(shù)的值借助L曲線法[16]計算獲得。

2.2 矩陣攝動理論

設(shè)不確定性結(jié)構(gòu)共含有k個不確定性參數(shù)，任意一個不確定性參數(shù)用bj(j=1,2,…,k)來描述，那么不確定性參數(shù)的集合向量b可以表示為

(22)

基于區(qū)間分析思想，將不確定性參數(shù)bj(j=1,2,…,k)用區(qū)間參數(shù)的形式表示，即

(j=1,2,…,k)

(23)

(j=1,2,…,k)

(24)

(j=1,2,…,k)

(25)

參數(shù)bj(j=1,2,…,k)的不確定性程度由不確定度uj來表征。uj越大，參數(shù)不確定性越嚴(yán)重，uj的計算公式為

(j=1,2,…,k)

(26)

bj(j=1,2,…,k)的攝動部分δj的范圍為

(j=1,2,…,k)

(27)

(j=1,2,…,k)

(28)

對于含有不確定性參數(shù)的振動系統(tǒng)來說，結(jié)構(gòu)特性矩陣M,C,K以及載荷向量F(t)都不再是定常值，而是關(guān)于不確定性參數(shù)集合b的函數(shù)，那么M,C,K和F(t)可改寫成M(b),C(b),K(b)和F(t,b)。

在矩陣攝動理論的推導(dǎo)過程中，首先需要計算結(jié)構(gòu)特性矩陣和載荷向量的中值和一階偏導(dǎo)。M(b),C(b),K(b)和F(t,b)在不確定性參數(shù)的中值bc處的值可表示為

M0=M(bc)，C0=C(bc)

(29，30)

K0=K(bc)，F(xiàn)0(t)=F(t,bc)

(31，32)

M(b)，C(b)，K(b)和F(t,b)關(guān)于任意一個不確定性參數(shù)bj(j=1,2,…,k)的一階偏導(dǎo)數(shù)在bc處的值可表示為

M1,j=?M(bc)/(?bj)

(j=1,2,…,k)

(33)

C1,j=?C(bc)/(?bj)

(j=1,2,…,k)

(34)

K1,j=?K(bc)/(?bj)

(j=1,2,…,k)

(35)

F1,j(t)=?F(t,bc)/(?bj)

(j=1,2,…,k)

(36)

在計算得到特征矩陣和載荷向量的中值和一階偏導(dǎo)之后，可以將M(b)，C(b)，K(b)和F(t,b)在bc的鄰域處近似表示為一階泰勒多項式展開式為

(37)

(38)

(39)

(40)

以上即為結(jié)構(gòu)動力學(xué)的矩陣攝動公式，將上述四個一階泰勒多項式代入振動微分方程(1)，可得

(41)

由于所有不確定性參數(shù)之間互不相關(guān)，所以式(41)成立的條件是每一個求和項均為0，即

(42)

(43)

式(42，43)與振動微分方程式(1)具有相似的形式，可借助2.1節(jié)的Newmark-β逐步積分的載荷幅值識別算法重構(gòu)F0(t)和F1,j(t)(j=1,2,…,k)的時間歷程。最后通過矩陣攝動公式將以上兩部分相加減，即可獲得動態(tài)激勵的上下邊界為

(44)

(45)

3 數(shù)值算例分析

為證明本文載荷識別算法的可行性和有效性，建立了2個含有不確定性參數(shù)的數(shù)值算例。每個算例中借助位移傳感器測量節(jié)點的位移響應(yīng)作為已知條件，然后采用本文算法分別對2個算例的動態(tài)激勵時間歷程上下界進行重構(gòu)，并和傳統(tǒng)的蒙特卡洛法進行對比。

3.1 漸變式載荷識別

算例1建立一個懸臂梁模型如圖1所示，載荷施加在節(jié)點8，位移傳感器分別位于節(jié)點7和節(jié)點10。

懸臂梁各項參數(shù)分別為

(1) 不確定性參數(shù)。彈性模量E=[68,72] GPa。

(2) 固定參數(shù)。泊松比0.33，密度2700×103kg/m3，長寬高尺寸為1000 mm×50 mm×20 mm。

對懸臂梁模型施加1個時間歷程變化比較平緩的時變載荷F1，F(xiàn)1的時間歷程如圖4所示。

圖1 受漸變式載荷激勵的懸臂梁模型Fig.1 Cantilever beam model with slowly changing excitation

在算法實施前需要進行兩步前處理，第一步通過仿真方式獲得節(jié)點動態(tài)響應(yīng)，第二步在動態(tài)響應(yīng)中加入5%的高斯白噪聲。前處理完成后采用本文算法識別動態(tài)激勵F1的中值以及關(guān)于不確定性參數(shù)E的一階偏導(dǎo)數(shù)，識別結(jié)果如圖2和圖3所示。

圖2 F1中值識別結(jié)果Fig.2 Midpoint value identification result of F1

圖3 ?F1/?E的識別結(jié)果Fig.3 Identification result of ?F1/?E

在獲得F1的中值以及關(guān)于E的一階偏導(dǎo)數(shù)的前提下，根據(jù)式(44,45)計算F1的上下邊界。同時，蒙特卡洛法也反演F1的上下邊界作為對照。兩種算法對F1的識別結(jié)果如圖4所示。同時，為了定量分析算法的精確度，選取了5個有代表性的時間點計算對應(yīng)的識別誤差，計算結(jié)果列入表1。

圖4 F1識別結(jié)果Fig.4 Identification result of F1

表1 算例1中F1的上下界偏移量Tab.1 Offset of identified load F1 in Example 1

從圖4可以看出，通過本文算法識別出的載荷上下邊界可以準(zhǔn)確地將F1的實際值包含在內(nèi)，體現(xiàn)了算法的有效性。同時，本文算法識別的載荷邊界與蒙特卡洛法的識別結(jié)果吻合度較好，說明算法的識別精度較高。由于受到5%噪聲的影響，識別的動態(tài)激勵在一些時間點上有輕微波動，但總體時間歷程保持平穩(wěn)，顯示出較強的抗噪性能。由表1可知，F(xiàn)1重建邊界的最大偏移量和最小偏移量分別為25.54%和10.21%。作為對照，蒙特卡洛法的識別誤差比矩陣攝動法稍高但相差不大。

3.2 突變式載荷識別

算例2建立一個桁架模型，載荷和位移傳感器的位置如圖5所示，其中3個位移傳感器所測信號均為節(jié)點的豎直位移信號。

圖5 受突變式載荷激勵的桁架模型Fig.5 Truss model with rapidly changing excitation

桁架的各項參數(shù)如下。

(1) 不確定性參數(shù)。4號桁架的長度l=[599,601] mm，11號桁架的直徑d=[19.5,20.5] mm。

(2) 固定參數(shù)。彈性模量200 GPa，泊松比 0.3，密度7800×103kg/m3，除4號桁架以外的所有水平和豎直桁架的長度為600 mm，除11號桁架以外的所有桁架直徑為20 mm。

對懸臂梁模型施加2個時間歷程變化比較迅速的突變式載荷F2和F3,F2和F3的時間歷程如圖6所示。

前處理過程首先通過仿真方式獲得節(jié)點動態(tài)響應(yīng)，然后在響應(yīng)信號中加入10%的高斯白噪聲。前處理完成后采用本文算法反演動態(tài)激勵F2和F3的中值以及關(guān)于不確定性參數(shù)l和d的一階偏導(dǎo)數(shù)，識別結(jié)果如圖6～圖8所示。

圖6 F2和F3的中值識別結(jié)果Fig.6 Midpoint value identification results of F2 and F3

圖7 ?F2/?d和?F3/?d的識別結(jié)果Fig.7 Identification results of ?F2/?d and ?F3/?d

圖8 ?F2/?l和?F3/?l的識別結(jié)果Fig.8 Identification results of ?F2/?l and ?F3/?l

在求得F2和F3的中值以及二階偏導(dǎo)數(shù)之后，分別采用本文算法推導(dǎo)的式(44，45)以及對照蒙特卡洛法反演F2和F3的上下邊界，兩種算法識別的F2和F3的上下邊界如圖9所示。同時，為了更加直觀地分析算法精度，選取了5個時間點計算識別結(jié)果和真實值的相對誤差，誤差的計算結(jié)果列入表2和表3。

算例2的突變式載荷在時間歷程上比算例1的漸變式載荷更復(fù)雜，不確定性參數(shù)增多，噪聲水平也更高。從圖9可以看出，識別的動態(tài)激勵F2和F3的上下界仍保持魯棒性并將實際載荷包絡(luò)在內(nèi)。同時，本文算法和蒙特卡洛法反演的載荷邊界基本保持一致，體現(xiàn)了較高的識別精度。由表2和表3可知，F(xiàn)2的最大和最小偏移量值分別為 23.09% 和 2.33%，F(xiàn)3的最大和最小偏移量值分別為34.29%和12.25%?？傮w來看，本文算法識別載荷邊界的偏移量值在各時間點上接近于蒙特卡洛法的偏移量值，證明本文方法重建的激勵上下界合理有效。

表3 算例2中F3的上下界偏移量Tab.3 Offset of identified load F3 in Example 2

圖9 F2和F3的識別結(jié)果Fig.9 Identification results of F2 and F3

4 結(jié) 論

針對含不確定性參數(shù)結(jié)構(gòu)的動載荷識別問題，本文在時域內(nèi)推導(dǎo)了一個將矩陣攝動理論和Newmark-β逐步積分法結(jié)合起來的綜合算法，該算法具有如下優(yōu)勢。

(1) 通過矩陣攝動理論將動載荷表示成中值和攝動量相加的一階泰勒展開式，避免了在不確定性參數(shù)的區(qū)間內(nèi)大量采樣計算，顯著提高了運算效率。

(2) 在計算載荷中值和攝動量時引入了Newmark-β逐步積分法并進行了改進，改進之后的Newmark-β逐步積分法本質(zhì)上是無條件穩(wěn)定而在形式上是一個簡單的線性離散方程，從而在不影響求解精度的同時又省略了復(fù)雜的時間步迭代。

數(shù)值算例結(jié)果證明，對于復(fù)雜程度不同的載荷以及不同的噪聲水平，該方法均可以高效準(zhǔn)確地反演出載荷上下邊界，顯示出了較高的識別精度，并具有較強的魯棒性和抗噪性。