單擺運(yùn)動(dòng)過程的MATLAB&Simulink建模與仿真

1、單擺運(yùn)動(dòng)過程的MATLA&Simulink建模與仿真作者：王軍本文章為一次學(xué)習(xí)總結(jié)，發(fā)到網(wǎng)上供大家參考，希望大家轉(zhuǎn)載的時(shí)候不要匿名篡改，保持良好的學(xué)術(shù)作風(fēng)。在高中物理學(xué)習(xí)過程中，我們接觸了單擺.當(dāng)時(shí)的單擺定義是：細(xì)線一端固定在懸點(diǎn)，另一端系一個(gè)小球，如果細(xì)線的質(zhì)量與小球相比可以忽略，球的直接與線的長(zhǎng)度相比也可以忽略，這樣的裝置就叫做單擺.根據(jù)上面的定義，得出以下試驗(yàn)結(jié)論.(i)當(dāng)擺角很小時(shí)，周期與振幅無(wú)關(guān)；(2)周期與擺球質(zhì)量無(wú)關(guān)；(3)單擺振動(dòng)的周期與擺長(zhǎng)有關(guān)；單擺周期的平方與擺長(zhǎng)成正比.以上結(jié)論是在理想條件下得到的結(jié)論，現(xiàn)對(duì)這個(gè)理想條件下的單擺進(jìn)行分析與仿真，將仿真結(jié)果與以上結(jié)論進(jìn)

2、行對(duì)比驗(yàn)證.1理想模式下單擺的數(shù)學(xué)模型.首先根據(jù)理想條件，擺線質(zhì)量忽略不計(jì)，空氣阻力忽略不計(jì).設(shè)擺線長(zhǎng)度為1,擺球質(zhì)量為m,重力加速度為g,系統(tǒng)的初始時(shí)刻為t=0,在任意t)0時(shí)刻擺球的線速度為v(t),角速度為3(t),角位移(t),以單擺的固定位置為坐標(biāo)原點(diǎn)建立直角坐標(biāo)系，水平方向?yàn)閤軸方向.示意圖如下所示：在t時(shí)刻，擺錘所受切向力ft(t)是重力mgi其運(yùn)動(dòng)圓弧切線方向上的分力，即ft(t)?mgsin(t)完全理想條件下，根據(jù)牛頓第二運(yùn)動(dòng)定律，切向加速度為:a(t)?gsin(t)因此得到單擺的運(yùn)動(dòng)微分方程組：dv(t)dt?g sin(t)dJt)迎dtl使用歐拉算法求解：將dv(t

3、)?v(t?dt)?v(t)和d(t)?(t?dt)?(t)代入式(1)及式(2)中，并以仿真步進(jìn)量作為dt的近似，得到基于時(shí)間的遞推方程：v(t+尸v(t)+gsin(t)(3)(t+尸(t)?v(t)(4)l注：本遞推方程僅適合于擺角,也就是要求無(wú)論初始速度多少，擺角的最大幅度不能超過90度，如果超過90度比并且初始速度為0時(shí)放手小球會(huì)自由下落一段時(shí)間才能擺動(dòng)，本遞推方程無(wú)法描述.據(jù)此編寫仿真程序：在MALA瑜令窗口輸入以下命令：dt=;%仿真步進(jìn)T=16;%仿真時(shí)間長(zhǎng)度t=0:dt:T;%仿真計(jì)算時(shí)間序列g(shù)=;L=;th0=;%初始擺角設(shè)置，不能超過無(wú)/2,即要求球擺動(dòng)開始時(shí)繩子就要有拉

4、力，如果初始擺角超過無(wú)/2,則球會(huì)經(jīng)過一陣自由落體后才能進(jìn)行擺動(dòng)，上面的遞推方程不能滿足該情形.v0=0;%初始擺速設(shè)置v=zeros(size(t);%程序存儲(chǔ)變量預(yù)先初始化，可提高執(zhí)行速度th=zeros(size(t);v(1)=v0;th(1)=th0;fori=1:length(t)%仿真求解開始v(i+1)=v(i)+g*sin(th(i).*dt;th(i+1)=th(i)-1./L.*v(i).*dt;end%使用雙坐標(biāo)系統(tǒng)來作圖AX,B1,B2=plotyy(t,v(1:length(t),t,th(1:length(t),'plot');set(B1,

5、9;LineStyle','-');%設(shè)置圖線型set(B2,'LineStyle',':');set(get(AX(1),'Ylabel'),'String','線速度v(t)m/s');%作標(biāo)注set(get(AX(2),'Ylabel'),'String','角位移th(t)/rad');xlabel('時(shí)間t/s');legend(B1,'線速度v(t)',2);legend(B2,'角位移th(

6、t)',1);在以上假設(shè)條件下得到仿真圖形如下F些N至第母*o-S能A-回TUL-在其他條件不變的情況下，僅僅改變細(xì)線的長(zhǎng)度L=3,再次進(jìn)彳f仿真，仿真圖象如下：6 1口12-2-嫉速度皿對(duì)比兩幅圖象可以看出，在理想條件下，同樣擺角下，單擺的擺臂變化，影響單擺的最大線速度以及單擺的周期，當(dāng)擺臂增加時(shí)，最大線速度增加，同時(shí)單擺的周期也增加此結(jié)論正好與最初單擺理想條件下的試驗(yàn)結(jié)論一致實(shí)際情況中，當(dāng)擺幅很小(<50),單擺的擺動(dòng)可以看作是簡(jiǎn)諧運(yùn)動(dòng)，現(xiàn)在更改單擺的初始擺角為50情況，也就是令th0=,L依然為，仿真結(jié)果如下：J 2 2 Q o4J.統(tǒng)速度M】)卜，.與任穆世陽(yáng)1：2現(xiàn)實(shí)情

7、況下單擺的數(shù)學(xué)模型現(xiàn)實(shí)情況下，繩子的質(zhì)量，擺球的半徑，空氣的阻力等等都對(duì)單擺的擺動(dòng)有影響，這些影響的主要作用就是阻止單擺的擺動(dòng)，為簡(jiǎn)單起見，可設(shè)單擺在擺動(dòng)中受到阻力fz,顯然阻力與擺錘的運(yùn)動(dòng)速度有關(guān)，即阻力是單擺線速度的函數(shù)：fz=f(v),fz(t)?kv(t)上式中，k>0為阻力比例系數(shù)，式中的負(fù)號(hào)表示阻力方向與擺錘運(yùn)動(dòng)方向相反。切向加速度由切向合力ft?fz產(chǎn)生，根據(jù)牛頓第二運(yùn)動(dòng)定律，有a(t)?gsin(t)?kv(t)m因此得到修正后的單擺運(yùn)動(dòng)微分方程組：dv(t)kv(t)/?gsin(t)?(5)dtm3迎(6)dtl仍然使用歐拉算法求解：將dv(t)?v(t?dt)?v(

8、t)和d(t)?(t?dt)?(t)代入式(5)及式(6)中，并以仿真步進(jìn)量作為dt的近似，得到基于時(shí)間的遞推方程:v(t+尸v(t)+(gsin(t)?kv(t)v(t)A(t+尸(t)?-pA據(jù)此編寫仿真程序：在MALA瑜令窗口輸入以下命令：subplot(2,1,1)dt=;%仿真步進(jìn)T=16;%仿真時(shí)間長(zhǎng)度t=0:dt:T;%仿真計(jì)算時(shí)間序列g(shù)=；L=；m=8;k=3;th0=;%初始擺角設(shè)置，不能超過無(wú)/2v0=0;%初始擺速設(shè)置v=zeros(size(t);%程序存儲(chǔ)變量預(yù)先初始化，可提高執(zhí)行速度th=zeros(size(t);v(1)=v0;th(1)=th0;fori=1:

9、length(t)%仿真求解開始v(i+1)=v(i)+(g*sin(th(i)-k./m.*v(i).*dt;th(i+1)=th(i)-1./L.*v(i).*dt;end%使用雙坐標(biāo)系統(tǒng)來作圖AX,B1,B2=plotyy(t,v(1:length(t),t,th(1:length(t),'plot');set(B1,'LineStyle','-');%設(shè)置圖線型set(B2,'LineStyle',':');set(get(AX(1),'Ylabel'),'String',&

10、#39;線速度v(t)m/s');%作標(biāo)注set(get(AX(2),'Ylabel'),'String','角位移th(t)/rad');xlabel('時(shí)間t/s');legend(B1,'線速度v(t)',2);legend(B2,'角位移th(t)',1);其它條件不變，令阻力比例系數(shù)k=,將兩次仿真的圖形放在一起，結(jié)果如下：S ID 12Rt間gb e 10其它條件不變，令阻力比例系數(shù)k=0,將將兩次仿真的圖形放在一起，結(jié)果如下：K. -3OIJ，:HiJ.li-2B B ID 1

11、2'14 USat rar 1/5日方同1八由上圖可以看出，在理想態(tài)下，當(dāng)線速度最大的時(shí)候，角位移為0,當(dāng)線速度為0時(shí)角位移最大，兩個(gè)同時(shí)發(fā)生沒有滯后.非理想狀態(tài)下，角位移與線速度0值和最大值發(fā)生會(huì)有一定的延遲，這是由于各方面阻力綜合原因造成的.3理想條件下特殊擺動(dòng)的數(shù)學(xué)模型.這里所說的特殊擺動(dòng)，指的是在理想條件下，=的情況下，擺球有一個(gè)很大的初始速度，可以使擺球圍繞懸點(diǎn)不停的旋轉(zhuǎn)（垂直方向，重力能夠起到作用），此時(shí)已經(jīng)不能稱之為單擺了，但是在理想條件下，小球的轉(zhuǎn)動(dòng)依然符合一定的規(guī)律，先對(duì)其進(jìn)行數(shù)學(xué)物理分析，再進(jìn)行MATLAB真.理力學(xué)原理的分析同理想狀態(tài)下的單擺分析完全一樣，只是單擺

12、是來回?cái)[動(dòng)，而此時(shí)是圍繞懸點(diǎn)來回旋轉(zhuǎn)，擺動(dòng)的時(shí)候角位移是在一個(gè)范圍內(nèi)變化，而旋轉(zhuǎn)的情況下角位移是不停的增加的，如果時(shí)間無(wú)限長(zhǎng)，則角位移無(wú)線長(zhǎng)，圖象顯示就是一條上升的曲線，沒有很直觀的感受，因此我們只考慮線速度的變化規(guī)律.基于時(shí)間的遞推方程依然如下：v(t+尸v(t)+(gsin(t)?kv(t)mv(t).(t+尸(t)?Al只是在MATLABB真中不需要圖象顯示角位移的圖象，僅僅顯示線速度v的圖象就好.MATLA項(xiàng)真程序如下：dt=;%仿真步進(jìn)T=16;%仿真時(shí)間長(zhǎng)度t=0:dt:T;%仿真計(jì)算時(shí)間序列g(shù)=；L=；m=8;k=0;%空氣阻力比例系數(shù)th0=pi;%初始擺角設(shè)置v0=10;%初

13、始擺速設(shè)置v=zeros(size(t);%程序存儲(chǔ)變量預(yù)先初始化，可提高執(zhí)行速度th=zeros(size(t);v(1)=v0;th(1)=th0;fori=1:length(t)%仿真求解開始v(i+1)=v(i)+(g*sin(th(i)-k./m.*v(i).*dt;th(i+1)=th(i)-1./L.*v(i).*dt;endplot(t,v(1:length(t)xlabel('時(shí)間t/s');ylabel('線速度v(t)');仿真圖象如下所示：I由上圖可以看出，單擺圍繞懸點(diǎn)做圓周運(yùn)動(dòng)，速度在最高點(diǎn)的時(shí)候最低，在最低點(diǎn)的時(shí)候最大，現(xiàn)在驗(yàn)證最大速

14、度值：Vm*max帶入仿真初始值：v0=10 g=1=從而得出vmax12.6，與仿真結(jié)果明顯一致。由此還可以引出許多特殊情形，比如用硬制木棍代替細(xì)擺線，力學(xué)原理就發(fā)生了很大的改變，仿真的曲線也會(huì)有所不同.4使用Simulink重新仿真理想模式下的單擺運(yùn)動(dòng).本次仿真以單變量建立微分方程，用示波器觀察的變化情況.由之前的分析易知道有：?jiǎn)螖[回復(fù)力Fmg sin ,而單擺的切向加速度22d (L ) .d_2 L 2dtdtmL4dt2mg sindj dt2g .一 sinL上式為一個(gè)二階非線性常微分方程.為了便于比較，將所有參數(shù)設(shè)置的與之前一致也就是g=,L=，初始(0)1.5,0,則£

15、;6.5dtLSimulink模型模塊如下所示模塊搭建好后，進(jìn)行參數(shù)設(shè)置，積分器Integrator為初始角速度設(shè)置為0,積分器Integrator為初始=，三角函數(shù)sin為默認(rèn)設(shè)置，兩個(gè)增益.gk16.5,k2L1.5,SimulationConfigurationParameters設(shè)置L為起始時(shí)間0.終止時(shí)間16,采樣周期auto,最大間隔，最小，解法為ode45（4階5階龍格庫(kù)塔算法）.全部設(shè)置好后，啟動(dòng)Simulink開始仿真，兩個(gè)示波器的輸出圖象如下TimoofTsetQn 門 |> r ' " 4 r.,- I1.U A與之前MATLAB!序仿真結(jié)果完全一致-10o時(shí)同g與 口 石 gl庭蜀至此，本次針對(duì)單擺運(yùn)動(dòng)的仿真全部結(jié)束.明顯可以看出Simulink在系統(tǒng)結(jié)構(gòu)明確的情況下，只需進(jìn)行簡(jiǎn)單的參數(shù)設(shè)置，就可以觀測(cè)任意環(huán)節(jié)的變量狀態(tài)，十分簡(jiǎn)便.關(guān)于本次仿真的一些感受：1. 仿真的條件要求很嚴(yán)格，不同條件（理