MATLAB在力學、機械中的應用舉例

MATLAB在力學、機械中的應用舉例

7.1理論力學7.2材料力學7.3機械振動7.1理論力學

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

解:◆建模本程序可用來對平面任意力系作簡化，得出一個合力。求合力的過程可分成兩步。第一步：向任意給定點o簡化，得到一個主矢和一個主矩Mo，即式中，是作用點的矢徑；是o點的矢徑。第二步：將此主矢和主矩向t點轉(zhuǎn)移，使其力矩Mt為零，成為一個合力Ft。令注意向量可以用數(shù)組表示，用1×2數(shù)組來表示平面向量，此式可化為用MATLAB的右除符號，可以得到合力作用點t的坐標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兩個分量[Fx(i)，F(xiàn)y(i)]=’);r(I,∶)=input(‘力F(i)的一個作用點的坐標r(i)=[rx,ry]=’);endro=input(′簡化中心ro的坐標ro=［xo，yo］=′); %輸入簡化中心的數(shù)據(jù)

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

fori=1∶N %計算各力對ro點的力矩

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 %求合力作用點的坐標◆程序運行結(jié)果最后一條語句從一個方程要求出兩個未知數(shù)rt(1)和rt(2)，這是一個欠定方程，事實上合力作用線將通過平面上的無數(shù)點，程序中用矩陣右除的方法只能給出無數(shù)個解中的一個解，即rt-ro中有一個分量是零的那個解。運行此程序，輸入

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］，又設簡化中心的坐標ro=［-1，-1］，答案為

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

rt=［-2.5000-1.0000］(合力作用線通過的某一點坐標)

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

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

解：

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

對桿2：

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

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

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

%設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′)◆程序運行結(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導彈攻擊目標的運動

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

解：

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

MATLAB數(shù)值積分要求把導數(shù)方程單獨列寫為一個函數(shù)程序，故其MATLAB程序由主程序和一個求導數(shù)的函數(shù)程序構(gòu)成。由于數(shù)值積分的步長是MATLAB按精度自動選取的，其間隔可變，因此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); %進行數(shù)值積分

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

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

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

vx=diff(z(∶，1))./dt;vy=diff(z(∶，2))./dt; %注意每差分一次序列長度減1wx=diff(vx)./dt(1∶Ldt-1);wy=diff(vy)./dt(1∶Ldt-1); %求二次導數(shù)［t(2∶Ldt)，x(2∶Ldt)，y(2∶Ldt)，wx，wy］%顯示數(shù)據(jù)下面是函數(shù)程序，寫成矩陣方程，存成一個文件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］◆程序運行結(jié)果把上面兩個程序均存到MATLAB的搜索路徑上。運行主程序并輸入以下參數(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ù)字積分會進入死循環(huán)而得不出結(jié)果。讀者可以思考，能否修改程序，使它能自動尋找到tfinal，并避免進入死循環(huán)。不過這就不能用現(xiàn)成的ode23函數(shù)，而要自己編寫數(shù)值積分子程序才行。將vm換成800m/s，并相應地把tfinal換成6，得到的軌跡位于圖7-3中原軌跡的左上方。圖7-3導彈跟蹤目標時的相對軌跡

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

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

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

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

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

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

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

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

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

1.主程序exn714a

globalL0L1L2L3th1

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

theta1=input(′當前角theta1=′);

theta3=input(′對應于theta1的theta3近似值=′);

th1=theta1;theta3=fzero(′exn714f′，theta3); %求當前輸出角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(‘對應于theta1最小值處的theta3(近似估計值)=’);

dt=2*pi/180/w;%桿1轉(zhuǎn)2°對應的時間增量

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的角速度，注意求導數(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是一個標量，而theta1在exn714b中是一個數(shù)組，因為函數(shù)exn714f中用到的是特定點的角度，是一個標量，所以不能用thetal，引入了th1作為其當前值。◆程序運行結(jié)果在L0=20，L1=8，L2=25，L3=20(或15)的條件下運行exn714b，根據(jù)提問，輸入w1=100，在thetal=0處，設theta3的近似初值為1弧度，所得的曲線如圖7-5(a)所示，相應的角速度變化規(guī)律如圖7-5(b)所示。若運行exn714a，其單點的數(shù)據(jù)與圖7-5一致，請讀者自行檢驗。利用MATLAB，可以用動畫來顯示四連桿的運動，運行程序集中的exn714d就可以看到它的結(jié)果。有興趣的讀者可以打開此程序并讀懂它的原理，它并不難懂。圖7-5四連桿機構(gòu)的輸入輸出角位置關(guān)系和輸出角速度(a)輸入輸出角位置關(guān)系；(b)輸出角速度

【例7-1-5】對于拋射體，設空氣阻力的方向與拋射體質(zhì)心，速度向量相反，大小與拋射體質(zhì)心速度的平方成正比。拋射體受力圖如圖7-6所示。考試空氣阻力，計算拋射體飛行的軌跡和距離。圖

7-6拋射體受力圖

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

來表示四個變量，方程組可寫成矩陣形式:為了調(diào)用MATLAB的數(shù)值積分法函數(shù)，要把其運動方程組寫成一個函數(shù)文件，取文件名為exn715f。該運動方程組是一個四行的向量方程組，表明系統(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]; %運動方程

2.主程序exn715

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

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

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

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

vx0=vMag*cos(vDir*(pi/180)); %計算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 %計算軌跡

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

%注意下一語句的意義:找y<0的下標所對應的x的最小值，以粗略計算射程

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

plot(［0，150］，［0，0］) %畫出x坐標線◆程序運行結(jié)果輸入初始速度(m/s):60輸入初速方向(度):45輸入飛行時間（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輸入飛行時間（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)過多少時間后，球?qū)o滑動地滾動，求此時球心的速度。

解：

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

v=v0-fgt

(7-7)

(7-８)將(7-7)式和(7-8)式聯(lián)立,可求得滿足v=rω的時刻為?！鬗ATLAB程序

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; %要計算的常數(shù)

wdot=F*r/J;

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

vdot=-F/(Q/g);

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

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

%求t1的方程

v=v0+vdot*t1

%求v的方程

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

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

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

解：

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

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

7-10靜不定三桿受力圖

此時應有五個方程如下:力平衡：

-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設X=［N1;N2;N3;Δx;Δy］，把上述五個線性方程組列成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);%計算桿長

K1=E*A1/L1;K2=E*A2/L2;K3=E*A3/L3; %計算剛度系數(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，實際上Δy不是零，這可從N2不等于零推知。在程序中用一個矩陣顯示數(shù)值相差很大的元素時，就得采用formatlong，以免丟失小的量。也可要求系統(tǒng)單獨顯示此元素的值，例如鍵入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懸臂梁受力圖

解：

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

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

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

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

Y=cumsum(A)*dx; %對轉(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%繪彎矩圖◆程序運行結(jié)果運行程序所得的結(jié)果如圖7-12所示。注意幾根曲線之間的積分關(guān)系。本題之所以簡單，是因為在x=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)角和撓度。設L=2m，q=1000N/m，M0＝900N·m，E=200×109N/m2，I=2×10-6m4。圖

7-13簡支梁受力圖

解：

◆建模此題解法基本上與例7-2-2相同，主要差別是要處理積分常數(shù)問題。支撐反力Na和Nb可由平衡方程求得，設，則各段彎矩方程為:對M/EI作積分，得轉(zhuǎn)角A，再作一次積分，得到撓度Y，每次積分都要出現(xiàn)一個待定積分常數(shù)此處設A0(x)=cumtrapz(M)*dx/EI。此處設Y0(x)=cumtrapz(A0)*dx。兩個待定積分常數(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的表達式

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)角（未計積分常數(shù)）

Y0=cumtrapz(A0)*dx; %由轉(zhuǎn)角積分求撓度（未計積分常數(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◆程序運行結(jié)果執(zhí)行本程序的結(jié)果如圖7-14所示。梯形積分累加函數(shù)cumtrapz與定積分函數(shù)trapz的不同在于cumtrapz類似于不定積分，逐點給出積分的值，因而得出一個數(shù)列，而trapz只給出積分到終點的一個值。這些函數(shù)都假定步長為1，因此累加的值必須乘以dx才與積分等價。用A=cumtrapz(M)來求面積，長度M為101，只能形成100個A。而cumsum則是把101個點逐個相加，相當于多算了一個點。準確地說，可以推導出

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

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

圖7-15鉆床受力圖

解：

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

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

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

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

r=roots(a);

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

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

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

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

d=0.1219m7.3機械振動

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

解：

◆建模單自由度阻尼系統(tǒng)振動方程如下(7-9)其中m為物體質(zhì)量，k為彈簧剛度，c為阻尼常數(shù)，f為所加的外力。為了研究它的固有振動,設外加力f=0。將方程整理后寫成(7-10)式中，為固有頻率，為阻尼系數(shù)。現(xiàn)取ζ=0.1~1，ωn=10，初始條件分別為x0=1、v0=0及x0=0、v0=1，進行討論。根據(jù)微分方程理論，這樣一個常系數(shù)二階方程，其通解的形式為，將它代入方程，得知p1、p2是特征方程的兩個根，而C1、C2則由初始條件決定。很多教材都用傳統(tǒng)的解析法解這個題目，在ζ<1時，p1、p2是一對共軛復根，p1=-ζωn+jωd，p2=-ζωn-jωd，其中，此時，解就可寫成正余弦函數(shù)的形式，常數(shù)C1、C2就轉(zhuǎn)化為A和φ:其中，。用MATLAB作為計算工具對這些公式進行計算和繪圖，比用手工計算和繪圖方便得多。先設x0=1，v0=1，時間區(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; %設定不同的ζ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ù)組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) %畫出三維圖形◆程序運行結(jié)果執(zhí)行此程序即可得到圖7-16所示結(jié)果。改變初始條件為x0=0，v0=1(程序exn731b)，可得到圖7-17所示結(jié)果。實際上后一組曲線就是系統(tǒng)的脈沖過渡函數(shù)。因為脈沖函數(shù)的幅度是無窮大，而持續(xù)時間卻是無限小，其面積為一個單位，所以脈沖激勵的最后效果(在t=+eps處)可形成一個單位的初速v0，由它產(chǎn)生的波形就是脈沖過渡函數(shù)的波形。圖

7-16初值x0=1、v0=0時的振動波形

圖

7-17初值x0=0、v0=1時的振動波形

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

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

上述公式之所以繁瑣，是因為要避免復數(shù)運算。而MATLAB本身就具備復數(shù)運算的功能，可使求解變得更為簡單。設原始微分方程為

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

這樣，p1、p2和C1、C2都已求出，只要給出t數(shù)組，就求出了x。要注意的是:雖然我們知道，在系統(tǒng)的實際運動中，x必然是實數(shù)，但在復數(shù)運算中，由于計算的誤差，難免會出現(xiàn)微小的虛部，使x變成復數(shù)。這會使繪圖語句無法執(zhí)行，因此要用real語句取出實部，才能繪圖。假如其他參數(shù)不變，要求出ζ=0.3的情況下此系統(tǒ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】設單自由度阻尼系統(tǒng)的質(zhì)量m=1kg，彈簧剛度系數(shù)K=100N/m，速度阻尼系數(shù)c=4N·s/m，求它在如下外力作用下的強迫振動，并畫出t≤1.2s的波形。

解:

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

t=0∶dt:1.2;

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

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

%給出外加力的采樣值

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

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

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

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