陳趙江　陳 敏

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

?

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

陳趙江陳 敏

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

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

達芬方程非線性振動非線性幅頻響應(yīng)跳躍現(xiàn)象MATLAB

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

1　達芬方程非線性幅頻響應(yīng)的理論分析

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

(1)

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

Ω=1+εσ

(2)

則式(1)可改寫為

(3)

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

(4)

(5)

(6)

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

(7)

(8)

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

(9)

(10)

把A表示成如下的級數(shù)形式

(11)

則式(10)改寫成

(12)

或

(13)

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

(14)

(15)

式中a和β都是實數(shù)，引進變換

(16)

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

φ′=σ-β′

(17)

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

(18)

(19)

由于我們考慮穩(wěn)態(tài)運動，即a和φ為常數(shù)值，設(shè)a′=0和φ′=0，可得到

(20)

(21)

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

(22)

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

(23)

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

(24)

其中，x1和x2為狀態(tài)量.則上述系統(tǒng)的雅可比矩陣為

(25)

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

(26)

因而可得系統(tǒng)的雅可比矩陣為

(27)

穩(wěn)態(tài)運動的穩(wěn)定性依賴于雅可比矩陣的特征值，式(27)所對應(yīng)的特征方程如下

(28)

展開此行列式，可得

(29)

本系統(tǒng)的本征值如下

(30)

(31)

其中

2　達芬振動方程非線性幅頻響應(yīng)分析的MATLAB實現(xiàn)

2.1MATLAB解析求解程序

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

%強迫振動達芬方程幅頻響應(yīng)曲線的解析求解

%參數(shù)設(shè)定

epsilon=0.2; %小參數(shù)

gamma1=4; %非線性參數(shù)

F1=0.5; %激勵幅值

zeta1=0.25;%阻尼項

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

forii=1:1:length(a)

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

%Omega1為較小的激勵頻率值

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

zeta1)^2));

%Omega1對應(yīng)的本征函數(shù)值

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對應(yīng)的本征函數(shù)值

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;% 對應(yīng)上文中的(31)式

ifFG<0

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

else

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

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;% 對應(yīng)上文中的(31)式

ifGF<0

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

else

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

holdon

end

xlabel('Omega') %x軸標題

ylabel('a')%y軸標題

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

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

end%滿足條件退出循環(huán)(計算到幅頻響應(yīng)曲線頂點)

end%結(jié)束對不同a值的計算

2.2MATLAB數(shù)值求解程序

(32)

在MATLAB中定義如下函數(shù)并保存為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中定義如下函數(shù)并保存為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

達芬方程幅頻響應(yīng)數(shù)值求解主程序main.m如下：

clc

%closeall

clearall

globalepsilongamma1F1zeta1Omega樣% 定義全局變量

%參數(shù)設(shè)定

epsilon=0.2; %小參數(shù)

gamma1 = 4; % 非線性參數(shù)

F1=0.5;%激勵幅值

zeta1=0.25; %阻尼項

np=400;% 計算的頻率數(shù)目

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));% 取穩(wěn)態(tài)幅值

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

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

yy=[yy;ymax];% 保存計算結(jié)果到y(tǒng)y向量

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　求解結(jié)果分析

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

(b)隨激勵幅值的變化

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

圖2　強迫振動達芬方程的幅頻響應(yīng)數(shù)值與解析解的比較

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

圖3　強迫振動達芬方程的幅頻響應(yīng)數(shù)值與解析解的比較

4　總結(jié)

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

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

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

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

2016-03-08)

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