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 k se][ được viết như sau:
[ ]
⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣
−
−
−
=
2 2
2
4
156 22
140 00
31304
54022156 13
007000140
420
Lk
DX Lkk
k
L LkLk
Lkk Lkk
k k
L
k
t
t t
a
t t
t t t t
a a
es
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 = ([ ]+[kk se]).
44 trang |
Chia sẻ: trungkhoi17 | Lượt xem: 456 | Lượt tải: 0
Bạn đang xem trước 20 trang tài liệu 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, để xem tài liệu hoàn chỉnh bạn click vào nút DOWNLOAD ở trên
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ứ
Các file đính kèm theo tài liệu này:
- giao_trinh_bao_toan_nang_luong_chuong_1_cac_phuong_phap_nang.pdf