Xác định tự động cơ cấu phá hủy sàn bêtông cốt thép bằng phương pháp đường xoay bất liên tục

48 KỸ THUẬT – CÔNG NGHỆ 
XÁC ĐỊNH TỰ ĐỘNG CƠ CẤU PHÁ HỦY SÀN BÊTÔNG CỐT THÉP 
BẰNG PHƯƠNG PHÁP ĐƯỜNG XOAY BẤT LIÊN TỤC 
NGUYỄN VĂN HIẾU 
Trường Đại học Kiến Trúc Thành phố Hồ Chí Minh - hieu.nguyenvan@uah.edu.vn 
NGUYỄN HUY GIA 
Trường Đại học Công Nghệ Sài Gòn - hgnguyen77@gmail.com 
ĐÀO ĐÌNH NHÂN 
Trường Đại học Kiến Trúc Thành phố Hồ Chí Minh - nhan.daodinh@uah.edu.vn 
(Ngày nhận: 9/9/2016; Ngày nhận lại: 11/11/16; Ngày duyệt đăng: 14/11/2016) 
TÓM TẮT 
Bài báo này giới thiệu một cách tiếp cận tính toán tự động cơ cấu phá hủy cũng như giá trị tối ưu cận trên của 
tải trọng tới hạn cho các loại kết cấu tấm chịu uốn. Phương pháp được sử dụng là phương pháp đường xoay bất liên 
tục dựa trên các cơ cấu chảy dẻo từ các điều kiện tối ưu để dự đoán tải trọng tới hạn. Kết cấu tấm được rời rạc thành 
các phần tử tam giác cứng dẻo tuyệt đối chỉ cho phép chảy dẻo và xoay quanh ba cạnh biên. Các ví dụ số minh họa 
cho phương pháp được đánh giá và so sánh với các lời giải chính xác hay thực nghiệm đã có trước đây. 
Từ khóa: Tải trọng tới hạn; Cơ cấu phá hủy; Phương pháp đường xoay bất liên tục. 
Automated determining the collapsed mechanism of reinforced concrete slabs by yield-
line method 
ABSTRACT 
This paper present an automatic computer program to find the yield-line solution of any given trial finite 
element geometry and the least load required to activate the mechanism of the plate bending. The yield-line method 
is used with the non-linear optimization techniques to predict the ultimate load corresponding to a critical yield-line 
mechanism. In this method the plate is subdivided into triangular mesh and the yield-lines are restricted to occur 
only on the element boundaries. Numerical examples are demonstrated by comparing the present results with 
analytical or experimental solutions available in the literature. 
Keywords: Ultimate load; Collapse mechanism; Yield-line method. 
1. Giới thiệu 
Lời giải chính xác trong các bài toán tấm 
khi kết cấu đạt tới trạng thái tới hạn rất phức 
tạp và chỉ có thể giải quyết được một số bài 
toán đơn giản bằng thủ công. Đây là một bài 
toán phân tích giới hạn của kết cấu đòi hỏi 
nhiều phép tính lặp chính xác. Ngoài ra, việc 
nghiên cứu các kết cấu bằng thực nghiệm rất 
khó khăn và tốn kém. Trong nhiều trường hợp 
việc tìm lời giải chính xác không thể tiến hành 
được, đặc biệt là đối với các kết cấu phức tạp 
như tấm/vỏ. Do đó, việc nghiên cứu kết cấu 
qua mô hình không chỉ đem lại hiệu quả kinh 
tế mà còn có ý nghĩa khoa học rất lớn. 
Các kết cấu tấm được sử dụng rất nhiều 
trong các ngành xây dựng, cơ khí, đóng tàu, 
hàng không,... Độ tin cậy và tuổi thọ của kết 
cấu tấm phụ thuộc chặt chẽ vào tính chất và 
cường độ của ngoại lực, vào vật liệu và sự 
chính xác của sơ đồ tính. Nếu sơ đồ tính càng 
chính xác thì việc tính toán càng phức tạp. 
Trong thiết kế kết cấu tấm, mục đích chính 
yếu của người kỹ sư là phải đảm bảo cho kết 
cấu có một hệ số an toàn thích hợp để kết cấu 
làm việc bình thường và không bị phá hoại 
dưới tải trọng thiết kế. Vì vậy, việc dự đoán 
tải trọng giới hạn mà kết cấu có khả năng chịu 
được cũng như cơ cấu phá hủy của kết cấu ở 
trạng thái tới hạn là cực kỳ quan trọng và cần 
thiết. Nó giúp cho người kỹ sư dự đoán được 
 TẠP CHÍ KHOA HỌC ĐẠI HỌC MỞ TP.HCM – SỐ 51 (6) 2016 49 
ứng xử của kết cấu, dự báo sự hình thành và 
phát triển vết nứt trong kết cấu cũng như đánh 
giá tuổi thọ của công trình. Sự gia tăng tải 
trọng ngoài vùng giới hạn đàn hồi của tấm 
dẫn đến sự hình thành các đường chảy dẻo 
(yield-line) và khi chúng phát triển lan tỏa 
hình thành cơ cấu chảy dẻo (yield-line 
mechanism) thì kết cấu tấm sẽ sụp đổ. Tải 
trọng tại thời điểm tấm bị sụp đổ được gọi là 
tải trọng tới hạn (critical collapsed load). Việc 
xác định chính xác tải trọng tới hạn này đóng 
vai trò rất quan trọng trong việc phân tích giới 
hạn của tấm. 
Phương pháp đường xoay bất liên tục 
trong phân tích giới hạn tấm dựa trên các cơ 
cấu chảy dẻo cho trước để dự đoán tải trọng 
tới hạn. Phương pháp này được đưa ra đầu 
tiên bởi Ingerslev (1923) và được tiếp tục phát 
triển bởi Johansen (1962), Wood (1961) và 
các nhà nghiên cứu khác như Mansfield 
(1957), Morley (1965) hay Johnson (1994, 
1995). Phương pháp này đưa ra trước một cơ 
cấu tương thích chuyển vị (hoặc vận tốc) bao 
gồm các miền tuyệt đối cứng giao nhau tại các 
đường chảy dẻo mà tại đó có xuất hiện sự 
xoay tương đối lẫn nhau. Ứng xử của vật liệu 
khi đó được xem như cứng-dẻo tuyệt đối. Từ 
đó việc xác định giá trị tới hạn của tải trọng 
được tiến hành dựa vào lý thuyết cận trên của 
phân tích giới hạn thông qua việc cân bằng 
năng lượng tiêu tán nội tại các đường chảy 
dẻo với năng lượng tiêu tán ngoại do tải trọng 
gây ra sự biến dạng của tấm theo cơ cấu phá 
hủy cho trước. Tuy phương pháp này đơn giản 
và hiệu quả nhưng việc áp dụng nó trong thực 
tế gặp nhiều khó khăn do sự phức tạp của các 
cơ cấu chảy dẻo cũng như việc xác định đâu 
là cơ cấu gãy đổ nguy hiểm nhất. Phương 
pháp này chỉ thích hợp cho việc tính toán thủ 
công một số bài toán đơn giản. Một số 
phương pháp tính mới gần đây như sử dụng 
phần tử xoay tự do kết hợp cận dưới của 
Salam Al-Sabah et al. (2013) hay phương 
pháp tối ưu lớp bất liên tục của Gilbert et al. 
(2014) rất hiệu quả trong việc xác định cơ cấu 
chảy dẻo của tấm nhưng số lượng nút hay 
phần tử rất lớn dẫn đến khối lượng tính toán 
khá lớn. 
Vì vậy trong bài báo này nhóm tác giả sẽ 
giới thiệu và áp dụng cải tiến một loại phần tử 
đặc biệt được đề nghị bởi Munro và Da 
Fonseca (1978) để có thể tự động hóa việc 
tính toán giá trị tải trọng tới hạn cũng như xác 
định cơ cấu phá hủy nguy hiểm của các kết 
cấu tấm chịu uốn một cách tự động thông qua 
công cụ máy tính với một số lượng rất ít phần 
tử mô phỏng. Việc tính toán này bao gồm: (1) 
rời rạc kết cấu khảo sát thành những phần tử 
tam giác Munro-Da Fonseca tuyệt đối cứng 
chỉ cho phép chảy dẻo, (2) xoay quanh ba 
cạnh biên và (3) kết hợp các điều kiện tối ưu 
để tìm lời giải tốt nhất. 
2. Giới thiệu sơ lược về phần tử 
Munro-Da Fonseca 
2.1. Quan hệ đối ngẫu tĩnh học và động học 
Theo phương pháp rời rạc hóa của 
Munro-Da Fonseca thì kết cấu tấm được rời 
rạc thành các phần tử tam giác kết nối nhau 
tại ba cạnh và ba điểm nút. Biến dạng ngang 
của tấm được diễn tả thông qua vector 
chuyển vị ngang (w) của các nút. Phần tử 
tam giác được giả định là tấm phẳng tuyệt 
đối cứng chỉ cho phép chảy dẻo và xoay 
quanh ba cạnh biên và các góc xoay dọc theo 
các cạnh phải là hằng số suốt chiều dài của 
cạnh tấm (xem Hình 1). 
+
1
3
+
+
f1
f2
f3
h2
h3h1
l2
b2
a2
Edge 21(x1,y1)
2
(x2,y2)
a3
l3
b3
Edge 3 Ed
ge
 1
b1
a1
l1
(x3,y3)
3
m3
m2
m1
2 (x2,y2)
(x1,y1)
(x3,y3)
 2
2
3
2
2
2
1
1
2l
lll
a
 1
23212321
1
))(())((
l
yyxxxxyy
h
 (a) (b) 
Hình 1. (a) Chi tiết hình học phần tử Munro-Da Fonseca; (b) Các biến nội lực của phần tử 
50 KỸ THUẬT – CÔNG NGHỆ 
Các góc xoay quanh cạnh tấm được lưu 
trong vector . Với các kích thước hình học 
như trong Hình 1, phương trình động học liên 
hệ giữa góc xoay i của cạnh thứ i của một 
phần tử tam giác với các thành phần chuyển vị 
của nó là 
1 1
1 1 1 1 1
1 1
2 2
2 2
2 2 2 2 2
33
3 3
3 3 3 3 3
1
1
1
e
e
e
b a
h l h l h
w
a b
w
l h h l h
w
b a
l h l h h



trong đó i được ký hiệu là góc xoay 
quanh cạnh của một phần tử đang xét. Phương 
trình động học liên hệ giữa vectơ góc xoay 
cạnh của toàn hệ có thể viết dưới dạng ma 
trận như sau: 
 = Ew (2) 
với E được định nghĩa là ma trận biến đổi 
động học 
Để phân biệt giữa góc xoay dương 
(sagging) và góc xoay âm (hogging), ta sử 
dụng hai vector không âm +, - sao cho điều 
kiện sau đây được thỏa mãn: 
  = + -- (3) 
Công thức đối ngẫu với biến số là 
moment chảy dẻo trên cạnh (m) và lực nút (f) 
có thể được thiết lập như sau: nếu trên Hình 
1b, moment trên một đơn vị dài được giả định 
là hằng số lần lượt là m1, m2, m3 dọc theo các 
cạnh của tam giác và lực nút tương ứng là f1
e
, 
f2
e
, f3
e
 thì điều kiện cân bằng cho phần tử tam 
giác là: 
3
2
1
322
2
11
1
33
3
211
1
33
3
22
2
1
3
2
1
1
1
1
m
m
m
hhl
b
hl
a
hl
a
hhl
b
hl
b
hl
a
h
f
f
f
e
e
e
Bằng cách lấy tổng các phân phối trên 
toàn miền ta sẽ thu được một hệ đầy đủ các 
điều kiện cân bằng như sau: 
 f = E
T
m (5) 
2.2. Quan hệ ứng xử 
Với điều kiện đồng nhất, tiêu chuẩn chảy 
dẻo cho tất cả moment uốn của các cạnh phần tử 
sẽ được giới hạn bởi moment kháng chảy dẻo 
dương m+ hay âm và m- tương ứng như sau 
+
*
-
*
mI
m
-I m
trong đó I là ma trận đơn vị; 
Sử dụng quan niệm thế năng chảy dẻo thu 
được 
T
* *π = N m-m 0 (7) 
 
+
+ -
*-
θ
θ = θ - θ = I -I = Nθ
θ
 (8) 
2.3. Công thức quy hoạch tuyến tính 
Các tải áp đặt tại nút được định nghĩa bởi 
vector f0 và công ngoại thực hiện trên chuyển 
vị ảo của cơ cấu được ràng buộc là một đơn 
vị: T
0f w 1 . Vì vậy các phương trình ràng 
buộc có thể được viết như sau 
+
+T
-o
-
θ
1 θ0 f
θ 0
0N -E θ
w
,
 (9) 
Do đó sẽ có ne+1 ràng buộc đẳng thức, 
với ne là số cạnh trong hệ lưới mà có thể chịu 
moment uốn. Hàm mục tiêu biểu diễn năng 
lượng tiêu tán trên các đường chảy dẻo là hàm 
cần phải cực tiểu hoá như sau 
Cực tiểu hàm số: + + - -* * * *z = m θ +m θ (10) 
trong đó + -* *m ,m là các vector mô tả các 
giá trị moment dẻo xuất hiện khi có phát sinh 
các góc xoay dương hoặc âm tương ứng. 
Bài toán có thể viết dưới dạng đối ngẫu 
như sau 
Cực đại hàm số  
λ
1
m
y .
 (1(1) 
với các điều kiện ràng buộc:
T +
*
T -
o
0 N λ m m
0
m 0f -E m
,
 (1(2) 
3. Công thức phi tuyến với các biến 
hình học là đường chảy dẻo 
Do phương pháp quy hoạch tuyến tính 
dựa trên phần tử Munro-Da Fonseca không 
xác định tự động được cơ cấu chảy dẻo tối ưu 
thật sự nếu không biết trước cơ cấu chảy dẻo 
cho trước dọc theo các cạnh của hệ lưới phần 
tử chọn trước nên nhu cầu khảo sát tìm kiếm 
một phương pháp tối ưu có thể hiệu chỉnh vị 
trí các nút trong hệ lưới phần tử để cực tiểu 
(1) 
(4) 
(6) 
 TẠP CHÍ KHOA HỌC ĐẠI HỌC MỞ TP.HCM – SỐ 51 (6) 2016 51 
hóa tải trọng tới hạn là thật sự cần thiết. Vì thế 
bài toán xác định cơ cấu và tải trọng tới hạn 
với ẩn số tọa độ nút có thể viết dưới dạng bài 
toán phân tích tối ưu phi tuyến như sau: 
Cực tiểu hàm số L() ;  Rq 
với điều kiện ràng buộc Ki() 0, 
i=1,2,...,nk 
(13) 
trong đó  là vector tọa độ các nút của hệ 
lưới. L() = z là hệ số tải trọng tới hạn định 
nghĩa ở phương trình (2) thông qua cực tiểu 
giá trị góc xoay +, - và độ võng w với giá trị 
 cho trước. q là số biến vị trí nút. Điều kiện 
ràng buộc (13) đảm bảo không có nút nào 
nằm giữa cạnh, không phần tử nào mất đi và 
không có phần tử mới nào tạo ra trong quá 
trình tối ưu hóa. Chi tiết về giải thuật tối ưu 
hóa có thể được tìm thấy trong các tài liệu 
tham khảo (Jennings (1996); Gill et al. 
(1981); McKeown et al. (1990); Thavalingam 
et al. (1998)). 
4. Ví dụ số 
Các ví dụ số trong phần này nhằm mô tả 
tính đơn giản và hiệu quả của phương pháp 
tính toán nêu trên trong việc tìm kiếm tự động 
cơ cấu gãy đổ và tải trọng giới hạn tương ứng 
gây ra sự sụp đổ đó. Các ví dụ được chọn để 
thuận lợi cho việc tính toán nhưng chúng cũng 
mô tả đầy đủ các điều kiện biên như cạnh 
ngàm, tựa đơn hay tự do và các dạng tải trọng 
như tải tập trung hay phân bố đều. Trong các 
tính toán thì giá trị moment chảy dẻo trên một 
đơn vị chiều dài được ký hiệu là m và giá trị 
tải trọng tới hạn dự đoán bằng tiến trình tối ưu 
hóa được đưa dưới dạng hệ số cơ cấu L()/m. 
Tỷ lệ moment kháng dẻo cốt thép 2 phương là 
=2 được chọn để mô tả sự trực hướng. 
4.1. Sàn bêtông cốt thép chữ nhật tựa 
đơn ở ba cạnh đặt thép trực hướng 
Xét một ô sàn bêtông cốt thép tựa đơn ở 
ba cạnh, chịu lực phân bố đều p và sàn được 
đặt thép trực hướng với các giá trị moment 
kháng dẻo của cốt thép 2 phương là mp, mp. 
Hai cơ cấu chảy dẻo xác định theo phương 
pháp giải tích cân bằng được vẽ như trên Hình 
2 phụ thuộc vào tỷ lệ hai cạnh của bản sàn. 
pm
m p
y
a-
2y
y
a
b
x
pm
mp
a
a/
2
a/
2
b
68233.0 
a
b

 3
222
2
b
a
b
aa
xoptmal
68233.0 
a
b
2
2
2
2
3
22
246



b
a
b
a
a
m
x
m
p
p
optimal
p
 

2
3
2
3
2
a
b
a
bb
y
 (a) (b) 
Hình 2. Cơ cấu gãy đổ và nghiệm giải tích với tỷ lệ hai cạnh 
 (a) lớn hơn 0.68233 và (b) nhỏ hơn 0.68233 
52 KỸ THUẬT – CÔNG NGHỆ 
Để khảo sát khả năng tự động tìm cơ cấu 
gãy đổ và tải trọng tới hạn bằng chương trình 
mô phỏng số, ta sẽ xét 2 trường hợp cụ thể 
như sau: 
Trường hợp (a): a=8; b=10; mp=100; =2: 
Hình 3. mô tả lưới ban đầu với các ký 
hiệu của cạnh, nút và phần tử được đánh nhãn 
lần lượt là L, N và T. Các biến số vị trí nút  
trong trường hợp này gồm có 1 là hoành độ x 
của nút N5; 2 là hoành độ x của nút N6; 3 là 
tung độ y của nút N6; 4 là tung độ y của nút 
N7; 5 là hoành độ x của nút N8; 6 là tung độ 
y của nút N9. 
Hình 3. Lưới khởi tạo 
Hình 4. mô tả kết quả lưới có hai phần tử 
T1 và T2 có xu hướng triệt tiêu và do đó các 
đường L2, L4 và L5 sẽ có xu hướng trùng 
nhau thành đường thẳng. Tổng góc xoay 4 và 
5 hầu như xấp xỉ bằng góc xoay 11. Điều này 
đã giúp cho quá trình tối ưu tự động dự đoán 
được sẽ có một đường chảy dẻo chảy thẳng 
vào góc của ô bản sàn đang khảo sát. Kết quả 
phân tích số với tối ưu hóa tự động được so 
sánh với lời giải chính xác ở Bảng 1. 
(a) Lưới khởi tạo (b) Lưới tối ưu 
Hình 4. Tiến trình tối ưu hóa tìm đường gãy đổ 
0 1 2 3 4 5 6 7 8 9 10
-2
-1
0
1
2
3
4
5
6
N1
N2
N3
N4N5
N6
N7
N8
N9
T1
T2 T3
T4
T5
T6
T7
T8
L1
L2L3 L4
L5
L6
L7 L8
L9L10
L11
L12
L13
L14
L15
L16
Initial Mesh: Nnode = 9; Nelem = 8; Nline =16; Ngeo = 6
Geometric variable  = [1; 2; 3; 4; 5; 6] = [x5; x6; y6; y7; x8; y9]
0 2 4 6 8 10
-2
-1
0
1
2
3
4
5
6
N1
N2 N3
N4N5
N6
N7
N8
N9
T1
T2 T3
T4
T5
T6
T7
T8
L2
L4
L5
L11
L15
First Iteration:  = [x5;x6;y6;y7;x8;y9] = [6;1.5;1;1.5;1;1.5]; L() = 29.2063
Edges Rotations: 2=-0.0162; 4=0.015; 5=0.0094; 11=0.0189; 15=0.0157
0 1 2 3 4 5 6 7 8 9 10
-2
-1
0
1
2
3
4
5
6
N1
N2 N3
N4N5
N6
N7
N8
N9
T3
T4
T5
T6T7
T8
L2
L4
L5
L11
L15
Optimum: 
 = [x5;x6;y6;y7;x8;y9]= [7.1123;5.9968;3.3726;0.0328;0.046671;3.7912];L()= 23.7552
Edge Rotations: 2=-0.0089; 4=0.0104; 5=0.0084; 11=0.0188; 15=0.0164
 TẠP CHÍ KHOA HỌC ĐẠI HỌC MỞ TP.HCM – SỐ 51 (6) 2016 53 
Bảng 1 
Tọa độ lưới tối ưu và giá trị tải tới hạn sàn tựa đơn với tỷ lệ cạnh lớn hơn 0.68233 
Lời giải 1 2 3 4 5 6 L () 
Phân tích số 7.1123 - - 0.0328 0.0467 - 23.7552 
Giải tích 7.1073 - - 0 0 - 23.7561 
Trường hợp (b): a=20; b=10; mp=100; =2: 
Hình 5. mô tả lưới khởi tạo ban đầu với 
các biến số vị trí nút  trong trường hợp này 
như sau 1 là hoành độ x của nút N5; 2 là 
hoành độ x của nút N6; 3 là tung độ y của nút 
N6; 4 là tung độ y của nút N7; 5 là hoành độ 
x của nút N8; 6 là tung độ y của nút N9. 
Hình 5. Lưới khởi tạo 
Hình 6. mô tả lưới và các đường chảy dẻo 
xuất hiện trong tiến trình tối ưu của thời điểm 
khởi tạo ban đầu và thời điểm kết thúc tìm ra 
cơ cấu gãy đổ và giá trị tải trọng tới hạn tương 
ứng cho ô sàn của trường hợp (b). 
(a) Lưới khởi tạo (b) Lưới tối ưu 
Hình 6. Tiến trình tối ưu hóa tìm đường gãy đổ 
-2 0 2 4 6 8 10 12
0
2
4
6
8
10
N1
N2 N3
N4N5
N6
N7
N8
N9
T1
T2
T3
T4
T5
T6
T7
T8
L1
L2L3 L4
L5
L6
L7
L8
L9L10
L11
L12
L13L14
L15
L16
Initial Mesh: Nnode = 9; Nelem = 8; Nline =16; Ngeo = 6
Geometric variable  = [1; 2; 3; 4; 5; 6] = [x5; x6; y6; y7; x8; y9]
-2 0 2 4 6 8 10 12
0
2
4
6
8
10
N1
N2 N3
N4N5
N6
N7
N8
N9
T1
T2
T3 T4
T5
T6
T7
T8
L2 L4
L5
L7
L16
First Iteration:  = [x5;x6;y6;y7;x8;y9] = [6;1.5;1;1.5;1;4]; L() = 8.5425
Edges Rotations: 2=-0.0036; 4=0.0018; 5=0.0021; 7=0.0033; 16=0.003
-2 0 2 4 6 8 10 12
0
2
4
6
8
10
N1
N2 N3
N4N5
N6
N7
N8
N9
T4
T5
T6
T7
T8
L2
L5
L7
L16
Optimum: 
= [x5;x6;y6;y7;x8;y9]= [9.6241;3.3073;1.8764;0.032627;0.05262;5.6082]; L()= 7.145
Edges Rotations: 2=-0.0025; 4=0; 5=0.0024; 7=0.0026; 16=0.0024
54 KỸ THUẬT – CÔNG NGHỆ 
x x
a
/2
a
/2
mp
pm
b
b-2x
a
Tương tự trường hợp (a), kết quả tối ưu 
hình trên hình 4.8.b cũng cho thấy phần tử T1, 
T2 và T3 có xu hướng triệt tiêu và các đường 
L2, L4, L5, L7 và L16 sẽ trở thành đường 
thẳng. Các giá trị góc xoay 5, 7 và 16 hầu 
như bằng nhau và quá trình tối ưu cho kết quả 
có đường chảy dẻo chảy thẳng vào góc của ô 
bản sàn đang khảo sát. 
Kết quả phân tích số với tối ưu hóa tự 
động cũng được so sánh với lời giải chính xác 
cho ở Bảng 2 bên dưới đây 
Bảng 2 
Tọa độ lưới tối ưu và giá trị tải tới hạn sàn tựa đơn với tỷ lệ cạnh nhỏ hơn 0.68233 
Lời giải 1 2 3 4 5 6 L () 
Phân tích số - - 0.0326 0.0526 5.6082 7.1450 
Giải tích - - - 0 0 5.5982 7.1452 
4.2. Sàn bêtông cốt thép chữ nhật ngàm 
bốn cạnh đặt thép trực hướng 
Xét một ô sàn bêtông cốt thép tựa ngàm ở 
bốn cạnh, chịu lực phân bố đều p và sàn được 
đặt thép trực hướng với các giá trị moment 
kháng dẻo của cốt thép 2 phương là mp, mp. 
Cơ cấu chảy dẻo với các giá trị tham biến xác 
định theo phương pháp giải tích cân bằng 
được vẽ trên Hình 7. 
Giá trị tới hạn của biến định vị x: 

 3
2
2
b
a
b
aa
x 
Giá trị tải trọng tới hạn tương ứng: 
2
)1(6
x
im
p
p 
2
2
2 3
)1(24


b
a
b
a
a
im
p
p 
(a) (b) 
Hình 7. Tấm chữ nhật ngàm chịu tải phân bố đều: (a) nghiệm giải tích; (b) cơ cấu gãy đổ 
Để mô phỏng số bài toán nêu trên, ta giả 
sử chọn a=10; b=20; mp=100; =2; i=1. 
Do tính đối xứng nên chỉ một nửa tấm sàn 
được mô hình và phân tích. Hình 7 mô tả lưới 
ban đầu với các ký hiệu của cạnh, nút và phần 
tử được đánh nhãn lần lượt là L, N và T. Các 
biến số vị trí nút  trong trường hợp này như 
sau: 1 là hoành độ x của nút N5; 2 là tung độ 
y của nút N5; 3 là hoành độ x của nút N6. 
Hình 8. Lưới khởi tạo 
-2 0 2 4 6 8 10 12
0
2
4
6
8
10
N1
N2 N3
N4
N5
N6
T1
T2
T3
T4
T5
L1
L2
L3
L4
L5
L6
L7
L8
L9
L10
Initial Mesh: Nnode = 6; Nelem = 5; Nline =10; Ngeo = 3
Geometric variable  = [
1
; 
2
; 
3
] = [x
5
; y
5
; y
6
]
 TẠP CHÍ KHOA HỌC ĐẠI HỌC MỞ TP.HCM – SỐ 51 (6) 2016 55 
Hình 9. mô tả lưới và các đường chảy 
dẻo xuất hiện trong tiến trình tối ưu của 
thời điểm khởi tạo ban đầu và thời điểm kết 
thúc tìm ra cơ cấu gãy đổ và giá trị tải 
trọng tới hạn tương ứng cho ô sàn là L()= 
35.4440. 
(a) Lưới khởi tạo (b) Lưới tối ưu 
Hình 9. Tiến trình tối ưu hóa tìm đường gãy đổ 
Cơ cấu gãy đổ cùng với giá trị tải trọng gây sụp đổ L()= 35.4440 này được so sánh với lời 
giải chính xác cho ở Hình 2a. Cho ở Bảng 3 sau đây: 
Bảng 3 
Tọa độ lưới tối ưu và giá trị tải trọng tới hạn của sàn chữ nhật ngàm bốn cạnh 
Lời giải 1 2 3 L () 
Phân tích số 8.2295 5.0000 5.0000 35.4440 
Giải tích 8.2295 5.0000 5.0000 35.4440 
4.3. Sàn bêtông cốt thép tam giác tựa 
đơn trên hai cạnh 
Xét 1 ô sàn bêtông cốt thép hình tam giác 
tựa đơn ở 2 cạnh, đặt thép đẳng hướng chịu tải 
tập trung tại trung điểm cạnh tự do như Hình 10. 
Lời giải chính xác cho hệ số tải trọng tới hạn 
o
p
c
P
m2
Hình 10. Tấm tam giác chịu tải tập trung 
Mô phỏng số ô sàn ở hình 4.9. với các 
thông số giả sử như sau: a=10; mp=100; Po=10. 
Hình 11. mô tả lưới khởi tạo ban đầu với 
các biến số vị trí nút  trong trường hợp này 
như sau 1 là hoành độ x của nút N5; 2 là 
hoành độ x của nút N6; 3 là tung độ y của nút 
N6; 4 là tung độ y của nút N7. 
Hình 11. Lưới khởi tạo 
-2 0 2 4 6 8 10 12
0
2
4
6
8
10
N1
N2 N3
N4
N5
N6
T1
T2
T3
T4
T5
L1
L2
L3
L4
L5
L7
L10
First Iteration:  = [x
5
;y
5
;y
6
] = [3;3.5;3]; P
critical
 = 51.0774
Edge Rotations: 
2
=0.01; 
3
=0.008; 
7
=0.011; 
8
=0
-2 0 2 4 6 8 10 12
0
2
4
6
8
10
N1
N2 N3
N4
N5 N6
T1
T2
T3
T4
T5
L1
L2
L3
L4
L7
L10
Optimum:  = [x
5
;y
5
;y
6
] = [8.2295;5;5]; P
critical
 = 35.444
Edge Rotations: 
2
=0.006; 
3
=0.006; 
7
=0.011
-2 0 2 4 6 8 10 12
0
2
4
6
8
10
N1
N2 N3
N4
N5
N6
N7
T1
T2
T3
T4
T5
T6
L1
L2
L3
L4
L5
L6 L7
L8
L9
L10
L11
L12
Initial Mesh: Nnode = 7; Nelem = 6; Nline =12; Ngeo = 4
Geometric variable  = [1; 2; 3; 4] = [x5; y5; y6; x7]
56 KỸ THUẬT – CÔNG NGHỆ 
Tiến trình tối ưu trên Hình 12b cho thấy 
kết quả phần tử T2 vàT3 có xu hướng triệt tiêu 
và các đường L2, L6, L7 và L11 có xu hướng 
thành thẳng hàng. Các giá trị góc xoay 5, 7 và 
16 cũng xấp xỉ như nhau. Điều này đã giúp cho 
quá trình tối ưu tự động dự đoán được đường 
chảy dẻo chảy thẳng vào góc của ô bản sàn 
đang khảo sát. Giá trị hệ số tải trọng tới hạn t m 
được là L()= 20.0005 so sánh với giá trị giải 
tích là 20.0000 đạt độ chính xác khá cao. 
(a) Lưới khởi tạo (b) Lưới tối ưu 
Hình 12. Tiến trình tối ưu hóa tìm đường gãy đổ 
5. Kết luận 
Bài báo này đã minh họa và đưa ra một 
chương trình tính toán tự động cơ cấu phá 
hủy cũng như giá trị tối ưu cận trên của tải 
trọng tới hạn cho các loại kết cấu tấm chịu 
uốn có các điều kiện biên khác nhau và tải 
trọng tác dụng bất kỳ dựa trên việc tối ưu 
hóa hệ lưới dùng phần tử Munro-Da 
Fonseca trong phân tích chảy dẻo. Các kết 
quả nghiên cứu số cũng đã được so sánh và 
kiểm chứng qua các lời giải giải tích cho 
thấy tính hiệu quả và độ chính xác tin cậy 
cao của phương pháp. Tuy nhiên kỹ thuật tối 
ưu đề cập trong bài báo này còn một số hạn 
chế nhất định trong việc đạt được sự hội tụ 
ổn định bởi sự hiện diện của sự bất liên tục 
về độ dốc trong hàm tối ưu. Điều này sẽ 
được tiếp tục nghiên cứu và cải thiện trong 
thời gian sau này để có thể khảo sát thêm 
nhiều tham số ảnh hưởng đến tải trọng tới 
hạn như bề dày tấm, các cách đặt lưới cốt 
thép theo phương bất kỳ 
Tài liệu tham khảo 
Ingerslev A (1923). The strength of rectangular slabs. The Structural Engineer, 1, 3-14. 
Johansen KW (1962). Yield line theory. London: Cement and Concrete Association. 
Wood RH (1961). Plastic and elastic design of slabs and plates. London: Thames & Hudson. 
Jones LL (1962). Ultimate load analysis of reinforced and prestressed concrete structures. London: Chatto and 
Windus. 
Mansfield EH (1957). Studies in collapse analysis of rigid plastic plates with a square yield diagram. Proceeding 
Royal Society, 241, 311-338. 
Morley CT (1965). Equilibrium methods for exact upper bounds of rigid plastic plates. In: Recent developments in 
yield line theory. London: Cement and Concrete Association, MCR Special Publication. 
-2 0 2 4 6 8 10 12
0
2
4
6
8
10
N1
N2 N3
N4
N5N6
N7
T1
T2
T3
T4
T5
T6
L2
L6
L7
L9
L11
First Iteration:  = [x4;x5;y5;y6;x7] = [5;2.5;3;3.5;2]; L() = 34.7455
Edges Rotations: 2=0.0132; 6=-0.026; 7=0.0229; 9=-0.006; 11=0.0285
-2 0 2 4 6 8 10 12
0
2
4
6
8
10
N1
N2 N3
N4
N5
N6
N7
T1
T3 T4
T5
T6
L2
L6
L7
L11
Optimum:  = [x5;y5;y6;x7] = [4.7972;4.7972;0.032595;0.03395]; L() = 20.0005
Edges Rotations: 2=0.0144; 6=-0.0142; 7=0.0139;11=0.0283
 TẠP CHÍ KHOA HỌC ĐẠI HỌC MỞ TP.HCM – SỐ 51 (6) 2016 57 
Johnson D (1994). Mechanism determination by automated yield line analysis. The Structural Engineer, 72 (19/4), 
323-327. 
Johnson D (1995). Yield-line analysis by sequential linear programming. International Journal Solids Structures, 
32, 1395-1404. 
Salam Al-Sabah, Abd; Falter, Holger (2013). Finite element lower bound "yield line" analysis of isotropic slabs 
using rotation-free elements. Engineering Structures, 53, 38-51. 
Gilbert, M.,He, L., Smith, C.C. & Le, C. (2014). Automatic Yield-Line Analysis of Slabs Using Discontinuity 
Layout Optimization. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 
470 (2168). 
Munro J, Da Fonseca AMA (1978). Yield line method by finite elements and linear programming. The Structural 
Engineer, 56 (2), 37-44. 
 Jennings A (1996). On the identification of yield-line collapse mechanisms. Engineering Structures, 18(4), 
332-337. 
 Gill PE, Murray W, Wright MH (1981). Practical optimisation. New York: Academic Press. 
McKeown JJ, Meegan D, Sprevak D (1990). An introduction to unconstrained optimisation. Bristol: Adam Hilger. 
Thavalingam, A., Jennings, A., McKeown, J.J., and Sloan.D (1998). A computerised method for rigid-plastic yield-
line analysis of slabs. Computers & Structures, 68(6), 601-612.