Áp dụng các mô hình thủy lực để mô phỏng và dự báo lũ trên hệ thống sông Hồng

Áp dụng các mô hình thủy lực để mô phỏng và dự báo lũ trên hệ thống sông Hồng. Nguyễn Tuấn Anh, Hoàng Văn Lai, Nguyễn Hồng Phong Viện Cơ học, Viện Khoa học và Công nghệ Việt Nam TÓM TẮT: Mô hình thủy lực đã và đang được áp dụng để mô phỏng và dự báo lũ trên các hệ thống sông. Trong các phần mềm thủy lực hiện đại, thí dụ như MIKE và HEC, thường sử dụng các hệ phương trình mô tả: sóng động lực học (dynamic wave model), mô hình sóng khuếch tán (diffusive wave model), mô hình sóng động học (kinematic wave model). Sóng khuếch tán và sóng động học là các mô hình được đơn giản hóa từ mô hình sóng động lực học và phù hợp hơn trong trường hợp số liệu không đầy đủ và chi tiết. Báo cáo này trình bày một số kết quả trong nghiên cứu và áp dụng các mô hình thủy lực trên để mô phỏng và dự báo lũ trên hệ thống sông Hồng. Báo cáo gồm 3 phần. Phần 1 giới thiệu tổng quan về các mô hình thủy lực, phương pháp giải số và các mô đun tính toán tương ứng. Phần 2 trình bày một số kết quả áp dụng các mô hình cho một số bài toán mẫu. Phần 3 giới thiệu việc thu thập và xử lý số liệu trên hệ thống sông Hồng phục vụ cho các mô hình và kết quả áp dụng để mô phỏng lũ năm 2008 và dự báo lũ năm 2009. Mở đầu. Sử dụng hợp lý tài nguyên nước nói chung, kiểm soát lũ lụt nói riêng, trên một lưu vực sông phức tạp là một vấn đề lớn phục vụ trực tiếp cho việc phát triển kinh tế - xã hội và bảo vệ môi trường. Mô hình hóa là một công cụ quan trọng trong sử dụng tài nguyên nước. Trong hệ thống mô hình cho một lưu vực sông, luôn luôn có mặt 2 loại mô hình là mô hình thủy văn và mô hình thủy lực. Trong trường hợp dòng chảy trên lưu vực sông chủ yếu được hạn chế trong lòng sông do có đê bao (thí dụ như lưu vực hệ thống sông Hồng - sông Thái Bình), thì cùng với mô hình thủy văn, mô hình thủy lực 1 chiều là các mô hình quan trọng nhất. Vì vậy, trong báo cáo này chúng tôi xin trình bày một số vấn đề liên quan đến các mô hình thủy lực 1 chiều trong mô phỏng và dự báo lũ cho hệ thống sông Hồng - sông Thái Bình. 1. Giới thiệu tổng quan về các mô hình thủy lực 1 chiều Mô hình thủy lực 1 chiều, mô tả chuyển động của nước trông sông hoặc kênh hở, được xây dựng dựa trên hệ phương trình Saint Venant ([1]) đầy đủ (mô hình sóng động lực học SDLH - dynamic wave model). Một trong những dạng của hệ phương trình này là: tc + x = q (1) t + x  A + gA x + Sf  = 0 (2) Trong hệ phương trình (1), (2): Q=Q(x,t) – lưu lượng của dòng chảy; Z=Z(x,t) – mực nước; q -lưu lượng phụ; Ac - Diện tích mặt cắt (kể cả vùng chứa); A - Diện tích chảy; Sf - Sức cản đáy n2QQ f A2R4/3 (3) R là bán kính thủy lực; n - Hệ số Manning. Phương trình (1) là phương trình bảo toàn khối lượng, phương trình (2) là phương trình bảo toàn động lượng. Nếu bỏ đi 2 số hạng đầu trong phương trình bảo toàn động lượng (2), ta thu được phương trình: x + Sf = 0 (4) Mô hình được xây dựng từ các phương trình (1) và (4) là mô hình sóng khuếch tán SKT (diffusive wave model, [2], [3]). Mô hình này có thể áp dụng để mô phỏng sóng chậm như sóng lũ trong tự nhiên ([4],[5]). Nếu bỏ đi tiếp thành phần h trong biểu thức Z = h - S0 , h là độ sâu, S0 là độ dốc của đáy sông, ta thu được phương trình: S0 = Sf (5) Mô hình được xây dựng từ các phương trình (1) và (5) là mô hình sóng động học SDH (kinematic wave model, [2], [3]). Mô hình này có thể áp dụng để mô phỏng sóng lũ chậm trên sông có độ dốc lớn ([5]). Khác với mô hình SDLH và SKT, trong mô hình SDH chi cần cho điều kiện biên lưu lượng ở thượng du. Mô hình SDH không mô tả được ảnh hưởng của mực nước ở hạ du đến dòng chảy. Trong áp dụng mô hình thủy lực 1 chiều, cùng với 3 mô hình SDLH, SKT, SDH còn 2 cách tiếp cận (approach) sau được nhiều người sử dụng, đó là cách tiếp cận kiểu Muskingum và Muskingum-Cunge. Cách tiếp cận kiểu Muskingum (MUS) là cách đơn giản nhất trong áp dụng mô hình thủy lực 1 chiều. Trong cách tiếp cận này, người ta chỉ sử dụng phương trình bảo toàn khối lương (1) và để giải một phương trình gồm 2 ẩn (A, Q), MUS đã sử dụng một giả thiết quan trọng là lưu lượng Q có quan hệ với mặt cắt A: Q = f(A) (6) Đặt c = dQ/dA (c có tên gọi là vận tốc truyền sóng) và ta có thể viết phương trình (1) trong trường hợp không có lưu lượng phụ về dạng sau: dA Q Q 1 Q Q dQ t x c t x Từ phương trình này ta thu được phương trình Q Q t x Trong MUS, lưu lượng Q trong phương trình (7) được tính tương đối đơn giản qua việc lựa chọn các hệ số cho từng đoạn sông. Vì vậy, chúng tôi không áp dụng cách tiệm cận này trong mô phỏng lũ trên hệ thống sông Hồng - sông Thái Bình. Cách tiếp cận kiểu Muskingum-Cunge (MUSC) được phát triển từ MUS qua phương trình khuếch tán sau ([6]) Q Q 1 2Q t x 2BS0 x2 (8) trong đó B là chiều rộng của sông. Phương trình (8) có thể thu nhận được từ các phương trình (1), (4) và (5). Tương tự như mô hình SDH, để tính lưu lượng Q theo MUSC, cũng chỉ cần cho điều kiên biên lưu lượng ở thượng du và MUSC cũng không mô tả được ảnh hưởng của mực nước ở hạ du đến dòng chảy. Phương pháp sai phân được sử dụng để tìm nghiệm số trị trong các mô hình thủy lực. 2. Kết quả kiểm định các mô hình thủy lực cho một số bài toán mẫu 2.1. Mô tả các bài toán mẫu. Có nhiều bài toán mẫu để kiểm tra khả năng mô phỏng của mô hình thủy lực 1 chiều. Trong trường hợp mô phỏng sóng lũ trong tự nhiên chúng tôi đã chọn các bài toán mẫu số 2 (dòng chảy dừng), số 6 (sóng khuếch tán) và số 7 (sóng động học) trong [5]. Bài toán mẫu số 2 có nghiệm giải tích (xem [5]) và được chọn để kiểm tra độ chính xác của nghiệm số trị thu được từ các loại mô hình. Bài toán mẫu số 6 là bài toán mẫu đặc trưng cho SKT. Bài toán mẫu số 7 là bài toán mẫu đặc trưng cho SDH. Bài toán mẫu số 2 - Dòng chảy dừng. Bài toán mẫu mô tả dòng chảy dừng trong kênh chữ nhật dài 10,000 m, có độ dốc dáy, biên lưu lượng trên thượng du, biên mực nước dưới hạ du không thay đổi theo thời gian. Bài toán mẫu số 6 - Sóng khuếch tán. Bài toán mô tả sóng lũ trên một kênh chữ nhật dài 50,000 m, có độ dốc là 0.0005 với lưu lượng vào đầu kênh biến đổi từ Q0 = 750 m3/giây lên Qm =2,000 m /giây trong khoảng thời gian là T = 1,800 giây, sau đó lưu lượng lại giảm trở lại Q0 cũng trong khoảng thời gian là T. Bài toán mẫu số 7 - Sóng động học. Tương tự như sóng khuếch tán, sóng động học mô tả sóng lũ trên một kênh chữ nhật dài 100,000 m, có độ dốc lớn là 0.005 với lưu lương vào đầu kênh biến đổi từ Q0 =1,000 m3/giây lên Qm = 2,000 m3/giây trong khoảng thời gian tương đối dài T = 43,200 giây (12 giờ), sau đó lưu lượng lại giảm trở lại Q0 cũng trong khoảng thời gian là 12 giờ. Tài liệu [5] chỉ mô tả nghiệm chính xác của một số bài toán mẫu (thí dụ số 1, số 2, số 9). Tài liệu [5] không mô tả nghiệm chính xác của các bài toán mẫu số 6 và số 7. Trong báo cáo này, chúng tôi xin mô tả nghiệm chính xác của bài toán mẫu số 7. Từ công thức (3) tính Sf, phương trình (5) và giả thiết rằng chu vi ước có thể lấy bằng chiều rộng B của sông, ta thu được các quan hệ sau: Q = AR2/3 1 A5/3 S0 0 n nB2/3 = k1A5/3,k1 = nB2/3 3/5 A =    1  (9) dQ = c = 5k1A2/3 = k2Q2/5 = f (Q),k2 = 5k13/5 . (10) Từ phương trình (7) và đẳng thức (10) ta thu được phương trình tải phi tuyến sau: Q + f (Q)Q = 0 (11) Nghiệm Q của phương trình (11) cần thỏa mãn điều kiện ban đầu và điều kiện biên của bài toán mẫu số 7. Sử dụng cách biều diễn nghiệm cho phương trình đạo hàm riêng phi tuyến trong [7], ta có thể biểu diễn nghiệm Q như sau: Q = Q(x,t) = Q0 = Q0 + (Qm - Q0)(-x1+t)/T = Qm - (Qm - Q0)[(-x1+t)/T-1] = Q0 , nếu t ≤ x1 , x1 = x1(x,Q) = , nếu x1≤ t ≤ x1+T , nếu x1+T ≤ t ≤ x1+2T , nếu x1+2T ≤ t x f (Q) (12.1) (12.2) (12.3) (12.4) Áp dụng công thức tính đạo hàm của hàm ẩn ta có thể kiểm tra được ngay Q là nghiệm chính xác của bài toán mẫu số 7. Xét điểm giữa x2 = 50,000 của kênh. Tại điểm x2 này nghiệm Q bắt đầu tăng khi t = x1/f(Q0) = x1(x2,Q0)/f(Q0) ≈ 2.072 giờ. Kết quả này phù hợp với kết quả đã công bố trong [5] (hình 1). Từ thời điểm này đến x1+2T, x1 trong (12.2) và (12.3) phụ thuộc phi tuyến vào nghiệm Q. Để tìm nghiệm Q cần giải phương trình phi tuyến (12.2) và (12.3). Có thể xấp xỉ nghiệm Q này bằng biểu thức tuyến tính sau Hình 1. NghiệmQ ([5]) tại x=0 và x=50,000 m Q’ = Q0 + (Qm - Q0)(- f (Q0 ) +t)/T (13) Sai số của phép xấp xỉ này là : Q’ - Q = (Qm - Q0)( f (Q) - f (Q0 ) )/T = (Qm −Q0 )x * f (Q0 )− f (Q) f (Q) f (Q0 ) 6 5 4 3 2 Với các giá trị đã cho trong bài toán mẫu số 7, ta có thể chứng minh được │Q’ - Q│≤ 0.07 . Tương tự ta có thể xấp xỉ nghiệm Q trong (12.3) bằng biểu thức tuyến tính sau H_TT H_CX Q’ = Qm - (Qm - Q0)[(- f (Q0 ) +t)/T – 1] (14) 1 Nghiệm xấp xỉ (13) và (14) là các tuyến tính theo t 0 1 8 15 22 29 36 43 50 57 64 71 78 85 92 99 và phù hợp với kết quả đã công bố trong [5] (hình Hình 2. Kết quả kiểm định mô hình SDLH 2.2. Kết quả kiểm định các bài toán mẫu. Trong bài toán mẫu số 2 SDLH cho kết quả rất chính xác (hình 2) với sai số tuyệt đối nhỏ hơn 0.015 m. SKT cho kết quả tương đối chính xác với sai số tuyệt đối nhỏ hơn 0.15 m. MUSC và SDH trong trường hợp này cho kết quả sai (vì không mô tả được ảnh hưởng của mực nước ở hạ du đến dòng chảy) Trong bài toán mẫu số 6 2100 SDLH và SKT cho các kết quả trùng với kết quả đã được công bố 1900 trong [5]. MUSC và SDH cho kết 1700 quả xấu hơn 2 mô hình trên. Trong bài toán mẫu số 7. 1500 x=0 Mô hình SDH (hình 3) và MUSC 1300 x=50,000 cho kết quả trùng với kết quả đã được công bố trong [5] trong hình 1100 1. Mô hình SDLH và mô hình SKT 900 cho các kết quả giống nhau và xấu 0 21600 43200 64800 86400 108000 129600 151200 hơn kết quả của mô hình SDH và Thoi gian (giay) MUSC. Hình 3. Kết quả kiểm định mô hình SDH 3. Kết quả áp dụng các mô hình thủy lực để mô phỏng và dự báo lũ trên hệ thống sông Hồng Theo SDLH phần mềm tính toán thủy 1000 lực 1 chiều IMECH_1D cho mạng sông 900 Hồng - sông Thái Bình đã được chúng tôi 700 xây dựng và áp dụng để dự báo lũ từ năm 600 2002 đến nay. Kết quả dự báo được đánh giá 500 Thuc do là tương đối tốt. Hình 4 là kết quả so sánh 400 48 gio mực nước thực đo với mực nước dự báo 200 trước 24 giờ và trước 48 giờ tại trạm thủy 100 văn Hà Nội. 0 Trên thượng du của hệ thống sông Thoi gian (gio) Hồng, việc đo đạc mặt cắt các sông hiện còn gặp một số khó khăn. Vì vậy, chúng tôi chỉ thu thập được thông tin về 2 đợt đo đạc trên thượng du sông Lô (hình 5) và các mặt cắt được trích từ bản đồ trên thượng du sông Đà (hình 6). Hình 5. Cao độ đáy sông Lô qua 2 lần đo Hình 6. Mặt cắt trên thượng du sông Đà Với các số liệu này, chúng tôi đã áp dụng SDH để mô phỏng lũ trên thượng du sông Lô năm 2008 (hình 7). Chúng tôi cũng đã kết nối mô hình thủy văn MARINE với MUSC và áp dụng áp dụng thử nghiệm để dự báo lưu lượng vào hồ Hòa Bình 2009 (hình 8). H. 7. Mô phỏng lũ trên thượng du sông Lô. H. 8. Dự báo lưu lượng vào hồ Hòa Bình, 5/7/09 ... - tailieumienphi.vn