Tấm và vỏ là các dạng kết cấu được sử dụng nhiều trong kỹ thuật
và chúng thường chịu biến dạng chịu uốn. Các phương trình PTHH đối
với các kết cấu tấm-vỏ thường phức tạp hơn nhiều so với các dạng kết
cấu khác. Chương 11 sẽ giới thiệu về hai lý thuyết tấm được sử dụng
phổ biến trong các bài toán kết cấu tấm-vỏ: lý thuyết tấm kinh điển của
Kirchoff (gọi tắt l à tấm Kirchoff) và lý thuyết tấm bậc nhất của Mindlin
(gọi tắt là tấm Mindlin).
Các thuật toán PTHH đối với tấm chịu uốn tương ứng với hai lý
thuy ết trên đã được thiết l ập chi tiết.
Phần tử vỏ được xem là tổ hợp của phần tử tấm chịu uốn và phần
tử tấm chịu trạng thái ứng suất phẳng.
299 trang |
Chia sẻ: banmai | Lượt xem: 2499 | Lượt tải: 1
Bạn đang xem trước 20 trang tài liệu Giáo trình kỹ thuật - Phương pháp phần tử hữu hạn, để xem tài liệu hoàn chỉnh bạn click vào nút DOWNLOAD ở trê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ú ý 0z 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ố 22
đ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 44 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 33; 44; 66 và 88. 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 2020 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 1010 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
Các file đính kèm theo tài liệu này:
- giaotrinh_phuong phap phan tu huu han.pdf