Giáo trình kỹ thuật - Phương pháp phần tử hữu hạn

h h xz yz x y AA AA dz Q Q Q k k     (12.24) với:   4,5 j i,);(' 1 1     kk k n k ijij hhCA (12.25) Trong tính toán, để có các thành phần ứng suất cắt đạt độ chính xác cao so với lời giải của lý thuyết đàn hồi, người ta phải sử dụng thêm các hệ số hiệu chỉnh ứng suất cắt. Mindlin đã đề xuất hàm hiệu chỉnh dưới dạng:               2 2/ 1 4 5)( h zzf (12.26) Khi ấy, từ (12.24) các thành phần lực cắt sẽ được tính theo biểu thức: 0 55 0 45 0 45 0 44 ; xzyzxxzyzy AAQAAQ   (12.27) trong đó:          n k kkkkkijij h hhhhCA 1 2 3 1 3 1 4,5 j i,; 1)( 3 4' 4 5 (12.28) Chú ý: Hàm hiệu chỉnh (hay hệ số hiệu chỉnh) được chọn tuỳ theo vật liệu và phụ thuộc vào tác giả. Chẳng hạn, Reissner chọn hệ số bằng 5/6 cho tấm trực hướng và đồng nhất. SinhVienKyThuat.Com 245 e. Phương trình ứng xử cơ học của tấm nhiều lớp. Phương trình quan hệ của tấm nhiều lớp được viết dưới dạng sau:                                                                                0 0 0 0 0 5545 4544 662616662616 262212262212 161211161211 662616662616 262212262212 161211161211 000000 000000 00 00 00 00 00 00 xz yz xy y x xy y x x y xy y x xy y x AA AA DDDBBB DDDBBB DDDBBB BBBAAA BBBAAA BBBAAA Q Q M M M N N N         (12.29) Hay có thể viết dưới dạng thu gọn:                                                             0'00 0 0   m A DB BA Q M N (12.30) trong đó: A =[Aij] (i, j =1, 2, 6) là ma trận độ cứng màng; D là ma trận độ cứng uốn; B là ma trận tương tác màng-uốn-xoắn; A’ =[Aij] (i, j =4, 5) là ma trận độ cứng cắt. Các phần tử của chúng được xác định theo các biểu thức:   )( 1 1 '    k n k kkijij hhQA (12.31a)   )( 2 1 2 1 2 1 '     kk n k kijij hhQB (12.31b)   )( 3 1 3 1 3 1 '     kk n k kijij hhQD (12.31c)    Txyyxm 000   là ma trận biến dạng màng,    Txyyx   là ma trận độ cong của tấm chịu uốn,    Txzyz 000   là ma trận biến dạng cắt của mặt trung bình. SinhVienKyThuat.Com 246 4.2. Mô hình hóa PTHH bài toán tấm composite lớp chịu uốn Dưới đây, chúng ta sẽ sử dụng phần tử tứ giác đẳng tham số, bốn nút ở đỉnh để giải bài toán uốn tấm composite. a. Véctơ chuyển vị nút phần tử Theo lý thuyết chuyển vị bậc nhất của Mindlin trên đây, tại mỗi nút của phần tử sẽ có 5 thành phần chuyển vị:  Tixiyiiii wvud 000 (12.32) Các thành phần chuyển vị này tương ứng với 5 bậc tự do tại mỗi nút:  Tiiiiii qqqqqq 4321  (12.33) b. Véctơ biến dạng Véctơ chuyển vị tại một điểm bất kỳ của phần tử tấm là:  Txywvud 000 (12.34) Véctơ này được xác định nhờ phép nội suy từ các chuyển vị nút và các hàm dạng: i i idNd    4 1 (12.35) Trong đó (nhắc lại về phần tử tứ giác bậc nhất đẳng tham số): các hàm dạng và đạo hàm của chúng được biểu diễn như sau:               11111111 4 1N (12.36a)                            1111 4 1 1111 4 1 ' ' N N (12.36b) Toạ độ hình học của một điểm bất kỳ trong phần tử cũng được nội suy từ toạ độ các điểm nút như sau: x = N1 x1 + N2 x2+ N3 x3+ N4 x4 y = N1 y1 + N2 y2 + N3 y3+ N4 x5 (12.37) Quan hệ giữa các đạo hàm riêng trong các hệ toạ độ thực và hệ toạ độ quy chiếu của một hàm f bất kỳ được biểu diễn như sau: SinhVienKyThuat.Com 247              y y fx x ff v y y fx x ff              (12.38) Hay                                  y f x f Jf f   (12.39) Với J là ma trận Jacobian được xác định bởi:                                    4 1 22 4 1 21 4 1 12 4 1 11 ; ; i i i i i i i i i i i i yNxJxNxJ yNyJxNxJ   (12.40) Thay (12.36b) vào (12.40) ta được:         432111 14 11 4 11 4 11 4 1 xxxxJ   (12.41a)         432112 14 11 4 11 4 11 4 1 yyyyJ   (12.41b)         432121 14 11 4 11 4 11 4 1 xxxxJ   (12.41c)         432122 14 11 4 11 4 11 4 1 yyyyJ   (12.41d) Và quan hệ ngược:                                                                f f JJ JJ Jf f J y f x f 1121 12221 det 1 (12.42) SinhVienKyThuat.Com 248 Khi đó, trường biến dạng được xác định qua véctơ chuyển vị dưới dạng: d L L Lm                      3 2 1 0   (12.43) Trong đó: Li là các ma trận toán tử đạo hàm được biểu diễn như sau:                            000 0000 0000 1 xy y x L ;                            xy y x L 000 0000 0000 2 ;                  0100 1000 3 x yL (12.44) Gọi a là véctơ chuyển vị nút của phần tử:    TTTTT dddda 4321 (12.45) Khi đó, các thành phần biến dạng được biểu diễn qua véctơ chuyển vị nút phần tử như sau:      4 1 111 i iim aBdNLdL (12.46a)      4 1 222 i ii aBdNLdL (12.46b)      4 1 333 0 i ii aBdNLqL (12.46c) Trong đó Bi được gọi là các ma trận biến dạng-chuyển vị và được xác định như sau:      413121111 NLNLNLNLB (12.47a)      423222122 NLNLNLNLB (12.47b)      433323133 NLNLNLNLB (12.47c) SinhVienKyThuat.Com 249 c. Ma trận độ cứng của phần tử của tấm Biểu thức của năng lượng biến dạng đàn hồi là: dVU V T e  2 1 (12.49) Khai triển (12.49), với chú ý 0z ta được:  dzdSU S h h xzxzyzyzxyxyyyxxe     2 2 2 1  (12.50) Thay các quan hệ nội lực ứng suất (12.14) và (12.19) vào (12.50), rồi lấy tích phân dọc theo z ta sẽ được biểu thức năng lượng biến dạng đàn hồi như sau: dS QQMkN MkNMkN U S xzxzyzyzxyxyxyxy yyyyxxxx e             000 00 2 1   (12.51) Thay biểu thức (12.30) vào (12.51) ta được:                          dS AkDkBk kBA U Se TT m T T mm T m e             0'02 1   (12.52) Thay các biểu thức (12.44a-c) vào (12.52) ta được: dS aBABaaDBBa aBBBaaBBBaaABBa U Se TTTT TTTTTT e             2322 122111 '2 1 (12.53) Cuối cùng, ta có thể viết (12.53) dưới dạng cô đọng sau:  akaU eTe 2 1  (12.54) Trong đó ek được gọi là ma trận độ cứng phần tử như sau:  dSBABDBBBBBBBBABBk Se TTTTTe   3322122111 ' (12.55) SinhVienKyThuat.Com 250 d. Véctơ lực nút phần tử Công do tải trọng phân bố p(x,y) tác dụng vuông góc với bề mặt tấm được biểu diễn bởi biểu thức:      ee S T P T S dSyxpBadSyxwyxp ),(),(),( (12.56) Với phần tử tứ giác bốn nút, ma trận Bp được xác định bởi:      4321 ppppp BBBBB (12.57) với:         0011 4 1001 pL (12.58a)         0011 4 1002 pL (12.58b)         0011 4 1003 pL (12.58c)         0011 4 1004 pL (12.58d) 5. CHƯƠNG TRÌNH TÍNH TẤM COMPOSITE LỚP CHỊU UỐN Xét tấm composite lớp chữ nhật với liên kết bản lề trên 4 cạnh chịu uốn bởi tải trọng phân bố đều như Hình 12.5. Tìm độ võng của tấm dựa trên lý thuyết tấm Mindlin. Kích thước của tấm chữ nhật: a=254 mm, b=508mm, h=12,7mm; cấu hình đối xứng, đúng trục:(900/00/00/900); tải trọng phân bố đều trên toàn bộ bề mặt tấm: p=0,6895 N/mm2. Cơ tính của các lớp như nhau: E1=144,8 gPa, E2=E3=9,65 gPa, G12=G13=4,14 gPa, G23=3,45 gPa; các hệ số Poatxong : 3,0231312   . Ở đây, với ma trận độ cứng uốn phần tử, ta lấy tích phân số 22 điểm Gauss; với ma trận độ cứng cắt, ta lấy tích phân số 1 điểm Gauss. Điều kiện biên sẽ được mô tả như sau: - Tại các nút nằm trên các cạnh biên song song với trục x: u, w và x bằng không. SinhVienKyThuat.Com 251 - Tại các nút nằm trên các cạnh biên song song với trục y: v, w và y bằng không. Chương trình nguồn %---------------------------------------------------------------------------- % Chuong trinh chuong 12 – Vi du 12.2 (P12_2) z y x 1 2 3 4 5 6 7 8 9 1 Hình 12.5. (a). Sơ đồ hoá tấm composite chữ nhật chịu uốn (b). Lưới 44 phần tử trên tấm chữ nhật (b) (a) p 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 y a b h x SinhVienKyThuat.Com 252 %---------------------------------------------------------------------------- % Mo ta bai toan % Tam composite, chiu lien ket tua ban le tren 4 canh; chiu tai trong deu. % Tim do vong cua tam, su dung phan tu tu giac bac nhat dang tham so % dua tren ly thuyet tam bac nhat Mindlin. (so do luoi nhu Hinh 12.5) % Mo ta cac bien % k = ma tran do cung phan tu % kb = ma tran do cung phan tu ung voi cac thanh phan uon % ks = ma tran do cung phan tu ung voi cac thanh phan cat % f = vecto luc nut phan tu % kk = ma trando cung tong the % ff = vecto luc nut phan tu % disp = vecto chuyen vi nut tong the % gcoord = toa do nut tong the % nodes = bang dinh vi nut phan tu % index = bang ghep noi phan tu % pointb = toa do cac diem Gauss cho tich phan so cac thanh phan uon % weightb = he so trong so cac diem Gauss cho cac thanh phan uon % points = toa do cac diem Gauss cho tich phan so cac thanh phan cat % weights = he so trong so cac diem Gauss cho cac thanh phan cat % bcdof = cac chuyen vi nut chiu rang buoc boi dieu kien bien % bcval = gia tri cac chuyen vi nut chiu rang buoc % B_b = ma tran bien dang - chuyen vi doi voi cac bien dang uon % D_b = ma tran do cung chong uon cua vat lieu % B_s = ma tran bien dang - chuyen vi doi voi cac bien dang cat % D_s = ma tran do cung chong cat cua vat lieu %---------------------------------------------------------------------------- %----------------------------------------- % cac tham so dieu khien luoi %----------------------------------------- SinhVienKyThuat.Com 253 clear noe_x=4; % so luong phan tu theo phuong x noe_y=4; % so luong phan tu theo phuong y noe=noe_x*noe_y; % tong so phan tu cua he nnel=4; % so luong nut cua phan tu ndof=5; % so bac tu do cua nut nnode=(noe_x+1)*(noe_y+1); % tong so nut cua he sdof=nnode*ndof; % tong so bac tu do cua ca he edof=nnel*ndof; % so bac tu do cua moi phan tu lop=254; % kich thuoc (chieu dai) tam (mm) ratio_b_a=2; wop=lop*ratio_b_a; % kich thuoc (chieu rong) tam (mm) ratio_h_a=0.05; nol_p=4; % so lop vat lieu top = lop*ratio_h_a; % chieu day tam aol_p=[0 90*pi/180 90*pi/180 0]; % cac goc dat cot (radial) nog_xb=2; nog_yb=2; % 2x2 diem Gauss-Legendre cho TP uon nog_b=nog_xb*nog_yb; % tong so diem Gauss nog_xs=1; nog_ys=1; % 1x1 Gauss-Legendre cho TP cat nog_s=nog_xs*nog_ys; % tong so diem Gauss %------------------------------------------ % du lieu tinh chat cau vat lieu lop %------------------------------------------ type_material=1; % chon 1 trong cac loai vat lieu switch (type_material) case 1 emodule_1=175.0e3; % modul dan hoi E_1 (N/mm^2) emodule_2=emodule_1/25; % modul dan hoi E_2 emodule_3=emodule_2; % modul dan hoi E_3 gmodule_12=0.5*emodule_2; % modul dan hoi truot G_12 gmodule_13=gmodule_12; % modul dan hoi truot G_13 gmodule_23=0.2*emodule_2; % modul dan hoi truot G_23 SinhVienKyThuat.Com 254 % He so Poisson nuy_12 = 0.25; nuy_13 = 0.25; nuy_23 = 0.45; case 2 emodule_1=144.8e3; emodule_2=9.65e3; emodule_3=emodule_2; gmodule_12=4.14e3; gmodule_13=gmodule_12; gmodule_23=3.45e3; nuy_12 = 0.3; nuy_13 = 0.3; nuy_23 = 0.49; case 3 emodule_1=19.2*psi_Pa; emodule_2=1.56*psi_Pa; emodule_3=emodule_2; gmodule_12=0.82*psi_Pa; gmodule_13=gmodule_12; gmodule_23=0.49*psi_Pa; nuy_12 = 0.24; nuy_13 = 0.24; nuy_23 = 0.49; end %---------------------------- % tai trong gay uon phan bo deu %---------------------------- p=13.8e-3; % tai trong phan bo deu (N/mm^2) %--------------------------------------------- % du lieu toa do nut gcoord(i,j) % trong do, i la chi so nut va j=1, chi toa do x va j=2 chi toa do y %--------------------------------------------- len_x_elm = lop/noe_x; len_y_elm = wop/noe_y; SinhVienKyThuat.Com 255 for row_index=1:noe_y+1 for col_index=1:noe_x+1 gcoord(((row_index-1)*(noe_x+1)+col_index),1) = … (col_index-1)*len_x_elm; gcoord(((row_index-1)*(noe_x+1)+col_index),2) = … (row_index-1)*len_y_elm; end end %--------------------------------------------------------- % Tinh do cao tam lop cua cac lop vat lieu %--------------------------------------------------------- z_p(1)= -top/2.0; for k=1:nol_p tol_p(k)=top/nol_p; % do day cua cac lop deu nhau (mm) z_p(k+1)=z_p(k)+tol_p(k); end %--------------------------------------------------------- % du lieu cua mang chi so nut tong the cua moi phan tu % nodes(i,j), trong do i la chi so phan tu, j la chi so nut %--------------------------------------------------------- % nodes=[1 2 5 4; 2 3 6 5; 4 5 8 7; 5 6 9 8]; % 4 phan tu 4 nut. for row=1:noe_y for col=1:noe_x elm=(row-1)*noe_x+col; nodes(elm,1)= (row-1)+elm; nodes(elm,2)= nodes(elm,1)+1; nodes(elm,4)= nodes(elm,1)+(noe_x+1); nodes(elm,3)= nodes(elm,4)+1; end end %-------------------------------------------- % dieu kien bien %-------------------------------------------- SinhVienKyThuat.Com 256 i=1; for node_indx=1:nnode temp_n=(node_indx-1)*5; % cac bien x=0 va x=a if ((gcoord(node_indx,1)= =0)||(gcoord(node_indx,1)= =lop)) bcdof(i)=temp_n+2; % v_i i=i+1; bcdof(i)=temp_n+3; % w_i i=i+1; bcdof(i)=temp_n+5; % teta_x_i i=i+1; end % cac bien y=0 va y=b if ((gcoord(node_indx,2)==0)||(gcoord(node_indx,2)==wop)) bcdof(i)=temp_n+1; % u_i =0 i=i+1; % tranh cap nhat cac nut goc if ((gcoord(node_indx,1)~= 0)&&(gcoord(node_indx,1)~=lop)) bcdof(i)=temp_n+3; i=i+1; end bcdof(i)=temp_n+4; i=i+1; end end bcval=zeros(size(bcdof)); % gia tri cac chuyen vi nut bi rang buoc = 0 %---------------------------------------------- % khoi tao cac vec to va ma tran %---------------------------------------------- ff=zeros(sdof,1); kk=zeros(sdof,sdof); disp=zeros(sdof,1); % vecto chuyen vi nut tong the SinhVienKyThuat.Com 257 index=zeros(edof,1); A=zeros(3,3); % ma tran do cung mang B=zeros(3,3); % ma tran tuong tac mang-uon-xoan D=zeros(3,3); % ma tran do cung uon A_t=zeros(2,2); % ma tran do cung cat %------------------------------------------------------------------------------- % tinh ma tran do cung phan tu, vecto luc nut va ghep noi phan tu %------------------------------------------------------------------------------- % % voi cac thanh phan gay uon % xac dinh toa do va trong so cac diem gauss [pointb,weightb]=Gauss_Point_2D(nog_xb,nog_yb); % voi cac thanh phan gay cat % xac dinh toa do va trong so cac diem gauss [points,weights]=Gauss_Point_2D(nog_xs,nog_ys); % tinh cac ma tran do cung vat lieu [A, B, D, A_t]=ABD_matrix(emodule_1,emodule_2,emodule_3,... gmodule_12,gmodule_13,gmodule_23,… nuy_12,nuy_13,nuy_23,nol_p,aol_p, z_p); for iel=1:noe % xet tung phan tu for i=1:nnel nd(i)=nodes(iel,i); % chi so tong nut the cua phan tu (iel) xcoord(i)=gcoord(nd(i),1); % toa do nut (phuong x) ycoord(i)=gcoord(nd(i),2); % toa do nut (phuong y) end k=zeros(edof,edof); f=zeros(edof,1); kb=zeros(edof,edof); ks=zeros(edof,edof); %------------------------------------------------------ % tinh tich phan so cac thanh phan uon %------------------------------------------------------ SinhVienKyThuat.Com 258 for intx=1:nog_xb x=pointb(intx,1); % toa do diem Gauss (phuong x) wtx=weightb(intx,1); % trong so diem Gauss (phuong x) for inty=1:nog_yb y=pointb(inty,2); % toa do diem Gauss (phuong y) wty=weightb(inty,2) ; % trong so diem Gauss (phuong y) % tinh cac ham dang va dao ham tai cac diem Gauss [shape,dhdr,dhds]=Shape_Func_4node(x,y); % tinh ma tran Jacobian jacob2=jacob_2D(nnel,dhdr,dhds,xcoord,ycoord); detjacob=det(jacob2); % dinh thuc ma tran Jacobian invjacob=inv(jacob2); % nghich dao ma tran Jacobian % xac dinh dao ham ham dang trong he toa do vat ly for i=1:nnel dhdx(i)=invjacob(1,1)*dhdr(i)+invjacob(1,2)*dhds(i); dhdy(i)=invjacob(2,1)*dhdr(i)+invjacob(2,2)*dhds(i); end % xac dinh ma tran chuyen vi – bien dang [B1, B2]=Bb_matrix_com1(nnel,dhdx,dhdy); B_p= Bp_matrix_com1(nnel,shape,edof); %-------------------------------------------- % tinh ma tran do cung ung voi thanh phan uon %-------------------------------------------- kb=kb+(B1'*A*B1+B1'*B*B2+B2'*B*B1+B2'*D*B2)… *wtx*wty*detjacob; %-------------------------------------------- % tinh vecto luc nut phan tu %-------------------------------------------- f=f+B_p'*wtx*wty*detjacob*p; end end % ket thuc doan chuong trinh tich phan so %------------------------------------------------------ % tinh tich phan so cac thanh phan cat SinhVienKyThuat.Com 259 %------------------------------------------------------ for intx=1:nog_xs x=points(intx,1); wtx=weights(intx,1); for inty=1:nog_ys y=points(inty,2); wty=weights(inty,2) ; % tinh cac ham dang va dao ham tai cac diem Gauss [shape,dhdr,dhds]=Shape_Func_4node(x,y); % tinh ma tran Jacobian jacob2=jacob_2D(nnel,dhdr,dhds,xcoord,ycoord); detjacob=det(jacob2); % dinh thuc ma tran Jacobian invjacob=inv(jacob2); % nghich dao ma tran Jacobian % xac dinh dao ham ham dang trong he toa do vat ly for i=1:nnel dhdx(i)=invjacob(1,1)*dhdr(i)+invjacob(1,2)*dhds(i); dhdy(i)=invjacob(2,1)*dhdr(i)+invjacob(2,2)*dhds(i); end % xac dinh ma tran chuyen vi – bien dang B_s=Bs_matrix_com1(nnel,dhdx,dhdy,shape); %---------------------------------------- % tinh ma tran do cung ung voi thanh phan cat %---------------------------------------- ks=ks+B_s'*A_t*B_s*wtx*wty*detjacob; end end % ket thuc doan chuong trinh tich phan so %-------------------------------- % tinh ma tran do cung phan tu %-------------------------------- k=kb+ks; % xay dung bang ghep noi phan tu index=sys_elm_dof_assoc(nd,nnel,ndof); kk=kk_build_2D(kk,k,index); % ghep noi ma tran do cunng tong SinhVienKyThuat.Com 260 the ff=ff_build_2D(ff,f,index); % ghep noi vecto luc nut tong the end %------------------------------------ % ap dat dieu kien bien %------------------------------------ [kk,ff]=boundary_aply_2D(kk,ff,bcdof,bcval); %--------------------------------- % giai he phuong trinh PTHH %--------------------------------- disp=kk\ff; for node=1:nnode displace(node)=disp((node-1)*5+3); end num=1:1:nnode; Result=[num' displace'] % in ket qua SinhVienKyThuat.Com 261 Các hàm sử dụng trong chương trình function [A, B, D, A_t]=ABD_matrix(emodule_1,emodule_2,emodule_3,... gmodule_12,gmodule_13,gmodule_23,… nuy_12,nuy_13,nuy_23,nol_p,aol_p, z_p) %------------------------------------------------------------------------ % Muc dich: % xac dinh cac ma tran do cung vat lieu composite lop % Cu phap: % [A, B, D]=ABD_matrix(emodule_1,emodule_2,emodule_3,... % gmodule_12,gmodule_13,gmodule_23,… % nuy_12,nuy_13,nuy_23,nol_p,aol_p) %------------------------------------------------------------------------ % Tinh toan cac ma tran do cung cua cac lop nuy_21 = nuy_12*emodule_2/emodule_1; nuy_31 = nuy_13*emodule_3/emodule_1; nuy_32 = nuy_23*emodule_3/emodule_2; temp_delta = 1-nuy_12*nuy_21; U1=(emodule_1+emodule_2+2*nuy_12*emodule_2)/(8*temp_delta); U2=(temp_delta*gmodule_12-nuy_12*emodule_2)/(2*temp_delta); U3=(emodule_1-emodule_2)/(2*temp_delta); U4=(emodule_1+emodule_2-2*nuy_12*emodule_2- 4*temp_delta*gmodule_12)/(8*temp_delta); for k=1:nol_p % C(k,i,j): the k layer, i row, j column of hardness coefficient matrix C(1,1,k)=3*U1+U2+U3*cos(2*aol_p(k))+U4*cos(4*aol_p(k)); C(2,2,k)=3*U1+U2-U3*cos(2*aol_p(k))+U4*cos(4*aol_p(k)); C(1,2,k)=U1-U2-U4*cos(4*aol_p(k)); C(2,1,k)=C(1,2,k); C(3,3,k)=U1+U2-U4*cos(4*aol_p(k)); C(1,3,k)=0.5*U3*sin(2*aol_p(k))+U4*sin(4*aol_p(k)); SinhVienKyThuat.Com 262 C(3,1,k)=C(1,3,k); C(2,3,k)=0.5*U3*sin(2*aol_p(k))-U4*sin(4*aol_p(k)); C(3,2,k)=C(2,3,k); C(4,4,k)=cos(aol_p(k))*cos(aol_p(k))*gmodule_13 + … sin(aol_p(k))*sin(aol_p(k))*gmodule_23; C(5,5,k)=sin(aol_p(k))*sin(aol_p(k))*gmodule_13 + … cos(aol_p(k))*cos(aol_p(k))*gmodule_23; C(4,5,k)=cos(aol_p(k))*sin(aol_p(k))*… (gmodule_13-gmodule_23); C(5,4,k)=C(4,5,k); end for i=1:5 for j=1:5 Temp1(i,j)=0; Temp2(i,j)=0; Temp3(i,j)=0; for k=1:nol_p Temp1(i,j) = Temp1(i,j) + C(i,j,k)*(z_p(k+1)-z_p(k)); Temp2(i,j) = Temp2(i,j) + C(i,j,k)*((z_p(k+1)^2)-(z_p(k)^2))/2; Temp3(i,j) = Temp3(i,j) + C(i,j,k)*((z_p(k+1)^3)-(z_p(k)^3))/3; end end end for i=1:3 for j=1:3 A(i,j) = Temp1(i,j); B(i,j) = Temp2(i,j); D(i,j) = Temp3(i,j); end end for i=1:2 for j=1:2 SinhVienKyThuat.Com 263 A_t(i,j) = Temp1(i+3,j+3)*5.0/6.0; % 5/6 He so hieu chinh end end %-------------------------------------------------------------------------- %-------------------------------------------------------------------------- function [B1,B2]=Bb_matrix_com1(nnel,dhdx,dhdy) %-------------------------------------------------------------------------- % Muc dich: % xac dinh cac ma tran chuyen vi – bien dang % cua cac thanh phan chuyen vi uon va do cong % Cu phap: % [B1, B2]=Bb_matrix_com1(nnel,dhdx,dhdy) % Mo ta cac bien: % nnel – so nut cua phan tu % dhdx – dao ham ham dang theo x % dhdy - dao ham ham dang theo y %-------------------------------------------------------------------------- for row =1:3 for col=1:nnel*5 B1(row,col)=0; B2(row,col)=0; end end for i=1:nnel i1=(i-1)*5+1; i2=i1+1; i3=i1+3; i4=i3+1; B1(1,i1)=dhdx(i); B1(2,i2)=dhdy(i); B1(3,i1)=dhdy(i); B1(3,i2)=dhdx(i); SinhVienKyThuat.Com 264 B2(1,i3)=dhdx(i); B2(2,i4)=dhdy(i); B2(3,i3)=dhdy(i); B2(3,i4)=dhdx(i); end %------------------------------------------------------------------------ function [B_p]=Bp_matrix_com1(nnel,shape,edof) %------------------------------------------------------------------------ % Muc dich: % xac dinh vecto dinh vi luc nut phan tu % Cu phap: % [B_p]=Bp_matrix_com1(nnel,shape, edof) % Mo ta cac bien: % nnel – so luong nut cua phan tu % shape – ham dang % edof - so bac tu do cua phan tu %------------------------------------------------------------------------ B_p=zeros(1,edof); % khoi tao vecto B_p for i=1:nnel ii=(i-1)*5+3; B_p(1,ii)=shape(i); end %------------------------------------------------------------------------ %------------------------------------------------------------------------ function [B_s]=Bs_matrix_com1(nnel,dhdx,dhdy,shape) %------------------------------------------------------------------------ % Muc dich: % xac dinh ma tran quan he chuyen vi – bien dang % cua cac thanh phan bien dang cat SinhVienKyThuat.Com 265 % Cu phap: % [B_s]=Bs_matrix_com1(nnel,dhdx,dhdy,shape) % Mo ta cac bien: % nnel – so nut cua phan tu % dhdx – dao ham ham dang theo x % dhdy - dao ham ham dang theo y % shape - ham dang %------------------------------------------------------------------------ for row =1:2 for col=1:nnel*5 B_s(row,col)=0; end end for i=1:nnel i1=(i-1)*5+3; i2=i1+2; i3=i1+1; B_s(1,i1)=dhdy(i); B_s(1,i2)=shape(i); B_s(2,i1)=dhdx(i); B_s(2,i3)=shape(i); end Kết quả số Result = Nút số Độ võng (mm) 1 0 2 0 3 0 4 0 5 0 SinhVienKyThuat.Com 266 6 0 7 0.0210 8 0.0291 9 0.0210 10 0 11 0 12 0.0256 13 0.0361 14 0.0256 15 0 16 0 17 0.0210 18 0.0291 19 0.0210 20 0 21 0 22 0 23 0 24 0 25 0 SinhVienKyThuat.Com 267 6. BÀI TẬP 12.1-12.3. Phát triển chương trình P12_2 cho bài toán đã cho trong Ví dụ 12.2, nhưng thay liên kết đơn trên 4 cạnh của tấm bằng các dạng sau đây: a. Liên kết ngàm (N) ở một biên y = 0 (TTNT); b. Liên kết ngàm ở biên y = 0 và liên kết tựa bản lề (B) ở biên y = b (TTNB); c. Liên kết ngàm tại 2 biên y = 0 và y =b (TTNN). 12.4-12.6. Phát triển chương trình P12_2 cho bài toán đã cho trong Ví dụ 12.2, với các yêu cầu của bài 12.1-12.3, nhưng thay tải trọng phân bố đều trên bề mặt tấm thành tải trọng phân bố sin như sau:  mm/Np; b ysin a xsinpp 100   12.7. Phát triển chương trình P12_2 để tính toán độ võng lớn nhất và các thành phần ứng suất tại điểm đó của tấm hình vuông có cấu hình lớp vuông xen lớp (0o/90o/0o/90o/...), số các cặp lớp (0o/90o) quy ước là n, lần lượt là n=1, n=5 và n=10. Vật liệu của các lớp là như nhau, có cơ tính là E1/E2 = 25 = 280 gPa, G12 = G13 = 0,5E2, G23 = 0,2 E2, 12=0,25. Xét bài toán với các điều kiện biên khác nhau và chịu tải trọng gây uốn phân bố sin như đã cho trong bài 12.3-12.6. Trong các trường hợp tính toán, các cạnh biên x = 0 và x = a chịu liên kết đơn (BB), hai cạnh còn lại sẽ là (NN), (BB) và (TT). Sử dụng phần tử tứ giác bậc nhất đẳng tham số và khảo sát bài toán với các mật độ lưới khác nhau 33; 44; 66 và 88. Chiều dày các lớp được xác định trong từng trường hợp cụ thể, với quy định tỷ lệ a/h = 100. SinhVienKyThuat.Com 268 Chương 13 PHẦN TỬ HỮU HẠN TRONG BÀI TOÁN ĐỘNG LỰC HỌC KẾT CẤU 1. GIỚI THIỆU Trong các chương trước, chúng ta đã xét các kết cấu chịu tác dụng của tải trọng tĩnh. Trong kỹ thuật, ta còn gặp các kết cấu chịu lực tác dụng tức thời của lực hoặc lực thay đổi theo thời gian, khi ấy khối lượng và gia tốc đóng vai trò quan trọng. Giả sử một kết cấu bị biến dạng đàn hồi, nếu ta loại bỏ lực tác dụng một cách tức thời, khi ấy kết cấu sẽ dao động xung quanh vị trí cân bằng. Chuyển động có chu kỳ này được gọi là dao động tự do. Số chu kỳ trong một đơn vị thời gian được gọi là tần số; chuyển vị lớn nhất từ vị trí cân bằng được gọi là biên độ. Dưới đây ta xét bài toán dao động tự do không có lực cản của vật (kết cấu). 2. MÔ TẢ BÀI TOÁN Người ta định nghĩa Lagrangean L là hiệu của động năng T và thế năng  của vật (hay hệ) như sau: L = T -  (13.1) Nguyên lý Haminton Trong một khoảng thời gian bất kỳ từ t1 đến t2, trạng thái chuyển động của vật sẽ làm cực trị phiếm hàm:  2 1 t t LdtI (13.2) Biểu diễn L theo các biến mở rộng  nn qqqqqq  ,,,,, 2121 ; trong đó t qq ii    . Khi đó, các phương trình chuyển động sẽ được viết dưới dạng: SinhVienKyThuat.Com 269 ni q L q L dt d ii ,,2,1;0              (13.3) Để minh hoạ, chúng ta xét ví dụ sau: Ví dụ Khảo sát hệ gồm hai khối lượng tập trung m1, m2 nối với nhau bởi các lò xo có độ cứng k1 và k2 như Hình 13.1. Động năng và thế năng của hệ được xác định bởi: 2 22 2 11 2 1 2 1 xmxmT    2122 2 11 2 1 2 1 xxkxk  Áp dụng L = T - , ta thu được phương trình chuyển động     0 0 12222 22 1221111 11                       xxkxm x L x L dt d xxkxkxm x L x L dt d     Hoặc dưới dạng ma trận:                                  0 0 0 0 2 1 22 221 2 1 2 1 x x kk kkk x x m m   Hay dưới dạng cô đọng hơn: 0 KXXM  (13.4) Trong đó M là ma trận của các khối lượng của các phần tử hữu hạn của hệ, K là ma trận độ cứng, X là ma trận chuyển vị, X là ma trận gia tốc. 11, xx  22 , xx  k1 k2 m2 m1 Hình 13.1. Hệ lò xo-khối lượng SinhVienKyThuat.Com 270 3. VẬT RẮN CÓ KHỐI LƯỢNG PHÂN BỐ Khảo sát một vật rắn với khối lượng phân bố (Hình 13.2). Biểu thức xác định thế năng đã trình bày trong (1.9) - Chương 1, còn động năng được xác định bởi:  V T dVuuT  2 1 (13.5) Trong đó  là khối lượng riêng của vật liệu; u là véctơ vận tốc tại điểm x với các thành phần u , v và w .  Twvuu   (13.6) Theo phương pháp PTHH, vật thể được chia thành các phần tử; trong mỗi phần tử, ta biểu diễn chuyển vị u theo các chuyển vị nút q nhờ các hàm dạng N: u = N q (13.7) Trong các bài toán động lực học, các thành phần chuyển vị q phụ thuộc vào thời gian. Vì vậy, véctơ vận tốc được xác định bởi: qNu   (13.8) Thay (13.8) vào (13.5) ta thu được động năng của phần tử: qdVNNqT e TT          2 1 (13.9) Biểu thức trong ngoặc được gọi là ma trận khối lượng của phần tử:  e Te dVNNm  (13.10) y x z uu , vv , ww , dV V Hình 13.2. Vật có khối lượng phân bố SinhVienKyThuat.Com 271 Động năng của cả vật (hệ) được xác định bởi: QMQqmqTT T e eT e e  2 1 2 1   (13.11) Mặt khác, thế năng của vật (hệ) được xác định bởi: FQKQQ TT  2 1 (13.12) Áp dụng hệ thức (13.3), ta thu được phương trình chuyển động FKQQM  (13.13) Trong đó M là ma trận khối lượng của các phần tử hữu hạn của hệ do các khối lượng phân bố qui đổi ra. Trong trường hợp vật dao động tự do, F = 0, do đó: 0 KQQM  (13.14) Khi dao động tự do, tất cả các điểm của hệ dao động cùng pha, vì vậy ta có thể viết: Q = U sin t (13.15) trong đó U là véctơ dao động nút và  (rad/s) là tần số góc. Thế phương trình (13.15) vào (13.14), ta được: KU = 2 M U (13.16) Đây là bài toán trị riêng tổng quát. Ta cũng có thể viết (13.16) dưới dạng: KU =  M U (13.17) Ở đây, U là véctơ riêng, biểu thị dạng dao động ứng với trị riêng . Trị riêng  là bình phương của tần số góc , còn tần số f = /2 được đo bằng Hz (số chu kỳ/giây). Mỗi giá trị tần số  ứng với một dạng dao động cụ thể. Có nhiều phương pháp giải hệ phương trình (13.17) để tìm các trị riêng và dạng dao động của hệ (xem thêm giáo trình phương pháp tính). Trong thực hành, người ta hay dùng các chương trình để giải bài toán trị riêng. SinhVienKyThuat.Com 272 4. MA TRẬN KHỐI LƯỢNG CỦA PHẦN TỬ CÓ KHỐI LƯỢNG PHÂN BỐ Xét một vật liệu có khối lượng riêng  bằng hằng số. Từ công thức (13.10) ta có:  e Te dVNNm  (13.18) 4.1. Phần tử một chiều Với phần tử một chiều như mô tả trong Hình 13.3. Ta có:  Tqqq 21        2 1 2 1 21 NNN (13.19) Khi đó:          21 12 62 1 1 eeTee e e Te lANdNlAdxANNm  (13.20) 4.2. Phần tử trong hệ thanh phẳng Phần tử giàn, như mô tả trong Hình 13.4. Ta có: 1 q1 q2 le Hình 13. 3. Phần tử 1 chiều dV=Adx 2 1 Hình 13.4. Phần tử giàn 2 q3 q4 q1 q2 u v SinhVienKyThuat.Com 273             21 21 4321 00 00 ; NN NN N qqqqqvuu TT (13.21) Trong đó N1, N2 cũng được xác định theo (13.19) Tương tự như trên, ma trận khối lượng của phần tử giàn cũng được xác định theo biểu thức:              2010 0201 1020 0102 6 eee lAm  (13.22) 4.3. Phần tử tam giác Với phần tử hữu hạn tam giác biến dạng hằng số, Hình 13.5. ta có:                311 321 321 654321 ;;1 000 000 ; NNN NNN NNN N qqqqqqq vuu T T Chú ý rằng: ; 12 1; 6 1 21 2 1 e e e e AdANNAdAN   Khi đó, ma trận khối lượng của phần tử được xác định bởi:                      201010 020101 102010 010201 101020 010102 12 eee tAm  (13.23) 1 q1 q2 2 q3 q4 3 q5 q6 u v x y Hình 13.5. Phần tử tam giác phẳng SinhVienKyThuat.Com 274 4.4. Phần tử tam giác đối xứng trục Với phần tử tam giác đối xứng trục, ta có  Twuu  Trong đó u và w là các chuyển vị hướng kính và hướng trục. Véctơ q và N được xác định tương tự như trường hợp phần tử tam giác ở trên. Khi ấy:   e T e Te dArNNdVNNm  2 (13.24) Vì: r = N1r1 + N2r2 + N3r3, do đó:    e Te dANNrNrNrNm 3322112 Chú ý ; 60 ; 30 ; 10 3212 2 1 3 1 e e e e e e AdANNNAdANNAdAN   Cuối cùng ta được                                      rrrrrr rrrrrr rrrrrr rrrrrr rrrrrr rrrrrr Am ee 2 3 40 3 20 3 20 02 3 40 3 20 3 2 3 202 3 40 3 20 0 3 202 3 40 3 2 3 20 3 202 3 40 0 3 20 3 202 3 4 10 3 12 3 12 1 2 3 1 2 3 23 1 23 1  (13.25) Trong đó: 3 321 rrrr  (13.26) SinhVienKyThuat.Com 275 4.5. Phần tử tứ giác Với phần tử tứ giác ở trạng thái ứng suất và biến dạng phẳng:             4321 4321 87654321 0000 0000 ; NNNN NNNN N qqqqqqqqqvuu TT (13.27) Ma trận khối lượng được xác định bởi:      1 1 1 1 det  ddJNNtm Te e (13.28) Để tính được ma trận khối lượng trên, ta sẽ áp dụng công thức tích phân số (xem Chương 8) 4.6. Phần tử dầm Với phần tử dầm Hermite, Hình 13.6, ta có: v = Hq (13.29) Khi ấy:    1 1 2  dlAHHm ee Te , sau khi lấy tích phân, ta được:                  eeee ee eeee ee eee llll ll llll ll lAm 422313 221561354 313422 135422156 420  (13.30) q2 q1 q3 q4 v le Hình 13.6. Phần tử dầm SinhVienKyThuat.Com 276 4.7. Phần tử khung Trong hệ toạ độ địa phương (x', y'), ma trận khối lượng của phần tử được xem như một tổ hợp của ma trận khối lượng của phần tử thanh và phần tử dầm. Khi ấy, ma trận khối lượng được xác định bởi:                         blblblbl blbblb aa blblblbl blbblb aa m eeee ee eeee ee e 22 222 2 42203130 22156013540 00200 31304220 13540221560 00002 ' (13.31) với ký hiệu: 420 ; 6 eeee lAblAa   Áp dụng ma trận chuyển vị L đã giới thiệu trong Chương 9, ta sẽ xác định được ma trận khối lượng me của phần tử khung trong hệ toạ độ chung bởi: LmLm eTe ' (13.32) 5. VÍ DỤ Cho một dầm ngàm hai đầu (Hình 13.7), chiều dài 2, khối lượng riêng , mặt cắt ngang tròn, diện tích A. Dầm chịu dao động uốn tự do. Hãy xác định các tần số dao động riêng. Lời giải Chia dầm ra hai phần tử có chiều dài bằng nhau. Từ ma trận độ cứng của phần tử chịu uốn (xem Chương 9), ta xác định được ma trận độ cứng chung K; chú ý đến điều kiện biên sẽ được: 1 2 1 2 3 l1=l l2=l A B Hình 13.7 SinhVienKyThuat.Com 277        23 80 024 ll EJK Áp dụng công thức (13.31), ta xác định được ma trận khối lượng của các phần tử, từ đó suy ra ma trận khối lượng chung M; chú ý đến điều kiện biên sẽ được:        280 0312 420 l AlM  Cuối cùng ta thiết lập được phương trình                          4 3 2 2 4 3 23 80 0312 42080 024 U U l Al U U ll EJ   do tính đối xứng, ta suy ra ngay: 4 2 24 2 1 420;31,32 Al EJ Al EJ      Nếu chia dầm ra nhiều phần tử, ta sẽ tính được các giá trị tần số góc chính xác hơn. 6. CHƯƠNG TRÌNH TÍNH TẦN SỐ DAO ĐỘNG TỰ DO CỦA DẦM VÀ KHUNG 6.1. Chương trình tính tần số dao động tự do của dầm Tính tần số dao động tự do của dầm công xôn như Hình 13.8. Dầm có chiều dài 1m, tiết diện mặt cắt ngang 2020 mm; khối lượng riêng vật liệu dầm là 1000Kg/m3; môđul đàn hồi kéo là 100 gPa, môđul đàn hồi trượt của vật liệu dầm là 40 gPa. Ở đây, ta sẽ minh hoạ bằng chương trình Matlab, với số phần tử là 4. Trong chương trình này, chúng ta sử dụng hàm  = eig(K,M) là một công cụ sẵn có của Matlab, cho phép tính nghiệm riêng của phương trình KU = MU. Ở đây, U là véctơ riêng, biểu thị dạng dao động ứng với trị riêng . Trị riêng  là bình phương của tần số góc , còn tần số f = /2 được đo bằng Hz (số chu kỳ/giây). Mỗi giá trị tần số  ứng với một dạng dao động cụ thể. SinhVienKyThuat.Com 278 Chương trình nguồn %---------------------------------------------------------------------------- % Chuong trinh so 1, chuong 13- Vi du 13.1 (P13_1) %---------------------------------------------------------------------------- % tinh tan so dao dong tu do cua dam, su dung phan tu dam voi 1 bac tu do % Mo ta bai toan % Tinh tan so dao dong tu do cua dam cong xon, su dung 4 phan tu. % dam co chieu dai 1 m, kich thuoc mat cat ngang 0.02 x0.02 (m) % va trong luong rieng la 1000 Kg/m^3. % mo dul dan hoi E = 100 GPa va mo dul dan hoi truot la 40 GPa % Mo ta cac bien % 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 chua cac bac tu do bi rang buoc theo dieu kien bien %----------------------------------------------------------------------------------- -- clear noe=4; % tong so phan tu nnel=2; % so nut cua phan tu ndof=3; % so bac tu do tai nut edof = nnel*ndof; % so bac tu do cua phan tu nnode=noe+1; % tong so nut cua he sdof=nnode*ndof; % so bac tu do cua he e_module=100*10^9; % modul dan hoi 1 2 3 4 1 2 3 4 5 1 m A-A A-A 0,02 Hình 13.8. Mô hình dầm dao động uốn tự do SinhVienKyThuat.Com 279 g_module=40*10^9; % modul dan hoi truot tleng=1; % chieu dai dam leng=tleng/noe; % chieu dai cac phan tu (deu nhau) h=0.02; % chieu cao mat cat dam b=0.02; % chieu rong mạt cat dam rho=1000; % trong luong rieng vat lieu dam bcdof(1)=1; % bac tu do thu 1 tai nut 1 bi rang buoc theo dieu kien bien bcdof(2)=2; % bac tu do thu 2 tai nut 1 bi rang buoc theo dieu kien bien bcdof(3)=3; % bac tu do thu 3 tai nut 1 bi rang buoc theo dieu kien bien kk=zeros(sdof,sdof); % khoi tao ma tran do cung tong the mm=zeros(sdof,sdof); % khoi tao ma tran khoi luong tong the index=zeros(edof,1); % khoi tao bang ghep noi for iel=1:noe % xet tuing phan tu trong he % xay dung bang ghep noi start = (iel-1)*(nnel-1)*ndof; for i=1:edof index(i)=start+i; end % tinh ma tran do cung phan tu va ma tran khoi luong phan tu [k,m]=beam_elm_3(e_module,g_module,leng,h,b,rho); % ghep noi ma tran do cung kk=kk_build_2D(kk,k,index); % ghep noi ma tran khoi luong mm=kk_build_2D(mm,m,index); end %---------------------------------- % ap dat dieu kien bien %---------------------------------- [kn,mn]=boundary_aply_beam(kk,mm,bcdof); fsol=eig(kn,mn); % giai he phuong trinh tri rieng va in ket qua fsol=sqrt(fsol) SinhVienKyThuat.Com 280 Các hàm sử dụng trong chương trình function [k,m]=beam_elm_3(e_module,g_module,leng,h,b,rho) %-------------------------------------------------------------- % Muc dich: % Xac dinh ma tran do cung va ma tran khoi luong phan tu % cua phan tu dam, voi cac chuyen vi nut chi la chuyen vi dai % vecto chuyen vi phan tu: {u_1_b u_1_t v_1 u_2_b u_2_t v_2} % Cu phap: % [k,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