一種求解非線性模態的改進Galerkin法

李 誠， 李鴻光

(1.上海衛星工程研究所，上海 201109；2.上海交通大學 機械系統與振動國家重點實驗室，上海 200240)

非線性模態理論自Ronsenberg提出開始，作為線性模態理論在包含某些剛度非線性(如幾何非線性)的振動系統中的一種自然拓展，一直以來受到國內外學者的廣泛關注。在該理論發展期間，關于非線性模態的定義、與其相關的非線性系統的動力學特性、以及非線性模態的數值計算及半解析求解等方面，不少學者都做出了重大的貢獻。Ronsenberg針對非線性自治保守系統，將系統各自由度之間的一種同步的周期運動模式首先定義為非線性模態[1]。該定義也適用于線性模態，所以可以認為非線性模態是線性模態的自然推廣。Shaw等[2]將系統相空間中二維不變流形上的運動與非線性模態運動相聯系，用不變流形的約束關系來描述非線性模態，且根據不變流形定理，不變流形與對應線性系統的不變子空間相切于平衡點，不變子空間即對應線性系統的線性模態。吳志強等[3]進一步推廣了該思想，將非線性模態定義為系統模態空間中偶數維不變流形上的運動，提出求解系統動力學方程的規范型來描述非線性模態上的動力學特性，將非線性模態進一步分為非耦合模態和內共振模態。在上述定義的基礎上，李欣業等[4]研究了內共振非線性系統的內共振模態的分岔特性。徐鑒等[5]研究了兩自由度非對稱系統在非奇異時的非線性模態，著重研究了疊加解的有效性與模態動力學方程靜態分岔間的關系。Cirillo等[6]揭示了非線性模態與Koopman算子譜特性之間的聯系。

關于非線性模態的計算，Nayfeh等[7-9]應用了多尺度法，King等[10]根據非線性模態的頻率與能量的關系，提出了一種基于能量的級數展開式求解方法，文獻[11-13]采用了規范型方法，求得級數展開形式的近似解。Shaw等[14-15]根據其關于非線性模態的不變流形的定義，在提出了相應級數展開式的漸近解法之后、發展出了基于Galerkin法的半解析方法[16]。一系列Galekin法的求解精度一般情況下均優于基于攝動法的漸近展開解，尤其在非線性系統響應幅值較大時，Galerkin法在精度上更具有顯著優勢。對于非線性模態所對應的一系列周期運動常采用打靶法等數值算法[17-18]，結合延拓連續法獲得數值解。

提高動力系統非線性模態的求解精度，其目的之一在于它可以描述系統接近非線性模態運動時各自由度之間的耦合動力學特性。保守系統非線性模態的一個重要特性，是當其對應的非保守系統發生主共振時，在相空間中該受迫振動的響應軌線發生在保守系統非線性模態的周圍。所以，非線性模態的曲面幾何結構，可用來描述非線性系統發生主共振時的動力學特性，以及系統各自由度之間的耦合關系[19-20]。本文以一種非線性雙級隔振器系統為例，說明在求解非線性模態時，經典的兩種Galerkin方法均存在一些限制和不足；提出了一種改進的Galerkin法，可在一定程度上彌補經典Galekin法的限制和不足，獲得整體更高精度的非線性模態曲面解。

1 非線性模態的控制方程

本文主要針對的是無阻尼的非線性自治系統，方程形如式(1)

(1)

式中：ηi(i=1,2,…,n)為該系統在平衡位置處的派生線性系統的廣義模態坐標；ωi為系統第i階線性模態的固有頻率；fi為作用在該階模態振子上的非線性回復力，假設其為廣義模態坐標位移向量的非線性函數。暫不考慮自由度振子之間出現內共振的情況，根據采用不變流形形式的非線性模態的定義，系統某一階非線性模態的運動中各自由度之間滿足如下的約束關系

(2)

(3)

主坐標的控制方程變為

(4)

待求的非線性模態約束關系也變為了

(5)

因該非線性模態的約束曲面一般是關于主坐標相角φ在[0,2π]上的周期連續函數，所以在后續的求解中其半解析解可以假設為關于φ的諧波函數疊加形式。將式(5)代入式(1)從自由度對應的微分方程，并采用鏈式法則，得

(6)

再將式(4)代入式(6)，兩邊同乘以a得到該階非線性模態的控制方程

(7)

2 Galerkin系列方法求解非線性模態

假設控制方程式(7)的解為展開式(8)的形式

(8)

式中：Pi(a,φ),Qi(a,φ)分別為關于a和φ的未知函數；C和D為展開式中各項的待求系數；Na和Nφ分別為展開式中關于a和φ的基函數的個數。在求解域{(a,φ)|a∈[0,a0],φ∈[0,2π]上應用Galerkin法，依次選取T(a,φ),U(a,φ)為試函數，得到式(9)、式(10)

(9)

i=1,2,…,n;i≠M;p=1,2,…,Na;q=1,2,…,Nφ

(10)

該式一般為待求系數C和D的非線性代數方程組，可用Powell hybrid算法[21]等非線性代數方程求根數值算法求解。

非線性模態假設解的展開式形式將決定其對應的展開式系數代數方程組的形式。因假設非線性模態的約束曲面是關于主坐標相角φ在[0,2π]上的周期連續函數，所以T(a,φ),U(a,φ)中關于φ的展開函數可選取一系列三角函數

Tk,l(a,φ)=Lk(a)cos[(l-1)φ],Uk,l(a,φ)=Lk(a)sin(lφ)

(11)

對于保守系統，因為在其非線性模態的周期運動中，各自由度的響應均保持同步，所以式(11)中關于位移響應的展開式選取余弦函數作為基函數，關于速度響應的展開式選取正弦函數作為基函數[16]。關于式(11)中的Lk(a)，一種方法(方法1)可選擇一組定義在[0,a0]上的多項式，滿足如下定義的正交性，可使積分式(9)、式(10)中相應的乘積項為0，簡化最終的代數方程組

(12)

同時，對于線性剛度項已完全解耦的式(1)，因其非線性模態曲面與派生線性系統的不變子空間相切于系統的平衡點，這里的不變子空間即為(a,φ)平面，所以上述的這組多項式的一次導數需在a=0處為0。圖1展示了Pesheck等研究中給出的前7個多項式函數。

圖1 非線性模態解展開式采用的正交多項式Fig.1 Orthogonal polynomials assumed in nonlinear normal mode (NNM) expansion

在該方法中將式(11)代入式(9)、式(10)，得到的代數方程組為

(13)

(14)

為了減小待求系數非線性方程組的復雜度，便于快速求解，另一種策略(方法2)是將不變流形曲面關于主坐標(a,φ)的求解區域沿a方向離散，劃分為幾個環狀區域，再逐個獨立求解。每個環狀求解域內，函數L(a)假設為線性函數。兩個相鄰的環域在分別求解后，在相接區域曲面的函數值一般并不相等，形成的不變流形曲面會出現階躍不連續，在接下來的例子中將體現這一點。一般可取相接處的節點在各自求解域中的函數值的均值，作為該節點的函數值，以保證解曲面的連續。

本文在這兩種方法的基礎上，采用線性的有限元格式作為假設函數L(a)，式(15)～式(17)示意圖，如圖2所示。

圖2 采用有限元格式的分段線性函數L(a)Fig.2 Piecewise linear function as finite element approximation L(a)

(15)

(16)

(17)

將式(15)～式(17)代入式(11)后再代入式(10)，采用Galerkin法求解。不同于方法2中每個積分方程的積分域為a方向上兩個節點之間的環域，此方法(方法3)中每個方程的積分域為包含對應節點的相鄰環域，得到的代數方程即同時考慮了在相關鄰域內非線性模態的控制方程對待求曲面的約束要求。因為節點處函數的連續性，各個子域上經過Galerkin積分得到的關于節點函數值的代數方程之間互相耦合，需要聯立求解。相比于方法2中求解一系列獨立的代數方程組再取平均的策略，本方法得到的節點值整體考慮了每個節點相關鄰域的影響，有可能獲得更精確的解，這將在后續的例子中展示。

采用方法3獲得的代數方程組由式(18)～式(23)給出

當p=1;q=1,2,…,Nφ;i,j=1,2,…,n;i,j≠M時

(18)

(19)

當p=2,3,…,Na-1;q=1,2,…,Nφ;i,j=1,2,…,n;i,j≠M時

(20)

(21)

當p=Na,q=1,2,…,Nφ;i,j=1,2,…,n;i,j≠M時

(22)

(23)

相比于方法1中式(13)、式(14)里的每個方程均耦合了全部的未知數，本方法生成的方程組中每個代數方程只含有與相關節點相鄰的節點值作為未知數，每個方程組的耦合項大大減少。因為代數方程組的這個特點，其對應的Jacobian矩陣也具有一定的稀疏性。后續的非線性代數方程組求根采用基于梯度的迭代數值算法，求解最小二乘意義下方程組的余量函數的最小值，其Jacobian矩陣的函數偏導數采用差分法進行近似。稀疏Jacobian矩陣示意圖，如圖3所示。

圖3(d)、圖3(e)中深色的網格代表非零元素。為了簡化近似Jacobian矩陣的計算量，減少調用待求方程組的計算次數，本方法根據生成的代數方程組、即式(18)～式(23)所具有的稀疏性，設計了特定的Jacobian矩陣計算流程(見附錄A)。設方程組式(18)～式(23)中從坐標的方程個數為Nsub=n-1，每個從坐標的位移約束函數P(a,φ)和速度函數Q(a,φ)中，關于a和φ的基函數的個數分別為Na和Nφ，總的未知數個數Nx=2(n-1)NaNφ。若不考慮Jacobian矩陣的稀疏性，采用向前差分計算待求函數在當前迭代點的近似Jacobian矩陣時，總共需調用計算函數值的次數為Nf=1+Nx次，即當前迭代點處的每個分量都需要單獨計算其獲得微小增量后對應的函數值以進行差分。參考方法3中Jacobian稀疏矩陣的示意圖3，并采用附錄A中特定的稀疏Jacobian矩陣計算流程后，總共需調用計算函數值Nf=1+2(n-1)·3·Nφ次。

圖3 稀疏Jacobian矩陣示意圖Fig.3 Schematic diagramof sparse Jacobian matrix

3 Galerkin系列方法的算例結果比較

下文以一種非線性雙級隔振器為例，研究其對應的動力系統所包含的非線性模態的動力學特性。

3.1 雙級非線性隔振器動力學模型

圖4給出了該隔振器的構型示意圖，相較于文獻[22]中引入非線性彈性元件的雙級隔振器，此例中的隔振器考慮了另一種相關構型，其附加質量只與負載所在的隔振平臺通過彈性及阻尼元件連接。

圖4 非線性雙級隔振器示意圖Fig.4 Schematic diagram of a nonlinear two-stage isolator

式(24)給出了該隔振器負載所在的隔振平臺受到諧波激勵時的控制方程。

(24)

式中：kv1和kv2分別為兩個縱向彈性元件的剛度；lo1和lo2分別為側向彈性元件kh1和kh2的原長。當兩個平臺的位移響應滿足x1<0.4l1和x2<0.4l2[23]時，方程中的剛度非線性項可以用泰勒展開式近似，式(24)簡化為式(25)

(25)

將系統方程投影到線性模態坐標空間，考慮其對應的保守系統，得到形如式(26)的控制方程，

(26)

圖5(b)同時給出了這一組主共振響應所對應的無阻尼式(26)的自由響應，該組自由響應即是發生在系統非線性模態上的周期運動響應，即圖5(c)中的軌線，均由打靶法數值求解出。圖5(c)中的各條強迫響應曲線均在非線性模態周期軌線的周圍，且兩者變化趨勢一致，說明此處無阻尼系統的非線性模態可以描述原系統的強迫響應在達到主共振峰時各自由度之間的耦合關系。圖5(c)、圖5(d)分別給出了式(26)在原物理坐標空間和廣義模態坐標空間內的非線性模態曲面，發生在曲面上的周期運動軌線與圖5(b)中的自由響應相對應。

圖5 主共振受迫響應與非線性模態自由周期響應Fig.5 Forced responses in primary resonance and the corresponding free periodic motions on the NNM

3.2 非線性模態半解析解的結果比較

為求解與原式(25)的主共振相對應的保守式(26)的非線性模態，采用以上的3種方法求解。表1給出了3種方法在不同求解域上選取的參數。展開式中關于a的基函數的個數Na，在方法1中指的是多項式基函數的最高階數；在方法2中指的是劃分的獨立子域的數量，也是獨立的代數方程組的個數；在方法3中指的是沿a方向劃分的網格數。對于所有的非線性代數方程組的數值求解，待求系數的迭代初值均取0，精度取1×10-6。注意到在方法1中，當Na≥5時，某些參數條件下求根的數值迭代中止(表1中用“無”表示)，無法求得滿足方程的數值解。

表1 求解參數和計算時間 {(a,φ)|a∈[0,0.5],φ∈[0,2π]}

由圖6(a)中方法2計算的解曲面的側視圖可以看出，在兩個相鄰的求解域上分別單獨計算非線性模態曲面解，在相接處不連續，非線性模態曲面存在間斷。對節點處的函數值取均值后，解曲面如圖6(a)中所示。由圖6(e)可以看出，在a∈[0,0.5]內，3種方法求得的解曲面在P1(a,0)處的投影差別不大，均接近于各個精確解的投影點。

圖6 非線性模態解曲面及其投影 {(a,φ)|a∈[0,0.5],φ∈[0,2π]}Fig.6 Approximation surfaces and profiles for NNM in {(a,φ)|a∈[0,0.5],φ∈[0,2π]}

圖7、圖8分別給出了該非線性模態關于位移約束曲面在a∈[0,1.0]，a∈[0,2.0]內采用上述3種方法求得的解曲面。圖7(a)中的解分別采用了：方法1，選取Na=5，Nφ=10的展開式求解；方法2，劃分了25個獨立的環域分別求解再在節點處取平均；方法3，選取a方向的網格數為25，作為采用了方法3求得的一系列解的代表。圖8(a)中的解曲面分別采用了：方法1，選取Na=5，Nφ=10的展開式；方法2，劃分了50個環域分別求解再在相接節點處取平均處理；方法3獲得的解曲面a方向劃分的網格數為50。圖7(b)～圖7(e)、圖8(b)～圖8(e)分別給出了在a∈[0,1.0]，a∈[0,2.0]內采用3種方法獲得的各個解曲面在P1(a,0)處的截面，同時給出對應求解域上部分精確周期解在該截面的投影點。

由圖7(a)、圖7(b)及圖8(a)、圖8(b)可以看出，當積分求解域增大時，方法1的部分假設解在值較大處與精確解存在偏差。如表2、表3所示，若展開式中的多項式選取更高的階數，有可能在未知系數的初始值為0時方程組非線性迭代求根不收斂；還有可能，解曲面雖然能在整個求解域內滿足加權積分平均意義上的控制方程，但實際得到的解在a=1處或a=2處已經遠遠偏離精確解。

表2 求解參數和計算時間 {(a,φ)|a∈[0,1.0],φ∈[0,2π]}Tab.2 Solution parameters and computing time {(a,φ)|a∈[0,1.0],φ∈[0,2π]}

表3 求解參數和計算時間 {(a,φ)|a∈[0,2.0],φ∈[0,2π]}

由圖7、圖8可以看出，在a∈[0.4,1.8]附近，方法2求得的解曲面在P(a,0)處的投影與各個精確解的投影點存在一定的偏差，網格密度加密后，該偏差依然存在。采用方法3，取不同網格密度均能獲得較好的解，解曲面在P(a,0)處的投影相比于方法1和方法2的解，更為接近各精確解的投影點。

圖7 非線性模態解曲面及其投影 {(a,φ)|a∈[0,1.0],φ∈[0,2π]}Fig.7 Approximation surfaces and profiles for NNM in {(a,φ)|a∈[0,1.0],φ∈[0,2π]}

圖8 非線性模態解曲面及其投影 {(a,φ)|a∈[0,2.0],φ∈[0,2π]}Fig.8 Approximation surfaces and profiles for NNM in {(a,φ)|a∈[0,2.0],φ∈[0,2π]}

3.3 一種加速計算策略

當網格數較大時，因為方法2中每個獨立子域求解的代數方程組未知數個數都遠小于對應方法3中未知數的個數，采用方法2的所有求解域的整體計算時間仍遠小于方法3的計算時間。而在方法3的上述應用中，未知數的迭代初值均取為0。為了利用方法3獲得更精確的解曲面，同時減少計算成本，我們進而提出一種綜合了方法2和方法3的策略，即先采用方法2求解，獲得解曲面的一個初步的全貌，再將方法2得到的解作為方法3中未知數的初值，方法2中的網格數、諧波基函數保持不變，采用方法3獲得相較于方法2更精確的解。本例中，方法3采用了方法2的解作為初值后的計算時間在表4中給出，即使加上方法2的求解時間，整體計算時間也大大減少。該解與初值為0時求得的解基本相同，如圖9所示。

表4 求解參數和計算時間 {(a,φ)|a∈[0,2.0],φ∈[0,2π]}Tab.4 Solution parameters and computing time {(a,φ)|a∈[0,2.0],φ∈[0,2π]}

圖9 非線性模態解曲面及其投影 {(a,φ)|a∈[0,2.0],φ∈[0,2π]}Fig.9 Approximation surfaces and profiles for NNM in {(a,φ)|a∈[0,2.0],φ∈[0,2π]}

4 結 論

本文提出的一種改進的 Galerkin方法可以較為準確地計算一類無阻尼非線性系統的非線性模態在不變流形定義下的解曲面。相比于已有的第一種非線性模態Galerkin求解法，該方法在假設函數中引入了一次多項式的有限元格式的展開式，使得待求系數的非線性代數方程組中各方程的耦合項減少。在非線性方程組迭代求根的算法中，根據Jacobian矩陣在該問題中特定的稀疏性，本文采用了一種特定的Jacobian矩陣計算流程，可以有效減少待求函數的計算次數，從而加速求解。相比于已有的兩種非線性模態Galerkin解法，當該方法在求解域較大時，仍能獲得較為準確的解。當該方法與采用分段獨立求解策略的第二種Galerkin法相結合，將其解作為本方法迭代計算的初值，不僅可以保證解的精度，更能夠進一步減少計算時間。

附錄A

非線性模態求解中稀疏Jacobian矩陣的計算流程JACdf_Sparse，Fortran Free format格式。

“!”后為注釋。

SUBROUTINE JACdf_Sparse(X, FCN, FJAC)

!‘X’輸入向量，計算‘X’處的Jacobian矩陣‘FJAC’作為輸出

!‘FCN’待求根的函數向量，返回‘X’處的函數值‘FF’

FJAC = 0 !Jacobian全部元素先賦0值

h = eps! 自變量增量h; eps為獲取浮點數精度數值的腳本

call FCN (X, FF, N) !調用FCN，計算X點處的函數值FF

Band: do i_Bnd = 1,3 !非零元素帶帶寬方向含有的分塊矩陣的塊數為3

Sub: do i_Sub = 1,j !從坐標的個數j

Harmonic: doi_Hrm = 1,N_phi !諧波基函數的總個數N_phi

Col_1: do i_X = i_Bnd, N_a, 3 !節點的總個數N_a；

每相隔3個節點所對應的自變量獲得微小增量

! 需獲得增量的自變量的編號kk和jj

kk=(i_Sub-1)*N_a*N_phi+(i_X-1)*N_phi+i_Hrm

jj = (i_Sub-1)*N_a*N_phi + (i_X-1)*N_phi + i_Hrm+ j*N_a*N_phi

! 將這兩個自變量添加微小的增量，并賦給中間變量：

Temp_X = X(kk) + h

Temp_Y = X(jj) + h

end do Col

!以上的‘Col循環’獲得了一組添加了增量的自變量，自變量之間的節點序號間隔為3

! 調用子函數FCN，輸入以上這組獲得增量的自變量向量，得到新的函數值向量Wa1, Wa2：

call FCN (Temp_X, Wa1, N)

call FCN (Temp_Y, Wa2, N)

Col_2: do i_X = i_Bnd, N_a，3

Row:do i_Row = 1 to N_a! 沿'行方向'遍歷，查找非零元素

if (i_Row >= i_X - 1 .and. iRow<= i_X + 1) then

! 判斷當前的元素是否在非零元素帶的分塊矩陣內

! 利用Wa1,Wa2,FF,h計算差商，近似當前元素的偏導FJAC(略)

end if

end do Row

end do Col_2

end do Harmonic

end do Sub

end do Band

end subroutine JACdf_Sparse