PHÂN TÍCH TĨNH VÀ DAO ĐỘNG TỰ DO CỦA TẤM COMPOSITE LAMINATE SỬ
DỤNG PHƯƠNG PHÁP PHẦN TỬ CHUYỂN ĐỘNG
STATIC AND FREE VIBRATION ANALYSIS OF LAMINATED COMPOSITE PLATES
USING MOVING ELEMENT METHOD
Lương Văn Hải, Phạm Hùng
laminate dưới tác động của các loại tải trọng trên các loại nền
thường được ứng dụng rộng rãi trong ngành xây dựng dân dụng
và cầu đường như: tính toán kết cấu móng, dầm, sàn của các tòa
nhà; thiết kế, tính toán kết cấu áo đường giao thông hay sân
bay,....
Trên thế giới có rất nhiều nghiên cứu về phân tích ứng xử
của kết cấu tấm composite laminate trong những thập kỷ gần
đây. Vào năm 1963, Thompson [1] đã nghiên cứu ứng xử của
kết cấu đường chịu tải trọng di chuyển bằng cách giả định
đường là tấm mỏng dài vô hạn tựa trên nền đàn hồi (resting
elastic foundation). Kim và Roesset (1998) [2] đã nghiên cứu
tấm vô hạn tựa trên nền đàn hồi Winkler chịu tải trọng di
chuyển sử dụng phép biến đổi Fourier. Wu và cộng sự (1987)
[3] đã sử dụng phương pháp phần tử hữu hạn để phân tích ứng
xử của tấm phẳng chịu tải trọng di chuyển. Pan và Atluri (1995)
[4] đã phân tích ứng xử của đường băng có kích thước hữu hạn
trên nền đàn hồi chịu tải trọng di chuyển sử dụng phương pháp
FEM/BEM kết hợp. Musharraf Zaman và cộng sự (1991) [5] đã
sử dụng phương pháp phần tử hữu hạn bốn nút để phân tích ứng
xử động của tấm dày trên nền đàn nhớt (viscoelastic foundation)
chịu tải trọng di chuyển. Như vậy, phương pháp phần tử hữu
hạn FEM truyền thống đã được sử dụng phổ biến để giải quyết
nhiều bài toán phức tạp. Trong phương pháp này, tất cả các ma
trận kết cấu sẽ được thực hiện trên một hệ trục tọa độ cố định.
Khi tải trọng di chuyển từ phần tử này sang phần tử khác thì
vectơ tải trọng phải được cập nhật sau mỗi bước thời gian. Do

(Moving Element Method) được phát triển nhằm phân tích ứng
xử tĩnh và dao động tự nhiên của tấm composite laminate dựa
trên lý thuyết biến dạng cắt bậc nhất. Tấm sẽ được chia nhỏ
thành những “phần tử chuyển động”. Những phần tử này không
phải chuyển động thật so với tấm đứng yên mà là chuyển động
giả tưởng cùng với lực tác động di chuyển trên kết cấu tấm. Do
đó, phương pháp này sẽ tránh được việc cập nhật vectơ tải trọng
so với mô hình tấm truyền thống. Các phương trình chuyển
động cũng như các ma trận của phần tử tấm được xây dựng trên
một hệ trục tọa độ chuyển động với vận tốc không đổi. Các kết
quả số trong phân tích tấm composite laminate cho thấy phương
pháp MEM đạt độ tin cậy cao so với các kết quả khác đã được
công bố trước đây.
Từ khóa: Phương pháp phần tử chuyển động, tấm composite
laminate, phân tích tĩnh, dao động tự do.
ABSTRACT
This paper develops a recently new method, namely moving
element method (MEM), for predicting the static and free
vibration response of laminated composite plates based on the
first shear deformation theory. In this method, the plate is
discretized into a number of elements called “the moving
elements”. These moving elements are not physical elements
fixed to the plate, but are conceptual elements that “flow” with
the moving load through the plate. Thus, the proposed method
eleminates the need of keeping track the location of moving
load relative to the element mesh. The governing equations of
motion as well as structural matrices of moving element are
formulated in a relative coordinate system travelling at a
constant speed. Numerical results for analysis of laminated
composite plates show that the MEM performs quite well

chuyển động của tấm, các ma trận kết cấu cũng được thiết lập
trong hệ tọa độ tương đối chuyển động cùng vận tốc của lực.
Mô hình tấm đặt trên nền đàn hồi và tấm không đặt trên nền
được khảo sát để xét ảnh hưởng của nền. Đồng thời, các thông
số như tỉ số module đàn hồi, tỉ số cạnh/bề dày và các dạng dao
động ảnh hưởng đến tần số dao động tự nhiên của kết cấu tấm
cũng được phân tích. Các kết quả thu được sẽ là tài liệu hữu ích
cho việc nghiên cứu và thiết kế các kết cấu tấm composite
laminate trong thực tiễn.
2. Cơ sở lý thuyết
Xét tấm composite laminate chịu biến dạng uốn bởi các lực
vuông góc với mặt phẳng tấm, hệ trục tọa độ Oxyz được chọn
sao cho mặt phẳng tọa độ Oxy trùng với mặt trung bình

về hệ tọa độ tự nhiên (ξ ,η ) như Hình 5.

Hình 5. Phần tử Q9 trong tọa độ tự nhiên
Các hàm dạng của phần tử Q9 trong hệ tọa độ tự nhiên
được cho bởi (Ghafoori, 2010 [12]):

Ω ⊂ R 2 và trục z vuông góc với mặt phẳng tấm. Theo lý
thuyết biến dạng cắt bậc nhất (Reddy, 1997 [11]), trường
chuyển vị tại một điểm bất kỳ trong tấm được biểu diễn thông
qua trường chuyển vị tại điểm tương ứng trên mặt trung hòa
như sau:
u (=
x , y , z ) uo ( x , y ) + z β x ( x , y )

 h h
=

N 6 = (1 − η 2 ) (ξ + 1) ξ
4
1
N8 = (1 − η 2 ) (ξ − 1) ξ
4

N1 =

(1)

trong đó u , v , w là các thành phần chuyển vị theo phương x, y,
z; β x , β y lần lượt là các góc xoay của pháp tuyến của mặt trung
hòa quanh trục Oy và Ox của hệ tọa độ địa phương với quy
ước chiều dương cho ở Hình 3, Ω là mặt trung hòa của tấm và
h là bề dày của tấm.

N2 =

(2)

Bằng cách sử dụng các hàm dạng, vectơ chuyển vị tại một
điểm bất kì u = u v w β x

dạng và chuyển vị của các nút

T

β y  sẽ được nội suy từ hàm

(1)


0

0

0

N1

0

0

0

...

0

N9

0

0

0

N1

0


0

0

0

N1

...

0

0

0

0

và d là vectơ chuyển vị nút:

và được đánh số từ 1-9 như được thể hiện trên Hình 4.

2


0 

0 



βx9

w09

T

(3)

Biến dạng của tấm bao gồm biến dạng màng, uốn và biến
dạng cắt. Các thành phần biến dạng này được cho bởi các công
thức sau:

ε

ε
ε


x

y

xy

  u
  β
 
 
=


yz

,y



ε + zκ
(4)
=



σ
σ

τ
τ

 τ









xx


0

61

Q62

Q66

0

0

0

0

Q55

Q54

0

0

Q45

Q44

n zk +1


Q16

0

0

21

Q22

Q26

0

0

61

Q62

Q66

0

0

0

0


yy

xy

xz

yz









N w = [ 0 0 N1 0 0 ... 0 0 N 9









∫

∫


12

2

12

)

∂ 2 u ( x, t )
∂ 2u ( r , t )
∂ 2u ( r , t ) ∂ 2u ( r , t )
 = 2
u
=
V2
−
2
V
+
∂t
∂r 2
∂r ∂t
∂t 2

(18)

∂w ( x, t )
∂w ( r , t ) ∂w ( r , t )
w ==

(16)

∂ 2 w ( x, t )
∂2 w ( r, t )
∂2 w ( r, t ) ∂2 w ( r, t )
(20)
 = 2
w
=
V2
− 2V
+
2
∂t
∂r
∂r ∂t
∂t 2
trong đó
T
 = u0 v0 w
0 βx βy 
u = u0 v0 w 0 βx β y  T , u
Tại thời điểm t , miền bài toán trong hệ tọa độ cố định là

T

2

(15)


(7)









(13)

giữa hai trục tọa độ được xác định như sau:

T

Ω′

(12)

Giả sử tải trọng di động theo phương x với vận tốc không
đổi V . Bằng cách sử dụng phương pháp MEM, một hệ tọa độ
( r , s ) gắn liền với tải trọng di động được thiết lập. Mối quan hệ

(6)

∫ δ u q d Ω′

=
+ δ w k f wd Ω′ + δ w c f wd


(k )

z

∫ δ ε Dε d Ω′ + ∫ δ γ D γ dΩ′ + ∫ δ u mud Ω′
p

(11)

n zk +1

1997 [11]).
Phương trình dạng yếu Galerkin cho phân tích ứng xử động
lực học của tấm composite laminate trên nền đàn nhớt có dạng:
p

)

Dmb =
∑

trong đó Qij là hằng số vật liệu biến đổi của lớp thứ k (Reddy,

T

6

k 1=
k 1

Q=
0
66
55
44
16
26
45

xx

=
=
; i , j 4,
5 κ

n k +1
(k ) 2
(i )
1 n 3
( zk +1 − zk3 )Q ij (i, j =
1, 2,6)
Db =
Q
z
dz
=
∑
∑
ij


Q11

∫ Q dz

(5)

y

Q16

(

2

s

(k )
(i )
1 n 2
( zk +1 − zk2 )Q ij (i, j =
1, 2,6)
Q
zdz
=
∑
ij
∫
2k1
k 1=

đưa ra từ định luật Hook:
(k )

D
 D

=
D

T

Ω
Ω
Ω
Ω

  T
T
T
dΩ
+ d  ∫ −2mVN N , r d Ω + ∫ N w c f N w d Ω  + d
=
∫ N mN
Ω
 Ω
Ω

(10)




T

Ω

Ce =

Ω

=
Ke

Ω

∫B

Ω

T
b

Db Bb d Ω + ∫ BTs Ds B s d Ω
Ω

(24)

T
∫ mV N N,rr d Ω − ∫ c f VN N w,r d Ω + ∫ k f N w N wd Ω
2


23
13

T

T
w

Ω

Ω

P=

Ω

∫ N q dΩ
T

(25)

Ω

trong đó

( ),r là đạo hàm bậc nhất theo

r và

( ),rr

∫qN

T

dΩ

(

w = 100 E2 wh / qa
*

3

Góc
hướng sợi

(28)

a/h

Ω

2

(29)

0

0



0.6564

0.6616

0.6693

0.6281

0.6458

20

0.4824

0.4852

0.4913

0.4516

0.4666

100

0.4248

0.4269

0.4337


0.4248

0.4269

0.4337

0.3675

0.4073

Tải phân bố (UDL)

00/900/00

3.1 Phân tích tĩnh
Trong phân tích này, mô hình tấm vuông composite
laminate tựa và không tựa trên nền đàn hồi với điều kiện biên là
tựa đơn 4 cạnh (S-S-S-S) và ngàm 4 cạnh (C-C-C-C) được khảo
sát.
3.1.1

)

Tải hình sin (SSL)

Tấm dao động với tần số ϖ thì phương trình cân bằng trở
thành
0
(K − ω M )u =


0.7572

0.6989

0.7195

100

0.6700

0.6608

0.6697

0.5744

0.6307

10

1.0253

1.0138

1.0250

0.9641

0.9883

(UDL) với q0 = 1N/m tại giữa tấm. Có 2 loại tấm được khảo sát
là tấm 3 lớp (00/900/00), 4 lớp (00/900/900/00). Kết cấu tấm sẽ
được chia thành phần tử có kích thước NxN với N = 16 .

Hình 7. So sánh độ võng không thứ nguyên w giữa các
phương pháp (tải phân bố hình sin SSL)
*

Tải phân bố hình sin (SSL)
Tải phân bố đều (UDL)
Hình 6. Mô hình tấm tựa đơn 4 cạnh

4


Hình 10. Chuyển vị của 3 loại tấm với mức lưới phần tử 20x20

Hình 8. So sánh độ võng không thứ nguyên w giữa các
phương pháp (tải phân bố đều UDL)
*

3.1.2

3.1.3

Tấm vuông chịu tải tập trung trên nền đàn hồi

Trong bài toán này, thông số tấm vẫn được sử dụng tương tự
Ví dụ 3.1.2 và có thêm hệ số nền đàn hồi k f = 107 N/m3 .


ngoài tác dụng. Kết quả này có thể áp dụng trong thực tế thiết
kế và thi công: muốn giảm ứng xử của tấm thì ta cần phải gia cố
nền với độ cứng thích hợp.

Hình 11. Chuyển vị tấm tại điểm đặt lực

Hình 9. Chuyển vị của tấm 5 lớp theo các mức lưới phần tử
Ngoài ra, dựa vào kết quả chuyển vị của từng dạng tấm
trong Hình 10 cho thấy: tấm có cùng bề dày h và cấu tạo bởi
cùng loại vật liệu nhưng tấm có nhiều lớp hơn thì độ võng cũng
như chuyển vị ở tâm tấm cũng nhỏ hơn. Nguyên nhân là khi các
lớp tấm có các hướng sợi khác nhau được ghép đan xen với
nhau thì sự liên kết các phần tử giữa các lớp tấm càng trở nên
chặt chẽ hơn và khăng khít hơn, tạo nên một kết cấu tổng thể
vững chắc và khó bị phá vỡ hơn.

Hình 12. So sánh giữa chuyển vị tại điểm đặt lực của tấm trên
nền đàn hồi và tấm không đặt trên nền đàn hồi
3.2 Phân tích dao động tự do
Trong mục này, việc khảo sát tính chính xác và sự hữu dụng
của phương pháp MEM trong việc phân tích tần số dao động tự
nhiên của kết cấu tấm được thực hiện. Hai điều kiện biên được
áp dụng cho tấm là bốn biên tựa đơn (S-S-S-S) và bốn biên
ngàm (C-C-C-C). Tấm được chia với nhiều hệ lưới khác nhau.
Kết quả của phương pháp MEM sẽ được so sánh với phương
5


pháp khác sử dụng các phần tử như FEM-9 và kết quả được
công bố trong bài báo của tác giả Reddy (1997) [11] và phương

phương pháp. Phương pháp MEM cho lời giải hoàn toàn giống
với FEM-9 và lời giải chính xác của Reddy (1997) [11]. Điều
này chứng tỏ rằng phương pháp MEM tỏ ra hiệu quả trong việc
phân tích tần số dao động tự nhiên của kết cấu tấm.

*

Phương
pháp
MEM9
ES-DSG3
MISQ20
HSDT
p-Ritz
RBF
HOIL
theory
Local
theory
Global
theory

Bảng 3. Bảng so sánh tần số dao động không thứ nguyên

(

ω = ωa / h
*

2


8.2985

9.5674

10.8542

11x11

8.2983

9.5672

10.8541

13x13

8.2982

9.5672

10.8540

15x15

8.2982

9.5671

10.8540


9.5613

10.8490

RBF

8.3101

9.5801

10.8640

Giải tích
(Reddy)

8.2982

9.5671

10.8540

(

động không thứ nguyên ω = ω a / h
2

)

ρ / E2

18.05
18.07
18.05

50
18.68
18.66
18.75
18.46
18.67
18.66

100
18.84
18.82
18.92
18.76
18.84
18.82

10.67

15.07

17.54

18.05

18.67


Hình 13. Sự ảnh hưởng tỉ số E 1 /E 2 trong phương pháp MEM so
với kết quả giải tích của Reddy
Hình 15. Sai số của tần số dao động không thứ nguyên ω các
phương pháp đối với phương pháp p-Ritz
Đồng thời, kết quả trên Hình 14 còn cho thấy nếu tấm càng
mỏng thì dao động tự nhiên của tấm càng lớn và điều này hoàn
toàn phù hợp với thực tế khảo sát.
*

Đồng thời từ Hình 13 cho thấy nếu tỉ số E 1 /E 2 càng lớn thì
sự hội tụ của tần số dao động tự nhiên càng chính xác và khi
lưới được chia càng mịn thì phương pháp MEM càng cho lời
giải tiến đến gần lời giải chính xác đưa ra trong bài báo của
Reddy (1997) [11]. Vì lưới càng mịn thì các phần tử sẽ được
6


3.2.3
Khảo sát các dạng dao động của tấm
Trong bài toán này, tấm composite laminate 3 lớp
(00/900/00) biên ngàm bốn cạnh (C-C-C-C) được khảo sát.
Thông số vật liệu tương tự Ví dụ 3.2.1, tỉ số module E 1 /E 2 =40,
tỉ số a/h =10. Tấm được khảo sát trong bài toán có hệ lưới phần
tử 15x15. Tần số các dạng dao dộng được so sánh với phương
pháp p-Ritz [13] để chứng tỏ tính đúng đắn của MEM.
Bảng 5 thể hiện tần số dao động tự nhiên ứng với 5 mode
dao động đầu tiên của phương pháp MEM và các phương pháp
đã công bố. Kết quả thu được cho thấy kết quả sử dụng phương
pháp MEM rất gần với kết quả của p-Ritz (sai số ở mode 5 là
0.47%). Điều này chứng tỏ tính chính xác của MEM trong phân

moving load on contiuum”, International Journal of
Mechanical Sciences, vol. 300, pp. 126-138, July.2006.
8. Chen YH, Huang YH, “Dynamic stiffness of infinite
Timoshenko beam on viscoelastic foundation in moving
coordinate”, International Journal for Numerical Methods
in Engineering 2000; 48:1-18.
9. Tran Minh Thi, Ang Kok Keng, and Luong Van Hai (2013).
Dynamic analysis of high-speed rail system on twoparameter elastic damped foundation. International
Conference on Advanced Computing and Applications
ACOMP, 23-25/10/2013, Ho Chi Minh City , Vietnam.
10. W.T. Xu, J.H. Lin, Y.H. Zhang, D. Kennedy and F.W.
Williams, “2D moving element method for random
vibration analysis of vehicles on Kirchho plate with Kelvin
foundation”, Latin American Journal of Solids and
Structures 6, pp. 169-183, 2009.
11. Reddy JN, “Mechanics of laminated composite plates –
Theory and Analysis”, New York: CRC Press, 1997.
12. E. Ghafoori, M. Asghari. “Dynamic analysis of laminated
composite plates traversed by a moving mass based on a
first order theory”, Composite structure, vol. 92, pp. 18651867, 2010.
13. H. H. Phan Dao, H. Nguyen Xuan, C. Thai Hoang and T.
Nguyen Thoi, “An Edge-Based Smoothed Finite Element
Method for Analysis of Laminated Composite Plates”,
International Journal of Computational Methods, vol. 10,
No.1, 2013.

Bảng 5. Các tần số dao động tự nhiên của tấm biên ngàm C-CC-C (a/h =10)
Modes
Phương pháp
1

Trong bài báo này việc phân tích ứng xử tĩnh và dao động tự
nhiên của tấm composite laminate sử dụng phương pháp phần
tử chuyển động MEM đã được thực hiện. Thông qua các kết
quả nghiên cứu, một số kết luận quan trọng có thể rút ra như
sau:
 Phương pháp MEM có tính khả thi cao trong việc phân tích
tĩnh và dao động tự do của kết cấu tấm composite laminate.
Các kết quả thu được từ MEM đã được kiểm chứng với
FEM và các kết quả đã công bố.
 Phương pháp MEM cho nghiệm chính xác và hội tụ nhanh
hơn so với phương pháp FEM truyền thống hay phương
pháp CS-DSG3.
7