回傳射線矩陣法在周期結(jié)構(gòu)振動分析中的應(yīng)用

陳 林，布英磊，葉 曦，王青山

（1.中國船舶及海洋工程設(shè)計研究院，上海200011；2.中南大學(xué)，機電工程學(xué)院，湖南 長沙 410083）

1 引言

周期結(jié)構(gòu)是一種新型的復(fù)合周期材料，其特有的振動帶隙特性在工程應(yīng)用上具有較為實際的應(yīng)用價值[1]。在周期結(jié)構(gòu)分類方面，通常根據(jù)周期結(jié)構(gòu)的周期維數(shù)進行劃分，將周期結(jié)構(gòu)劃分為一維，二維和三維單元。因一維周期結(jié)構(gòu)較為簡單，雖然較二維，三維的研究較少，但是與工程應(yīng)用領(lǐng)域較為貼合。目前，計算周期結(jié)構(gòu)振動帶隙的計算方法主要有[2]傳遞矩陣法，平面波展開法，集中質(zhì)量法等，上述算法有著個各自的特點以及計算的領(lǐng)域范圍。例如，傳遞矩陣法是一維周期結(jié)構(gòu)的帶隙特性計算中較為常見的計算方法，其計算量較小，但是較難處理二維，三維的周期結(jié)構(gòu)模型[3]；平面波展開法是周期結(jié)構(gòu)問題研究中較為常見的一種方法，將彈性模量等參數(shù)按照傅里葉級數(shù)進行展開，結(jié)合布魯赫定理，將彈性波的波動方程轉(zhuǎn)化為本征值進行求解，進而得到周期結(jié)構(gòu)的能帶結(jié)構(gòu)[4]；集中質(zhì)量法是利用離散化的物理思想，將集中質(zhì)量引入到二維，三維周期結(jié)構(gòu)帶隙的計算過程之中[5]。上述三種方法均存在著較為明顯的特點以及計算的領(lǐng)域范圍。

對回傳射線矩陣法的原理進行闡述，對其理論思想以及求解過程進行描述，并且給出以一維等截面周期結(jié)構(gòu)對于軸向波的傳輸特性的計算作為算例。

2 算法原理闡述

回傳射線矩陣法（MRRM）是1998年由Howard和Pao公開提出的，隨著時代的進步和社會的發(fā)展，人們對回傳射線矩陣法的研究逐漸深入，并且對其計算領(lǐng)域的范圍逐漸進行深入的擴展研究[6]。近年來，隨著對周期結(jié)構(gòu)理論的逐漸研究，復(fù)合周期結(jié)構(gòu)梁中彈性波的傳輸特性逐漸得到重視，引入回傳射線矩陣法計算周期結(jié)構(gòu)的振動特性曲線。

回傳射線矩陣法基本思想是[7]：將整個結(jié)構(gòu)劃分為若干單元，在各個單元內(nèi)建立對偶坐標(biāo)系，以控制微分方程為基礎(chǔ)，結(jié)合傅里葉變化，將時域范圍內(nèi)的偏微分方程轉(zhuǎn)化為頻域范圍內(nèi)的常微分方程，將彈性波的波動方程的解表達成簡諧波的形式，根據(jù)在各個局域坐標(biāo)內(nèi)的出入關(guān)系，將簡諧波設(shè)置成為入射波和出射波，將入射波的波幅和出射波的波幅作為未知量進行表達[8]。根據(jù)各自單元內(nèi)物理量之間的關(guān)系得出相應(yīng)的相位關(guān)系，根據(jù)節(jié)點處的力平衡方程和位移協(xié)調(diào)方程得出相應(yīng)的散射關(guān)系，將各自節(jié)點的散射關(guān)系以及相位關(guān)系按照一定的關(guān)系進行組合排列，形成整體相位矩陣與散射矩陣，進而得到整體結(jié)構(gòu)的回傳矩陣，計算穩(wěn)態(tài)響應(yīng)，得出對應(yīng)的頻率響應(yīng)函數(shù)曲線。

3 理論分析

3.1 單元劃分與對偶變換

回傳射線矩陣法是基于單元劃分的思想，又因為一維周期結(jié)構(gòu)是在x軸方向上周期排列的復(fù)合結(jié)構(gòu)，二者共同擁有的單元劃分思想為回傳射線矩陣法的計算提供了較為便利的基礎(chǔ)。

采用的局部坐標(biāo)系是根據(jù)右手螺旋定則進行設(shè)定[9]。對于單元ij，引入兩組對偶坐標(biāo)系，如圖1所示。分別定義為（x，y，z）ij和（x，y，z）ji，其中，（x，y，z）ij表示由由節(jié)點i沿著梁單元的軸線指向節(jié)點j的局部坐標(biāo)，對應(yīng)的，xij為由節(jié)點i到節(jié)點j的方向為正方向，則xji為節(jié)點j到節(jié)點i的方向為正方向。通過右手螺旋法則，可以得到相對應(yīng)的yij和yji，zij和zji，可以看出yij和yji反向，zij和zji同向。同時定義，在對偶坐標(biāo)系之中，xij=lij-xji。

圖1 局部坐標(biāo)定義Fig.1 Local Coordinate Systems for Beam Members

在上述定義的對偶坐標(biāo)之中，可以確定在每個對偶坐標(biāo)之中，軸向力N，剪切力Q，彎矩M等力向量以及軸向位移u，縱向位移v以及轉(zhuǎn)角位移φ之間的關(guān)系，如式（1）所示。

對于單元ij而言，可以定義廣義力向量Fij和廣義位移向量Vij，二者的具體形式為：

因回傳射線矩陣法的推導(dǎo)列式以及算法思想是在局部坐標(biāo)系之中進行建立的，需要將局部坐標(biāo)系之中的物理量轉(zhuǎn)化為整體坐標(biāo)系之中，因此需要進行對偶變換。對于單元ij內(nèi)的廣義力向量和廣義位移向量而言，其對偶變換關(guān)系為：

式中：θ—局部坐標(biāo)系中的x軸與整體坐標(biāo)系中x軸的夾角；Fi

Oj，Vi

O

j

—整體坐標(biāo)系下的廣義力向量和廣義位移向量。

3.2 相位關(guān)系矩陣

為了將波動方程解寫成簡諧波的形式，將簡諧波定義為入射波和出射波，以單元12為例，入射波向量定義為：

對應(yīng)的出射波向量定義為：

根據(jù)單元內(nèi)的對偶關(guān)系，結(jié)合廣義力向量與廣義位移向量的簡諧波形式，可以得出單元范圍內(nèi)的相位關(guān)系：

可以得到單一單元范圍內(nèi)的相位關(guān)系：

根據(jù)結(jié)構(gòu)的節(jié)點定義，按照一定的排列順序可以將單元相位關(guān)系矩陣組裝成總體相位關(guān)系矩陣：

其中，P為6n×6n階對角陣。通過節(jié)點的定義與單元的劃分關(guān)系，可以得出總體結(jié)構(gòu)的相位關(guān)系：

3.3 散射關(guān)系矩陣

回傳射線矩陣法的散射關(guān)系矩陣是基于節(jié)點處的力平衡關(guān)系以及位移協(xié)調(diào)關(guān)系進行計算。假設(shè)節(jié)點i處有nj個單元，則作用在節(jié)點i出的外力與單元內(nèi)力之間的力平衡關(guān)系表示如下：

式中：Tij—上述描述的轉(zhuǎn)換矩陣；fi—作用在節(jié)點i處的外力源向量；Mi，Ki—節(jié)點i處的集中質(zhì)量矩陣以及彈簧剛度矩陣。

結(jié)合廣義力向量的簡諧波形式，可以將上述力平衡公式寫成：

式中：Ai

Fj

O（x，ω），Di

Fj

O（x，ω）—6×6階廣義力傳播矩陣，可以根據(jù)單元控制微分方程推到得出。

根據(jù)節(jié)點的位移協(xié)調(diào)條件，可以得出：

結(jié)合廣義位移向量的簡諧波形式，可以得出：

式中：Ai

Vj

O（0，ω），Di

Vj

O（0，ω）—6×6階廣義位移傳播矩陣，可以對單元的控制微分方程進行推導(dǎo)得出。

根據(jù)上述力平衡關(guān)系矩陣以及位移協(xié)調(diào)關(guān)系矩陣，聯(lián)立式（10）與式（11）消去位移向量ui可以得出節(jié)點i的散射關(guān)系：

式中：Si—局部關(guān)系散射矩陣；si—節(jié)點處的外力源向量，根據(jù)節(jié)點所受外力以及按照整個結(jié)構(gòu)的節(jié)點排列順序，可以得出整體結(jié)構(gòu)的散射關(guān)系矩陣S：

式中：S—6n×6n階對角陣。通過節(jié)點與單元的劃分定義，按照節(jié)點的排列順序，則可以得到總體散射關(guān)系：

3.4 回傳矩陣與穩(wěn)態(tài)響應(yīng)

為了將總體結(jié)構(gòu)的關(guān)系矩陣進行組裝，需要引入相位關(guān)系轉(zhuǎn)換矩陣。我們可以了解到，出入波向量d和之間的具有相同的元素，但是二者之間各元素的排列順序不同，因此需要引入轉(zhuǎn)換矩陣U，可以得出：則有a=PUD

將上述總體相位關(guān)系矩陣以及散射關(guān)系矩陣進行聯(lián)立，可以得出：

其中，定義R=SPU為系統(tǒng)的回傳射線矩陣。

將相應(yīng)的動力響應(yīng)寫成穩(wěn)態(tài)響應(yīng)的模式：

將相應(yīng)的頻率代入廣義力與廣義位移的表達式，可以得出時域范圍內(nèi)的額廣義力與廣義位移向量，如下所示，進而可以求出頻率范圍內(nèi)的穩(wěn)態(tài)響應(yīng)，得出對應(yīng)的頻率響應(yīng)函數(shù)曲線[6]。

經(jīng)過數(shù)據(jù)處理得出相應(yīng)的頻率響應(yīng)函數(shù)曲線，研究分析其振動特性。

4 數(shù)值算例分析

本節(jié)基于Euler-Bernoulli梁理論，以一維等截面周期結(jié)構(gòu)結(jié)構(gòu)的穩(wěn)態(tài)響應(yīng)最為算例進行計算，本節(jié)采取三周期的一維等截面周期結(jié)構(gòu)模型，如圖2所示。周期結(jié)構(gòu)晶格，如圖3所示。具體參數(shù)如下：鋁：密度ρ1=2799kg/m3，泊松比υ1=0.33，彈性模量E1=7.21×1010。有機玻璃：密度ρ2=1142kg/m3，泊松比υ2=0.36，彈性模量E1=0.32×1010。結(jié)構(gòu)物理參數(shù)：l1=l2=50mm，橫截面為圓形r1=r2=2mm定義晶格常數(shù)N=l1+l2=100mm，組分u=l1/（l1+l2）。

圖2 周期結(jié)構(gòu)模型Fig.2 Structure of Periodic Structures Beam

圖3 晶格簡圖Fig.3 The Schematic of the Unit Cell

基于回傳射線矩陣法的單元劃分思想，將一維等截面周期結(jié)構(gòu)模型按照節(jié)點與單元的思想進行劃分，得到的節(jié)點定義與單元劃分，如圖4所示。

圖4 單元與節(jié)點定義Fig.4 The Definition of Unit and Node

在每一段梁單元之中建立對偶坐標(biāo)系，以梁單元12以及單元23之中為例進行建立，得到的對偶坐標(biāo)系，如圖5所示。

圖5 對偶坐標(biāo)建立Fig.5 The Establish of Dual Coordinate

基于Euler-Bernoulli梁理論，推導(dǎo)單元的控制微分方程，得出控制微分方程的解，并且寫成簡諧波的形式，得出對應(yīng)的廣義力向量與廣義位移向量：

通過前文所述的結(jié)構(gòu)的相位關(guān)系，將每個梁單元之中的物理量進行相應(yīng)的對偶變換，可以得出單元的相位關(guān)系矩陣，按照節(jié)點的定義順序，將節(jié)點1至節(jié)點6的單元相位關(guān)系矩陣進行組裝，得出總體結(jié)構(gòu)的相位關(guān)系矩陣，其中P1-6為6×6階矩陣

通過各個節(jié)點處的力平衡關(guān)系以及位移協(xié)調(diào)關(guān)系，得到整體結(jié)構(gòu)中各個節(jié)點的散射關(guān)系矩陣，其中，節(jié)點0與節(jié)點7的散射矩陣為3×3階，其余散射矩陣為6×6階矩陣。

通過引入轉(zhuǎn)換關(guān)系矩陣U，轉(zhuǎn)換矩陣U的特點為：其每一行以及每一列有且僅有一個元素為1，其余元素為零，并且其逆矩陣為其自身。

對于整個結(jié)構(gòu)，通過引入邊界條件，得出總體結(jié)構(gòu)對應(yīng)的回傳射線矩陣，求解相應(yīng)的穩(wěn)態(tài)響應(yīng)，得出對應(yīng)頻率范圍內(nèi)的振動特性曲線。

本節(jié)算例之中，為了減小與理想周期結(jié)構(gòu)結(jié)構(gòu)之間的差距，采取自由邊界條件進行計算[10]，取一維等截面周期結(jié)構(gòu)的左端面為激勵輸入端，取右端面為激勵響應(yīng)拾取端；本節(jié)算例施加簡諧位移軸向激勵，其振幅為1mm。

為了驗證本節(jié)算例計算一維等截面周期結(jié)構(gòu)對軸向波的頻率響應(yīng)函數(shù)曲線的正確性，因此采取有限元分析法對上述算例進行有限元分析計算[11]。

采用AnsysWorkbench進行有限元分析，具體設(shè)置參數(shù)如下：

材料參數(shù)與周期結(jié)構(gòu)物理參數(shù)：與前文結(jié)構(gòu)描述相同。

網(wǎng)格劃分：因模型較為基本，且較為簡單，因此采用自由網(wǎng)格劃分的方法，定義六面體單元的基本尺寸為2mm。

施加載荷：施加軸向簡諧位移載荷，位移的幅值為1mm，激振頻率范圍為（0～40）kHz。

分析參數(shù)設(shè)置：因施加的載荷為簡諧位移載荷，采用諧響應(yīng)分析模塊，采用完全算法，即：Full算法。

圖6 有限元分析示意圖Fig.6 The Schematic of Finite Element Analysis

為了比較分析研究本節(jié)算例的數(shù)值計算結(jié)果以及有限元分析方法之間的關(guān)系，由圖7可以看出，由有限元分析方法計算得出的曲線與由回傳射線矩陣法計算得出的頻率響應(yīng)函數(shù)曲線基本保持吻合，因此可以證明回傳射線矩陣法計算一維等截面周期結(jié)構(gòu)對于軸向波的頻率響應(yīng)函數(shù)曲線的正確性。

圖7 有限元與數(shù)值計算對比Fig7 The Comparison of FEM and Numerical

5 結(jié)論

依據(jù)回傳射線矩陣法的理論思想，依據(jù)其單元劃分以及節(jié)點定義思想，結(jié)合單元的相位關(guān)系矩陣以及節(jié)點的散射關(guān)系矩陣，通過引入轉(zhuǎn)換矩陣以及總體關(guān)系矩陣的組裝，結(jié)合邊界條件，得出總體結(jié)構(gòu)的回傳射線矩陣，求得結(jié)構(gòu)中的穩(wěn)態(tài)響應(yīng)，得出頻率范圍內(nèi)的振動特性曲線。以一維等截面周期結(jié)構(gòu)結(jié)構(gòu)對與軸向波的傳輸特性曲線進行了計算，并且與有限元分析方法進行了對比，證明回傳射線矩陣法對周期結(jié)構(gòu)結(jié)構(gòu)振動特性計算的正確性。