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

毛虎平 王偉能 續(xù)彥芳 張艷崗 董小瑞

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

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

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

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

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

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

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

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

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

1 切比雪夫譜元法

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

1.1 振動問題及其積分形式

考慮振動問題的一般形式

其中Ar為關(guān)聯(lián)矩陣，關(guān)聯(lián)著質(zhì)量、阻尼和剛度，并假設(shè)與時間t無關(guān)，x，f是時間t的函數(shù)。

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

其中T表示時間域。

1.2 時間單元劃分

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

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

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

1.3 振動微分方程離散

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

應(yīng)用重心插值公式，節(jié)點基函數(shù)可以表示為

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

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

為了獲得一般單元Ωe的節(jié)點基函數(shù)，需要進行節(jié)點坐標轉(zhuǎn)化。那么節(jié)點基函數(shù)在標準單元Ωst和一般單元Ωe中的關(guān)系可以表示為

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

在本研究中，一維坐標轉(zhuǎn)化可以表示為

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

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

1.4 邊界條件施加并求解

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

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

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

對于非線性方程組

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

考慮函數(shù)F:?n→?n，其中

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

Newton-Raphson迭代公式表示為

其中ΔX=Xi+1-Xi。

3 Duffing型非線性振動方程

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

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

近似解析解為

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

4 單擺的非線性振動

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

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

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

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

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

圖4 非線性振動單擺的響應(yīng)

求出位移響應(yīng)θ(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 結(jié)語

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

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

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

[1]聞邦椿，李以農(nóng)，徐培民，等.工程非線性振動[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]毛虎平，吳義忠，李建軍，等.時間譜元法在動態(tài)響應(yīng)優(yōu)化中的應(yīng)用[J].振動工程學報，2013，26(3)：395-403.

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

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

[15]張艷崗，蘇鐵熊，毛虎平，等.動態(tài)應(yīng)力解空間譜元離散的關(guān)鍵時間點識別方法[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]彭海闊，孟光.基于譜元法的梁結(jié)構(gòu)中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].北京：機械工業(yè)出版社，1986.

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

[25]周凱紅王元勛李春植.微分求積法在單擺非線性振動分析中的應(yīng)用[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－），男，副教授，碩士生導師，主要從事振動理論與工程數(shù)值分析方法，結(jié)構(gòu)動態(tài)響應(yīng)優(yōu)化方法。E-mail:maohp@nuc.edu.cn

續(xù)彥芳（1968－），女，副教授，碩士生導師，主要從事武器系統(tǒng)設(shè)計與應(yīng)用研究，結(jié)構(gòu)動力學仿真分析。E-mail:xuyanfang1968@163.com