李 誠(chéng)，李鴻光

（上海交通大學(xué) 機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240）

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

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

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

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

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

第M階模態(tài)坐標(biāo)uM、為該階非線性模態(tài)的“主自由度”，其余坐標(biāo)ui、(i≠M(fèi))為對(duì)應(yīng)的“從自由度”，各個(gè)從自由度和主自由度之間在非線性模態(tài)運(yùn)動(dòng)中滿足式(2)的約束關(guān)系。文獻(xiàn)[12]對(duì)主坐標(biāo)對(duì)引入一種極坐標(biāo)變換，將主坐標(biāo)對(duì)(uM,) 轉(zhuǎn)換為(a,φ)：

主坐標(biāo)對(duì)的控制方程化為：

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

之所以先引入該極坐標(biāo)變換，是因?yàn)榉蔷€性模態(tài)的不變流形曲面一般是關(guān)于主坐標(biāo)中的相角φ在[0,2π]上的周期連續(xù)函數(shù)，其半解析解便可假設(shè)為關(guān)于φ的諧波函數(shù)的展開式。將式(5)分別代入方程組(1)中各個(gè)從自由度振子的微分方程對(duì)應(yīng)的1階微分形式，應(yīng)用鏈?zhǔn)椒▌t得：

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

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

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

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

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

因?yàn)椴蛔兞餍味x下非線性模態(tài)的約束曲面是關(guān)于主坐標(biāo)相角φ在[0,2π] 上的周期連續(xù)函數(shù)，所以對(duì)于T(a,φ)、U(a,φ)關(guān)于φ的展開函數(shù)可選取一系列諧波函數(shù)作為基函數(shù)：

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

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

2.1 Galerkin分片求解法

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

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

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

2.2 基于譜單元的Galerkin法

在幅值a維度上采用譜單元進(jìn)行離散近似。應(yīng)用Patera提出的譜元法[18]，在將求解域沿a方向劃分為相接的各個(gè)單元后，在每個(gè)單元內(nèi)采用Lagrange插值函數(shù)作為基函數(shù)，插值點(diǎn)選擇為第二類Chebyshev 多項(xiàng)式在所在單元定義域上對(duì)應(yīng)的各個(gè)零點(diǎn)。表1給出了定義域?yàn)閇-1，1]上第m階的第二類Chebyshev多項(xiàng)式的零點(diǎn)。

表1 第二類Chebyshev多項(xiàng)式的零點(diǎn)

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

圖2 譜單元上的基函數(shù)示意圖

將[a0,aend] 劃分為N個(gè)單元，設(shè)第r個(gè)單元含有包括單元兩個(gè)端點(diǎn)在內(nèi)的共mr個(gè)節(jié)點(diǎn)相鄰單元r與單元r+1 的共同節(jié)點(diǎn)坐標(biāo)求解域內(nèi)共有個(gè)插值節(jié)點(diǎn)，其中單元內(nèi)的節(jié)點(diǎn)個(gè)，單元兩端的節(jié)點(diǎn)共N+1 個(gè)。第r個(gè)單元內(nèi)、除兩個(gè)端點(diǎn)以外的插值節(jié)點(diǎn)函數(shù)值對(duì)應(yīng)的基函數(shù)分別為：

當(dāng)r=1,…,N,s=2,…,mr-1

在各個(gè)單元的端點(diǎn)，如第r個(gè)端點(diǎn)(r=1,…,N+1)，節(jié)點(diǎn)函數(shù)值對(duì)應(yīng)的基函數(shù)分別為：

當(dāng)r=1時(shí)，

當(dāng)r=2,…,N時(shí)，

當(dāng)r=N+1時(shí)，

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

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

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

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

圖3 稀疏Jacobian矩陣示意圖

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

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

3.1 動(dòng)力學(xué)模型

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

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

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

3.2 結(jié)果比較

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

表2 應(yīng)用譜單元Galerkin法時(shí)選取的求解參數(shù)

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

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

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

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

圖6給出了采用兩類方法求解的第2 階非線性模態(tài)的各個(gè)解曲面及相應(yīng)的一組周期軌線，同時(shí)給出了初值為= (1.8,0)和=(2.0,0)時(shí)依據(jù)各近似解曲面得到的時(shí)域響應(yīng)及對(duì)應(yīng)的參考解。采用兩類Galerkin方法時(shí)選取的參數(shù)與求解第1階非線性模態(tài)時(shí)選取的相應(yīng)參數(shù)基本相同。系統(tǒng)(19)的第2 階非線性模態(tài)中，u2、為主自由度，u1為從自由度。由圖6可以看出，采用一系列譜單元Galerkin 法仍可獲得較為準(zhǔn)確的第2 階非線性模態(tài)解曲面及其曲面上發(fā)生的周期響應(yīng)。圖5(a)和圖6(a)同時(shí)給出了采用分片Galerkin 法在獨(dú)立子環(huán)域片數(shù)為25 時(shí)的解?？梢钥闯黾词共捎昧烁?xì)的網(wǎng)格，分片Galerkin 法在響應(yīng)較大區(qū)域依然無法獲得較為精確的解。

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

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

4 結(jié)語

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