Bài giảng Phương pháp phần tử hữu hạn

,m]=beam_elm_3(e_module,g_module,leng,h,b,rho) 
% Mo ta cac bien: 
% k – ma tran do cung phan tu (kich thuoc 6x6) 
% m – ma tran khoi luong phan tu (kich thuoc 6x6) 
% e_module – modul dan hoi 
% g_module – modul cat 
% leng – chieu dai phan tu 
% h, b – chieu cao va chieu rong mat cat ngang cua dam 
% rho – trong luong rieng vat lieu dam 
%--------------------------------------------------------------- 
% ma tran do cung 
 a1=(g_module*leng*b)/(4*h); 
 a2=(g_module*h*b)/leng; 
 a3=(e_module*h*b)/(6*leng); 
 a4=g_module*b/2; 
 k= [ a1+2*a3 -a1+a3 a4 a1-2*a3 -a1-a3 -a4;... 
 -a1+a3 a1+2*a3 -a4 -a1-a3 a1-2*a3 a4;... 
 a4 -a4 a2 a4 -a4 -a2;... 
 a1-2*a3 -a1-a3 a4 a1+2*a3 -a1+a3 -a4;... 
 -a1-a3 a1-2*a3 -a4 -a1+a3 a1+2*a3 a4;... 
 -a4 a4 -a2 -a4 a4 a2]; 
% ma tran khoi luong 
 m=zeros(6,6); 
 mass=rho*h*b*leng/4; 
 m=mass*diag([1 1 2 1 1 2]); 
%---------------------------------------------------------------------- 
SinhVienKyThuat.Com
281 
function [kk,mm]=boundary_aply_beam(kk,mm,bcdof) 
%---------------------------------------------------------------------- 
% Muc dich: 
% Ap dat cac dieu kien bien len he phuong trinh tri rieng 
% [kk]{x}=lamda[mm]{x} 
% Cu phap: 
% [kk,mm]=boundary_aply_beam(kk,mm,bcdof) 
% Mo ta cac bien: 
% kk – ma tran do cung tong the truoc khi ap dat dieu kien bien 
% mm - ma tran khoi luong tong the truoc khi ap dat dieu kien bien 
% bcdof – vecto cac bac tu do chiu rang buoc theo dieu kien bien 
%--------------------------------------------------------------------- 
 n=length(bcdof); 
 sdof=size(kk); 
 for i=1:n 
 c=bcdof(i); 
 for j=1:sdof 
 kk(c,j)=0; 
 kk(j,c)=0; 
 mm(c,j)=0; 
 mm(j,c)=0; 
 end 
 mm(c,c)=1; 
 end 
Kết quả số 
fsol = 
Mode Tần số (Hz) 
1 200 
2 1260 
3 4040 
SinhVienKyThuat.Com
282 
6.2. Chương trình tính tần số dao động tự do của khung 
Ví dụ 13.2. 
Tính tần số dao động tự do của khung công xôn như Hình 13.9. 
Tiết diện mặt cắt ngang 1010 mm; khối lượng riêng vật liệu khung là 
1000Kg/m3; môđun đàn hồi kéo nén vật liệu khung là 100gPa. Ở đây ta 
sẽ xây dựng chương trình tính với lưới gồm 10 phần tử có kích thước 
đều nhau, được mô tả như Hình 13.9. 
1 
1 
1 m 
1m 
Hình 13.9. Dao động tự do của khung phẳng 
2 
3 
4 
5 
6 
7 
x 
y 
0.01m 
0,01 
A-A 
A-A 
10 8 9 10 11 
SinhVienKyThuat.Com
283 
Chương trình nguồn 
%---------------------------------------------------------------------------- 
% Chuong trinh so 2, chuong 13 – Vi du 13.2 (P13_2) 
%---------------------------------------------------------------------------- 
% Mo ta bai toan 
% Tim tan so dao dong rieng cua khung hinh chu L duoc cau tao tu 2 
thanh 
% moi thanh co do dai 1 m. Ca 2 thanh co cung tiet dien ngang 
0.01x0.01 m. 
% Mo dul dan hoi E=100 gPa; khoi luong rieng vat lieu thanh 1000 
Kg/m^3. 
% Chuong trinh dung luoi 10 phan tu. 
% Mo ta cac bien 
% x va y = cac toa do nut toan cuc 
% k = ma tran do cung phan tu 
% kk = ma tran do cung tong the 
% m = ma tran khoi luong phan tu 
% mm = ma tran khoi luong tong the 
% index = bang ghep noi phan tu 
% bcdof = vecto chuyen vi nut chiu rang buoc theo dieu kien bien 
%---------------------------------------------------------------------------- 
clear 
b=0.01; % chieu rong mat cat thanh (mm) 
h=0.01; % chieu cao mat cat thanh (mm) 
noe=10; % so luong phan tu 
nnel=2; % so luong nut cua phan tu 
ndof=3; % so luong bac tu do cua moi nut 
nnode=(nnel-1)*nel+1; % tong so nut trong he 
sdof=nnode*ndof; % tong so bac tu do cua he 
% toa do x, y cua cac nut trong he truc chung 
x(1)=0; y(1)=0; 
x(2)=0; y(2)=0.2; 
x(3)=0; y(3)=0.4; 
x(4)=0; y(4)=0.6; 
x(5)=0; y(5)=0.8; 
x(6)=0; y(6)=1; 
x(7)=0.2; y(7)=1; 
SinhVienKyThuat.Com
284 
x(8)=0.4; y(8)=1; 
x(9)=0.6; y(9)=1; 
x(10)=0.8; y(10)=1; 
x(11)=1; y(11)=1; 
e_module=100*10^9; % modul dan hoi 
area=b*h; % dien tich mat cat ngang 
xi=(b*h^3)/12; % momen quan tinh mat cat ngang 
rho=1000; % khoi luong rieng vat lieu khung 
bcdof(1)=1; % thanh phan u tai nut 1chiu rang buoc boi dieu kien 
bien 
bcdof(2)=2; % thanh phan v tai nut 1chiu rang buoc boi dieu kien 
bien 
bcdof(3)=3; % goc xoay tai nut 1chiu rang buoc boi dieu kien bien 
kk=zeros(sdof,sdof); 
mm=zeros(sdof,sdof); 
index=zeros(nel*ndof,1); 
for iel=1:noe % xet tung phan tu 
index=feeldof1(iel,nnel,ndof); % xay dung bang ghep noi phan tu 
node1=iel; % chi so nut tong the cua nut thu 1 phan tu 'iel' 
node2=iel+1; % chi so nut tong the cua nut thu 2 cua phan tu 'iel' 
x1=x(node1); y1=y(node1); % toa do x, y cua nut thu 1 
x2=x(node2); y2=y(node2); % toa do x, y cua nut thu 2 
leng=sqrt((x2-x1)^2+(y2-y1)^2); % chieu dai phan tu 'iel' 
if (x2-x1)==0; % tinh goc giua truc dia phuong x va truc chung X 
 beta=pi/2; 
else 
 beta=atan((y2-y1)/(x2-x1)); 
end 
% tinh ma tran do cung phan tu va ma tran khoi luong phan tu 
[k,m]=frame_element_2(e_module,xi,leng,area,rho,beta,1); 
kk=kk_build_2D(kk,k,index); % ghep noi ma tran do cung tong the 
mm=kk_build_2D(mm,m,index); % ghep noi ma tran khoi luong tong 
the 
end 
[kn,mn]=boundary_aply_beam(kk,mm,bcdof); % ap dat dieu kien bien 
fsol=eig(kn,mn); % giai he phuong trinh tri rieng 
fsol=sqrt(fsol) % in ket qua 
SinhVienKyThuat.Com
285 
Các hàm sử dụng trong chương trình 
%-------------------------------------------------------------- 
function 
[k,m]=frame_element_2(e_module,xi,leng,area,rho,beta,ipt) 
%-------------------------------------------------------------- 
% Muc dich: 
% xac dinh ma tran do cung phan tu va ma tran khoi luong phan tu 
% cua phan tu khung 2 chieu 
% voi vecto chuyen vi {u_1 v_1 theta_1 u_2 v_2 theta_2} 
% Cu phap: 
% [k,m]=frame_element_2(e_module,xi,leng,area,rho,beta,ipt) 
% Mo ta cac bien: 
% k – ma tran do cung phan tu (kich thuoc 6x6) 
% m – ma tran khoi luong phan tu (kich thuoc 6x6) 
% e_module – modul dan hoi 
% xi – mo men quan tinh cua mat cat ngang 
% leng – chieu dai phan tu 
% area – dien tich mat cat ngang cua khung 
% rho – khoi luong rieng (kg/m^3) 
% beta – goc nghieng giua truc dia phuong x va truc chung X 
%-------------------------------------------------------------------------- 
% ma tran do cung trong he truc dia phuong 
 a=e_module*area/leng; 
 c=e_module*xi/(leng^3); 
 kl=[a 0 0 -a 0 0;... 
 0 12*c 6*leng*c 0 -12* c 6*leng*c;... 
 0 6*leng*c 4*leng^2*c 0 -6*leng*c 2*leng^2*c;... 
 -a 0 0 a 0 0;... 
 0 -12*c -6*leng*c 0 12*c -6*leng*c;... 
 0 6*leng*c 2*leng^2*c 0 -6*leng*c 4*leng^2*c]; 
% ma tran quay (bien doi he toa do) 
 r=[ cos(beta) sin(beta) 0 0 0 0;... 
SinhVienKyThuat.Com
286 
 -sin(beta) cos(beta) 0 0 0 0;... 
 0 0 1 0 0 0;... 
 0 0 0 cos(beta) sin(beta) 0;... 
 0 0 0 -sin(beta) cos(beta) 0;... 
 0 0 0 0 0 1]; 
% ma tran do cung phan tu tinh trong he truc chung 
 k=r'*kl*r; 
% consistent mass matrix 
 mm=rho*area*leng/420; 
 ma=rho*area*leng/6; 
 ml=[2*ma 0 0 ma 0 0;... 
 0 156*mm 22*leng*mm 0 54*mm -13*leng*mm;... 
 0 22*leng*mm 4*leng^2*mm 0 13*leng*mm … 
 -3*leng^2*mm;... 
 ma 0 0 2*ma 0 0;... 
 0 54*mm 13*leng*mm 0 156*mm -22*leng*mm;... 
 0 -13*leng*mm -3*leng^2*mm 0… 
 -22*leng*mm 4*leng^2*mm]; 
% ma tran khoi luong trong he toa do chung 
 m=r'*ml*r; 
Kết quả số 
Mode Tần số (Hz) 
1 34 
2 92 
3 455 
4 667 
SinhVienKyThuat.Com
287 
7. BÀI TẬP 
13.1. Cho kết cấu dầm như Hình 13.7.1. 
a. Hãy xây dựng ma trận độ cứng tổng thể cho kết cấu và ma trận 
khối lượng hệ; 
b. Thực hiện tính toán bằng tay, xác định tần số dao động tự do nhỏ 
nhất của dầm; 
c. Phát triển chương trình P13_1 để thực hiện theo các yêu cầu ở ý 
a và b trên đây; và kiểm tra, so sánh kết quả. 
13.2. Phát triển chương trình P13.1, xác định các tần số dao 
động tự do của kết cấu dầm như Hình 13.7.2. So sánh kết quả khi tính ở 
hai trường hợp: sử dụng lưới 2 phần tử và lưới 4 phần tử. 
13.3. Bằng cách tính tay và phát triển chương trình P13_1 để 
xác định hai tần số dao động tự do thấp nhất của hệ thống trục thép 
mang các bánh răng (coi như khối lượng tập trung) như chỉ ra trên Hình 
13.7.3. Ở hai trường hợp như sau: 
a. Coi cả 3 ổ bi như các gối đơn 
b. Mỗi ổ bi được coi là các gối đỡ mềm, độ cứng là 45kN/mm. 
800 mm 
x 
Hình 13.7.2 
75 mm 
25 mm 
300 mm 400 mm 
A1=1200 mm2 A2=900 mm2 
x 
Hình 13.7.1 
SinhVienKyThuat.Com
288 
 13.4. Phát triển chương trình P13_2 để xác định hai tần số dao 
động tự do thấp nhất của khung thép như chỉ ra trên Hình 13.7.4. 
600 
600 
300 
15 
30 
0,5 
1 
Khung thép Mặt cắt ngang 
Hình 13.4 
10 kN 
5 kN 
75mm 75mm 45mm 45mm 
Hình 11.7.3 
SinhVienKyThuat.Com
289 
TÀI LIỆU THAM KHẢO 
1. Trần Ích Thịnh - Trần Đức Trung - Nguyễn Việt Hùng. 
Phương pháp phần tử hữu hạn trong kỹ thuật. Đại học Bách Khoa 
– Hà Nội, 2000. 
2. Tirupathi R. Chandrupatla – Ashok D. Belegundu. Introduction 
Finite Elements in Engineering. Third Edition. 
3. Young W. Hwon - Hyochoong Bang. The Finite Element Method 
Using MATLAB. Second Editor. CRC Press, 2000. 
4. J. N. Reddy. An Introduction To The Finite Element Method. 
Third Edition. Tata McGraw-Hill, 2005. 
5. Klaus – Jürgen Bathe. Finite Element Procedures. Prentice-Hall 
of India, New Delhi, 2005. 
6. K Chandrashekhara. Theory of Plates. Universities Press, 2001. 
7. O. C. Zienkiewicz and R. L. Taylor. The Finite Element Method, 
Fifth Edition. Volume 2, Solid Mechanics. Butterworth Heinemann, 
2000. 
8. O. C. Zienkiewicz and K. Morgan. Finite Element and 
Approximation. New York: Wiley – Iterscience, 1982. 
9. Akin J. E. Finite Element for Analysis and Design. Academic 
Press Limited, London, 1994. 
10. Batoz J. L. Et Dhatt DG. Modélesation des structues par 
élements finis.Vol. 1, 2, 3. Ed. Hermès. Paris, 1995. 
11. Dhatt G. Et Touzot G. Une présentation de la méthode des 
élements finis. Maloine S.A. Editeur, 1981. 
12. Ochoa O. O, Readdy, J. N. Finite Element Analysis of Composite 
Laminates. Klwer Academic Publisher, 1992. 
SinhVienKyThuat.Com