Đề tài Matlab - Tìm hiểu và ứng dụng giải một số bài toán kĩ thuật

c tiết diện nguy hiểm là A và B dựa vào các ph−ơng 
trình cân bằng sau: 
 MB + lBC x RC = 0 
(5) 
 MA+(-lA x R0) = 0 
(6) 
Để giải các ph−ơng trình véc tơ (3), (4), (5) và (6) ta đ−a ra các véc tơ 
đơn vị của hệ toa độ Đề các là i, j và k và dùng hàm cross(...,...) để nhân 
có h−ớng hai véc tơ. Cách làm đ−ợc thể hiện rõ trong file plkd.m d−ới 
đây: 
% File plkd.m 
% Tinh phan luc tai khop dong va mo men uon 
% Cac Vec to don vi cua he toa do De cac 
I=[1 0 0]'; 
J=[0 1 0]'; 
K=[0 0 1]'; 
% Vec to vi tri 
la=100; lb=250; lc=350; lbc=100; 
La=100*I; Lb=250*I; Lc=350*I; Lbc=lbc*I; 
% Luc tac dung 
Pay=-1.4*J; Paz=.5*K; 
Pby=-.9*J; Pbz=1.4*K; 
% Tinh phan luc goi C 
LcRc=-cross(La,Pay)-cross(La,Paz)-cross(Lb,Pby)-
cross(Lb,Pbz); 
Rc0=LcRc/lc; % Vec to phan luc tai 
goi C(chua quay) 
Q=[0 0 0; 0 0 -1; 0 1 0]; % Ma tran chuyen 
(quay) 
Rc=-Q*Rc0 % Vec to phan luc tai 
goi C 
% Tinh phan luc tai goi O 
 R0=-(Pay+Paz+Pby+Pbz+Rc) 
% Tinh mo men uon tai tiet dien B va A 
 Mb=-cross(Lbc,Rc) 
 Ma=-cross(-La,R0) 
% --------------------------------------------- 
Chạy ch−ơng trình plkd.m trên màn hình Command Window, kết quả 
của các phản lực khớp động tại O , C và các mô men uốn tại mặt cắt A và 
B nhận đ−ợc nh− sau: 
Rc = 
 0 
 1.0429 
 -1.1429 
R0 = 
 0 
 1.2571 
 -0.7571 
Mb = 
 0 
 -114.2857 
 -104.2857 
Ma = 
 0 
 -75.7143 
 -125.7143 
ằ plkd 
Rc = 
 0 
 1.0429 
 -1.1429 
R0 = 
 0 
 1.2571 
 -0.7571 
Mb = 
 0 
 -114.2857 
 -104.2857 
Ma = 
 0 
 75.7143 
 125.7143 
6. Tính dao động 
a. Tính dao động của hệ một bậc tự do 
Giả sử có hệ dao động một bậc tự do nh− hình vẽ: m kà khối l−ợng của 
vật, c là độ cứng lò xo (N/m), k là hệ số cản của giảm chấn thuỷ lực 
(Ns/m), còn F là lực kich thích với F=F0coswt ( F0 là biên độ của lực – 
N, w là tần số góc – rad/s). 
Viết ph−ơng trình định luật 2 cho vật ; 
 mx’’= -kx’-cx +F0coswt (1) 
trong đó x, x’, x’’ lần l−ợt là chuyển vị, x 
vận tốc và gia tốc của vật. k F 
ta biến đổi (1) nh− sau: 
 c m 
 x’’= -kx’/m - cx/m + +F0coswt/m (2) 
Đặt x1=x’, x2= x thì x1’= x’’, x2’= x1, ph−ơng trình vi phân cấp 2 trên sẽ 
đ−ợc viết thành hệ hai ph−ơng trình bậc nhất với hai ẩn là x1 và x2 nh− 
d−ới đây: 
 x1’= - kx1/m - cx2m + F0coswt/m (3)
 x2’= x1 
Ch−ơng trình để giải hệ (3) đ−ợc viết trong các file ml00.m và ml0.m 
nh− d−ới đây: 
% File ml00.m 
% Chuongtrinh tinh dao dong he 1 bac tu do co 
giam chan 
% thuy luc chui luc kich dong dang F=F0cos(wt) 
global m k c w F0 
m=100; % Khoi luong vat - kg 
c=10000; % Do cung lo xo - N/m 
k=80; % He so can nhot Ns/m 
w1=sqrt(c/m) % Tan so dao dong rieng 
w=6; % Tan so luc kich dong 
tc=15; % Thoi gian dao dong 
x0=[0 .01 ]' ; % Dieu kien ban dau (van 
toc, toa do) 
F0=20 ; % Bien do luc kich dong - N 
 [t,x]=ode45('ml0',tc,x0); 
subplot(211) 
plot(t,x(:,2)) 
title(' TINH DAO DONG HE MOT BAC TU DO') 
ylabel('Chuyen vi cua vat - cm') 
grid 
subplot(212) 
plot(t,x(:,1)) 
xlabel('thoi gian -s') 
ylabel('Van toc cua vat - m/s') 
grid 
pause 
% File ml0.m 
function xc=ml0(t,x) 
 global m k c w F0 
 xc=zeros(2,1); 
 xc(1)=-k/m*x(1)-c/m*x(2)+F01*cos(w*t)/m; 
 xc(2)=x(1); 
Kết quả chạy ch−ơng trình cho d−ới dạng đồ thị d−ới đây: 
0 5 1 0 1 5
-1
-0 .5
0
0 .5
1
 TIN H D A O D O N G H E M O T B A C TU D O
C
hu
ye
n 
vi
 c
ua
 v
at
- 
cm
0 5 1 0 1 5
-0 .1
-0 .0 5
0
0 .0 5
0 .1
tho i g ia n -s
V
an
 t
oc
 c
ua
 v
at
 -
 m
/s
b. Tính dao động hệ hai bậc tự do 
Có hệ dao dộng hai bậc tự do nh− hình vẽ: khối l−ợng 
của các vật là m1, m2, độ cứng các lò xo là c1và c2, hệ x2 
F1
số cản nhớt là k1 và k2, lực kích thich đặt lên hệ là c2 k2
F1=F10sinw1t, F2=F20sinw2t. 
 x1 
F2 
 c1 
k1
Viết hệ ph−ơng trình định luật 2 Niu tơn cho 
từng vật: 
m1x1’’=- k1x1’- c1x1+k2(x2’-x1’)+c2(x2-x1) +F01sinw1t 
m2x2’’= - k2 (x2’-x1’)- c2(x2 -x1)+F02sinw2t 
x1’’ = - (k1+k2)x1’/m1- (c1+c2)x1/ m1+k2x2’/ m1+c2x2 / m1+ F01sinw1t/ m1
 (1) 
x2’’ = k2x1’/ m2 + c2x1 / m2 - k2x2’ / m2 - c2x2/ m2 + F02sinw2t/ m2
T−ơng tự nh− khi giải ph−ơng trình vi phân câp hai ở phần trên, ta cũng 
đặt: 
y1 = x1’, y2 = x1, y3 = x2’, y4 = x2. 
 (2) 
Rút ra: 
 y1’=x1’’ 
 (3) 
 y2’=y1
 y3’= x2’’ 
 y4’=y3 
Thay (1) vào (3) và chú ý đến (2) ta lập đ−ợc hệ gồm bốn ph−ơng trình 
vi phân tuyến tính bậc nhất. Hệ ph−ơng trình này và cách giải nó đ−ợc 
cho trong file ml11.m và ml1.m ở d−ới đây. 
Nhân tiện đây ta cũng thử tác dụng của giảm chấn động lực: theo lí 
thuyết, nếu lực kích thích vào vật 2 bằng 0 ( F2=0), và các thông số của 
bộ giảm chấn động lực gồm khối l−ợng m2 và lò xo có độ cứng c2 thoã 
mãn điều kiện: c2/m2 = w1
2 ( w1 là tần số lực kích thích của lực F1 lên vật 
1) thì dao động của vật chính 1 có khối l−ợng m1 đó sẽ bằng 0. Trong 
ch−ơng trình ta cho F20=0, m2=10 kg, c2=510 N/m và w1= 22 /mc . 
% File ml11.m 
% Tinh dao dong he hai bac tu do 
global m1 m2 k1 k2 c1 c2 w F10 F20 
m1=100; % Khoi luong vat thu nhat - kg 
c1=10000; % Do cung lo xo thu nhat - N/m 
m2=10; % Khoi luong vat thu 2 - kg 
c2=510; % Do cung lo xo thu 2 - N/m 
k1=220; k2=0; % He so can nhot cua giam chan 
- Ns/m 
w20=sqrt(c2/m2) 
x0=[0 .0 0 .0]'; % Dieu kien ban dau 
F10=100; % Bien do luc tac dung vao 
vat 1 – N 
F20=0; % -- 2 
w1=w20; % Tan so goc cua luc k.thich 
- rad/s 
 [t,y]=ode45('ml1',tc,x0); 
subplot(211) 
plot(t,y(:,2)*100) 
title('TINH DAO DONG HE HAI BAC TU DO') 
ylabel('Chuyen vi cua vat 1 - cm') 
grid 
subplot(212) 
plot(t,y(:,4)*100) 
xlabel(' Thoi gian - s') 
ylabel('Chuyen vi cua vat 2 - cm') 
grid 
pause 
% File ml1.m 
% Chuong trinh con tinh dao dong he hai bac tu 
do 
function yc=ml1(t,y) 
 global m1 m2 k1 k2 c1 c2 w F10 F20 
 yc=zeros(4,1); 
 yc(1)=-(k1+k2)/m1*y(1)-(c1+c2)/m1*y(2)+... 
k2/m1*y(3)+c2/m1*y(4)+F10*sin(w1*t)/m1; 
 yc(2)= y(1); 
 yc(3)= k2/m2*y(1)+c2/m2*y(2)-k2/m2*y(3)-
c2/m2*y(4)+... 
 F20*sin(w2*t)/m2; 
 yc(4)= y(3); 
Chạy ch−ơng trình sẽ cho kết quả d−ới dạng đồ thị sau: 
0 2 4 6 8 1 0 1 2 1 4 1 6 1 8 2 0
-2
-1
0
1
2
TIN H D A O D O N G H E H A I B A C TU D O
C
hu
ye
n 
vi
 c
ua
 v
at
 1
 -
 c
m
0 2 4 6 8 1 0 1 2 1 4 1 6 1 8 2 0
-4 0
-2 0
0
2 0
4 0
 Th o i g ia n - s
C
hu
ye
n 
vi
 c
ua
 v
at
 2
 -
 c
m
Ta thấy một điều thú vị là : đúng nh− dự đoán, biên độ dao động của vật 
chính 1 giảm dần theo thời gian và sẽ tắt hẳn sau ít giây nữa. 
c. Xác định tần số dao động riêng của hệ nhiều bậc tự do 
Với hệ cơ học nhiều bậc tự do việc xác định các tần số dao động riêng 
của hệ có một ý nghĩa quan trọng. Matlab giúp ta giải bài toán này không 
mấy khó khăn. 
Ví dụ: Có hệ gồm bốn toa xe móc nối với nhau, hãy xác định tần số dao 
động riêng của hệ theo chiều dọc. Biết khối l−ợng các toa xe là m1=45T, 
m2=35T, m3=40T, m4=50T, độ cứng các lò xo là c1=c2=800KNm, 
c3=1000KN/m. 
 m1 x1 m2 x2 m3 x3 m4 x4
Để giải bài toán tr−ớc hết ta phải lập các ma trận khối l−ợng và ma trân 
độ cứng nh− sau: 
- Ma trận khối l−ợng: 
 M=diag(m1, m2, m3, m4); ( Ma trận chéo có các phần tử là các 
khối l−ợng) 
-Ma trận độ cứng: 
 c1 -c1 0 0 
 C = -c1 c1+c2 -c2 0 
 0 -c2 c2+c3 -c3
 0 0 -c3 c3
- Dùng hàm eig để xác định các trị riêng và véc tơ riêng: 
 [ X, L] = eig( C, M); 
trong đó L là ma trận đ−ờng chéo mà mỗi phần tử trên đ−ờng chéo là một 
trị riêng (eigenvalue) là w1
2, w2
2, w3
2
 và w4
2 . Còn w1 w2, w3 và w4 chính là 
các tần số dao động riêng của hệ (rad/s). 
Ma trận (hàng) của các tần số dao động riêng đ−ợc xác định nh− sau: 
 w = diag(sqrt(L)); 
Kết quả tính bằng số cụ thể nhận đ−ợc nhờ chạy ch−ơng trình trong file 
tansd.m d−ới đây: 
% File tánd.m 
% CHUONG TRINH TINH TAN SO DAO DONG RIENG 
% CUA HE NHIEU BAC TU DO 
m1=45000; m2=35000; 
 m3=40000; m4=50000; % Khoi luong cac toa xe 
c1=800000;c2=c1; % Do cung lo xo N/m. 
c3=1000000; 
% Ma tran khoi luong va ma tran do cung 
m=[m1 m2 m3 m4]; 
M=diag(m); 
C=[c1 -c1 0 0; -c1 c1+c2 -c2 0 ; 0 -c2 c2+c3 
-c3;... 
 0 0 -c3 c3]; 
% Tinh tri rieng va vec to rieng 
 ,L]=eig(C,M); 
wb=diag(L); % Cac tri rieng 
w=sqrt(wb) % Cac tan so dao dong rieng - 
rad/s 
Các giá trị tần số dao động riêng đ−ợc ch−ơng trình cho d−ới đây: 
 w = 
 8.6524 
 6.5472 
 0 + 0.0000i 
 3.2807 
Ta thấy có một tần số dao động riêng bằng 0 : w3=0. Điều này về mặt lí 
thuyết dao động thì hoàn toàn hợp lí vì nó t−ơng ứng vối một dạng dao 
động riêng của hệ trong đó toàn hệ chuyển động nh− một vật rắn. 
Kết luận 
Qua hai phần nghiên cứu tìm hiểu phần mềm lập trình Matlab và ứng 
dụng nó để giải một số bài toán kĩ thuật, chúng tôi thấy rằng đây là một 
phần mềm rất mạnh, dễ sử dụng và rất có hiệu quả trong tính toán, đặc 
biệt thích hợp cho việc nghiên cứu trong lĩnh vực khoa học kĩ thuật. 
 Đồng thời chúng tôi cũng thấy rõ rằng: để hiểu cặn kẽ cũng nh− sử dụng 
đ−ợc hết khả năng của Matlab thì đó cũng không phải là chuyện dễ, vì 
bản thân nó là một phần mềm lớn. Tuy nhiên, trong phạm vi của đề tài 
nghiên cứu này, những vấn đề cơ bản và những ứng dụng phổ biến của 
phần mềm đã đ−ợc trình bày và đã có những ví dụ minh hoạ cụ thể. Từ 
những phần cơ bản này, ng−ời đọc qua sử dụng thực tế phần mềm sẽ có 
khả năng tự tìm hiểu sâu hơn nhờ tra cứu phần Help trong máy. 
Cách trình bày chúng tôi đã cố gắng sao cho ngắn gọn, dễ hiểu và thiên 
về thực hành. Một trong những mục đích chính của nhóm nghiên cứu là 
biến kết quả nghiên cứu này thành tài liệu để tìm hiểu và sử dụng Matlab 
cho những ng−ời quan tâm, đặc biệt là cho sinh viên. Hi vọng là với tài 
liệu này, ng−ời học có thể tự thực hành trên máy tính để lập các ch−ơng 
trình tính toán cho bài toán của mình. 
Do khuôn khổ thời gian có hạn, còn một phần thứ hai nữa của Matlab mà 
chúng tôi ch−a kịp nghiên cứu: đó là phần SIMULINK. Đây là phần 
mềm mô phỏng rất mạnh và hiệu quả của Matlab, cho phép giảiđ−ợc 
nhiều loại bài toán phức tạp. Chúng tôi hi vọng sẽ có điều kiện nghiên 
cứu tiếp trong thời gian tới.