采用譜單元Galerkin法求解非線性模態(tài)

李 誠，李鴻光

（上海交通大學 機械系統(tǒng)與振動國家重點實驗室，上海 200240）

非線性模態(tài)理論作為非線性振動研究的一個重要分支，是線性模態(tài)理論在某些非線性振動系統(tǒng)中的一種對應拓展。自Ronsenberg 引入非線性模態(tài)NNM(Nonlinear Normal Modes)的概念以來，國內外不少學者在非線性模態(tài)的定義、動力學特性以及非線性模態(tài)的求解方法等方面都做出了重要貢獻。Ronsenberg首先將非線性自治保守系統(tǒng)中各個自由度的同步周期運動模式定義為非線性模態(tài)[1]。Shaw和Pierre 引入不變流形的概念定義非線性模態(tài)。不變流形定理表明，非線性系統(tǒng)的不變流形與派生線性系統(tǒng)的不變子空間相切于系統(tǒng)平衡點，該不變子空間為派生線性系統(tǒng)的線性模態(tài)。系統(tǒng)發(fā)生在不變流形上的運動即為非線性模態(tài)運動[2]。吳志強和陳予恕發(fā)展了該思想[3]，定義非線性模態(tài)為系統(tǒng)模態(tài)空間中偶數維不變流形上的運動，通過求解系統(tǒng)動力學方程的規(guī)范型(Normal Form)來構造非線性模態(tài)的控制方程，根據方程將非線性模態(tài)分為非耦合模態(tài)和內共振模態(tài)。在該定義的基礎上，李欣業(yè)等[4]對多自由度內共振非線性系統(tǒng)中內共振模態(tài)的分岔特性進行了研究。徐鑒等[5]研究了系統(tǒng)在非奇異條件下非線性模態(tài)疊加解的有效性與模態(tài)動力學方程靜態(tài)分岔的關系。Cirillo 等[6]發(fā)現了非線性模態(tài)與Koopman 算子譜特性之間的關系。非線性模態(tài)的求解方法包括Nayfeh 等采用的多尺度法[7]、Vakakis等[8]提出的一種基于能量的冪級數展開漸近解法，文獻[9－10]中采用的規(guī)范型方法等。Shaw和Pesheck等應用非線性模態(tài)的不變流形的定義，采用了相應的基于冪級數展開的漸近解法[11]，更進一步提出了基于Galerkin法的半解析方法[12]。Galekin法的求解精度一般均優(yōu)于基于攝動法的漸近展開解，尤其在非線性系統(tǒng)響應幅值較大時，Galerkin法在精度上更具有顯著優(yōu)勢。對于保守系統(tǒng)非線性模態(tài)所對應的一系列周期運動可采用打靶法等相關改進算法[13－15]，結合延拓連續(xù)法獲得一系列數值解。

非線性模態(tài)的曲面幾何結構可用來描述非線性系統(tǒng)發(fā)生模態(tài)運動時各自由度之間的相對關系[16]。對于動力系統(tǒng)的非線性模態(tài)，更高的求解精度可以更準確地揭示系統(tǒng)在非線性模態(tài)運動時各部分之間的耦合動力學特性。本文提出了一種基于譜單元的Galerkin 法以求解不變流形定義下的非線性模態(tài)曲面。以一個兩自由度非線性系統(tǒng)為例，指出在求解非線性模態(tài)時已有的Galerkin分片方法所求解區(qū)域較大處解的精度仍可能存在不足。基于譜單元的Galerkin 法可獲得整體上更高精度的非線性模態(tài)曲面解。

1 不變流形定義下的非線性模態(tài)微分方程

考慮保守非線性系統(tǒng)在其廣義模態(tài)坐標下的控制方程：

方程組中ui(i=1,2,…,n)是系統(tǒng)投影到派生線性系統(tǒng)模態(tài)空間上的廣義模態(tài)坐標，ωi即為系統(tǒng)第i階線性模態(tài)的固有頻率，設作用在該階模態(tài)振子上的回復力fi為廣義模態(tài)坐標位移向量ui的非線性光滑函數。假設各模態(tài)振子之間未出現內共振，根據非線性模態(tài)所采用的不變流形的定義，系統(tǒng)第M階非線性模態(tài)的運動中各自由度之間滿足如下的約束關系：

第M階模態(tài)坐標uM、為該階非線性模態(tài)的“主自由度”，其余坐標ui、(i≠M)為對應的“從自由度”，各個從自由度和主自由度之間在非線性模態(tài)運動中滿足式(2)的約束關系。文獻[12]對主坐標對引入一種極坐標變換，將主坐標對(uM,) 轉換為(a,φ)：

主坐標對的控制方程化為：

待求的非線性模態(tài)的約束關系式(2)變?yōu)椋?/p>

之所以先引入該極坐標變換，是因為非線性模態(tài)的不變流形曲面一般是關于主坐標中的相角φ在[0,2π]上的周期連續(xù)函數，其半解析解便可假設為關于φ的諧波函數的展開式。將式(5)分別代入方程組(1)中各個從自由度振子的微分方程對應的1階微分形式，應用鏈式法則得：

將式(4)代入式(6)，等式兩邊同乘以a以避免解在a=0處的奇異性，最終得到不變流形定義下第M階非線性模態(tài)的控制方程：

2 求解非線性模態(tài)的Galerkin法

針對非線性模態(tài)不變流形的微分方程式(7)，假設解為展開式(8)：

Pi(a,φ)、Qi(a,φ)為關于主坐標對a和φ的未知函數，展開式中各項基函數的待求系數為C和D，關于a和φ的基函數的個數分別為Na和Nφ。在求解域{(a,φ)|a∈[0,a0],φ∈[0,2π]}上應用Galerkin 法，將式(8)作為假設解代入后將該余量式選取為T(a,φ)，U(a,φ)作為試函數，令積分式為0，得式(9)、式(10)：

該方程組經過加權積分最終轉化為各展開項待求系數C和D的非線性代數方程組，可應用各種數值求根算法如Powell混合算法[17]等求解。而展開式系數代數方程組的形式由對應的非線性模態(tài)假設解的形式決定，并將影響后續(xù)方程組數值迭代求根過程中Jacobian矩陣的形式。

因為不變流形定義下非線性模態(tài)的約束曲面是關于主坐標相角φ在[0,2π] 上的周期連續(xù)函數，所以對于T(a,φ)、U(a,φ)關于φ的展開函數可選取一系列諧波函數作為基函數：

因為在保守系統(tǒng)非線性模態(tài)的周期運動中，各自由度的響應均保持同步，所以式(11)中選取余弦函數作為位移響應的展開式基函數，選取正弦函數作為速度響應的展開式基函數[12]。將式(11)代入式(9)、式(10)，得到如下代數方程組：

式中：p=1,2,…,Na,q=1,2,…,Nφ;i,j=1,2,…,n,i,j≠M。關于系數C和D的代數方程gp,qi=0和hp,qi=0 分別來自于積分式(9)和式(10)。因方程組(1)中的各非線性剛度項fM和fi(i≠M)耦合了所求的其他位移響應ui，所以式(12)、式(13)中的各代數方程可能含有相應的待求系數C。

2.1 Galerkin分片求解法

根據上述的求解策略，文獻[12]提出了一種Galerkin 分片求解法，將不變流形曲面關于主坐標(a,φ)的求解區(qū)域沿a方向離散，劃分為若干個環(huán)狀域，分別獨立求解各個環(huán)狀域內的不變流形子曲面。每個環(huán)狀求解域內，將函數L(a) 假設為線性函數。分別在兩個相鄰的環(huán)域內進行求解，所得的在相接處兩個子曲面的函數值一般不同，整體的不變流形曲面會出現階躍不連續(xù)，一般可取相接處的節(jié)點在各自求解域中的函數值的均值作為該節(jié)點的函數值，以保證解曲面的連續(xù)。如圖1所示，圖1(a)和圖1(b)分別為某一不變流形曲面uSub=P(uM,)由原點向外第1、第3、第5個和第2、第4個環(huán)狀域上求解出的子曲面，各子曲面均在各自子域內獨立求解。圖1(c)給出了整體不變流形曲面的側視圖，可以看出在相鄰曲面的相接處均存在曲面階躍不連續(xù)。圖1(d)顯示了在各個相鄰子曲面相接處函數值取均值后的整體非線性模態(tài)曲面，子曲面相接處函數值連續(xù)。

圖1 采用分片Galerkin法求解出的非線性模態(tài)曲面

若采用更多更細的環(huán)狀子域分別求解非線性模態(tài)曲面，相鄰子曲面相接處的階躍不連續(xù)將減小，但合成后的解曲面的精度仍存在不足，在后續(xù)算例中將給出說明。

2.2 基于譜單元的Galerkin法

在幅值a維度上采用譜單元進行離散近似。應用Patera提出的譜元法[18]，在將求解域沿a方向劃分為相接的各個單元后，在每個單元內采用Lagrange插值函數作為基函數，插值點選擇為第二類Chebyshev 多項式在所在單元定義域上對應的各個零點。表1給出了定義域為[-1，1]上第m階的第二類Chebyshev多項式的零點。

表1 第二類Chebyshev多項式的零點

當L(a) 采用一系列譜單元進行插值近似后，不同于分片法，相鄰單元的共用節(jié)點處的函數值應相等。圖2給出了譜單元上各插值基函數的示意圖，設第r個單元具有4 個插值點，第r+1 個單元具有5個插值點，將圖中各插值函數結合插值節(jié)點的函數值以近似某解曲線。

圖2 譜單元上的基函數示意圖

將[a0,aend] 劃分為N個單元，設第r個單元含有包括單元兩個端點在內的共mr個節(jié)點相鄰單元r與單元r+1 的共同節(jié)點坐標求解域內共有個插值節(jié)點，其中單元內的節(jié)點個，單元兩端的節(jié)點共N+1 個。第r個單元內、除兩個端點以外的插值節(jié)點函數值對應的基函數分別為：

當r=1,…,N,s=2,…,mr-1

在各個單元的端點，如第r個端點(r=1,…,N+1)，節(jié)點函數值對應的基函數分別為：

當r=1時，

當r=2,…,N時，

當r=N+1時，

將式(14)至式(17)代入式(11)后再代入式(9)、式(10)，采用Galerkin 法求解。不同于分片Galerkin法，其積分方程的積分域為a方向上兩個節(jié)點之間的環(huán)域，而此方法中每個方程的積分域包含了對應節(jié)點的相鄰環(huán)域，最終得到的代數方程同時考慮了相關鄰域內控制方程式(9)式(10)對非線性模態(tài)曲面的約束。由于在節(jié)點處函數具有連續(xù)性，在各個子域上經Galerkin積分獲得的有關節(jié)點函數值的代數方程之間相互耦合，需聯立求解。相比于分片Galerkin 法中求解一系列獨立的代數方程組再對節(jié)點處函數值取平均的策略，此方法得到的節(jié)點值從整體上考慮了每個節(jié)點相鄰鄰域的影響，有可能獲得更精確的解，這將在算例研究中進行驗證。

后續(xù)的非線性代數方程組求根采用基于梯度的迭代數值算法，求解最小二乘意義下方程組的余量函數的最小值。迭代過程中需要計算待求函數的Jacobian 矩陣。采用本方法求解非線性模態(tài)曲面展開系數時，其相關的Jacobian 矩陣含有式(18)中的4個分塊矩陣。進行數值迭代求根時，矩陣中的函數偏導數采用差分法進行近似。

其中：p,k=1,2,…,Na,q,l=1,2,…,Nφ;i,j=1,2,…,n,i,j≠M數值迭代求根中的Jacobian 矩陣將由于在譜單元中選取不同階數的插值多項式而呈現出其對應特定的稀疏性。設從坐標個數j=1，諧波基函數個數Nφ=2，a方向劃分N=2個單元，第一個單元采用3個節(jié)點的Lagrange插值多項式，第二個單元采用4個節(jié)點的插值多項式。兩個單元內的插值節(jié)點分別對應第二類Chebyshev多項式的第3階、第4階多項式的根。以此為例，圖3給出了系數C和D的代數方程組在迭代求根中的Jacobian矩陣結構示意圖。

圖3中具有陰影底紋的網格表示矩陣中的非零元素。在系數方程組的數值求根迭代過程中，為了簡化Jacobian 矩陣的計算，根據其特定的稀疏性只計算矩陣的非零元素。

圖3 稀疏Jacobian矩陣示意圖

3 兩種Galerkin方法求解非線性模態(tài)的比較

本節(jié)以一個具有強二次及三次非線性的兩自由度振動系統(tǒng)為例，采用上述的兩種方法計算系統(tǒng)所包含的各階非線性模態(tài)，并比較其精確性。

3.1 動力學模型

以文獻[15]給出的非線性質量彈簧振動系統(tǒng)為研究對象，如圖4所示。集中質量mo可以在u1、u2方向上運動，由兩個原長均為L的彈簧連接。沿用文獻中采用的Green-Lagrange應變來定義大變形下彈簧的應變：由此推導出的系統(tǒng)動力學方程式(19)中含有二次及三次的光滑非線性剛度項。

圖4 兩自由度非線性振動系統(tǒng)示意圖

設mo=1kg，L=1m，得到系統(tǒng)的動力學方程。采用文獻[15]中的參數設置，考慮k1=1N、k2=2N時的情況。

3.2 結果比較

采用以上兩類方法求解振動系統(tǒng)(19)的各階非線性模態(tài)。首先應用分片Galerkin 法，將求解域劃分為12 個獨立子環(huán)域，令各子環(huán)域a方向的寬度相等，取諧波基函數的個數為6。再應用譜單元Galerkin法求解非線性模態(tài)，分別采用了4種網格劃分及譜單元階數的組合，如表2所示。

表2 應用譜單元Galerkin法時選取的求解參數

各種組合策略中，設置相鄰插值節(jié)點間的距離與分片Galerkin 法中獨立子環(huán)域片數為12 時環(huán)域a方向的寬度接近，以比較兩類方法在近似的離散程度下所求得的解的精度差異。本例中各種譜單元的組合策略中均令單元的長度相等，取諧波基函數的個數為6。應用兩類方法時對所生成的非線性代數方程組進行數值求解，待求系數的迭代初值均取為0，精度取為1×10-6。圖5給出了運用兩類Galerkin方法求得的第1階非線性模態(tài)的近似曲面。

圖5(a)同時給出了發(fā)生在系統(tǒng)該階非線性模態(tài)曲面上的一系列周期運動響應，即圖中實線所示的一組封閉相軌線。該組周期軌線均通過打靶法進行數值求解得到。取某一初值為(u10,u˙10,u20,u˙20)=(1.7,0,0.009 617,0)的周期響應，代入原方程式(13)，通過Runge-Kutta法進行數值積分，可得參考周期解的時域響應u1(t)和u2(t)。再取另一組初值可得周期解時域響應u1(t)和u2(t)，見圖5(b)，以作參考比較。

當采用兩類方法分別求得非線性模態(tài)曲面近似展開式中的待求系數后，即得到了式(2)形式的非線性模態(tài)曲面顯式表達式。將其代入方程組(1)中主坐標uM所對應的控制方程中的非線性回復力項，可得關于uM的單自由度的微分方程，方程中只含有uM一個未知函數。取初值對該方程采用Runge-Kutta 法進行數值積分，得到響應解uM(t)。再將該響應序列uM(t)代入式(2)，可得對應發(fā)生在非線性模態(tài)不變流形曲面上的從坐標響應序列ui(t)。uM(t)和ui(t)描述了非線性模態(tài)近似解曲面上的一個響應軌線，因為其是根據近似展開解(2)得出，所以可通過它與對應的參考周期解比較來判斷非線性模態(tài)近似解的精確性。對于第1階非線性模態(tài)，uM(t)=u1(t)，ui(t)=u2(t)，圖5(b)中分別給出了初值為和=(1.9,0)時依據各近似解曲面得出的時域響應。

由圖5可以看出，采用各種譜單元的Galerkin法求得的非線性模態(tài)曲面，相比于分片Galerkin 法的結果，均具有較高的精度。即使計算幅值較大的非線性模態(tài)曲面上的周期軌，根據由譜單元Galerkin法獲得的近似解曲面，進而求得的曲面上的周期響應也與參考解較為接近。

圖6給出了采用兩類方法求解的第2 階非線性模態(tài)的各個解曲面及相應的一組周期軌線，同時給出了初值為= (1.8,0)和=(2.0,0)時依據各近似解曲面得到的時域響應及對應的參考解。采用兩類Galerkin方法時選取的參數與求解第1階非線性模態(tài)時選取的相應參數基本相同。系統(tǒng)(19)的第2 階非線性模態(tài)中，u2、為主自由度，u1為從自由度。由圖6可以看出，采用一系列譜單元Galerkin 法仍可獲得較為準確的第2 階非線性模態(tài)解曲面及其曲面上發(fā)生的周期響應。圖5(a)和圖6(a)同時給出了采用分片Galerkin 法在獨立子環(huán)域片數為25 時的解。可以看出即使采用了更細的網格，分片Galerkin 法在響應較大區(qū)域依然無法獲得較為精確的解。

圖5 第1階非線性模態(tài)解曲面及其周期響應

圖6 第2階非線性模態(tài)解曲面及其周期響應

4 結語

本文采用一種基于譜單元的Galerkin 法，可以較為準確地計算一類無阻尼非線性振動系統(tǒng)在不變流形定義下的非線性模態(tài)。相比于已有的非線性模態(tài)Galerkin 分片求解法，該方法選取第二類Chebyshev多項式的零點構造每個譜單元的Lagrange插值函數，對求解域進行整體求解。在展開系數的迭代求解中，Jacobian由于譜單元中選取不同階數的多項式而呈現特定的稀疏性。相較于分片求解法，該方法在求解域較大時仍能獲得較為準確的解。