Bài giảng Mô hình nước ngầm - Phần 1: Cơ sở lý thuyết - Nguyễn Mai Đăng

Nội dung text: Bài giảng Mô hình nước ngầm - Phần 1: Cơ sở lý thuyết - Nguyễn Mai Đăng

1. 10/6/2013 MÔ HÌNH NƯỚC NGẦM Phần 1: Cơ sở lý thuyết Nggyuyễn Mai Đăng Bộ môn Thủy văn & Tài nguyên nước dang@wru.vn 0989.551.699 Quá trình ứng dụng mô hình Xác định và mô tả vấn đề 1. Xác định vấn đề cầnmôphỏng – Các yếutố cầnmôhìnhhóa – Mối quan hệ và tương tác giữa chúng Số liệu Cơ sở lý thuyết mô hình – Độ chính xác 2. Xác định cơ sở lý thuyết và phát triểnmôhình – Mô tả về mặttoánhọc Phát triển mô hình – Loạimôhình – Phương pháp số –lậptrìnhtrênmáytính Hiệu chỉnh mô hình & – Lướitínhtoán, điềukiệnbiên, đ/k ban đầu Đánh giá thông số 3. Hiệuchỉnh thông số mô hình – Tạmthờixácđịnh các thông số củamôhình Kiểm định mô hình & – So sánh kếtquả tính toán từ mô hình vớisố liệuthực đo Phân tích độ nhạy – Hiệuchỉnh các thông số cho đếnkhikếtquả mô hình đạt được độ chính xác đặt ra. 4. Kiểm định Tổng hợp và hoàn – Thiếtlậpbộ số liệu độclập(so vớigiaiđoạnhiệuchỉnh) chỉnh mô hình – Chạymôhình(giữ nguyên bộ thông số xác định đượctừ giai đoạnhiệuchỉnh) – So sánh, đánh giá sự sai khác giữakếtquả tính toán và thực Ứng dụng mô hình đo Trình diễn các kết quả 1
2. 10/6/2013 Các công cụ để giải quyết các vấn đề của nước ngầm • Các phương pháp tương tự và vậtlý – Some of the first methods used. • Các phương pháp phân tích – Những gì chúng ta đã thảo luận cho đến thời điểm này? www.epa.state.oh.us – Gặpkhókhănkhi: các biên không đều, điềukiện biên khác nhau, tính chất không đồng nhất, không đẳng hướng, nhiềugiai đoạn, phi tuyến • Các phương pháp số – Chuyển các phương trình vi phân từng phần(đạohàmriêng) củanướcngầm www.isws.illinois.edu sang hệ phương trình vi phân thường (đạo hàm toàn phần) hoặc chuyển sang các phương trình đạisốđểgiảiranghiệm bài toán. Mô hình khái niệm • Mô tả đại diện của hệ thống nước ngầm kết hợp với giải thích các điều kiện địa chất và thủy văn. • Quá trình nào là quan trọng đối với mô hình? • Các biên là gì? • Các giá trị thông số có tồn tại không? • Các giá trị thông số phải được thu thập là gì? 2
3. 10/6/2013 Chúng ta thực sự muốn tính toán cái gì? • Dòng chảy hướng ngang trong tầng ngậm nước có áp và thấm yếu ∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ K' W ∂h ⎜Tx ⎟ + ⎜Ty ⎟ + (h0 − h) ± ∑Qwδ (x − xw ) = S ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ b' w=1 ∂t Dòng chảyThấmNguồn chảy Lượng trữ Bề mặt đất vào/chảy ra • Các phương trình chủđạủ đạo Cộtnt nướctc tầng có áp • Các điều kiện biên Lớp có áp (không thấm) • Các điều kiện ban đầu Qx Tầng có áp h b z y K x Tầng đá gốc Phương pháp sai phân hữu hạn (Finite Difference Method) • Phương pháp sai phân hữu hạn – Thay thế các đạo hàm trong phương trình nước ngầm bằng các ‘xấp xỉ chuỗi Taylor’ – Thiết lập hệ phương trình đại số để giải tìm nghiệm 3
4. 10/6/2013 Chuỗi Taylor • Khai triển chuỗi Taylor đối với h(x) tại điểm (x+∆x) lân cận điểm x (∆x)2 (∆x)n h(x+∆x) = h(x)+∆xh′(x)+ h′′(x)+ + hn (x) + 2! L n! L • Nếu cắt chuỗi số từ các trị số n+1 trở đi, thì sai số của chuỗi số là: (∆x) n+1 error < h n+1 (x) = ϑ (∆x n +1 ) (n+1)! max Đạo hàm bậc 1 “tiến” (+∆x) • Phát triển chuỗi Taylor “tiến” cho hàm mực nước h(x) tại vị trí lân cận điểm x là: (∆x)2 h(x+∆x) = h(x)+∆x ⋅ h′(x)+ h′′(x)+ 2! L Trong đó đạo hàm bậc 1 được tính xấp xỉ như sau: ∂h h(x+∆x)-h(x) h′(x) = ≈ ∂x ∆x ∆x ∆x x x − ∆x x x + ∆x 4
5. 10/6/2013 Đạo hàm bậc 1 “lùi” (‐∆x) • Phát triển chuỗi Taylor “lùi” cho hàm mực nước h(x) tại vị trí lân cận điểm x là: (∆x)2 h(x-∆x) = h(x)-∆x ⋅ h′(x)+ h′′(x)- 2! L Trong đó đạo hàm bậc 1 được tính xấp xỉ như sau: ∂h h(x)-h(x-∆x) h′(x) = ≈ ∂x ∆x ∆x ∆x x x − ∆x x x + ∆x Đạo hàm bậc 2 (∆x)2 h(x+∆x) = h(x)+∆x ⋅ h′(x)+ h′′(x)+ 2! L (∆x)2 h(x+∆x) = h(x) − ∆x ⋅ h′(x)+ h′′(x) − 2! L Trong đó h”(x) là đạo hàm bậc 2 được tính xấp xỉ là: ∂2h h(x+∆x)-2h(x)+h(x-∆x) h′′(x) = ≈ ∂x2 ∆x2 5
6. 10/6/2013 Xấp xỉ sai phân hữu hạn h(x) = hi h(x + ∆x) = hi+1 h(x − ∆x) = hi−1 ∆x ∆x h′(x) = hi′ x h -h hi -h i −1 i +1 i hi+1-2hi+hi−1 h′ = hi′ = h′′= i ∆x ∆x i ∆x 2 Đạo hàm Đạo hàm Đạo hàm lùi bậc 1 tiến bậc 1 trung tâm bậc 2 Lưới và rời rạc (Grids and Discretization) • Quá trình rời rạc y, j • Lưới được thiết lập để phủ kín Lưới vùng tính toán • Mục tiêu –dự báo mực nước tại Vùng tính toán các điểm nút của lưới i,j+1 – Xác định ảnh hưởng của bơm ∆ y i-1,j i,j i+1,j – Dòng chảy từ sông, • Phương pháp sai phân hữu hạn i,j-1 Điểm nút – Được sử dụng phổ biến do tính đơn giản của phương pháp, dễ hiểu, dễ sử dụng và lại đạt kết x, i quả tốt ∆x – Sử dụng tốt cho vùng địa hình đơn giản Ô lưới 6
7. 10/6/2013 Lưới không gian 3 chiều (Three‐Dimensional Grids) • Một hệ thống tầng ngậm nước được chia thành các khối chữ nhật bởi các lưới. • Lưới được tổ chức bởi các hàng (i), cột (j), lớp (k), mỗi một khối (block) được gọi là một “ô lưới“ (cell) • Các loại lớp địa tầng – Có áp – Không áp – Chuyển đổi Có thể mô phỏng cho các lớp có vật liệu khác nhau Dòng chảy tầng ngậm nước có áp 1 chiều (1‐D Confined Aquifer Flow) • Dòng chảy đồng nhất, đẳng hướng, 1 chiều, và Bề mặt đất (Ground surface) có áp. • Phương trình mô tả cho Lớp không thấm h dòng chảy xác định ở trên: A Nút (Node) Tầng ngậm nước (Aquifer) ∆x ∂ ⎛ ∂h ⎞ ∂h hB ⎜T ⎟ = S ∂x ∂x ∂t b ⎝ ⎠ z y x • Điều kiện ban đầu (initial i = 0 1 2 3 4 5 6 7 8 9 10 condition ): Ô lưới (Grid Cell) h ( x ,0 ) = 6 .1 m • Điều kiện biên (boundary ∆x= 1 m, L = 10 m, b = 1.5 m conditions): hA = 6.1 m, hB = 1.5 m, h(0, t) = 6.1 m K = 0.5 m/d, S = 0.02 h(L, t) = 1.5 m 7
8. 10/6/2013 Tính toán xấp xỉ cho các đạo hàm (Derivative Approximations) • Phương trình nước ngầm 1 chiều, đồng nhất, đẳng hướng, có áp: ∂2h S ∂h t, l = ∂x2 T ∂t i, l+1 • Đạo hàm bậc 2 đối với (WRT with respect ∆t i−1, l i+1,l to) x i, l 2 l l l l x, i ∂ h hi+1-2hi +hi−1 ∆x 2 ≈ 2 ∂x i ∆x i, l−1 • Đạo hàm bậc 1 đối với (WRT) t l l +1 l l l l −1 ∂h hi -hi ∂h h -h ≈ ≈ i i ∂t ∆t i ∂t i ∆l Tiếnlùi Phương pháp sai phân sơ đồ hiện (Explicit Method) • Sử dụng thông tin (H, Q) tại các nút ở bước thời gian trước để tính cho nút ở t, l bước thời gian hiện tại (sơ đồ hiện). i, l +1 • Tính toán cho từng nút một (tất cả các điểm không gian và thời gian) cho đến ∆t i −1, l i, l i +1,l hết miền tính toán. x, i • Phải cẩn thận khi chọn bước thời gian ∆x (∆t), nếu lớn quá sẽ không ổn định trong tính toán, nếu nhỏ quá thì tính toán lâu, i, l −1 do vậy chọn vừa đủ nhỏỏể để bài toán ổn định và cũng đáp ứng được sai số cho phép. 2 l l l l+1 l ∂ h S ∂h h − 2h + h S h − h = i-1 i i+1 = i i ∂x2 T ∂t ∆x2 T ∆t Từ pt vi phân chuyển sang sai phân (xấp xỉ) 8
9. 10/6/2013 Phương pháp sai phân sơ đồ hiện (Explicit Method) hl − 2hl + hl S hl+1 − hl i-1 i i+1 = i i ∆x2 T ∆t T ∆t (h l − 2h l + h l ) = h l +1 − h l S ∆x 2 i-1 i i+1 i i T ∆t Đặt: r = S ∆x2 l+1 l l l l hi = hi +r(hi-1 −2hi +hi+1) Bước thời gian l+1 Bước thời gian l Chưa biết Đã biết Phương pháp sai phân sơ đồ hiện (Explicit Method) Bề mặt đất (Ground surface) l +1 l l l l hi = hi + r (hi-1 − 2hi + hi+1 ) Lớp không thấm hA Nút tính toán Tầng ngậm nước ∆x T ∆t h r = B S ∆x2 b S ∆t = r∆x 2 T i= 0 1 2 3 4 5 6 7 8 9 10 Xem xét: ÔlÔ lưới S r = 0.48 ∆t = r∆x 2 = 0.0128 d ≈ 18.5 min T r = 0.52 ∆ t = 0.0139 d ≈ 20 min ∆x= 1m, L = 10m, b = 1.5m hA = 6.1m, hB = 1.5m, K = 0.5m/d, S = 0.02 9
10. 10/6/2013 Kết quả tính theo sơ đồ hiện (∆t = 18.5 min; r = 0.48 0.5) Kết quả đường mực nước ko ổn định, do ∆t lớn 10
11. 10/6/2013 Phân tích quá trình không ổn định đó diễn ra như thế nào? • Tại thời gian t = 0 Æ ko có dòng Bề mặt đất chảy • Khi t > 0 Æ có dòng chảy • Lượng nước có thể thoát ra từ Lớp ko thấm h lượng trữ trong một ô lưới trong A khoảng thời gian ∆t: Tầng ngậm nước ∆x h ∆ V = S ⋅ ∆ x (1) ⋅ ∆ h B 1 b • Nước chảy ra khỏi một ô lưới ∆x trong thời khoảng ∆t là: ∆h i= 0 1 2 i-1 i i+1 8 9 10 ∆V2 = T ∆t ∆x / 2 Ô lưới i • Theo nguyên lý bảo toàn khối lượng: ∆V2 ≤ ∆V1 Vậy nếu r > 0.5 Æ khoảng thời ∆h T ⋅ ⋅∆t ≤ S ⋅∆x ⋅∆h gian qua lớn Æ ô lưới nước ngầm ∆x / 2 không thể chứa đủ nước Æ không T ∆t 1 r = ≤ ổn định S ∆x2 2 Phương pháp sai phân sơ đồ ẩn (Implicit Method) • Sử dụng thông tin (H, Q) ở một nút thời gian trước tính toán t, l cho 3 nút thời gian sau (sơ đồ i −1, l +1 i, l +1 i +1,l +1 ẩn). ∆t • Giải đồng thời cho tất cả các i −1, l i, l i +1,l x, i nút trong miền tính toán cùng ∆x một thời gian (khác với sơ đồ i −1, l −1 i, l −1 i +1,l −1 hiện là giải từng bước một). • PP này có bản chất luôn luôn ổn định. ∂2h S ∂h hl+1 − 2hl+1 + hl+1 S hl+1 − hl = i-1 i i+1 = i i ∂x2 T ∂t ∆x 2 T ∆t Chuyển từ pt vi phân sang pt sai phân hữu hạn 11
12. 10/6/2013 Phương pháp sơ đồ ẩn (Implicit Method) hl+1 − 2hl+1 + hl+1 S hl+1 − hl i-1 i i+1 = i i ∆x 2 T ∆t T ∆t (h l +1 − 2h l +1 + h l +1 )= h l +1 − h l S ∆x 2 i-1 i i+1 i i T ∆t r = S ∆x2 l +1 l +1 l +1 l − rhi-1 + (1 + 2r)hi − rhi+1 = hi Bước thời gian l+1 Bước thời gian l chưa biết Đã biết Dòng chảy ổn định 2 chiều (2‐D Steady‐State Flow) y, j Node No. (1,5) • Phương trình cơ bản Unknown heads (5,4) Known heads ∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ W ∂h (1,5) (2,5) (3,5) (4,5) ⎜Tx ⎟ + ⎜Ty ⎟ ± Qw = S ⎜ ⎟ ∑ (0,4) (1,4) (2,4) ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ w=1 ∂t (3,4) (4,4) (5,4) ∆y (0,3) (1,3) (2,3) (3,3) (4,3) (5,3) • Dòng chảy đồng nhất, đẳng hướng, ko có giếng (0,2) (1,2) (2,2) (3,2) (4,2) (5,2) ∂2h ∂2h + = 0 (0,1) (1,1) (2,1) (3,1) (4,1) ∂x2 ∂y2 (5,1) (1,0) (2,0) (3,0) (4,0) h − 2h + h h − 2h + h x,i i-1,j i,j i+1,j + i,j-1 i,j i,j+1 = 0 ∆x ∆x 2 ∆y 2 • Bình quân theo không gian (bình quân trong ô lưới) h + h + h + h h = i-1,j i+1,j i,j- 1 i,j+ 1 i,j 4 12
13. 10/6/2013 Dòng chảy 2 chiều, ko đồng nhất, ko đẳng hướng (2‐D Heterogeneous Anisotropic Flow) y i-1/2 i+1/2 node (i,j) ∂ ⎛ x ∂h ⎞ ∂ ⎛ y ∂h ⎞ ∆x ⎜T ⎟ + ⎜T ⎟= 0 ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ j+1 cell Qy,j+1/2 j+1/2 Tx và Ty là hệ số chuyển nước theo phương x và y Qx,i-1/2 Qx,i+1/2 ∆y j ⎛ x ∂h ⎞ ⎛ x ∂h ⎞ ⎛ y ∂h ⎞ ⎛ y ∂h ⎞ ⎜T ⎟ − ⎜T ⎟ ⎜T ⎟ − ⎜T ⎟ ⎝ ∂x ⎠ ⎝ ∂x ⎠ ⎝ ∂y ⎠ ⎝ ∂y ⎠ i+1/ 2.j i−1/ 2.j + i.j +1/ 2 i.j −1/ 2 = 0 ∆x ∆y Dòng chảy 2 chiều, ko đồng nhất, ko đẳng hướng (2‐D Heterogeneous Anisotropic Flow) • Hệ số chuyển nước trung bình, điều hòa ⎛ x h i+ 1 ,j − h i,j ⎞ ⎛ x h i,j − h i- 1 ,j ⎞ ⎜ T i+ 1 / 2 ,j ⎟ − ⎜ T i- 1 / 2 ,j ⎟ ⎝ ∆ x ⎠ ⎝ ∆ x ⎠ ∆ x ⎛ h − h ⎞ ⎛ h − h ⎞ ⎜ T y i,j+ 1 i,j ⎟ − ⎜ T y i,j i,j- 1 ⎟ ⎜ i,j+ 1 / 2 ∆ y ⎟ ⎜ i,j- 1 / 2 ∆ y ⎟ + ⎝ ⎠ ⎝ ⎠ = 0 ∆ y x x x Ti+1,jTi,j Ti+1/ 2,j = 2 x x Ti+1,j+Ti,j A i,j h i + 1 ,j + B i,j h i −1 ,j + C i,j h i,j + 1 + D i,j h i,j −1 + E i,j h i,j = 0 13
14. 10/6/2013 Các vấn đề chuyển đổi (Transient Problems) ∂ ⎛ x ∂h ⎞ ∂ ⎛ y ∂h ⎞ ∂h ⎜T ⎟ + ⎜T ⎟ = S ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂t (h l +1 − hl ) A h l +1 + B h l +1 + C h l +1 + D h l +1 + E h l +1 = S i,j i,j i,j i+1,j i,j i-1,j i,j i,j+1 i,j i,j-1 i,j i,j i,j ∆t l +1 l +1 l +1 l +1 l +1 l Ai,j hi+1,j + Bi,j hi-1,j + C i,j hi,j+ 1 + Di,j hi,j- 1 + E i,j′ hi,j = Fi,j hi,j S E ′ = E − i,j i, j i,j ∆t MODFLOW • Là mô hình toán được phát triển bởi Cục địa chất Mỹ (USGS – US Geological Survey) • Sử dụng phương pháp sai phân hữu hạn • Một số Version đã phát triển tham khảo: • Gia diện đồ họa MODFLOW: – GWV ‐ Groundwater Vitas: ‐ vistas. com/ : – GMS – Ground Modeling System: – PMWIN –Modflow for Window: 14
15. 10/6/2013 MODFLOW có thể mô phỏng được cái gì? Mô hình MODFLOW có thể mô phỏng được những yếu tố sau: 1. Tầng ngậm nước không áp và có áp 2. Đứt gãy và các biên không thấm 3. Các địa tầng có áp hạn mịn và các lớp đáy liên kết 4. Đại tầng có áp ‐ dòng chảy ngầm và sự thay đổi lượng trữ 5. Sông –trao đổi nước ở các tầng ngậm nước 6. Lưu lượng nước từ các kênh tiêu và suối 7. Dòng chảy phù du –trao đổi nước ở các tầng ngậm nước 8. Hồ chứa –trao đổi nước ở các tầng ngậm nước 9. Bổ cập nước ngầm từ mưa và tưới 10. Bốc thoát hơi nước 11. Các giếng bơm hút và bổ cập nước ngầm 12. Xâm nhập mặn 15