陳 睿, 劉 杰, 張 正, 張 偉

(1. 湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410082；2. 吉首大學(xué)物理與機(jī)電工程學(xué)院， 湖南 吉首 416000 )

引 言

在機(jī)械、結(jié)構(gòu)系統(tǒng)的力學(xué)計(jì)算、結(jié)構(gòu)設(shè)計(jì)、故障診斷中，載荷的確定是一個(gè)非常重要的問(wèn)題，它為結(jié)構(gòu)的設(shè)計(jì)、計(jì)算以及分析提供可靠的載荷依據(jù)，為減小振動(dòng)、提高結(jié)構(gòu)的可靠性、安全性，提供確切的環(huán)境條件。然而在很多實(shí)際工程中，如導(dǎo)彈及航天器在飛行時(shí)的載荷、發(fā)電站工作時(shí)機(jī)組及脈動(dòng)水壓力、橋梁所受運(yùn)動(dòng)態(tài)載荷等，由于技術(shù)條件或經(jīng)濟(jì)條件，對(duì)結(jié)構(gòu)所受外載荷難以直接測(cè)量甚至根本無(wú)法實(shí)測(cè)。因此，近年來(lái)基于結(jié)構(gòu)測(cè)取的響應(yīng)識(shí)別結(jié)構(gòu)所受外載荷的技術(shù)取得了很大的進(jìn)步。

載荷識(shí)別在結(jié)構(gòu)動(dòng)力學(xué)中屬于第二類反問(wèn)題，它根據(jù)已知系統(tǒng)的動(dòng)態(tài)特性和實(shí)測(cè)的動(dòng)力響應(yīng)反求出結(jié)構(gòu)所受的動(dòng)態(tài)激勵(lì)。動(dòng)態(tài)載荷識(shí)別的研究最早可追溯到20世紀(jì)70年代末，從軍事用途發(fā)展起來(lái)的[1]。目前動(dòng)態(tài)載荷識(shí)別的方法主要有頻域法和時(shí)域法兩大類方法[2,3]。頻域法發(fā)展較早，是比較成熟的識(shí)別方法，其應(yīng)用較為廣泛。該方法動(dòng)態(tài)標(biāo)定簡(jiǎn)單、識(shí)別精度較高，但要求測(cè)試數(shù)據(jù)的樣本具有一定的長(zhǎng)度，一般只適用于穩(wěn)態(tài)或平穩(wěn)隨機(jī)載荷識(shí)別，對(duì)沖擊型瞬態(tài)載荷的識(shí)別則受到限制。時(shí)域法提出較晚，它從系統(tǒng)動(dòng)力學(xué)方程出發(fā)，根據(jù)響應(yīng)的時(shí)間歷程直接確定動(dòng)態(tài)載荷的時(shí)間歷程。時(shí)域法直觀，廣泛應(yīng)用于工程中，其優(yōu)點(diǎn)在于時(shí)域法可對(duì)非平穩(wěn)動(dòng)態(tài)載荷識(shí)別，對(duì)沖擊型載荷識(shí)別更是具有實(shí)際工程價(jià)值。目前國(guó)內(nèi)外已開(kāi)展了廣泛研究，并取得了一定的成果[4～6]。

本文從系統(tǒng)動(dòng)力學(xué)方程出發(fā)，基于最優(yōu)輸出跟蹤的基本思想，在載荷位置確定的基礎(chǔ)上，通過(guò)實(shí)測(cè)的動(dòng)力響應(yīng)反求出結(jié)構(gòu)所受的動(dòng)態(tài)載荷。該方法針對(duì)集中力載荷，實(shí)現(xiàn)了多源動(dòng)態(tài)載荷識(shí)別，并采用L曲線法確定了最優(yōu)輸出跟蹤中性能指標(biāo)的關(guān)鍵參數(shù)。數(shù)值算例結(jié)果表明了本文提出的基于最優(yōu)輸出跟蹤的動(dòng)態(tài)載荷識(shí)別方法的有效性和可行性。該方法在載荷作用的全部時(shí)間域內(nèi)整體進(jìn)行反求，應(yīng)用方便，計(jì)算量較小且效率較高，并具有較強(qiáng)的抑制噪聲能力。

1 載荷識(shí)別的正問(wèn)題

結(jié)構(gòu)動(dòng)力響應(yīng)的有限元離散形式為[7]

(1)

可將式(1)改寫為

(2)

假設(shè)一空間向量

(3)

建立相應(yīng)的系統(tǒng)方程

(4)

該系統(tǒng)結(jié)構(gòu)響應(yīng)的輸出為

y=Dw

(5)

若實(shí)際測(cè)點(diǎn)響應(yīng)量為z，(z=[z1z2…zn]，n為測(cè)點(diǎn)個(gè)數(shù))。該響應(yīng)可以為位移響應(yīng)或速度響應(yīng)。如果實(shí)測(cè)響應(yīng)量z已知，需求結(jié)構(gòu)的外載荷，就轉(zhuǎn)化為最優(yōu)跟蹤問(wèn)題，即尋求最優(yōu)控制u，使得結(jié)構(gòu)響應(yīng)輸出y(t)盡量跟蹤上已知實(shí)測(cè)的響應(yīng)輸出z(t)。

記e為跟蹤誤差，即

e(t)=y(t)-z(t)

(6)

2 基于最優(yōu)輸出跟蹤的載荷識(shí)別

最優(yōu)輸出跟蹤問(wèn)題是尋求一種最優(yōu)控制，既能保證跟蹤誤差盡可能小，又要避免使用取值“過(guò)大”的控制，故選用下式的最優(yōu)化性能指標(biāo)[8]

(7)

式中t0為起始時(shí)間，tf為終止時(shí)間，f表示末值型性能指標(biāo)參數(shù)，Q表示平均誤差度量性能指標(biāo)參數(shù)，R表示最小能量控制參數(shù)，且?t∈[t0,tf],f≥0,Q≥0,R>0。本文假設(shè)這些性能指標(biāo)參數(shù)都為常數(shù)。

由于載荷識(shí)別問(wèn)題中結(jié)構(gòu)的外載荷是肯定存在的，所以肯定存在最優(yōu)控制。本文的載荷識(shí)別方法基于載荷位置已知，并通過(guò)實(shí)測(cè)獲得測(cè)點(diǎn)的響應(yīng)來(lái)反求載荷。此外，該載荷識(shí)別問(wèn)題注重時(shí)間歷程，對(duì)終端誤差不關(guān)注，進(jìn)而可忽略終端響應(yīng)量的誤差，即令f=0。因此，最優(yōu)性能指標(biāo)式(7)可改寫為

(8)

式中r=R/Q。

從式(4)和(8)易知哈密爾頓函數(shù)為

λTAw+λTBu

(9)

u*(t)=R-1BTλ

(10)

由極大值原理知共軛方程和橫截條件分別為

λ(tf)=0

(11)

將式(10)代入式(4)并結(jié)合式(11)得正則方程和端點(diǎn)條件

(12)

根據(jù)式(12)，可以假設(shè)

λ(t)=-p(t)w(t)+g(t)

(13)

微分式(13)并結(jié)合式(12)得

(14)

將式(13)代入式(11)，并聯(lián)合式(14)得

p(t)BR-1BTp(t)]w(t)+[p(t)BR-1BT-

(15)

由于式(15)對(duì)任意w(t)均成立，利用w(t)的任意性可知，p(t)和g(t)分別滿足

(16)

由式(11)及w(tf)的任意性得

p(tf)=0，g(tf)=0

(17)

從而可以得到最優(yōu)跟蹤器為

(18)

進(jìn)而得到最優(yōu)跟蹤器相應(yīng)的響應(yīng)

(19)

因此，要最終求得最優(yōu)跟蹤器u*(t)，即結(jié)構(gòu)所受的外載荷，就要求解式(16)和(17)，即帶終端條件的Riccati微分方程以及帶初始條件的微分方程(19)。

通過(guò)上述方法求解微分方程求得p,g,x，即可求出最優(yōu)跟蹤器，即結(jié)構(gòu)所受的外載荷u*(t)。

3 L曲線法確定關(guān)鍵性能指標(biāo)參數(shù)

利用上述方法進(jìn)行載荷識(shí)別過(guò)程中，需要確定最優(yōu)化性能指標(biāo)中的指標(biāo)參數(shù)，即式(8)中的權(quán)重系數(shù)r。當(dāng)測(cè)量響應(yīng)中存在誤差時(shí)，跟蹤誤差和所獲得的取值都會(huì)不同程度的變化，因此該權(quán)重系數(shù)r=R/Q的選取至關(guān)重要，既要保證跟蹤誤差盡可能小，又要避免取值“過(guò)大”。為了能有效地確定r值，本文采用適應(yīng)性強(qiáng)的L曲線法以獲取最佳的r值。

L曲線法是指用對(duì)數(shù)尺度來(lái)描述殘差范數(shù)和解的范數(shù)的曲線對(duì)比[10]，該方法的典型特征是尺度圖形中出現(xiàn)一條明顯的L曲線，曲線的拐點(diǎn)所對(duì)應(yīng)的參數(shù)當(dāng)作優(yōu)化參數(shù)。式(8)中eT(t)e(t)和uT(t)u(t)都是權(quán)重系數(shù)r的函數(shù)，選擇不同的r值，以eT(t)e(t)的對(duì)數(shù)為橫坐標(biāo)，uT(t)u(t)的對(duì)數(shù)為縱坐標(biāo)，得到許多點(diǎn)(lg(eT(t)e(t))，lg(uT(t)u(t)))，經(jīng)過(guò)曲線擬合得到一條曲線。再利用L曲線法選擇最佳的r值，其中本文采用最短距離來(lái)定位L曲線的最優(yōu)點(diǎn)。

4 算 例

為了驗(yàn)證上述基于最優(yōu)輸出跟蹤的載荷識(shí)別方法的正確性和有效性，下面給出3個(gè)多源載荷識(shí)別的數(shù)值算例。在數(shù)值算例中，測(cè)點(diǎn)的位移響應(yīng)通過(guò)有限元數(shù)值仿真計(jì)算得到。此外，測(cè)點(diǎn)的選擇至關(guān)重要，該測(cè)點(diǎn)的響應(yīng)要對(duì)外載荷具有較強(qiáng)的敏感性，則一般通過(guò)敏感性分析來(lái)選取合適的測(cè)點(diǎn)，且測(cè)點(diǎn)的個(gè)數(shù)要大于等于外載荷的數(shù)目，以避免載荷識(shí)別時(shí)出現(xiàn)不適定性問(wèn)題。載荷識(shí)別前，在仿真計(jì)算得到的位移響應(yīng)中加入一定水平的隨機(jī)噪聲以模擬實(shí)驗(yàn)測(cè)量的響應(yīng)，此時(shí)帶噪聲的位移響應(yīng)可用下式表示

Yerr=Ycal+lnoise·std(Ycal)·rand(-1,1)

(20)

式中Ycal為計(jì)算得到的位移響應(yīng)；std(Ycal)為計(jì)算的位移響應(yīng)的標(biāo)準(zhǔn)差；lnoise為百分?jǐn)?shù)，代表噪聲水平；rand(-1,1)為[-1,1]之間的隨機(jī)數(shù)。

4.1 算例 1

如圖1所示一均勻的四層剪切樓板結(jié)構(gòu)[11]，每一層樓板有m=175.126 kg，k=57.328 kN/m，并在第4層樓板和第2層樓板分別施加指數(shù)衰減載荷F1和三角形載荷F2的剪切力模擬瞬態(tài)沖擊載荷，以及在第1層和第3層分別布置測(cè)點(diǎn)。

圖1 結(jié)構(gòu)簡(jiǎn)圖

按結(jié)構(gòu)動(dòng)力學(xué)的方法可寫出結(jié)構(gòu)對(duì)應(yīng)的質(zhì)量陣和剛度陣：

當(dāng)位移響應(yīng)中未加入隨機(jī)噪聲時(shí)，利用上述載荷識(shí)別的方法進(jìn)行載荷反求重構(gòu)，其結(jié)果如圖2所示。從圖中可以看出，在測(cè)量響應(yīng)沒(méi)有噪聲污染的情況下，本方法可以很好地實(shí)現(xiàn)動(dòng)態(tài)載荷時(shí)程的重構(gòu)。同時(shí)，確定最佳r值的L曲線如圖3所示。

當(dāng)在測(cè)點(diǎn)位移響應(yīng)加入5%水平的隨機(jī)噪聲時(shí)，利用本文所述方法進(jìn)行載荷反求重構(gòu)的結(jié)果如圖4所示。從圖中可以看出，在響應(yīng)數(shù)據(jù)受噪聲影響的情況下，本方法仍可較好地進(jìn)行動(dòng)態(tài)載荷識(shí)別。圖5給出了確定最佳r值的L曲線。

表1給出了當(dāng)測(cè)點(diǎn)位移含有5%水平的隨機(jī)噪聲時(shí)，9個(gè)不同時(shí)間點(diǎn)下識(shí)別的載荷值及其相對(duì)誤差。從表中可以看出，本方法在測(cè)量響應(yīng)受到較低水平噪聲影響下識(shí)別載荷結(jié)果的相對(duì)誤差均在10%以內(nèi)。從上述結(jié)果可看出，本文所述的載荷識(shí)別方法能有效地在時(shí)域上重構(gòu)出施加在結(jié)構(gòu)上的多源載荷。

圖2 多源載荷識(shí)別的結(jié)果(零噪聲水平)

圖3 確定最佳r值的L曲線(0%噪聲水平)

圖4 多源載荷識(shí)別的結(jié)果(5%噪聲水平)

圖5 確定最佳r值的L曲線(5%噪聲水平)

表1 不同時(shí)刻識(shí)別的動(dòng)態(tài)載荷及其相對(duì)誤差(5%噪聲水平)

4.2 算例2

圖6所示為一個(gè)高空索道塔架結(jié)構(gòu)有限元模型。該結(jié)構(gòu)立柱下端材料為等邊角鋼∠200×16，上端為等邊角鋼∠200×14，橫桿、斜桿、撐桿下端為不等邊角鋼∠160×100×12，上端為不等邊角鋼∠140×90×10，其他用不等邊角鋼∠100×80×8，材料特性參數(shù)為密度7 860 kg/m3，彈性模量210 GPa，泊松比0.33，阻尼假設(shè)為比例阻尼，與質(zhì)量陣相關(guān)的系數(shù)為10，與剛度陣相關(guān)的系數(shù)為2×10-4。整個(gè)塔架結(jié)構(gòu)采用梁?jiǎn)卧?，塔架的四個(gè)支腳固支。在圖6中所示的箭頭處分別施變頻正弦載荷F1和隨機(jī)動(dòng)態(tài)載荷F2，這些載荷模擬了高空索道在工作過(guò)程中鋼索對(duì)塔架的垂向作用力。

圖6 高空索道塔架結(jié)構(gòu)的有限元模型

圖6中圓點(diǎn)處標(biāo)出了響應(yīng)測(cè)量的位置，通過(guò)仿真計(jì)算，可以得到測(cè)點(diǎn)的位移響應(yīng)，在仿真得到的位移響應(yīng)中加入10%水平的隨機(jī)噪聲來(lái)模擬實(shí)驗(yàn)測(cè)量的響應(yīng)。利用本文提出的載荷識(shí)別方法進(jìn)行載荷時(shí)間歷程重構(gòu)，其結(jié)果如圖7所示。從圖中可以看出，在測(cè)量響應(yīng)數(shù)據(jù)受到較高噪聲水平影響下，本文方法能夠正確有效地實(shí)現(xiàn)動(dòng)態(tài)載荷的識(shí)別。圖8給出了確定最佳r值的L曲線。

圖7 多源載荷識(shí)別的結(jié)果(10%噪聲水平)

圖8 確定最佳r值的L曲線(10%噪聲水平)

4.3 算例3

圖9所示一個(gè)25桿桁架結(jié)構(gòu)，其中桿(1)～(4)有相同的橫截面積A1=400 mm2，桿(16)～(25)、桿(11)～(15)和桿(5)～(10)的橫截面積分別為A2=500 mm2，A3=600 mm2和A4=800 mm2。橫向和縱向桿的長(zhǎng)度L=15.24 m，桿的彈性模量為200 GPa，泊松比為0.33，密度為2 800 kg/m3。阻尼假設(shè)為比例阻尼，與質(zhì)量陣相關(guān)的系數(shù)為0，與剛度陣相關(guān)的系數(shù)為0.003。連接點(diǎn)12為餃接支座，6，8和10為滾動(dòng)支座。在3個(gè)不同的節(jié)點(diǎn)位置分別施加不同時(shí)間歷程的載荷，其中一種為兩個(gè)周期的正弦載荷，另外兩種分別為一個(gè)周期的三角載荷和半個(gè)周期的三角載荷，其作用位置分別如圖9中箭頭F1，F(xiàn)2和F3所示。

建立結(jié)構(gòu)的有限元模型，采用桿單元?jiǎng)澐志W(wǎng)格，共25個(gè)單元和12個(gè)節(jié)點(diǎn)。選擇測(cè)點(diǎn)響應(yīng)為連接點(diǎn)2處的橫向位移、連接點(diǎn)3處的縱向位移和連接點(diǎn)7的橫向位移。

圖9 25桿桁架結(jié)構(gòu)

為了測(cè)試該算法對(duì)噪聲的適應(yīng)能力，在位移響應(yīng)中加入15%水平的隨機(jī)噪聲作為測(cè)量響應(yīng)，該動(dòng)態(tài)載荷識(shí)別的結(jié)果如圖10所示。從圖中可以看出，在測(cè)量響應(yīng)數(shù)據(jù)受到較高噪聲水平影響下，本文載荷識(shí)別的方法可以較準(zhǔn)確地實(shí)現(xiàn)動(dòng)態(tài)載荷時(shí)程的重構(gòu)。圖11給出了確定最佳r值的L曲線。

圖10 多源載荷識(shí)別的結(jié)果(15%噪聲水平)

圖11 確定最佳r值的L曲線(15%噪聲水平)

表2給出了在15%噪聲水平下3個(gè)不同時(shí)刻識(shí)別載荷的相對(duì)誤差。識(shí)別載荷的相對(duì)誤差都在12%以內(nèi)，這說(shuō)明本文所述的多源載荷識(shí)別方法能有效地抑制噪聲對(duì)識(shí)別結(jié)果的影響，具有較好的穩(wěn)健性。

表2 15%噪聲水平下不同時(shí)刻識(shí)別載荷的相對(duì)誤差

5 結(jié) 論

工程實(shí)際中動(dòng)態(tài)載荷的識(shí)別是一個(gè)難度較高且較為復(fù)雜的反問(wèn)題，有許多共性問(wèn)題需要亟待解決。本文基于最優(yōu)輸出跟蹤的思想，提出了一種多源動(dòng)態(tài)載荷識(shí)別的方法，并采用L曲線法確定了性能指標(biāo)中的關(guān)鍵參數(shù)。數(shù)值算例表明該方法能有效抑制噪聲對(duì)載荷識(shí)別結(jié)果的影響，能較為準(zhǔn)確、穩(wěn)健地實(shí)現(xiàn)載荷的反演識(shí)別。該方法計(jì)算簡(jiǎn)單快捷、抗噪聲能力較強(qiáng)，具有一定的工程實(shí)際實(shí)用價(jià)值。

參考文獻(xiàn)：

[1] Bartlett F D, Flannelly W D. Modal verification of force determination for measuring vibration loads[J]. Journal of the American Helicopter Society， 1979: 10—18.

[2] Liu Y, Shepard WS. Dynamic force identification based on enhanced least squares and total-squares schemes in the frequency domain[J]. Journal of Sound and Vibration, 2005, 282: 37—60.

[3] 劉杰. 動(dòng)態(tài)載荷識(shí)別的計(jì)算反求技術(shù)研究 [D]. 長(zhǎng)沙：湖南大學(xué), 2011.Liu Jie. Research on computational inverse techniques in dynamic load identification [D]. Changsha:Hunan University, 2011.

[4] 韓旭, 劉杰, 李偉杰, 等. 時(shí)域內(nèi)多源動(dòng)態(tài)載荷的一種計(jì)算反求技術(shù) [J]. 力學(xué)學(xué)報(bào), 2009, 41(4): 598—605.Han Xu, Liu Jie, Li Weijie, et al. A computational inverse technique for reconstruction of multisource loads in time domain [J]. Chinese Journal of Theoretical and Applied Mechanics, 2009, 41(4): 598—605.

[5] 伍乾坤, 韓旭, 劉杰, 等. 一種直接求解的動(dòng)態(tài)載荷識(shí)別方法[J]. 應(yīng)用力學(xué)學(xué)報(bào), 2011, 28(2): 201—205.Wu Qiankun, Han Xu, Liu Jie, et al. A dynamic load identification method based on direct solution [J]. Chinese Journal of Applied Mechanics, 2011, 28(2): 201—205.

[6] 章繼峰, 張博明, 杜善義. 最優(yōu)化方法在動(dòng)態(tài)載荷時(shí)程重構(gòu)中的應(yīng)用研究[J]. 振動(dòng)與沖擊, 2006, 25(4): 5—8.Zhang Jifeng, Zhang Boming, Du Shanyi. Application of optimization method to kinetic loading time history reconstruction [J]. Journal of Vibration and Shock, 2006, 25(4): 5—8.

[7] Paultre P. Dynamics of Structures [M]. Wiley-ISTE, 2010.

[8] Anderson B D O, Moore J B. Linear Optimal Control[M]. Prentice-Hall Inc, 1971.

[9] 李榮華, 劉播. 微分方程數(shù)值解法[M]. 北京: 高等教育出版社, 2009.Li Ronghua, Liu Bo. Numerical Solution of Differential Equations [M]. Beijing: Higher Education Press, 2009.

[10] Hansen P C, Leary D P. The use of the L-curve in the regularization of discrete ill-posed problems [J]. SIAM Journal on Scientific and Statistical Computing, 1993, 14:1 487—1 503.

[11] 帕茲 M. (美)著．結(jié)構(gòu)動(dòng)力學(xué)——理論與計(jì)算[M]. 李裕澈, 劉勇生，譯. 北京: 地震出版社, 1993.Paz M. Structural Dynamics-Theory and Computation [M]. Beijing: Publish House of Seismological Bureau, 1993.