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

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

1 引 言

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

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

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

2 載荷識(shí)別算法

2.1 Newmark-β逐步積分算法

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

(1)

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

C=c1M+c2K

(2)

(0≤γ≤1)

(3)

(0≤β≤0.5)

(4)

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

(5)

(6)

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

(7)

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

(8)

對(duì)式(8)變形可得

(9)

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

(10)

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

(11)

(12)

(13)

Cv=-C0(C+ΔtK)

(14)

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

(15)

(16)

令

(17)

式(16)可以簡(jiǎn)化成振動(dòng)離散方程為

Y=HF

(18)

式中

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

e=Y-HF

(19)

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

J=(eHe)+λ(FHF)

(20)

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

F=(HHH+λI)-1HHY

(21)

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

2.2 矩陣攝動(dòng)理論

設(shè)不確定性結(jié)構(gòu)共含有k個(gè)不確定性參數(shù)，任意一個(gè)不確定性參數(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ì)算公式為

(j=1,2,…,k)

(26)

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

(j=1,2,…,k)

(27)

(j=1,2,…,k)

(28)

對(duì)于含有不確定性參數(shù)的振動(dòng)系統(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òng)理論的推導(dǎo)過程中，首先需要計(jì)算結(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)于任意一個(gè)不確定性參數(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)

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

(37)

(38)

(39)

(40)

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

(41)

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

(42)

(43)

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

(44)

(45)

3 數(shù)值算例分析

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

3.1 漸變式載荷識(shí)別

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 突變式載荷識(shí)別

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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