Áp dụng phương pháp không lưới cho tính toán cọc đơn trong môi trường đất đàn hồi ba chiều

 ĐỊA KỸ THUẬT SỐ 3-2015 46 
ÁP DỤNG PHƯƠNG PHÁP KHÔNG LƯỚI CHO TÍNH TOÁN 
CỌC ĐƠN TRONG MÔI TRƯỜNG ĐẤT ĐÀN HỒI BA CHIỀU 
LÊ ĐỖ KIÊN, VƢƠNG VĂN THÀNH 
NGHIÊM MẠNH HIẾN* 
The meshless method for single pile behavior in tri-dimentioned 
elastic medium 
Abstract: This paper presents a novel method to analyze the behavior of the 
pile-soil system in a linear elastic soil medium based on the meshless method. 
The meshless method used in this study is Moving Least Square (MLS). 
Results of an analysis of single pile under vertical load using meshless 
method are good agreement with the results from finite element analysis. 
1. ĐẶT VẤN ĐỀ * 
Phƣơng pháp phần tử hữu hạn (PTHH) là 
một phƣơng pháp phổ biến áp dụng vào cơ học 
tính toán trong nhiều thập kỷ qua, phƣơng pháp 
này đã có những đóng góp đáng kể cho sự phát 
triển của khoa học và kỹ thuật. Tuy nhiên, 
phƣơng pháp PTHH không hoàn toàn phù hợp 
với các vấn đề có lƣới biến dạng phức tạp của 
vật liệu hay trƣờng hợp xuất hiện những biến 
dạng không liên tục nhƣ lan truyền vết nứt dọc 
theo đƣờng bất kỳ và các vết nứt phức tạp. Bên 
cạnh đó, phƣơng pháp PTHH cũng gặp khó 
khăn liên quan đến việc chia lƣới và chia lại 
lƣới trong vấn đề tối ƣu hóa lƣới phần tử hoặc 
trong phân tích ảnh hƣởng của vật liệu đa miền. 
Khác với các phƣơng pháp PTHH, phƣơng 
pháp không lƣới chỉ sử dụng một tập hợp các 
điểm nút, các xấp xỉ và hàm dạng đƣợc xây dựng 
hoàn toàn dựa trên các nút, không sử dụng lƣới 
hoặc các phần tử trong phƣơng pháp này. Điều 
này hạn chế đƣợc những khó khăn liên quan đến 
hệ lƣới và đƣa ra một cách tiếp cận linh hoạt hơn 
trong các ứng dụng vào tính toán cơ học. 
Phƣơng pháp không lƣới bắt đầu đƣợc phát 
triển từ những năm 1980, đến nay đã có rất nhiều 
*
 Trường Đại học Kiến trúc Hà Nội 
 DĐ: 0972056219 
 Email: kienlicogi86@gmail.com 
phƣơng pháp không lƣới khác nhau đƣợc xây 
dựng nhƣ: Phƣơng pháp không lƣới bình phƣơng 
di chuyển nhỏ nhất MLS, phƣơng pháp không 
lƣới cục bộ Petrov-Galerkin (MLPG), phƣơng 
pháp không lƣới sử dụng tích phân điểm PIM [3] 
 Trong bài báo, tác giả trình bày quy trình cụ 
thể của phƣơng pháp không lƣới bình phƣơng di 
chuyển nhỏ nhất (MLS) áp dụng cho bài toán địa 
kỹ thuật, đồng thời phát triển phƣơng pháp này 
để phân tích bài toán tƣơng tác của cọc đơn và 
nền đất đàn hồi ba chiều. 
2. PHƢƠNG PHÁP KHÔNG LƢỚI ÁP 
DỤNG TRONG BÀI TOÁN BA CHIỀU 
Quy trình của phƣơng pháp không lƣới bình 
phƣơng di chuyển nhỏ nhất MLS cũng giống 
nhƣ các phƣơng pháp không lƣới khác, đều bao 
gồm 2 bƣớc [2],[4]: 
- Bƣớc 1: Lập hàm dạng. 
- Bƣớc 2: Phân tích không lƣới. 
2.1. Lập hàm dạng 
Hàm dạng không lƣới đƣợc xây dựng thông 
qua các hàm xấp xỉ hoàn toàn dựa trên các nút. 
Xét một hàm vô hƣớng chƣa xác định của một 
biến trƣờng u(x) trong miền . Các xấp xỉ bình 
phƣơng di chuyển nhỏ nhất MLS của u(x) đƣợc 
xác nghĩa nhƣ sau [4]: 
1
( ) ( ).a ( ) ( ).a( )
m
h T
i i
i
u x p x x p x x
  (1) 
ĐỊA KỸ THUẬT SỐ 3-2015 47 
Trong đó: 
- p(x): là hàm cơ sở của các không gian tọa 
độ x, các hàm cơ sở p(x) đƣợc xây dựng từ tam 
giác Pascal. 
- p
t
 : là hàm chuyển của p. 
Trong không gian ba chiều, hàm cơ sở 
( )Tp x bậc 2 đƣợc định nghĩa nhƣ sau [2]: 
2 2 2( ) 1, , , , , , , , ,Tp x x y z x y z xy yz zx . 
- m: là số lƣợng các hàm cơ sở. 
- a(x) : là các hệ số tƣơng ứng và là hàm của 
tọa độ không gian x. Số lƣợng các hệ số a(x) 
phụ thuộc vào bậc và kích thƣớc của hàm cơ sở. 
Hệ số a(x) đƣợc xác định theo phƣơng trình 
tuyến tính sau [2],[4]: 
A(x)a(x)=B(x)U 
hay a(x) = A
-1
(x).B(x)U (2) 
trong đó 1 2{u ,u ,...,u }
T
n
U (3) 
1
W ( )p( )p ( )
n
t
I I I
I
A x x x
 (4) 
1 1 2 2
{ ( )p( ), ( )p( ),.., ( )p( )}
n n
B W x x W x x W x x (5) 
W( )
i
x : là hàm trọng số tại nút thứ I. 
Tác giả lựa chọn các miền hỗ trợ có dạng 
hình chữ nhật, kích thƣớc của miền hỗ trợ theo 
các hƣớng x, y và z tƣơng ứng là dsx , dsy và dsz. 
Hàm trọng số tƣơng ứng với miền hỗ trợ hình 
chữ nhật đƣợc xác định nhƣ sau: 
Wi(x)= W ix(x). W iy(x). 
W iz(x) = W rx. W ry. W rz (6) 
với W ix(x), W iy(x) và W iz(x) là hàm trọng 
số tiêu chuẩn theo hƣớng x, y và z. Các hàm 
trọng số có dạng đƣờng cong bậc 4 đƣợc xác 
định theo GR Liu và Liu [2],[4]: 
2 3 4
ix ix ix ix
ix
ix
1 6r 8r 3r víi 0 r 1
W(r )
0 víi r 1
2 3 4
iy iy iy iy
iy
iy
1 6r 8r 3r víi 0 r 1
W(r )
0 víi r 1
 (7) 
2 3 4
iz iz iz iz
iz
iz
1 6r 8r 3r víi 0 r 1
W(r )
0 víi r 1
Với 
 i
ix
sx
x x
r
d
; 
 i
iy
sy
y y
r
d
và 
 i
iz
sz
z z
r
d
Kết hợp (1) và (2), xấp xỉ u(x) có thể đƣợc 
biểu diễn nhƣ sau: 
' 1( ) ( )A ( ) B(x) U ( ) Uhu x p x x x  (8) 
trong đó: 
' 1( ) ( )A ( ) B(x)x p x x  là hàm dạng. 
2.2. Phân tích không lƣới 
Xét vấn đề của cơ học vật rắn đàn hồi tuyến 
tính trong một miền Ω đƣợc giới hạn bởi biên Γ. 
Hệ phƣơng trình vi phân từng phần và điều kiện 
biên đƣợc viết dƣới dạng sau [1],[2]: 
- Phƣơng trình cân bằng: 
0TL b trong Ω (9). 
- Điều kiện biên tự nhiên: 
n t trên Γt (10) 
- Điều kiện biên cần thiết: 
u u trên Γu (11) 
Trong đó: 
-  : Véc tơ ứng suất. 
- u: Véc tơ chuyển vị, đối với vấn đề 3 chiều, 
x
y
z
u
u v

 
 
 
- b: Véc tơ lực khối. 
- t : Lực kéo quy ƣớc trên lực kéo biên (biên 
tự nhiên). 
- u : Chuyển vị quy ƣớc trên chuyển vị biên 
(biên cần thiết). 
- n: Các véc tơ đơn vị tại một điểm trên biên 
tự nhiên. 
- L: Toán tử khác biệt, đối với vấn đề 3 chiều 
 ĐỊA KỸ THUẬT SỐ 3-2015 48 
0 0
0 0
0 0
0
0
0
x
y
z
L
z y
z x
y x
  
  
  
    
    
    
Các biến phân tiêu chuẩn hình thức dạng yếu của phƣơng trình (9) có dạng sau [2]: 
 ( ) ( ) 0
t
T T TL u DLu d u bd u td  
  
   (12). 
D là ma trận ứng suất – biến dạng, đối với vật liệu đẳng hƣớng: 
2 2
2 2
2 2
1
3
3
3
1 0 0 0
1 0 0 0
1 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
D
với 1
(1 )
(1 )(1 2 )
E 
 
; 
2
1


 và 3
1 2
2(1 )


Sử dụng các hàm dạng không lƣới MLS trên n nút trong các miền hỗ trợ cục bộ: 
( ) ( )
n
h
I I
I
u x x u  hoặc 
0 0
0 0
0 0
II
h
II
n n
h
II I I
I I
I I
u
uu
vu uv


 

   
    
   
  (13) 
Sử dụng phƣơng trình (13), Luh trở thành: 
n
h
I I
I
Lu L u  = 
0 0
0 0
0 0
0 0
0 0 .
0
0 0
0
0
I
n
I I
I
I
x
y
z
u
z y
z x
y x



  
  
  
        
    
 
,
,y
,z
,z ,y
,z ,
,y ,
0 0
0 0
0 0
0
0
0
I x
I
n
I
I I I
II I
I I x
I I x
u B u



 
 
 
 (14) 
Trong đó ,I x , ,yI và ,zI là các đạo hàm của hàm dạng MLS đối với x, y và z. 
ĐỊA KỸ THUẬT SỐ 3-2015 49 
BI là ma trận biến dạng của nút I. 
Thay phƣơng trình (13) và (14) vào phƣơng trình (12) trở thành: 
d d d 0
t
T T T
n n n n
I I J I I I I I
I J I I
B u D B u u b u t  
  
     
    (15) 
- Xét thành phần thứ nhất trong phƣơng trình (15) 
d d
T
n n n n
T T
I I J I I I J I
I J I J
B u D B u u B D B u 
 
  
    
= 
IJ
IJ
B DB d .
n nn n
T T T
I I J J I J
I J I J
K
u u u K u 

   
 = 1 11 1 2 12 2 1...
T T T
N N N
u K u u K u u K u   + 2 21 1 2 22 2 2 2...
T T T
N N
u K u u K u u K u   
+.+ 1 1 2 2 ...
T T T
N N N N N NN N
u K u u K u u K u   = TU KU (16) 
Với K: là ma trận độ cứng tổng thể đƣợc xây dựng từ ma trận độ cứng của các nút. 
 U: là véc tơ chuyển vị tổng thể đƣợc xây dựng từ véc tơ chuyển vị của các nút. 
- Tiếp theo, xét thành phần thứ hai trong phƣơng trình (15): 
Tu bd

 = bd d
b
I
T
n n n
T T T b
I I I I I I
I I I
F
u u b u F  
 
    
   (17) 
với F1
b
 là véc tơ lực khối của nút, b T
I I
F bd

   
Vế phải của phƣơng trình (17) : 
n
T b
I I
I
u F = 1 1 2 1 ...
T b T b T b
N N
u F u F u F   
= 
1
1 (1 2 )
(2 1)
... ...
b
T T
N x N
b
N Nx
F
u u
F
 
 
 
 
= T bU F (18) 
F
b
 là véc tơ lực khối tổng thể đƣợc tập hợp từ 
các vectơ lực khối của tất cả các nút trong toàn 
bộ miền tính toán. 
Thực hiện tƣơng tự với thành phần thứ 3 của 
phƣơng trình (15), véc tơ lực khối đƣợc thay thế 
bởi các véc tơ lực kéo trên biên tự nhiên và tích 
phân trên miền biên tự nhiên Γ. Các véc tơ lực 
kéo tại nút là: dt TI IF t

   (19) 
Kết hợp các phƣơng trình (16), (18) và (19), 
phƣơng trình (15) trở thành: 
( ) ( ) 0T T b T tU KU U F U F   
Hoặc ( ) ( ) 0T b tU KU F F (20) 
Do U là bất kỳ, phƣơng trình (20) thỏa 
mãn chỉ khi: 
( ) ( ) 0b tKU F F hoặc KU F (21) 
với F là véc tơ lực khối tổng thể: 
( ) ( )b tF F F 
Các chuyển vị nút thu đƣợc bằng cách giải 
phƣơng trình (20), sau đó thông qua mối quan 
hệ tuyến tính giữa ứng suất – biến dạng có thể 
xác định đƣợc trạng thái ứng suất tại các điểm 
trong môi trƣờng đất đàn hồi. 
Tác giả đã xây dựng các chƣơng trình con 
tính hàm dạng và các đạo hàm của hàm dạng, bổ 
sung vào phần mềm SSI3D để tính toán chuyển 
vị và ứng suất tại các điểm của hệ cọc – đất. 
 ĐỊA KỸ THUẬT SỐ 3-2015 50 
3. VÍ DỤ MINH HỌA 
Tính toán cọc đơn có đƣờng kính 1,0m ; 
chiều dài 30m chịu tải trọng đứng tại đỉnh cọc 
P= 1000 tấn. Các đặc trƣng của vật liệu cọc và 
môi trƣờng mà cọc nằm trong đƣợc trình bày 
trong bảng 1 và bảng 2. Do tính đối xứng nên 
chỉ 1/4 mô hình thực tế đƣợc xây dựng để giảm 
thời gian tính toán, mô hình tính toán đƣợc trình 
bày trong hình 1. 
Bảng 1: Đặc trƣng vật liệu làm cọc 
Đặc trƣng Đơn vị Giá trị 
Mô đun đàn hồi T/m2 2700000 
Hệ số Poisson - 0.2 
Trọng lƣợng riêng T/m3 2.5 
Bảng 2: Đặc trƣng đất nền 
Đặc trƣng Đơn vị Giá trị 
Mô đun đàn hồi T/m2 4000 
Hệ số Poisson - 0.3 
Trọng lƣợng riêng T/m3 1.9 
a) b) c) 
Hình 1: Vị trí các nút 
a) không gian b) mặt bằng c) mặt đứng 
Kết quả tính toán chuyển vị của cọc theo độ 
sâu đƣợc trình bày trong hình 2. Mô hình tƣơng 
tự đƣợc xây dựng trên phần mềm Plaxis 2D theo 
bài toán đối xứng trục. Kết quả chuyển vị tính 
toán thu đƣợc tại đỉnh cọc là 0.017 m và tại mũi 
cọc là 0.01 m, kết quả này phù hợp với kết quả 
tính toán theo phƣơng pháp không lƣới bằng 
phần mềm SSI3D. 
STT Điểm 
Tọa độ 
điểm (X) 
Tọa độ 
điểm (Y) 
Chuyển vị thẳng 
đứng UY (m) 
STT Điểm 
Tọa độ 
điểm 
(X) 
Tọa độ 
điểm 
(Y) 
Chuyển vị 
thẳng đứng 
UY (m) 
1 16178 0 0 0,017144468 16 10163 0 -15 0,01186453 
2 15777 0 -1 0,016591175 17 9762 0 -16 0,011636637 
3 15376 0 -2 0,016167579 18 9361 0 -17 0,011422201 
4 14975 0 -3 0,015735774 19 8960 0 -18 0,011221083 
5 14574 0 -4 0,015327275 20 8559 0 -19 0,011033172 
6 14173 0 -5 0,014934254 21 8158 0 -20 0,010858398 
ĐỊA KỸ THUẬT SỐ 3-2015 51 
STT Điểm 
Tọa độ 
điểm (X) 
Tọa độ 
điểm (Y) 
Chuyển vị thẳng 
đứng UY (m) 
STT Điểm 
Tọa độ 
điểm 
(X) 
Tọa độ 
điểm 
(Y) 
Chuyển vị 
thẳng đứng 
UY (m) 
7 13772 0 -6 0,014558229 22 7757 0 -21 0,010696728 
8 13371 0 -7 0,014198372 23 7356 0 -22 0,010548174 
9 12970 0 -8 0,013854434 24 6955 0 -23 0,010412801 
10 12569 0 -9 0,013526035 25 6554 0 -24 0,010290737 
11 12168 0 -10 0,013212843 26 6153 0 -25 0,010182188 
12 11767 0 -11 0,012914536 27 5752 0 -26 0,01008744 
13 11366 0 -12 0,012630816 28 5351 0 -27 0,010007087 
14 10965 0 -13 0,012361406 29 4950 0 -28 0,009941257 
15 10564 0 -14 0,012106053 30 4549 0 -29 0,009892469 
16 10163 0 -15 0,01186453 31 4035 0 -30 0,00986384 
Hình 2: Chuyển vị nút theo độ sâu 
4. KẾT LUẬN 
Trong bài báo, tác giả đã trình bày quy trình 
cụ thể của phƣơng pháp không lƣới bình 
phƣơng di chuyển nhỏ nhất và vận dụng phƣơng 
pháp này cho bài toán tính toán cọc đơn trong 
môi trƣờng đất nền đàn hồi tuyến tính. Ví dụ 
tính toán đối với cọc đơn chịu tải trọng đứng 
cho kết quả tính toán phù hợp với kết quả tính 
toán theo phần mềm Plaxis 2D. 
TÀI LIỆU THAM KHẢO 
1. G.R. Liu (2003); “Meshfree Method: 
Moving beyond the finite element Method”. 
National University of Singapore, Singapore. 
2. G.R. Liu and Y.T. Gu, (2003); “An 
Introduction to Meshfree Methods and Their 
Programming”. National University of 
Singapore, Singapore. 
3. Huafeng Liu and Pengcheng Shi, (2003); 
“Meshfree Particle Method”. Department of 
Electrical and Electronic Engineering Hong 
Kong University of Science and Technology, 
Hong Kong. 
4. Youping Chen, James D. Lee and Azim 
Eskandarian, (2006); “Meshless Methods in Solid 
Mechanics”. Springer Science+Business Media, 
Inc., 233 Spring Street, New York, USA. 
Người phản biện: GS.TS. ĐỖ NHƢ TRÁNG