Mô phỏng dòng chảy trong sông bằng sóng động học một chiều phi tuyến

Nội dung text: Mô phỏng dòng chảy trong sông bằng sóng động học một chiều phi tuyến

1. Tạp chí Khoa học ĐHQGHN: Các Khoa học Trái đất và Môi trường, Tập 32, Số 3S (2016) 14-19 Mô phỏng dòng chảy trong sông bằng sóng động học một chiều phi tuyến Bùi Văn Chanh1,*, Trần Ngọc Anh2,3, Lương Tuấn Anh4 1Đài Khí tượng Thủy văn khu vực Nam Trung Bộ, Trung tâm KTTV Quốc gia, Bộ TNMT, 22 Pasteur, Nha Trang, Khánh Hòa, Việt Nam 2Khoa Khí tượng Thủy văn và Hải dương học, Trường Đại học Khoa học Tự nhiên, ĐHQGHN, 334 Nguyễn Trãi, Hà Nội, Việt Nam 3Trung tâm Động lực học Thủy khí Môi trường, Trường Đại học Khoa học Tự nhiên, ĐHQGHN, 334 Nguyễn Trãi, Hà Nội, Việt Nam 4Viện Khoa học Khí tượng Thủy văn và Biến đổi Khí hậu, Bộ Tài nguyên và Môi trường, 62 Nguyễn Chí Thanh, Đống Đa, Hà Nội, Việt Nam Nhận ngày 08 tháng 8 năm 2016 Chỉnh sửa ngày 26 tháng 8 năm 2016; Chấp nhận đăng ngày 16 tháng 12 năm 2016 Tóm tắt: Mô phỏng dòng chảy thượng nguồn các con sông là rất quan trọng và cần thiết, do hạn chế về số liệu nên việc mô phỏng gặp nhiều khó khăn. Trong nghiên cứu này trình bày phương pháp mô phỏng dòng chảy phân bố bằng mô hình sóng động học phi tuyến, vừa giải quyết hạn chế vấn đề số liệu vừa đáp ứng yêu cầu mô phỏng và cho kết quả nhanh hơn. Mô hình sóng động học phi tuyến được xây dựng từ hệ phương trình Saint Venant, trong đó gồm một chương trình sóng động học phi tuyến giải hệ phương trình bằng phương pháp lặp Newton và một chương trình sóng động học tuyến tính phục vụ tính toán giá trị lưu lượng ban đầu. Chương trình sóng động học tuyến tính được xây dựng từ hệ phương trình Saint Venant và được giải bằng sơ đồ sai phân ẩn 4 điểm. Chương trình sóng động học tuyến tính sau khi lập trình cho kết quả trùng khớp với kết quả tính toán trong giáo trình Thủy văn ứng dụng của Vante Chow [1, 1988]. Mô hình sóng động học phi tuyến gồm 2 phần, trong đó phần đầu là chương trình sóng động học phi tuyến, phần sau là chương trình sóng động học phi tuyến. Trong nghiên cứu đã xây dựng phương pháp, sơ đồ giải, lập trình chương trình tính toán cho mô hình. Mô hình được lập trình bằng ngôn ngữ lâp trình Fortran 90 và kiểm tra chất lượng mô phỏng tại trạm thủy văn Tà Pao, Võ Xu trên sông La Ngà tỉnh Bình Thuận. Kết quả mô phỏng của mô hình khá tốt, tuy nhiên mô hình có nhược điểm là không mô phỏng được cho đoạn sông có ảnh hưởng triều, nước vật. Từ khóa: Sóng động học, Saint Venant, phương pháp lặp Newton. 1. Mở đầu * bố một chiều. Các phương trình liên tục và động lượng bảo toàn và không bảo toàn bỏ qua Phương trình Saint Venant có nhiều dạng dòng bên, lực cản của gió và các tổn thất xoáy giản hóa khác nhau, mỗi dạng xác định một mô được dùng để định nghĩa các loại mô hình khác hình diễn toán dòng chảy không ổn định phân nhau về diễn toán dòng chảy không ổn định phân bố một chiều. ___ Phương trình động lượng bao gồm các * Tác giả liên hệ. ĐT.: 84-915620289 Email: buivanchanh@gmail.com thành phần thuộc các quá trình vật lý điều khiển 14
2. B.V. Chanh và nnk. / Tạp chí Khoa học ĐHQGHN: Các Khoa học Trái đất và Môi trường, Tập 32, Số 3S (2016) 14-19 15 dòng động lượng. Các thành phần này là: thành A = αQβ (7) phần gia tốc địa phương mô tả sự thay đổi của Trong phương trình Manning với So = Sf và động lượng do thay đổi của vận tốc trong thời R = A/P ta có: gian, thành phần gia tốc đối lưu mô tả sự thay 1 đổi của động lượng gây ra bởi sự thay đổi của 2 5 1.49S0 3 vận tốc dọc theo kênh, thành phần áp lực tỉ lệ QA 2 (8) với sự thay đổi độ sâu của nước dọc theo kênh, nP 3 thành phần trọng lực tỉ lệ với độ dốc đáy S và 0 Viết lại phương trình (8) cho A từ đó tìm thành phần ma sát tỉ lệ với độ dốc ma sát Sf. Thành phần gia tốc địa phương và gia tốc đối lưu được α và β = 0.6 như sau: đại biểu cho tác động của các lực quán tính lên 3 2 5 dòng chảy. nP 3 3 AQ 5 (9) + Phương trình liên tục: 1.49 S0 QA  0 (1) 0.6 x  t 2 + Phương trình động lượng: nP 3 A (10) 2 1Q 1  Q  y 1.49 S0 g g S0 S f 0 (3) A t A  x A  x Sóng động học chi phối dòng chảy khi các Phương trình (1) chỉ phụ thuộc vào A và Q, lực quán tính và áp lực có thể bỏ qua, trong trong đó A được xác định trong phương trình sóng động học lực ma sát và trọng lực cân bằng (7). Đạo hàm riêng phương trình (7) của biến A nhau nên dòng nước chảy không có gia tốc. Đối và Q theo t và thế vào phương trình (1) được với sóng động học, đường năng song song với phương trình (11). Thế phương trình (11) vào đáy kênh và dòng chảy trong đoạn nguyên tố là phương trình (5) được phương trình (12). một dòng đều ổn định (vì S0=Sf). Phương trình (12) được sai phân theo sơ đồ Sóng động học tạo nên do sự thay đổi trong tuyến tính theo phương trình (19), sai phân theo dòng chảy như thay đổi về lưu lượng nước hoặc sơ đồ phi tuyến theo phương trình (24). tốc độ sóng là vận tốc truyền thay đổi dọc theo AQ 1  kênh dẫn. Tốc độ sóng phụ thuộc vào loại sóng Q (11) đang xét và có thể hoàn toàn khác biệt với vận t  t tốc dòng nước. Đối với sóng động học, các QQ 1  thành phần gia tốc và áp suất trong phương Q q (12) trình động lượng đã bị bỏ qua nên chuyển động x  t của sóng được mô tả chủ yếu bằng phương trình liên tục. Do đó sóng đã mang tên sóng động học và động học nghiên cứu chuyển động 2. Xây dựng mô hình trong đó không xét đến ảnh hưởng của khối lượng và lực. Mô hình sóng động học được xác Mô hình được xây dựng trên ngôn ngữ lập định bằng các phương trình như sau: trình Fortran 90 gồm hai phần chính là mô hình - Phương trình liên tục: sóng động học tuyến tính và phi tuyến. Trong QA  đó mô hình tuyến tính được sử dụng để làm q (5) nghiệm thứ nhất của mô hình phi tuyến. Mô x  t hình tuyến tính được giải bằng sơ đồ sai phân - Phương trình động lượng: ẩn, mô hình phi tuyến được giải bằng phương So = Sf (6) pháp lặp Newton.
3. 16 B.V. Chanh và nnk. / Tạp chí Khoa học ĐHQGHN: Các Khoa học Trái đất và Môi trường, Tập 32, Số 3S (2016) 14-19 2.1. Sơ đồ sóng động học tuyến tính Áp dụng sơ đồ sai phân ẩn: uj 1 u j 1 u j 1 i 1 i 1 i (13) x x uj 1 u j 1 u j i 1 i 1 i 1 (14) t t QQQj 1 j 1 j 1 i 1 i 1 i (15) x x QQQj 1 j 1 j i 1 i 1 i 1 (16) t t QQj 1 j Q i i 1 (17) 2 qj 1 q j q i 1 i 1 (18) 2 Thay thế các phương trình từ (15) đến (18) vào phương trình (12) được phương trình sai Hình 2. Sơ đồ khối tính toán phân sóng động học tuyến tính cho sơ đồ ẩn sóng động học tuyến tính. như sau:  1 Qj 1 Q j 1 Q j Q j 1 Q j 1 Q j q j 1 q j 2.2. Sơ đồ sóng động học phi tuyến i 1 i  i 1 i i 1 i 1 i 1 i 1 (19) x 2 t 2 Sai phân phương trình (12) như sau: j i 1  1 j 1 j 1 j 1 j j 1 j t QQ qj 1 q j Qj 1  Q j i 1 j t 2 i 1 i 1 Qi 1 Q i A i 1 A i 1 q i 1 q i 1 i i 1 (21) x 2 j 1 x t 2 Qi 1 (20) j j 1  1 t QQi 1 i   j 1 j 1 x 2 AQi 1 i 1 (22) g j j  AQi 1 i 1 (23) Thế phương trình (22) và (23) vào (21) ta được: j 1 j   tj 1 j 1 t j 1 j qi 1 q i 1 Qi 1 Q i 1 Q i Q i 1 t (24) x x 2 Phương trình này đã được sắp xếp cho lưu j+1 lượng chưa biết Qi+1 nằm ở vế trái và các đại lượng đã biết nằm ở vế phải. Đây là phương j 1 trình phi tuyến đối với Qi 1 do đó cần được giải bằng phương pháp số, trong chương trình lập trình và sơ đồ khối dưới đây áp dụng phương pháp lặp Newton. Mô hình tuyến tính Hình 1. Sơ đồ sai phân ẩn giải phương trình sóng. xây dựng trong mô hình phi tuyến được thể hiện động học tuyến tính. trong khối ước lượng ban đầu bằng cách sử dụng ước lượng tuyến tính 20 như hình dưới đây.
4. B.V. Chanh và nnk. / Tạp chí Khoa học ĐHQGHN: Các Khoa học Trái đất và Môi trường, Tập 32, Số 3S (2016) 14-19 17 j 1 j j 1 t  q q j 1 j i 1 i 1 Ước lượng giá trị khởi đầu của Qi 1 trong C Qi Q i 1 t (25) x 2 mỗi quá trình lặp có ảnh hưởng quan trọng đến sự hội tụ của sơ đồ. Một cách tiếp cận là sử Từ đó một sai số dư f(Qj 1 ) được xác định i 1 dụng nghiệm của sơ đồ tuyến tính, phương trình bằng phương trình (24). (20) như là nghiệm gần đúng thứ nhất của sơ đồ j 1 t j 1 j 1  phi tuyến. Li Simons và Stevens (1975) [1] sau f( Qi 1 ) Q i 1 ( Q i 1 ) C (26) khi tiến hành các phân tích về tính ổn định đã x chỉ ra sơ đồ sử dụng phương trình (24) là một j+1 Đạo hàm bậc nhất của f(Qi+1 ) như sau: sơ đồ ổn định không điều kiện và có thể sử t dụng các trị của Δt/Δx trong một phạm vi khá f'( Qj 1 )  ( Q j 1 ) 1 (27) i 1 x i 1 rộng mà không tạo ra sai số lớn trong hình dạng của đường quá trình lưu lượng. j 1 j 1 Mục tiêu là tìm Qi 1 để buộc f(Qi 1 ) bằng không. Sử dụng phương pháp lặp Newton và các bước lặp k = 1, 2, 3, 3. Mô phỏng thử nghiệm mô hình f() Qj 1 k j 1 j 1 i 1 Mô hình sau khi lập trình được kiểm tra với (QQi 1 ) k 1 ( i 1 ) k 'j 1 (28) f() Qi 1 k ví dụ 9.6.1 trong giáo trình Thủy văn Ứng dụng Tiêu chuẩn hội tụ cho quá trình lặp là: của Vente Chow [1]. Kết quả của mô hình trùng khớp với kết quả tính toán của ví dụ trên, tương j 1) f( Qi 1 ) k 1  (29) quan kết quả tính được thể hiện trong hình 4. Mô hình sóng động học một chiều phi tuyến sau khi xây dựng được mô phỏng thử nghiệm mô phỏng dòng chảy trên sông La Ngà đoạn từ đập thủy điện Đa Mi đến trạm thủy văn Phú Hiệp. Điều kiện ban đầu xác định như sau: độ rộng trung bình của sông là 95m, độ dốc sông trung bình 1%, chiều dài sông 115310m, số đoạn sông là 8 đoạn, bước thời gian mô phỏng là 360 phút, lưu lượng ban đầu 40m3/s, hệ số nhám Manning được xác định từ bảng tra thủy lực của M.F. Xripnut cho sông La Ngà 0.032. Hình 3. Sơ đồ khối tính toán Hình 4. Tương quan kết quả tính bằng chương trình sóng động học phi tuyến. và của Vente Chow.
5. 18 B.V. Chanh và nnk. / Tạp chí Khoa học ĐHQGHN: Các Khoa học Trái đất và Môi trường, Tập 32, Số 3S (2016) 14-19 vực sông La Ngà của Đài Khí tượng Thủy văn khu vực Nam Trung Bộ. Mô hình được thiết lập cho khu vực hạ lưu hồ thủy điện Đa Mi với các trạm mưa Tà Pao, Võ Xu, Suối Kiết, Đông Giang, Mê Pu, bốc hơi tính theo trạm khí tượng Phan Thiết. Bộ thông số MIKE NAM được hiệu chỉnh và kiểm định tại trạm thủy văn Đại Nga trên sông La Ngà, đánh giá chất lượng hiệu chỉnh theo chỉ tiêu Nash là 90%, kiểm định là 87%, đạt loại tốt theo tiêu chí của Tổ chức Khí Hình 5. Đường tính toán và thực đo trạm Tà Pao. tượng Thế giới (WMO). Từ số liệu thực đo tại tạm Tà Pao và Phú Hiệp tiến hành hiệu chỉnh theo phương pháp thử sai đã tìm được bộ thông số hiệu chỉnh với độ dốc sông trung bình là 1,2%, hệ số nhám là 0.0357. Kết quả hiệu chỉnh tại trạm thủy văn Tà Pao theo chỉ tiêu Nash là 93,5% đạt loại tốt theo tiêu chí của WMO; trạm thủy văn Võ Xu là 82,0%, đạt loại khá theo tiêu chí của WMO. Tài liệu tham khảo Hình 6. Đường tính toán và thực đo trạm Võ Xu. [1] Ven Techow, David R.Maidment, Larry Mô hình mô phỏng dòng chảy từ dữ liệu xả W.Mays (1988), Applied Hydrology, New York, hồ Đa Mi và gia nhập khu giữa được tính từ mô McGraw-Hill, 1988. hình MIKE NAM của Đan Mạch với thời đoạn [2] Nguyễn Hữu Khải, Nguyễn Thanh Sơn, Mô 6 giờ, thời gian mô phỏng từ 01h/01/6 đến hình toán thủy văn, Nxb Đại học Quốc gia Hà Nội, 2003. 07h/14/8 năm 2016. Mô hình MIKE NAM được kế thừa từ phương án dự báo thủy văn lưu Somulation River Discharge by Non-linear One Dimension Kinematic Wave Bui Van Chanh1, Tran Ngoc Anh2,3, Luong Tuan Anh4 1South Center Regional Hydro - Meteorologial Center, NHMS, MONRE 2Faculty of Hydro-Meteorology & Oceanography, VNU University of Science, 334 Nguyen Trai, Hanoi, Vietnam 3Center for Environmental Fluid Dynamic, VNU University of Science, 334 Nguyen Trai, Hanoi, Vietnam 4Vietnam Institute of Meteorology Hydrology and Climate Change, MONRE, 62 Nguyen Tri Thanh, Dong Da, Hanoi, Vietnam Abstract: Simulating discharge on upstream of rivers is very important and necessary, by restrict about data so the simulating is very difficult. This research that is presented distributed discharge
6. B.V. Chanh và nnk. / Tạp chí Khoa học ĐHQGHN: Các Khoa học Trái đất và Môi trường, Tập 32, Số 3S (2016) 14-19 19 simulating method by by non-linear one dimension kinematic wave model. This model overcome restrict by lack data and satisfy demand by simulating, calculate more quickly. Non-linear one dimension kinematic wave model is buit by Saint Venant system equation that include non-linear one dimension kinematic wave program and linear one dimension kinematic wave program. Saint Venant system equation in the non-linear and linear program are solved correlatively by Newton iteraction method and 4 points latent defference schedule. The linear program is used to calculate initial dischare for the non-linear program and result of the linear program is same result in Applied Hydrology syllabus of Vante Chow [1, 1988]. The linear program is first section and the non- linear program is second section in the model. This research buit method, solving schedule, programming calculating for the model. Fortran 90 programming languge is use to buid the model and checked simulating quality at Ta Pao, Vo Xu hydrology stations on La Nga river Binh Thuan province. Simulating result of the model is good, but the model can’t simulate on river section that is effected by tide and backwater and that is restrict of the model. Keywords: Kinematic wave, Saint Venant, Newton iteraction method.