Tính toán công trình ngầm trong nền san hô chịu tác dụng đồng thời của động đất và lực thuỷ tĩnh

1 Tạp chí Khoa học và Công nghệ biển T10 (2010). Số 3. Tr 01 - 13 TÍNH TOÁN CÔNG TRÌNH NGẦM TRONG NỀN SAN HÔ CHỊU TÁC DỤNG ðỒNG THỜI CỦA ðỘNG ðẤT VÀ LỰC THUỶ TĨNH NGUYỄN TẤT NGÂN Học viện Kỹ thuật quân sự Tóm tắt. Bài báo trình bày phương pháp tính tương tác giữa công trình ngầm trong nền san hô chịu tác dụng của tải trọng do ñộng ñất và lực thuỷ tĩnh nên. Bài toán ñược giải quyết bằng mô hình phẳng, trên cơ sở phương pháp phần tử hữu hạn (PTHH), kết hợp phần biến dạng phẳng dạng tam giác 3 ñiểm nút và phần tử tiếp xúc hai chiều (2D). Với thuật toán kết hợp tích phân trực tiếp Newmark và phương pháp lặp Newton-Raphson tác giả lập trình tính toán trong môi trường Matlab ñể giải phương trình phi tuyến liên kết. Kết quả tính toán số phù hợp quy luật cơ học và toán học cho thấy khả năng của chương trình ñã lập. Nội dung của bài báo có thể làm tài liệu tham khảo cho việc nghiên cứu và tính toán, thiết kế các kết cấu công trình ngầm trong nền san hô. I. MỞ ðẦU ðộng ñất là hiện tương thiên nhiên ñã và sẽ gây ra nhiều thảm hoạ ñối với con người và các công trình xây dựng. Từ trước tới nay, con người ñã có rất nhiều nỗ lực trong việc nghiên cứu phòng, chống ñộng ñất ñể bảo vệ sinh mạng của mình và tài sản vật chất xã hội. Một trong những hướng nghiên cứu hiện nay là ñưa ra các phương pháp tính toán, thiết kế công trình ñảm bảo giảm thiểu tác ñộng của ñộng ñất, nâng cao hiệu quả sử dụng, tăng tuổi thọ của công trình. Mặc dù ñã có những thành công ñáng kể trong lĩnh vực này, song thực tế vẫn còn nhiều thảm hoạ do ñộng ñất gây ra trên Thế giới hiện nay. Khi ñộng ñất xảy ra, do hiện tượng dao ñộng mạnh của nền ñất, nên kết cấu cũng bị dao ñộng, do ñó tự bản thân kết cấu (hệ kết cấu - nền) cũng chịu tác dụng của lực quán tính do gia tốc chuyển ñộng của nền (gia tốc nền) gây ra. Do vậy, trong tính toán công trình chịu tác dụng của ñộng ñất, gia tốc nền là một trong những yếu tố rất quan trọng. Các công trình trên các ñảo san hô xa bờ cũng có khả năng chịu tác ñộng của ñộng ñất, khi các công trình này hư hỏng, ngoài sự thiệt hại về vật chất, con người như các công trình trong ñất liền thì vấn ñề bảo ñảm ñiều kiện tác chiến trên ñảo không còn, ảnh hưởng lớn ñến khả năng giữ vững chủ quyền trên biển - ñiều này là một trong những ñặc thù của 2 các công trình biển, ñảo. Vì vậy, nghiên cứu tính toán các công trình ngầm trong nền san hô chịu tác dụng của tải trọng ñộng ñất là vấn ñề cần thiết. Hình 1: Gia tốc nền trong trận ñộng ñất năm 1995 tại Kobe - Nhật Bản [5] II. CÁC GIẢ THIẾT, MÔ HÌNH VÀ PHƯƠNG PHÁP TÍNH 1. Các giả thiết và mô hình tính của bài toán Bài toán ñược giải quyết dựa vào các giả thiết sau: • Kết cấu công trình ñàn hồi, biến dạng tuyến tính. • Mỗi lớp nền san hô là vật liệu ñồng nhất, ñẳng hướng, ñàn hồi tuyến tính. Quá trình kết cấu làm việc, không có hiện tượng tách, trượt giữa các lớp nền với nhau. • Hệ kết cấu công trình và nền làm việc trong ñiều kiện biến dạng phẳng. Liên kết giữa kết cấu và nền san hô ñược thay thế bằng liên kết nút giữa các phần tử biến dạng phẳng 2D thông qua liên kết với phần tử tiếp xúc Goodman [1, 3]. Liên kết tiếp xúc giữa kết cấu và nền san hô là liên kết một chiều. • Khi tính toán, ñối với vật liệu san hô, bỏ qua lực dính giữa kết cấu và nền. Theo hướng này, tách từ hệ thực bán vô hạn ra một miền hữu hạn bao gồm kết cấu và một phần nền gọi là miền nghiên cứu, trên biên miền nghiên cứu ñược ñặt các liên kết, việc tính toán ñược thực hiện trên miền nghiên cứu ñã xác ñịnh. Kích thước biên của miền nghiên cứu ñược xác ñịnh theo phương pháp lặp [1,7] (hình 2). 3 Hình 2: Mô hình thực của bài toán 2. Phương pháp tính và các loại PTHH ðể tính toán kết cấu theo mô hình nói trên, tác giả sử dụng phương pháp PTHH, trong ñó các loại phần tử sử dụng là: ðối với nền san hô và kết cấu sử dụng loại phần tử tam giác phẳng 3 ñiểm nút. ðối với lớp tiếp xúc giữa kết cấu và nền san hô sử dụng loại phần tử tiếp xúc tứ giác 4 ñiểm nút (four nodes slip element). III. THIẾT LẬP CÁC PHƯƠNG TRÌNH CƠ BẢN BẰNG PHƯƠNG PHÁP PTHH 1. Các quan hệ ñối với phần tử thuộc nền và kết cấu Xét phần tử tam giác phẳng 3 ñiểm nút, việc ñánh số thứ tự nút của phần tử và lựa chọn hệ toạ ñộ cục bộ của phần tử ñược quy ước như trên hình 3. Chuyển vị tại một ñiểm bất kỳ trong phạm vi phần tử ñược nội suy từ véc tơ chuyển vị nút của phần tử theo biểu thức [2]: { } [ ] { }mm UNu = (1) trong ñó: { } { }Tyx uuu = - véc tơ chuyển vị tại ñiểm (x,y) thuộc phần tử, [ ]mN - ma trận hàm dạng của phần tử thứ m, 4 { }mU - véc tơ chuyển vị nút của phần tử thứ m, nó có dạng: { } { }Ty3x3y2x2y1x1m UUUUUUU = (2) Ma trận hàm dạng của phần tử có dạng: [ ]       = 322 321 m N0N0N0 0N0N0N N (3) r s HÖ täa ®é côc bé (r, s) HÖ täa ®é chung (x, y) 1 (0, 0) 3 (0, 1) x 2 (1, 0) y Hình 3: Số hiệu nút của phần tử, hệ toạ ñộ chung và toạ ñộ cục bộ của phần tử Biểu thức hàm dạng của phần tử tam giác viết trong toạ ñộ cục bộ và tọa ñộ chung ñược cho trong bảng 1. Bảng 1: Hàm dạng của phần tử tam giác Toạ ñộ cục bộ Toạ ñộ chung N1 = 1- r - s N1 = (a1 + b1x + c1y)/2∆ N2 = r N2 = (a2 + b2x + c2y)/2∆ N3 = s N3 = (a3 + b3x + c3y)/2∆ trong ñó 33 22 11 yx1 yx1 yx1 2 1 =∆ - là diện tích phần tử, với: 123312231 213132321 122133113223321 xxc;xxc;xxc ;yyb;yyb;yyb ;yxyxa;yxyxa;yxyxa −=−=−= −=−=−= −=−=−= . 5 Véc tơ biến dạng của phần tử{ }ε ñược viết như sau: { } [ ] { }mm UB=ε (4) Ma trận biến dạng-chuyển vị của phần tử [ ]mB ñược xác ñịnh: [ ]                   ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = x N y N x N y N x N y N y N0 y N0 y N0 0 x N0 x N0 x N B 332211 321 321 m (5) Ma trận ñộ cứng phần tử ñược xác ñịnh theo biểu thức sau: [ ] [ ] [ ][ ]∫= mA mm T mm dABDBhK (6) trong ñó: [ ]D là ma trận vật liệu phần tử, phụ thuộc vào modul Young E, hệ số Possion ν. Sử dụng phương phương pháp tích phân số cầu phương Gauss-Legendre, biểu thức (6) ñược viết lại: [ ] [ ] [ ][ ] i2i1ii,mn 1i T i,mm WWJdetBDBhK ∑ = = (7) trong ñó: i - ñiểm tích phân thứ i có toạ ñộ (ri, si), n - tổng số ñiểm tích phân, detJ - ñịnh thức của ma trận Jacobi, W1i, W2i - trọng số của phép tích phân cầu phương, chỉ số i biểu diễn giá trị của hàm tại ñiểm Gauss thứ i. Ma trận khối lượng của phần tử ñược xác ñịnh theo công thức sau: [ ] [ ] [ ]∫ρ= mA mm T mm dANNhM (8) và cũng ñược tính theo phương pháp tích phân số cầu phương Gauss-Legendre. trong ñó: ρ - khối lượng riêng của vật liệu của phần tử thứ m; Am - diện tích phần tử thứ m; h - chiều dày phần tử, 2. Các quan hệ ñối với phần tử thuộc lớp tiếp xúc (PTTX) Phần tử tiếp xúc 2 chiều ñược bố trí vào lớp tiếp xúc giữa kết cấu công trình và nền san hô thể hiện tính chất liên kết một chiều của nền. Mô hình phần tử tiếp xúc có chiều dày “bằng không” ñược chỉ ra như trên hình 4, do chiều dày t của phần tử gần như bằng 0 nên các cặp nút 1 và 4, 2 và 3, 5 và 6 có cùng toạ ñộ. 6 (a) (b) (c) Hình 4: Mô hình phần tử tiếp xúc Goodman (1÷4 – các nút). a) Sơ ñồ hình học của phần tử tiếp xúc; b) Quan hệ ứng suất pháp tuyến với biến dạng pháp tuyến; c) Quan hệ ứng suất tiếp tuyến với biến dạng tiếp tuyến Quan hệ giữa số gia ứng suất và số gia biến dạng: [ ]       γ∆ ε∆ =       τ∆ σ∆ seD (9) trong ñó ma trận [Dse] ñược xác ñịnh theo biểu thức: [ ]       = ξ η k0 0k Dse (10) với: kη và kξ tương ứng là ñộ cứng pháp tuyến và ñộ cứng tiếp tuyến của phần tử. Các hệ số này có thể ñược xác ñịnh bằng thực nghiệm hoặc xác ñịnh thông qua các ñặc trưng của nền [1,3]: ( )( )( )ν−ν+ ν− =η 211 1Ek ; ( )ν+==ξ 12 EGk (11) với: ν và E tương ứng là hệ số Poisson và mô ñun Young của vật liệu nền. Chuyển vị một ñiểm bất kỳ của phần tử trong hệ trục toạ ñộ ñịa phương ñược xác ñịnh thông qua chuyển vị của nó trong hệ toạ ñộ tổng thể và góc lệch trục α:               ββ− ββ =       v u cossin sincos v u 1 1 (12) Ma trận ñộ cứng trong hệ toạ ñộ chung: [ ] [ ] [ ][ ]dxdyBDBK sese T sese ∫∫= (13) trong ñó: [ ] [ ]localse B cossin sincos B         ββ β−β = (14) ðối với PTTX 4 ñiểm nút, biến dạng trong phần tử là ñều, ta có: 7 [ ]         −− −− = 0N0N0N0N N0N0N0N0 B ' 4 ' 3 ' 2 ' 1 ' 4 ' 3 ' 2 ' 1 local với: ξ= d dNN i'i , ξ (Từ -1 ñến +1) dọc theo chiều dài phần tử, Ni là hàm dạng của PTTX [1, 3]. Sử dụng phép biến ñổi ∫∫ ∫ ξ= − 1 1 dJdetdxdy , ma trận ñộ cứng PTTX ñược xác ñịnh: [ ] [ ] [ ][ ]∫ ξ= − 1 1 sese T sese dJdetBDBK (15) trong ñó [J] là ma trận Jacobin [1,3]. Ma trận ñộ cứng của PTTX trong công thức (15) ñược xác ñịnh bằng phương pháp tích phân số cầu phương Gauss-Legendre và sử dụng sơ ñồ 1 ñiểm tích phân (ξ = 0, w = 2). 4. Phương trình chuyển ñộng của kết cấu Khi ñộng ñất, gia tốc nền sẽ gây ra gia tốc chuyển vị của các khối lượng tập trung tại các nút của kết cấu { }gUɺɺ , ñược biểu diễn [4,5,6,7]: { } { } gg urU ɺɺɺɺ = . (16) trong ñó: { }gUɺɺ là véc tơ có số chiều bằng bằng số chiều véc tơ nút { }U của hệ, { }r là véc tơ hệ số ảnh hưởng chuyển vị nút biểu thị chuyển vị tại các nút do chuyển vị ngang của nền bằng ñơn vị gây ra, cấu trúc véc tơ gồm các phần tử 1 và 0 - tuỳ vào bậc tự do của mỗi nút và số chiều của gia tốc nền. Gọi { }eR là véc tơ tải trọng do ñộng ñất tác dụng lên các nút (là lực quán tính do gia tốc nền gây ra ñối với ñối với các khối lượng quy nút) của hệ. Lúc này ta có: { } [ ]{ } [ ]{ } gge urMUMR ɺɺɺɺ −=−= . (17) Phương trình chuyển ñộng của hệ kết cấu chịu tác dụng ñồng thời của tải trọng ñã cho và ñộng ñất, có dạng [4,5,6,7]: [ ]{ } [ ]{ } [ ]{ } { } { }ed RRUKUCUM +=++ ɺɺɺ , (18) trong ñó: { }dR là véc tơ tải trọng nút do tải trọng ñã cho gây ra. 8 IV. THUẬT TOÁN GIẢI PHƯƠNG TRÌNH CHUYỂN ðỘNG CỦA KẾT CẤU Sau khi ñưa vào các ñiều kiện biên, (18) có dạng tổng quát như sau [2, 4, 5, 6]: [ ]{ } [ ]{ } [ ]{ } { } { }ed RRUKUCUM +=++ ɺɺɺ , (19) Như ñã biết, do kể ñến sự tách, trượt cục bộ giữa kết cấu và nền san hô, phương trình chuyển ñộng (19) là phương trình phi tuyến hình học và ñược viết lại như sau: [ ]{ } { }( )[ ]{ } { }( )[ ]{ } { } { } { }RRRUUKUUCUM ed =+=++ ɺɺɺ (20) Ma trận cản tổng thể ñược xác ñịnh: [ ] [ ] [ ] { }( )[ ]UCKMC =β+α= , với: 21 21 ; 2 ωβω=α ω+ω ξ =β ; (ω1, ω2 là 2 tần số ñầu tiên trong dao ñộng tự do, ξ là tỷ số cản). ðể giải phương trình chuyển ñộng (20) tác giả sử dụng phương pháp tích phân trực tiếp Newmark kết hợp với phương pháp lặp Newton – Raphson, với các bước tóm tắt thuật giải như trong bảng 2. Bảng 2: Tóm tắt thuật toán giải bài toán tương tác kết cấu công trình và nền san hô chịu tác dụng ñộng ñất theo phương pháp Newmark kết hợp với lặp Newton-Raphson I. Nhập dữ liệu ban ñầu: 1. Sơ ñồ, kích thước kết cấu 2. ðặc trưng vật liệu 3. Nhập tải trọng (nếu có) 4. Nhập giản ñồ gia tốc nền )t(u gɺɺ II. Tính toán sơ bộ ban ñầu: 1. Tính các ma trận cơ bản: [K], [M], [C], véc tơ ảnh hưởng { }r . 2. Xác ñịnh bước tích phân ∆t và các tham số tích phân α, δ. Tính các hệ số của phương pháp: a0, a1, a2, a3, a4, ..., a7. 3. Xác ñịnh các véc tơ ñiều kiện ban ñầu: { } { }{ }000 U,U,U ɺɺɺ . 4. Tính ma trận ñộ cứng hiệu quả [ ]∗K . 9 II. Tính lặp cho các bước tích phân tại các thời ñiểm: t = ∆t, t = 2∆t, t=3∆t, ... theo các bước sau: 1. Tính véc tơ tải trọng hiệu quả { }∗ ∆+ ttR tại thời ñiểm t + ∆t 2. Dùng phương pháp lặp xác ñịnh véc tơ chuyển vị nút{ }ttU ∆+ tại thời ñiểm t + ∆t: [ ]{ } { } { } )1i(tt)1i(tt)i(tt FRUK −∆+−∗ ∆+∗ ∆+ −=∆ , { } { } { } )i()1i(tt)i(tt UUU ∆+= −∆+∆+ , Kiểm tra ñiều kiện hội tụ: { } { } { } Dt)i(tt )i( UU U ε≤ − ∆ ∆+ 3. Xác ñịnh các véc tơ gia tốc và vận tốc tại thời ñiểm t + ∆t: { } { } { }( ) { } { } { } { } { } { }tt7t6ttt t3t2ttt0tt UaUaUU UaUaUUaU ∆+∆+ ∆+∆+ ++= −−−= ɺɺɺɺɺɺ ɺɺɺɺɺ 4. Tính ứng suất, biến dạng và nội lực tại thời ñiểm t + ∆t. 5. Kiểm tra ñiều kiện tách, trượt, cập nhật ma trận phần tử. 6. Tăng bước thời gian t = t + ∆t và tiến hành lặp lại từ bước 1. V. VÍ DỤ SỐ Trên cơ sở các thuật toán ñã thiết lập, tác giả ñã lập chương trình tính có tên là Eartquake_2010. Chương trình ñược viết bằng ngôn ngữ MATLAB chạy trên máy tính sử dụng hệ ñiều hành Windows xp. Tính kết cấu công sự trong nền san hô chịu ñồng thời của lực thuỷ tĩnh tác dụng trực tiếp lên kết cấu và tải trọng ñộng ñất theo phương ngang, giản ñồ gia tốc nền như hình 5. Sơ ñồ bài toán như hình 2, với các số liệu cụ thể như sau: - Lực thuỷ tĩnh: ( )       +−= 2 1 H yptp 00 , với p0 = 200N/cm 2 . - Kích thước hình bao của kết cấu: chiều cao H = 4 m, bề rộng B = 5 m, chiều dày t = 0,5 m. Vật liệu có E = 0,34×107N/cm2, ν = 0,3, ρ = 2,5×10-3 kg/cm3. 10 Hình 5: Biểu ñồ gia tốc nền với giá trị cực ñại amax = 170,85cm/s2 [6] Sau khi lặp ứng với gia tốc nền lớn nhất, xác ñịnh ñược kích thước miền tính toán: B0 = 20 m, chiều cao H0 = 14 m, với sai số chuyển vị ε = 2,16%, tương ứng với 5 lớp nền (thứ tự từ 1 ñến 5) có ñặc trưng cơ lý như bảng 3. Bảng 3: Thuộc tính vật liệu nền san hô [1] Lớp ðộ sâu (m) E (N/cm2) ν ρ (kg/cm3) Hệ số ma sát f Hệ số giảm chấn ξ 1 1 2,8.104 0,22 2,5.10-3 0,21 2 4 2,1.105 0,25 2,8.10-3 0,32 3 8 2,0.106 0,25 2,9.10-3 0,41 4 10 2,6.105 0,25 2,0.10-3 0,47 5 14 2,3.106 0,25 3,0.10-3 0,38 0,05 Mô hình PTHH của bài toán như hình 6. 11 Hình 6: Mô hình PTHH của bài toán Số lượng phần tử tổng cộng là 186, bao gồm: 32 phần tử kết cấu, 20 phần tử tiếp xúc và 134 phần tử nền. Thời gian tính toán tcal = 20s, bước thời gian ∆t = 0,02s. Kết quả tính: Hình 7 và 8 là ví dụ các kết quả ñáp ứng chuyển vị ngang U và ứng suất σx theo thời gian của ñiểm C (ðiểm bên trái công sự). Hình 7: Sự thay ñổi chuyển vị ngang tại ñiểm A thuộc ñỉnh công sự 12 Hình 8: Sự thay ñổi ứng suất σy tại ñiểm C thuộc ñỉnh công sự VI. KẾT LUẬN Nội dung bài báo ñã ñạt ñược: - Xây dựng mô hình, thuật toán PTHH và chương trình phần mềm trong môi trường Matlab giải bài toán tương tác giữa kết cấu công trình và nền san hô chịu tác dụng ñồng thời của lực thuỷ tĩnh và tải trọng do ñộng ñất gây nên. - Tính toán trên ví dụ cụ thể cho thấy kết quả tính phù hợp quy luật cơ học. - Kết quả bài báo có thể làm tài liệu tham khảo và là công cụ tính toán tương tác giữa kết công trình và nền san hô chịu tải trọng ñộng ñất phục vụ an ninh quốc phòng. TÀI LIỆU THAM KHẢO 1. Nguyễn Thái Chung, 2006. Nền san hô và sự làm việc của cọc trong nền san hô, Luận án TSKT, Học viện KTQS. 2. Bathe K.J and Wilson E.L., 1978. Numerical Method in Finite Method Analyis Prentice, Hall of India Private Limited, New Delhi. 13 3. Pizhong Qiao, Jialai Wang, 2004. Mechanics and fracture of crack tip deformable bi-material interface, International Journal of Solids and Structures 41 (2004) 7423 – 7444. 4. Giovanni Li Destri Nicosia, 2008. On seismic design and advanced numerical modelling of diaphragm walls under earthquake loading, Istituto Universittario di studi superiori di pavia universita degli studi di pavia. 5. Chen Guoxing và Zhuang Haiyang, 2009. Analysis on the erathquake response of subway station based on the substructuring suntraction method, Institute of geotechnical engineering, Nanjing university of technology, China. 6. Balideh S.,... 2009. Seismic analysis of underground spaces to propagation of seismic waves, Journal of applied sciences 9 (9), pp.1615-1627. 7. Jun Seong Choi, Jong Seh Lee, Jae Min Kim, 2009. Nonlinear earthquake respomse analysis of 2-D underground structures with soil–structures interaction including separation and sliding at interface,15thASCE engineering mechanics conference, Columbia University, New York. ANALYSIS OF UNDERGROUND STRUCTURE IN CORAL FOUNDATION UNDER EARTHQUAKE AND HYDROSTATIC LOAD NGUYEN TAT NGAN Summary: In this study, the author present the calculating method interaction between underground structure and coral foundation under earthquake and hydrostatic load. The problem analysis with finite element method and using plane strain three nodes triangular element and 2D slip element (for slip layer between underground structure and coral foundation). The paper results can be used in design as reference for the structures in the coral foundation under earthquake loads. Ngày nhận bài: 16 - 7 - 2010 Người nhận xét: PGS. TS. Phạm Tiến ðạt