Ứng dụng FEM nghiên cứu bước thời gian thích hợp trong quá trình cố kết - Nguyễn Quang Hùng

ỨNG DỤNG FEM NGHIÊN CỨU BƯỚC THỜI GIAN THÍCH HỢP<br /> TRONG QUÁ TRÌNH CỐ KẾT<br /> Nguyễn Quang Hùng<br /> NCS đại học Vũ Hán, Hồ Bắc Vũ Hán 430072<br /> <br /> <br /> Tóm tắt: Hiện tại vấn đề nâng cao độ chính xác và tính ổn định của phương pháp FEM đang được nhiều nhà khoa học<br /> quan tâm. Bài báo này dựa trên lý thuyết cố kết cơ bản Biot tiến hành nghiên cứu bước thời gian tự thích ứng trong<br /> quá trình cố kết nhằm nâng cao tính chính xác cũng như độ ổn định của FEM. Từ những nguyên lý cơ bản suy luận và<br /> đưa ra công thức khống chế bước thời gian. Từ đây có thể khống chế toàn bộ trường ứng suất cũng như áp lực lỗ rỗng.<br /> Kết quả nghiên cứu nhằm chính xác hóa quá trình phát triển ứng suất – biến dạng , điều này hết sức có ý nghĩa trong<br /> quá trình thiết kế cũng như thi công đập vật liệu địa phương , đặc biệt là đập có tường lõi bằng đất.<br /> <br /> <br /> <br /> 1. Mở đầu d2  Es<br /> 1 <br /> Nghiên cứu trạng thái ứng suất của các công<br /> trình đất đặc biệt là các đập cao có xét đến quá E<br /> d3 <br /> trình cố kết của đất là một nội dung quan trọng 2(1  2  )<br /> trong quá trình phân tích ổn định công trình [16]. E (1   )<br /> Phương pháp phần tử hữu hạn là một công cụ ES <br /> (1   )(1  2 )<br /> mạnh để giải bài toán này. Nội dung của bài báo<br /> trình bày đóng góp mới nhằm nâng cao độ chính<br /> u,v,w lần lượt là chuyển vị theo phương các<br /> xác lời giải phương trình cố kết cơ bản Biot.<br /> trục x,y,z.<br /> 2. Phương trình cố kết cơ bản Biot [2,3,11]<br /> p: áp lực lỗ rỗng.<br /> Các giả thiết cơ bản của lý thuyêt cố kết<br /> 2.2 Phương trình liên tục<br /> Biot:<br /> Dựa vào định luật Dacxi có:<br /> Ngoài tính thấm ra, tính chất của đất là đồng<br /> nhất đẳng hướng. Tính chất thấm của đất không  k h p<br /> v x  <br /> thay đổi theo không gian và thời gian, không xét   w x<br /> đến lực quán tính.  k h p<br /> 2.1 Phương trình cân bằng v y   (2)<br /> Đối với đất đồng nhất đẳng hướng kết hợp   w y<br /> phương trình cân bằng và nguyên lý ứng suất  k p<br /> hiệu quả thu được: v z   h<br />   w z<br />   2u  2u  2u Trong đó<br />  d1 2  d 3 ( 2  2 ) <br />  x y z w=wg : Trọng lượng nước trong lỗ rỗng<br />  2<br />  v  2w p kh : Hệ số thấm của đất.<br />  ( d 2  d 3 )(  ) 0<br />   x  y  x  z x Phương trình liên tục:<br />   2v  2v  2v 1 2 p 2 p 2 p  u v w<br />  d1 2  d 3 ( 2  2 )  [kh (   )]  (   ) (3)<br />  y x z (1) w x 2 y 2 z 2 t x y z<br />  2 2<br />  ( d  d )(  u   w )   p  0 3. Rời rạc hóa phương trình không gian cơ bản<br /> 2 3<br />   y x  y z y 3.1 Rời rạc hóa phươg trình cân bằng<br />  2 2 2<br />   w  w  w [ K ]{ }  [C ]{ p}  {R F } (4)<br />  d1 z 2  d 3 ( y 2  x 2 )  Trong đó<br /> <br />   2u  2v p [ K ]  [ B]T [ D][B]dv là ma trận cứng tổng thể.<br /> (<br />  2 d  d 3 )(  ) 0 V<br />  zx zy z<br /> [C]  [ B]T [M ][N ]dv là ma trận ngẫu hợp.<br /> V<br /> Trong đó<br /> d1=Es {RF }  [ N ]T {F}ds là vector tải trọng ngoài.<br /> S<br />  {R q }  t ({R q }  [ H ]{ p} n ) là vector gia tăng<br /> {F } là vector lực biên. lưu lượng trong thời đoạn t .<br /> [ M ]  [1 1 1 0 0 0]T  là tham số tích phân có thể lựa chọn trong<br /> 3.2 Rời rạc hóa phương trình liên tục phạm vi 0.5~1.0<br /> [C ]T {}  [ H ]{ p}  {R q } (5) 5. Khống chế sai số [5,7,8]<br /> { }<br /> Lựa chọn<br /> Trong đó {}  u nInt1  u n  tu n   1tu n<br /> t<br />  Int (8)<br /> [ H ]   [ B s ]T [k ][ B s ]dv là ma trận thấm u n 1  u n  u n<br /> V Int<br />   p n 1  p n  tp n   2 tp n<br />  Int (9)<br />    p n 1  p n  p n<br />  x <br /> <br /> [ B s ]   N 1 N2 ... N 8 <br /> y Sau khi biến đổi và áp dụng khai triển Taylor<br />  <br />   thu được kết quả như sau :<br />  z <br />  Int E 1 2 3<br /> u n 1  u n 1  1tu n  2 t un  O ( t )<br /> k h 0 0 (10)<br /> 1  <br /> [k ]  0 kh 0  là ma trận hệ số thấm  p Int  p E   tp  1 t 2 p  O ( t 3 )<br /> w   n 1 n 1 2 n<br /> 2<br /> n<br />  0 0 k v <br /> {Rq }   [ N ]T v n ds là ma trận lưu lượng. Sai số cục bố được tính theo công thức:<br /> S<br /> <br />  u 1 Int<br /> 4. Rời rạc hóa phương trình không gian cơ  e   t 2 ( u n  1  u n ) (11)<br /> bản [3,6,12] <br /> Xét miền tính toán tại hai thời điểm t n và  e p   t 1 ( p Int  p )<br />  2<br /> n 1 n<br /> <br /> t n 1 tương ứng có các biến lượng:<br /> { } n , { }n 1 , { p} n , { p} n 1 . Trong thời đoạn<br /> 6. Khống chế bước thời gian[3,10]<br /> Xét tới tác dụng của ứng suất và áp lực lỗ rỗng,<br /> t  t n 1  t n có thể lấy {}  { } / t , sai số Euclid được chọn để khống chế bước thời<br /> { p }  {p} / t . Lựa chọn hàm số tuyến tính: gian trong quá trình tính toán.<br /> {}  {}n  ({}n 1  {}n ) / t<br /> Sai số áp lực lỗ rỗng:<br />  ( : 0  t ) (6) 1<br /> { p}  { p}n  ({ p}n 1  { p}n ) / t  p  e p  t ( p nInt1  p n )<br /> 2<br /> Hoặc e p<br /> <br />  N   F ds<br /> T<br /> R F  <br /> S p max<br /> <br /> [ K ]{ }  [C]{p}  {RF } Sai số ứng suất<br /> (7) 1 Int<br />  T<br /> [C] { } t[ H ]{p}  {Rq }  u  e u  t (u n 1  u n )<br /> Trong đó: 2<br /> Hoặc<br /> Vector gia tăng tải trọng trong thời đoạn t .<br /> eu<br /> là vector gia <br /> Fx  u<br />   tăng tải trọng max<br /> {F}  Fy  biên trong thời<br />   Bước thời gian tính toán trong thời đoạn tiếp<br /> Fz  đoạn t .<br /> theo được tính theo công thức:<br /> p<br /> t new<br />  tp p<br /> <br />  t old<br />  u<br /> t new<br />  tu u<br /> <br />  t old<br />  1 t   p   2 t<br />  1 t   u   2 t<br /> ηt là sai số khống chế định trước Hình 2 Bảy điểm tích phân số<br />  1  (0,1.0) ,  2  1.0 Hình 1 và hình 2 thể hiện phần tử và hình ảnh<br /> Bước thời khống chế phải đồng thời thỏa điểm tích phân dùng trong chương trình . Trong<br /> mãn 2 điều kiện về sai số áp lực lỗ rỗng và sai<br /> đó O thể hiện bậc tự do chuyển vị ▲ thể hiện độ<br /> số ứng suất<br /> 7. Chương trình ứng dụng và tính chính xác bậc do áp lực lỗ rỗng.<br /> của nó<br /> Dùng ngôn ngữ lập trình Fortran, dựa theo 7.1 Ví dụ 1<br /> Dùng chương trình SCRIP 2004 để tính toán<br /> lý thuyết đã trình bày ở trên lập chương trình<br /> ví dụ bài toán cố kết 1 chiều nhằm kiểm tra tính<br /> tính toán SCRIP 2004.<br /> chính xác của chương trình.<br /> Chương trình SCRIP- 2004 có một vài đặc<br /> Ví dụ dùng ở đây được thể hiện trên hình 3 có<br /> điểm chính như sau :<br /> các thông số chính sau: µ=0.301, E=3 Mpa,<br /> - Phạm vi tính toán: có thoát nước và không<br /> kv=10-6 cm/s , H= 10m ,qo=100kPa, t0=70, mặt<br /> thoát nước , bài toán cố kết Biot 2 chiều hoặc<br /> trên thoát nước và mặt dưới không thoát nước.<br /> đối xứng trục.<br /> - Vật liệu vùng tính toán: đàn hồi đồng nhất<br /> hoặc không đồng nhất đẳng hướng , trạng thái<br /> giới hạn của đất (Cam-clay, modified Cam-clay)<br /> [9]<br /> <br /> - Loại phần tử: phần tử tam giác, tứ giác bậc<br /> 1 hoặc bậc 2, bài toán cố kết được tăng thêm độ<br /> tự do áp lực lỗ rỗng.<br /> - Kỹ thuật phi tuyến: Chương trình SCRIP-<br /> 2004 dùng phương pháp tích phân từng bước<br /> giải bài toán phi tuyến (increment approach).<br /> Trong bài toán này chọn phần tử loại tam<br /> giác 6 điểm nút với 7 điểm tích phân. Đối với<br /> bài toán 2 chiều, phần tử nêu trên có 15 bậc tự<br /> do trong đó 12 bậc tự do chuyển vị và 3 bậc tự Hình 3 Ví dụ điển hình bài toán cố kết 1 chiều<br /> do áp lực lỗ rỗng.<br /> <br /> <br /> <br /> <br /> Hình 1 Phần tử cố kết 6 nút dùng trong chương<br /> trình SCRIP 2004 Hình 4 Lưới phần tử bài toán cố kết 1 chiều<br /> Lưới phần tử được thể hiện trong hình 4 bao gồm<br /> 14 phần tử, 45 nút, tổng cộng có 106 bậc tự do.<br /> Như thể hiện trong hình 1 , điều kiện biên bao quả tính toán bằng chương trình được lấy từ nút 8<br /> của lưới phần tử( nút trên cùng phía bên trái của<br /> gồm : 2 mặt trái , phải hạn chế chuyển vị theo<br /> lưới)<br /> chiều y, mặt dưới có hạn chế chuyển vị theo Ví dụ này đã chứng minh được tính chính xác<br /> của chương trình SCRIP 2004.<br /> chiều z, không thoát nước.<br /> <br /> (1)<br /> (2)<br /> 7.2 Ví dụ 2<br /> Giải bài toán tính toán nền bán vô hạn đối<br /> xứng chịu tác dụng của tải trọng phân bố đều.<br /> Các thông số đưa vào tính toán:<br /> Ex=Ey=3Mpa;µ=0.25;H=10m;qo=30kPa;<br /> kx=ky=10-6cm/s;T=694 ngày<br /> <br /> <br /> <br /> <br /> Áp lực lỗ rỗng<br /> Hình 5<br /> So sánh kết quả tính toán áp lực lỗ rỗng<br /> <br /> Hình 7 Kết quả tính toán lún tại điểm A<br /> <br /> Tiến hành tính toán phân tích với nhiều giá trị<br /> (1) sai số khống chế cho phép η để phân tích hiệu quả<br /> (2)<br /> của bước thời gian tự thích ứng. Từ kết quả trên<br /> hình 8 cho thấy toàn bộ quá trình lún tại điểm A<br /> với các giá trị sai số η khác nhau và hoàn toàn có<br /> thể dựa vào η để khống chế độ chính xác của tính<br /> toán. Với η=5% kết quả giải bằng phương pháp số<br /> đã rất sát với kết quả giải bằng phương pháp lý<br /> luận.<br /> 8. Kết luận<br /> Hình 6 Bài báo này dựa trên lý thuyết cơ bản cố kết<br /> So sánh kết quả tính lún và độ cố kết Biot, nhằm nâng cao độ chính xác của tính toán,<br /> bài báo đã tiến hành nghiên cứu bước thời gian tự<br /> Trong hai hình 5 và 6 : thích ứng trong quá trình cố kết. Dùng những<br /> : Lời giải lý thuyết nguyên lý cơ bản để đưa ra công thức khống chế<br /> : Lời giải bằng phương pháp số sai số , từ đó có thể khống chế toàn bộ sai số của<br /> trường ứng suất và trường áp lực lỗ rỗng. Kết quả<br /> Kết quả tính toán so sánh áp lực lỗ rỗng giái nghiên cứu này nhằm chính xác hóa quá trình phát<br /> theo phương pháp lý thuyết và phương pháp số triển ứng suất biến dạng. Những kết quả đạt được<br /> được thực hiện bằng chương trình SCRIP 2004 hết sức có ý nghĩa trong quá trình thiết kế và thi<br /> thể hiện trên hình 5 tại các thời điểm t=14, 42, công công trình vật liệu địa phương.<br /> 70, 140 ngày tương ứng với các độ cố kết Tv = Vấn đề lựa chọn giá trị sai số cho phép [η] vẫn<br /> 0.05, 0.15, 0.25, 0.5. Kết quả phân tích trên hình còn nhiều tồn tại . Hiện nay chỉ hoàn toàn dựa vào<br /> 5 và hình 6 cho thấy rõ ràng kết quả tính toán kinh nghiệm của người tính toán để quyết định giá<br /> bằng chương trình SCRIP 2004 với phương trị này.<br /> pháp số hoàn toàn phù hợp với lý thuyết. Kết<br />  Bài báo này dùng 2 sai số phân biệt η=η(u) [7] Chen ShengHong Ứng dụng FEM phân tích dàn dẻo<br /> và η=η(p) để khống chế độ sai số của trường tự thích ứng mái dốc âu thuyền công trình Tam<br /> Hiệp . Báo cơ học đất Trung Quốc.1998,19(1) 13-19<br /> ứng suất và trường áp lực lỗ rỗng là chưa hoàn (Bản tiếng Trung Quốc)<br /> thiện. Trong thực tế 2 trường này có những tác [8] Deng JianHui Xiong WenLin Phương pháp tự động<br /> động, quan hệ tương hỗ lẫn nhau. Do vậy chia lưới phần tử tam giác tự thích ứng trong miền<br /> nghiên cứu tiếp theo sẽ kết hợp 2 sai số này và tính toán phức tạp. Báo cơ học đất Trung Quốc<br /> đưa ra được một hàm sai số η=η(u,p) duy nhất . 1994/15/2 43-54 (Bản tiếng Trung Quốc)<br /> [9] A.M.Britto & M.J.Gunn. Critical state soil mechanics<br /> via finite elements. University of Bristol.<br /> TÀI LIỆU THAM KHẢO [10] Zienkiewicz O C, Zhu J Z. A simple error estimator<br /> [1] Xie YongLi, Pan QiuYuan. Phương trình vật lý and adaptive procedure for practical enginerring<br /> cố kết biến dạng lớn và hình thức ma trận của nó. analysis[J] . Int.J.Num.Meth.Eng.1897,24(2):337-357.<br /> Tạp chí giao thông học viện giao thông Trùng [11] Biot M.A.General theory of Three Dimensional<br /> Khánh. 1994, 13(2):77-81. (bản tiếng Trung Quốc) Consolidation. J.App.Physis, Vol. 12 ,155. 1941<br /> [2] Gong XiaoNan. Tính toán phân tích côngtình đất. [12] Geudehus, G., Finite Elements in Geomechanics,<br /> Nhà xuất bản xậy dựng Trung quốc, 1999. (Bản Jonhn Wiley&Sons, 1977<br /> tiếng Trung Quốc) [13] Shoji, M.& Matsumoto, T.Consolidation of<br /> [3] Chen ShengHong. Phân tích ổn định Mái dốc cao Embankment Foundation, Soil and Foundations,<br /> và đập vật liệu địa phương. Nhà xuất bản thủy lợi Vol.16, No.1, 59~74, 176<br /> thủy điện Trung quốc.2001. (Bản tiếng Trung [14] Murry, R.T. Development in Two- and Three-<br /> Quốc) Dimentional Consolidation Theory, Development in<br /> [4] Chen JianRong. Lý thuyết ,ứng dụng tính toán lún Soilmechanics, 1978<br /> và cố kết không thoát nước của nền. Báo Kỹ thuật [15] Olson,R.E.Consolidation under Time-Dependent<br /> công trình cảng 1994, 4 : 51-56. (Bản tiếng Trung Loading,J.Geotech.Eng.Div.,ASCE,Vol.103,No.GT1,<br /> Quốc) 55~60, 1977<br /> [5] Chen ShengHong Wang JingSong Zhang JunLu.<br /> Ứng dụng FEM phân tích đàn dẻo tự thích ứng [16] Nguyễn Quang Hùng, Chen ShengHong. Phân Tích<br /> trong công trình thủy công. Báo thủy lợi Trung động lực. Tạp chí Người Xây Dựng, ISSN 0866 8531,<br /> quốc.1996(2) 68-75. (Bản tiếng Trung Quốc) Tháng 7/2004<br /> [6] Ku Ke. Lý thuyết phân tích FEM và ứng dụng. Xuất<br /> bản lần thứ 2. Bắc kinh Nhà xuất bản khoa học.<br /> 1989. (Bản tiếng Trung Quốc)<br /> <br /> <br /> <br /> SUMMARY<br /> Study on adaptive time step of consolidation geotechnical problems<br /> by finite element method<br /> NGUYEN QUANG HUNG<br /> (School of Water Resources and Hydropower, Wuhan University, Wuhan 430072, China)<br /> <br /> The accuracy and stabilility of a problem solved by FE method is concerned and taken in to account by<br /> many researchers. In this paper in order to improve the numerical accuracy and stabilility, based on Biot’s<br /> consolidation theory the adapted time step algorithm is introduced for solving consolidation problems. The method<br /> is to estimate local and global erros of stress state and pore pressure. By the link between errors and time step<br /> length the time step adaptivity strategy has been proposed.<br /> This study can provide well understanding of development of the relation between stress and strain by time.<br /> Furthermore the study will be very useful for the design and construction of dams with earth core wall.<br /> Keywords: consolidation; finite element method; adaptive time step; error norm.<br />