基于MATLAB的強迫振動達芬方程的非線性幅頻響應分析*③

陳趙江　陳 敏

(浙江師范大學數理與信息工程學院　浙江 金華　321004)

?

基于MATLAB的強迫振動達芬方程的非線性幅頻響應分析*③

陳趙江陳 敏

(浙江師范大學數理與信息工程學院浙江 金華321004)

利用MATLAB研究了強迫振動達芬(Duffing)方程的非線性幅頻響應特性，分析了達芬方程非線性幅頻響應近似解析求解和數值求解的方法和步驟，給出了相應的MATLAB求解程序，并將解析解與數值解結果進行了比較.仿真程序和結果能夠加深學生對非線性振動相關知識的理解，提高大學物理及相關力學課程的課堂教學效果.

達芬方程非線性振動非線性幅頻響應跳躍現象MATLAB

強迫振動達芬方程(Duffingequation)在非線性振動領域是一個非常經典的實例，被廣泛用來說明跳躍現象(Jumpingphenomena)、頻響曲線彎曲和其他非線性行為.理解這個簡單的低階非線性方程的振動特性是學習更為復雜的非線性系統的基礎.達芬方程解的跳躍現象與其非線性頻響特性有直接的關系，因此對頻響特性的研究非常重要.在目前大學物理以及相關力學課程的教學中，學生對非線性振動知識了解甚少，一知半解.如果能利用軟件程序幫助學生了解掌握晦澀難懂的非線性振動知識，這將大大有益于相關教學工作的開展.目前已往文獻對達芬方程非線性幅頻響應的分析大多利用MAPLE和MATHEMATICA數學軟件且缺少數值解和解析解的比較分析.本文給出了達芬方程非線性幅頻響應近似解析求解和數值求解的方法和步驟，利用MATLAB的強大計算功能編寫了相應的求解程序，并將解析解與數值解結果進行了比較分析，以滿足學生對非線性振動知識的進一步理解.

1　達芬方程非線性幅頻響應的理論分析

含阻尼和外部驅動力的無量綱達芬方程一般具有如下形式[1]

(1)

其中u是位移，t是時間，ζ是阻尼比，γ是立方剛度系數，F和Ω分別為激勵幅值和激勵頻率.需注意的是，式(1)中所有物理量均為無量綱量，因此Ω其實是激勵頻率和系統固有頻率之比,即

Ω=1+εσ

(2)

則式(1)可改寫為

(3)

我們采用多尺度法[1]對上式進行近似解析求解，引入尺度T0=t和T1=εt

(4)

(5)

(6)

的式(3)的解.把式(4)～(6)代入式(3)，并令ε的同次冪的系數相等，我們得到

(7)

(8)

式(7)的通解可以表示為

(9)

(10)

把A表示成如下的級數形式

(11)

則式(10)改寫成

(12)

或

(13)

將式(13)分成實部和虛部，得出

(14)

(15)

式中a和β都是實數，引進變換

(16)

把式(14)和(15)化為自治系統，由式(16)

φ′=σ-β′

(17)

把式(16)代入式(14)和(15)，并利用式(17)可得

(18)

(19)

由于我們考慮穩態運動，即a和φ為常數值，設a′=0和φ′=0，可得到

(20)

(21)

對式(20)和(21)取平方，將結果相加，得出

(22)

通常稱之為頻率響應方程.利用式(2)我們將上式改寫為Ω關于a的形式， 得到

(23)

從式(23)可知，對于某一幅值響應a，其對應兩個激勵頻率Ω1,2，但這兩個頻率所對應的響應a可能是不穩定的.為了表征解的穩定性，考慮如下的動力學系統方程

(24)

其中，x1和x2為狀態量.則上述系統的雅可比矩陣為

(25)

根據式(18)和式(19)，在本系統中狀態向量分別為a和φ，f1和f2的表達式如下

(26)

因而可得系統的雅可比矩陣為

(27)

穩態運動的穩定性依賴于雅可比矩陣的特征值，式(27)所對應的特征方程如下

(28)

展開此行列式，可得

(29)

本系統的本征值如下

(30)

(31)

其中

2　達芬振動方程非線性幅頻響應分析的MATLAB實現

2.1MATLAB解析求解程序

從式(22)可以看出為了得到頻率響應曲線，可解出a2關于σ的函數，或者可以解出σ關于a的函數.顯然，后一種方法比較簡單.因此，在編程時根據式(23)進行計算.以下為MATLAB求解程序：

%強迫振動達芬方程幅頻響應曲線的解析求解

%參數設定

epsilon=0.2; %小參數

gamma1=4; %非線性參數

F1=0.5; %激勵幅值

zeta1=0.25;%阻尼項

a=linspace(0.1,2.0,100); % 計算響應幅值范圍

forii=1:1:length(a)

%根據(23)式，對應一個響應幅值，可能存在兩個激勵頻率值；

%Omega1為較小的激勵頻率值

Omega1(ii)=(1+3*epsilon*gamma1/8*a(ii)^2-sqrt((epsilon*F1/(2*a(ii)))^2-(epsilon*

zeta1)^2));

%Omega1對應的本征函數值

lamOmega11=sqrt(-((Omega1(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega1(ii)-1)/epsilon-9*gamma1/8*a(ii)^2))-zeta1;

lamOmega12=-sqrt(-((Omega1(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega1(ii)-1)/epsilon-9*gamma1/8*a(ii)^2))-zeta1;

%Omega2為較大的激勵頻率值

Omega2(ii)=(1+3*epsilon*gamma1/

8*a(ii)^2+sqrt((epsilon*F1/(2*a(ii)))^2-

(epsilon*zeta1)^2));

%Omega2對應的本征函數值

lamOmega21=sqrt(-((Omega2(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega2(ii)-1)/epsilon-9*gamma1/8*a(ii)^2))-zeta1;

lamOmega22=-sqrt(-((Omega2(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega2(ii)-1)/epsilon-9*gamma1/8*a(ii)^2))-zeta1;

iflamOmega11==conj(lamOmega21)

%如果兩個本征值相等退出計算(即此時計算到頻響曲線最高點)

plot(Omega1(ii),a(ii),'k.','MarkerSize',15')

holdon

break

else

FG=((Omega1(ii)-1)/epsilon-3*gamma1/8*a(ii)^2)*((Omega1(ii)-1)/epsilon-9*

gamma1/8*a(ii)^2)+zeta1^2;% 對應上文中的(31)式

ifFG<0

plot(Omega1(ii),a(ii),'color',[.5 .5 .5],'MarkerSize',15) %不穩定點

else

plot(Omega1(ii),a(ii),'k.','MarkerSize',15)%穩定點

holdon

end

GF=((Omega2(ii)-1)/epsilon-3*gamma1/

8*a(ii)^2)*((Omega2(ii)-1)/epsilon-9*

gamma1/8*a(ii)^2)+zeta1^2;% 對應上文中的(31)式

ifGF<0

plot(Omega2(ii),a(ii),'.','color',[.5 .5 .5],'MarkerSize',15) %不穩定點

else

plot(Omega2(ii),a(ii),'k.','MarkerSize',15)%穩定點

holdon

end

xlabel('Omega') %x軸標題

ylabel('a')%y軸標題

xlim([0.5,1.5]); %x軸范圍

ylim([0,2]);%y軸范圍

end%滿足條件退出循環(計算到幅頻響應曲線頂點)

end%結束對不同a值的計算

2.2MATLAB數值求解程序

(32)

在MATLAB中定義如下函數并保存為duffing.m

functiondy=duffing(t,y)

globalepsilongamma1F1zeta1Omega

dy=zeros(2,1);

dy(1)=-zeta1*y(1)+F1/2*sin(y(2));

dy(2)=(Omega-1)/epsilon-3*gamma1/

8*y(1)^2+F1/(2*y(1))*cos(y(2));

end

(33)

在MATLAB中定義如下函數并保存為duffing1.m

functiondydt=duffing1(t,y)

globalepsilongamma1F1zeta1Omega

dydt=zeros(2,1);

dydt(1)=y(2);

dydt(2)=-2*epsilon*zeta1*y(2)-y(1)-epsilon*gamma1*y(1)^3+epsilon*F1*cos

(Omega*t);

end

達芬方程幅頻響應數值求解主程序main.m如下：

clc

%closeall

clearall

globalepsilongamma1F1zeta1Omega樣% 定義全局變量

%參數設定

epsilon=0.2; %小參數

gamma1 = 4; % 非線性參數

F1=0.5;%激勵幅值

zeta1=0.25; %阻尼項

np=400;% 計算的頻率數目

Omega1=linspace(.5,1.5,np); % 正向掃頻時的頻率

Omega2=linspace(1.5,.5,np); % 反向掃頻時的頻率

yy=[];

Y0=[0.1 0.1];% 初始值

fori=1:1:length(Omega1)

Omega=Omega1(i);

[T,Y] =ode45(@duffing,[0 400],Y0);

%[T,Y] =ode45(@duffing1,[0 400],Y0);

nn=length(Y(:,1));

ymax=max(Y(nn-round(nn/2):nn,1));% 取穩態幅值

ymax1=max(Y(nn-round(nn/2):nn,2));

Y0=[ymaxymax1];% 下一次計算的初始值

yy=[yy;ymax];% 保存計算結果到yy向量

end

plot(Omega1,yy,'k','LineWidth',1.5)% 正向掃頻曲線用黑色表示

holdon

% 反向掃頻

yy1=[];

ymax1=max(Y(nn-round(nn/2):nn,1));% 反向掃頻的初始值

ymax2=max(Y(nn-round(nn/2):nn,2));% 反向掃頻的初始值

forj=1:1:length(Omega2)

Omega=Omega2(j)

[T,Y1] =ode45(@duffing,[0 400],[ymax1ymax2]);

%[T,Y] =ode45(@duffing1,[0 400],Y0);

nn1=length(Y1(:,1));

ymax1=max(Y1(nn1-round(nn1/2):nn1,1));

ymax2=max(Y1(nn1-round(nn1/2):nn1,2));

yy1=[yy1;ymax1];

end

plot(Omega2,yy1,'color',[.5 .5 .5],

'LineWidth',1.5)% 反向掃頻曲線用灰色表示

xlabel('Omega') %x軸標題

ylabel('a')%y軸標題

xlim([0.5,1.5]); %x軸范圍

ylim([0,2]); %y軸范圍

3　求解結果分析

(a)隨非線性參數的變化

(b)隨激勵幅值的變化

其次，我們利用2.2節的幅頻響應數值求解程序對強迫振動達芬方程的一階近似幅頻響應進行求解分析，其中使用的MATLAB子函數為duffing.m.圖2中黑色和灰色曲線分別為正向和反向掃頻計算結果.從圖2可以發現當正向掃頻(即頻率比Ω增大)時，響應幅值一直增加.當進一步增加時，就發生從A點到B點的跳躍并伴隨著響應幅值的突然減小.隨后，響應幅值隨Ω的增大而減小.對應著A點的最大振幅只有從較低的頻率開始增加才可能達到.在另一方面，反向掃頻(即頻率比Ω減小)時，響應幅值也增加.當Ω進一步減小時就發生從C點到D點的跳躍并伴隨著響應幅值的突然增加.隨后，響應幅值隨Ω的減小而減小.從圖中也可發現解析解中的不穩定區域(從A點到C點)在數值求解時是無法得到的，幅頻響應一階近似解析解和一階近似數值解結果符合很好.

圖2　強迫振動達芬方程的幅頻響應數值與解析解的比較

最后，我們對強迫振動達芬方程的幅頻響應做進一步的數值分析，其中使用的MATLAB子函數為duffing1.m.計算結果如圖3所示，從圖中可以發現數值解與近似解析解存在一定差異，且偏離Ω=1 越大兩者的差別越大.這是因為這里數值求解的是最初的不包含近似的微分方程式(3)，而解析解是滿足Ω≈1的一階近似解.因此，當Ω越接近1時，兩者解的結果越接近.從上述分析也可知近似解析解式(22)是在弱阻尼、弱非線性系數并接近共振時得出的結果，并不適用于所有的情況，這一點在學習的時候要特別注意.

圖3　強迫振動達芬方程的幅頻響應數值與解析解的比較

4　總結

本文對強迫振動達芬方程的非線性幅頻響應特性進行了研究，分析了達芬方程非線性幅頻響應解析求解和數值求解的方法和步驟，給出了相應的MATLAB求解程序，并將解析解與數值解結果進行了比較分析.仿真程序和結果能夠加深學生對非線性振動的理解，提高大學物理以及相關力學課程的教學效果.此外，雖然本文只討論了達芬方程的非線性幅頻響應特性，但仿照本文的求解步驟和程序也容易對達芬方程的非線性相頻特性進行分析.

1A·H·奈弗,D·T·穆克著. 非線性振動. 宋家骕, 羅惟德, 陳守吉譯. 北京：高等教育出版社, 1980

2張志涌. 精通MATLABR2011a. 北京：北京航空航天大學出版社, 2011

陳趙江(1980-)，男，博士，主要從事大學物理教學和研究.

2016-03-08)

③浙江師范大學“十二五”省級實驗教學示范中心重點建設項目和2015年度浙江師范大學教改項目資助.