Á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