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

1. Công thức tích phân 1- điểm Khi n = 1:    11 1 1  fwdfI    (8.33) Có 2 tham số w1 và 1 cần xác định, với yêu cầu là (8.33) cho kết quả chính xác khi f() là một đa thức bậc nhất. Điều này có nghĩa là nếu f() là một đa thức bậc nhất, thì sai số:     011 1 1 10     fwdaa (8.34a) hoặc   02 11010  aawa (8.34b)   02 11110  awwa (8.34c) Từ (8.34c), suy ra  = 0 khi:      0 2 1 1  w (8.35) Với hàm f tổng quát và tuỳ ý, ta có:    02 1 1 fdfI     (8.36) Kết quả tính giống như phương pháp trung điểm (phương pháp hình chữ nhật), như mô tả ở Hình 8.2 sau. 6.1.2. Công thức tích phân 2- điểm  f(0) f -1 0 1 Phần diện tích xấp xỷ = 2f(0) Phần diện tích chính xác f() Hình 8.2. Cầu phương 1 điểm Gauss SinhVienKyThuat.Com 134 Khi n = 2      2211 1 1  fwfwdfI    (8.37) Ở đây, có 4 thông số cần xác định là: w1; w2 ; 1 và 2. Chọn   332210  aaaaf  , ta có:        02211 1 1 3 3 2 210     fwfwdaaaa (8.38) Để thoả mãn (8.38):  = 0, khi đó ta phải có:             0 3 2 0 2 3 22 3 11 2 22 2 11 2211 21    ww ww ww ww (8.39) Hệ phương trình phi tuyến (8.39) có nghiệm duy nhất là:       57735,0 3 1 1 21 21  ww (8.40) Từ đây, chúng ta có thể kết luận rằng phương pháp cầu phương Gauss n-điểm sẽ cho kết quả chính xác nếu f là một đa thức bậc (2n-1) hoặc nhỏ hơn. Bảng dưới cho các giá trị của wi và i theo công thức Gauss với n =1,...,6. Bảng 8.1. Điểm Gauss và hàm trọng lượng n i wi n i wi 1 0 2 5 0,9061798459 0,2369268851 2 0,5773502692 1 0,5384693101 0,4786286705 3 0,7745966692 0,5555555556 0 0,5688888889 0 0, 888888889 6 0,9324695142 0,1713244924 4 0,8611363116 0,3478548451 0,6612093865 0,3607615730 0,3399810436 0,6521451549 0,2386191861 0,4679139346 6.2. Tích phân số hai biến SinhVienKyThuat.Com 135 Xét tích phân    ddfI      1 1 1 1 , (8.41) Ta có           n i n j jiji n n jiij n ii fwwfwwdfwI 1 11 1 1 1 1 ,,,  (8.42) 6.3. Tích phân ma trận độ cứng Từ công thức (8.26):      1 1 1 1 det  ddJDBBtk Te e Chú ý ke là ma trận độ cứng (8x8) đối xứng, do đó ta chỉ cần lấy tích phân các số hạng ở phía trên đường chéo chính là đủ. Ký hiệu    ijTe JDBBt det,   (8.43) và áp dụng công thức (8.42), ta được        22 2 21212212111 2 1 ,,,,   wwwwwwkij (8.44) Trong đó:      57735,0;57735,0 1 2211 21  ww Các điểm Gauss theo công thức tích phân 2- điểm ở trên được mô tả trên Hình 8.3. Nếu gọi các điểm Gauss lần lượt là 1, 2, 3, 4 thì kij trong (8.44) được viết dưới dạng:    4 1K Kkji Wk (8.45) Trong đó K là giá trị của  và WK là hàm trọng số tại điểm tích phân k. Chú ý rằng WK = (1)(1) = 1. Công thức (8.45) rất thuận lợi cho chúng ta khi lập trình tính toán ma trận độ cứng. SinhVienKyThuat.Com 136 7. TÍNH ỨNG SUẤT Không giống như phần tử tam giác có biến dạng hằng số; với phần tử tứ giác, ứng suất  = DBq không phải là hằng số trong phạm vi một phần tử, mà chúng thay đổi theo  và . Trong thực tế, ứng suất được tính tại các điểm Gauss, đây cũng là những điểm được sử dụng trong khi tính ma trận độ cứng ke. Như vậy, với phần tử tứ giác, ta sẽ có 4 giá trị ứng suất. Để tránh sinh ra nhiều số liệu, ta chỉ cần tính ứng suất tại một điểm cho mỗi phần tử, chẳng hạn điểm ( = 0, = 0). 8. VÍ DỤ Khảo sát một phần tử chữ nhật như hình vẽ dưới đây. Giả thiết phần tử chịu trạng thái ứng suất phẳng và có các thành phần chuyển vị:    mmq T00081,0152,0076,0051,000 . Biết E = 206 gPa,  = 0,3. Hãy xác định các đại lượng J, B và  tại điểm ( = 0,  = 0). 1 1 2 3 4   3 1  3 1  3 1  3 1  Hình 8.3. Điểm Gauss theo qui tắc tích phân 2 điểm 3 (2,1) y x 2 (2,0) 1 (0,0) 4 (0,1) q7 q8 q5 q6 q3 q4 q1 q2 SinhVienKyThuat.Com 137 Bài giải Theo công thức (8.10), ta có:                                  2 10 01 111212 111212 4 1   J Áp dụng công thức (8.19), ta được:                0 2 110 1000 000 2 1 21 1A Áp dụng công thức (8.21) và (8.23) với  = 0,  = 0 ta sẽ được:                     4 1 2 1 4 1 2 1 4 1 2 1 4 1 2 1 2 10 2 10 2 10 2 10 0 4 10 4 10 4 10 4 1 B Ma trận D được xác định bởi (6.8)              35,000 013,0 03,01 09,01 10206 3D Cuối cùng ta tính được  tại điểm  = 0,  = 0:    MPaDBq T282159461 SinhVienKyThuat.Com 138 9. CHƯƠNG TRÌNH Trở lại với ví dụ ở mục 7, chương 6. Ở đây, ta sẽ sử dụng phần tử tứ giác bậc nhất đẳng tham số với 4 phần tử để xây dựng chương trình tính cho kết cấu, như mô tả trong Hình 8.5. Trong chương trình, việc lấy tích phân số sẽ được thực hiện với 22 điểm Gauss. Theo sơ đồ lưới PTHH, điều kiện biên sẽ được mô tả như sau: Tại nút 1, các thành phần chuyển vị u và v bằng không, tại nút 2 u bằng không; như vậy điều kiện biên tương ứng là Q1 = Q2 = Q3 = 0 Chương trình nguồn. %---------------------------------------------------------------------------- % Chuong trinh so 1, chuong 8 - (P8_1) %---------------------------------------------------------------------------- % Phan tich ung suat cua tam chiu keo su dung phan tu tu giac bac 1 % (Hinh. 8.5 mo ta luoi phan tu) % Mo ta cac bien % k = ma tran do cung phan tu % f = vecto luc nut phan tu % kk = ma tran do cung tong the % ff = vecto luc nut tong the % disp = vecto chuyen vi nut tong the % eldisp = vecto chuyen vi nut phan tu % stress = ma tran cac thanh phan ung suat % strain = ma tran cac thanh phan bien dang % gcoord = toa do nut Hình 8.5. Tấm chữ nhật chịu kéo trong mặt phẳng 1 2 2 1 3 4 4 3 6 5 8 7 10 9 p SinhVienKyThuat.Com 139 % nodes = ma tran chi so nut phan tu % index = vecto chua cac bac tu do (chuyen vi) tong the tai moi phan tu % bcdof = vecto chua cac bac tu do (chuyen vi nut) tai bien cua ket cau % bcval = vecto chi dieu kien bien (gia tri chuyen vi nut tai bien) %---------------------------------------------------------------------------- %------------------------------------ % Cac tham so dau vao %------------------------------------ clear noe_x=4; % so luong cac phan tu theo phuong x noe_y=1; % so luong cac phan tu theo phuong y noe=noe_x*noe_y; % tong so phan tu nnel=4; % tong so nut cua moi phan tu ndof=2; % tong so bac tu do tai moi nut nnode=(noe_x+1)*(noe_y+1); %tong so nut cua he sdof=nnode*ndof; % tong so bac tu do cua he edof=nnel*ndof; % so bac tu cua phan tu E_module=200e3; % modul dan hoi khi keo nen poisson=0.3; % he so Poisson x_length=400; % chieu dai tam (mm) y_length=100; % chieu rong tam (mm) t=10; % chieu day tam (mm) p=200; % tai trong phan bo deu (MPa=N/mm) nog_x=2; nog_y=2; % 2x2 Gauss-Legendre quadrature nog_xy=nog_x*nog_y; % number of sampling points per element %--------------------------------------------- % nhap gia tri toa do cac nut: gcoord(i,j), % trong do: i la chi so nut, j =1,2 chi toa do x hay y %--------------------------------------------- len_x_elm = x_length/noe_x; len_y_elm = y_length/noe_y; for col_index=1:(noe_x+1) for row_index=1:(noe_y+1) SinhVienKyThuat.Com 140 gcoord((row_index+(col_index-1)*(noe_y+1)),1)=... (col_index-1)*len_x_elm; gcoord((row_index+(col_index-1)*(noe_y+1)),2)=... (row_index-1)*len_y_elm; end end gcoord %--------------------------------------------------------- % Nhap chi so nut o moi phan tu: nodes(i,j) % trong do: i la chi so cua phan tu va j la chi so nut %--------------------------------------------------------- nodes=[1 3 4 2; 3 5 6 4; 5 7 8 6; 7 9 10 8]; %------------------------------------- % Nhap cac dieu kien bien %------------------------------------- bcdof=[1 2 3]; % bon bac tu do dau tien chiu rang buoc bcval=[0 0 0]; % voi gia tri bang khong %----------------------------------------- % khoi tao cac ma tran va cac vecto %----------------------------------------- ff=zeros(sdof,1); % vecto luc nut tong the kk=zeros(sdof,sdof); % ma tran do cung tong the disp=zeros(sdof,1); % vecto chuyen vi nut tong the eldisp=zeros(edof,1); % vecto chuyen vi nut phan tu stress=zeros(noe,3); % ma tran ung suat strain=zeros(noe,3); % ma tran bien dang index=zeros(edof,1); % veto chi so B_matrix=zeros(3,edof); % ma tran chuyen vi-bien dang D_matrix=zeros(3,3); % ma tran do cung vat lieu %---------------------------- % vecto luc nut %---------------------------- ff(17)=(y_length/noe_y)*p/2; % luc dat tai nut so 9 theo phuong x ff(19)=(y_length/noe_y)*p/2; % luc dat tai nut so 10 theo phuong x %----------------------------------------------------------------- SinhVienKyThuat.Com 141 % tinh toan cac ma tran do cung phan tu, vec to luc nut phan tu va % ghep noi thanh ma tran do cung tong the, vecto luc nut tong the %----------------------------------------------------------------- % tinh toa do va trong so cac diem Gauss [point,weight]=Gauss_Point_2D(nog_x,nog_y); % tinh ma tran do cung vat lieu D_matrix=D_matrix_2D(1,E_module,poisson); for iel=1:noe % xet tung phan tu for i=1:nnel nd(i)=nodes(iel,i); % chi so nut tong the cua nut thu i phan tu iel xcoord(i)=gcoord(nd(i),1); % toa do x cua nut thu i ycoord(i)=gcoord(nd(i),2); % toa do y cua nut thu i end k=zeros(edof,edof); % khoi tao ma tran do cung phan tu iel % tich phan so for intx=1:nog_x x=point(intx,1); % toa do x cua diem Gauss wtx=weight(intx,1); % trong so diem gauss for inty=1:nog_y y=point(inty,2); % toa do y cua diem Gauss wty=weight(inty,2) ; % trong so diem gauss % tinh cac ham dang va cac dao ham tai cac diem Gauss [shape,dhdr,dhds]=Shape_Func_4node(x,y); % 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 cac 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 SinhVienKyThuat.Com 142 % ----------------------------------- % tinh ma tran chuyen vi-bien dang B % ----------------------------------- B_matrix=B_matrix_2D(nnel,dhdx,dhdy); %------------------------------ % tinh ma tran do cung phan tu %------------------------------ k=k+t*B_matrix'*D_matrix*B_matrix*wtx*wty*detjacob; end end % ket thuc doan chuong trich tinh phan so % xay dung vecto chi so bac tu do tong the gan ket voi tung phan tu index=sys_elm_dof_assoc(nd,nnel,ndof); % ghep noi ma tran do cung tong the kk=kk_build_2D(kk,k,index); end %----------------------------- % ap dat dieu kien bien %----------------------------- [kk,ff]=boundary_aply_2D(kk,ff,bcdof,bcval); %------------------------------ % tinh vecto chuyen vi nut tong the %------------------------------ disp=kk\ff; %------------------------------ % in ket qua chuyen vi %------------------------------ num=1:1:sdof; displace=[num' disp] %------------------------------ % tinh ung suat phan tu %------------------------------ for ielp=1:noe % xet tung phan tu for i=1:nnel nd(i)=nodes(ielp,i); % chi so nut tong the cua nut thu i phan tu iel xcoord(i)=gcoord(nd(i),1); % toa do x cua nut thu i SinhVienKyThuat.Com 143 ycoord(i)=gcoord(nd(i),2); % toa do y cua nut thu i end %-------------------------------- % tich phan so %-------------------------------- intp=0; for intx=1:nog_x x=point(intx,1); % toa do x diem gauss wtx=weight(intx,1); % trong so diem gauss for inty=1:nog_y y=point(inty,2); % stoa do y diem gauss wty=weight(inty,2) ; % trong so diem gauss intp=intp+1; % tinh cac ham dang va cac dao ham tai cac diem Gauss [shape,dhdr,dhds]=Shape_Func_4node(x,y); % ma tran Jacobian jacob2=jacob_2D(nnel,dhdr,dhds,xcoord,ycoord); detjacob=det(jacob2); % determinant of Jacobian invjacob=inv(jacob2); % inverse of Jacobian matrix % xac dinh cac 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 % ----------------------------------- % tinh ma tran chuyen vi-bien dang B_p % ----------------------------------- B_matrix=B_matrix_2D(nnel,dhdx,dhdy); % xay dung vecto chi so bac tu do tong the gan ket voi tung phan tu index=sys_elm_dof_assoc(nd,nnel,ndof); %------------------------------------------------------- % xac dinh vec to chuyen vi nut phan tu %------------------------------------------------------- for i=1:edof eldisp(i)=disp(index(i)); end SinhVienKyThuat.Com 144 estrain=B_matrix*eldisp; % tinh bien dang estress=D_matrix*estrain; % tinh ung suat for i=1:3 strain(intp,i)=estrain(i); % luu kq cua moi phan tu stress(intp,i)=estress(i); % luu kq cua moi phan tu end location=[ielp,intx,inty] % in vi tri bieu dien US stress(intp,:) % gia tri ung suat end end % ket thuc phan tinh tich phan so end SinhVienKyThuat.Com 145 Các hàm sử dụng trong chương trình function [point_1D,weight_1D]=Gauss_Point_1D(nog_1D) %------------------------------------------------------------------- % Muc dich: % xac dinh các diem tich phan va ham trong so % cua cac diem Gauss-Legendre voi tich phan so 1 bien % Cu phap: % [point1,weight1]=1_D_Gauss_Point(nog_1D) % Mo ta cac bien: % nog_1D – so luong cac diem tich phan % point_1D - vector chua tao do cac diem Gauss % weight_1D - vector chua cac he so trong so cua diem Gauss %------------------------------------------------------------------- % khoi tao point_1D=zeros(nog_1D,1); weight_1D=zeros(nog_1D,1); % tim cac diem tich phan va he so trong so tuong ung if nog_1D==1 % tich phan 1 diem point_1D(1)=0.0; weight_1D(1)=2.0; elseif nog_1D==2 % tich phan 2 diem point_1D(1)=-0.577350269189626; point_1D(2)=-point_1D(1); weight_1D(1)=1.0; weight_1D(2)=weight_1D(1); elseif nog_1D==3 % tich phan 3 diem point_1D(1)=-0.774596669241483; point_1D(2)=0.0; point_1D(3)=-point_1D(1); weight_1D(1)=0.555555555555556; weight_1D(2)=0.888888888888889; weight_1D(3)=weight_1D(1); SinhVienKyThuat.Com 146 elseif nog_1D==4 % tich phan 4 diem point_1D(1)=-0.861136311594053; point_1D(2)=-0.339981043584856; point_1D(3)=-point_1D(2); point_1D(4)=-point_1D(1); weight_1D(1)=0.347854845137454; weight_1D(2)=0.652145154862546; weight_1D(3)=weight_1D(2); weight_1D(4)=weight_1D(1); else % tich phan 5 diem point_1D(1)=-0.906179845938664; point_1D(2)=-0.538469310105683; point_1D(3)=0.0; point_1D(4)=-point_1D(2); point_1D(5)=-point_1D(1); weight_1D(1)=0.236926885056189; weight_1D(2)=0.478628670499366; weight_1D(3)=0.568888888888889; weight_1D(4)=weight_1D(2); weight_1D(5)=weight_1D(1); end %------------------------------------------------------------------- function [point_2D,weight_2D]=Gauss_Point_2D(nog_xa,nog_ya) %------------------------------------------------------------------- % Muc dich: % xac dinh các diem tich phan va ham trong so % cua cac diem Gauss-Legendre voi tich phan so 2 bien % Cu phap: % [point_2D,weight_2D]=2_D_Gauss_Point(nog_xa,nog_ya) % Mo ta cac bien: % nog_xa - so diem tich phan theo phuong x % nog_ya - so diem tich phan theo phuong y % point_2D - vector chua cac diem Gauss SinhVienKyThuat.Com 147 % weight_2D - vector cac he so trong so cua diem Gauss %------------------------------------------------------------------- % xac dinh gia tri max giua nog_xa va nog_ya if nog_xa > nog_ya nog=nog_xa; else nog=nog_ya; end % khoi tao point_2D=zeros(nog,2); weight_2D=zeros(nog,2); % tim cac diem tich phan va he so trong so tuong ung % tinh cac diem Gauss theo phuong x [point_xa,weight_xa]=Gauss_Point_1D(nog_xa); % tinh cac diem Gauss theo phuong x [point_ya,weight_ya]=Gauss_Point_1D(nog_ya); % tinh cac diem Gauss theo 2 phuong for intx=1:nog_xa % quadrature in x-axis point_2D(intx,1)=point_xa(intx); weight_2D(intx,1)=weight_xa(intx); end for inty=1:nog_ya % quadrature in y-axis point_2D(inty,2)=point_ya(inty); weight_2D(inty,2)=weight_ya(inty); end %------------------------------------------------------------------------ function [shape,dos_r,dos_s]=Shape_Func_4node(rvalue,svalue) %------------------------------------------------------------------------ % Muc dich: % tinh ham dang cua phan tu tu giac bac nhat dang tham so % va cac dao ham cua no tai cac diem Gauss trong he truc tu nhien % Cu phap: % [shape,dos_r,dos_s]=4node_Shape_Func(rvalue,svalue) % Mo ta cac bien: SinhVienKyThuat.Com 148 % shape – ham dang cua phan tu tu giac bac nhat % dos_r – dao ham ham dang theo r (Coxi) % dos_s - dao ham ham dang theo s (Nheta) % rvalue – gia tri toa do r cua diem Gauss % svalue - gia tri toa do s cua diem Gauss % Chu y: % toa do nut 1 la (-1,-1), toa do nut 2 la (1,-1) % toa do nut 3 la (1,1), toa do nut 4 la (-1,1) %------------------------------------------------------------------------ % Ham dang shape(1)=0.25*(1-rvalue)*(1-svalue); shape(2)=0.25*(1+rvalue)*(1-svalue); shape(3)=0.25*(1+rvalue)*(1+svalue); shape(4)=0.25*(1-rvalue)*(1+svalue); % cac dao ham dos_r(1)=-0.25*(1-svalue); dos_r(2)=0.25*(1-svalue); dos_r(3)=0.25*(1+svalue); dos_r(4)=-0.25*(1+svalue); dos_s(1)=-0.25*(1-rvalue); dos_s(2)=-0.25*(1+rvalue); dos_s(3)=0.25*(1+rvalue); dos_s(4)=0.25*(1-rvalue); %-------------------------------------------------------------------------- SinhVienKyThuat.Com 149 Kết quả số - Chuyển vị nút: Bậc tự do (q) Giá trị Bậc tự do (q) Giá trị Bậc tự do (q) Giá trị Bậc tự do (q) Giá trị 1 0 6 0.0000 11 0.0200 16 -0.0030 2 0 7 0.0100 12 -0.0030 17 0.0400 3 0 8 -0.0030 13 0.0300 18 0.0000 4 - 0.0030 9 0.0200 14 0.0000 19 0.0400 5 0.0100 10 0.0000 15 0.0300 20 -0.0030 - Ứng suất tại các điểm Gauss: Điểm Gauss Phần tử 1, 1 1, 2 2, 1 2, 2 Phần tử 1 20.0000 20.0000 20.0000 20.0000 Phần tử 2 20.0000 20.0000 20.0000 20.0000 Phần tử 3 20.0000 20.0000 20.0000 20.0000 Phần tử 4 20.0000 20.0000 20.0000 20.0000 Chú ý: Thành phần ứng suất trong bảng là x, các thành phần ứng suất còn lại y và xy đều bằng không. Kết quả này hoàn toàn tương tự như đã nhận được trong Ví dụ 6.8.1 ở Chương 6. SinhVienKyThuat.Com 150 10. BÀI TẬP 8.1. Trở lại kết cấu trong bài toán 6.1 trên hình 8.10.1. Ta vẫn đi xác định chuyển vị tại hai điểm A và B, và phân bố của ứng suất trong tấm, nhưng ở đây ta sử dụng các phần tử chữ nhật. Vật liệu của tấm là thép có môđun đàn hồi E=210 gPa và hệ số Poisson = 0,22. Hãy chia tấm ra thành hai phần tử chữ nhật theo hai cách được mô tả trong các hình: Hình 8.10.1a,b. So sánh kết quả trong hai trường hợp trên và so sánh với kết quả tìm được bằng các phần tử tam giác như trong bài 6.1. Biết: a = 0,2 m; b = 0,3 m; t = 5mm và q = 8 kN/m. 8.2. Một kết cấu được mô tả như Hình 8.10.2. Bằng các phần tử chữ nhật tuyến tính, hãy giải lại bài toán theo yêu cầu như bài toán 6.4. c 2b 2b a a b P 2b b a a Hình 8.10.2. Kết cấu lắp ghép a = 50mm; b = 75mm; c = 500mm P = 8kN. 4 3 2 1 b a B A q a b q B 2 1 4 3 A 1 2 1 2 Hình 8.10.1. Sơ đồ hoá PTHH với phần tử tứ giác SinhVienKyThuat.Com 151 8.3. Xác định ma trận độ cứng tổng thể các lưới phần tử như trên hình vẽ 8.10.3. Giả thiết mỗi nút có một bậc tự do và gọi [K(e)] là ma trận độ cứng của phần tử thứ e. (Kết quả có thể biểu diễn thông qua các số hạng Kij(e)). 8.4. Trở lại kết cấu trong bài tập 6.5 (Hình 6.9.5). Bằng cách kết hợp hợp lý các phần tử tam giác tuyến tính với các phần tử chữ nhật tuyến tính, hãy tính lại hệ số tập trung ứng suất trong các kết cấu. 1 2 3 4 5 6 7 8 1 2 3 4 1 2 3 4 5 1 2 3 4 5 6 7 8 9 10 Hình 8.10.3. Sơ đồ lưới phần tử SinhVienKyThuat.Com 152 Chương 9 PHẦN TỬ HỮU HẠN TRONG TÍNH TOÁN KẾT CẤU DẦM VÀ KHUNG 1. GIỚI THIỆU Dầm và khung được ứng dụng rộng rãi trong kỹ thuật. Trong chương này chúng ta sẽ áp dụng phương pháp phần tử hữu hạn để tính dầm sau đó mở rộng cho kết cấu khung hai chiều. Ta chỉ xét dầm có mặt cắt ngang đối xứng so với mặt phẳng tải trọng. Sơ đồ hoá dầm chịu uốn và biến dạng (độ võng) của trục dầm được minh hoạ như Hình 9.1 dưới đây. Trong trường hợp biến dạng nhỏ, ta đã có kết quả quen biết sau: y J M  (9.1) E    (9.2) Phương trình độ võng: J2 2 E M dx vd  (9.3) Mkm Pkp p y x Hình 9.1. (a) Sơ đồ hoá dầm chịu uốn; (b). Biến dạng của trục dầm x v (a) (b) kp km SinhVienKyThuat.Com 153 Trong đó:  là ứng suất pháp,  là biến dạng dài, M là mômen uốn nội lực trên mặt cắt ngang, v là độ võng của trục x và J là mômen quán tính của mặt cắt ngang đối với trục trung hoà (trục z đi qua trọng tâm mặt cắt ngang). 2. THẾ NĂNG Năng lượng biến dạng trong một phân tố có chiều dài dx được xác định bởi: dAdxdu  A2 1  (9.4) Thay (9.1) và (9.2) vào (9.4) ta được: dx EJ MdxdAy EJ MdAdxdu A 2 2 2 2 2 A 2 1 2 1 2 1          (9.5) Từ (9.5) ta tính được năng lượng biến dạng toàn phần trong dầm: dx dx vdEJU         L 2 2 2 1 (9.6) Thế năng  của dầm được xác định bởi:         km kmkm kp kpkp LL MvPvdxdx dx v  0 2 0 2 2 pdEJ 2 1 (9.7) Trong đó: p là tải trọng phân bố trên một đơn vị dài; Pkp là lực tập trung tại điểm kp; Mkm là mômen của ngẫu lực tại điểm km; vkp là độ võng tại điểm kp và km là góc xoay của mặt cắt ngang tại km. 3. HÀM DẠNG HERMITE Giả sử ta chia dầm thành bốn phần tử như Hình 9.2, mỗi phần tử có 2 nút; mỗi nút có 2 bậc tự do. Hai bậc tự do của nút i được ký hiệu là Q2i-1 và Q2i. Trong đó Q2i-1 là độ võng; Q2i là góc xoay. Q là véctơ chuyển vị chung: Q = [ Q1, Q2 , Q3, ..., Q10 ]T (9.8) Bốn bậc tự do địa phương của mỗi phần tử được ký hiệu bởi: q = [q1, q2, q3, q4]T. trong đó: q1, q3 là độ võng và q2, q4 là góc xoay. SinhVienKyThuat.Com 154 Các hàm dạng để nội suy chuyển vị trên một phần tử dầm sẽ được xác định theo ; trong đó -1   1. Các hàm dạng đối với phần tử dầm khác với các hàm dạng mà ta đã biết trong các chương trước. Vì có cả chuyển vị ngang (độ võng) và góc xoay do đó chúng ta đưa vào hàm dạng Hermite như sau: Hi = a i + b i  + c i2 + di 3; với (i = 1, 2, 3, 4) (9.9) Các hàm dạng Hermite thoả mãn các điều kiện: H1 H1’ H2 H2’ H3 H3’ H4 H4’  = -1 1 0 0 1 0 0 0 0  = 1 0 0 0 0 1 0 0 1 Các hệ số ai, bi, ci và di được xác định nhờ các điều kiện trong bảng trên. Kết quả cho ta: H1 = ( 1 - )2(2 + )/4 = (2 - 3 + 3)/4 H2 = (1 - )2 (1 + )/4 = (1 -  - 2 + 3)/4 H3 = (1 + )2(2 - )/4 = (2 + 3 - 3)/4 (9.10) H4 = (1 + )0( - 1)/4 = (-1 -  + 2 + 3)/4 Biểu diễn hình học các hàm dạng như Hình 9.3: 1 2 3 4 e Q10 Q9 Q6 Q8 Q5 Q6 Q3 Q4 Q1 Q2 q4 q3 q1 q2 2 1 v1 v’2 v2 v’1 (a) (b) Hình 9.2. Rời rạc dầm bằng phần tử hữu hạn SinhVienKyThuat.Com 155 4. MA TRẬN ĐỘ CỨNG CỦA PHẦN TỬ DẦM Trước hết, ta sử dụng các hàm dạng Hermite trên để xây dựng biểu thức chuyển vị v(): v() = H1v1 + H2 (dv/d)1 + H3v2 + H4 (dv/d)2 (9.11) Mặt khác:  222 1 2 1 1221 21 xxxxxxx  (9.12) Vì (x2 - x1 = le) là độ dài phần tử, do đó: d ldx e 2  ; và  d dx dx dv d dv  , nên ta có quan hệ đạo hàm sau: dx dvl d dv e 2   (9.13) Chú ý rằng (dv/dx) tại nút 1 và nút 2 chính là q2 và q4, suy ra:   44332211 22 qHlqHqHlqHv ee  (9.14) hay viết dưới dạng cô đọng: Hình 9.3. Hàm dạng Hermite SinhVienKyThuat.Com 156 v = H q (9.15) Trong đó:      4321 22 HlHHlHH ee (9.16) Năng lượng biến dạng của phần tử được xác định bởi: dx dx vdEJU e e 2 2 2 2 1         Từ (9.13) suy ra: d dv ldx dv e 2  và 2 2 22 2 4 d vd ldx vd e  Thay v = Hq, ta được:                               22 31 2 3 22 31 2 3 16 2