基于子結構解耦的連接特性識別新方法1)

范新亮 王 彤 夏遵平

(南京航空航天大學振動工程研究所,機械結構力學及控制國家重點實驗室,南京 210016)

引言

機械結構中廣泛存在著各種連接形式,如螺栓連接、鉚接、焊接以及膠接等.連接結構的動力學特性決定著裝配體的動響應預測精度,因此準確地對其進行辨識顯得尤為重要.而由于連接結構的多樣性和復雜性[1-3],難以建立其一般的解析模型[4],因此常基于子結構解耦技術并利用模態或頻響函數等測試數據對連接特性進行辨識.該技術從裝配體和殘余結構的動力學行為中提取連接結構的動態特性[5],其中基于模態試驗數據的方法中振型的測試誤差極易導致識別結果不理想[6],因此更適于有限元環境下應用;而基于頻響函數試驗數據的方法則由于測試簡單、數據信息充足且精度較高,得到廣泛關注與研究[7-12].

Okubo 和Miyazaki[13]最早將基于頻響函數的子結構解耦技術應用于連接特性識別,經過不斷地發展與改善,許多方法被提出,并在一些實際結構中得到了應用[14-17].雖然形式不盡相同,但幾乎所有方法的辨識公式都與Ren 和Beards[18]提出的基本公式相似或等價[19],其思路均為利用實測的裝配體頻響函數與實測或有限元模型計算的殘余結構頻響函數來獲取連接特性.連接特性識別所面臨的主要難題之一是裝配體結構界面自由度頻響函數的測量[19].對此,?eli?等[4-5]在Ren 等所提基本公式上進行了改進,詳細討論了忽略部分界面自由度及內部自由度對辨識結果的影響,并對比了逐頻求逆的直接法與對各個頻率點進行平均的最小二乘法,然而其在忽略部分界面自由度時的辨識結果不佳.此外,采用旋轉自由度的測試技術[16]及未測頻響函數估計技術[19-20]是解決該問題的可行方案,但卻增加了測試或計算的成本,且額外地引入了測試誤差或估計誤差.連接特性識別的另一難題是辨識方程對測試噪聲的高度敏感性,即測試數據微小的誤差也可能導致結果與真實特性完全無關,且不同形式的求解方程在辨識精度上有著較大差異[18],這是由過多的矩陣求逆以及系統方程本身的嚴重病態造成的.對含噪數據進行降噪是減小噪聲干擾最直接的方法,Peeters 等[7]采用模態模型對測試頻響函數進行濾波并應用于子結構解耦,盡管測試噪聲得以濾除,但濾波后的頻響函數對真實值的擬合程度仍難以控制.替代降噪技術的有Zhang 等[21-22]提出的一種噪聲自校準參數辨識方法,以諧波基對噪聲校準量參數化并與系統參數通過誤差函數共同進行估計,有效提高了參數辨識的準確性.在優化辨識公式的形式方面,Wang 等[23]減少了對噪聲有放大作用的矩陣求逆運算,使得方法對實測數據的適用性增強,但由于公式中的簡化使得在高頻時辨識效果不佳.Wang等[19,24]根據基本公式推導了統一形式的辨識方程,并利用該方程將實測與估算的未測頻響函數一起引入辨識過程抑制噪聲,提高了識別結果的準確性,然而由于采用了逐頻求逆的直接法,在由動剛度矩陣序列估計連接結構物理參數的過程中將產生偏差.

為此,本文在子結構解耦基本方程基礎上,提取其測試自由度分量并推導了顯含連接動剛度矩陣的辨識方程,該方程直接利用可測的裝配體頻響函數進行辨識,而無需針對其界面自由度信息進行未測頻響函數的估計或旋轉自由度及部分不可測界面自由度的測量.而后將該方程寫為增量格式,得到了連接特性識別的迭代形式,其殘差相對全量格式顯著減小,因而數值穩定性更好[25],且通過迭代可使得由連接動剛度與殘余結構頻響函數重構的裝配體頻響函數擬合值不斷逼近于測量值.最后,將動剛度曲線以多項式進行擬合,得到了擬合系數的估計方程,并篩選合適的頻率點聯立方程組求解,減少了病態矩陣的求逆運算且改善了抗噪性.彈簧質量系統數值算例和雙梁連接結構試驗表明該方法辨識效果較好,能有效地對實際連接特性進行識別.

1 理論背景

1.1 連接結構動力學特性識別方程

連接結構、殘余結構及其組成的裝配體結構如圖1 所示,令I 代表殘余結構或裝配體結構部分內部節點的自由度,J 代表殘余結構或連接結構全部界面節點的自由度,A 和R 分別代表裝配體結構和殘余結構,則有輸入-輸出關系

圖1 連接結構,殘余結構與裝配體結構Fig.1 Joint,residual and assembly structure

其中XI,XJ與FI,FJ分別為裝配體結構在內部自由度上、界面自由度上的位移矢量與的力矢量,而X與F為由其組合而成的位移矢量與力矢量,頂標-與^則分別代表殘余子結構與連接結構的相應物理量.HA,HR和HJ分別為裝配體結構、殘余結構和連接結構的頻響函數矩陣.

由界面自由度上位移與力的協調關系有

可得基于子結構解耦識別連接特性的基本公式[19]

其中PR,J為篩選HR中界面自由度所在列的布爾矩陣,即PR,JTHR=[HRJIHRJJ]

連接結構動力學特性識別問題一般忽略連接結構的內部自由度(需滿足連接件內部自由度上無外力的假設),且由XJ所指為完備界面自由度知

設實際對裝配體結構進行振動測試的自由度上位移矢量為X*,且X*=P*TX,其中P*亦為布爾矩陣.取P*T左乘式(5)得

式(8)經移項后左乘廣義逆(P*THRPR,J)+得

注意此處求取廣義逆必須滿足行大于列的條件,即要求測試自由度的數目應大于界面自由度數目.上式兩端左乘HC并將式(6)和式(7)代入得

移項并在等式兩端左乘DJ得

而由AX=B的極小二乘解為X=(A)+B知上式等價于

將式(9)記為

式(10)與式(8)形式一致,其中左、右等效頻響為

式(10)即為顯含連接動剛度DJ的辨識方程,是后文推導增量格式以及引入多項式擬合的基礎.其中HR通過校準的殘余子結構有限元模型計算得到,因此對式(10)兩端右乘裝配體結構振動測試的激勵自由度布爾矩陣PI后,辨識過程所需測試數據僅為裝配體在所選測試自由度上的頻響函數P*THAPI,而無需測量或估計其界面自由度的量.且式(10)由式(5)提取測試自由度上的平衡方程后經矩陣變換得到,界面協調條件仍成立,因此仍是精確的.此外,所選測試自由度能準確辨識DJ的必要條件為擁有足夠多有效信息.

1.2 識別方程增量形式

連接結構界面自由度數目較多時,由于待求未知量的數目增加而使得辨識精度明顯下降,且噪聲對辨識結果的擾動顯著增強.為提高方程數值穩定性并通過迭代逼近最佳辨識值,將式(10)變換為增量形式.若給定連接結構動剛度初值DJ0,欲識別其真實值DJ,由式(10)分別有

又由式(11)知右等效頻響函數的差量為

將其代入式(14)可得

其中I為單位矩陣,將上式記為

因此,由式(5),式(16)和式(17)可得識別連接特性的迭代方程

其中各迭代量為

式(18)和式(19)右乘PI后,給定有限元模型計算數據HR、裝配體測試數據P*THAPI與初值后,即可由迭代更新連接結構動剛度.

全量形式的辨識公式(10)以殘余結構與裝配體頻響函數的殘差為極小化目標函數,增量式(18)則以范數遠小于前者的裝配體頻響函數測試值與初始擬合值的殘差為目標函數,參數識別的數值穩定性得到提高[25],且減小了噪聲對辨識結果的干擾.

1.3 動剛度曲線多項式擬合

若采用直接法[5]逐個頻率點依次求解式(18),對含噪聲矩陣(且在界面自由度較多時常因測試自由度信息不足而呈現病態)求逆后將使誤差放大.且連接結構對裝配體結構動力學特性的影響在部分頻率點幾乎可以忽略,使相應頻率點的辨識結果對噪聲極為敏感[16].為此,以多項式函數擬合動剛度曲線,代入式(18)得到關于擬合系數的方程,并選擇合適的頻率點聯立求解,則原問題轉化為頻域內極小化殘差函數的參數估計問題.

連接結構動剛度矩陣第m行第n列(m,n=1,2,···,NJ,NJ為界面自由度數目)對應元素DJmn(ωk)可用多項式函數表示為(1+iωkβ)α0mn-ωk2α2mn,其中β為剛度比例阻尼系數,ωk為第k個頻率點,而擬合系數α0mn,α2mn實際上為連接結構剛度矩陣、質量矩陣的對應分量,則DJ的增量可寫為

取ej為振動測試時篩選位移矢量X中第j個激勵自由度所在行的位置向量,在式(18)兩端右乘ej得

將式(20)代入(21)得

對上式兩端進行矩陣的列拉直得

式中V(·) 表示矩陣的列拉直即為第g個迭代步的殘差向量,而為殘差關于擬合系數的靈敏度矩陣,γ(g)則為擬合系數增量.

式(22)對應于頻率點ωk以及第j個激勵自由度,對k(k=1,2,···,Nf,Nf為頻率點數目)及j(j=1,2,···,NI,NI為激勵自由度數目)組合可得頻域上多輸入多輸出的參數估計方程組f(g)=S(g)γ(g).當所聯立的頻率點及激勵自由度包含足夠多系統有效信息時,可有效減小待估計參數較多(即連接結構界面自由度較多)時靈敏度矩陣的病態程度,且通過聯合全頻域上所有測點與激勵點的信息可提高辨識的準確性.

方法迭代收斂性的證明[21]過程為,若所選測試自由度集、激勵自由度集、頻率點集使得形成的S(g)滿秩,記且

αT為真實擬合系數時,相應殘差f(αT)為0,即αT是迭代方程(23)的不動點.而

且由S(g)滿秩即S(g)TS(g)可逆及有界性知存在常數M1,M2使得

故||α(g+1)-αT||≤M1M2||α(g)-αT||2=d(g)||α(g)-αT||,當任意α(g)滿足||α(g)-αT||＜1/M1M2時,有d(g)＜1,迭代方程(23)為壓縮映射,α(g)將收斂至不動點αT.由此表明給定了合適的連接特性初值后,方法將收斂于真實值.

1.4 算法實現

1.4.1 方程歸一化及頻率點的選擇

式(22)中系數矩陣Sj的各列之間存在顯著的幅值差異,導致其條件數較大,因此按式(24)(25)對其各列進行歸一化

其中[],r與[]r,分別代表矩陣第r列與第r行,norm(.)代表對列取范數.

實際測試試驗中,共振區受噪聲污染小,且對參數變化靈敏度更高[26],故信號信噪比較低時可僅選擇共振區的頻率點參與辨識過程,當信噪比較高時則可選擇給定帶寬的共振區鄰域.此外,由于式(22)以裝配體頻響函數測量值與擬合值的殘差范數為極小化目標,因此可將其幅值相關性系數[27]作為篩選頻率點的準則: 相關性較好的頻率點,其殘差較小且可信度較高;反之,相關性較低則說明該頻率點可能受噪聲污染較為嚴重,且殘差較大不利于準確辨識.幅值相關性系數αa計算式為

其中,? 代表取實部,σmin為所允許的相關性最小值,實踐中發現取σmin＞0 時識別效果較好.

1.4.2 連接結構動力學特性識別流程

連接特性辨識流程如圖2 所示,具體過程如下.

圖2 連接結構動力學特性辨識流程Fig.2 Workflow of the joint dynamic properties identification

步驟1: 確定完備界面自由度及其對應的布爾矩陣PR,J,裝配體結構的測試自由度及其對應的布爾矩陣P*,激勵自由度及其對應的位置向量ej(j=1,2,···,NI).

步驟2: 根據激勵自由度ej及測試自由度P*對裝配體結構進行振動測試,得到測試頻響函數P*THAej.選擇共振區(或其鄰域)內的頻率點序列L1并利用校準后的殘余子結構有限元模型計算相應頻響函數HRP*,HRPR,J.置g=0.

根據式(26)在 L1中選擇參與辨識過程的頻率點 L2.

步驟4: 由式(11)、式(19)計算所選頻率點序列 L2對應的左、右等效頻響函數HEl*,HEr*ej及轉換矩陣TE*(g) .

步驟5: 計算頻率點序列 L2中ωk處的靈敏度矩陣及殘差向量代入式(22),對k和j組合形成超定方程組f(g)=S(g)γ(g),并進行歸一化.求解擬合系數增量γ(g)后得到ΔDJ(g),更新連接動剛度矩陣DJ(g +1)=DJ(g)+ΔDJ(g),置g=g+1.

步驟6: 重復步驟3～5,直至滿足停止準則,即||f(g)||≤ε,輸出連接結構動剛度矩陣的擬合系數.

當殘余結構有限元模型規模較大時,其頻響函數HRP*,HRPR,J的計算量較大且是決定算法效率的關鍵一環,因此需在迭代前計算頻率點序列 L1對應的HR.且該模型須經過精細的校準,否則HR將包含較大的誤差,導致錯誤的識別結果.

2 算例分析

2.1 仿真算例

為驗證所提方法的正確性及抗噪性,以具有10 個平動自由度的簡單線性多自由度系統(圖3)進行數值仿真.將其視為由連接結構與兩個殘余子結構組成的裝配體(殘余子結構1、2 分別在自由度處與連接結構耦合),且該連接結構同時包含質量、剛度及阻尼效應.系統各參數為1.0×107N/m,殘余結構與連接結構的剛度比例阻尼系數β均取1.0×10-7.

圖3 連接結構與殘余結構耦合Fig.3 Coupling of joint and residual structure

裝配體結構計算的頻響函數添加噪聲后作為測試值HAej,殘余結構計算的頻響函數作為所需HR.連接動剛度矩陣所取初值DJ0為

在無噪聲情形下,由本文方法利用測試頻響函數P*THAej(P*不包含界面自由度按1.4 中所述流程進行辨識,并與利用完備頻響函數HAej由逐頻求逆的直接法[5]進行對比.由圖4(a)知兩種方法辨識的DJ與HR重構的裝配體頻響函數均與真實值完全一致,驗證了所提方法僅利用所選測試自由度(無需包含完備的界面自由度)上的裝配體頻響函數識別連接特性的正確性.

為檢驗所提方法的抗噪性,在前述算例的頻響函數真實值中加入10%的白色噪聲與10%的有色噪聲后再次進行識別,由圖4(b)知兩種方法在測試噪聲較小的共振區識別效果均較好,而在其余噪聲較大的頻率點,所提方法的準確度高于直接方法.對比兩者所識別DJ的對角元素與非對角元素可知(圖5),本文方法辨識結果受噪聲的影響遠小于直接法.

圖4 裝配體結構頻響函數重構值與真實值對比Fig.4 Comparison of reconstructed and accurate assembly FRF

圖5 所辨識連接結構動剛度實部的對比Fig.5 Comparison of real part of the identified joint dynamic stiffness

2.2 實驗驗證

2.2.1 膠接

為驗證所提方法識別實際連接結構的有效性,首先對一個T 形結構進行了實驗.該裝配體結構由兩根相同的長為500 mm、寬為50 mm 的梁(殘余結構)在與一鐵塊在界面處膠接而成,且忽略連接界面的非線性行為.實驗中采用了N-Modal 模態測試系統,PCB 力錘以及PCB 三軸傳感器.測試頻率范圍設置為1000 Hz,取1600 條譜線.如圖6 所示,裝配體結構采用彈性繩懸掛的方式模擬自由-自由邊界條件,傳感器位于梁I 中心位置,在均勻分布的126 個測試點上進行錘擊,獲得其頻響函數P*THAPI.辨識公式中(P*THRPR,J)+要求測試自由度數目大于連接結構自由度數目,然而為抑制噪聲[24]并避免P*THRPR,J出現零列或呈現病態,應布置更多包含結構有效信息的測試自由度.

圖6 含一處連接的耦合結構及其測試裝置Fig.6 The coupled structure with one joint and its test setup

采用Shell63 板單元對殘余結構建模,并按材料參數劃分為兩個組.利用自由狀態的測試數據對單根梁進行模型修正[27-31]后,有限元模型計算的頻響函數與其測試值匹配程度較高(圖7),由此得到校準的殘余結構有限元模型.其材料參數如表1 所示,其中D與E代表密度與彈性模量,且殘余結構與連接結構的剛度比例阻尼系數均取為1.0 × 10-7.由該校準模型計算得殘余結構頻響HRP*和HRPR,J.

表1 梁I 和II 的模型材料參數Table 1 The material parameters of Beam I and Beam II

圖7 校準的殘余子結構有限元模型計算頻響與測試值對比Fig.7 Comparison of the residual structure's FRF: calibrated FEM and experimental model

殘余結構與連接結構的耦合如圖8 所示.連接結構兩端的虛擬節點分別通過MPC 約束依賴于梁I、梁II 被膠接區域內的節點,使得虛擬節點的平移及轉動自由度上的位移代表了膠接區域內節點的位移均值,且虛擬節點上的力被均勻分配至膠接區域內的節點上[1].初值DJ0的取法為: 以密度為90 000 kg/m3、彈性模量為3.0×107Pa 的梁單元連接兩虛擬節點,其形成的動剛度矩陣作為與HR所重構的裝配體頻響函數初始擬合值與測試值對比如圖9 所示.

圖8 T 形梁的連接結構與殘余結構耦合Fig.8 Coupling of Joint and residual structure of T-Beam

圖9 迭代前裝配體頻響函數重構值與測試值對比Fig.9 Comparison of reconstructed and experimental assembly FRF before iteration

按1.4 中所述流程識別得到的DJ擬合系數常數項、二次項作為連接結構的剛度、質量矩陣,與殘余結構有限元模型重組可得到修正的裝配體有限元模型,其頻響函數結果與測試值匹配較好(圖10),表明本文方法能有效識別實際連接結構動態特性.值得注意的是,由于實際連接結構剛度、質量矩陣的對稱性、稀疏性及其他特性將對DJ的擬合系數產生一系列約束,需在辨識方程中予以考慮.迭代過程中共振區序列 L1中各頻率點的裝配體重構頻響與測試頻響的幅值相關性如圖11 所示,可見其隨著迭代步的增加而不斷提高,即殘差范數||f(g)||逐漸減小,驗證了算法的收斂性.

圖10 迭代后裝配體頻響函數重構值與測試值對比Fig.10 Comparison of reconstructed and experimental assembly FRF after iteration

圖11 頻響函數相關性的迭代歷程Fig.11 Iterative process of FRF correlation

2.2.2 螺栓連接

為進一步驗證所提方法對較復雜連接動力學特性的識別能力,對一個包含兩處螺栓連接的雙梁結構進行了實驗.該裝配體結構由兩根相同的長度為500 mm 的L 形梁(殘余結構)通過螺栓連接而成,測試中使用的數據采集系統如圖12 所示.測試頻率范圍設置為2500 Hz,取1600 條譜線.裝配體結構以海綿墊支撐以模擬自由-自由邊界條件,三軸傳感器位于螺栓附近,在均布的138 個測試點上進行錘擊,獲得其頻響函數P*THAPI.采用Shell63 板單元對單根L 形梁建模,并按材料參數劃分為3 個組.經模型修正后得到校準的殘余結構有限元模型,其計算的頻響函數與測試值吻合較好(圖13).材料參數如表2 所示,且殘余結構與連接結構的剛度比例阻尼系數均取為1.0 × 10-7.利用該模型可較準確地計算出包含完備界面自由度信息的殘余結構頻響函數.

表2 梁I 和II 的模型材料參數Table 2 The material parameters of Beam I and Beam II

圖13 殘余結構校準后模型頻響與測試值對比Fig.13 Comparison of the residual structure's FRF: calibrated FEM and experimental model

殘余結構有限元模型與連接結構的耦合如圖14所示,與實驗一類似,通過MPC 約束將梁I、梁II 的螺栓壓緊區域內的節點與連接結構耦合,其中壓緊區域取螺栓直徑的兩倍范圍.該結構包含兩處相同的連接,取約束方程使得兩者連接參數相同,即由于兩處連接無共同自由度,因此總體連接動剛度矩陣為初值的取法為:以密度為7900 kg/m3、彈性模量為1.8×1011Pa的梁單元連接虛擬節點,其動剛度矩陣作為與HR重構的裝配體頻響函數初始擬合值與測量值對比如圖15 所示.

圖14 L 形梁的連接結構與殘余結構耦合Fig.14 Coupling of Joint and residual structure of L-Beam

圖15 迭代前裝配體頻響函數重構值與測試值對比Fig.15 Comparison of reconstructed and experimental assembly FRF before iteration

同樣地,所辨識連接結構與殘余結構重組得到裝配體結構有限元模型,其頻響函數結果與測量值對比如圖16 所示,盡管在550～ 650 Hz 的高頻段有所差異,但150～ 550 Hz 的頻率范圍內均匹配較好.由此證明了所提方法能有效地對較復雜的多處連接特性進行識別.圖17 為動剛度曲線相對初值歸一化的擬合系數隨迭代的變化,可以發現部分參數振蕩較為明顯,裝配體結構動響應對這些參數通常有低敏感性的特點因而易受噪聲干擾且難以準確辨識,也正因該特點,其辨識誤差對對連接結構預測裝配體動力學特性的準確性影響較小[16].

圖16 迭代后裝配體頻響函數重構值與測試值對比Fig.16 Comparison of reconstructed and experimental assembly FRF after iteration

圖17 擬合系數的迭代歷程Fig.17 Iterative process of fitting coefficient

3 結論

本文提出一種新的連接特性識別方法,直接利用可測的裝配體頻響函數進行辨識,避免了旋轉自由度及部分不可測界面自由度的測量或估計.通過具有收斂性質的增量方程增加了界面自由度數目較多時辨識的數值穩定性;并以多項式擬合動剛度曲線,在所選頻率點上聯立方程得到參數估計公式,減少了對病態矩陣的求逆,提高了辨識結果的準確性.通過數值仿真算例了驗證所提方法的正確性及抗噪性,通過含膠接、螺栓連接等結構的試驗驗證了方法識別較復雜的實際連接結構特性的有效性.