Giáo trình Bảo toàn năng lượng - Chương 1: Các phương pháp năng lượng - Biến dạng phẳng và ứng suất phẳng

− (b) Ứng suất phẳng xác định theo quan hệ: σ = D ε Matrận cứng phần tử tính theo côngthức: [k] = = [B]dVBCB V T ]][[][∫ T[D] [B] t = tA[B]∫ A dA T[D] [B] (c) Matrận [k] = [kn] + [ks], bậc 6x6 có dạng: 154 ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ −−− −− −− − − −= 2 2121212131213121322132 2 212131213121322132 2 31313131323132 2 3131323132 2 323232 2 32 2 )1(4 ][ xxyyxxyxxxy yyxyyyxyy xxyxxyy yyxyy DXxxy y A Etkn ννν νν νν ν ν υ (d) [ ] ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ −−− −− −− − − += 2 2121212131213121322132 2 212131213131322132 2 31313131323132 2 3131323132 2 323232 2 32 )1(8 yyxyyyxyyyx xxyyxxyxx yyxyyyx xyyxx yyx x A Etks υ (e) Nếu ký hiệu: 3 2 1 3 2 1 33 22 1 1 11 1 111 x x A c y y A b yx yx A a =−== , các hệ số mang chỉ số j = 2 và 3 tính theo qui luật tương tự, có thể viết công thức cho chuyển vị dạng sau: ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎦ ⎤ ++++ ++ ⎢⎣ ⎡ ++ ++++=⎭⎬ ⎫ ⎩⎨ ⎧ 3 2 1 333222 333 111 222111 0 00 00 01 ),( ),( U U U ycxbaycxba ycxba ycxba ycxbaycxba Ayx yxu v (f) Vector ứng suất mang dạng: ( ) ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ −−−−−−− −−− −−− ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ −−−−−−−− −−− −−− −=⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ 3 2 1 121213 122131 132131 132323 132332 132332 2 )( 2 1)( 2 1)( 2 1 )()()( )()()( )( 2 1)( 2 1)( 2 1 )()()( )()()( 1 U U U yyxxyy xxyyxx yyyyxx xxyyxx yyxxyy yyxxyy A E xy y x υυυ υ υυ υυυ υυ υ υτ σ σ (g) Những phần tử tam giác dùng trong bài toán phẳng đa dạng. Một số trong nhóm có thể thấy như sau đây. 155 Hình Hàm plante giúp xác định ma trận [k], ký hiệu Ke và vector lực fe phần tử bốn cạnh trong bài toán phẳng. Qui ước rằng nếu người dùng ghi ptype = 1 chúng ta sẽ xem xét bài toán trạng thái ứng suất phẳng, còn ptype = 2 sẽ bàn đến trạng thái biến dạng phẳng. Giá trị của ptype đăng ký tại vector ep = [ptype t ]. Toạ độ xi, i=1,2,3 , yi, i=1,2,3 ghi tại ex[]và ey[]. function [Ke,fe]=plante(ex,ey,ep,D,eq) % Ke=plante(ex,ey,ep,D) % [Ke,fe]=plante(ex,ey,ep,D,eq) ptype=ep(1); t=ep(2); bx=0.; by=0.; if nargin==5; bx=eq(1); by=eq(2); end C=[ 1 ex(1) ey(1) 0 0 0 0 0 0 1 ex(1) ey(1) 1 ex(2) ey(2) 0 0 0 0 0 0 1 ex(2) ey(2) 1 ex(3) ey(3) 0 0 0 0 0 0 1 ex(3) ey(3)]; A=1/2*det([ones(3,1) ex' ey']); %--------- plane stress -------------------------------------- if ptype==1 B=[0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0]*inv(C); colD=size(D,2); if colD>3 Cm=inv(D); Dm=inv(Cm([1 2 4],[1 2 4])); else Dm=D; end Ke=B'*Dm*B*A*t; fe=A/3*[bx by bx by bx by]'*t; %--------- plane strain -------------------------------------- 156 elseif ptype==2 B=[0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0]*inv(C); colD=size(D,2); if colD>3 Dm=D([1 2 4],[1 2 4]); else Dm=D; end Ke=B'*Dm*B*A*t; fe=A/3*[bx by bx by bx by]'*t; else error('Error ! Check first argument, ptype=1 or 2 allowed') return end Tính ứng suất và biến dạng trong phần tử tam giác nhờ hàm plants(ex,ey,ep,D,ed). function [es,et]=plants(ex,ey,ep,D,ed) % [es,et]=plants(ex,ey,ep,D,ed) ptype=ep(1); rowed=size(ed,1); rowex=size(ex,1); %--------- plane stress -------------------------------------- if ptype==1 colD=size(D,2); if colD>3 Cm=inv(D); Dm=inv(Cm([1 2 4],[1 2 4])); else Dm=D; end if rowex==1 incie=0; else incie=1; end et=[];es=[];ie=1; for i=1:rowed C=[ 1 ex(ie,1) ey(ie,1) 0 0 0 0 0 0 1 ex(ie,1) ey(ie,1) 1 ex(ie,2) ey(ie,2) 0 0 0 0 0 0 1 ex(ie,2) ey(ie,2) 1 ex(ie,3) ey(ie,3) 0 0 0 0 0 0 1 ex(ie,3) ey(ie,3)]; B=[0 1 0 0 0 0; 0 0 0 0 0 1; 0 0 1 0 1 0]*inv(C); 157 ee=B*ed(i,:)'; if colD>3 ss=zeros(colD,1); ss([1 2 4])=Dm*ee; ee=Cm*ss; else ss=Dm*ee; end et=[et;ee']; es=[es;ss']; ie=ie+incie; end %--------- plane strain -------------------------------------- elseif ptype==2 colD=size(D,2); if rowex==1 incie=0; else incie=1; end et=[];es=[];ie=1;ee=zeros(colD,1); for i=1:rowed C=[ 1 ex(ie,1) ey(ie,1) 0 0 0 0 0 0 1 ex(ie,1) ey(ie,1) 1 ex(ie,2) ey(ie,2) 0 0 0 0 0 0 1 ex(ie,2) ey(ie,2) 1 ex(ie,3) ey(ie,3) 0 0 0 0 0 0 1 ex(ie,3) ey(ie,3)]; B=[0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0]*inv(C); e=B*ed(i,:)'; if colD>3 ee([1 2 4])=e; else ee=e; end et=[et;ee']; es=[es;(D*ee)']; ie=ie+incie; end else error('Error ! Check first argument, ptype=1 or 2 allowed') return end Ví dụ 1: Sử dụng công thức có sẵn của phần tử tam giác làm việc ở trạng thái ứng suất phẳng tính ứng suất và biến dạng tấm hình vuông cạnh l. Chiều dày tấm t = const. Mô đun đàn hồi vật liệu E, hệ số Poisson ν = 1/6. Nhờ tính đối xứng cấu hình và đối xứng tải, đưa một nửa tấm vào tính toán. Phần dùng tính toán của tấm chia làm 4 phần tử tam giác, đánh số (1), (2), (3), (4) 158 như giới thiệu tại hình. Số nút dùng trong tính toán 6. Căn cứ cách đánh số này xác định tọa độ các nút cho mỗi phần tử. 1 2 3 4 5 6 7 8 9 (3) (4) (2)(1) L L L/4 L/ 4 x y Hình 8 Phần tử (1). x1 = 0 ; y1 = 0; x2 = 0; y2 = - ½ l x3 = ½ l; y3 = - ½l Phần tử (2). x1 = 0 ; y1 = 0; x3 = ½ l ; y3 = - ½l x4 = ½l ; y4 = 0 Phần tử (3). x5 = 0 ; y5 = 0; x1 = 0; y1 = - ½ l x4 = ½ l; y4 = - ½ l Phần tử (4). x5 = 0 ; y5 = 0; x4 = ½ l; y4 = - ½ l x6 = ½ l; y6 = 0 Vector tải chỉ có mặt trong thành phần {p} của phần tử (1), mang giá trị P/2. Chuyển vị ngang nút số 1, số 2, số 5 tại trục đối xứng bị hạn chế, nút 6 bị hạn chế theo chiều đứng, v6 = 0. Sau khi tính các hệ số cứng phần tử tiến hành tổng hợp ma trận cứng toàn hệ thống và xử lý điều kiện biên, sẽ nhận được hệ phương trình sau. 159 ⎪⎪ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ = ⎪⎪ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − − − −− −− −− −−−− 0 0 0 0 0 0 0 728,0 214,0728,0 086,00456,1 214,03,03,0456,1 00514,0086,0728,0 00214,0214,00728,0 0000214,0086,0728,0 0514,0428,03,003,0514,0456,1 2 1 6 5 4 4 3 3 2 1 P u u u DX Et v v v v v Kết quả giải hệ phương trình đại số có dạng: ⎪⎪ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ − = ⎪⎪ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ 3909,0 1090,1 7772,0 0066,0 1090,1 0692,0 9033,1 2728,1 6 5 4 4 3 3 2 1 Et P u u u v v v v v Ứng suất tấm tính tại các nút thực hiện theo (g). Ví dụ tính cho phần tử (4) chúng ta nhận được ứng suất sau. ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ − −×= ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ × ⎥⎥ ⎥ ⎦ ⎤ − − ⎢⎢ ⎢ ⎣ ⎡ −− − = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ 536,0 465,1 583,0 0 3909,0 7772,0 0066,0 0109,1 0 208,0208,00 2/2/2/ 2/2/2/ 208,0208,00 00.2/ 002/ 114,4 2 tl P ll lll lll ll l l tl E xy y x υ υυ υ τ σ σ Một trong những ứng dụng hữu ích phần tử trạng thái ứng suất phẳng là nghiên vứu tập trung ứng suất. Ví dụ tiếp theo giới thiệu kết quả tính ứng suất gần lỗ khoét hình tròn của cơ cấu chịu ứng suất nén đơn vị, hình 9. 160 Hình 9 Phần tử ba cạnh, trong vật thể tròn xoay Chuyển vị trong phần tử được trình bày trong hệ tọa độ Ozr, trục 0z hướng lên trên, trục 0r hướng tâm, bắt đầu từ tâm tại trục tròn xoay. u = = [N]{δ} (a) ⎭⎬ ⎫ ⎩⎨ ⎧ ),( ),( zrw zru [N] = (b) ⎥⎦ ⎤⎢⎣ ⎡ kji kji NNN NNN 000 000 {δ} = , còn {U U U U i j k ⎧ ⎨⎪ ⎩⎪ ⎫ ⎬⎪ ⎭⎪ i} = , tương tự vậy chúng ta viết U⎭⎬ ⎫ ⎩⎨ ⎧ i i w u j, Uk. Hàm = N N N j k ⎧ ⎨⎪ ⎩⎪ ⎫ ⎬⎪ ⎭⎪ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ++ ++ ++ zcrba zcrbaj zcrba A kkk jj iii 2 1 , với A - diện tích tam giác. (c) A= ( )jkijkiikkjji zrzrzrzrzrzr +++++21 161 Hình 10 Quan hệ biến dạng-chuyển vị trong phần tử thể hiện qua biểu thức: {ε} = = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ rz zz rr ε ε ε ε θθ ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ + r w z u z w r u r u ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = [B]{δ} (d) trong đó: [B] = ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ kkjjii kji kji kji bcbcbc ccc rNrNrN bbb A 000 0)/(0)/(0)/( 000 2 1 (e) Ma trận cứng tính theo thứ tự: [k] = [B]T[D][B] dV, trong đó dV = 2πr Ve ∫ Ve ∫ m A (f) [k] = [B]T[D][B]. 2πrm A, (g) với rm = ( ri + rj + rk)/3; (h) [D] = E( ) ( )( ) 1 1 1 2 − + − ν ν ν 1 1 0 1 1 1 0 1 0 0 1 2 2 1 0 1 1 0 1 ν ν ν ν ν ν ν ν ν ν ν ν ν ν / ( ) / ( ) / ( ) / ( ) ( ) / [ ( )] / ( ) / ( ) − − − − − − − − ⎡ ⎣ ⎢⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥⎥ Vector tải trên biên xác định từ biểu thức: 162 { } rds f f L L L L L L dSfNp z r S k k j j i i S s T s ij π2][ 1 1 ⎭⎬ ⎫ ⎩⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ == ∫∫ (i) Vector lực do nhiệt độ thay đổi, gây ra tính theo cách chúng ta vẫn dùng cho các phần tử trong nhóm phần tử phẳng. [ ] ( ) [ ] ArB TETEDBp III V T t πυ αα 2 0 1 1 1 21 0 1 1 1 ][}{ ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ −= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = ∫ (j) Biểu thức xác định lực khối được viết dạng sau: { } rdA f f L L L L L L dVfNp z r A k k j j i i V b T b π2][ ⎭⎬ ⎫ ⎩⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ == ∫∫ (k) trong đó r = ri Li + rj Lj + rk Lk Phần tử bốn cạnh vật tròn xoay có dạng sau: Hình 11 163 Tấm mỏng Tấm mỏng dùng trong tài liệu theo định nghĩa Timoshenko, thỏa mãn các giả thiết, trong đó có giả thuyết của Kirchhoff như đã ghi trong lý thuyết tấm. Hình 12 Dưới tác động của lực và moment ( tập trung hoặc phân bố) tác động theo phương pháp tuyến tấm bị uốn. Ứng suất trong tấm được ký hiệu như trên hình 11. Hình 13. 164 Hình 14. Ký hiệu lực, momen của tấm bị uốn. Ứng suất này gồm ứng suất cắt σyz, σxz, σxy và ứng suất tại các mặt cắt ngang σxx, σyy. Moment và lực được định nghĩa như sau: Mx = ; Mzdz t t x∫+− 2/2/ σ y = ; Mzdztt y∫+− 2/2/ σ xy = zdztt xy∫+− 2/2/ σ qx = σ− +∫ tt//22 xz dz; qy = σ−+∫ tt//22 yz dz Quan hệ giữa chúng thể hiện trong hệ phương trình cân bằng lực. 0=++ p y q x q yx ∂ ∂ ∂ ∂ x xyx q y M x M =+ ∂ ∂ ∂ ∂ y yxy q y M x M =+ ∂ ∂ ∂ ∂ trong đó p - lực phân bố, tác động theo phương pháp tuyến với mặt tấm. Biến dạng tấm thể hiện bằng: {ε} = = -z ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ x y x ε ε ε ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ yx w y w x w ∂∂ ∂ ∂ ∂ ∂ ∂ 2 2 2 2 2 2 (a) 165 Nếu biểu diễn chuyển vị tấm theo cách quen thuộc: w = [N]{δ} (b) còn {δi } = = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ yi xy iw θ θ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ⎟⎠ ⎞⎜⎝ ⎛ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛− i i i x w y w w ∂ ∂ ∂ ∂ (c) Vector biến dạng: {ε} = = [B]{δ} (d) ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ x y x ε ε ε trong đó: [Bi] = ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ − − yx Ni y Ni x Ni ∂∂ ∂ ∂ ∂ ∂ ∂ ][2 ][ ][ 2 2 2 2 2 (e) Vector ứng suất tấm được xét trong bài toán này: {σ} ≡ {M} = [D]( {ε} - {ε0 } ) (f) Matrận [D] được ghi tại phần đầu tài liệu. Với vật liệu đẳng hướng, [D] trong bài toán phẳng có dạng sau. [D] = ( )2 3 112 υ− Et ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − 2/)1(00 01 01 υ υ υ Phương trình chuyển vị của tấm dưới tác động lực phân bố p(x,y) có dạng: D pw yyxx =⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ++ 4 4 22 4 4 4 2 ∂ ∂ ∂∂ ∂ ∂ ∂ D = )1(12 2 3 υ− Et Ứng suất trong phần tử tấm, tính tại các mặt cắt qua tâm phần tử có dạng: σx = Nx / t + 12Mx.z /t3 σy = Ny /t + 12My.z / t3 τxy = Nxy / t ; τyx = Nyx / t 166 Ứng suất chính tính bằng công thức: σ1 = 12 (σx + σy) + 1 2 22 4)( xyxy τσσ ++ σ2 = 12 (σx + σy) - 1 2 22 4)( xyxy τσσ +− Phần tử tấm PLATE có khả năng tiếp nhận tải trọng tập trung, tải trọng phân bố đều trên bề mặt phần tử. Với tấm hình tam giác sẽ có 3 nút đặt tại ba đỉnh, số bậc tự do 18. Tấm phẳng 4 nút có 24 bậc tự do. Phiếm hàm áp dụng cho tấm bị uốn có dạng chung: Π = 1 2 V ∫ {ε}T {σ}dV - p(x,y)wdV (g) V ∫ Từ đó matrận cứng phần tử tấm được tính theo công thức quen thuộc: [k] = [B] A ∫∫ T [D] [B] dxdy (h) Vector lực tính theo: {F} = - [N] A ∫∫ T p(x,y)dxdy (i) Tấm mỏng hình chữ nhật, cạnh axb, 12 bậc tự do. Phương trình chuyển vị: w(x,y) = [P]{a} = [ 1 x y x2 y2 xy x2y xy2 x3 y3 x3y xy3 ] {a} (j) Phương trình trên đây thỏa mãn điều kiện liên tục của bài toán uốn tấm: ∇2 ∇2 w = 0 (k) Aùp dụng phương trình trên cho chuyển vị các nút sẽ nhận được: {δ} = [C]{a} (l) Matrận [C] có kích thước 12x12, phụ thuộc vào toạ độ các nút i, i= 1, 2, 3, 4. {a} = [C]-1 {δ} (m) Từ đó: w = [N] {δ} = [P] [C]-1 {δ} (o) [N] = [P] [C]-1 (p) Biểu thức Ni, i=1, 2, 3, 4 sau khai triển: [Ni] = 1 2 [ ( ξ0 + 1)(η0 + 1) (2 + ξ0 + η0 - ξ2 - η2), 167 aξi(ξ0 +1)2(ξ0 -1)( η0 +1), bηi(ξ0 +1)( η0 +1)2(η0 -1) ], (q) trong đó: ξ = x x a c− η = y y b c− ξ0 = ξξi η0 = ηηi Thay các giá trị trên vào biểu thức tính {ε} sẽ nhận được: {ε} = [Q]{a} = [Q]{C]-1 {δ} và từ đó: [B] = [Q]{C]-1 (r) [Q] = 0 0 0 2 0 0 6 2 0 0 6 0 0 0 0 0 0 2 0 0 2 6 0 6 0 0 0 0 2 0 0 4 4 0 6 62 2 − − − − − − − ⎡ ⎣ ⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥ x y xy x y xy x y x y − Matrận cứng phần tử tính theo (h). Nếu đưa biểu thức ( r) vào vị trí [B] tại (h) sẽ nhận được biểu thức tính [k] như sau: [k] = {[C]-1 }T ( ∫∫ [Q]T [D] [Q] dxdy ) [C]-1. Bằng cách tương tự biểu thức tính vector lực giờ có thể mang dạng: {pe} = { -[C]-1 }T ∫∫ [P]T pdxdy. (s) Hàm platre.m phục vụ tính ma trận cứng [k], vector tải fe phần tử tấm chữ nhật như sau: Dữ liệu đầu vào gồm: ex= [x1 x2 x3 x4 ] ey= [y1 y2 y3 y4 ] ep = [t] Hình 14 2a 2b t eq = [qz] matrận C hiểu như sau đây: 168 ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ −−−−−− −−−− −−−−−− −−−−−−−− −−−− −− −−−−−− −−−− −−−− −−−−−− = 3223 3322 33322322 3223 3322 33322322 3223 3322 33322322 3222 3322 33322322 302302010 332020100 1 302302010 332020100 1 302302010 332020100 1 302302010 332020100 1 bbabababa ababababa abbababbaabababa bbabababa ababababa abbababbaabababa bbabababa ababababa abbababbaabababa bbabababa ababababa abbababbaabababa C function [Ke,fe]=platre(ex,ey,ep,D,eq) % Ke=platre(ex,ey,ep,D) % [Ke,fe]=platre(ex,ey,ep,D,eq) Lx=ex(3)-ex(1); Ly=ey(3)-ey(1); t=ep(1,1); % D=t^3/12*D; % A1=Ly/(Lx^3); A2=Lx/(Ly^3); A3=1/Lx/Ly; A4=Ly/(Lx^2); A5=Lx/(Ly^2); A6=1/Lx; A7=1/Ly; A8=Ly/Lx; A9=Lx/Ly; % C1= 4*A1*D(1,1)+4*A2*D(2,2)+2*A3*D(1,2)+5.6*A3*D(3,3); C2=-4*A1*D(1,1)+2*A2*D(2,2)-2*A3*D(1,2)-5.6*A3*D(3,3); C3= 2*A1*D(1,1)-4*A2*D(2,2)-2*A3*D(1,2)-5.6*A3*D(3,3); C4=-2*A1*D(1,1)-2*A2*D(2,2)+2*A3*D(1,2)+5.6*A3*D(3,3); C5=2*A5*D(2,2)+A6*D(1,2)+0.4*A6*D(3,3); C6=2*A4*D(1,1)+A7*D(1,2)+0.4*A7*D(3,3); C7=2*A5*D(2,2)+0.4*A6*D(3,3); C8=2*A4*D(1,1)+0.4*A7*D(3,3); C9= A5*D(2,2)-A6*D(1,2)-0.4*A6*D(3,3); C10=A4*D(1,1)-A7*D(1,2)-0.4*A7*D(3,3); C11=A5*D(2,2)-0.4*A6*D(3,3); C12=A4*D(1,1)-0.4*A7*D(3,3); C13=4/3*A9*D(2,2)+8/15*A8*D(3,3); C14=4/3*A8*D(1,1)+8/15*A9*D(3,3); C15=2/3*A9*D(2,2)-8/15*A8*D(3,3); C16=2/3*A8*D(1,1)-8/15*A9*D(3,3); C17=2/3*A9*D(2,2)-2/15*A8*D(3,3); C18=2/3*A8*D(1,1)-2/15*A9*D(3,3); C19=1/3*A9*D(2,2)+2/15*A8*D(3,3); C20=1/3*A8*D(1,1)+2/15*A9*D(3,3); C21=D(1,2); % Keq=zeros(12,12); 169 Keq(1,1:12)=[C1 C5 -C6 C2 C9 -C8 C4 C11 -C12 C3 C7 -C10]; Keq(2,2:12)=[C13 -C21 C9 C15 0 -C11 C19 0 -C7 C17 0]; Keq(3,3:12)=[C14 C8 0 C18 C12 0 C20 -C10 0 C16]; Keq(4,4:12)=[C1 C5 C6 C3 C7 C10 C4 C11 C12]; Keq(5,5:12)=[C13 C21 -C7 C17 0 -C11 C19 0]; Keq(6,6:12)=[C14 C10 0 C16 -C12 0 C20]; Keq(7,7:12)=[C1 -C5 C6 C2 -C9 C8]; Keq(8,8:12)=[C13 -C21 -C9 C15 0]; Keq(9,9:12)=[C14 -C8 0 C18]; Keq(10,10:12)=[C1 -C5 -C6]; Keq(11,11:12)=[C13 C21]; Keq(12,12)=[C14]; Keq=Keq'+Keq-diag(diag(Keq)); % if nargin==5 R1=eq*Lx*Ly/4; R2=eq*Lx*Ly^2/24; R3=eq*Ly*Lx^2/24; % feq(:,1)=[R1 R2 -R3 R1 R2 R3 R1 -R2 R3 R1 -R2 -R3]'; end Ke=Keq; fe=feq; Tính ứng suất và biến dạng điểm trong tấm chịu uốn thực hiện hàm platrs(ex,ey,ep,D,ed). Vector ep = [t], ghi chiều dày tấm. Vector ed chứa chuyển vị nút tính theo hàm extract: ed = [u1 u2 . . . u12]; Kết quả tính ghi trong các ma trận sau: es = [M V ]; et = κT. function [es,et]=platrs(ex,ey,ep,D,ed) % [es,et]=platrs(ex,ey,ep,D,ed) Lx=ex(3)-ex(1); Ly=ey(3)-ey(1);t=ep(1,1); % D=t^3/12*D; % A1=D(2,2)/2/Ly; A2=D(1,1)/2/Lx; A3=D(1,2)/2/Ly; A4=D(1,2)/2/Lx; A5=D(3,3)/2/Ly; A6=D(3,3)/2/Lx; A7=4*D(3,3)/Lx/Ly; B1=6*D(2,2)/Ly/Ly/Ly; B2=6*D(1,1)/Lx/Lx/Lx; B3=-3*D(2,2)/Ly/Ly; B4=3*D(1,1)/Lx/Lx; B5=D(1,2)/Lx/Ly; % mx=A3*(-ed(:,2)-ed(:,5)+ed(:,8)+ed(:,11))+A2*(ed(:,3)-ed(:,6)- ed(:,9)+ed(:,12)); 170 my=A1*(-ed(:,2)-ed(:,5)+ed(:,8)+ed(:,11))+A4*(ed(:,3)-ed(:,6)- ed(:,9)+ed(:,12)); mxy=A6*(ed(:,2)-ed(:,5)-ed(:,8)+ed(:,11))+A5*(-ed(:,3)- ed(:,6)+ed(:,9)+ed(:,12))... +A7*(ed(:,1)-ed(:,4)+ed(:,7)-ed(:,10)); m1=0.5.*(mx+my)+sqrt(0.25.*(mx-my).^2+mxy.^2); m2=0.5.*(mx+my)-sqrt(0.25.*(mx-my).^2+mxy.^2); alfa=0.5*180/pi*atan2(mxy,(mx-my)/2); vx=B5*(-ed(:,2)+ed(:,5)- ed(:,8)+ed(:,11))+B4*(ed(:,3)+ed(:,6)+ed(:,9)+ed(:,12))... +B2*(-ed(:,1)+ed(:,4)+ed(:,7)-ed(:,10)); vy=B3*(ed(:,2)+ed(:,5)+ed(:,8)+ed(:,11))+B5*(ed(:,3)-ed(:,6)+ed(:,9)- ed(:,12))... +B1*(-ed(:,1)-ed(:,4)+ed(:,7)+ed(:,10)); es=[mx my mxy vx vy]; et=-inv(D)*[mx;my;mxy]; et=et'; Hàm extract như sau đây, dùng chung cho tòn bộ chương trình tính. function [ed]=extract(edof,a) % ed=extract(edof,a) [nie,n]=size(edof); % t=edof(:,2:n); % for i = 1:nie ed(i,1:(n-1))=a(t(i,:))'; end Tấm tam giác ba nút, 9 bậc tự do. Hình 15 Phương trình chuyển vị: w(x,y) = [P]{a}= [ 1 x y x2 y2 x3 (x2y + xy2 ) y3 ] {a} (a) Phương trình trên đây thỏa mãn điều kiện liên tục của bài toán uốn tấm: ∇2 ∇2 w = 0 (b) Áp dụng phương trình trên cho chuyển vị các nút sẽ nhận được: {δ} = [C]{a} (c) Matrận [C] có kích thước 9x9, phụ thuộc vào toạ độ các nút i, i= 1, 2, 3. 171 [C] = 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 2 0 0 3 0 1 0 0 0 0 0 1 0 0 1 0 2 0 2 3 0 1 0 2 0 3 2 0 2 2 2 2 3 2 2 2 2 2 2 3 3 3 2 3 3 3 2 3 3 3 2 3 3 3 2 3 3 3 3 3 3 3 2 3 2 3 3 3 2 3 2 3 3 − − − − + + − − − − − + ⎡ ⎣ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥⎥⎥y y y y y y x y x x y y x x y x y y x y x y x y x y x y x y ( ) ( ) ( ) ⎥⎥⎥⎥⎥⎥⎥ y {a} = [C]-1 {δ} (d) Từ đó: w = [N] {δ} = [P] [C]-1 {δ} (e) [N] = [P] [C]-1 (f) {ε} = [B] {δ} = [β]{a} = [β] [C]-1 {δ} từ đó [B] = [β] [C]-1. trong đó [β] = -z 0 0 0 2 0 0 6 2 0 0 0 0 0 0 2 0 2 6 0 0 0 0 2 0 0 4 0 x y x y x y( )+ ⎡ ⎣ ⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥ Matrận cứng của phần tử tính theo (h) phần trên:[k] = t [B] A ∫∫ T [D] [B] dxdy hoặc [k] = {[C]-1 }T ( ∫∫ [β]T [D] [β] dxdy ) [C]-1. Tích phân nằm trong ngoặc đơn có dạng:[k’] = ∫∫ [β]T [D] [β] dxdy = ∫∫− A dxdy Et )1(12( 2 3 ν ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ + −− +−+++−+ − 2 2 2 36)(123612012000 ])1(8 ))(812[( )(12)(4))(1(4)(4000 3612012000 404000 )1(20000 4000 000 00 0 yyxxyyy xy yx yxxyxyxyx xxx DX υυυ υ υυυυυ υ υ υ (g) Tích phân , trong hệ tọa độ chung hoặc hệ tọa độ riêng đều thực hiện theo qui tắc chung nêu sau đây. ∫ A dxdy 232 1 yxAdxdy A ==∫ 172 2 2 36 1 yxAXxdxdy C A ==∫ ( )32236 1 yyyxAYydxdy C A +==∫ ( ) ( ) ( )[ ] 23322222 12112 yxXXXXXXAAXdxdyx CkCjCiCA =−+−+−+=∫ ( )( ) ( )( ) ( )([ ]) ( )32223 224 1 12 yyyx YYXXYYXXYYXXAAYXxydxdy CkCkCjCjCiCiCC A += −−+−−+−−+=∫ ( ) ( ) ( )[ ] ( )2332222322222 12112 yyyyyxYYYYYYAAYdxdyy CkCjCiCA ++=−+−+−+=∫ Trong các công thức XC và YC tính theo cách sau: ( )kjiC XXXX ++= 31 ( )kjiC YYYY ++= 31 Trường hợp tấm không nằm trùng mặt 0xy của hệ tọa độ chung, cần tính chuyển các đặc trưng nhờ ma trận chuyển, như đã đề cập phần trên. Ma trận chuyển [λ] có dạng: [Λ] = , trong đó: [T] = T T T 0 0 0 0 0 ⎡ ⎣ ⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥ 0 1 0 0 0 0 0 0 0 l m l y m x x oy ⎡ ⎣ ⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥ Matrận [Λ], 9x9 có dạng: [Λ] = (h) 9 8 7 6 5 4 3 2 1 0000000 0000000 001000000 0000000 0000000 000001000 0000000 0000000 000000001 00 00 00 00 00 00 ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ yy xx yy xx yy xx ml ml ml ml ml ml 173 Matrận [k] và vector lực được tính chuyển sang hệ tọa độ chung theo công thức: [k]chung = [Λ]T [k] [Λ], {p}chung = [Λ]T{p}. Vector tải do ảnh hưởng nhiệt xem xét theo các cách tương tự. Nếu phân bố nhiệt trong tấm ghi bằng công thức: 0),,( Tt zTzyxT += (i) Trong đó T và T0 là hằng số chỉ nhiệt độ trung bình và chênh lệch nhiệt độ giữa mặt trên z = t/2 và mặt dưới z = -t/2 của tấm. Vector tải do ảnh hưởng nhiệt trong trường hợp này mang dạng: [ ] [ ]{ }dVDBp V T∫= 0}{ ε (j) Công thức này còn được hiểu trong khuôn khổ phần tử tam giác: { } [ ]( ) [ ] [ ] [ ]( ) ∫ ∫ ∫ ∫ −= − −= − ⎪⎪ ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎪⎪ ⎨ ⎧ + ⎟⎠ ⎞⎜⎝ ⎛ += = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎟⎠ ⎞⎜⎝ ⎛ += 2/ 2/ 3 0 1 2/ 2/ 0 1 6 )(2 6 0 1 0 0 0 . 0 1 1 t tz A T t tz A TT dA y yx x x dzT t zTzC dAT t zTDBdzCp α α Sau tích phân công thức cuối được viết lại dạng sau đây: { } [ ]( ) ( ) ( ) ⎪⎪ ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎪⎪ ⎨ ⎧ + ++ ×−−= − 32 323 3 1 2 0 3/ 1 0 1 0 0 0 )1(6 yy yyx x CAtTEp T υ α (k) 174 Vec to này khi chuyển sang hệ tọa độ chung sẽ là: [ ]{ }pT (l) Những phần tử tấm tam giác uốn, thường gặp trong tính toán: Hình 16 Ví dụ 1: Sử dụng phần tử tấm PLATE được trình bày dưới đây cho trường hợp tấm bị ngàm chiều dài cạnh l, chịu tác động áp lực phân bố đềàu p, theo phương pháp tuyến. Nhờ tính đối xứng của tấm qua hai trục và tính đối xứng của bản thân tải trọng chỉ sử dụng 1/4 tấm trong tính toán. Xét hai trường hợp điều kiện biên khác nhau như sau đây. a) Tấm ngàm cả bốn cạnh Nếu chia tấm vuông cạnh L làm bốn phần tử, độ cứng mỗi phần tử tính theo (s) dạng ma trận cỡ 12x12 và vec to tải cỡ 12x1 có dạng như nêu tại (t), giành cho phần tử PLATE. Nhờ tính đối xứng, như vừa nêu, chỉ sử dụng 1 phần tử trong tính toán tiếp theo. Tại biên qua trục đối xứng x = L/2 chuyển vị uj = 0, θx = 0, j số thứ tự nút trên đường đối xứng. Tại biên qua trục đối xứng y = L/2 chuyển vị vj = 0, θy = 0. 175 Nút chính giữa tấm chỉ tham gia chuyển động theo hướng Oz, w ≠ 0, các thành phần khác đều bằng 0. Điều kiện biên bài toán đặt ra rằng chuyển vị tất cả các nút tại ngàm trượt tiêu. Và như vậy sau khi xử lý điều kiện biên phương trình cân bằng chỉ còn chứa một ẩn là chuyển vị theo hướng Oz của nút tại giữa tấm. k.w = Q Kết quả giải phương trình cho phép xác định w = 1,14.10-3 4qL D và tiếp đó xác định các thành phần lực theo w. b) Tấm ngàm tại hai cạnh đối diện, hai cạnh còn lại tựa trên gối. Có thể tăng độ chính xác phép tính bằng biện pháp tăng số phần tử trong quá trình mô hình hóa tấm. Phương án tiếp theo, chia ¼ tấm vừa đề cập làm 4 phần tử bằng biện pháp chia đôi hai cạnh của hình. Phân bố nút, phần tử được trình bày tại hình. Điều quan tâm bây giờ là xác định chuyển vị nút 1, momen uốn Mx = My tính tại nút này. 1 2 3 4 5 6 7 8 9 (3) (4) (2)(1) L L L/4 L/ 4 x y Hình 17 Xác lập ma trận cứng tiến hành theo các bước: mỗi phần tử có ma trận cứng xác lập trong hệ tọa độ cục bộ, cỡ 12x12. Ví dụ, phần tử (1) sẽ có ma trận dạng sau: 176 [ ] 6 5 2 1 4 4 4 3 3 3 2 2 2 1 1 1 1 12,12 1 11,12 1 1,12 1 12,11 1 1,11