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

陳趙江　陳 敏

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

?

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

陳趙江陳 敏

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

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

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

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

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

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

(1)

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

Ω=1+εσ

(2)

則式(1)可改寫(xiě)為

(3)

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

(4)

(5)

(6)

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

(7)

(8)

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

(9)

(10)

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

(11)

則式(10)改寫(xiě)成

(12)

或

(13)

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

(14)

(15)

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

(16)

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

φ′=σ-β′

(17)

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

(18)

(19)

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

(20)

(21)

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

(22)

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

(23)

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

(24)

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

(25)

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

(26)

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

(27)

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

(28)

展開(kāi)此行列式，可得

(29)

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

(30)

(31)

其中

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

2.1MATLAB解析求解程序

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

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

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

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

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

F1=0.5; %激勵(lì)幅值

zeta1=0.25;%阻尼項(xiàng)

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

forii=1:1:length(a)

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

%Omega1為較小的激勵(lì)頻率值

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

zeta1)^2));

%Omega1對(duì)應(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為較大的激勵(lì)頻率值

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

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

(epsilon*zeta1)^2));

%Omega2對(duì)應(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)

%如果兩個(gè)本征值相等退出計(jì)算(即此時(shí)計(jì)算到頻響曲線最高點(diǎn))

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

ifFG<0

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

else

plot(Omega1(ii),a(ii),'k.','MarkerSize',15)%穩(wěn)定點(diǎ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;% 對(duì)應(yīng)上文中的(31)式

ifGF<0

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

else

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

holdon

end

xlabel('Omega') %x軸標(biāo)題

ylabel('a')%y軸標(biāo)題

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

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

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

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

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

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

clc

%closeall

clearall

globalepsilongamma1F1zeta1Omega樣% 定義全局變量

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

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

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

F1=0.5;%激勵(lì)幅值

zeta1=0.25; %阻尼項(xiàng)

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

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

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

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];% 下一次計(jì)算的初始值

yy=[yy;ymax];% 保存計(jì)算結(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軸標(biāo)題

ylabel('a')%y軸標(biāo)題

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

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

3　求解結(jié)果分析

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

(b)隨激勵(lì)幅值的變化

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

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

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

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

4　總結(jié)

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

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

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

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

2016-03-08)

③浙江師范大學(xué)“十二五”省級(jí)實(shí)驗(yàn)教學(xué)示范中心重點(diǎn)建設(shè)項(xiàng)目和2015年度浙江師范大學(xué)教改項(xiàng)目資助.