強非線性振動系統的通用化求解程序及應用

張琪昌，趙薔薇，王 煒

(1.天津大學 機械學院，天津市非線性動力學與混濁控制重點實驗室，天津 300072;2.天津大學 內燃機燃燒學國家重點實驗室，天津 300072)

就強非線性振動系統而言，相關定量理論發展至今已經產生了一系列有價值的研究方法，內容涉及穩態響應計算以及分岔分析兩個主要方面。這其中，張佑啟、陳樹輝等［1-2］引入變換 α =α(ε，ω0，ω1)，提出改進的LP方法(MLP)，利用小參數α將給定的強非線性系統轉化為弱非線性系統;唐駕時［3-4］在MLP方法中引入新的參數變換 α=α(ε，nω0/m，ω1)，解決了強非線性系統次諧共振解問題;張琪昌等基于Nayfeh提出的復規范形理論［5］，提出了待定固有頻率法［6-7］，該方法通過在系統中引入新的未知頻率ω1取代復規范形理論中的原有基頻ω0，將規范形理論的適用范圍由弱非線性振動系統拓展到了強非線性振動系統領域，獲得了系統的高精度穩態漸近解，并在后續工作中開展了有關低維系統的靜態與動態分岔問題分析。在此基礎上，郝淑英、王煒等［8-9］采用該方法計算了兩自由度強非線性振動系統的穩態響應，張琪昌、萬浩川等［10］則將這種方法進一步推廣到三自由度系統。

考慮到高維系統多頻激勵共振關系的確定比較復雜，同時規范形求解過程中涉及近恒同變換系數計算的重復性工作較多，因此時至今日有關待定固有頻率法在強非線性振動系統中的應用仍以低維、經典系統為主，且不同模態振型間僅維持簡單的內共振關系。鑒于此，本文針對多自由度強非線性振動系統，采用待定固有頻率法，運用符號運算語言Mathematica編寫了高階穩態響應計算的通用化求解程序。首先應用模塊化思想對已有理論進行優化，將整體分析流程劃分為相對獨立的計算單元。此后按照邏輯關系編制運算程序，并且設計了相對獨立單元之間的程序接口。利用該程序化方法，研究了一類多自由度碳納米管強非線性振動系統的穩態響應，考察了具體參數取值情況下阻尼系數和激振力對于振動的影幅值響。最后將本文得到的解析解與數值結果和傳統規范形結果進行了對比，不但驗證了漸近解結果的求解精度，也充分體現了通用化程序設計思想在提高待定固有頻率法分析復雜問題計算效率方面的有效性。

1 多自由度強非線性振動系統待定固有頻率法

考慮n自由度強非線性振動系統:

其中:kji為線性項系數，fj為非線性項表達式，Fjcos(Ωjt)為外加激勵，j=1，…，n(n為系統自由度)。首先要對式(1)進行解耦。通過線性變換，方程(1)簡化為如下形式:

考慮到強非線性振動系統振動中非充分小擾動量對于振動基頻的影響，引入新的待定固有頻率ωj0替代傳統復規范形算法中采用的系統原有基頻ωj，即=，…，n(n為系統自由度)。由此原方程轉化為:

引入復變量 ξ1，…，ξn和z1，…，zn，令:

求解方程組(3)的前兩式，得到:

將式(4)中第一式對t求導，并將式(2)、式(3)代入式(4)，得到:

依據規范形理論，引入如下形式的近恒同非線性變換:

其中:

其中，k為非線性項階次，k=1，…，s，Γj，k為引入的非線性變換系數。以僅含有三階非線性項的系統為例，將式(6)代入式(5)，得到:

2 多自由度強非線性振動系統的規范形

為對式(7)進行簡化，需要確定系統的共振項。本文考慮到的共振情形主要包含以下幾類:

(1)非共振情形，即N0Ωj≠N1ω10…≠Njωj0(N0-Nj為沒有公約數的正整數，j=1，…，n);

(3)次諧共振，即 Ωj≈sωj0;

(5)內共振情形，包含:mrωr0≈mtωt0;或msωs0≈mrωr0+mtωt0，其中ms，mr，mt為互質正整數;

(6)組合共振，即上述共振形式的組合情況，可以有幾種共振形式同時存在。

由于多自由度系統往往體現出復雜的共振關系，因此針對具體問題逐一計算規范形的環節，成為制約強非線性振動系統穩態響應計算和影響運算效率的主要因素。為此，本文程序將待定固有頻率法的運算步驟嵌入到已有半單及非半單分岔系統規范形的Mathematica運算程序之中，實現了由計算機程序自動篩選相應的共振項，完成有限階截斷規范形的計算功能。

將其代入式(7)，分離實部與虛部，根據系統具有定常解條件=0，可以得到待定瞬時頻率ωj0、振幅aj，以及漸近解vj的表達式。

3 程序化方法思路

在本文的求解程序中，將整體分析流程劃分為:多自由度系統的簡化解耦、共振項篩選、規范形計算、平均方程獲取、穩態響應計算(待定固有頻率、振幅，漸近解)等相對獨立的計算單元;同時在各單元之間設置了數據流傳輸接口。共振條件、階次u作為初始條件存入文檔CHZ.txt中。程序運行過程中，根據消除vj中非共振項的條件，由hj1至hju依次求出各階非線性變換項系數。該程序適用于具有常規非線性項的多自由度系統，在輸入初始條件(階次與共振條件)之后，就可以直接得到系統的規范形和穩態漸近解，也可以由得到的平均方程開展靜態分岔分析，具有通用性。

程序框圖如圖1。

圖1 程序框圖Fig.1 Procedure structure diagram

4 研究碳納米管的強非線性振動規律

碳納米管復合材料集小尺度、低密度、高強度以及優良的力電性質于一身，其制造和應用作為一個全新的科學領域，是目前研究的熱點問題之一。作為一種性能優異的納米材料，其適用于納米尺度下，制作微機械器件、塔接電路等。另外，由于小直徑、高強度等特點，碳納米管探針在高分辨率原子力顯微鏡、DNA序列測定等領域中發揮了重大作用［11-12］。

碳納米管復合材料所具有的特殊的力電性質和其本身的結構特性與受力變形形態相關，因此，研究在外部激振力作用下碳納米管的力學行為是當前熱門問題，為進一步發展和完善納米材料的應用、優化設計等方面提供非常重要的理論基礎和參考價值。

圖2 雙壁碳納米管受力圖Fig.2 Vibrating load on double-walled carbon nanotube

4.1 建立方程

由于較大的長徑比，碳納米管的力學分析模型通常等效為橫向受力梁。本文以耦合雙壁碳納米管為例［13］，外加橫向簡諧激振力。D0代表外徑，Di代表內徑，如圖2。

其內外管壁耦合振動方程如下:

其中，w1，w2為內外壁彎曲撓度，E1，E2為彈性模量，I1為慣性矩，A1，A2為截面積，ρ1，ρ2為管密度，c1，c2為范德華作用系數，ξ為粘滯阻尼系數。

令wi=Xi(x)ui(t)，Xi(x)=(1-cos(2πx/L))，i=1，2。對上述方程組進行伽遼金離散，引入內管半徑R及納米管橫向振動基頻［14］w0，得到無量綱化系統:

方程組系數見附錄。

4.2 計算系統的穩態響應

運行求解程序，首先對(9)式線性部分進行解偶，得到兩自由度控制方程:

相對線性力而言，方程中非線性力處于更高量級(由(10)式系數計算結果可以看出)，系統為強非線性振動系統。內管壁振幅很小，以下只分析外層管壁的運動行為。

系統穩態漸近解表達式及相圖如下。由圖3可見，程序化方法具有很高的精度。

粘滯阻尼系數是影響納米管振幅的重要因素，在穩態漸近解表達式中，非線性項中包含阻尼系數，初始階段，阻尼系數增長會使納米管的振幅增大，隨著阻尼系數的繼續增加，納米管振幅會有明顯減小，體現了阻尼系數對振動的抑制程度，如圖4所示。實際應用中，復合材料的粘滯阻尼系數與多種因素相關，包含溫度、組成材料含量比、空氣環境等，通常可以調節這些因素以增大或減小振動幅值。

圖3 碳納米管-外壁運動相圖Fig.3 Phase pattern of outer wall

外加激勵幅值對納米管的振動具有很大影響，在周期解存在的情況下，隨著振幅增加，納米管振動幅值增加，如圖5所示。激勵因素可以人為實施控制，以達到抑振目的。

圖4 粘滯阻尼系數對外壁振幅的影響Fig.4 Influence of viscous damping coefficient on amplitude of outer wall

圖5 激振力幅值對外壁振幅的影響Fig.5 Influence of exciting force amplitude on amplitude of outer wall

5 結論

本文進行了如下研究:

(1)采用待定瞬時固有頻率法，運用Mathematica程序語言編寫了計算強非線性振動系統規范形的通用程序，拓展了多自由度系統規范形計算的范圍，使其不再局限于某一類方程形式，并且考慮了復雜的共振情況，由計算機代替人工計算，提高了運算效率。

(2)采用雙壁碳納米管具體實例，應用連續梁理論，將其橫向振動表示為內外管壁耦合非線性微分方程，用程序進行分析，得到了系統穩態解表達式及相圖，與Nayfeh解相比較，具有更高的精度，證明了程序化方法的有效性。

(3)得到了粘滯阻尼系數、激振力幅值與納米管外壁振動振幅的關系曲線，經過分析每個參數對于振動影響，得到結論:激振力幅值增大會使納米管外壁振幅增加，隨著粘滯阻尼系數的增大，振幅先有小幅度的增加后迅速減小。此結論可用于指導具有工程背景的實際問題及微納米構件的優化設計，同時證明了本文程序的通用性與實用性。

［1］陳樹輝.強非線性振動系統的定量分析方法［M］.北京:科學出版社，2007.

［2］ Yang C H，Zhu S M，Chen S H.A modified elliptic lindstedt-poincaré method forcertain strongly non-linear oscillators［J］.Journal of Sound and Vibration，2004，273(4-5):921-932.

［3］唐駕時.求強非線性系統次諧共振解的MLP方法［J］.應用數學和力學，2000，21(10):1039-1045.

［4］Hu H，Tang J H.A convolution integral method for certain strongly nonlinear oscillations［J］.Journal of Sound and Vibration，2005，285(4-5):1235-1241.

［5］Nayfeh A H.Method of normal forms［M］.New York，John Wiley＆ Sons，1993.

［6］張琪昌，郝淑英，陳予恕.用范式理論研究強非線性振動問題［J］.振動工程學報，2000，13(3):481-486.

［7］張琪昌，王 煒，郝淑英.研究強非線性振動問題的最簡規范形方法［J］.應用力學學報，2008，25(4):698-702.

［8］郝淑英，王 煒，張琪昌.研究兩自由度強非線性振動系統的規范形方法［J］.振動工程學報，2007，20(4):422-426.

［9］王 煒，張琪昌，田瑞蘭.兩自由度強非線性振動系統的漸近解及分岔分析［J］.振動與沖擊，2008，27(5):130-133.

［10］萬浩川，張琪昌，王 煒.復規范形理論在研究三自由度強非線性振動問題中的應用［J］.振動與沖擊，2010，29(8):189-194.

［11］ Nessim G D.Properties，synthesis，and growth mechanisms of carbon nanotubes with special focus on thermal chemical vapor deposition［J］.Nanoscale，2010，2(8):1306-1323.

［12］ Hoseinzadeh M S，Khadem S E.Thermoelastic vibration and demping analysis of double-walled carbon nanotubes based on shell theory［J］.Physica E，2011，48(6):1146-1154.

［13］ Hawwa M A，Al-Qahtani H M.Nonlinear oscillations of a double-walled carbon nanotube［J］.Computational Materials Science，2010，48(1):140-143.

［14］Natsuki T，Lei X W，Ni Q Q，et al.Vibrational analysis of double-walled carbon nanotubes with inner and outer nanotubes of different lengths［J］.Physics Letters A，2010，374(46):4684-4689.

附錄

參數取值:ρ1=ρ2=1.13 ×103kg/m3，E1=0.926 TPa，E2=0.923 TPa，Di=5 nm，t=0.34 nm，c1=71.11 GPa，c2=2.57 ×104GPa/nm2，L=60 nm，ξ=0.1，F=4 N，Ω =1.8 rad/s。