Xây dựng ma trận độ cứng phần tử tấm – Gân ứng dụng trong tính toán kết cấu tấm composite lớp có gân tăng cứng bằng phương pháp phần tử hữu hạn

T¹p chÝ Khoa häc & C«ng nghÖ - Sè 2(42)/N¨m 2007 
 29
XÂY DỰNG MA TRẬN ĐỘ CỨNG PHẦN TỬ TẤM – GÂN ỨNG DỤNG 
TRONG TÍNH TOÁN KẾT CẤU TẤM COMPOSITE LỚP 
CÓ GÂN TĂNG CỨNG BẰNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN 
 Ngô Như Khoa (Trường Đại học Kỹ thuật Công nghiệp - ĐH Thái Nguyên) - 
Đỗ Tiến Dũng (Trường Cao đẳng Công nghiệp Việt - Hung) 
1. Giới thiệu 
Nhờ có ưu điểm nổi trội về khả năng chịu lực trong khi chi phí về vật liệu và trọng lượng 
kết cấu được giảm ở mức đáng kế, mà các kết cấu tấm-vỏ có gân tăng cứng đã được sử dụng rất 
phổ biến ở hầu hết các ứng dụng trong kỹ thuật và đời sống, cho dù là các kết cấu chế tạo từ các 
loại vật liệu kinh điển hay các kết cấu được chế tạo từ vật liệu composite lớp. Tuy nhiên, trong 
thực tế của ngành cơ học kỹ thuật, việc tính toán cơ học đối với các kết cấu tấm-vỏ có gân tăng 
cứng luôn được xem là rất phức tạp và cho đến nay vẫn chưa có được lời giải tổng quát, đặc biệt 
là các kết cấu bằng vật liệu có tính dị hướng cao như composite lớp. Vì vậy, vấn đề này đã và 
đang được sự quan tâm của nhiều nhà nghiên cứu cơ học trong và ngoài nước. Ví dụ như, gần 
đây Kolli và Chandrashekhara [3] sử dụng phần tử đẳng tham số với các hàm nội suy khác nhau 
cho tấm và dầm để phân tích ứng xử phi tuyến của tấm gân Composite bằng việc sử dụng phần 
tử tứ giác 9 nút và phần tử gân 3 nút dựa trên lý thuyết tấm của Mindlin. Các tác giả Y.V.Satish 
Kumar, Madhujit Mukhopadhyay[4] sử dụng một phần tử tấm gân mới để phân tích ổn định cho 
kết cấu tấm có gân tăng cứng bằng vật liệu composite lớp, phần tử này là một sự tổ hợp của phần 
tử tam giác ứng suất phẳng của Allman và một phần tử uốn Mindlin –Kirchhoff rời rạc; mô hình 
này cũng có khả năng áp dụng đối với bài toán có số gân bất kỳ và hướng tuỳ ý. Nhóm tác giả 
Guanghui Qing, Jiajun Qiu, Yanhong Liu [5] dựa trên nghiệm bán giải tích của lý thuyết phương 
trình véctơ trạng thái, một mô hình toán học mới để phân tích dao động tự do của tấm gân nhiều 
lớp đã được phát triển bằng cách xem xét riêng biệt các phần tử tấm và gân; phương pháp này 
dựa trên điều kiện tương thích về ứng suất và biến dạng tại các điểm nút giao tiếp giữa tấm và 
gân; các tác giả cũng sử dụng phần tử tứ giác bậc nhất 4 nút và với phạm vi nghiên cứu giới hạn 
trong các kết cấu có gân bố trí dọc theo các cạnh của tấm. Bên cạnh một số công trình quốc tế đã 
công bố trên như đã liệt kê trên, gần đây cũng đã có một số công trình trong nước như: Nhóm 
các tác giả Trần Ích Thịnh, Trần hữu Quốc [2] đã nghiên cứu bài toán dao động của các kết cấu 
tấm có gân tăng cứng bằng vật liệu Composite; ở nghiên cứu này, các tác giả đã sử dụng mô 
hình phần tử tứ giác đẳng tham số 9 nút và phần tử dầm 3 nút độc lập với cùng hàm nội suy; ma 
trận độ cứng của phần tử dầm được xây dựng dựa trên điều kiện tương thích về chuyển vị tại mặt 
liên kết giữa tấm và gân trên cơ sở lý thuyết tấm của Mindlin; tương tự như các nghiên cứu của 
các tác giả khác, nghiên cứu này cũng chỉ khảo sát với các kết cấu có gân bố trí dọc theo các 
cạnh, hay việc chia lưới phải phụ thuộc vào sơ đồ bố trí của gân. 
Mục đích được đặt ra trong báo cáo này là xây dựng được mô hình phần tử có thể áp 
dụng cho bài toán kết cấu tấm có gân tăng cứng ở dạng tổng quát (kết cấu có số lượng gân bất 
kỳ, hướng gân không nhất thiết phải song song với các cạnh bên của tấm). Tư tưởng chính để 
thực hiện trong báo cáo là rời rạc hoá kết cấu bởi các phần tử dạng tam giác bậc hai, trong đó 
bao gồm các phần tử đơn thuần là phần tử tấm chịu uốn truyền thống và các phần tử có sự tổ hợp 
với thành phần gân. Ở các phần tử tổ hợp ta xem phần tử tấm và phần gân là hai thành phần của 
T¹p chÝ Khoa häc & C«ng nghÖ - Sè 2(42)/N¨m 2007 
 30
một thể thống nhất, như vậy ma trận độ cứng của phần tử tổ hợp sẽ là tổng ma trận độ cứng của 
các thành phần. Tuy nhiên, điểm mấu chốt ở đây đó là ma trận độ cứng của thành phần gân được 
xây dựng trên cơ sở của việc biểu diễn trường biến dạng trong gân thông qua một trường chuyển 
vị trung gian lấy trên phần tử tấm, và trường chuyển vị này được xác định nhờ việc nội suy từ 
các thành phần chuyển vị nút của phần tử tấm. 
2. Ứng xử cơ học của tấm và dầm composite lớp 
* Mô hình bài toán: 
Mô hình bài toán tấm có 
gân tăng cứng có thể được biểu 
diễn như hình bên (Hình 1). Để 
đơn giản hoá việc biểu diễn, 
chúng tôi sử dụng cách biểu diễn 
với vật liệu đơn lớp và không mô 
tả điều kiện liên kết, tuy nhiên liên 
kết của kết cấu cũng sẽ được khảo 
sát ở dạng tổng quát. 
Trong mô hình tổng quát, 
hệ trục chung lấy theo hệ quy 
chiếu của tấm (x,y,z) như hình vẽ, hệ trục địa phương của gân là (x’,y’,z’). Trong đó trục z’≡ z, 
x’≡ phương của gân và hệ trục quy đổi của gân ( )zyx ,, là hệ trục (x,y,z) tịnh tiến đi 1 khoảng 





 +
2
,0,0 Hh , với h và H tương ứng là chiều dày của tấm và chiều cao của gân. 
* Ứng xử cơ học của tấm composite lớp 
Quan hệ ứng suất-biến dạng trong tấm đối với lớp thứ k được biểu diễn như sau: 
kxz
yz
xy
y
x
kkxz
yz
xy
y
x
QQ
QQ
QQQ
QQQ
QQQ
































=
















γ
γ
γ
ε
ε
τ
τ
τ
σ
σ
5545
4544
662616
262212
161211
000
000
00
00
00
 (1) 
Trong đó, Qij là các hệ số độ cứng thu gọn, được xác định từ các hằng số đàn hồi của vật 
liệu [1]; và các thành phần biến dạng trong (1) được xác định qua chuyển vị như sau: 
xxx zk+=
0εε ; yyy zk+=
0εε ; xyxyxy zk+=
0γγ ; 0yzyz γ=γ ; 0xzx γ=γ ; (2) 
trong đó: 
x
u
x ∂
∂
=
0
0ε ;
x
v
y ∂
∂
=
0
0ε ; 
x
k xx ∂
∂
=
θ
; 
y
k yy ∂
∂
=
θ
; 
xy
k yxxy ∂
∂
+
∂
∂
=
θθ
; 
x 
pxy 
z 
y bg2 
Hg2 
Hg1 
bg1 
Hình 1. Mô hình kết cấu tấm chịu uốn 
có gân tăng cứng 
T¹p chÝ Khoa häc & C«ng nghÖ - Sè 2(42)/N¨m 2007 
 31
x
v
y
u
xy ∂
∂
+
∂
∂
=
00
0γ ; 
y
w
yyz ∂
∂
+= θγ ; 
x
w
xxz ∂
∂
+= θγ (3) 
Trong (3), u, v và w là các chuyển vị ở mặt giữa của tấm và xθ , yθ là các góc xoay. Tích 
phân (1) theo chiều dầy tấm, ta thu được phương trình quan hệ ứng suất - biến dạng: 
 { } { }ttt DN 0ε





= (4) 
Trong đó { }tN là vectơ lực màng và mô men uốn, xoắn; 




 tD là ma trận độ cứng vật 
liệu, được xác định từ các hằng số đàn hồi của vật liệu [1]: 
{ } { }Txzyzxyyxxyyxt kkk 000000 ,,,,,,, γγγεεε = (5) 
[ ] ( )Txyyxxzyzxyyxt MMMQQNNNN , , , , , , ,= (6) 
11 12 16 11 12 16
12 22 26 12 22 26
16 26 66 16 26 66
11 12 16 11 12 16
12 22 26 12 22 26
16 26 66 16 26 66
44 45
45 55
0 0
0 0
0 0
0 0
0 0
0 0
0 0 0 0 0 0
0 0 0 0 0 0
t t t t t t
t t t t t t
t t t t t t
t t t t t t
t
t t t t t t
t t t t t t
t t
t t
A A A B B B
A A A B B B
A A A B B B
B B B D D D
D
B B B D D D
B B B D D D
A A
A A





  =  










 
 
 
 
 (7) 
trong đó: 
( ) ( ) ( )∑∫
=
−
=
L
k
h
h k
t
ij
t
ij
t
ij
t
ij
k
k
dzzzQDBA
1
2
1
,,1,, với i=1,2,6 (8) 
( ) ( ) ( ) ( )[ ]∑∫
=
−
=
L
k
h
h k
t
k
t
k
tttt k
k
dzQfQffQfAAA
1
55
2
2452144
2
1554544
1
,,,, (9) 
Trong đó tijQ là các hệ số vật liệu của lớp thứ k trong hệ toạ độ tấm, f1 và f2 là các hệ số 
hiệu chỉnh cắt và L là số lớp vật liệu trong tấm. 
* Quan hệ ứng suất-biến dạng của dầm composite lớp 
Các gân gia cường được mô hình bằng các dầm nhiều lớp và đặt theo phương bất kỳ đối 
với cạnh tấm. Dầm chỉ chịu uốn, khi đó thành phần chuyển vị v(x,z) là vô cùng bé nên bỏ qua, 
trong tính toán chỉ còn 2 thành phần chuyển vị u(x,z) và w(x) [6]. 
Trường chuyển vị theo Mindlin được biểu diễn như sau: 
( ) ( ) ( )xzxuz,xu xθ+= 0 (10.a) 
( ) ( )xwxw 0= (10.b) 
T¹p chÝ Khoa häc & C«ng nghÖ - Sè 2(42)/N¨m 2007 
 32
Trường biến dạng của dầm: 
xxx zk+= 0εε ; 0xzxz γγ = (11) 
Với: 
x
u
x ∂
∂
=
0
0ε ;
x
k xx ∂
∂
=
θ
; xxz
x
w
θγ +
∂
∂
=
0
0
 (12) 
Trường ứng suất: 
kxz
x
k
'
'
'
'
'
kxy
xz
yz
y
x
C
C
C
C
C




























=


















γ
ε
σ
σ
σ
σ
σ
0
0
0
0
0
16
55
45
12
11
 (13) 
3. Hệ phương trình phần tử hữu hạn 
Sử dụng phần tử tam giác 6 nút, mỗi nút 5 bậc tự do để mô hình tấm, ta có thể viết các 
chuyển vị tại điểm bất kỳ trong tấm như sau: 
 { } ( )[ ] ( ){ }ty,xN ei
6
1i
e
i uu
e ∑
=
= (14) 
trong đó: { } { }Tyixi0i0i0iei ,,w,v,u θθ=u ; [ ] [ ]Ieiei NN = 
và eiN là các hàm nội suy Lagrange, [I] là ma trận đơn vị 5×5 và ui, vi,... là các chuyển vị nút. 
Ta viết lại (14) dưới dạng: 
 { } { }N qe e ei =  u (15) 
trong đó: [ ] [ ] [ ] [ ]e6e2e1ei N,...,N,NN = 
 { } { } { } { }( )TeTeTee 921 ,...,,q uuu= 
Thay (15) vào (3) ta được 
 { } [ ]{ }eete qBε = (16) 
Để mô hình cho gân, ta xét trường hợp tổng quát, kết cấu tấm có gân tăng cứng hợp với 
trục x góc ϕ , có giá trị giới hạn là 01800 ≤≤ ϕ . 
Đối với kết cấu tấm có gân tăng cứng, khi ta rời rạc hoá bởi các phần tử, sẽ tồn tại các 
phần tử có chứa đoạn gân tăng cứng như hình 2. 
T¹p chÝ Khoa häc & C«ng nghÖ - Sè 2(42)/N¨m 2007 
 33
 Véc tơ chuyển vị của điểm M bất kỳ trên mặt trung bình của gân sẽ được xác định thông 
qua véc tơ chuyển vị Mt (M là hình chiếu của Mt theo phương vuông góc với mặt phẳng tấm) 
Khi đó, véc tơ chuyển vị tại Mt : 
{ } { }Ttytxtttet wvuu θθ ,,,, 000= (17) 
được nội suy bởi các thành phần chuyển vị nút của phần tử tấm như sau: 
{ } { } [ ]{ }eeet qNq
NN
NN
NN
NN
NN
u =
















=
61
61
61
61
61
00000000
00000000
00000000
00000000
00000000
⋯
⋯
⋯
⋯
⋯
 (18) 
Các chuyển vị của gân trong hệ trục (x’, y’, z’) được xác định theo các chuyển vị của tấm 
theo hệ trục (x, y, z), nhờ đó mà điều kiện tương thích về chuyển vị của phần tử thanh và phần tử 
tấm tự động được thoả mãn [4]. 
{ } [ ][ ]et
e
t
y
t
x
t
t
t
e
'y
'x
'ex
uEw
v
u
'w
'v
'u
=
































=


















=
φ
θ
δ
δ
θ
θ
0
0
0
0
0
0
10000
01000
00
0
1
0
0
1
0
0
0001
u
 (19) 
trong đó: ( ) 2hH +=δ là khoảng cách giữa mặt phẳng trung bình tấm với mặt phẳng trung bình gân. 
Biến dạng của gân trong hệ trục toạ độ địa phương ( )zyx ,, là: 
 { } { }
zxx
Tg γεε ,= (20) 
Các thành phần biến dạng trên gồm có: { } { }0
,
00
,,
zxxxx
Tg γθεε = 
1 2 
3 
4 
5 6 
Mt 
Mt 
M 
Mặt trung bình của tấm 
Mặt trung bình của gân 
Hình 2. Mô hình phần tử tấm chứa một đoạn gân 
x 
x’ 
y 
 ϕ 
T¹p chÝ Khoa häc & C«ng nghÖ - Sè 2(42)/N¨m 2007 
 34
Các thành phần biến dạng trong hệ trục ( )zyx ,, được chuyển đổi sang hệ trục (x’, y’, z’) như 
sau [1]: 
{ } [ ]{ }'00 gxgg T εε = (21) 
trong đó: { }'0gxε là các thành phần biến dạng của { }g0ε được biểu diễn trong hệ toạ độ (x’, y’, z’) 
theo biểu thức: 
{ } { }0
''
0
''','','','','
0
''
0
'
0
'
'0
,,,,,,,, zxzyyxxyyyxxyxyx
Tgx γγθθθθγεεε = (22) 
và ma trận chuyển đổi hệ trục biến dạng : 
{ }


















=
ϕϕ
ϕϕϕϕ
ϕϕϕ
cossin
sinsinsincos
sinsincos
T g
0000000
002
2
12
2
1000
0000002
2
1
22
22
 (23) 
Biểu thức (22) biểu diễn qua chuyển vị { }'uex dưới dạng: 
{ } [ ] 'xg'gx dL=0ε
 (24) 
Trong đó các ma trận toán tử Lg được xác định như sau: 
T
g
xy
yx
xy
xy
yx
L






























∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
=
010
''
0000
10
'
00
'
000
''
0000000
000000
''
0
000000
'
0
'
 (25) 
Cuối cùng ta có các thành phần biến dạng được viết lại dưới dạng: 
{ } [ ]{ } [ ][ ]{ } [ ]{ }aBaNLdL gg'xg'gx ===0ε
 (26) 
Trong đó các Bg là ma trận liên hệ biến dạng chuyển vị được xác định như sau: 
[ ] [ ][ ]NLB gg = (27) 
Biểu thức thế năng của vật thể được xác định là: 
WU +=Π (28) 
T¹p chÝ Khoa häc & C«ng nghÖ - Sè 2(42)/N¨m 2007 
 35
trong đó: 
dV
2
1U
V
Tσε∫= là năng lượng biến dạng đàn hồi, và pdSwpdSdW
SS
T
∫∫ −=−= là công của 
lực phân bố. Áp dụng nguyên lý cực tiểu hoá thế năng ta thiết lập được hệ phương trình PTHH: 
0=Πδ
 (30) 
Năng lượng biến dạng đàn hồi của tấm được xác định bởi: 
{ } { }1 12 2
Tt t t T
t tA
U dxdy = = ∫ ε D ε q K q (31) 
Năng lượng biến dạng đàn hồi của gân được xác định bởi: 
{ } [ ] [ ] [ ][ ][ ]{ } qKqxdqBTDTBqbU geT
g
gggTgTgTgg
e 2
1
2
=





= ∫
ℓ
 (32) 
Trong đó Dg là ma trận độ cứng của phần tử gân được biểu diễn như là: 
[ ]












=
g
gg
gg
g
A
DB
BA
D
11
1111
1111
'00
0
0
 (33) 
Ma trận độ cứng của tấm và gân được xác định là: 
[ ] [ ]dSB'ABBDBBBBBBBBABK
St
TTTTTt
e ∫ ++++= 3322122111 (34) 
trong đó: 
[ ] [ ] [ ] [ ][ ]NENL.......NLNLB 121111 = (35.a) 
[ ] [ ] [ ] [ ][ ]NENL.......NLNLB 222122 = (35.b) 
[ ] [ ] [ ] [ ][ ]NENL.......NLNLB 323133 = (35.c) 
[ ] [ ] ( )∑
=
−
−==
n
k
kkkijij hh'QAA
1
1
 (36.a) 
[ ] [ ] ( )∑
=
−
−==
n
k
kkkijij hh'QBB
1
2
1
2
2
1
 (36.b) 
[ ] [ ] ( )∑
=
−
−==
n
k
kkkijij hh'QDD
1
3
1
3
3
1
 với: i, j = 1, 2, 6. (36.c) 
[ ] [ ] ( )∑
=
−
−==
n
k
kkkijij hh'Q'A'A
1
1 với: i, j = 4, 5. (36.d) 
T¹p chÝ Khoa häc & C«ng nghÖ - Sè 2(42)/N¨m 2007 
 36
[ ] [ ] [ ] [ ] [ ] xdBTDTBbK
g
gggTgTg
g
g
e ∫=
ℓ
 (37) 
Ma trận độ cứng của phần tử tổ hợp tấm-gân được xác định bởi: 
g
e
t
ee KKK += (38) 
Nếu như phần tử gân chỉ thuộc 1 phần tử tấm thì ma trận độ cứng của phần tử tổ hợp thực 
hiện như biểu thức (38), nếu thuộc cả hai phần tử tấm thì coi như chỉ thuộc 1 phần tử và phần tử 
còn lại coi như không chứa gân. 
Công do tải trọng phân bố đều được xác định bởi: 
 [ ] { } FqdSyxpBqdSyxwyxp T
S
T
P
T
S ee
== ∫∫∫∫ ),(),(),( (39) 
trong đó: [ ]PB là ma trận chuyển đổi lực bề mặt, nó được xác định từ các ma trận định vị và các 
hàm dạng như sau: 
[ ] [ ][ ]621 .........NNNLB PP = (40) 
Trong trường hợp tổng quát, nếu coi các tải trọng tác dụng lên mặt trung bình của tấm, 
bao gồm các tải trọng tác dụng trong mặt phẳng tấm và kể cả tải trọng tác dụng vuông góc với bề 
mặt, thì Lp là ma trận đơn vị có kích thước (5×5). 
Thay (31), (32) và (39) vào (30) ta thu được: [ ]{ } { }FqK = , là hệ phương trình PTHH cần 
thiết lập. 
4. Kết luận 
Với mục đích xây dựng mô hình phần tử phục vụ cho bài toán kết cấu tấm có gân tăng 
cứng ở dạng tổng quát, khắc phục những hạn chế của những mô hình phần tử truyền thống, báo 
cáo này đã thiết lập được mô hình phần tử tổ hợp giữa phần tử tấm và thành phần gân chịu uốn. 
Mô hình này cho phép áp dụng đối với các kết cấu ở dạng tổng quát, cụ thể là có thể áp dụng đối 
với các kết cấu có số lượng gân bất kỳ, hướng gân không nhất thiết phải song song với các cạnh 
bên của tấm. Đây chính là tính ưu việt nổi bật của mô hình phần tử mới được xây dựng, có khả năng 
khắc phục những hạn chế của các mô hình PTHH dạng tấm và dạng gân truyền thống vẫn được sử 
dụng trong hầu hết các nghiên cứu về cơ học kết cấu composite trong và ngoài nước hiện nay. Ma 
trận độ cứng phần tử tổ hợp đã xây dựng đối với phần tử tam giác bậc 2, dựa trên lý thuyết tấm 
Mindlin. 
Việc xây dựng thuật toán PTHH và chương trình tính toán số để khảo sát các bài toán 
tĩnh và động các kết cấu tấm có gân tăng cứng bằng vật liệu composite lớp sẽ được nhóm tác giả 
tiếp tục trình bày trong các báo cáo khoa học tiếp theo
 
T¹p chÝ Khoa häc & C«ng nghÖ - Sè 2(42)/N¨m 2007 
 37
Tóm tắt 
 Bài báo này tập trung vào việc nghiên cứu xây dựng mô hình phần tử hữu hạn tổ hợp 
giữa phần tử tấm với phần tử gân truyền thống. Nhờ mô hình này, vấn đề tính toán và phân tích 
cơ học các kết cấu tấm composite lớp có gân tăng cứng bằng phương pháp phần tử hữu hạn 
(PTHH) có thể được giải quyết ở dạng tổng quát so với các dạng phần tử truyền thống khác. Nội 
dung chính được trình bày trong báo cáo là xây dựng ma trận độ cứng phần tử cho phần tử tấm 
tam giác 6 nút mà được tổ hợp với phần tử dầm, dựa trên tư tưởng chuyển vị của phần tử dầm 
được nội suy từ trực tiếp từ các nút của phần tử tấm. Trong đó, lý thuyết tấm được sử dụng là lý 
thuyết chuyển vị bậc nhất của Mindlin. 
Summary 
Building stiffness matrix of stiffened-plate element to applyi in laminated composite 
stiffened-plate analysis using finite element method 
 This paper concentrates on building finite element modeling, which is the combination 
of traditional plate element and beam plate. Based on this element, the calculation and analysis 
of laminated composite stiffened-plate mechanical by using finite element method (FEM) should 
be resolved in general forms. Main content of the paper is steps to build of stiffned matrix for 6 
node triangular elements; the new element is combinated of plate and beam elements. The 
solution is based on beam element’s displacement field is interpolated from plate element’s 
nodes. In which, the Mindlin plate theory was applied. 
Tài liệu tham khảo 
[1]. Trần Ích Thịnh (1994). Vật liệu composite- cơ học và tính toán kết cấu. Nxb Giáo dục, Hà 
Nội. 
[2]. Trần Ích Thịnh - Trần Hữu Quốc (2006), “Phân tích cơ học kết cấu tấm Composite lớp có 
gân gia cường bằng phương pháp phần tử hữu hạn”, Tuyển tập công trình khoa học hội nghị 
cơ học vật rắn biến dạng toàn quốc lần thứ 8, Thái Nguyên. 
[3]. M. Kolli - K. Chandrashekhara, Finite element analysis of siffened laminated plates under 
transverse loading, University of Missouri-Rolla, Rolla, Missouri 65409, USA. 
[4]. Y.V Satish Kumar - Madhujit Mukhopadhyay (1999), A new triangular stffened plate 
element for laminate analysis. Indian Institute of Technology, Kharagpur 721 302, India. 
[5]. Guanghui, Jiajun Qiu- Yanhong Liu,(2001) Free vibration analysis of stffened laminated 
plates. Tianjin University ,92 Weijin Road, Tianjin 300072, People’s Republic of China. 
[6]. Valery V. Vasiliev & Evgeny V. Morozov(2001), Mechanics and Analysis of Composite 
Materials. 
[7]. A.H Sheikha, A. Chakrabartib (2002), A new plate bendung element based on higher-order 
shear deformation theory for the analysis of composite plates. Engineering College, 
Jalpaiguri-735102, West Bengal, India.