對稱間隙單自由度振動系統的Simulink仿真

崔錦濤，汪 諍

（蘭州交通大學 機電工程學院，蘭州 730070）

實際機械系統中廣泛存在著各種非線性因素，因此工程實際中的振動系統絕大多數屬于非線性系統。近些年來，人們對于實際工程中含間隙振動系統的問題越來越重視[1-4]。

對于非線性振動問題的求解，多借助于數值方法[5-8]，但是這種方法往往需要編制復雜的計算程序，也會借助多提動力學軟件進行仿真，然而軟件仿真往往針對特定的模型進行的（例如齒輪系統、連桿機構、鏈傳動等）[9-11]。對于一些特定的模型，大型通用多體動力學軟件往往沒有相應的模型。因此，要想通過多體動力學軟件進行非線性問題的仿真并得到較為精準的結果，還是需要與其他方法（有限元分析、過程控制）相結合進行聯合仿真。

軟件本文將在對稱間隙碰撞系統模型基礎上，實現了一種基于Simulink的含間隙非線性動力學仿真模型。與傳統的數值仿真算法程序相比，Simulink是一種框圖式的動態系統建模仿真工具，可以實現多個變量同時輸出顯示，并且在模型運行過程中動態顯示變化過程。并且，此方法對于實現Simulink與多體動力學軟件進行含間隙類非線性動力學聯合仿真提供了參考意義。

1 物理模型和數學模型

典型的對稱間隙單自由度碰撞模型如圖1所示。質量塊M通過彈簧（K0）阻尼（C）系統連接，振子M在簡諧波激勵的作用下沿著水平光滑平面運動。設線性彈簧在完全釋放時的位置為零點，以零點建立向右的一維坐標系統。當振子的位移＞B（或＜-B）時，受到剛度為K1的彈性約束，在簡諧力的作用下，如此往復發生碰撞。

圖1 對稱間隙單自由度碰撞模型Fig.1 Single degree of freedom collision model with symmetrical gap

1.1 數學模型

圖1所示對稱間隙單自由度振動系統為典型的非線性振動系統。其主要含有分段線性彈性力F（X），其數學模型為

分段線性彈性力為

取無量綱化參數為

則無量綱化形式為

其中

1.2 分叉圖

取式（2）中的參數 μk=32，ξ=0.23，δ=0.02。 采用龍格庫塔法對系統進行數值積分，得到系統的分岔圖如圖2所示，圖中橫坐標為激振頻率ω，縱坐標為振子運動到坐標零點時的瞬時速度。由圖可見，系統對ω的變化極為敏感，并且隨著激振頻率的改變，系統出現了周期三分岔、混沌、倍化分岔、周期運動等特性。

圖2 分岔圖Fig.2 Bifurcation diagram

為了更好地描述分岔圖中這幾種典型的運動特性，取 ω 分別為 2.1，2.115，2.166，2.188；對系統進行數值仿真，得到了周期三分岔、混沌、倍化分岔、周期運動相對應的相圖，如圖3所示。

圖3 系統在不同激振頻率下的相圖Fig.3 Phase diagram of system on different excitation frequency

2 Simulink建模及仿真

Simulink可以實現建模、仿真、動態分析的可視化仿真工具，以各種類型的子模塊可以搭建出線性系統、非線性系統，以及數字信號處理系統。一般來說，1個動態Simulink系統的構成分為3個部分，即信號的輸入、主系統、系統的響應。文中使用連續系統模塊Continuous，數學運算模塊Math Operations，信號源模塊Source等，搭建出對稱間隙單自由度碰撞系統的模型，并使用信號流路模塊Signal Routing確定系統的穩態輸出；利用雙標量圖形模塊XY Graph繪制系統相圖[12-14]。

2.1 積分模塊設計

在Simulink中，由微分方程建立仿真模型通常有3種方法：①將微分方程轉換成傳遞函數Discrete TransterFcn；②將微分方程轉換成空間狀態方程State-Space；③利用過程信號連續時間積分模塊Integrator直接建模。在此采用方法③。

其中

圖4 積分模塊Fig.4 Integral module

2.2 非線性模塊設計

非線性部分即式（3）中的 f（y1），對其進行轉換，得

式（4）中的數學運算可采用Math Operation模塊庫中的Product數據相乘模塊、Add信號相加模塊實現；分段邏輯關系可采用Commonly Used Blocks模塊庫中的Switch模塊實現。Switch模塊基于第2輸入端（中間）的值，判斷輸出第1輸入端（上部）還是第3輸入端（底部）。當第2輸入端的值大于模塊閾值（P值）時，第1輸入端作為輸出；反之第3輸入端作為輸出。

由于Switch模塊的通過判斷標準Criteria for passing first input，僅包含中間輸入是否大于等于閾值（U2≥Threshold），因此對式（4）中的定義域做如下轉變：

其中，采用Gain信號增益模塊，對輸入信號x取負增益可以得到-x。

3 Simulink仿真

Simulink中模塊的仿真過程可以解釋為模塊不斷地接受信號、更新狀態和計算輸出過程。同樣，模型的仿真是系統在特定的仿真時間段內不斷根據源信號計算狀態和輸出的過程。連續時間系統的仿真過程是將連續的時域離散成相應的采樣時間點，然后通過數值積分的方式求解響應，Simulink將采樣時間點分為主時間步和微時間步2類，后者是對前者的進一步細分。

連續時間系統中的離散狀態更新僅發生在主時間步上，而在2個主時間步之間保持不變。同時，在仿真結果輸出中，Simulink只顯示主時間步的狀態和輸出（實際上主時間步的結果是微時間步上數值積分累積而成的）。微時間步上進行的工作，則是計算連續狀態的微分以更新連續狀態，同時計算每個模塊的輸出。

系統的仿真一般經歷模型的初始化和仿真執行2個階段。

初始化階段在此階段Simulink將進行以下操作：

1）模塊的參數估值。模型中很多模塊都具有參數（可能是變量或表達式）其具體值由MatLab確定。

2）模型的各個子系統被展開，取而代之的是一個大系統，然后模型中的模塊按更新次序排序（即仿真執行時的順序）。模塊的更新順序存在2個原則：①每個模塊必須在它所驅動的任何一個模塊更新之前更新；②非直接饋入模塊可以按任何次序更新，只要它們在其所驅動的直接饋入模塊之前更新。

3）確定每個模塊輸入輸出信號的屬性，諸如數據類型、數值類型、信號寬度（標量、向量或矩陣）。若模塊沒有明確的信號屬性，Simulink采用屬性傳遞的方法確定待定模塊的特性，即未知模塊的屬性由驅動它的模塊屬性唯一決定。

4）確定模型中所有模塊的采樣時間，分配和初始化用于存儲每個模塊狀態和輸出當前值的存儲空間。

模型執行階段在模型的執行階段Simulink通過數值積分實現仿真。所運用的仿真解法器取決于模型提供它的連續狀態的微分能力。計算微分可以分成2步：首先，按照排序所確定的次序計算每個模塊的輸出；然后再根據當前時刻其輸入狀態來確定狀態的微分，得到微分向量后再將其返回解法器，后者用它計算下一采樣時間點的狀態向量。一旦新的狀態向量計算完畢，被采樣的數據源模塊和接收模塊才被更新。

Simulink提供給用戶一系列的數值積分器。對于連續時間系統，根據步長的變化情況分為固定步長和可變步長兩大類：1）固定步長求解器 ODE1（一階 Euler法）、ODE2 （二階改進 Euler法）、ODE4（四階Runge-Kutta法）；2）可變步長求解器 ODE45（四五階 Runge-Kutta法）、ODE23（二三階 Runge-Kutta法）、ODE113（二三階 Adams-Bashforth 法）。

選用ODE45作為積分求解器，系統的信號輸出采用雙標量圖形輸出模塊XY Graph，第1坐標量輸入y1，第2坐標量輸入y2；為了得到系統穩態響應，以系統仿真時間為信號，采用輸出選擇開關（Switch）控制閾值的方法控制信號輸出。根據數學模型，將模塊輸入輸出參數之間正確連接，最終搭建完成對稱間隙單自由度振動系統的Simulink仿真模型，設計如圖5所示。設置參數U=32，S=0.23，P=0.02；仿真參數設置如下：

圖5 Simulink中的仿真模型Fig.5 Simulation model use Simulink

仿真開始時間為0/s；

仿真停止時間為200/s；

取樣開始時間為100/s；

取樣停止時間為200/s；

最大步長為0.02；

最小步長為0.001；

求解器為Ode45（Dormand-Prince）；

相對誤差為10-6；

絕對誤差為10-6。

取正弦信號（Sine Wave1）頻率分別為2.1，2.115，2.166，2.188；運行模型得到不同激振頻率下的相圖，如圖6所示。

圖6 Simulink模型中仿真的相圖Fig.6 Simulation phase diagram use Simulink model

4 結語

采用Simulink可以建立與對稱間隙單自由度振動系統數學模型完全等效的仿真模型，建模過程簡單，不需要復雜繁瑣的編程；以分模塊的方式建立系統（包含積分模塊、非線性模塊、時間控制模塊），使得龐大的動力模型構建過程變得簡單易行。Simulink具備優秀的積分算法和豐富的求解器類型，可以達到數值計算的精度水準，能夠滿足工程實際應用。Simulink的示波器可以實現動態實時顯示圖像的特點，并且在運行過程中可以調整模型參數進行假設仿真分析，監視參數變化對運動特性的影響。這種方法的效率要高于程序計算，且操作簡便。