變截面懸臂梁結構動載荷辨識方法

吳肖，曾捷,*，胡子康，李明，胡錫濤

1. 南京航空航天大學 機械結構力學及控制國家重點實驗室，南京 210016 2. 中國航空綜合技術研究所 裝備服務產品部，北京 100028

載荷識別主要用于計算惡劣服役環境下結構主承力構件所受載荷或者無法通過傳感器獲取的外部激勵[1]。動載荷識別一般分為直接測量法和間接識別法2種[2]，直接測量法是通過在飛行器表面設置大量測壓孔感知整個表面的壓力分布。此方式特點是測量精度高、計算量小，但缺點是打孔給飛行器結構引入損傷，一定程度上影響了飛行器結構力學性能。

間接識別算法是在結構表面布置傳感器網絡直接采集載荷響應信息來反求載荷，此種方法不需破壞結構、操作較簡單。國內外許多學者對此進行了大量深入的研究，Hong等提出了一種倒譜域的機械動態載荷譜識別方法，基于高階譜的諧波恢復，建立了ARMA(Auto Regvessive Moving Avevage Model)模型，確定格林函數的系數，求出動載荷，但只是進行了仿真驗證，缺乏工程驗證[3]。Maes等基于桿件的模態特性，假設了Timoshenko梁理論，在考慮桿的轉動慣量和傳感器質量的基礎上，提出了一種估算組合結構桿件軸向力的方法，但不適用于變截面等復雜結構的動載荷識[4]。Ma等提出了一種基于卡爾曼濾波和遞推最小二乘算法的在線遞推動載荷方法。此算法對高斯白噪聲具有很好的抑制作用，通過正弦、三角脈沖、矩形脈沖、一系列脈沖和隨機脈沖仿真，驗證了該方法對由噪聲測量值估計梁結構外部激勵具有良好的性能，但未實現載荷作用位置辨識[5]。間接識別法屬于數學上的逆問題，因而存在矩陣病態導致的識別精度不高和計算量偏大的問題[6]。目前關于間接識別法的研究中，除了矩陣病態問題，還有傳感器的安裝和布局等難點[7-8]。

由于光纖傳感技術能夠實現針對應變、溫度、振動等多種物理量的直接測量，具有便于分布式組網、靈敏度高、芯徑細、重量輕、柔韌性好、適于強電磁干擾以及腐蝕性危險環境使用等獨特優點[9]。基于此，本文研究了一種基于光纖光柵傳感器與卡爾曼濾波器的變截面懸臂梁結構動載荷識別算法。該方法基于自適應卡爾曼濾波器反饋修正的思想進行載荷識別，一定程度上解決了間接識別法中矩陣求逆所導致的病態問題，給出了載荷所用位置判斷依據，并通過仿真與實驗驗證了此算法具有較好的抗噪能力、較高的載荷辨識精度和工程適用性。

1 光纖光柵應變傳感原理

光纖光柵傳感器(Fiber Bragg Grating，FBG)基于光纖周期或纖芯折射率的變化，當光纖Bragg光柵感受到外界溫度、應變變化時，會引起光柵柵區折射率的變化，從而使得其對應中心波長發生改變[10]。光纖光柵傳感原理，如圖1所示。

由光纖光柵傳感器原理可知，假定環境溫度恒定，此時FBG僅受軸向應變因素影響，當中心波長變化量為ΔλB時，其應變表達式為

(1)

式中：Pe為光纖光柵有效彈光系數，取為常數0.22;λB為FBG反射光譜中心波長; ΔλB為中心波長偏移量[11]。

圖1 光纖光柵傳感器原理Fig.1 Principle of fiber Bragg grating sensor

2 動載荷識別原理

2.1 變截面梁有限元分析

針對如圖2所示長度為l、厚度為h、上底邊寬度為b0、下底邊寬度為b1的錐形變截面梁單元形式，其寬度b(x)、截面積A(x)和慣性矩I(x)可以表示為[12]

圖2 錐形變截面梁單元示意圖Fig.2 Schematic diagram of tapered variable section beam element

(2)

根據結構力學的有限元分析可知，等截面梁單元剛度矩陣Ke可以表示為

(3)

對于等截面梁單元的一致質量矩陣Me可以表示為

(4)

式中：EI為梁單元抗彎剛度；A為梁的截面面積；b為梁截面寬度[13]；N為梁單元的形函數矩陣；ρ為材料密度；B為幾何矩陣。將式(2)代入式(3)、式(4)，可以得到變截面梁單元剛度矩陣與質量矩陣。

2.2 單元節點劃分

通過有限元法將梁離散成n個單元，在每個單元內選取2個應變測量點，共計2n個應變測量點。由此，可通過形函數分別建立起每個單元內所測應變與單元節點的位移、轉角之間的映射關系:

(5)

式中：ε1～ε2n為n個單元對應的2n個應變值，實驗中可由光纖光柵傳感器測得；ω1～ωn為n個單元節點的位移；θ1～θn為n個單元節點的轉角；B(i,ξi)=[-6+12ξi,l(-4+6ξi), 6-12ξi,l(-2+6ξi)]×y/l2,ξ為光纖光柵傳感器粘貼位置，y=h/2為梁中性軸距離其上表面的距離。

2.3 建立狀態方程

離散化后的梁結構為二階振動系統，其2n個自由度系統的振動方程為

MX″(t)+CX′(t)+KX(t)=F(t)

(6)

式中：M、C和K分別為2n×2n的質量矩陣、阻尼矩陣和剛度矩陣；F(t)為2n×1的載荷向量；X(t)為系統的位移；X′(t)為系統的速度；X″(t)為系統的加速度。

對于阻尼矩陣在實模態空間內不能解耦的情況，引入狀態坐標變換[14]:

Y(t)=[X(t)X′(t)]

(7)

將式(7)代入式(6)，可把上述矩陣形式的微分方程轉換為狀態空間方程，表示為

(8)

Y′(t)=DY(t)+GF(t)

(9)

Z(t)=HY(t)

(10)

式中：H為觀測矩陣；Z(t)為觀測序列，即結構的應變響應[15]。

(11)

采樣時間間隔設置為ΔT，將式(9)和式(10)離散化，得到線性離散系統的狀態方程和觀測方程：

Y(k+1)=ΦY(k)+Ψ(F(k)+w(k))

(12)

Z(k)=ΗY(k)+v(k)

(13)

Φ=exp(D*ΔT)

(14)

(15)

式中：Y(k)為狀態向量；Φ為狀態轉移矩陣(描述動態系統從時間k到時間k+1狀態的轉移)；Ψ為驅動矩陣；w(k)和v(k)分別為互不相干的系統白噪聲和觀測白噪聲，其系統特性為[16]：數學期望(系統白噪聲)為0，即E[w(k)]=0；系統白噪聲自適應函數為E[w(k)wT(k)]=Qkdkl；數學期望(觀測白噪聲)為0，即E[v(k)]=0；觀測白噪聲自適應函數為E[v(k)vT(k)]=Rkdkl,dkl為自相關函數；Q、R為協方差矩陣,Q=QwI2n×2n，R=RvI2n×2n,Qw和Rv為w(k)白噪聲給定的協方差初始值。

2.4 卡爾曼濾波器原理

卡爾曼濾波器可分為2個部分：時間更新預測和測量更新修正[17-18]。

時間更新：將光纖光柵傳感器所測應變響應作為輸入變量進行時間方程更新。

(16)

P(k/(k-1))=ΦP(k-1/(k-1))ΦT+ΨQΨT

(17)

根據上一狀態的協方差P(k-1/(k-1))和系統過程中的協方差Q，求出當前狀態的協方差預測值P(k/(k-1))。

測量更新：

S(k)=HP(k/(k-1))HT+R

(18)

利用式(17)得到當前時刻的協方差預測值P(k/(k-1))，再結合觀測矩陣H，求得新息協方差S(k)。

Ka(k)=P(k/(k-1))HTS-1(k)

(19)

利用新息協方差S(k)以及觀測矩陣H，得到卡爾曼增益矩陣Ka(k)。

P(k/k)=[I-Ka(k)H]P(k/(k-1))

(20)

(21)

利用觀測方程Z(k)、觀測矩陣H以及當前時間的狀態估計值Y(k/(k-1))，求得當前時間的新息序列Z(k)。

(22)

2.5 載荷識別器原理

基于卡爾曼濾波器的載荷估計器是以光纖光柵傳感器測得的應變響應作為觀測信號，通過卡爾曼濾波器生成的增益矩陣、新息序列以及協方差矩陣，得到靈敏度矩陣和估計力的增益矩陣。在此基礎上，利用廣義回歸模型及其最小二乘算法估計載荷特征。

基于卡爾曼濾波器的動載荷識別器方程如下所示[20]：

Bs(k)=H[ΦMs(k-1)+I]Ψ

(23)

Ms(k)=[I-Ka(k)H][ΦMs(k-1)+I]

(24)

(25)

Pb(k)=[I-Kb(k)Bs(k)]γ-1Pb(k-1)

(26)

(27)

基于光纖光柵傳感器和卡爾曼濾波器的動載荷識別流程，如圖3所示。

圖3 基于光纖光柵傳感器與卡爾曼濾波器的動載荷識別流程Fig.3 Dynamic load identification process based on FBG sensor and Kalman filter

3 有限元數值仿真

3.1 變截面懸臂梁結構建模

利用Ansys Workbench對錐形梁單元形式的變截面梁結構進行有限元數值仿真，材料密度為2 770 kg/m3，彈性模量為71 GPa，泊松比為0.3，上底邊長為100 mm，下底邊長250 mm，模型長度為1 200 mm，厚度為6 mm。

將變截面懸臂梁結構模型平均劃分為3個單元，從左到右依次定義相應3個單元節點：A～C。在每個單元內的1/3處和2/3處設置應變監測提取點，則共有6個應變提取點，從固支端開始依次將其定義為1～6。單元節點劃分與應變提取點位置，如圖4所示，其中矩形標記為單元內應變提取點，圓形標記為單元節點。

為驗證此載荷辨識算法的抗噪能力，在提取應變響應后加上σ=1×10-9，Qw=1×10-15的高斯白噪聲作為系統的觀測序列代入相關算法。設置衰減因子γ=0.65，采樣頻率為500 Hz。

假設變截面懸臂梁結構系統阻尼為比例阻尼，阻尼矩陣是質量矩陣和剛度矩陣的線性組合，可以表示為

C=αM+βK

(28)

式中:α為質量比例常數，β為剛度比例常數，設置α=0.01，β=0.02。

為進一步描述動載荷識別效果，量化動載荷識別誤差，采用平均絕對誤差MAE和均方根誤差RMSE來描述動載荷識別誤差。

(29)

(30)

圖4 單元節點劃分與測點位置圖Fig.4 Unit node division and measuring point location

3.2 有限元仿真結果驗證

3.2.1 單點加載下動載荷辨識效果

1) 動載荷大小識別

依次在單元節點A、B和C施加正弦激勵：F=5sin(3πt)，識別效果如圖5所示。

根據圖5可知，單點載荷作用下，單元節點識別載荷效果良好，載荷識別曲線與仿真施加載荷曲線基本吻合，各單元節點載荷識別誤差相差不大。單元節點A處的動載荷識別精度較好，單元節點B和單元節點C處的載荷識別曲線在峰值處略微偏小，這主要是因為懸臂梁模型在外加載荷激勵過程中，距離固支端越遠的測點，提取的應變響應受到外部激勵擾動越明顯，使得載荷辨識精度相對于固支端有所降低。

單點激勵下各單元節點載荷識別誤差，如表1 所示。

2) 動載荷激勵位置識別

根據式(5)函數關系，本文提出通過提取變截面模型每個單元內2個位置的應變響應來計算相應單元節點處的位移和轉角作為載荷識別算法的輸入變量，而輸出值是每個單元節點對應的載荷與彎矩。因此，算法識別出來的載荷與其相應單元節點呈一一對應關系，也即在載荷激勵位置，將會存在顯著的載荷響應，據此可以判斷出動載荷施加位置。

圖5 載荷識別效果Fig.5 Load identification effect

表1 單點激勵下各單元節點載荷識別誤差

以在單元節點B處施加激勵載荷為例，單元節點B識別的有效載荷和單元節點A、單元節點C識別到的擾動載荷響應曲線，如圖6所示。可以發現，在單元節點B(距離固支端0.8 m)處識別的載荷顯著大于單元節點A(距離固支端0.4 m)、單元節點C(距離固支端1.2 m)，且單元節點A、單元節點C處識別的載荷波動很小，無明顯變化規律，因此可以判斷載荷施加位置為單元節點B處，這與模擬仿真所加載荷位置相符。

圖6 各單元節點載荷識別曲線Fig.6 Load identification curves of each element node

3.2.2 三單元節點同時加載下動載荷辨識效果

為進一步驗證此算法在復合工況下動載荷識別效果，選擇在單元節點A、單元節點B和單元節點C處同時施加3個正弦激勵，具體加載情況如表2所示。

3個單元節點同時加載條件下動載荷識別效果，如圖7所示。從圖7可以看出，3個單元節點同時加載工況下，單元節點A、B、C處對應的載荷識別曲線與有限元仿真結果基本吻合，僅在峰值處均與存在一定偏差。單元節點C處識別的載荷受其他激勵載荷影響，識別曲線出現一定的波動。相對于單點激勵工況，由于加載幅值的增加，3個單元節點的載荷識別誤差略有增大。

表2 三單元節點同時加載工況Table 2 Three-unit joint simultaneous loading conditions

圖7 3個單元節點同時加載下動載荷識別效果Fig.7 Dynamic load identification effect under simultaneous loading of three element nodes

3個單元節點同時加載條件下動載荷識別誤差，如表3所示。

表3 3個單元節點同時加載下動載荷識別誤差

4 實驗驗證

4.1 變截面懸臂梁結構動載荷識別實驗系統

變截面懸臂梁結構單元節點劃分與傳感器布局示意圖，如圖8所示。變截面懸臂梁結構分布式光纖動態載荷識別實驗系統，如圖9所示。

圖8 單元節點劃分與傳感器布局Fig.8 Unit node division and sensor layout

圖9 分布式光纖動載荷監測與識別實驗系統Fig.9 Distributed optical fiber dynamic load monitoring and identification experimental system

實驗平臺主要由光纖光柵傳感器、光纖光柵解調儀、計算機、33500B型信號發生器、HEAS-50型功率放大器、力顯示器、CZL-1012型力傳感器、NI USB-6341數據采集卡、HEV-50型激振器等組成。

4.2 單點加載工況下動載荷辨識效果

在單元節點C和單元節點B進行單點激勵，通過算法識別出來的動載荷與力傳感器直接測得的真實激勵進行數據分析與對比，激勵工況及動載荷識別效果如下。

4.2.1 正弦激勵

單元節點B正弦激勵(7 Hz)動載荷識別效果，如圖10所示。由圖10(b)可以判別載荷激勵位置為單元節點B。

圖10 單元節點B正弦激勵(7 Hz)動載荷識別效果Fig.10 Dynamic load identification effect of unit node B with sinusoidal excitation (7 Hz)

4.2.2 方波激勵

單元節點C處方波激勵(1 Hz)動載荷識別效果，如圖11所示。由圖11(b)可以判別載荷激勵位置為單元節點C。

圖11 單元節點C方波激勵(1 Hz)動載荷識別效果Fig.11 Identification effect of dynamic load excited by C square wave (1 Hz) of unit node

4.2.3 鋸齒波激勵

單元節點C處鋸齒波激勵(1 Hz)動載荷識別效果，如圖12所示。由圖12(b)可以判別載荷激勵位置為單元節點C。

根據圖10～圖12可知，在單點加載工況下，載荷識別曲線與真實激勵曲線吻合度較好，只是在最大幅值處出現一定偏差。由于在施加方波、鋸齒波時激振器會產生一定的頓挫，在頓挫瞬間，激勵本身出現振蕩，使得載荷識別曲線相較于實際施加載荷曲線存在一定波動。

通過各單元節點載荷識別曲線，根據每個單元節點對應的識別載荷可以判斷出真實激勵單元節點。

單點激勵下，不同頻率與波形激勵載荷對應的動載荷辨識誤差，如表4所示。

圖12 單元節點C鋸齒波激勵(1 Hz)動載荷識別效果Fig.12 Identification effect of dynamic load excited by sawtooth wave (1 Hz) of element node C

表4 不同頻率與波形激勵載荷對應的載荷辨識誤差

4.3 兩單元節點同時加載工況下動載荷辨識效果

在單元節點B(激勵頻率2 Hz)和單元節點C(激勵頻率4 Hz)同時施加外部激勵荷載，載荷識別效果如圖13所示。

根據圖13所示，在兩點同時激勵的復合動載荷加載工況中，載荷辨識曲線與實際加載曲線輪廓基本吻合，算法能夠較好地識別出單元節點B和單元節點C處對應的外部激勵載荷特征，具有良好的識別效果。根據圖14可以判斷出單元節點載荷激勵位置，分別為單元節點B和單元節點C。

圖13 單元節點B、C同時激勵下動載荷識別效果Fig.13 Identification effect of dynamic load under simul taneous excitation of B and C nodes

圖14 單元節點載荷激勵位置識別Fig.14 Location identification of load excitation of element node

動載荷識別結果中呈現2個激勵頻率波形相互疊加的情況，這是因為在單元節點B和單元節點C處同時激勵的過程中，2個激振器通過變截面懸臂梁結構相互作用，激勵本身出現耦合所致。這表明該算法針對復合載荷同時加載工況，具有較好的動載荷識別能力。兩單元節點同時激勵實驗誤差，如表5所示。

表5 兩單元節點同時激勵對應的載荷辨識誤差

綜合上述分析，動載荷識別結果也存在一定偏差，分析其原因可能包括：一是系統結構離散的單元數較少；二是光纖光柵傳感器與變截面懸臂梁結構上表面未均勻膠結或粘貼方向有偏差，引起應變測量誤差；三是光纖光柵傳感器中心波長偏移量受結構振動擾動影響較大，尤其是在最大振幅時識別精度有所下降。

5 結 論

本文針對飛行器結構載荷監測需求，推導了變截面梁模型質量矩陣、剛度矩陣等參數信息，提出了基于分布式光纖傳感器與卡爾曼濾波器的變截面懸臂梁動載荷辨識方法。

1) 基于自適應卡爾曼濾波器反饋修正的思想進行載荷識別，一定程度上解決了常規間接識別法中矩陣求逆所導致的病態問題。

2) 借助ANSYS有限元分析軟件，開展針對結構在不同載荷激勵狀態下的數值仿真，驗證了該算法對于多種不同載荷形式辨識的可行性以及良好的噪聲抑制能力。

3) 基于形函數建立的變截面懸臂梁單元所測應變與單元節點的位移、轉角之間映射關系，提出了動載荷激勵位置判別方法。

4) 構建了基于分布式光纖傳感器的變截面懸臂梁應變監測與動載荷識別系統，試驗結果表明，該算法對于復合載荷工況具有較好的識別能力和工程適用性。

5) 本文所述方法采用全分布式光纖傳感網絡，相較于電學式測量方法，具有抗電磁干擾、質量輕不改變被測結構動力學特征、適用于狹小空間機載環境的非視覺測量等特點，能夠為風洞和飛行試驗中飛行器載荷辨識提供技術支撐。