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 - Những phần tử thường dùng

g hợp các nút của phần tử còn thêm một bậc tự do, ví dụ còn dịch chuyển dọc trục, ma trận cứng phần tử phải cộng thêm các thành phần của phần tử BAR nêu tại dòng đầu phần đang đề cập. Đây sẽ là phần tử khung phẳng. 101 [ ] ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − − − − = 2 2 22 22 3 4 612 00 2604 6120612 0000 L LDX J AL LLL LL J AL J AL L EJk Ma trận chuyển áp dụng cho trường hợp này có dạng: [ ] ; 0 0 ⎥⎦ ⎤⎢⎣ ⎡= C C T với [ ] trong đó φ - góc giữa hướng trục Ox hệ tọa độ chính với O’x’ của hệ tọa độ cục bộ, gắn liền với phần tử. ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − = 100 0cossin 0sincos φφ φφ C Vector lực tác động lên phần tử Trường hợp tổng quát có thể xác định công ngọai lực tác động lên phần tử dưới dạng: W = ∫ p(x).u(x) dx (a) Thay u = [N] {δ } vào (a) có thể viết: W = ∫ p[N]T dx {δ }, tích phân tiến hành dọc chiều dài L dầm. (b) Biểu thức lực tổng quát có dạng: {F} = p {N]∫L 0 T dx Trường hợp p(x) = const vector p có dạng sau, nếu không tính thành phần dọc trục: { } ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ − = 12 2 12 2 2 2 pL pL pL pL F (c) Trường hợp phân bố tải trọng hình thang p(x) = L ab )( − x + a, trong đó a, b - giá trị tải trọng phân bố tại nút 1 và nút 2 của dầm. 102 { } ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ +− + + + = )20/30/( 35,015,0 )30/20/( 15,035,0 baL ba baL ba F Trường hợp dầm chịu tác động lực tập trung P, đặt tại nút cách nút 1 đoạn ξ: {F} = ∫ P.δ( ξ - x) [N]T dx, trong đó δ(ξ - x) - hàm Dirac. Dưới đây chúng ta cùng xem xét những hàm viết bằng ngôn ngữ Matlab phục vụ tính các thành phần ma trận cứng phần tử ma trận chuyển {T}, vector tải fe. Chương trình viết bằng Matlab sau đây trích từ tài liệu dạy phương pháp PTHH khoa Cơ học kết cấu và Cơ học vật rắn trường đại học Lund soạn. Chương trình tính [k] trong hệ tọa độ riêng, ký hiệu Ke, và ma trận chuyển ký hiệu G với tên gọi Beam2e: function [Ke,fe]=beam2e(ex,ey,ep,eq); % Ke=beam2e(ex,ey,ep) % [Ke,fe]=beam2e(ex,ey,ep,eq) ---------------------- b=[ ex(2)-ex(1); ey(2)-ey(1) ]; L=sqrt(b'*b); n=b/L; E=ep(1); A=ep(2); I=ep(3); qx=0; qy=0; if nargin>3; qx=eq(1); qy=eq(2); end Kle=[E*A/L 0 0 -E*A/L 0 0 ; 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2; 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L; -E*A/L 0 0 E*A/L 0 0 ; 0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2; 0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L]; fle=L*[qx/2 qy/2 qy*L/12 qx/2 qy/2 -qy*L/12]'; G=[n(1) n(2) 0 0 0 0; -n(2) n(1) 0 0 0 0; 0 0 1 0 0 0; 0 0 0 n(1) n(2) 0; 0 0 0 -n(2) n(1) 0; 0 0 0 0 0 1]; Ke=G'*Kle*G; fe=G'*fle; Chương trình tính [k] trong hệ tọa độ riêng, ký hiệu Ke, và ma trận chuyển ký hiệu G với tên gọi Beam2t dùng cho dầm kiểu Timoshenko, trong đó ảnh hưởng ứng suất cắt thể hiện trong hệ số ks. 103 Hình 4. Dầm Timoshenko Công thức tính ma trận [k] có dạng: [ ] ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ + + +−+ + − +−+ + ++−++ − = μ μ μμ μ μ μμ μ μμμμ 1 2/14 1 6 1 12 00 1 2/12 1 60 1 2/14 1 6 1 120 1 6 1 12 0000 23 2 2323 L EJ L EJ L EJDX L AE L EJ L EJ L EJ L EJ L EJ L EJ L EJ L AE L AE k với sGAkL EJ 2 12=μ function [Ke,fe]=beam2t(ex,ey,ep,eq) % function [Ke,fe]=beam2t(ex,ey,ep,eq) % ep = [E G A I ks ] element properties, % E: Young's modulus % G: shear modulus % A: Cross section area % I: Moment of inertia % ks: shear correction factor if nargin==3; eq=[0 0]; end b=[ex(2)-ex(1);ey(2)-ey(1)]; L=sqrt(b'*b); n=b/L; E=ep(1); Gm=ep(2); A=ep(3); I=ep(4); ks=ep(5); % m=(12/L^2)*(E*I/(Gm*A*ks)); % Kle=E/(1+m)*[A*(1+m)/L 0 0 -A*(1+m)/L 0 0; 0 12*I/L^3 6*I/L^2 0 -12*I/L^3 6*I/L^2; 0 6*I/L^2 4*I*(1+m/4)/L 0 -6*I/L^2 2*I*(1-m/2)/L; 104 -A*(1+m)/L 0 0 A*(1+m)/L 0 0; 0 -12*I/L^3 -6*I/L^2 0 12*I/L^3 - 6*I/L^2; 0 6*I/L^2 2*I*(1-m/2)/L 0 -6*I/L^2 4*I*(1+m/4)/L]; % fle=L*[eq(1)/2 eq(2)/2 eq(2)*L/12 eq(1)/2 eq(2)/2 -eq(2)*L/12]'; % G=[n(1) n(2) 0 0 0 0; -n(2) n(1) 0 0 0 0; 0 0 1 0 0 0; 0 0 0 n(1) n(2) 0; 0 0 0 -n(2) n(1) 0; 0 0 0 0 0 1]; % Ke=G'*Kle*G; fe=G'*fle; Chương trình tính [k] trong hệ tọa độ riêng, ký hiệu Ke, và ma trận chuyển ký hiệu G với tên gọi Beam2w dùng cho dầm trên nền đàn hồi. Trong công thức tính ma trận cứng hệ số tải nền đàn hồi theo hướng pháp tuyến ký hiệu kt, theo hướng dọc trục ký hiệu ka. Ma trận cứng thứ nhất dầm trên nền đàn hồi với hệ số nền đàn hồi kt và ka, ký hiệu được viết như sau: ][ esk [ ] ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − = 2 22 4 22156 00140 31304 1354022156 007000140 420 Lk LkkDX k LkLkL LkkLkk kk Lk t tt a tt tttt aa e s Ma trận cứng của dầm trên nền là tổng của [k], thành phần thứ nhất như đã tính tại beam2e cùng [kse] như đang đề cập. Công thức xác định Ke - ma trận cứng dầm trên nền đàn hồi = [ ] [ ]( )eskk + . function [Ke,fe]=beam2w(ex,ey,ep,eq) % Ke=beam2w(ex,ey,ep) % [Ke,fe]=beam2w(ex,ey,ep,eq) % b=[ ex(2)-ex(1); ey(2)-ey(1) ]; L=sqrt(b'*b); n=b/L; % E=ep(1); A=ep(2); I=ep(3); ka=ep(4); kt=ep(5); % qx=0; qy=0; if nargin>3; qx=eq(1); qy=eq(2); end % K1 =[E*A/L 0 0 -E*A/L 0 0 ; 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2; 105 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L; -E*A/L 0 0 E*A/L 0 0 ; 0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2; 0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L]; % K2=L/420*[140*ka 0 0 70*ka 0 0 ; 0 156*kt 22*kt*L 0 54*kt -13*kt*L ; 0 22*kt*L 4*kt*L^2 0 13*kt*L -3*kt*L^2; 70*ka 0 0 140*ka 0 0 ; 0 54*kt 13*kt*L 0 156*kt -22*kt*L ; 0 -13*kt*L -3*kt*L^2 0 -22*kt*L 4*kt*L^2]; % Kle=K1+K2; % fle=L*[qx/2 qy/2 qy*L/12 qx/2 qy/2 -qy*L/12]'; % G=[n(1) n(2) 0 0 0 0; -n(2) n(1) 0 0 0 0; 0 0 1 0 0 0; 0 0 0 n(1) n(2) 0; 0 0 0 -n(2) n(1) 0; 0 0 0 0 0 1]; % Ke=G'*Kle*G; fe=G'*fle; Những hàm giúp xác định lực và chuyển vị trong hệ tọa độ cục bộ gắn liền phần tử đi theo ba hàm trên có tên tương ứng beam2s, beam2ts, beam2ws giới thiệu tiếp dưới đây. Dữ liệu đầu ra của beam2s và beam2ts: es = [N V M]; edi = [u v], của hàm beam2ws là es = . ⎥⎦ ⎤⎢⎣ ⎡ 222 111 MVN MVN Hàm beam2s % es=beam2s(ex,ey,ep,ed) % es=beam2s(ex,ey,ep,ed,eq) % [es,edi,eci]=beam2s(ex,ey,ep,ed,eq,n) EA=ep(1)*ep(2); EI=ep(1)*ep(3); b=[ ex(2)-ex(1); ey(2)-ey(1) ]; L=sqrt(b'*b); if length(ed(:,1)) > 1 disp('Only one row is allowed in the ed matrix !!!') return end qx=0; qy=0; if nargin>4; qx=eq(1); qy=eq(2); end ne=2; if nargin>5; ne=n; end; C=[0 0 0 1 0 0; 0 0 0 0 0 1; 0 0 0 0 1 0; L 0 0 1 0 0; 0 L^3 L^2 0 L 1; 0 3*L^2 2*L 0 1 0]; 106 n=b/L; G=[n(1) n(2) 0 0 0 0; -n(2) n(1) 0 0 0 0; 0 0 1 0 0 0; 0 0 0 n(1) n(2) 0; 0 0 0 -n(2) n(1) 0; 0 0 0 0 0 1]; M=inv(C)*(G*ed'-[0 0 0 -qx*L^2/(2*EA) qy*L^4/(24*EI) qy*L^3/(6*EI)]' ); A=[M(1) M(4)]'; B=[M(2) M(3) M(5) M(6)]'; x=[0:L/(ne-1):L]'; zero=zeros(size(x)); one=ones(size(x)); u=[x one]*A-(x.^2)*qx/(2*EA); du=[one zero]*A-x*qx/EA; v=[x.^3 x.^2 x one]*B+(x.^4)*qy/(24*EI); % dv=[3*x.^2 2*x one zero]*B+(x.^3)*qy/(6*EI); d2v=[6*x 2*one zero zero]*B+(x.^2)*qy/(2*EI); d3v=[6*one zero zero zero]*B+x*qy/EI; N=EA*du; M=EI*d2v; V=-EI*d3v; es=[N V M]; edi=[u v]; eci=x; Hàm beam2ts function [es,edi,eci]=beam2ts(ex,ey,ep,ed,eq,n) % [es,edi,eci]=beam2ts(ex,ey,ep,ed,eq,n) if nargin==5; n=2; end; if nargin==4; n=2; eq=[0 0]; end; ne=n; % EA=ep(1)*ep(3); EI=ep(1)*ep(4); GAK=ep(2)*ep(3)*ep(5); alfa=EI/GAK; % b=[ex(2)-ex(1);ey(2)-ey(1)]; L=sqrt(b'*b); n=b/L; % qx=eq(1);qy=eq(2); % C=[0 0 0 1 0 0; 0 0 0 0 0 1; 0 6*alfa 0 0 1 0; L 0 0 1 0 0; 0 L^3 L^2 0 L 1; 0 3*(L^2+2*alfa) 2*L 0 1 0]; % G=[n(1) n(2) 0 0 0 0; -n(2) n(1) 0 0 0 0; 0 0 1 0 0 0; 0 0 0 n(1) n(2) 0; 107 0 0 0 -n(2) n(1) 0; 0 0 0 0 0 1]; % M=inv(C)*(G*ed'-[0 0 0 -qx*L^2/(2*EA) qy*L^4/(24*EI)-qy*L^2/(2*GAK) qy*L^3/(6*EI)]' ); C2=[M(1) M(4)]'; C4=[M(2) M(3) M(5) M(6)]'; % x(1:ne,1)=[0:(L/(ne-1)):L]'; % one=ones(size(x)); zero=zeros(size(x)); % u=[x one]*C2-qx/(2*EA)*x.^2; du=[one zero]*C2-qx*x/EA; % v=[x.^3 x.^2 x one]*C4+qy/(24*EI)*(x.^4)-qy/(2*GAK)*(x.^2); dv=[3*x.^2 2*x one zero]*C4+qy*(x.^3)/(6*EI)-qy*x/GAK; % teta=[3*(x.^2+2*alfa*one) 2*x one zero]*C4+qy*(x.^3)/(6*EI); dteta=[6*x 2*one zero zero]*C4+qy*(x.^2)/(2*EI); % N=EA*du; M=EI*dteta; V=GAK*(dv-teta); % es=[N V M]; edi=[u v teta]; eci=[x]; Hàm beam2ws function es=beam2ws(ex,ey,ep,ed,eq) % es=beam2ws(ex,ey,ep,ed,eq) if length(ed(:,1)) > 1 disp('Only one row is allowed in the ed matrix !!!') return end b=[ ex(2)-ex(1); ey(2)-ey(1) ]; L=sqrt(b'*b); n=b/L; if nargin==4; qx=0; qy=0; end if nargin==5; qx=eq(1); qy=eq(2); end E=ep(1); A=ep(2); I=ep(3); ka=ep(4); kt=ep(5); K1=[ E*A/L 0 0 -E*A/L 0 0 ; 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2; 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L ; -E*A/L 0 0 E*A/L 0 0 ; 0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2; 0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L ]; K2=L/420*[140*ka 0 0 70*ka 0 0 ; 0 156*kt 22*kt*L 0 54*kt -13*kt*L ; 0 22*kt*L 4*kt*L^2 0 13*kt*L -3*kt*L^2; 70*ka 0 0 140*ka 0 0 ; 108 0 54*kt 13*kt*L 0 156*kt -22*kt*L ; 0 -13*kt*L -3*kt*L^2 0 -22*kt*L 4*kt*L^2]; Kle=K1+K2; fle=L*[qx/2 qy/2 qy*L/12 qx/2 qy/2 -qy*L/12]'; G=[n(1) n(2) 0 0 0 0; -n(2) n(1) 0 0 0 0; 0 0 1 0 0 0; 0 0 0 n(1) n(2) 0; 0 0 0 -n(2) n(1) 0; 0 0 0 0 0 1]; P=(Kle*G*ed'-fle); es=[-P(1) -P(2) -P(3) P(4) P(5) P(6)]; Ví dụ sử dụng phần tử BEAM 2D. Ví dụ 1: Ứng dụng giải bài toán dầm liên tục hai nhịp, ngàm đầu phía phải, tựa trên hai gối, một gối tại đầu phía trái, gối thứ hai chính giữa dầm, hình . Dầm được chia làm ba phần tử, PT 1 dài l, phần tử tiếp theo PT 2 dài l/2, và phần tử cuối dài l/2. Phần tử thích hợp dùng trong trường hợp này là phần tử dầm trong mặt phẳng, thường ký hiệu bằng 2D BEAM. Phần tử chịu uốn, mỗi nút có hai bậc tự do. Hình 6 Xác lập ma trận cứng các phần tử. [ ] ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − − − = )1( 22 )1( 21 )1( 12 )1( 11 2 22 )1( 4 612 264 612612 kk kk DX ll l llll l EJk 109 [ ] [ ] ⎥⎦ ⎤⎢⎣ ⎡=⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − − − == )3( 44 )3( 43 () 34 )2( 33 )2( 33 )2( 32 )2( 23 )2( 22 2 22 )3()2( 8 2496 2248 24962496 kk kk kk kk DX ll l llll l EJkk Tập họp ma trận cứng [K] = ∑k(e): [ ] ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ + += )3( 44 )3( 43 )3( 34 )3( 33 )2( 33 )2( 32 )2( 23 )2( 22 )1( 22 )1( 21 )1( 12 )1( 11 kk kkkk kkkk kk K Vector chuyển vị nút các phần tử xác định như sau: { } ;0 11 1)1( ⎭⎬ ⎫ ⎩⎨ ⎧= ⎭⎬ ⎫ ⎩⎨ ⎧= θθδ v { } ;0 22 2)2( ⎭⎬ ⎫ ⎩⎨ ⎧= ⎭⎬ ⎫ ⎩⎨ ⎧= θθδ v { } ; 3 3)3( ⎭⎬ ⎫ ⎩⎨ ⎧= θδ v { } ; 0 0)4( ⎭⎬ ⎫ ⎩⎨ ⎧=δ Vector ngoại lực xác định dưới dạng: { } ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − − − = 12/ 2/ 12/ 2/ 2 2 )1( ql ql ql ql P Lực tác động tại nút số 3: = P = ql/2. Tập họp các vecto lực tính cho tất cả 3 phần tử cấu thành dầm sẽ được vector lực toàn hệ {P}. Phương trình cân bằng giờ đây có dạng: ⎪⎪ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ − + + = ⎪⎪ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − −−− − −−− − −−− − − 4 4 12 12 22 12 12 4 4 3 3 2 2 1 1 2424 24962496 2424 24961922496 24186 2418108612 66 612612 0 84 41604 0 4122 2 2 2 2 22 222 22 22 M R R R l EJ ql ql ql ql ql ll llll ll lllll lll lllll ll llll θ θ θ θ v v v v Sau xử lý điều kiện biên, phương trình trở thành: 110 ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ −= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − 0 2 12 12 v 1604 019224 424122 24 2 2 3 3 2 1 2 ql ql ql ll l l EJ θ θ θ Từ đó có thể tính: EJ ql l 3 3 3 2 1 0022,0 0015,0 0089,0 0253,0 v ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − − = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ θ θ θ Vector nội lực tính cho các nút mỗi phần tử có dạng: Phần tử thứ nhất. ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − − = ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ − − − + ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− − − = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ l ql ql ql ql ql ll llll ll llll l EJ M N M N 983,0 5984,0 0 4016,0 12 2 12 2 0 0 4626 612612 2646 612612 2 2 2 1 22 22 2 2 1 1 θ θ Phần tử thứ hai: ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − − − = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− − − = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ l l ql ll llll ll llll l EJ M N M N 054,0 3048,0 0984,0 3048,0 v 0 824424 24962496 424824 24962496 3 3 2 22 22 2 2 1 1 θ θ Phần tử thứ 3: 111 ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ −= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− − − = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ l l ql v ll llll ll llll l EJ M N M N 0448,0 1968,0 0536,0 1968,0 0 0 824424 24962496 424824 24962496 3 3 22 22 2 2 1 1 θ Ví dụ 2: Sử dụng phần tử BEAM 2D, chịu uốn và chịu cắt dùng tính khung phẳng thường gặp trong kết cấu phương tiện vận tải. Kết cấu cứng trên tàu có kích thước như sau: l12 = 3m = l; l 23 = 6,6m = 2,2 l; l24 = 9m = 3l. Momen quán tính mặt cắt ngang các dầm, tính bằng cm4: J12 = 0,6. 104 = J; J23 = 0,9.104 = 1,5J; J24 = 2,4.104 = 4J; Tải tác động lên kết cấu: q2 = 22,6 kN/m = q; q4 = 90,4 kN/m = 4q. Hình 7 Mô hình hoá kết cấu được trình bày tại hình 7. Ma trận cứng các phần tử, tính trong hệ tọa độ gắn liền với phần tử như sau. [ ] ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − − = 22 1211 2 10 22 1010 )1( 4 612 0010 2604 120612 00100010 k kk ll DX l lll l EJk 112 [ ] ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − − = 33 2322 2 10 22 1010 )2( 73,2 86,169,1 0010 36,186,1073,2 86,169,1086,169,1 00100010 k kk ll DX l llll l EJk Ma trận chuyển áp dụng cho phần tử số 2, góc nghiêng giữa hệ tọa độ cục bộ và hệ tọa độ chung α = π/2, có dạng: [ ] [ ] [ ] ;0 0 ⎥⎦ ⎤⎢⎣ ⎡= R R T với [ ] ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − = 100 001 010 R Ma trận cứng phần tử số 2 trong hệ tọa độ chung tính bằng biểu thức [T][k(2)][T]T: ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− 3332 2322 10 2 1010 22 73,2 010 86,1069,1 36,1086,173,2 0100010 86,1069,186,1069,1 kk kk DX ll l llll l EJ Với phần tử số 3: [ ] ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − − = 4442 2422 2 10 22 1010 )3( 33,5 67,278,1 0010 67,267,2033,5 67,278,1067,278,1 00100010 kk kk ll DX l llll l EJk Ma trận cứng toàn hệ tính trong hệ tọa độ chung: 113 [ ] [ ] ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ ++= )3( 44 )2( 33 )3( 24 )2( 23 )3( 22 )2( 22 )1( 22 )1( 12 )1( 11 0 00 kDX k kkkkk kk K Vector lực tính trong hệ tọa độ cục bộ, gắn liền với phần tử, cho các phần tử số 1 và 3, chịu tải phân bố theo luật tuyến tính mang dạng sau: { } ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ − = l l qlP 05,0 35,0 0 033,0 15,0 0 )1( ; { } ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ − = l l qlP 1,2 65,4 0 65,1 85,2 0 )3( Hệ phương trình {K].{δ}= {P}, trong đó bậc của [K] 8x8, hai vector còn lại bậc 8x1. Điều kiện biên: u1 = v1 = u2 = v2 = u3 = v3 = u4 = v4 = θz4 = 0. Áp đặt điều kiện biên cho hệ phương trình vừa nêu có thể viết lại phương trình ở dạng sau cùng: ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ 0 60,1 033,0 73,236,10 36,106,122 024 3 2 1 l l ql l EJ z z z θ θ θ Sau khi thay q = 22,6 kN/m; J = 6000cm4, l = 3m và E=2,058.105 MPa (2,058.1011 N/mm2), sẽ tìm được các giá trị góc xoay, tính bằng rad. 2 3 2 1 10. 375,0 754,0 336,0 − ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ − − = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ z z z θ θ θ Biểu đồ phân bố lực cắt tính bằng kN và momen uốn, tính bằng kN.m như sau. Phần tử 1: [ ] ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ −= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − + ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ 65,106 14,58 038,0 4,24 05,0 35,0 033,0 15,0 0 0 2 1)1( 2 2 1 1 l l qlk M N M N z z θ θ Phần tử 2: 114 [ ] ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ −= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ 0696,0 672,9 73,63 672,9 0 0 3 2)2( 2 2 1 1 z zk M N M N θ θ Phần tử 3: [ ] ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − − − = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − + ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ 00,510 89,342 20,170 61,165 1,2 65,4 65,1 85,2 0 0 0 2)3( 2 2 1 1 l l qlk M N M N zθ Ví dụ sử dụng Matlab giải bài toán uốn dầm trong khuôn khổ phương pháp PTHH. Chương trình beam.m sau đây là một trong nhiều cách thể hiện. echo on %----- Topology ------------------------------------------------- Edof=[1 1 2 3 4 5 6; 2 4 5 6 7 8 9; 3 7 8 9 10 11 12]; %----- Stiffness matrix K and load vector f --------------------- K=zeros(12); f=zeros(12,1); f(5)=-10000; %----- Element stiffness matrices ------------------------------ E=2.1e11; A=45.3e-4; I=2510e-8; ep=[E A I]; ex=[0 3]; ey=[0 0]; Ke=beam2e(ex,ey,ep) %----- Assemble Ke into K --------------------------------------- K=assem(Edof,K,Ke); %----- Solve the system of equations and compute reactions ------ bc=[1 0; 2 0; 11 0]; [a,Q]=solveq(K,f,bc); %----- Section forces ------------------------------------------- Ed=extract(Edof,a); es1=beam2s(ex,ey,ep,Ed(1,:)); es2=beam2s(ex,ey,ep,Ed(2,:)); es3=beam2s(ex,ey,ep,Ed(3,:)); %----- Results ----------------Hình 8---------------------------- a ; Q; es1; es2; es3 %------------------------ end ----------------------------------- 4m 6m E E I A A 75kN/m I C 1kN 1 1 1 2 A 1E 1 1I 22 echo off 115 es1=beam2s(ex1,ey1,ep1,Ed(1,:),eq1, 20) Ỵ es1 = 1.0e+05 * -2,2467 0,1513 0,1981 -2,2467 0,1513 0,1662 . . -2,2467 0,1513 -0,4070 Ví dụ 3: Sử dụng phần tử BEAM 2D giải bài toán ổn định dầm. Công biến dạng dầm trường hợp bị uốn và nén gồm hai thành phần: công biến dạng do uốn và do uốn - nén. Lực dọc trục tại phần này ký hiệu bằng T. ( )∫+= L dxwEJAELTU 0 2 2 " 2 1 2 1 (a) Công biến dạng do tác động lực nén dọc trục tính theo công thức: ( )∫−≈ LN dxwTAELTW 0 2 2 ' 2 1 2 1 (b) Thế năng của dầm trongtrường hợp này được hiểu là: ( ) ( )∫∫ +=−=Π LLN dxwTdxwEJWU 0 2 0 2 ' 2 1" 2 1 (c) Từ biểu thức thứ nhất phía phải phương trình chúng ta đã xây dựng công thức tính ma trận cứng [k]e, dạng: [ ] ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − = 2 22 3 4 612 264 612612 LDX L LlL LL L EJke (*) Ma trận “cứng” xây dựng từ biểu thức sau vế phải công thức được tiến hành theo cách tương tự. Ma trận thứ hai này ký hiệu [kg]e. Nếu ký hiệu [N] – hàm hình dáng, {u} – vector chuyển vị nút của phần tử, công thức tính đạo hàm hàm w = [N]{u} sẽ là: }]{[' uN dx dw = Sử dụng hàm Hermite trong xây dựng hàm hình dáng: 5 4 6 8 9 11 12 101 2 3 1 2 3 116 Hình 9 [ ] ⎥⎦ ⎤⎢⎣ ⎡ +−−+−+−= 2 2 3 3 22 2 3 3 2 326634166 L x L x L x L x L x L x L x L xN dx d (d) Từ đó có thể thấy: ( ) { } { }uNu dx Ndw ]'[][' == ; ( ) { } {uNNuw TT ]'[]'[' 2 = } Biểu thức tính công theo biểu thức thứ hai, vế phải công thức (c) được viết như sau: { } [ ] [ ] { }udxNNTu L TT ⎥⎦ ⎤⎢⎣ ⎡ ∫ ''21 Ma trận [kg ]e sẽ là: [ ] [ ] ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −− − ⎟⎠ ⎞⎜⎝ ⎛== ∫ 2 22 4 336 34 336336 30 ]'[' LDX L LLL LL L TdxNNTk L T eg (e) Ví dụ: Xác định tải giới hạn lên dầm độ cứng EJ = 1,2kN.m2, L = 750mm, và a = 500mm. C PA B Hình 10 Dầm được chia làm 2 phần tử, ghi số 1 – đoạn AC, 2 – đoạn CB. Số thứ tự các nút tại A, C, B sẽ là 1, 2, 3. Lực T trong công thức gán bằng –P như nêu tại hình 10. Độ cứng các phần tử tính theo công thức vừa nêu: 117 [k]1 = ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −− − − ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − 1 5,136 25,05,11 5,1365,136 15 1 312 5,031 312312 5,0 102,1 3 3 DX P DX x [k]2 = ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −− − − ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − 25,0 75,036 0625,075,025,0 75,03675,036 15 2 25,0 5,112 125,05,125,0 5,1125,112 25,0 102,1 3 3 DX P DX x Tập họp hai ma trận cứng từ hai phần tử nhận được ma trận cứng toàn hệ: [ ] [ ]( ) 6 5 4 3 2 1 1,08,4 1,01,0 1,08,402,7 001,0 001,04,21,04,2 10 2,19 2,1156,921 6,92,1158,28 2,1156,9214,868,1036 008,48,286,9 008,282,1158,282,115 ][ 30 1 120 1 60 1 15 1 3 2 1 ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − − − − ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − − − =+= − − = ∑ DX P DX kkK e ege Điều kiện biên bài toán xác định rằng w1 = w3 = 0. Sau xử lý điều kiện biên, trong trường hợp này các dòng số 1, số 3 và các cột 1, 3 bị loại khỏi ma trận cứng, dạng cuối cùng của [K] sẽ là: 2 4 5 6 2 4 5 6 [ ] 6 5 4 2 1,08,4 1,01,0 00 10 2,19 2,1156,921 6,92,1158,28 008,46,9 30 1 120 1 60 1 15 1 3 ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −− ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −= − − DX P DX K Có thể mượn biến sau đây thay vào vị trí của P nhằm tạo thuận lợi lúc giải bài toán và tăng độ chính xác phép tính: P = 12x104xp (**) Sau thay thế, định thức từ ma trận độ cứng có dạng: 118 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )ppp ppp pppp pp K 42,18122,1156,90 122,1155766,921122,1150 6,9122,115128,2828,4 0028,486,9 10)det( 3 −+−+ +−−+− ++−−+ +− = Trị riêng của bài toán được xác định như sau: p = 0,1314; 0,86549; 2,18863; 6,56217 Giá trị nhỏ nhất của p thay thế vào biểu thức (**) giúp xác định tải giới hạn: Pcr = 12x104x0,1314 = 15,77x103N = 15,77kN. Ví dụ 4: Ứng dụng phần tử chịu uốn và xoắn đồng thời Ma trận cứng phần tử 2D này tập họp từ BEAM, biểu thức (f) và (j) của TOR. [ ] ⎥⎦ ⎤⎢⎣ ⎡ − −= 11 11 L GIk T [ ] ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ −− − = EJ EJLEJLDX LGJ EJEJLJ EJLEJLEJLEJL EJLLGJ L k T T 12 64 00 126012 62064 0000 1 2 2 22 22 3 Phần tử nhóm này cần thiết khi xử lý các bài toán kết cấu giàn. Trong các kết cấu thuộc thân tàu thủy, ô tô, máy bay giàn gồm các dầm hoặc cơ cấu tương đương dầm, bố trí dọc và ngang. Thông thường hai hệ thống dầm vuông góc với nhau, trong trường hợp ấy chúng ta gặp hệ dầm trực giao. Giàn2 trong thân tàu thủy có dạng tiêu biểu, giới thiệu tại hình 11. Hình 11. 2 Giàn, còn có thể viết “dàn” trong tiếng Việt. Từ kỹ thuật tương ứng dùng trong sách bằng tiếng Anh là grillage, trong tiếng Pháp grigalle plan. Sách xuất bản tại Nga viết là перекрытие, còn người Trung Hoa gọi đây là 扳架. 119 Lực tác động lên giàn theo phương pháp tuyến với mặt giàn. Phần tử 2D beam trong trường hợp tính giàn phải chịu uốn và xoắn. Ví dụ sau đây trình bày cách sử dụng phần tử 2D BEAM tính kết cấu giàn đối xứng. Tại kết cấu trình bày tại hình có thể nhận xét, cấu hình giàn đối xứng qua trục Ox và cả trục Oy. Trong trường hợp này nên sử dụng ¼ giàn vào tính toán nhằm tiết kiệm công sứ