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

ĐỊ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