Xác định tải trọng giới hạn trong bài toán Prandtl bằng phương pháp sai phân hữu hạn thông qua Matlab và FLAC 2D

ĐỊA KỸ THUẬT SỐ 4-2018 11 
XÁC ĐỊNH TẢI TRỌNG GIỚI HẠN TRONG BÀI TOÁN 
PRANDTL BẰNG PHƯƠNG PHÁP SAI PHÂN HỮU HẠN 
THÔNG QUA MATLAB VÀ FLAC2D 
NGÔ THỊ THANH HƢƠNG*, VŨ BÁ THAO** 
Determining ultimate bearing capacity of soils by the finite defference 
method using Matlab and Flac 3D sofwares in the solution Prandtl 
Abstract: The the Prandlt problem has determined the value of the 
soil-bearing capacity of earth foundation which affected by a vertical 
strip load on a half plane. The author proposes a new solution to 
determine the soil-bearing capacity based on the finite difference 
method. This problem was solved by computation program which coded 
by author on Matlab software. Besides, the FLAC 2D, a finite 
difference software, is also used for simulation to solve the Prandtl 
problem. The results shows that the value of the soil-bearing capacity 
which determined by proposed method is compatible with results from 
Prandlt problem and FLAC 2D software. 
1. XÂY DỰNG BÀI TOÁN XÁC ĐỊNH 
TẢI TRỌNG GIỚI HẠN CỦA NỀN ĐẤT* 
Xét bài toán phẳng, phân tố đất chịu tác dụng 
của các ứng suất z , x , xz và trọng lƣợng bản 
thân  nhƣ Hình 1. 
Dƣới tác dụng của tải trọng tĩnh, phân tố 
đất bất kỳ ở trong trạng thái cân bằng, ổn định 
và đủ sức chịu tải nếu nó thỏa mãn các điều 
kiện dƣới đây. 
 Phƣơng trình cân bằng tĩnh học: 
 (1) 0 




zx
zx
x

(2) 0 






xz
xz
z 
trong đó: 
x, z ,xz - các ứng suất có hiệu tƣơng 
ứng với các phƣơng của các trục tọa độ. 
* Khoa Công trình, Tr ng Đ i h c Công nghệ Giao 
thông Vận tải 
** 
Phòng Nghiên cứu Địa kỹ thuật, Viện Thủy Công, 
Viện Khoa h c Thủy lợi Việt Nam 
* Email: vubathao@gmail.com 
z
dz
o x
z
xz

x

z
zx

x
xz
zx
z
 x
x
x
dx
xz
dzz
zx 
 z
z dz








x dx
Hình 1: Phân t đất chịu tác dụng của ứng suất 
 Điều kiện về ứng suất trong đất: đất phải 
luôn luôn chịu nén, hay: 
 0 x và 0 z . (3) 
 Điều kiện cân bằng bền và ổn định: 
Nếu đất nằm trong trạng thái cân bằng bền và 
ổn định thì ứng suất tiếp lớn nhất trong đất max 
phải có giá trị nhỏ nhất (max min), điều kiện 
đó có thể suy ra từ nguyên lý Castiliano (1847-
1884), nó đƣợc phát biểu nhƣ sau: “Trong tất cả 
các tr ng thái cân bằng lực có thể thì tr ng thái 
cân bằng thực xảy ra khi thế năng biến d ng là 
cực tiểu”. Nguyên lý trên có thể đƣợc viết dƣới 
 ĐỊA KỸ THUẬT SỐ 4-2018 12 
dạng công thức nhƣ sau: 3 
 min
1
2
1 2
max dVG
Z
V
 , (4) 
trong đó: Z - thế năng biến dạng do ứng suất 
tiếp  gây ra; 
max- ứng suất tiếp lớn nhất tại điểm đang xét; 
G - môđun trƣợt của đất; max
1

G
- biến dạng trƣợt; 
V - miền lấy tích phân - thể tích khối đất 
đƣợc xét. 
 Điều kiện bền Mohr-Coulomb: 
Ngoài thỏa mãn các phƣơng trình cân bằng, 
điều kiện bền và ổn định nói trên, đất nền là đủ 
sức chịu tải thì các ứng suất phát sinh trong đất 
dƣới tác dụng của tải trọng phải đáp ứng đƣợc 
điều kiện bền f(k) của Mohr-Coulomb: 
0. f(k) ctg  . (5) 
Nhƣ vậy, biểu thức (4) chính là Hàm mục 
tiêu của bài toán xác định tải tr ng giới h n của 
nền đất với các ràng buộc (1), (2), (3) và điều 
kiện bền Mohr-Coulomb (5). Bài toán trên là 
bài toán quy hoạch phi tuyến (do điều kiện bền 
Mohr-Coulomb là phi tuyến), giải đƣợc nó, ta 
xác định đƣợc tải trọng giới hạn của nền đất. 
2. ĐỀ XUẤT PHƯƠNG PHÁP XÁC ĐỊNH 
TẢI TRỌNG GIỚI HẠN TRONG BÀI TOÁN 
PLANDTL DỰA TRÊN PHƯƠNG PHÁP SAI 
PHÂN HỮU HẠN 
2.1. Đề xuất phƣơng pháp 
MÆt tho¸ng ®Êt n»m ngang
(a)
(b)
(c)
p
1 2 i1
i,j
j
m
n
0
x
z
i,j+1
i+1,j
 x
 Z
i,j
i+1,j+1
p
Hình 2. Sơ đồ sai phân 
a - Mô hình kh i đất tính toán; b và c - L ới sai phân và kích th ớc ô l ới 
Sơ đồ tính dƣới dạng đối xứng và sơ đồ lƣới 
sai phân có dạng nhƣ trên hình 2. Số thứ tự nút 
lƣới sai phân theo trục oz thay đổi trong khoảng 
1m và theo trục ox là 1n ( 0n - nút giữa trong 
khoảng n1 ). Điều kiện biên của bài toán là 
trạng thái ứng suất các nút ở cạnh dƣới, cạnh 
trên và hai bên của lƣới, cụ thể nhƣ sau: 
+ Tại mặt trên của khối đất: các nút không 
chịu tác dụng của tải trọng: chỉ có ẩn số là x 
còn z=0 và xz=zx=0. 
+ Các nút chịu tác dụng của tải trọng: chỉ có 
ẩn số là x còn z= p và xz=zx = 0. 
+ Ứng suất tại các nút biên còn lại là 
chƣa biết. 
Nhƣ vậy, để giải đƣợc bài toán thì phải xác 
định sơ đồ sai phân sao cho ứng suất của các nút 
trên biên đều nằm trong các phƣơng trình cân 
bằng và trong hàm mục tiêu. Trên nguyên tắc 
đó, đối với mỗi ô lƣới ta viết phƣơng trình cân 
bằng cho điểm giữa của ô lƣới (điểm nằm giữa 4 
ĐỊA KỸ THUẬT SỐ 4-2018 13 
nút của mỗi ô lƣới) và sơ đồ sai phân đƣợc sử 
dụng trong bài toán là sai phân trung tâm. 
Do các điều kiện biên của bài toán là ứng suất 
ở các nút trên biên và với cách chọn sơ đồ sai 
phân nhƣ trên, ta có thể mở rộng khối đất theo 
phƣơng ngang và theo chiều sâu đến vô hạn. 
Phƣơng trình cân bằng và hàm mục tiêu: 
Xét ô lƣới đƣợc xác định bởi 4 nút: (i,j), 
(i,j+1), (i+1,j), (i+1,j+1). Phƣơng trình cân 
bằng đƣợc viết cho điểm nằm giữa 4 nút của 
ô lƣới đang nên phƣơng trình (1) và (2) 
có dạng: 
 ;0
1
22
1
22
),()1,(),1()1,1(
),(),1()1,()1,1(
zx
x
zx
z
ji
xz
ji
xz
ji
xz
ji
xz
ji
x
ji
x
ji
x
ji
x


 (6) 
.0
1
22
1
22
),(),1()1,()1,1(
),()1,(),1()1,1(
zxzx
x
zx
z
ji
xz
ji
xz
ji
xz
ji
xz
ji
z
ji
z
ji
z
ji
z



 (7) 
Điều kiện đất luôn chịu nén (3) đối với mỗi nút lƣới sẽ là: 
 0),( jix và 0
, jiz . (8) 
Để có hàm mục tiêu (4) dƣới dạng sai phân, ta lƣu ý ứng suất tiếp max và ứng suất pháp  
tƣơng ứng tại điểm đang xét đƣợc xác định đối với bài toán phẳng nhƣ sau: 
2
31
max


 ; 
2
31 
 (9) 
và 
22
2
2
3,1 xz
zxzx 

 
 (10) 
Khi đó, biểu thức của hàm mục tiêu (4) đƣợc biểu diễn theo sai phân có dạng: 
 .min zx
2G
1
.
1 2),(
2
j)(i,
z
j)(i,
x2
max 
  
ji
xz
i jS
dzdx
G


 (11) 
Mô đun trƣợt G chỉ có ở trong hàm mục tiêu (11) (không có trong điều kiện ràng buộc). Về mặt 
toán học G hệ số bình quân gia quyền. Về mặt cơ học, trong tính toán, tác giả xét hai trƣờng hợp là 
G=const. 
Với điều kiện bền Mohr-Coulomb (5), thay (10) vào (9), sau đó vào (5), ta đƣợc: 
 0cossin
2
)(
4
)(
)( 2
2



ckf zxxz
zx (12) 
Nhƣ vậy, theo sai phân hữu hạn, bài toán tìm 
tải trọng giới hạn đối với nền đất là bài toán quy 
hoạch phi tuyến với hàm mục tiêu (11), các ràng 
buộc (6), (7), (8) và điều kiện bền Mohr-
Coulomb (12). Để giải bài toán trên, có thể sử 
dụng các phƣơng pháp giải bài toán quy hoạch 
phi tuyến khác nhau. Ở đây, tác giả sử dụng 
hàm fmincon có sẵn của phần mềm MATLAB 
và chọn phƣơng pháp thử dần. 
2.2. Lập trình xác định tải trọng giới hạn 
bằng MATLAB 
Sử dụng bài toán đặt ra ở trên, xét trƣờng 
 ĐỊA KỸ THUẬT SỐ 4-2018 14 
hợp c≠0; =0, =0, tác giả đã tiến hành giải 
với lƣới sai phân có kích thƣớc m=8, n=15 và 
nhận đƣợc kết quả ứng với các trƣờng hợp tác 
dụng tải trong khác nhau trên nền đất nhƣ 
dƣới đây. 
 Với giá trị áp lực của tải trọng p<4c, ứng 
suất trong tất cả các điểm nút tính toán của 
lƣới sai phân đều đáp ứng điều kiện bất đẳng 
thức f(k)<0 trong biểu thức (12). Đất ở trạng 
thái ổn định. 
 Với giá trị áp lực của tải trọng p =4c, xuất 
hiện chảy dẻo (đạt điều kiện cân bằng trong biểu 
thức (12), f(k)=0) tại hai điểm nút ứng với (i=8, 
j=1) và (i=8, j=2) nhƣ trên hình 3a. Các đƣờng 
đẳng bền f(k) của điều kiện bền Mohr-Coulomb 
nhƣ trên hình 3b. 
p -0.09
-0
.0
9
-0.
09
-0
.0
9
-0
.0
9
-0.0
9 -0.09
-0
.0
9
-0
.0
8
-0
.0
8
-0.
08
-0.08
-0
.0
8
-0
.0
80
2 4 6 8 10 12 14
2
4
6
8
Hình 3a: Các điểm chảy dẻo Hình 3b: Các đ ng đẳng bền f(k) 
 Với giá trị áp lực của tải trọng p =4.3c, xuất 
hiện chảy dẻo tại các điểm nút nhƣ trƣờng hợp 
p =4.0c. 
 Với giá trị áp lực của tải trọng p =4.6c, xuất 
hiện chảy dẻo (đạt điều kiện cân bằng trong biểu 
thức (12), f(k) = 0) tại ba điểm nút ứng với (i=8, 
j=1); (i=8, j=2) và (i=7, j=2) nhƣ trên hình 4a. 
Các đƣờng đẳng bền f(k) của điều kiện bền 
Mohr-Coulomb nhƣ hình 4b. 
p -0.09 -0.0
9
-0
.0
9
-0
.09
-0.
09 -0.0
9
-0
.0
9
-0.08
-0
.0
8
-0.0
8 -0.08
-0
.0
8
-0
.0
8
-0
.07
-0
.0
7
-0.07
-0
.0
7
0 0
2 4 6 8 10 12 14
2
4
6
8
Hình 4a: Các điểm chảy dẻo Hình 4b: Các đ ng đẳng bền f(k) 
 Với giá trị áp lực của tải trọng p =5c, xuất 
hiện chảy dẻo (đạt điều kiện cân bằng trong biểu 
thức (12), f(k) = 0) tại bốn điểm nút ứng với 
(i=8, j=1); (i=8, j=2); (i=7, j=2) và (i=7, j=1) 
nhƣ trên hình 5a. Các đƣờng đẳng bền f(k) của 
điều kiện bền Mohr-Coulomb nhƣ hình 5b. 
ĐỊA KỸ THUẬT SỐ 4-2018 15 
p
-0
.0
9
-0.09
-0.
09 -0
.09
-0
.0
9
-0
.0
8
-0
.0
8
-0.08 -0.0
8
-0.0
8
-0.08
-0
.0
7
-0
.0
7
-0.07 -0
.0
7
-0.0
7
-0
.0
6
-0
.0
6
-0.06
-0
.06
0 0
2 4 6 8 10 12 14
2
4
6
8
Hình 5a: Các điểm chảy dẻo Hình 5b: Các đ ng đẳng bền f(k) 
 Với giá trị áp lực của tải trọng p =5.23c, 
xuất hiện chảy dẻo (đạt điều kiện cân bằng trong 
biểu thức (12), f(k) = 0) tại bốn điểm nút ứng 
với (i=8, j=1); (i=8, j=2); (i=7, j=2) và (i=7, j=1) 
nhƣ trên hình 6a. Các đƣờng đẳng bền f(k) của 
điều kiện bền Mohr-Coulomb nhƣ hình 6b. 
p
-0
.0
9
-0
.0
9
-0.09
-0.09
-0.09
-0
.0
9
-0.09 -0.09
-0.08
-0.08
-0.08
-0.08
-0
.08
-0.08
-0
.0
8
-0.07
-0
.0
7
-0.07
-0
.0
7
-0.07
0 0
2 4 6 8 10 12 14
2
4
6
8
10
Hình 6a: Các điểm chảy dẻo Hình 6b: Các đ ng đẳng bền f(k) 
+ Trƣờng hợp áp lực của tải trọng p=5.24c 
hoặc lớn hơn - bài toán không có nghiệm. 
3. XÁC ĐỊNH TẢI TRỌNG GIỚI HẠN 
TRONG BÀI TOÁN PRANDTL BẰNG 
PHẦN MỀM FLAC 2D 
3.1. Giới thiệu phần mềm Flac 2D 
Phần mềm Flac 2D (hoặc Flac 3D) dựa trên 
phƣơng pháp sai phân hữu hạn, áp dụng cho 
môi trƣờng liên tục nhằm mô phỏng các bài toán 
phẳng (hoặc bài toán không gian) trong môi 
trƣờng đất/đá. Sử dụng các phần tử tiếp xúc 
(interface) để biểu diễn các mặt không liên tục: 
khe nứt, đứt gẫy, phân lớp. Phần mềm có ƣu 
điểm nổi bật trong các bài toán biến dạng lớn, 
ứng xử phi tuyến. Là phần mềm mở, tích hợp 
sẵn 12 mô hình ứng xử vật liệu và các phần tử 
kết cấu: dầm, neo, vỏ  Ngƣời sử dụng có thể 
tích hợp thêm các mô hình vật liệu và các phần 
tử kết cấu khác. Phân tích mô phỏng đƣợc các 
bài toán tĩnh, bài toán động, bài toán nhiệt và từ 
biến 6. 
3.2. Xác định tải trọng giới hạn bằng 
Flac 2D 
Nhằm đánh giá kết quả tính toán bằng 
MATLAP, phần mềm sai phân hữu hạn FLAC 
2D đƣợc sử dụng để mô phỏng bài toán tải trọng 
phá hủy nền Prandtl với các điều kiện biên, tải 
trọng, tính chất cơ lý đất nền tƣơng tự nhƣ phân 
tích bằng MATLAP. 
Sơ đồ tính toán FLAC 2D biểu diễn trên 
 ĐỊA KỸ THUẬT SỐ 4-2018 16 
Hình 7, với điều kiện biên là biên mặt trên tự 
do, biên đáy và biên phải giới hạn chuyển vị cả 
hai phƣơng ngang và đứng, biên trái chỉ giới 
hạn chuyển vị ngang. Tải trọng tác dụng lên 
móng đƣợc mô phỏng thông qua vận tốc của 
một điểm nút tại góc trên phía trái mô hình, 
tƣơng tự nhƣ cách gia tải trong MATLAB. Theo 
hƣớng dẫn của phần mềm Flac 2D, giá trị vận 
tốc lấy là 2.5x10-5m/bƣớc tính toán, nhằm đảm 
bảo đủ nhỏ để không ảnh hƣởng đến kết quả 
tính toán trong nội mô hình. Kích thƣớc lƣới mô 
hình đƣợc chia đúng theo mô hình phân tích 
bằng MATLAB, 8 phần tử theo phƣơng đứng và 
16 phần tử theo phƣơng ngang (Hình 8). 
Hình 7: Sơ đồ tính toán trong FLAC 2D 
Hình 8: Phân chia l ới phần tử 
Kết quả tính toán vùng biến dạng dẻo xuất hiện 
ở các giai đoạn tính toán ứng với các giá trị tải 
trọng p khác nhau đƣợc thể hiện trên Hình 9. Có 
thể thấy rằng, khi giá trị p nhỏ hơn 4c mô hình 
không xuất hiện vùng biến dạng dẻo, khi p có trị 
tƣơng đƣơng với các giá trị p trong tính toán bằng 
MATLAB thì vị trí và số lƣợng phần từ xuất hiện 
biến dạng dẻo là hoàn toàn trùng khớp với kết quả 
tính toán MATLAB. Điều đó khẳng định kết quả 
phân tích tính toán theo phƣơng pháp tác giả đề 
xuất là chính xác. Kết quả tính toán FLAC cũng 
cho thấy sức chịu tải giới hạn của móng p = 5.21c 
(Hình 10), xấp xỉ bằng kết quả tính toán MATLAB 
(p = 5.23c) và của Prandlt (p=5.14c) 1. 
Hình 9a: Tải tr ng tác dụng p = 4.0c Hình 9b: Tải tr ng tác dụng p = 4.3c 
Hình 9c: Tải tr ng tác dụng p = 4.6c Hình 9d: Tải tr ng tác dụng p = 5.3c 
Hình 9: Vùng biến d ng dẻo ở các giai đo n tính toán ứng với tải tr ng p khác nhau 
ĐỊA KỸ THUẬT SỐ 4-2018 17 
Hình 10. Các đ ng quan hệ giữa tải tr ng và chuyển vị 
theo ph ơng đứng t i điểm đặt móng 
Giữa kết quả FLAC 2D và MATLAB có 
sự sai khác về khả năng mô phỏng bài toán ở 
trạng thái phá hủy. Theo đó, sau khi mô hình 
xuất hiện vùng biến dạng dẻo lớn thì 
MATLAB dừng lại, trong khi đó, FLAC 2D 
mô phỏng đƣợc biến dạng lớn nên mặc dù tải 
trọng p đã vƣợt giá trị sức chịu tải giới hạn 
nhƣng phần mềm vẫn tiếp tục phân tích đƣợc 
(Hình 11). Hình 11b và 11c cho thấy, xu 
hƣớng chuyển vị và dịch chuyển của khối phá 
hủy phù hợp với mô hình phá hủy của Prandlt 
(Hình 12). 
Hình 11a: Vùng phá hủy dẻo 
Hình 11c: Hình đẳng h ớng véc tơ chuyển vị 
Hình 11c: Hình đẳng h ớng véc tơ chuyển vị 
Hình 11. Kết quả tính toán FLAC sau khi tải tr ng đ t sức chịu tải p=5.21c 
Hình 12. Mô hình phá hủy của Prandlt đ i với nền đất không ma sát d ới tác dụng 
 tải tr ng móng băng 
 ĐỊA KỸ THUẬT SỐ 4-2018 18 
4. KẾT LUẬN 
Tác giả đƣa ra phƣơng pháp mới xác định tải 
trọng giới hạn dựa trên phƣơng pháp sai phân 
hữu hạn và đƣợc giải thông qua viết chƣơng 
trình trên nền tảng MATLAB. Kết quả tính toán 
đƣợc so sánh với phần mềm sai phân hữu hạn 
FLAC 2D. Một số kết luận đạt đƣợc nhƣ sau: 
- Giá trị tải trọng tác dụng p =5.23c là giá trị 
tải trọng giới hạn tính toán ttghp - khi trong đất 
xuất hiện cơ cấu phá hủy. Tải trọng giới hạn xác 
định đƣợc tt
ghp =5.23c so với tải trọng giới hạn 
theo Prandtl Prandtlghp = 5.14c có sai số là 1,75%. 
- Kết quả tính toán tải trọng giới hạn bằng 
phần mềm FLAC 2D phù hợp với phần mềm 
MATLAB và Prandlt. Mô hình phá hủy nền 
phân tích bằng FLAC 2D phù hợp với 
Prandlt. 
TÀI LIỆU THAM KHẢO 
1. Ngô Thị Thanh Hƣơng và Hồ Sĩ Lành, Cơ 
h c đất, Nhà xuất bản Xây dựng, 2017. 
2. Ngô thị Thanh Hƣơng, Nghiên cứu tính 
toán tải tr ng giới h n của nền đất, Tạp chí Địa 
kỹ thuật, số 2 năm 2011, trang 56-61. 
3. Ngô thị Thanh Hƣơng, Nghiên cứu tính 
toán ứng suất trong nền đất các công trình giao 
thông, Luận án Tiến sĩ, năm 2012. 
4. Bùi Minh Trí (2001), Quy ho ch toán h c, 
Nhà xuất bản khoa học và kỹ thuật. 
5. Jeffery Cooper (2001), A matlab 
companion for multivariable calculus, 
University of Maryland. 
6. ITASCA Consulting Group, H ớng dẫn 
phần mềm sai phân hữu h n FLAC 5.0 - Fast 
Lagrangian Analysis of Continua. 
Ng i phản biện: PGS.TS HOÀNG VĨNH AN