崔錦濤，汪 諍

（蘭州交通大學(xué) 機電工程學(xué)院，蘭州 730070）

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

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

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

1 物理模型和數(shù)學(xué)模型

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

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

1.1 數(shù)學(xué)模型

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

分段線性彈性力為

取無量綱化參數(shù)為

則無量綱化形式為

其中

1.2 分叉圖

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

圖2 分岔圖Fig.2 Bifurcation diagram

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

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

2 Simulink建模及仿真

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

2.1 積分模塊設(shè)計

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

其中

圖4 積分模塊Fig.4 Integral module

2.2 非線性模塊設(shè)計

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

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

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

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

3 Simulink仿真

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

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

系統(tǒng)的仿真一般經(jīng)歷模型的初始化和仿真執(zhí)行2個階段。

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

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

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

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

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

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

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

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

圖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 結(jié)語

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