MATLAB在力學(xué)、機(jī)械中的應(yīng)用舉例

MATLAB在力學(xué)、機(jī)械中的應(yīng)用舉例

7.1理論力學(xué)7.2材料力學(xué)7.3機(jī)械振動(dòng)7.1理論力學(xué)

【例7-1-1】給定由N個(gè)力(i=1，2，…，N)組成的平面任意力系，求其合力。

解:◆建模本程序可用來對(duì)平面任意力系作簡化，得出一個(gè)合力。求合力的過程可分成兩步。第一步：向任意給定點(diǎn)o簡化，得到一個(gè)主矢和一個(gè)主矩Mo，即式中，是作用點(diǎn)的矢徑；是o點(diǎn)的矢徑。第二步：將此主矢和主矩向t點(diǎn)轉(zhuǎn)移，使其力矩Mt為零，成為一個(gè)合力Ft。令注意向量可以用數(shù)組表示，用1×2數(shù)組來表示平面向量，此式可化為用MATLAB的右除符號(hào)，可以得到合力作用點(diǎn)t的坐標(biāo)rt為式中，rt和ro都是1×2的數(shù)組，可由此式解出rt?！鬗ATLAB程序

clear，N=input(′輸入力的數(shù)目N=′) %輸入力系中各力的數(shù)據(jù)

fori=1∶N

I,F(i,∶)=input(‘力F(i)的x，y兩個(gè)分量[Fx(i)，F(xiàn)y(i)]=’);r(I,∶)=input(‘力F(i)的一個(gè)作用點(diǎn)的坐標(biāo)r(i)=[rx,ry]=’);endro=input(′簡化中心ro的坐標(biāo)ro=［xo，yo］=′); %輸入簡化中心的數(shù)據(jù)

Fo=sum(F)， %求主矢Fo=［Fox，F(xiàn)oy］

fori=1∶N %計(jì)算各力對(duì)ro點(diǎn)的力矩

M(i)=F(i，2)*(r(I,1)-ro(1))-F(i，1)*(r(i，2)-ro(2));end Mo=sum(M) %相加求主矩

rt=Mo/［Fo(2);-Fo(1)］+ro %求合力作用點(diǎn)的坐標(biāo)◆程序運(yùn)行結(jié)果最后一條語句從一個(gè)方程要求出兩個(gè)未知數(shù)rt(1)和rt(2)，這是一個(gè)欠定方程，事實(shí)上合力作用線將通過平面上的無數(shù)點(diǎn)，程序中用矩陣右除的方法只能給出無數(shù)個(gè)解中的一個(gè)解，即rt-ro中有一個(gè)分量是零的那個(gè)解。運(yùn)行此程序，輸入

N=3，F(xiàn)(1，∶)=［2，3］，r(1，∶)=［-1，0］，

F(2，∶)=［-4，7］，r(2，∶)=［1，-2］，

F(3，∶)=［3，-4］，r(3，∶)=［1，2］，又設(shè)簡化中心的坐標(biāo)ro=［-1，-1］，答案為

Fo=［16］，Mo=-9(即x方向分力為1，y方向分力為6)

rt=［-2.5000-1.0000］(合力作用線通過的某一點(diǎn)坐標(biāo))

【例7-1-2】求圖7-1(a)所示桿系的支撐反力Na、Nb、Nc。設(shè)已知：G1=200;G2=100;L1=2;

;θ1=30°;θ2=45°。圖7-1桿系結(jié)構(gòu)及受力圖

解：

◆建模畫出桿1和桿2的受力圖(如圖7-1(b)、(c)所示)，列出方程。對(duì)桿1：

對(duì)桿2：

這是包含六個(gè)未知數(shù)Nax、Nay、Nbx、Nby、Ncx和Ncy的六個(gè)線性代數(shù)方程，要解這個(gè)方程組，通常是要尋找簡化的步驟，但用了MATLAB工具，也可以不簡化，把方程組寫成矩陣形式AX=B，用矩陣除法X=A＼B直接來解。在本題中，X和B都是6×1列向量，而A是6×6階方陣。在編寫程序時(shí)，盡量用文字變量，先輸入已知條件，在程序開始處給它們賦值，這樣得出的程序具有一定的普遍性，若要修改參數(shù)，只需修改頭幾行的數(shù)據(jù)即可?！鬗ATLAB程序

G1=200;G2=100;L1=2;L2=sqrt(2);%給原始參數(shù)賦值

theta1=30*pi/180;theta2=45*pi/180;%將度化為弧度

%設(shè)X=［Nax；Nay；Nbx；Nby；Ncx；Ncy］，則系數(shù)矩陣A、B可寫成下式

A=［1,0,0,0,1,0;0,1,0,0,0,1;0,0,0,0,-sin(theta1),

cos(theta1);...0,0,1,0,-1,0;0,0,0,1,0,-1;0,0,0,0,

sin(theta2),

cos(theta2)］

B=[0;G1;G1/2*cos(theta1);0;G2;-G2/2*cos(theta2)]

X=A＼B;%用左除求解線性方程組

disp(′Nax，Nay，Nbx，Nby，Ncx，Ncy′) %顯示結(jié)果

disp(X′)◆程序運(yùn)行結(jié)果

A=1.0000 0 001.00000 01.00000001.0000 0 0 00-0.50000.8660 0 01.00000-1.00000 0 00 1.00000-1.0000 0 0000.70710.7071B=[0200.000086.60250100.0000-35.3553]′

NaxNayNbxNbyNcxNcy

95.0962154.9038-95.0962145.0962-95.096245.0962圖7-2導(dǎo)彈攻擊目標(biāo)的運(yùn)動(dòng)

【例7-1-3】設(shè)導(dǎo)彈M速度分別為vm=1000m/s和800m/s，其速度向量始終對(duì)準(zhǔn)速度為vt=500m/s的直線飛行目標(biāo)T，發(fā)射點(diǎn)在目標(biāo)運(yùn)動(dòng)方向的左(4000m)前(3000m)方，如圖7-2所示，試求導(dǎo)彈軌跡及其加速度。

解：

◆建模在與目標(biāo)固連的等速直線運(yùn)動(dòng)坐標(biāo)(慣性坐標(biāo)系)中列寫動(dòng)點(diǎn)M的方程。因動(dòng)點(diǎn)坐標(biāo)與目標(biāo)T固連，牽連速度，動(dòng)點(diǎn)為M，它的絕對(duì)速度。由速度合成定理得相對(duì)速度，列出其在x、y兩方向的投影方程，得求其積分，即可求得其軌跡x=x(t)，y=y(t)?！鬗ATLAB程序

MATLAB數(shù)值積分要求把導(dǎo)數(shù)方程單獨(dú)列寫為一個(gè)函數(shù)程序，故其MATLAB程序由主程序和一個(gè)求導(dǎo)數(shù)的函數(shù)程序構(gòu)成。由于數(shù)值積分的步長是MATLAB按精度自動(dòng)選取的，其間隔可變，因此dt要用數(shù)組表示。主程序exn713: vt=input(′vt=′);vm=input(′vm=′); %輸入主程序及函數(shù)程序共用的參數(shù)

z0=input(′[x0;y0]=′);%輸入數(shù)值積分函數(shù)需要的參數(shù)

tspan=input(′tspan=［t0，tfinal］=′); %輸入數(shù)值積分函數(shù)需要的參數(shù) ［t，z］=ode23(′exn713f′，tspan，z0); %進(jìn)行數(shù)值積分

plot(z(∶，1)，z(∶，2)); %繪圖%在慣性坐標(biāo)中，M點(diǎn)位置的導(dǎo)數(shù)是相對(duì)速度，而其二次導(dǎo)數(shù)則為絕對(duì)加速度

dt=diff(t);Ldt=length(dt); %為了求導(dǎo)數(shù)，先求各時(shí)刻處t的增量

x=z(∶，1);y=z(∶，2); %把z寫成x、y兩個(gè)分量形式

vx=diff(z(∶，1))./dt;vy=diff(z(∶，2))./dt; %注意每差分一次序列長度減1wx=diff(vx)./dt(1∶Ldt-1);wy=diff(vy)./dt(1∶Ldt-1); %求二次導(dǎo)數(shù)［t(2∶Ldt)，x(2∶Ldt)，y(2∶Ldt)，wx，wy］%顯示數(shù)據(jù)下面是函數(shù)程序，寫成矩陣方程，存成一個(gè)文件exn713f。

functionzprime=exn713f(t，z，vt，vm) globalvtvm r=sqrt(z(1)^2+z(2)^2); zprime=［-vt-vm*z(1)/r;-vm*z(2)/r］◆程序運(yùn)行結(jié)果把上面兩個(gè)程序均存到MATLAB的搜索路徑上。運(yùn)行主程序并輸入以下參數(shù)：

vt=500;vm=1000

［x0;y0］=［3000;4000］

tspan=［t0，tfinal］=［0，4.5］得出圖形如圖7-3所示，數(shù)據(jù)如表7-1所示，為節(jié)省篇幅，表中省略了一些數(shù)據(jù)。注意：在給定tfinal時(shí)，必須使它小于遭遇點(diǎn)的值，否則數(shù)字積分會(huì)進(jìn)入死循環(huán)而得不出結(jié)果。讀者可以思考，能否修改程序，使它能自動(dòng)尋找到tfinal，并避免進(jìn)入死循環(huán)。不過這就不能用現(xiàn)成的ode23函數(shù)，而要自己編寫數(shù)值積分子程序才行。將vm換成800m/s，并相應(yīng)地把tfinal換成6，得到的軌跡位于圖7-3中原軌跡的左上方。圖7-3導(dǎo)彈跟蹤目標(biāo)時(shí)的相對(duì)軌跡

【例7-1-4】四連桿機(jī)構(gòu)如圖7-4所示，輸入桿L1的轉(zhuǎn)角θ1=ω1t，ω1=100rad/s，求輸出桿L3的轉(zhuǎn)角θ3隨時(shí)間的變化規(guī)律，并求其角速度和角加速度。解：

◆建模四連桿機(jī)構(gòu)的運(yùn)動(dòng)方程由X和Y方向的長度關(guān)系確定為

L1cosθ1+L2cosθ2-L3cosθ3-L0=0

(7-1)L1sinθ1+L2sinθ2-L3sinθ3=0

(7-2)從上述兩個(gè)方程中消去θ2，便可化成一個(gè)只包括θ1和θ3的方程，給定θ1，可求出滿足此方程的θ3。圖7-4四連桿機(jī)構(gòu)的幾何關(guān)系由（7-2）式得

sinθ2＝(L3sinθ3-L1sinθ1)/L2

(7-3)將(7-1)式中的cosθ2代以，得出(7-4)在θ1給定時(shí)，求能使f(θ3)=0的θ3值，然后，θ2就可由(7-3)式求得。為了求能使f(θ3)=0的θ3值，可調(diào)用MATLAB中的fzero函數(shù)。為此，要把f=f(θ3)單獨(dú)定義為一個(gè)MATLAB函數(shù)exn714f，在主程序中要調(diào)用它。為了把長度參數(shù)傳給子程序，在主程序和子程序中都加了全局變量語句(global)，但全局變量容易造成程序的混亂，要特別小心，在復(fù)雜的程序中應(yīng)盡量避免使用。求得θ1，θ2和θ3后，就不難根據(jù)桿1的角速度求出桿3的角速度，其方法有以下兩種:

(1)求瞬時(shí)速度，這是通常理論力學(xué)的解法，其依據(jù)就是桿2兩端點(diǎn)a和b的速度沿桿長方向的分量相等，通過三角關(guān)系，有從而由桿2兩端點(diǎn)a和b的速度沿桿長垂直方向的分量之差，可以求出桿2的角速度

(2)求運(yùn)動(dòng)全過程的角位置、角速度和角加速度曲線，這只有借助于計(jì)算工具才能做到，因?yàn)橛檬止に阋粋€(gè)點(diǎn)就不勝其煩，算幾十個(gè)點(diǎn)更是很難想象。而由MATLAB編程調(diào)用fzero函數(shù)時(shí)，要求給出一個(gè)近似猜測(cè)值，若連續(xù)算幾十點(diǎn)，前一個(gè)解就可作為后一個(gè)解的猜測(cè)值，從而帶來了方便。本書將提供exn714a和exn714b兩個(gè)程序來分別表述這兩種方法，它們所要調(diào)用的函數(shù)程序命名為exn714f?！鬗ATLAB程序

1.主程序exn714a

globalL0L1L2L3th1

L0=20;L1=8;L2=25;L3=20; %輸入基線及三根桿的長度L1、L2、L3

theta1=input(′當(dāng)前角theta1=′);

theta3=input(′對(duì)應(yīng)于theta1的theta3近似值=′);

th1=theta1;theta3=fzero(′exn714f′，theta3); %求當(dāng)前輸出角theta3

theta2=asin((L3*sin(theta3)-L1*sin(theta1))/L2);

w1=input(′w1=′);

w3=L1*w1*cos(pi/2-theta1+theta2)/(L3*cos(theta3-pi/2-theta2))

2.主程序exn714b

globalL0L1L2L3th1

L0=20;L1=8;L2=25;L3=20; %輸入基線及三根桿的長度L1、L2、L3

w1=input(′桿1角速度w1=′);

theta1=linspace(0，2*pi，181); %把桿1每圈分為180份，間隔2°

theta3=input(‘對(duì)應(yīng)于theta1最小值處的theta3(近似估計(jì)值)=’);

dt=2*pi/180/w;%桿1轉(zhuǎn)2°對(duì)應(yīng)的時(shí)間增量

th1=theta1(1);theta3(1)=fzero(′exn714f′，theta3);

%求初始輸出theta3

fori=2∶181th1=theta1(i);theta3(i)=fzero(′exn714f′，theta3(i-1)); %調(diào)用fzero函數(shù)逐次求theta3

endsubplot(1,2,1),plot(theta1,theta3),ylabel(′theta3′),grid

%畫曲線

w3=diff(theta3)/dt; %求桿3的角速度，注意求導(dǎo)數(shù)后數(shù)組長度小于1subplot(1,2,2),plot(theta1(2:length(theta1))，w3);grid %畫角速度曲線

3.子程序(函數(shù)程序)exn714f

functiony=exn714f(x)

globalL0L1L2L3th1

y=L1.*cos(th1)+L2*sqrt(1-(L3*sin(x)-L1*sin(th1)).

^2/L2/L2)…-L3*cos(x)-L0;在上述程序中注意th1是一個(gè)標(biāo)量，而theta1在exn714b中是一個(gè)數(shù)組，因?yàn)楹瘮?shù)exn714f中用到的是特定點(diǎn)的角度，是一個(gè)標(biāo)量，所以不能用thetal，引入了th1作為其當(dāng)前值?！舫绦蜻\(yùn)行結(jié)果在L0=20，L1=8，L2=25，L3=20(或15)的條件下運(yùn)行exn714b，根據(jù)提問，輸入w1=100，在thetal=0處，設(shè)theta3的近似初值為1弧度，所得的曲線如圖7-5(a)所示，相應(yīng)的角速度變化規(guī)律如圖7-5(b)所示。若運(yùn)行exn714a，其單點(diǎn)的數(shù)據(jù)與圖7-5一致，請(qǐng)讀者自行檢驗(yàn)。利用MATLAB，可以用動(dòng)畫來顯示四連桿的運(yùn)動(dòng)，運(yùn)行程序集中的exn714d就可以看到它的結(jié)果。有興趣的讀者可以打開此程序并讀懂它的原理，它并不難懂。圖7-5四連桿機(jī)構(gòu)的輸入輸出角位置關(guān)系和輸出角速度(a)輸入輸出角位置關(guān)系；(b)輸出角速度

【例7-1-5】對(duì)于拋射體，設(shè)空氣阻力的方向與拋射體質(zhì)心，速度向量相反，大小與拋射體質(zhì)心速度的平方成正比。拋射體受力圖如圖7-6所示?？荚嚳諝庾枇Γ?jì)算拋射體飛行的軌跡和距離。圖

7-6拋射體受力圖

解：◆建模在例6-2-1中，研究過不計(jì)空氣阻力的拋射體飛行軌跡，考慮空氣阻力后，程序要復(fù)雜一些。因?yàn)閤和y兩個(gè)方向的方程會(huì)通過空氣阻力互相耦合，必須聯(lián)立起來求解。根據(jù)受力圖7-6可列寫其運(yùn)動(dòng)方程:式中c為拋射體的空氣阻力系數(shù)。本來可以把兩個(gè)速度導(dǎo)數(shù)的方程分別聯(lián)立起來求解，先求出速度，再積分而求出位置。但在MATLAB中，階次高并不造成困難，分成兩步反而增加編程的工作量，所以往往選擇一次解出這個(gè)四階方程。這里設(shè)了一個(gè)四維列向量

來表示四個(gè)變量，方程組可寫成矩陣形式:為了調(diào)用MATLAB的數(shù)值積分法函數(shù)，要把其運(yùn)動(dòng)方程組寫成一個(gè)函數(shù)文件，取文件名為exn715f。該運(yùn)動(dòng)方程組是一個(gè)四行的向量方程組，表明系統(tǒng)為四階?！鬗ATLAB程序

1.函數(shù)程序exn715f

functionrdot=exn715f(t，r)

c=0.01;g=9.81;m=1; %給出空氣阻力系數(shù)及重力加速度(m/s^2)

vm=sqrt(r(3)^2+r(4)^2);%速度大小

rdot=[r(3);r(4);-c*vm*r(3)/m;-c*vm*r(4)/m-g]; %運(yùn)動(dòng)方程

2.主程序exn715

clear;y0=0;x0=0; %初始位置

vMag=input(′輸入初始速度(m/s):′); %輸入初始速度

vDir=input(′輸入初速方向(度):′);

tf=input(′輸入飛行時(shí)間（s）:′); %輸入飛行時(shí)間

vx0=vMag*cos(vDir*(pi/180)); %計(jì)算x，y方向的初始速度

vy0=vMag*sin(vDir*(pi/180));

0=［0;0;vx0;vy0］;

［t，r］=ode45(′exn715f′，［0，tf］，r0)，

%數(shù)值積分(調(diào)用函數(shù)程序exn715f)

plot(r(∶，1)，r(∶，2))，holdon %計(jì)算軌跡

%ode45規(guī)定返回的結(jié)果中：t是列向量，各時(shí)刻的r成為四列向量

%注意下一語句的意義:找y<0的下標(biāo)所對(duì)應(yīng)的x的最小值，以粗略計(jì)算射程

xmax=min(r(find(r(∶，2)<0)，1))

plot(［0，150］，［0，0］) %畫出x坐標(biāo)線◆程序運(yùn)行結(jié)果輸入初始速度(m/s):60輸入初速方向(度):45輸入飛行時(shí)間（s）:6.2 t

x y vx vy 0

0 0 42.426442.4264 0.3002

11.7236

11.3046

36.0492

33.3239 1.164637.864032.179325.785716.5362

…

…

…

…

…5.6509111.85784.571410.0477-22.17166.2000117.0149-8.21908.7531-24.3438換新的參數(shù)：輸入初始速度(m/s):60輸入初速方向(度):35輸入飛行時(shí)間（s）:6得到近似射程xmax=123.1946。其軌跡如圖7-7所示。讀者可思考如何能求出射程的精確值。圖

7-7考慮空氣阻力后的拋射體軌跡

【例7-1-6】給定半徑r=0.1m、重量Q=2kg的均質(zhì)保齡球，球的初始速度v0=3m/s，初始角速度ω0=0，地面的摩擦系數(shù)f=0.05，問經(jīng)過多少時(shí)間后，球?qū)o滑動(dòng)地滾動(dòng)，求此時(shí)球心的速度。

解：

◆建模保齡球受力情況如圖7-8所示，接觸面之間打滑時(shí)，摩擦力使圓柱質(zhì)心減速，而使其轉(zhuǎn)動(dòng)加速。當(dāng)圓柱觸地點(diǎn)C的線速度達(dá)到0，即v=ω*r時(shí)，進(jìn)入純滾動(dòng)狀態(tài)。圖7-8圓柱運(yùn)動(dòng)受力圖已知球的質(zhì)量為，轉(zhuǎn)動(dòng)慣量為，可列出動(dòng)力學(xué)方程積分可得：

v=v0-fgt

(7-7)

(7-８)將(7-7)式和(7-8)式聯(lián)立,可求得滿足v=rω的時(shí)刻為。◆MATLAB程序

r=0.1;Q=2;g=9.81; %輸入常數(shù)

f=0.05;v0=3;w0=0;

J=Q*r^2/2.5/g;F=f*Q; %要計(jì)算的常數(shù)

wdot=F*r/J;

%繞質(zhì)心轉(zhuǎn)動(dòng)加速度方程

vdot=-F/(Q/g);

%質(zhì)心線加速度方程

t1=(v0-w0*r)/(wdot*r-vdot)

%求t1的方程

v=v0+vdot*t1

%求v的方程

◆程序運(yùn)行結(jié)果運(yùn)行此程序的結(jié)果為:t1=0.8737，v=2.1429即經(jīng)過0.8737s后，保齡球進(jìn)入純滾動(dòng)狀態(tài)，此時(shí)質(zhì)心速度為2.1429m/s。7.2材料力學(xué)

【例7-2-1】拉壓桿系的靜不定問題。由n根桿組成的桁架結(jié)構(gòu)如圖7-9所示，受力P的作用，各桿截面面積分別為Ai，材料彈性模量為E，求各桿的軸力Ni及節(jié)點(diǎn)A的位移。圖

7-9任意靜不定桿系受力圖

解：

◆建模先列寫具有普遍意義的方程，設(shè)各桿均受拉力，A點(diǎn)因各桿變形而引起的x方向位移為Δx，y方向位移為Δy，由幾何關(guān)系得變形方程即其中，為桿i的剛度系數(shù)。再加上兩個(gè)力平衡方程

共有n+2個(gè)方程，其中包含n個(gè)未知力和兩個(gè)待求位移Δx和Δy，方程組可解。因?yàn)檫@又是一個(gè)線性方程組,可寫成D*X=B的標(biāo)準(zhǔn)形式，所以可由MATLAB的矩陣除法X=D＼B解出。算例:設(shè)三根桿組成的桁架如圖7-10所示，掛一重物P=3000N，設(shè)L=2m，各桿的截面積分別為A1=200×10-6m2，A2=300×10-6m2，A3=400×10-6m2，材料的彈性模量E=200×109N/m2，求各桿受力的大小。圖

7-10靜不定三桿受力圖

此時(shí)應(yīng)有五個(gè)方程如下:力平衡：

-N1cosα1-N2-N3cosα3=0 N1sinα1-N3sinα3=０位移協(xié)調(diào)：

N1/K1=Δxcosα1+Δysinα1 N2/K2=Δy N3/K3=Δxcosα3-Δysinα3設(shè)X=［N1;N2;N3;Δx;Δy］，把上述五個(gè)線性方程組列成D*X=B的矩陣形式?！鬗ATLAB程序

P=3000;E=200e9;L=2

A1=200e-6;A2=300e-6;A3=400e-6;

a=pi/3;a2=pi/2;a3=3*pi/4;

L1=L/sin(a1);L2=L/sin(a2);L3=L/sin(a3);%計(jì)算桿長

K1=E*A1/L1;K2=E*A2/L2;K3=E*A3/L3; %計(jì)算剛度系數(shù)

%為避免語句太長，給系數(shù)矩陣按行賦值

D(1，∶)=［cos(a1)，cos(a2)，cos(a3)，0，0］;D(2，∶)=［sin(a1)，sin(a2)，sin(a3)，0，0］;D(3，∶)=［1/K1，0，0，-cos(a1)，-sin(a1)］;D(4，∶)=［0，1/K2，0，-cos(a2)，-sin(a2)］;D(5，∶)=［0，0，1/K3，-cos(a3)，-sin(a3)］;B=［P;0;0;0;0］;formatlong，X=D＼B %求解線性方程組，用長格式顯示結(jié)果◆程序執(zhí)行結(jié)果執(zhí)行此程序，用formatlong顯示的結(jié)果為

若用普通格式顯示，將得出Δy=0.0000，實(shí)際上Δy不是零，這可從N2不等于零推知。在程序中用一個(gè)矩陣顯示數(shù)值相差很大的元素時(shí)，就得采用formatlong，以免丟失小的量。也可要求系統(tǒng)單獨(dú)顯示此元素的值，例如鍵入x(5)，系統(tǒng)將給出ans=1.970475034321116e-005。讀者還可改變幾根桿的剛度系數(shù)，看它們?nèi)绾斡绊懜鳁U受力的分布。

【例7-2-2】長為L的懸臂梁如圖7-11所示，左端固定，在離固定端L1處施加力P，求它的轉(zhuǎn)角和撓度。已知梁的彈性模量E=200×109N/m2和截面慣性矩I=2×10-5m4

。圖

7-11懸臂梁受力圖

解：

◆建模材料力學(xué)中從彎矩求轉(zhuǎn)角要經(jīng)過一次積分，而從轉(zhuǎn)角求撓度又要經(jīng)過一次積分，這不僅很麻煩而且容易出錯(cuò)。在MATLAB中，可用cumsum函數(shù)或更精確的cumtrapz函數(shù)作近似的不定積分，只要x取得足夠密，其結(jié)果是相當(dāng)精確的，且程序非常簡單。本題采用cumsum函數(shù)。解題的關(guān)鍵還在于正確地列寫彎矩方程，請(qǐng)讀者注意程序中的這部分。本題的彎矩方程為轉(zhuǎn)角

撓度◆MATLAB程序

clear

L=2;P=2000;L1=1.5; %給出已知常數(shù)

E=200e9;I=2e-5;

x=linspace(0，L，101);dx=L/100; %將x分成100段，步長為L/100

n1=L1/dx+1; %確定x=L1處對(duì)應(yīng)的下標(biāo)

M1=-P*(L1-x(1:n1));%第一段彎矩賦值

M2=zeros(1，101-n1);%第二段彎矩賦值(全為零)

M=［M1，M2］;%全梁的彎矩A=cumsum(M)*dx/(E*I); %對(duì)彎矩積分求轉(zhuǎn)角

Y=cumsum(A)*dx; %對(duì)轉(zhuǎn)角積分求撓度

subplot(3，1，1)，plot(x，M)，grid%繪彎矩圖

subplot(3，1，2)，plot(x，A)，grid%繪彎矩圖

subplot(3，1，3)，plot(x，Y)，grid%繪彎矩圖◆程序運(yùn)行結(jié)果運(yùn)行程序所得的結(jié)果如圖7-12所示。注意幾根曲線之間的積分關(guān)系。本題之所以簡單，是因?yàn)樵趚=0處，轉(zhuǎn)角和撓度都為零，因此兩次積分的積分常數(shù)恰好都為零。如果它們不為零，程序中就得有確定積分常數(shù)的語句，這可在例7-2-3中看到。圖

7-12懸臂梁彎矩(M)、轉(zhuǎn)角（A）和撓度（Y）曲線

【例7-2-3】簡支梁左半部分受均勻分布載荷q作用，右邊L/4處受集中力偶M0作用(如圖7-13所示)，求其彎矩、轉(zhuǎn)角和撓度。設(shè)L=2m，q=1000N/m，M0＝900N·m，E=200×109N/m2，I=2×10-6m4。圖

7-13簡支梁受力圖

解：

◆建模此題解法基本上與例7-2-2相同，主要差別是要處理積分常數(shù)問題。支撐反力Na和Nb可由平衡方程求得，設(shè)，則各段彎矩方程為:對(duì)M/EI作積分，得轉(zhuǎn)角A，再作一次積分，得到撓度Y，每次積分都要出現(xiàn)一個(gè)待定積分常數(shù)此處設(shè)A0(x)=cumtrapz(M)*dx/EI。此處設(shè)Y0(x)=cumtrapz(A0)*dx。兩個(gè)待定積分常數(shù)Ca和Cy可由邊界條件Y(0)=0及Y(L)=0確定:Y(0)=Y0(0)+Cy=0Y(L)=Y0(L)+Ca·L+Cy=0于是可得即◆MATLAB程序

%輸入已知參數(shù)L，q，M0，E，I后，先求兩鉸鏈的支撐反力Na和Nb

L=2;q=1000;M0=900;E=200e9;I=2e-6;

Na=(3*q*L^2/8-M0)/L;Nb=(q*L^2/8+M0)/L;%求支撐反力

x=linspace(0，L，101);dx=L/100;%將x分為100小段

M1=Na*x(1∶51)-q*x(1∶51).^2/2; %分三段用數(shù)組列出M的表達(dá)式

M2=Nb*(L-x(52∶76))-M0;

M3=Nb*(L-x(77∶101));M=［M1，M2，M3］; %列寫完整的M數(shù)組

A0=cumtrapz(M)*dx/(E*J); %由M積分求轉(zhuǎn)角（未計(jì)積分常數(shù)）

Y0=cumtrapz(A0)*dx; %由轉(zhuǎn)角積分求撓度（未計(jì)積分常數(shù)）

C=［0，1;L，1］＼［-Y0(1);-Y0(101)］; %由邊界條件求積分常數(shù)Ca，Cy

Ca=C(1)，Cy=C(2)，

A=A0+Ca;Y=Y0+Ca*x+Cy;%求出轉(zhuǎn)角與撓度的完整值

subplot(3，1，1)，plot(x，M)，grid %繪圖

subplot(3，1，2)，plot(x，A)，grid

subplot(3，1，3)，plot(x，Y)，grid◆程序運(yùn)行結(jié)果執(zhí)行本程序的結(jié)果如圖7-14所示。梯形積分累加函數(shù)cumtrapz與定積分函數(shù)trapz的不同在于cumtrapz類似于不定積分，逐點(diǎn)給出積分的值，因而得出一個(gè)數(shù)列，而trapz只給出積分到終點(diǎn)的一個(gè)值。這些函數(shù)都假定步長為1，因此累加的值必須乘以dx才與積分等價(jià)。用A=cumtrapz(M)來求面積，長度M為101，只能形成100個(gè)A。而cumsum則是把101個(gè)點(diǎn)逐個(gè)相加，相當(dāng)于多算了一個(gè)點(diǎn)。準(zhǔn)確地說，可以推導(dǎo)出

cumtrapz(M)=cumsum(M)-M(1)/2-M/2實(shí)際上只要點(diǎn)取得足夠多，直接用cumsum(M)代替cumtrapz(M)，在工程上也是可以接受的。圖7-14例7-2-3的彎矩、轉(zhuǎn)角和撓度曲線(a)彎矩曲線;(b)轉(zhuǎn)角曲線;(c)撓度曲線

【例7-2-4】拉彎合成部件的截面設(shè)計(jì)。這一設(shè)計(jì)計(jì)算將歸結(jié)為解一個(gè)三次代數(shù)方程,過去要用試湊法反復(fù)運(yùn)算,本例顯示了用MATLAB求解的簡潔。鉆床立柱如圖7-15所示。設(shè)P＝15kN，許用拉應(yīng)力[σ]=35MPa，鉆頭軸與立柱軸距離為0.4m，試求立柱直徑。

圖7-15鉆床受力圖

解：

◆建模立柱受到拉力P和彎矩Pl作用，兩者產(chǎn)生的拉應(yīng)力之和為最大拉應(yīng)力，令它小于［σ］，即把代入上式后，得到求直徑d的方程這個(gè)三次代數(shù)方程可用MATLAB多項(xiàng)式求根的roots函數(shù)求解。◆MATLAB程序

P=input(′P=′)，l=input(′l=′)，%輸入力和偏心距

asigma=input(′［σ］=′)，%輸入許用拉應(yīng)力

a=［asigma*pi/32，0，-P/8，-P*l］; %求三次代數(shù)方程的系數(shù)向量

r=roots(a);

%求代數(shù)方程的根

d=r(find(imag(r)==0))

%只取實(shí)根◆程序運(yùn)行結(jié)果運(yùn)行此程序，按提示輸入以下條件

P=15000，l=0.4，［σ］=35e6得到的解為

d=0.1219m7.3機(jī)械振動(dòng)

【例7-3-1】分析單自由度阻尼系統(tǒng)的阻尼系數(shù)對(duì)其固有振動(dòng)模態(tài)的影響。

解：

◆建模單自由度阻尼系統(tǒng)振動(dòng)方程如下(7-9)其中m為物體質(zhì)量，k為彈簧剛度，c為阻尼常數(shù)，f為所加的外力。為了研究它的固有振動(dòng),設(shè)外加力f=0。將方程整理后寫成(7-10)式中，為固有頻率，為阻尼系數(shù)?，F(xiàn)取ζ=0.1~1，ωn=10，初始條件分別為x0=1、v0=0及x0=0、v0=1，進(jìn)行討論。根據(jù)微分方程理論，這樣一個(gè)常系數(shù)二階方程，其通解的形式為，將它代入方程，得知p1、p2是特征方程的兩個(gè)根，而C1、C2則由初始條件決定。很多教材都用傳統(tǒng)的解析法解這個(gè)題目，在ζ<1時(shí)，p1、p2是一對(duì)共軛復(fù)根，p1=-ζωn+jωd，p2=-ζωn-jωd，其中，此時(shí)，解就可寫成正余弦函數(shù)的形式，常數(shù)C1、C2就轉(zhuǎn)化為A和φ:其中，。用MATLAB作為計(jì)算工具對(duì)這些公式進(jìn)行計(jì)算和繪圖，比用手工計(jì)算和繪圖方便得多。先設(shè)x0=1，v0=1，時(shí)間區(qū)間為t=0～2s，ζ按步長0.1由0增加到1，可以得到如下的MATLAB程序?！鬗ATLAB程序(exn731a)

clear，wn=10;tf=2;x0=1;v0=0;

forj=1∶10zeta(j)=0.1*j; %設(shè)定不同的ζwd(j)=wn*sqrt(1-zeta(j)^2); %求ωda=sqrt((wn*x0*zeta(j)+v0)^2+(x0*wd(j))^2)/wd(j); %求振幅Aphi=atan2(wd(j)*x0,v0+zeta(j)*wn*x0); %用atan2是為了求四象限相角

t=0∶tf/1000∶tf; %設(shè)定自變量數(shù)組x(j，∶)=a*exp(-zata(j)*wn*t).*sin(wd(j)*t+phi); %求過渡過程

end

plot(t,x(1,∶),t,x(2,∶),t,x(3,∶),t,x(4,∶),t,x(5,∶),%繪圖

t，x(6,∶)，t，x(7,∶)，t，x(8,∶)，t，x(9,∶)，t，x(10,∶))grid，figure，mesh(x) %畫出三維圖形◆程序運(yùn)行結(jié)果執(zhí)行此程序即可得到圖7-16所示結(jié)果。改變初始條件為x0=0，v0=1(程序exn731b)，可得到圖7-17所示結(jié)果。實(shí)際上后一組曲線就是系統(tǒng)的脈沖過渡函數(shù)。因?yàn)槊}沖函數(shù)的幅度是無窮大，而持續(xù)時(shí)間卻是無限小，其面積為一個(gè)單位，所以脈沖激勵(lì)的最后效果(在t=+eps處)可形成一個(gè)單位的初速v0，由它產(chǎn)生的波形就是脈沖過渡函數(shù)的波形。圖

7-16初值x0=1、v0=0時(shí)的振動(dòng)波形

圖

7-17初值x0=0、v0=1時(shí)的振動(dòng)波形

鍵入mesh(t，zeta，x)，可以得到以ζ為參數(shù)的脈沖響應(yīng)。三維圖形如圖7-18所示，從中可以更形象地看出ζ對(duì)固有振動(dòng)模態(tài)的影響，因?yàn)闆]有彩色，所以視覺效果要差一些，讀者在計(jì)算機(jī)屏幕上看要好得多，并且可以再鍵入rotate3d命令，以便用鼠標(biāo)拖動(dòng)三維圖形旋轉(zhuǎn)，獲得更清晰的概念。圖

7-18不同ζ對(duì)固有振動(dòng)模態(tài)影響的三維圖

上述公式之所以繁瑣，是因?yàn)橐苊鈴?fù)數(shù)運(yùn)算。而MATLAB本身就具備復(fù)數(shù)運(yùn)算的功能，可使求解變得更為簡單。設(shè)原始微分方程為

先用roots函數(shù)求p1、p2，語句為:p=roots(［a2，a1，1］)，然后不必管它們是否為復(fù)數(shù)，把對(duì)t求導(dǎo)，得代入初始條件可得x0=C1+C2，v0=C1p1+C2p2將這兩個(gè)線性方程聯(lián)立，解出C1、C2。

這樣，p1、p2和C1、C2都已求出，只要給出t數(shù)組，就求出了x。要注意的是:雖然我們知道，在系統(tǒng)的實(shí)際運(yùn)動(dòng)中，x必然是實(shí)數(shù)，但在復(fù)數(shù)運(yùn)算中，由于計(jì)算的誤差，難免會(huì)出現(xiàn)微小的虛部，使x變成復(fù)數(shù)。這會(huì)使繪圖語句無法執(zhí)行，因此要用real語句取出實(shí)部，才能繪圖。假如其他參數(shù)不變，要求出ζ=0.3的情況下此系統(tǒng)的脈沖響應(yīng),則可編出程序exn731c如下，其核心語句只有中間的三句，語句簡單多了。

wn=10;tf=2;x0=1;v0=0;zeta=0.3;t=0∶tf/1000∶tf; %輸入?yún)?shù)和自變量數(shù)組

p=roots(［1，2*zeta*wn，wn^2］) %求特征方程的根

C=inv(［1，1;p(1)，p(2)］)*［x0;v0］

%求各暫態(tài)分量常數(shù)

x=real(C(1)*exp(p(1)*t)+C(2)*exp(p(2)*t)); %用real消除虛數(shù)

plot(t，x)，gridon %繪圖

【例7-3-2】設(shè)單自由度阻尼系統(tǒng)的質(zhì)量m=1kg，彈簧剛度系數(shù)K=100N/m，速度阻尼系數(shù)c=4N·s/m，求它在如下外力作用下的強(qiáng)迫振動(dòng)，并畫出t≤1.2s的波形。

解:

◆建模首先求出系統(tǒng)的脈沖過渡函數(shù)h(t)，則強(qiáng)迫振動(dòng)的波形就等于h(t)和外加力f(t)作卷積的結(jié)果，即在MATLAB中，h(t)和f(t)都可用數(shù)組h和f表示，其取樣間隔dt應(yīng)相同，便有x=conv(h，f)*dt卷積計(jì)算很繁，通常講振動(dòng)時(shí)只討論一些標(biāo)準(zhǔn)的有解析表達(dá)式的激勵(lì)信號(hào)，如方波、正弦波等，而用了MATLAB就不必受限制。在本題中脈沖過渡函數(shù)用極點(diǎn)留數(shù)函數(shù)residue求得，然后用卷積函數(shù)conv根據(jù)輸入函數(shù)和脈沖過渡函數(shù)求輸出。極點(diǎn)留數(shù)法是以拉普拉斯變換為基礎(chǔ)的，設(shè)m=1，對(duì)振動(dòng)方程(7-9)的兩端作拉氏變換，得設(shè)輸入為單位脈沖，其拉普拉斯變換F(s)=1，把X(s)用極點(diǎn)留數(shù)表示，有p1、p2是X的兩個(gè)極點(diǎn)，r1、r2則是對(duì)應(yīng)的留數(shù)，X(s)的拉普拉斯反變換為除了將C換成了r，這個(gè)公式和上題完全相仿。好處是MATLAB中有專門的函數(shù)來同時(shí)求出p和r，其調(diào)用方式為 ［r，p］=residue(b，a)其中b、a分別是X(s)分子、分母多項(xiàng)式的系數(shù)向量。于是可編出本題的程序exn732?！鬗ATLAB程序

m=1;c=4;K=100;dt=0.015; %輸入給定的參數(shù)

w0=sqrt(K/m); %求系統(tǒng)固有頻率

zeta=c/sqrt(m*K)/2; %求系統(tǒng)固有阻尼系數(shù)

a=［1，2*zeta*w0，w0^2］;b=1; %求分母、分子的系數(shù)［r，p］=residue(b，a); %求極點(diǎn)、留數(shù)

t=0∶dt:1.2;

h=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t); %求出系統(tǒng)的脈沖響應(yīng)

f=［1∶10，10*ones(1，70)］;

%給出外加力的采樣值

x=conv(h，f)*dt; %把脈沖響應(yīng)和外加力作卷積

plot(t(1∶80)，x(1∶80))%繪圖

v1=diff(x)/dt; %求導(dǎo)得出速度，注意求導(dǎo)后數(shù)組長度少1

[t(1∶80)′，f(1∶80)′，x(1∶80)′，［0，v1(1∶79)］′] %列出結(jié)果◆程序運(yùn)行結(jié)果執(zhí)行此程序的結(jié)果如圖7-19所示。其大量的數(shù)值結(jié)果予以刪略。讀者不妨把程序中的h改用例7-3-