非線性振動分析的切比雪夫譜元法

毛虎平 王偉能 續彥芳 張艷崗 董小瑞

(1.中北大學 機械與動力工程學院，太原 030051；2.煤科集團杭州環保研究院有限公司，杭州 311201）

非線性振動分析的切比雪夫譜元法

毛虎平1王偉能2續彥芳1張艷崗1董小瑞1

(1.中北大學 機械與動力工程學院，太原 030051；2.煤科集團杭州環保研究院有限公司，杭州 311201）

為了進一步探索Chebyshev時間譜元法求解非線性的振動問題，從Bubnov-Galerkin方法出發，在第二類Chebyshev正交多項式極點處；用重心Lagrange插值來構造節點基函數及其特性，推導了非線性振動問題的伽遼金譜元離散方案，借助Newton-Raphson法求解非線性方程組。對于非線性單擺，還需要將二分法和重心Lagrange插值結合求解角頻率。以Duffing型非線性振動和非線性單擺振動問題為例，驗證了此方法具有現實可行和高精度的優點。

振動與波；非線性振動；切比雪夫正交多項式；譜元法；牛頓—拉夫遜方法

盡管許多工程問題可以用線性振動近似，但還是有很多工程振動需要考慮非線性。例如，大角度單擺、振動輸送機、換熱器直管、葉輪機葉片、高彈聯軸節軸系、高速列車行駛時氣體的阻力及材料產生彈塑性變形構成的振動系統等[1—3]，均需通過非線性微分方程進行分析。非線性振動不符合疊加原理，通常應用數值方法進行分析。

Steven Orszag[4]于1969年提出了譜方法[5,6]之后，給研究者所關注的高精度數值分析帶來希望，然而其不能處理復雜設計域、不能近似非光滑函數等缺點[7]限制了其發展。考慮到譜方法的高精度以及指數收斂和有限元方法處理邊界靈活的特性，學者Patera于1984年提出了譜元法，通過在Gauss-Lobatto-Legendre（GLL）點處Lagrange插值來構造節點基函數，并應用于流體動力學數值分析[8]。30多年來，由于譜元法的高精度和快速收斂的特點得到了極大關注，并被成功應用于科學和工程的很多領域[9-11]。在動態響應優化中，譜元法精確求解動力學控制方程結合高斯—勒讓德—羅巴托（GLL）點以滿足動態約束條件，獲得更好優化的解[12]。在機械故障診斷中，用譜元法模擬帶裂紋的三維板結構的導波激勵與接受以及波的傳播[13]。將仿真時間分為若干步，采用逐步時間譜元法[14]仿真三維懸臂梁，獲得與ANSYS仿真一致的結果，而效率高于ANSYS。文獻[15]將譜元離散方案應用于結構動態應力關鍵時間點識別。Zhao J M[16]采用Chebyshev最小二乘譜元法詳細分析并求解了半透明介質的輻射傳熱。林偉軍[17]應用Modal basis譜元法詳細闡述了彈性波傳播模擬的理論公式，并應用Chebyshev正交多項式展開。彭海闊等[18]通過Legengre譜元法模擬結構彈性波的傳播。秦國良等[19]提出了時空藕合譜元方法，并將其用于帶第一類邊界條件的非齊次一維、二維、三維波動方程的求解。Bar-Yoseph P Z等對非線性一維對流問題、非線性Euler-Bernoulli梁從時—空耦合以及對非線性動力系統應用譜元法進行分析[20—22]。

本文通過在Chebyshev正交多項式極點處重心Lagrange插值構造節點基函數，提出求解非線性振動問題的Chebyshev譜元法。

1 切比雪夫譜元法

譜元近似融合了譜近似和有限元近似的優點，譜元近似可以自由選擇插值次數，獲得p收斂，而有限元近似可以柔性地處理復雜設計域并自由地選擇單元尺寸，獲得h收斂。所謂譜元是正交多項式光滑函數的有限級數。由于數值求解非線性振動問題是以線性振動問題為基礎。因此，首先對線性振動問題進行分析。

1.1 振動問題及其積分形式

考慮振動問題的一般形式

其中Ar為關聯矩陣，關聯著質量、阻尼和剛度，并假設與時間t無關，x，f是時間t的函數。

在切比雪夫譜元法中，為了得到振動問題的數值解，運用Bubnov-Galerkin法，引入一個權函數W，與方程（1）兩邊同時相乘并在時間域上積分，得到了振動問題的積分形式

其中T表示時間域。

1.2 時間單元劃分

作為一種有限元方法，解空間Ω被劃分為Ne個互相不重疊的單元空間，即

譜元法通過在每一個單元Ωe中進行譜擴展來近似一個函數。

將單元節點基函數作為形函數，在單元Ωe上，振動位移可以近似為

1.3 振動微分方程離散

本研究中，采用切比雪夫第二類多項式來構造節點基函數。在標準區間[-1，1]上，N階節點基函數可以表示為拉格朗日插值多項式，其通過N+1個Chebyshev-Gauss-Lobatto點，也就是

應用重心插值公式，節點基函數可以表示為

圖1 6階切比雪夫Lagrange插值多項式

圖1中顯示了節點基函數的科羅尼克δ的特性，這就保證了公式（4）中擴展系數與節點值一致，并且保證施加邊界條件。

為了獲得一般單元Ωe的節點基函數，需要進行節點坐標轉化。那么節點基函數在標準單元Ωst和一般單元Ωe中的關系可以表示為

其中t=t(ξ)定義坐標從一般單元Ωe到標準單元Ωst的轉化，?t是關于t的梯度操作算子，?ξ是關于ξ的梯度操作算子，J是雅可比矩陣是定義在標準單元Ωst上的節點基函數。

在本研究中，一維坐標轉化可以表示為

微分方程左面的表達式，稱為最小二乘譜元法。本研究采用伽遼金譜元法。將φj作為權函數代入式（10），轉化為線性方程組，獲得

矩陣A、B、D通過Gauss-Chebyshev-Lobatto求積公式獲得。

1.4 邊界條件施加并求解

速度初始條件，將K的第一行和第一列除了第一個元素外都強制等于零，對應的F中第一個元素強制等于速度初值；位移初始條件，將K的第(N+1)行和第(N+1)列除了第(N+1,N+1)個元素外都強制等于零，對應F中第(N+1)個元素強制等于位移初值。線性方程組式（11）可以直接求解。

2 非線性振動問題的牛頓—拉夫遜方法

對于非線性振動問題中的非線性項，先直接求微分，再加入到線性振動問題的離散公式中，將其轉化為牛頓—拉夫遜迭代格式進行迭代求解。

對于非線性方程組

其中X(t)—n維解向量，—n維函數向量。

考慮函數F:?n→?n，其中

那么F(x1,x2,…,xn)的雅可比矩陣為

Newton-Raphson迭代公式表示為

其中ΔX=Xi+1-Xi。

3 Duffing型非線性振動方程

Duffing型非線性振動方程可以寫為

其中ε,F是給定的常數，ω是外載荷的頻率，也是常數。

近似解析解為

從圖2、圖3可以看出，本文方法獲得解與近似精確解非常吻合。

4 單擺的非線性振動

單擺的非線性振動方程可以表達為

其中g為重力加速度，l為擺長，θ為擺角。初始條件為

采用單元數10，插值次數6，通過伽遼金離散方案得到非線性方程組，利用Newton-Raphson法求解，當初始擺角時，獲得如圖4所示的擺角、角速度和角加速度，并且與ODE 45求解器計算結果比較，很好的吻合。

圖2 Duffing型振動問題的響應（第一種初始條件）

圖3 Duffing型振動問題的響應（第二種初始條件）

圖4 非線性振動單擺的響應

求出位移響應θ(t)，可以獲得兩個時間點ti，tj，滿足θ(ti)>0，θ(tj)<0且ti

從表1可看出，初始擺角θ0<135°時，本文方法可以獲得最大的絕對誤差0.01%，而2階攝動解最大的絕對誤差為6.1%，DQ法最大的絕對誤差為0.02%。當θ0=150°時，本文方法獲得最大的絕對誤差為1.16%，而2階攝動解最大的絕對誤差為15.85%，DQ法最大的絕對誤差為1.25%。

表1 非線性單擺振動的初始擺角和固有頻率的比值。

5 結語

(1)采用重心Lagrange插值近似單元未知函數，可以獲得精確單元插值微分矩陣，通過有限元節點共享特性可以獲得全局插值微分矩陣，最后獲得非線性代數方程組；

(2)結合Newton-Raphson法，可以同時獲得非線性振動問題的位移和速度，進而通過微分方程中加速度與位移和速度的關系求出加速度；

(3)對于非線性單擺振動，求出角位移后，結合二分法可以精確求出不同初始擺角時的角頻率，并與其他方法比較，說明本文方法精度最高。

[1]聞邦椿，李以農，徐培民，等.工程非線性振動[M].北京：科學出版社，2007.

[2]李安軍，邢桂菊，周麗雯.換熱器直管非線性振動分析與控制[J].噪聲與振動控制，2007，27(5)：50-53.

[3]唐駕時，彭海.阻尼對葉片非線性振動的影響[J].噪聲與振動控制，2013，33(5)：15-18.

[4]Orszag S A.Numerical methods for the simulation of turbulence[J].Physics of Fluids(1958-1988),2004,12(12): II-250-II-257.

[5]Guo B.Spectral methods and their applications[M].World Scientific,1998.

[6]Boyd J P.Chebyshev and Fourier spectral methods[M]. Courier Dover Publications,2013.

[7]Valenciano J,Chaplain M A J.A laguerre-legendre spectral-element method for the solution of partial differential equations on infinite domains:Application to the diffusion of tumour angiogenesis factors[J].Mathematical and Computer Modelling,2005,41(10):1171-1192.

[8]Patera A T.A spectral element method for fluid dynamics: laminar flow in a channel expansion[J].Journal of Computational Physics,1984,54(3):468-488.

[9]High-order methods for incompressible fluid flow[M]. Cambridge University Press,2002.

[10]Zhu W,Kopriva D A.A spectral element approximation to price European options with one asset and stochastic volatility[J].Journal of Scientific Computing,2010,42(3): 426-446.

[11]Zhu W,Kopriva D A.A spectral element approximation to price European options.II.The Black-Scholes model with two underlying assets[J].Journal of Scientific Computing,2009,39(3):323-339.

[12]毛虎平，吳義忠，李建軍，等.時間譜元法在動態響應優化中的應用[J].振動工程學報，2013，26(3)：395-403.

[13]李富才，彭海闊，孫學偉，等.基于譜元法的板結構中導波傳播機理與損傷識別[J].機械工程學報，2013，48 (21)：57-66.

[14]毛虎平，蘇鐵熊，李建軍.基于逐步時間譜元法的結構動態響應仿真[J].中北大學學報:自然科學版，2013，34 (004)：424-430.

[15]張艷崗，蘇鐵熊，毛虎平，等.動態應力解空間譜元離散的關鍵時間點識別方法[J].機械工程學報，2014，50(5)：82-87.

[16]Zhao J M,Liu L H.Least-squares spectral element method for radiative heat transfer in semitransparent media[J].Numerical Heat Transfer,Part B:Fundamentals, 2006,50(5):473-489.

[17]林偉軍.彈性波傳播模擬的Chebyshev譜元法[J].聲學學報，2007，32(6)：525-533.

[18]彭海闊，孟光.基于譜元法的梁結構中Lamb波傳播特性研究[J].噪聲與振動控制，2009，29(6)：62-66.

[19]耿艷輝，秦國良，王陽，等.Galerkin時空耦合譜元法求解聲波動方程[J].聲學學報，2013，38(3)：306-318.

[20]Bar-Yoseph P,Moses E,Zrahia U,et al.Space-time spectral element methods for one-dimensional nonlinear advection-diffusion problems[J].Journal of Computational Physics,1995,119(1):62-74.

[21]Bar-Yoseph P Z,Fisher D,Gottlieb O.Spectral element methods for nonlinear temporal dynamical systems[J].Computational Mechanics,1996,18(4):302-313.

[22]Bar-Yoseph P Z,Fisher D,Gottlieb O.Spectral element methods for nonlinear spatio-temporal dynamics of an Euler-Bernoulli beam[J].Computational Mechanics, 1996,19(1):136-151.

[23]谷口修著.尹傳家等譯.振動工程大全(下)[M].北京：機械工業出版社，1986.

[24]呂中榮，劉濟科.擺的振動分析[J].暨南大學學報:自然科學與醫學版，1999，20(1)：42-45.

[25]周凱紅王元勛李春植.微分求積法在單擺非線性振動分析中的應用[J].力學與實踐,2003,25(3):50-52.

Chebyshev Spectral Element Method forAnalysis of Nonlinear Vibration Problems

MAO Hu-ping1,WANG Wei-neng2,XU Yan-fang1, ZHANG Yan-gang1,DONG Xiao-rui1
（1.College of Mechanical and Power Engineering,North University of China,Taiyuan 030051,China; 2.CCTEG Hangzhou Environmental Research Institute,Hangzhou 311201,China）

The solution of nonlinear vibration problems was studied by using Chebyshev spectral elements method. The node-based functions were constructed by barycentric Lagrange interpolation at the pole points of Chebyshev orthogonal polynomials of the 2nd kind which characteristics were analyzed by using Bubnov-Galerkin method.Galerkin discretization scheme for the nonlinear vibration problems was derived.Finally,the nonlinear equations were solved by Newton-Raphson method.For nonlinear single pendulums,the angular frequencies were solved using the combination of the dichotomy with the barycentric Lagrange interpolation.Two examples of Duffing-type vibration equations and nonlinear vibration of pendulums were employed to illustrate the feasibility and advantages of high-precision of the proposed method.

vibration and wave;nonlinear vibration;Chebyshev orthogonal polynomials;spectral element method; Newton-Raphson method

TB53；TH113.l

：A

：10.3969/j.issn.1006-1335.2015.01.015

1006-1355(2015)01-0073-05

2014－06－03

國家自然科學基金項目資助(51275489)

毛虎平（1974－），男，副教授，碩士生導師，主要從事振動理論與工程數值分析方法，結構動態響應優化方法。E-mail:maohp@nuc.edu.cn

續彥芳（1968－），女，副教授，碩士生導師，主要從事武器系統設計與應用研究，結構動力學仿真分析。E-mail:xuyanfang1968@163.com