嚴(yán)國軍, 楊 彥, 顧建華, 姚月琴，祁 淼,4

(1. 鹽城工業(yè)職業(yè)技術(shù)學(xué)院 機(jī)電工程學(xué)院，江蘇 鹽城 224000; 2. 東南大學(xué) 機(jī)械工程學(xué)院，江蘇 南京 210000； 3. 南京理工大學(xué) 控制工程學(xué)院，江蘇 南京 210000；4. 江蘇大學(xué) 汽車工程學(xué)院，江蘇 鎮(zhèn)江 212000)

0 引 言

在我國，電動(dòng)農(nóng)業(yè)機(jī)械還處于研究階段。無論從市場(chǎng)需求，還是能源發(fā)展方向等方面考慮，電動(dòng)農(nóng)業(yè)機(jī)械都是我國新農(nóng)村建設(shè)和現(xiàn)代農(nóng)業(yè)發(fā)展的需要。我國“十二五”規(guī)劃曾明確提出：“提高設(shè)施裝備水平，提升產(chǎn)業(yè)發(fā)展質(zhì)量，促進(jìn)設(shè)施農(nóng)業(yè)可持續(xù)發(fā)展”。隨著新能源車輛的發(fā)展，作為農(nóng)業(yè)大國，純電動(dòng)拖拉機(jī)正在慢慢替代傳統(tǒng)燃油拖拉機(jī)，并且應(yīng)用在農(nóng)業(yè)的各個(gè)方面。

純電動(dòng)車輛動(dòng)力傳動(dòng)系統(tǒng)中普遍存在機(jī)電耦合現(xiàn)象，研究機(jī)電系統(tǒng)耦合振動(dòng)，可減小系統(tǒng)非線性特性的影響。純電動(dòng)車輛動(dòng)力傳動(dòng)系統(tǒng)的扭轉(zhuǎn)振動(dòng)問題在電磁激勵(lì)和外部激勵(lì)的作用下變得復(fù)雜與突出。在機(jī)電耦合方面，林利紅等人對(duì)純電動(dòng)車輛動(dòng)力傳動(dòng)系統(tǒng)進(jìn)行機(jī)電耦合動(dòng)力學(xué)分析和實(shí)驗(yàn)研究，建立了系統(tǒng)的機(jī)電耦合振動(dòng)仿真模型，仿真分析了由于電流調(diào)節(jié)器參數(shù)、阻尼、諧波擾動(dòng)、間隙以及負(fù)載擾動(dòng)等因素引起的機(jī)電耦合振動(dòng)動(dòng)態(tài)過程，并將物理實(shí)驗(yàn)研究與理論分析結(jié)果進(jìn)行對(duì)比和驗(yàn)證[1]。在車用電機(jī)非線性扭轉(zhuǎn)振動(dòng)方面，陳星等人采用多尺度法研究了純電動(dòng)車輛動(dòng)力傳動(dòng)系統(tǒng)的扭轉(zhuǎn)振動(dòng)系統(tǒng)特性，仿真結(jié)果表明合理的設(shè)計(jì)和控制電機(jī)的電磁參數(shù)和機(jī)械參數(shù)，可以解決電機(jī)扭振系統(tǒng)出現(xiàn)的跳躍分岔等非線性特性現(xiàn)象[2-3]。

針對(duì)低速重載起步工況下純電動(dòng)拖拉機(jī)機(jī)電耦合傳動(dòng)系統(tǒng)受到非周期沖擊干擾力作用，分析扭轉(zhuǎn)角與電磁轉(zhuǎn)矩的耦合影響。為此筆者通過建立純電動(dòng)拖拉機(jī)機(jī)電耦合傳動(dòng)系統(tǒng)的非線性扭轉(zhuǎn)振動(dòng)模型，求解無擾動(dòng)Hamilton系統(tǒng)的平衡點(diǎn)；利用Melnikov法求解分岔閾值曲線，分析電機(jī)控制參數(shù)對(duì)純電動(dòng)拖拉機(jī)機(jī)電耦合傳動(dòng)系統(tǒng)非線性動(dòng)力學(xué)行為的影響；基于Melnikov理論分析研究傳動(dòng)系統(tǒng)的分岔和混沌特性，從而求解系統(tǒng)穩(wěn)定性運(yùn)行邊界條件。

1 機(jī)電耦合扭轉(zhuǎn)振動(dòng)模型

1.1 電機(jī)模型的建立

不同情況下，永磁同步電機(jī)的穩(wěn)定運(yùn)行時(shí)的相量如圖1。

圖1 永磁同步電機(jī)相量Fig. 1 Phasor diagram of permanent magnet synchronous motor

由圖1得到如下關(guān)系式：

(1)

通過求解式(1)可得到電機(jī)定子電流的直、交軸分量：

(2)

定子相電流：

(3)

而電機(jī)的輸入功率：

P1=mUI1cosφ

(4)

當(dāng)忽略定子電阻，由式(4)可得到電機(jī)的電磁功率：

(5)

則電機(jī)的電磁轉(zhuǎn)矩：

(6)

式中：m為電機(jī)相數(shù)；Ω為機(jī)械角速度，且ω=pΩ，其中，ω為電機(jī)的電角速度，p為電機(jī)的極對(duì)數(shù)；E0為空載反電動(dòng)勢(shì)，可由式(7)計(jì)算得到。

(7)

式中：f為電源頻率；Kdp為基波繞組系數(shù)；N為每相串聯(lián)匝數(shù)；Φ10為永磁體基波磁通；Am為永磁體面積；Br為永磁體剩磁；σ0為空載漏磁系數(shù)；bm0為永磁體空載工作點(diǎn)假設(shè)值；KΦ為氣隙磁通波形系數(shù)；αi為計(jì)極弧系數(shù)；αp為極弧系數(shù)，δ為氣隙長度，D為定子內(nèi)徑。

Te=b1sinθ+b2sin 2θ

(8)

由圖1可知，θ=φ+ψ，因此可得到sinθ和sin 2θ三角函數(shù)展開式。

利用sinx≈x-1/6x3，cosx≈1-1/2x2，泰勒展開的近似式將含φ的三角函數(shù)展開并化簡，即可得到作用于轉(zhuǎn)子的電磁轉(zhuǎn)矩：

Te=k0+k1φ-k2φ2-k3φ3

(9)

式中：k0=b1sinψ+b2cos 2ψ；k1=b1cosψ+2b2sinψ；k2=1/2b1sinψ+2b2sin 2ψ；k3=1/6b1cosψ+4/3b2cos 2ψ。

1.2 機(jī)電耦合扭轉(zhuǎn)振動(dòng)模型

通過分析車載永磁同步電機(jī)機(jī)電耦合扭轉(zhuǎn)振動(dòng)特性，將轉(zhuǎn)子系統(tǒng)和傳動(dòng)系統(tǒng)進(jìn)行簡化，可以等效為二質(zhì)體機(jī)電耦合轉(zhuǎn)子模型[4-5]，如圖2。圖中永磁同步電機(jī)轉(zhuǎn)子與機(jī)械轉(zhuǎn)子之間彈性聯(lián)軸器聯(lián)接，電機(jī)機(jī)座可視為剛性。圖中Te是電磁轉(zhuǎn)矩，TL是負(fù)載轉(zhuǎn)矩。

圖2 二質(zhì)體機(jī)電耦合轉(zhuǎn)子系統(tǒng)Fig. 2 Two plastid electromechanical coupled rotor system

根據(jù)牛頓定律，二質(zhì)體轉(zhuǎn)子系統(tǒng)扭轉(zhuǎn)振動(dòng)動(dòng)力學(xué)方程為：

(10)

(11)

當(dāng)轉(zhuǎn)子系統(tǒng)處于穩(wěn)態(tài)運(yùn)行時(shí)，式(10)通過方程變換并消除得到公式中的靜態(tài)項(xiàng)：

(12)

則機(jī)電耦合扭轉(zhuǎn)振動(dòng)動(dòng)力學(xué)方程可轉(zhuǎn)化為無量綱形式：

(13)

式中：v=J1/J2；r=1/(1+v)；μ=(1/J1+1/J2)C；β=k2r2/J1；γ=k3r3/J1；η2=ω0-k1r/J1，ω0為不考慮機(jī)電耦合時(shí)轉(zhuǎn)子系統(tǒng)的機(jī)械固有頻率，可由式(14)進(jìn)行計(jì)算：

(14)

從式(14)可看出，動(dòng)力學(xué)方程具有平方項(xiàng)和立方項(xiàng)，方程為典型的非線性動(dòng)力學(xué)方程。因此，考慮到機(jī)電耦合效應(yīng)，傳動(dòng)系統(tǒng)轉(zhuǎn)子扭轉(zhuǎn)振動(dòng)方程將有明顯的非線性特征。

2 系統(tǒng)平衡點(diǎn)分析

選取的拖拉機(jī)電機(jī)額定功率為7.5 kW，額定電壓為96 V，其他相關(guān)電機(jī)參數(shù)如表1。

表1 電機(jī)參數(shù)表Table 1 Motor parameter list

式中：ε為小參數(shù)，且0<ε?1。

無擾動(dòng)系統(tǒng)的勢(shì)函數(shù)在研究非線性動(dòng)力學(xué)行為中具有十分重要的作用，因此需要討論無擾動(dòng)系統(tǒng)的函數(shù)。當(dāng)ε=0時(shí)，排除阻尼項(xiàng)和擾動(dòng)項(xiàng)，函數(shù)方程變?yōu)闊o擾動(dòng)系統(tǒng)：

(16)

在排除阻尼項(xiàng)和擾動(dòng)項(xiàng)后，轉(zhuǎn)子動(dòng)力學(xué)系統(tǒng)為Hamilton系統(tǒng)，其Hamilton方程為：

(17)

式(17)對(duì)應(yīng)的勢(shì)函數(shù)為：

(18)

Hamilton系統(tǒng)平衡點(diǎn)穩(wěn)定性分為如下幾種情況：

1)情況A

當(dāng)β2≤4γ(1-κ)，Hamilton系統(tǒng)只有一個(gè)平衡點(diǎn)P0(0，0)：

①情況A1：當(dāng)κ<1時(shí)，P0是中心點(diǎn)。

②情況A2：當(dāng)κ>1時(shí)，P0是鞍點(diǎn)。

2)情況B

當(dāng)β2>4γ(1-κ)和4γ(1-κ)>0時(shí)，有3個(gè)平衡點(diǎn)，分別為P0(0，0)，P1(x1，0)，P2(x2，0)：

①情況B1：當(dāng)γ>0時(shí)，P1(x1，0)是鞍點(diǎn)，P0(0，0)和P2(x2，0)是中心點(diǎn)。當(dāng)β=0.8，κ=0.35，γ=0.2時(shí)，可以得到勢(shì)函數(shù)和相軌跡如圖3(a)和圖3(b)，從圖3(a)和圖3(b)中可看出存在橫截同宿軌道。

②情況B2：當(dāng)γ<0時(shí)，P1(x1，0)是中心點(diǎn)，P0(0，0)和P2(x2，0)是鞍點(diǎn)。當(dāng)β=0.8，κ=1.65，γ=-0.2時(shí)，可得到勢(shì)函數(shù)和相軌跡如圖3(c)和3圖(d)，從圖3(c)和3圖(d)中可以看出存在橫截異宿軌道。

圖3 B1和B2情況下Hamilton系統(tǒng)勢(shì)函數(shù)和相軌跡Fig. 3 Potential functions and phase trajectories of Hamilton systems in the case of B1 and B2

3)情況C

當(dāng)β2>4γ(1-κ)和4γ(1-κ)<0時(shí)，有3個(gè)平衡0點(diǎn)，分別為P0(0，0)，P1(x1，0)，P2(x2，0)：

①情況C1：當(dāng)γ>0時(shí)，P0(0，0)是鞍點(diǎn)，P1(x1，0)和P2(x2，0)是中心點(diǎn)。當(dāng)β=0，κ=1.5，γ=0.2時(shí)，可得到勢(shì)函數(shù)和相軌跡如圖4(a)和圖4(b)，從圖4(a)和圖4(b)中可以看出存在橫截同宿軌道。

②情況C2：當(dāng)γ<0時(shí)，P0(0，0)是中心點(diǎn)，P1(x1，0)和P2(x2，0)是鞍點(diǎn)。當(dāng)β=0.8，κ=1.65，γ=-0.2時(shí)，可得到勢(shì)函數(shù)和相軌跡如圖4(c)和圖4(d)，從圖中可以看出存在橫截異宿軌道。

圖4 C1和C2情況下Hamilton系統(tǒng)勢(shì)函數(shù)和相軌跡Fig. 4 Potential functions and phase trajectories of Hamilton systems in the case of C1 and C2

綜合以上分析可知：情況A中只存在一個(gè)平衡點(diǎn)，將不做討論；情況B和C均具有3個(gè)平衡點(diǎn)，而在情況B1和C1中，對(duì)較小的擾動(dòng)幅值，Hamilton系統(tǒng)具有封閉軌道，此時(shí)將產(chǎn)生橫截同宿軌道，可能導(dǎo)致同宿混沌；在情況B2和C2中，增加激勵(lì)幅值，將會(huì)產(chǎn)生橫截異宿軌道，可能導(dǎo)致異宿混沌。

從圖3和圖4中可以看出，考慮機(jī)電耦合作用的Hamilton系統(tǒng)是雙勢(shì)對(duì)稱或不對(duì)稱波動(dòng)系統(tǒng)原型，同時(shí)由于機(jī)電耦合動(dòng)力學(xué)方程系數(shù)的復(fù)雜多變，所以機(jī)電系統(tǒng)扭轉(zhuǎn)振動(dòng)將會(huì)表現(xiàn)出復(fù)雜的非線性動(dòng)力學(xué)行為。車用永磁同步電機(jī)傳動(dòng)系統(tǒng)為獲得更高的傳動(dòng)范圍和高功率密度，通常設(shè)計(jì)和控制于弱磁效應(yīng)，因此，筆者只討論ψ>0的情況。當(dāng)內(nèi)功率因數(shù)角ψ在[0，π]內(nèi)變化時(shí)，Hamilton系統(tǒng)的情況如圖5。分析圖可知，系統(tǒng)只存在B1和C1兩種情況，筆者將以C1為例進(jìn)行進(jìn)一步討論(B1可采用相似的過程進(jìn)行計(jì)算討論)。

圖5 機(jī)電傳動(dòng)系統(tǒng)的平衡點(diǎn)的穩(wěn)定性分類Fig. 5 Stability classification of the equilibrium points ofelectromechanical transmission system

已知例C1中，P0(0，0)是鞍點(diǎn)，P1(x1，0)和P2(x2，0)是中心點(diǎn)。系統(tǒng)Hamilton量為：

(19)

當(dāng)H(0，0)=0時(shí),存在兩條連接鞍點(diǎn)(0，0)的同宿軌道{q0(t)}。同宿軌道求解過程如式(20)：

(20)

由式(20)可得：

(21)

積分得：

(22)

當(dāng)t=0，x0=±1時(shí)，c=0。整理后得兩條同宿軌的參數(shù)方程為：

(23)

3 混沌閾值

當(dāng)擾動(dòng)項(xiàng)增加到Hamilton系統(tǒng)時(shí)，封閉的同宿軌將被破壞，此時(shí)可能會(huì)導(dǎo)致分岔，全局同宿分岔預(yù)示著混沌運(yùn)動(dòng)。Melnikov方法可以估算存在同宿分岔的參數(shù)空間，并可以轉(zhuǎn)化為混沌，從而確定機(jī)電耦合扭轉(zhuǎn)振動(dòng)混沌運(yùn)動(dòng)的閾值[6-7]。因此，在分析同宿分岔時(shí)，將系統(tǒng)(15)寫成狀態(tài)方程的形式為：

(24)

Melnikov方法實(shí)際上是通過計(jì)算雙曲不動(dòng)點(diǎn)的穩(wěn)定和不穩(wěn)定流形之間的距離來確定系統(tǒng)是否存在橫截同宿點(diǎn)。根據(jù)Smale-birkhoff同宿理論，當(dāng)μ≠0，f≠0，且足夠小時(shí)，系統(tǒng)(15)可能存在Smale馬蹄存在意義下的混沌。根據(jù)Melinkov理論，Melnikov方程為：

(25)

式中：q0(τ)=(x±，y±)，代表Hamilton系統(tǒng)的同宿軌道；F∧G=f1g2-f2g1。

將式(24)代入式(25)中，初始條件為t=0，x0=±1，y0=0時(shí)，y±(τ)是時(shí)間τ的奇函數(shù)。則可得到：

(26)

式中，I1和I2分別為部分積分：

所以，對(duì)于一定頻率，當(dāng)滿足式(27)時(shí)，系統(tǒng)可能產(chǎn)生Smale馬蹄變換意義下的混沌。

(27)

式(26)第一項(xiàng)積分I1可通過數(shù)學(xué)軟件計(jì)算直接求得I1=0.372 7；第二項(xiàng)積分I2在求解過程相當(dāng)復(fù)雜，所以采用數(shù)值積分法，通過采用數(shù)值積分法畫出混沌閾值曲線，然后根據(jù)閾值曲線同樣能夠?qū)ο到y(tǒng)做定性和定量分析。Melnikov函數(shù)數(shù)值積分法求解過程如下：

(28)

對(duì)式(28)分離變量，并對(duì)兩邊同時(shí)積分，得到：

(29)

將式(28)和式(29)代入I2中得到：

(30)

式(30)可以計(jì)算得到對(duì)應(yīng)ω值的I2值，代入式(31)即可得出隨ω變化的Melnikov閾值曲線。以f作為控制參數(shù)時(shí)，Melnikov閾值曲線如圖6。在一定條件下，在f-ω平面上存在混沌區(qū)域。當(dāng)f在該區(qū)域取值時(shí)，系統(tǒng)存在Smale馬蹄變換意義下的混沌。

圖6 同宿分岔的閾值曲線(μ=2.5)Fig. 6 Threshold curve of homoclinic bifurcation (μ=2.5)

4 數(shù)值計(jì)算分析

利用Runge-Kutta法對(duì)擾動(dòng)系統(tǒng)(15)進(jìn)行數(shù)值仿真，來進(jìn)行驗(yàn)證上述用Melnikov方法做出的判斷是否正確，從而對(duì)得到對(duì)應(yīng)的時(shí)域波形、相軌跡、龐加萊截面圖和幅值譜進(jìn)行分析。圖7為μ=2.5，ω=1.0，ε=0.1，時(shí)，F(xiàn)=εf，以f為控制參數(shù)時(shí)的系統(tǒng)分岔圖，初始條件為x=0，y=0。

圖7 f為控制參數(shù)時(shí)的系統(tǒng)分叉Fig. 7 Bifurcation diagram of the system when f is thecontrol parameter

從圖7中可以看出，當(dāng)f為控制參數(shù)時(shí)，機(jī)電耦合非線性扭振系統(tǒng)表現(xiàn)出復(fù)雜的非線性振動(dòng)現(xiàn)象[8]：當(dāng)激勵(lì)值f較小時(shí)(f<3.6)，系統(tǒng)為1周期振蕩；當(dāng)激勵(lì)值f增大到3.6～3.8區(qū)間，2周期振蕩取代1周期振蕩；當(dāng)f>3.8時(shí)，控制參數(shù)f進(jìn)一步增加，系統(tǒng)在倍周期分岔后進(jìn)入混沌運(yùn)動(dòng)；當(dāng)f>3.8時(shí)，機(jī)電耦合系統(tǒng)表現(xiàn)為周期運(yùn)動(dòng)與混沌運(yùn)動(dòng)交替存在。

便于更好的理解，選取3組數(shù)值對(duì)周期運(yùn)動(dòng)變?yōu)榛煦邕\(yùn)動(dòng)的過程進(jìn)行具體分析系統(tǒng)，分別選取控制參數(shù)f=2，f=3.2，f=4，如圖6(b)。當(dāng)ω=1.0時(shí)，控制參數(shù)f由0逐漸增加，橫截同宿點(diǎn)橫截交叉發(fā)生于f≥3.6處。圖8～圖10表示了3種情況下時(shí)域波形、相軌跡、龐加萊截面和幅值譜。f=2對(duì)應(yīng)的點(diǎn)處于閾值曲線之下，沒有達(dá)到混沌運(yùn)動(dòng)的前提條件，如圖6。

圖8 f=2時(shí)1周期運(yùn)動(dòng)Fig. 8 1 period motion graph when f=2

圖9 f=3.2時(shí)2周期運(yùn)動(dòng)Fig. 9 2 period motion graph when f=3.2

圖10 f=4時(shí)混沌態(tài)Fig. 10 Chaotic state graph when f=4

由圖8可以看出,f=2時(shí)，機(jī)電耦合扭轉(zhuǎn)系統(tǒng)為1周期運(yùn)動(dòng)；當(dāng)f開始增大并超過臨界閾值，系統(tǒng)開始失去穩(wěn)定性，并存在發(fā)生混沌運(yùn)動(dòng)可能性。從圖9可以看出，f=3.2時(shí)，此時(shí)運(yùn)行工況點(diǎn)處于閾值曲線之上，系統(tǒng)運(yùn)動(dòng)為2周期運(yùn)動(dòng)，系統(tǒng)發(fā)生倍周期分岔，但并未發(fā)生混沌運(yùn)動(dòng)。從圖10可以看出：f=4時(shí)，運(yùn)行工況點(diǎn)處于閾值曲線之上，進(jìn)入如圖7中的混沌區(qū)域，其計(jì)算結(jié)果如圖10；時(shí)域波形并無周期性，相軌跡為多個(gè)交叉的環(huán)形圈，龐加萊截面圖具有雜亂無章的多個(gè)點(diǎn)，頻譜分析結(jié)果具有能量集中區(qū)域，這種現(xiàn)象表明該情況下系統(tǒng)發(fā)生混沌運(yùn)動(dòng)[9-10]，這與我們之前Melnikov方法解析得到的結(jié)果是一致的。

5 結(jié) 論

以永磁同步電機(jī)的機(jī)電復(fù)合傳動(dòng)系統(tǒng)為研究對(duì)象，考慮機(jī)電耦合關(guān)系，建立永磁同步電機(jī)的電磁轉(zhuǎn)矩模型，并通過非線性動(dòng)力學(xué)理論得到電磁參數(shù)作用下的機(jī)電耦合參數(shù)振動(dòng)界限以及規(guī)律。分析得到如下結(jié)論：

1)產(chǎn)生非線性振動(dòng)的主要原因是由于車用永磁同步電機(jī)扭轉(zhuǎn)角引起的電磁激勵(lì)。

2)考慮機(jī)電耦合作用，扭轉(zhuǎn)振動(dòng)系統(tǒng)的固有特性與電機(jī)的運(yùn)行狀態(tài)以及結(jié)構(gòu)參數(shù)相關(guān)，且在某些情況下，固有頻率有較大的降幅。

3)考慮到車用電機(jī)采用弱磁控制，并且固有頻率受到內(nèi)功率因數(shù)角的影響，因此車用永磁同步電機(jī)需要合理控制內(nèi)功率因數(shù)角的范圍。

4)純電動(dòng)拖拉機(jī)機(jī)電耦合傳動(dòng)系統(tǒng)的控制參數(shù)激勵(lì)值f逐漸增大時(shí)，系統(tǒng)具有周期、擬周期和混沌等復(fù)雜的動(dòng)力學(xué)行為。當(dāng)出現(xiàn)非線性振動(dòng)時(shí)，機(jī)電耦合傳動(dòng)系統(tǒng)受到非周期沖擊干擾力作用，產(chǎn)生較大強(qiáng)度動(dòng)載荷并引起非線性扭振，嚴(yán)重影響純電動(dòng)拖拉機(jī)機(jī)電耦合傳動(dòng)系統(tǒng)可靠性和耐久性。因此，應(yīng)避免純電動(dòng)拖拉機(jī)在低速重載工況下長時(shí)間運(yùn)行。