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

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

非線性振動(dòng)分析的切比雪夫譜元法

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

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

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

振動(dòng)與波；非線性振動(dòng)；切比雪夫正交多項(xiàng)式；譜元法；牛頓—拉夫遜方法

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

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

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

1 切比雪夫譜元法

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

1.1 振動(dòng)問(wèn)題及其積分形式

考慮振動(dòng)問(wèn)題的一般形式

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

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

其中T表示時(shí)間域。

1.2 時(shí)間單元?jiǎng)澐?/p>

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

譜元法通過(guò)在每一個(gè)單元Ωe中進(jìn)行譜擴(kuò)展來(lái)近似一個(gè)函數(shù)。

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

1.3 振動(dòng)微分方程離散

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

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

圖1 6階切比雪夫Lagrange插值多項(xiàng)式

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

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

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

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

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

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

1.4 邊界條件施加并求解

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

2 非線性振動(dòng)問(wèn)題的牛頓—拉夫遜方法

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

對(duì)于非線性方程組

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

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

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

Newton-Raphson迭代公式表示為

其中ΔX=Xi+1-Xi。

3 Duffing型非線性振動(dòng)方程

Duffing型非線性振動(dòng)方程可以寫(xiě)為

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

近似解析解為

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

4 單擺的非線性振動(dòng)

單擺的非線性振動(dòng)方程可以表達(dá)為

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

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

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

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

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

求出位移響應(yīng)θ(t)，可以獲得兩個(gè)時(shí)間點(diǎn)ti，tj，滿足θ(ti)>0，θ(tj)<0且ti

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

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

5 結(jié)語(yǔ)

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

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

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

[1]聞邦椿，李以農(nóng)，徐培民，等.工程非線性振動(dòng)[M].北京：科學(xué)出版社，2007.

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

[3]唐駕時(shí)，彭海.阻尼對(duì)葉片非線性振動(dòng)的影響[J].噪聲與振動(dòng)控制，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]毛虎平，吳義忠，李建軍，等.時(shí)間譜元法在動(dòng)態(tài)響應(yīng)優(yōu)化中的應(yīng)用[J].振動(dòng)工程學(xué)報(bào)，2013，26(3)：395-403.

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

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

[15]張艷崗，蘇鐵熊，毛虎平，等.動(dòng)態(tài)應(yīng)力解空間譜元離散的關(guān)鍵時(shí)間點(diǎn)識(shí)別方法[J].機(jī)械工程學(xué)報(bào)，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].聲學(xué)學(xué)報(bào)，2007，32(6)：525-533.

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

[19]耿艷輝，秦國(guó)良，王陽(yáng)，等.Galerkin時(shí)空耦合譜元法求解聲波動(dòng)方程[J].聲學(xué)學(xué)報(bào)，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]谷口修著.尹傳家等譯.振動(dòng)工程大全(下)[M].北京：機(jī)械工業(yè)出版社，1986.

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

[25]周凱紅王元?jiǎng)桌畲褐?微分求積法在單擺非線性振動(dòng)分析中的應(yīng)用[J].力學(xué)與實(shí)踐,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

國(guó)家自然科學(xué)基金項(xiàng)目資助(51275489)

毛虎平（1974－），男，副教授，碩士生導(dǎo)師，主要從事振動(dòng)理論與工程數(shù)值分析方法，結(jié)構(gòu)動(dòng)態(tài)響應(yīng)優(yōu)化方法。E-mail:maohp@nuc.edu.cn

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