基于擬靜態初值的載荷識別數值修正算法

徐　菁， 張　方， 姜金輝， 浦玉學， 蔣　祺

(南京航空航天大學 航空宇航學院 機械結構力學及控制國家重點實驗室，南京　210016)

?

基于擬靜態初值的載荷識別數值修正算法

徐菁， 張方， 姜金輝， 浦玉學， 蔣祺

(南京航空航天大學 航空宇航學院 機械結構力學及控制國家重點實驗室，南京210016)

摘要：推導擬靜態法獲得載荷初值及穩定不發散計算結果，分析獲得初值不準確原因，推導獲得新的載荷識別方法。該方法基于數值原理，利用步步修正達到減少累積誤差效果，可較大程度提高計算穩定性，并將數值迭代修正方法用于有限元模型載荷識別。結果表明，該方法能較準確識別出工程中常見的多種載荷，并具有一定抗噪能力。

關鍵詞：載荷識別；時域；擬靜態法；數值修正

動載荷識別技術為結構動態設計關鍵，屬于結構動力學第二類逆問題。其原理為據已知系統特性及動響應求得結構所受動載荷。研究內容為加載位置識別及動載荷量值識別。目前，動載荷識別技術主要有頻、時域兩類方法，其中時域法可直觀表現載荷隨時間歷程的變化規律。

張方等[1-2]利用廣義正交多項式導出基于多節點時間有限元的載荷識別模型，將廣義正交域線性算子逆運算替代時域的復雜逆卷積關系，通過計算正交多項式系數求解載荷，使復雜結構分布動態載荷識別問題得以解決。Nordstrǒm等[3]采用時延方法將載荷識別過程中存在的不適定問題轉化為適定問題，研究時間延遲選取對激勵力識別精度影響。Hashemi等[4]為識別施加于簡支梁的動載荷引入遺傳算法，仿真結果表明該方法的正確性及對測量噪聲、傳感器位置的魯棒性。Gunawan等[5]采用正則化二次樣條函數擬合系統沖擊載荷，利用基于L曲線的TSVD方法求解動態載荷。該方法雖能提供較精確的載荷識別結果，但不能準確識別沖擊載荷卸載過程。Kazemi等[6]利用應變響應時間積分識別動載荷。結果表明，用應變識別的載荷波動較大，需用正則化方法修正，且結果對測試噪聲較敏感。韓旭等[7]將動態載荷表示為一系列脈沖或階躍函數的疊加，并用零相位濾波器、正則化技術及優化策略實現載荷的穩定重構。郭杏林等[8]針對結構受多點分布載荷識別問題提出精細正則化算法，對載荷識別求解過程中的方程病態問題用正則化奇異值分解技術進行處理。姜金輝等[9-10]研究基于Wilson-θ反分析法的動載荷時域識別方法及θ取值、時間間隔、阻尼比、測量點位置數目及采樣中噪聲干擾等因素對動載荷識別精度影響，并在Wilson-θ反分析法基礎上提出多種修正方案，使計算結果趨于收斂。朱濤等[11]引入Newmark積分獲得系統離散運動方程，并基于動態規劃法及Bellman最優化原理，推導出動態優化載荷識別算法公式，該方法抗噪性能較好。王靜等[12]利用系統響應構造狀態空間方程組，建立精確非遞推連鎖計算格式的時域內動載荷反演模型，通過正則化技術克服模型算子的病態特性及響應測量中噪聲不利影響。

本文通過擬靜態法獲得不發散的識別力初值，并說明改初值偏離真實載荷原因。在此基礎上推導獲得基于數值迭代的修正算法并推廣到多輸入模型；用仿真算例及實驗考察該算法的收斂性及抗噪能力。

1擬靜態法獲得初值

用數值迭代算法前提之一需獲得可靠初值范圍，此范圍須含真實值且大小合適。本文所用擬靜態載荷法即為能方便獲得初值方法，計算過程簡單，結果不發散。

加速度信號為工程中載荷識別常見輸入形式，Wilson-θ法計算由加載力得到加速度響應時具有良好的數值穩定性。取θ>1.37即可保證結果無條件收斂，本文取θ=1.4，結合Wilson-θ法，可推導出擬靜態法求載荷公式。對離散多自由度系統，動力學平衡方程為

(1)

在Wilson-θ法計算中[13]，設外力在t～θΔt時間內線性變化，則有

f(t～θΔt)=f(t)+θ(f(t+Δt)-f(t))

(2)

由此，可得t+θΔt時刻系統等效靜力方程為

(3)

(i=1,2,…,n)

(4)

(5)

求解上式可得di(t+θΔt)，再據式(4)，式(3)的解滿足線性關系，即

(6)

設有m個測點，若已測出加速度響應，據Wilson-θ法假設，有

(7)

(8)

(9)

(10)

該方程組的解可分為3種情況：①m≤n即測點數小于未知載荷數時不能得到唯一解，無法識別；②m=n時可求得唯一解；③m≥n時對位移矩陣求廣義逆，可得一組最小二乘解。

求出{λ}t+θΔt后，據式(4)可求出t+θΔt時刻動載荷值fi(t+θΔt)，由于上一步載荷fi(t)已知，故可識別出下一步載荷為

fi(t+Δt)=

[fi(t+θΔt)-fi(t)]/[θ+fi(t)]

(11)

2數值迭代修正算法

擬靜態法計算簡便，雖對已知加速度響應求解所得載荷計算結果不會發散，但假定結構處于靜載荷狀態，獨立計算每個時間間隔點，未考慮一個時間步長內前一時間點載荷、響應對后一時間點響應影響，故不適用于動載荷識別問題。

為進一步獲得較精確結果，可用數值迭代修正擬靜態法求得載荷初值，使每個時間間隔點載荷更接近真實載荷。其中，二分法與黃金分割法為常用數值迭代算法。由Wilson-θ法知

(12)

(12)化簡為

(13)

若在第b1,…,bi,…bn自由度上施加載荷，在第a1,…,aj,…,am自由度上測量響應加速度，則有

f(t+Δt)={0,…,0,f(t+Δt)b1,0,…,0,f(t+Δt)bj,

(15)

即

?

?

簡寫為

(16)

在該方程組中取n個方程，得

(17)

解方程，消去未知量，只保留f1得

(18)

式中：Ai1為系數矩陣去掉第i行第1列代數余子式；其余未知量用f1表示，則有

(19)

式中：Ain,1j為系數矩陣去掉第i、n行第1、j列代數余子式。令

(20)

對g(f1)求導，即

(21)

由擬靜態法獲得識別載荷的初始時間序列f(t)。設時間點tm對應的力為f(tm)，令a0=-rf(tm)，b0=rf(tm)，r為區間放大倍數，且能滿足g(a0)g(b0)<0，則[a0,b0]為初始含根區間(設f(tm)>0)。據單調函數性質知，此區間內g(f1)零點即為真實載荷。一般，若已獲得含根區間[ak,bk](k=1,2,…)，則令

(22)

用二分法[14]計算時取q=0.5；用黃金分割法計算則q=0.618，若

(23)

(24)

實際計算中所用終止原則為bk-ak≤10-p，p為正整數，有

ε≤(1-q)(bk-ak)=10-p(1-q)

(25)

式(24)可化為

(26)

求解以上不等式得

(27)

所需迭代步數上限為

(28)

式中：[ ]為取整。

可見，步步迭代算法增加的計算量主要取決于每次K*的大小。因此確定所用迭代法后，計算耗時與初始含根區間邊界值a0,b0及精度要求p密切相關，而結果的精確性主要由p決定。

3仿真計算

建立自由-自由梁矩形截面梁計算模型，梁長2 m，截面0.08 m×0.02 m，各階阻尼比0.02，彈性模量72 GPa，材料密度2 700kg/m3，劃分為20個有限元單元，每節點有平動、轉動2個自由度，共42個自由度。結構單元節點示意見圖1。

圖1　結構單元節點示意圖Fig.1 Diagram of structure elements and nodes

3.1算例1

在第29、33自由度分別加載f1=10sin(4πt)N，f2=15sin(8πt)N的正弦力。已知數據為第15、21、35自由度的加速度響應。無噪聲下，取Δt=0.001 s，識別結果見圖2。圖2(a)為靜態擬合算法所得識別力，可見，此擬合算法雖不能獲得準確載荷，但結果未現發散且能穩定在一范圍內，符合數值迭代修正的初值要求。圖2(b)為修正后的識別力，故在精度要求合適前提下，數值迭代算法可獲得較準確的計算載荷。

圖2 二輸入載荷識別對比圖Fig.2Two-inputloadidentification圖3 三輸入載荷識別對比圖Fig.3Three-inputloadidentification

3.2算例2

在上算例基礎上增加一輸入載荷，所測加速度響應增加到5個，即第15、29、33自由度上分別加載f1=8sin(4πt)N，f2=10sin(6πt)N，f3=15sin(8πt)N的正弦力，已知數據為第11、17、21、35、39自由度的加速度響應。3個輸入載荷對比見圖3。

實際問題中往往所遇為非規則正弦加載。對以上算例3個自由度的工程常見集中載荷，考察數值迭代修正法的識別效果。圖4(a)為在第15、29、33自由度上加載半正弦、三角波及方波載荷，圖4(b)為在第15、29、33自由度上加載鋸齒波、組合正弦及沖擊載荷。由圖4知，該加載下結果均未發散，說明此算法對各種載荷具有良好的適應性及收斂性。

工程中噪聲影響不可忽略，具不具備良好抗噪性能是判斷算法穩定與否的重要依據。為驗證數值修正算法的抗噪性能，在兩算例正弦輸入下，對已知加速度響應加入5%高斯白噪聲。圖5(a)、(b)分別為二、三輸入下數值迭代修正算法識別結果。由圖5看出，修正后識別曲線雖與真實結果存在一定偏離，但整體并未出現發散，能基本體現真實加載。

4實驗驗證

實驗模型采用仿真計算的自由-自由梁，將其劃分為20個有限元單元。

主要儀器設備有：NI PXI-4496機箱，NI USB-4431數據采集卡，PCB模態力錘，ICP型加速度傳感器及連接線等。將實驗梁用橡皮繩懸掛于支架上，實驗裝置見圖6。

利用m+p SmartOffice軟件進行模態實驗，獲得梁前幾階固有頻率及阻尼比。建立梁仿真模型、計算固有頻率，并用實驗數據對仿真模型進行修正，使兩者盡量吻合，彈性模量為68 GPa，材料密度2 700kg/m3，模型參數見表1。

用力錘快速敲擊實驗梁，利用Labview Signal Express軟件采集數據，獲得固定位置加速度響應。用其作為修正模型的加載響應，利用數值修正算法回推計算輸入載荷，對比實驗所測得真實輸入載荷，分析算法的可行性。

圖4　多形式加載識別對比圖Fig.4 The identification in multiple types of load

圖5　加噪后的識別結果圖Fig.5 The identification after adding noise

參數第1階第2階第3階第4階第5階第6階第7階第8階第9階固有頻率/Hz實驗值25.771.7141.4235.1348.3486.1645.7828.81032.2計算值25.771.1139.4230.5344.3481.1640.8823.71031.1 模態阻尼比/%2.5011.2870.4552.0470.2210.1680.1580.1590.164

圖6　實驗裝置Fig.6 Experimental setup figure

取第7、15號節點為錘擊加載點，第5、13號節點為響應采集點，采樣率20 kHz。據加速度與模型參數，用數值修正算法計算獲得沖擊載荷曲線，見圖7。

設fpeak為加載峰值，進一步分析實驗結果數據，引入峰值誤差參考量Epeak

(29)

實驗中，Epeak,1=4.12%,Epeak,2=3.24%。

由圖7及誤差看出，識別曲線基本能反映真實加載趨勢。沖擊載荷加、卸載過程雖曲線平直段誤差較大，但總體未現發散趨勢。實驗中造成誤差原因較多，如梁的有限元模型與實際模型并不完全一致。實驗儀器設備亦存在不穩定及固有誤差，且錘擊點亦會有偏離。

圖7　實驗結果曲線Fig.7 The identified curve of experimental result

5結論

載荷識別計算過程中，誤差累積導致的發散難以克服。對此，本文提出載荷識別的新思路，在每個時間步長內利用數值算法進行迭代修正。仿真、實驗結果表明，該修正算法可有效獲得穩定的識別結果，基本收斂于真實加載。

參 考 文 獻

[1] 張方,朱德懋. 基于廣義正交域的一種動載荷識別方法研究[J].南京航空航天大學學報, 1996, 9(1):1-8.

ZHANG Fang, ZHU De-mao. A new theoretical study of dynamic load identification based on generalized polynomial expansion[J]. Journal of Nanjing University of Aeronautics & Astronautics, 1996, 9(1):1-8.

[2] 張方,秦遠田,鄧吉宏. 復雜動態載荷識別技術研究[J].振動工程學報, 2006, 19 (1): 81-85.

ZHANG Fang, QIN Yuan-tian, DENG Ji-hong. Research of identification technology of dynamic load distributed on the structure[J]. Journal of Vibration Engineering, 2006,19(1): 81-85.

[3] Nordstr?m L J L, Nordberg T P. A time delay method to solve non-collocated input estimation problems[J]. Mechanical Systems and Signal Processing, 2004,18:1469-1483.

[4] Hashemi R,Kargarnovin M H.Vibration base identification of impact force using genetic algorithm[J]. International Journal of Mechanical Systems Science and Engineering, 2007,1(4): 204-210.

[5] Gunawan F E, Homma H, Morisawa Y. Impact force estimation by quadratic spline approximation[J]. Journal of Solid Mechanics and Materials Engineering, 2008, 2(8): 1092-1103.

[6] Kazemi M, Hematiyan M R, Ghavami K. An efficient method for dynamic load identification based on structural response[C]//EngOpt 2008 International Conference on Engineering Optimization,Rio de Janeiro, Brazil, 2008: 1-5.

[7] 韓旭,劉杰,李偉杰,等. 時域內多源動態載荷的一種計算反求技術[J]. 力學學報, 2009, 41(4): 595-602.

HAN Xu, LIU Jie, LI Wei-jie, 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): 595-602.

[8] 郭杏林,毛玉明,趙巖,等. 基于Markov參數精細積分法的載荷識別研究[J]. 振動與沖擊, 2009, 28(3):27-30.

GUO Xing-lin, MAO Yu-ming, ZHAO Yan, et al. Load identification based on precise time-step integration for markov parameters[J]. Journal of Vibration and Shock,2009,28(3):27-30.

[9] 徐菁,張方,馬小艷. 時域動載荷識別的若干修正算法[C]//第十屆全國振動理論及應用學術會議論文集,南京, 2011: 796-802.

[10] 姜金輝,徐菁,張方,等. Wilson-θ反分析法的動載荷識別精度的若干問題[J]. 振動、測試與診斷, 2013,33(5):783-788.

JIANG Jin-hui, XU Jing, ZHANG Fang,et al. Accuracy of the load identification based on inverse analysis of Wilson-θmethod[J]. Journal of Vibration, Measurement & Diagnosis,2013,33(5):783-788.

[11] 朱濤,肖守訥,陽光武. 一種新的時域動載荷識別方法[J]. 西安交通大學學報, 2012, 47(6): 968-973.

ZHU Tao, XIAO Shou-ne, YANG Guang-wu. A new time-domain method for force identification[J]. Journal of Southwest Jiaotong University, 2012, 47(6): 968-973.

[12] 王靜, 陳海波,王婧. 基于精細積分的沖擊載荷時域識別方法研究[J]. 振動與沖擊, 2013, 32(20):81-85.

WANG Jing, CHEN Hai-bo, WANG Jing. Impulsive load identification in time domain based on precise time-integration method[J]. Journal of Vibration and Shock, 2013, 32(20):81-85.

[13] 張雄,王天舒. 計算動力學[M].北京:清華大學出版社, 2007.

[14] 藺小林. 計算方法[M]. 西安:西安電子科技大學出版, 2009.

基金項目：航空科學基金(2012ZA52001)；高等學校博士學科點專項科研基金(20123218120005)；中央高校基本科研業務費專項資金(NS2012080)；國家自然科學基金資助項目(51305197);機械結構力學及控制國家重點實驗室自主研究課題(0115K01)

收稿日期：2014-06-30修改稿收到日期：2015-01-22

通信作者張方 男，教授，博士生導師，1962年8月生

中圖分類號:TH 113; O 327

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.02.007

Numerical correcting algorithm for load identification based on quasi-static initial value

XU Jing, ZHANG Fang, JIANG Jin-hui, PU Yu-xue, JIANG Qi

(State Key Laboratory of Mechanics and Control of Mechanical Structures,College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

Abstract:A quasi-static method was presented to obtain the initial value of loading, which is a simple and quick algorithm to get stable and convergent results. The reason why the initial value of the algorithm would be not accurate was analysed, and a new load identification algorithm was deduced. The method is based on the principle of numerical algorithm and uses the idea of step-by-step correction, which can reduce the effect of cumulative error and greatly improve the stability of calculation. This numerical iteration method was applied to the load identification of a finite element model. The results of simulation and experiment show that the method can accurately identify the common variety of loads in engineering, and it has a certain ability to resist noise.

Key words:load identification; time domain; quasi-static method; numerical correction

第一作者 徐菁 女，博士生，1988年9月生