›
›

Báo cáo Nghiên cứu, ứng dụng mô hình kết nối Marine và IMech1D dự báo lưu lượng vào hồ Hòa Bình

pdf 90 trang phuongnguyen 1860
Download
Bạn đang xem 20 trang mẫu của tài liệu "Báo cáo Nghiên cứu, ứng dụng mô hình kết nối Marine và IMech1D dự báo lưu lượng vào hồ Hòa Bình", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

  • pdfnghien_cuu_ung_dung_mo_hinh_ket_noi_marine_va_imech1d_du_bao.pdf

Nội dung text: Báo cáo Nghiên cứu, ứng dụng mô hình kết nối Marine và IMech1D dự báo lưu lượng vào hồ Hòa Bình

  1. 1 TRƯỜNG . KHOA .  Báo cáo tốt nghiệp Đề tài: NGHIÊN CỨU, ỨNG DỤNG MÔ HÌNH KẾT NOOIS MARINE VÀ I,ECH1D DỰ BÁO LƯU LƯỢNG VÀO HỒ HÒA BÌNH
  2. 2 MỤC LỤC MỤC LỤC 1 MỞ ĐẦU 6 Chương 1. ĐẶC ĐIỂM CỦA LƯU VỰC VÀ CỦA CÁC MÔ HÌNH LỰA CHỌN NGHIÊN CỨU 8 1.1. Đặc điểm địa lý tự nhiên, quy luật dòng chảy lũ trên lưu vực sông Đà và vai trò của hồ Hòa Bình [4]-[6]: 8 1.1.1. Đặc điểm mưa gây lũ [4]: 8 1.1.2 Đặc điểm dòng chảy lũ sông Đà [4]: 9 1.1.3. Vai trò của hồ Hòa Bình [4]: 10 1.2. Tổng quan về mô hình thủy văn [3]-[5]: 11 1.3. Tổng quan về mô hình thủy lực [3]; [6]: 12 Chương 2. PHẦN MỀM THỦY VĂN THAM SỐ PHÂN BỐ MARINE 15 2.1. Cơ sở khoa học của phần mềm thủy văn tham số phân bố Marine: 15 2.1.1. Mô hình dòng chảy trên bề mặt lưu vực [3]; [5]; [6]; [8]: 15 2.1.2. Mô hình thấm Green Ampt [12]; [16]: 16 2.2. Cấu trúc dữ liệu trong Marine [16]: 18 Chương 3. PHẦN MỀM THỦY LỰC MỘT CHIỀU IMECH1D 21 3.1. Các thành phần của hệ thống [5]; [6]; [9]: 21 3.1.1. Mạng sông: 21 3.1.1.1. Nút sông: 21 3.1.1.2. Đoạn sông: 22 3.1.2. Ô ruộng (Ô chứa): 22 3.2. Mô hình toán học [5]; [6]; [9]: 22 3.2.1. Mô hình toán học một đoạn sông: 22 3.2.2. Mô hình toán học của một ô ruộng [5]; [6]; [9]: 23 3.3. Lược đồ sai phân [1]; [5]; [6]: 23 3.4. Tuyến tính hóa hệ phương trình (3.5), (3.7), (3.8): 25 3.4.1. Tuyến tính hoá các biểu thức đơn giản [10]: 26 3.4.2. Tuyến tính hoá biểu thức có lực cản đáy: 26
  3. 3 3.4.3. Tuyến tính hoá biểu thức trao đổi nước qua đê: 27 3.5. Thuật giải hệ phương trình đại số tuyến tính [1]; [5]; [6]: 27 3.6. Các thuật toán phụ trợ sử dụng trong xây dựng bộ chương trình tính toán thủy lực một chiều IMech1D [5]; [6]: 32 3.6.1. Khái toán mặt cắt: 32 3.6.2. Tạo giá trị mực nước và lưu lượng làm điều kiện ban đầu: . 33 3.6.3. Vấn đề xác định hệ số nhám và chỉnh kết quả: 34 Chương 4. ỨNG DỤNG MÔ HÌNH KẾT NỐI MARINE-IMECH1D CHO LƯU VỰC SÔNG ĐÀ PHẦN TRÊN LÃNH THỔ VIỆT NAM 35 4.1. Mô hình kết nối Marine và Imech1D [5]: 35 4.2. Xử lý số liệu cho mô hình kết nối Marine-IMech1D: 36 4.2.1. Xử lý bản đồ địa hình: 36 4.2.1.1. Xác định hướng của dòng chảy và độ tích tụ của dòng chảy trên DEM: 37 4.2.1.2. Tạo mạng sông suối từ DEM: 39 4.2.1.3. Phân chia lưu vực trên nền DEM: 40 4.2.2. Xử lý bản đồ phân loại đất: 40 4.2.3. Xử lý bản đồ hiện trạng sử dụng đất: 42 4.2.4. Xây dựng bản đồ phân bố mưa trong lưu vực: 42 4.2.5. Tích hợp các mặt cắt sông vào lớp sông suối trên nền DEM: 44 4.2.6. Xử lý các số liệu khác: 45 Chương 5. NÂNG CAO CHẤT LƯỢNG TÍNH TOÁN CỦA MÔ HÌNH BẰNG KỸ THUẬT LỌC KALMAN 46 5.1. Quá trình cần đánh giá (ước lượng) [13]: 46 5.2. Các vấn đề tính toán (bản chất tính toán) của lọc Kalman [13]; [17]: 48 5.2.1. Định nghĩa về các ước lượng tiên nghiệm và hậu nghiệm: 48 5.2.2. Bước dự báo – cập nhật (ước lượng tiên nghiệm): 48 5.2.3. Bước hiệu chỉnh (ước lượng hậu nghiệm): 50 5.2.4. Tìm Kalman gain (blending factor) K: 51 5.3. Thuật toán lọc Kalman rời rạc: 54
  4. 4 5.3.1.Cập nhật theo thời gian – dự báo (ước lượng tiên nghiệm) (predict): 54 5.3.1.1. Phép tính s 1: 54 5.3.1.2. Phép tính số 2: 54 5.3.2. Cập nhật theo đo đạc – chỉnh sửa (ước lượng hậu nghiệm) (correct): 55 5.3.2.1. Phép tính số 3: 55 5.3.2.2. Phép tính số 4: 55 5.3.2.3. Phép tính số 5: 55 Chương 6. CÁC KẾT QUẢ ỨNG DỤNG MÔ HÌNH KẾT NỐI MARINE VÀ IMECH1D ĐỂ DỰ BÁO LƯU LƯỢNG VÀO HỒ HÒA BÌNH 56 6.1. Kết quả kiểm tra bài toán mẫu cho 10 lưu vực bộ phận: 56 6.2. Sử dụng mô hình kết nối Marine-IMech1D để dự báo lại trận lũ năm 2006 và hiệu chỉnh các tham số của mô hình: 57 6.2.1. Nhận định chung tình hình lũ sông Đà năm 2006: 57 6.2.2. Kết quả tính toán dự báo lại cho trận lũ năm 2006 bằng mô hình kết nối MARINE-IMECH1D: 57 6.3. Kết quả sử dụng mô hình kết nối Marine-IMech1D tác nghiệp cho mùa lũ năm 2009: 60 KẾT LUẬN 63 DANH MỤC CÔNG TRÌNH CỦA TÁC GIẢ 64 TÀI LIỆU THAM KHẢO 66 Tiếng Việt 66 Tiếng Anh 66 Tiếng Pháp 67 PHỤ LỤC 68 1. Kết quả kiểm tra bài toán mẫu cho 10 lưu vực bộ phận của lưu vực sông Đà bằng MARINE: 68 2. Kết quả kiểm định bộ chương trình tính toán thủy lực một chiều IMech1D bằng các bài toán kiểm định mẫu (Test Cases). 74 2.1. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 1: SÓNG XẢ TRONG KÊNH CHỮ NHẬT 74
  5. 5 2.2. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 2: DÒNG CHẢY ÊM, ĐỀU TRONG KÊNH HÌNH CHỮ NHẬT 75 2.3. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 3: DÒNG CHẢY ĐỀU CÓ LƯU LƯỢNG PHỤ 76 2.4. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 4: DÒNG CHẢY ĐỀU CÓ CÔNG TRÌNH 77 2.5. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 5: SÓNG ĐỘNG LỰC HỌC 78 2.6. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 6: SÓNG KHUẾCH TÁN 79 2.7. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 7: SÓNG ĐỘNG HỌC 80 2.8. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 8: SÓNG LŨ QUA HỒ CHỨA 81 2.9. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 9: NHIỄU ĐỊA PHƯƠNG TRONG DÒNG CHẢY DỪNG 82 2.10. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 10: HÌNH HỌC KHÔNG ĐỀU TRONG DÒNG CHẢY DỪNG 83 2.11. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 11: DÒNG CHẢY KHÔNG DỪNG TRONG KÊNH CÓ LÒNG DẪN PHỨC HỢP 85 2.12. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu 12: PHÂN LƯU 88
  6. 6 MỞ ĐẦU Do ảnh hưởng của tình trạng biến đổi khí hậu toàn cầu và hiện trạng khai thác, sử dụng đất trên bề mặt lưu vực bị thay đổi nhiều nên lũ lụt trên hệ thống sông Đà đang có chiều hướng ngày một diễn biến phức tạp hơn. Chính điều này đã ngày một gây thêm nhiều khó khăn cho công tác dự báo lưu lượng vào hồ Hòa Bình cũng như công tác phòng chống lụt bão và điều hành hồ chứa thủy điện Hòa Bình. Sông Đà là một hệ thống sông lớn và là chi lưu lớn nhất trong ba chi lưu của hệ thống sông Hồng. Lưu lượng nước đổ về sông Hồng phần lớn là từ sông Đà chảy về. Chính vì vậy sông Đà có ảnh hưởng rất lớn đến tình trạng lũ lụt trên khu vực đồng bằng châu thổ sông Hồng. Để phòng chống lũ cho khu vực đồng bằng châu thổ sông Hồng và đặc biệt là chống lũ cho thành phố Hà Nội đòi hỏi phải kiểm soát được lũ sông Đà và vận hành công trình chống lũ hồ thủy điện Hòa Bình một cách hợp lý. Dự báo trước lưu lượng vào hồ Hòa Bình là một nhiệm vụ quan trọng. Bởi vì, để vận hành dược hồ Hòa Bình phục vụ đa mục tiêu cần phải biết trước được lưu lượng vào hồ. Luận văn thạc sỹ này được đặt ra trong hoàn cảnh thực tế là cần phải xây dựng một công cụ cho phép dự báo trước lưu lượng vào hồ Hòa Bình từ các số liệu đầu vào đã biết bao gồm: - Các thông tin về lưu vực sông Đà bị giới hạn nằm trên phần lãnh thổ Việt Nam. - Số liệu mưa thực đo đã biết của các trạm đo trong lưu vực. - Số liệu mưa dự báo tại các trạm đo trong lưu vực bởi các mô mình dự báo mưa. - Lưu lượng chảy vào từ phần lưu vực thuộc lãnh thổ Trung Quốc. - Số liệu về lưu lượng và mực nước tại một số trạm đo trên hệ thống sông Đà và của hồ Hòa Bình. - Các thông số của hồ Hòa Bình, công trình thủy điện Hòa Bình - Thông tin về địa hình, hiện trạng sử dụng đất, của lưu vực sông Đà. - Các thông tin phụ trợ khác. Trên cơ sở đó nội dung của luận văn được đặt ra với mục tiêu là khai thác, sử dụng các mô hình toán tiên tiến mà thế giới hiện đang nghiên cứu phát triển
  7. 7 để xây dựng công cụ dự báo trước lưu lượng vào hồ Hòa Bình 48 giờ. Qua sự nghiên cứu, phân tích nhiều mô hình thủy văn và thủy lực khác nhau, cuối cùng mô hình thủy văn Marine và mô hình thủy lực IMech1D được lựa chọn để phát triển và kết nối thành mô hình kết nối Marine và IMech1D phục vụ bài toán dự báo lưu lượng vào hồ Hòa Bình. Trên cơ sở đó nội dung của luận văn bao gồm 6 chương chính là: Chương 1: trình bày các thông tin tổng quan về đề tài bao gồm, thông tin về lưu vực nghiên cứu, các thông tin về mô hình thủy văn, thủy lực được lựa chọn nghiên cứu. Chương 2: Nghiên cứu cơ sở khoa học và phát triển mô hình thủy văn tham số phân bố Marine. Chương 3: Nghiên cứu cơ sở khoa học và khai thác mô hình thủy lực một chiều IMech1D. Chương 4: Kết nối mô hình thủy văn Marine với mô hình thủy lực một chiều IMech1D thành một mô hình thống nhất. Ứng dụng mô hình này cho lưu vực sông Đà để dự báo lưu lượng vào hồ Hòa Bình Chương 5: Nghiên cứu ứng dụng kỹ thuật lọc Kalman để nâng cao độ chính xác của mô hình đã kết nối Marine và IMech1D. Chương 6: Trình bày các kết quả sử dụng mô hình kết nối Marine- IMech1D dự báo lưu lượng vào hồ Hòa Bình trước 48 giờ. Cuối cùng là phần kết luận và một số phụ lục.
  8. 8 Chương 1. ĐẶC ĐIỂM CỦA LƯU VỰC VÀ CỦA CÁC MÔ HÌNH LỰA CHỌN NGHIÊN CỨU 1.1. Đặc điểm địa lý tự nhiên, quy luật dòng chảy lũ trên lưu vực sông Đà và vai trò của hồ Hòa Bình [4]-[6]: Sông Đà có diện tích 52.900 km2, lưu vực hẹp, kéo dài theo hướng tây bắc- đông nam tới 380km, rộng trung bình 80km, phần thuộc địa phận Việt Nam có diện tích là 26800km2, chiếm khoảng 50,7% diện tích toàn lưu vực. Lòng chính ở thượng lưu hẹp, nhiều thác ghềnh; ở hạ lưu, lòng sông mở rộng, độ dốc sông trung bình là 3,58%. Độ cao bình quân lưu vực là 965m, độ dốc bình quân là 36,8%. Đường phân thủy phía đông của lưu vực là dãy núi Hoàng Liên Sơn, Pu Luông với đỉnh cao từ 2500-3000m. Phía tây có dãy núi cao Phu Huổi Long, Phutama, Phu Tung và Phu Sang. Phía bắc có dãy núi cao Pusi-Lung và Ngũ Đài Sơn, phía đông nam là vùng núi thấp Ba Vì, Viên Nam và Đối Thôi. Địa hình lưu vực là dạng núi và cao nguyên đều cao, chia cắt mạnh theo chiều thẳng đứng; các dãy núi, cao nguyên và thung lũng xếp song song theo hướng tây bắc- đông nam. 1.1.1. Đặc điểm mưa gây lũ [4]: Sự sắp xếp song song của địa hình núi, cao nguyên và thung lũng sông có tác động rõ rệt tới khí hậu trên lưu vực. Dãy núi cao Hoàng Liên Sơn - Puluông như một bức tường tự nhiên ngăn cản và làm suy yếu ảnh hưởng của gió đông bắc. Các dãy núi cao ở biên giới Việt-Lào lại tạo ra hiệu ứng fơn đối với gió mùa tây nam. Điều kiện địa hình và vị trí của lưu vực đã qui định khí hậu với hai mùa: mùa đông khô lạnh, mùa hè nhiều mưa ở vùng cao và khô nóng ở vùng thấp. Mưa lớn trên lưu vực thường bắt đầu sớm, vào khoảng tháng VI, tháng VII. Vùng bắc và tây bắc là vùng núi cao có khí hậu ẩm ướt đến rất ẩm, lượng mưa trung bình nhiều năm từ 1500 đến 2700mm, lượng mưa mùa hè (tháng V- IX) chiếm tới trên 70% tổng lượng mưa năm. Vùng núi thấp Sơn La-Mộc Châu, mùa hè chịu ảnh hưởng của gió mùa tây nam, lượng mưa trung bình năm thấp, chỉ 1100 đến 1500mm, trong đó lượng mưa mùa hè dưới 1000mm. Trên lưu vực sông Đà tồn tại những trung tâm mưa lớn như trung tâm mưa ở sườn tây dãy Hoàng Liên Sơn thuộc các lưu vực sông nhánh Nậm Na, Nậm Mu, lượng mưa trung bình năm khoảng 2500mm (trên lưu vực Nậm Na-
  9. 9 mưa trung bình năm tới trên 2000mm : tại Phong Thổ lượng mưa trung bình năm là 2202mm, PaTần 2997mm, Sình Hồ 2682mm; trên lưu vực Nậm Mu lượng mưa trung bình năm tới 2454mm, ở thượng lưu lên tới 2700-2800mm). Tại vùng phía tây dãy Hoàng Liên Sơn thấy rõ qui luật lượng mưa tăng theo độ cao lưu vực, mưa tập trung vào các tháng V-X, đặc biệt là các tháng VI-VIII; lượng mưa mùa hè chiếm trên 90%, lượng mưa các tháng VI-VIII chiếm 50- 60% lượng mưa năm. Trung tâm mưa lớn tại phần lưu vực thuộc địa phận Việt Nam gần biên giới Việt-Trung là tâm mưa lớn nhất, lượng mưa năm thay đổi tùy từng vị trí từ 2400 đến 3000mm, mưa tập trung nhiều nhất vào các tháng VI-VIII. 1.1.2 Đặc điểm dòng chảy lũ sông Đà [4]: Địa hình núi cao, chia cắt mạnh, độ dốc lớn, thung lũng sâu, hẹp với lượng mưa lớn lại tập trung vào một vài tháng trong năm nên tạo điều kiện hình thành mạng lưới sông dày đặc, ít sông lớn, hướng của các dòng sông suối trùng với hướng của lưu vực. Mật độ sông suối lớn nhất ở vùng núi phía tây Hoàng Liên Sơn lên tới 1,5-1,7km/km2. Phía hữu ngạn sông Đà, do có lượng mưa ít hơn đáng kể so với các vùng khác nên sông suối thưa hơn, chỉ từ 0,5 đến 1,5km/km2, thường dưới 1,0km/km2. Trên sông Đà và các dòng sông nhánh như Nậm Na, Nậm Mu, các sông suối nhỏ đổ vào dòng chính thường phân bố đều dọc sông. Vùng cao nguyên đá vôi mưa ít, sông suối thưa, dòng chảy nhỏ hơn. Nguồn sinh dòng chảy quan trọng nhất trên sông Đà nằm ở phần lưu vực thuộc vùng biên giới Việt - Trung và vùng sườn phía tây dãy Hoàng Liên Sơn, nơi có môdun dòng chảy năm từ 30-40 l/s/km2 và hơn nữa. Ở các nơi khác trên lưu vực, lượng dòng chảy thường không vượt quá 20 l/s/km2 (biểu 1.3). Dòng chảy sông tập trung vào các tháng mùa lũ, chiếm tới 69-78% tổng lượng dòng chảy năm. Mùa lũ trên sông Đà thường bắt đầu vào tháng V, kết thúc vào cuối tháng IX đầu tháng X. Lũ lớn nhất thường xảy ra vào cuối tháng VII, nửa đầu tháng VIII. Dòng chảy lũ trên sông Đà lớn, tập trung nhanh và không đồng bộ ở các phần khác nhau của lưu vực là một đặc điểm nổi bật nhất của dòng chảy sông Đà. Trong điều kiện địa lý tự nhiên thuận lợi cho dòng chảy lũ hình thành trên các phụ lưu sông Đà, nhất là lưu vực Nậm Na, Nậm Mu hai phụ lưu lớn nhất bên tả ngạn, thường xuất hiện những trận lũ đặc biệt lớn gây tác hại nghiêm
  10. 10 trọng. Mô đun dòng chảy lũ lớn nhất đạt tới 2000-3000 l/s/km2 - thuộc loại lớn nhất ở Việt Nam. Trên dòng chính, lượng dòng chảy lũ chiếm bình quân từ 77,6 đến 78,5% dòng chảy năm, dòng chảy tháng VIII- tháng có dòng chảy lớn nhất năm - chiếm tới 23,7% dòng chảy năm. Dòng chảy lũ sông Đà thuộc loại lớn nhất trên hệ thống sông Hồng. Mô đun đỉnh lũ tại Lai Châu là 324 l/s/km2 xảy ra vào các tháng VII năm 1966 và 428 l/s/km2 vào tháng VIII năm 1945. Mô đun đỉnh lũ tại Hòa Bình lên tới 454 l/s/km2 vào tháng VII năm 1964. Nhìn chung, trên đoạn sông từ Lai Châu về Hòa Bình thấy rõ quy luật tăng dần môdun dòng chảy cực đại khi diện tích lưu vực tăng. Điều này chứng tỏ rằng lượng gia nhập đáng kể ở phần lưu vực thuộc địa phận Việt Nam. Tại Lai Châu, biên độ lũ lớn nhất đạt tới 25 mét, cao nhất ở Việt Nam, với cường suất lũ lên lớn nhất tới 77,4 cm/h. Dòng chảy lũ tập trung nhanh như vậy nên công tác dự báo thủy văn gặp khó khăn lớn, mà để giải quyết vấn đề này đòi hỏi phải có một mô hình tương đối nhạy với qúa trình thay đổi dòng chảy trong sông. Trên cơ sở xác định thời gian truyền lũ trung bình ở các đoạn sông chính từ Mường Tè về tới Hòa bình và trên các phụ lưu chính Nậm Na, Nậm Mu, thấy rằng, trong mùa lũ, thời gian truyền lưu lượng có ít nhiều khác nhau khi lũ lên và lũ xuống, tuy nhiên trong tính toán và dự báo có thể lấy thời gian trung bình truyền lũ từ Lai Châu về tới Tạ Bú là 12-18 giờ, từ Tạ Bú về tới Hòa Bình là 12-24 giờ trong tự nhiên, hiện nay khi có hồ chứa Hòa Bình thời gian truyền lũ rút ngắn còn 6 - 12h tuỳ theo mực nước hồ. Lưu ý rằng, thời gian truyền lũ trên các đoạn sông chính và trên các phụ lưu còn phụ thuộc vào vị trí tâm mưa trên lưu vực. Như vậy, với điều kiện kỹ thuật thủy văn và điều kiện thông tin khí tượng thủy văn hiện có thì thời gian dự kiến thực tế của dự báo lưu lượng và mực nước tại trạm Tạ Bú trên sông Đà không thể vượt quá 18 giờ, tại Hòa Bình - không vượt quá 36 giờ. Để kéo dài thời gian dự kiến của dự báo có thể sử dụng các thông tin về lượng mưa dự báo trong thời gian dự kiến ở các phần lưu vực khác nhau. Dự báo lượng mưa trong 24 giờ và 48 giờ tới, trong điều kiện hiện nay, thường kém chính xác cả về định tính và định lượng. Do vậy, những trị số dự báo lưu lượng và mực nước trên sông Đà với thời gian dự kiến trên 36 giờ chỉ nên dùng để tham khảo hoặc chỉ nên xem như những nhận định khả năng. 1.1.3. Vai trò của hồ Hòa Bình [4]: Hồ Hòa Bình là hồ thủy điện lớn được xây dựng với nhiều mục đích sử dụng khác nhau bao gồm: - Mục đích chống lũ
  11. 11 - Mục đích phát điện - Mục đích tưới tiêu - Phục vụ phát triển giao thông đường thủy - Chính vì hồ có nhiều vai trò trong các lĩnh vực khác nhau nên việc vận hành hồ sao cho đảm bảo được lợi ích tối đa của các mục tiêu là công việc hết sức khó khăn. Để hỗ trợ công việc điều hành hồ cần phải nhiều thông tin hỗ trợ. Một trong các thông tin hỗ trợ vô cùng quan trọng đó là dự báo trước lưu lượng vào hồ. Đây là một trong những lý do mà tác giả chọn đề tài cho luận văn thạc sỹ của mình. 1.2. Tổng quan về mô hình thủy văn [3]-[5]: Hiện nay trên thế giới có rất nhiều mô hình thủy văn khác nhau. Các mô hình này có thể được phân thành hai loại chính đó là: - Mô hình thủy văn dạng hộp đen: Đây là các mô hình được xây dựng và phát triển từ rất lâu rồi. Các mô hình này chủ yếu được xây dựng trên cơ sở các công thức kinh nghiệm. Chính vì vậy các mô hình này đòi hỏi rất ít dữ liệu (thông tin đầu vào), nhưng ngược lại, việc sử dụng các mô hình này đòi hỏi phải có nhiều kinh nghiệm mới đạt được kết quả tốt. Loại các mô hình này rất thích hợp với điều kiện thiếu các thông tin đầu vào của lưu vực. Ngày nay trên thế giới các mô hình loại này đã và đang được sử dụng rất hạn chế. Thế giới đang có su hướng chuyển sang sử dụng các mô hình hiện đại hơn đó là các mô hình thủy văn tham số phân bố (trình bày ở phần dưới đây). - Mô hình thủy văn tham số phân bố: Đây là loại mô hình thủy văn được xây dựng trên cơ sở giải hệ phương trình Saint Venant hai chiều. Loại mô hình này hiện đang được nghiên cứu phát triển và ứng dụng ở các nước tiên tiến. Đây chính là su hướng phát triển của mô hình thủy văn. Mô hình thủy văn loại này đòi hỏi nhiều thông tin đầu vào, các thông tin phải chi tiết, chính xác. Kết quả tính của các mô hình này cũng ít bị phụ thuộc vào kinh nghiệm của người sử dụng. Trong thời kỳ phát triển của khoa học kỹ thuật, đặc biệt là sự phát triển vượt bậc của công nghệ GIS thì mô hình thủy văn tham số phân bố là loại mô hình có ưu thế sử dụng nhiều hơn. Các mô hình loại này đã và đang được phát triển có thể kể đến các mô hình sau: - GBHM (của Nhật) - Casd2D-SET (của Mỹ)
  12. 12 - Marine (của Pháp) - TopModel (của Mỹ) - SWAT (của Mỹ) - HydroTel (của Canada) - Và nhiều mô hình khác Trong các mô hình đó, mỗi mô hình có một số dặc điểm riêng khác nhau tùy từng mục tiêu của bài toán mà mô hình đó đặt ra. Marine là mô hình được xây dựng để giải quyết bài toán lũ lớn trên các lưu vực có độ dốc cao. Đặc điểm nầy rất phù hợp với lưu vực sông Đà, vì lưu vực sông Đà có độ dốc lớn, lũ sông đà cũng lớn (trong mùa lũ) do lượng mưa lớn. Hơn thế Viện Cơ học có toàn bộ mã nguồn của Marine nên có thể tiếp tục phát triển, tùy biến sao cho phù hợp với điều kiện của Việt Nam. Đây chính là lý do Marine được chọn để kết nối với IMech1D trong luận văn này. 1.3. Tổng quan về mô hình thủy lực [3]; [6]: Với mục tiêu là xây dựng được một bộ chương trình tính toán thuỷ lực một chiều có khả năng mô phỏng được dòng chảy trong hệ thống sông phức tạp với các đặc điểm riêng biệt như hệ thống sông ở miền Bắc Việt Nam, cụ thể là hệ thống sông Hồng và sông Thái Bình. Trong các năm 70, 80, 90 các nhà khoa học Việt Nam đã có nhiều nỗ lực để thực hiện việc này. Bộ chương trình VRSAP (Vietnam River System And Plains) của Giáo sư, Anh hùng Lao động Nguyễn Như Khuê và cộng sự được đánh giá như là một thành tựu lớn trong việc xây dựng mô hình tính toán thuỷ lực một chiều cho hệ thống sông Hồng - Thái Bình. Tuy vậy, bộ chương trình này đã được xây dựng trên nền của trình độ toán học, tin học của nước ta ở những năm 80, còn phải chấp nhận nhiều giả thiết nên khả năng mô phỏng các đặc thù của dòng chảy trên hệ thống sông Hồng - Thái Bình của VRSAP còn hạn chế. Những năm qua, trong khuôn khổ của các dự án quốc tế, chúng ta có được một số chương trình tính toán thuỷ lực một chiều đã được thương mại hoá như WENDY (Hà Lan), MIKE (Đan Mạch). Vì các bộ chương trình này được chuyển giao cho chúng ta không có mã nguồn và đi kèm một số lượng khóa cứng hạn chế nên việc sử dụng nhất là việc kết nối với các bộ chương trình khác sẽ gặp nhiều vấn đề khó vượt qua. Gần đây nhiều chuyên gia thuỷ lực của ta đã khai thác bộ chương trình tính toán thuỷ lực của Hoa Kỳ FLDWAV. Có thể nói bộ chương trình FLDWAY đã
  13. 13 sử dụng 1 số thành tựu mới trong lĩnh vực toán học, thuỷ lực, tin học, nhưng bộ chương trình này chỉ mô phỏng cho hệ thống sông hình cây, trong khi đó hệ thống sông Hồng - Thái Bình lại có nhiều đoạn vòng. Các bộ chương trình khác lấy tự do từ internet cũng đã được sử dụng như HEC, HEC-RAS. Trong khuôn khổ của đề tài "Nghiên cứu cơ sở khoa học cho các giải pháp tổng thể dự báo phòng tránh lũ lụt ở đồng bằng sông Hồng" nhóm nghiên cứu về lũ lụt của Viện Cơ học đặt ra nhiệm vụ là kế thừa các thành tựu của các tác giả khác và của nhóm trong thời gian trước để xây dựng được một bộ chương trình tính toán thuỷ lực một chiều IMech1D có khả năng mô phỏng được các đặc thù của dòng chảy trên hệ thống sông Hồng - Thái Bình: a) Quá trình lưu lượng ở thượng lưu được điều tiết bằng các hồ chứa lớn như Thác Bà, Hoà Bình, Sơn La, Tuyên Quang. b) Khi cần thiết phải vận hành công trình phân lũ đập Đáy, vận hành các khu chậm lũ Tam Thanh, Lập Thạch, Lương Phú. c) Hệ thống sông Hồng - Thái Bình có nhiều đoạn vòng và nhiều khu bối. d) Ở các vùng gần cửa sông dòng chảy bị ảnh hưởng bởi chế độ triều Biển Đông. Bộ chương trình do nhóm đề tài xây dựng đã vượt qua được các bài toán kiểm định mẫu của các Phòng Thí nghiệm thuỷ lực lớn ở Châu Âu đề xuất và được hiệu chỉnh, kiểm định bằng số liệu của các trận lũ đã xảy ra cho các năm: 1996, 1999 và 2000. Trong mùa lũ từ năm 2002 đến nay, nhóm đề tài đã sử dụng bộ chương trình để dự báo quá trình lũ trên hệ thống sông Hồng - Thái Bình. Kết quả dự báo lũ của nhóm đề tài đã được Uỷ ban Phòng chống lụt bão Trung ương sử dụng để tham khảo khi ra các quyết định nhằm kiểm soát lũ trên hệ thống sông Hồng, Thái Bình. IMech1D đã được chính thức chuyển giao cho một số cơ quan để phục vụ công tac phòng chống lụt bão cho hệ thống sông Hồng sông Thái Bình và hệ thống sông Hương ở Thừa Thiên Huế. Cụ thể là đã chuyển giao cho: - Cục phòng chống lụt bão và Quản lý Đê điều của Bộ Nông nghiệp và Phát triển Nông thôn. - Trung tâm Dự báo Khí tượng Thủy văn trung ương của Bộ Tài nguyên và Môi trường. - Trung tâm dự báo Khí tượng Thủy văn tỉnh Thừa Thiên Huế
  14. 14 - Chi cục Phòng chống lụt bão và Quản lý Đê điều tỉnh Thừa Thiên Huế Từ năm 2004 đến nay (2009) IMech1D đã được chính thức sử dụng trong công tác dự báo mực nước các sông trong hệ thống sông Hồng và sông Thái Bình trong mùa lũ (tác nghiệp) tại Viện Cơ học và Trung tâm Dự báo Khí tượng Thủy văn Trung ương để hỗ trợ công tác điều hành các hồ chứa và phòng chống lụt bão.
  15. 15 Chương 2. PHẦN MỀM THỦY VĂN THAM SỐ PHÂN BỐ MARINE Marine là phần mềm thủy văn tham số phân bố, được nghiên cứu và phát triển bởi Viện Cơ học Chất lỏng Toulouse, Cộng hòa Pháp. Trong khuôn khổ dự án quốc tế FLOCODS, Marine được chuyển giao cho Viện Cơ học năm 2004. Từ đó tới nay Marine luôn được nghiên cứu cải tiến và phát triển cho phù hợp với điều kiện ứng dụng ở Việt Nam. Marine đã được sử dụng để nghiên cứu tính toán cho nhiều lưu vực sông ở Việt Nam. Trên thế giới Marine được đánh giá cao và được khai thác sử dụng ở nhiều nước như Pháp, Hà Lan, Brasin, 2.1. Cơ sở khoa học của phần mềm thủy văn tham số phân bố Marine: 2.1.1. Mô hình dòng chảy trên bề mặt lưu vực [3]; [5]; [6]; [8]: MARINE mô phỏng quá trình hình thành dòng chảy sinh ra bởi mưa trên lưu vực dựa trên phương trình bảo toàn khối lượng: V u.grad(V) P0 t (2.1) Trong đó: V: là thể tích khối chất lỏng xét. U: là vận tốc của dòng chảy giữa các ô lưới. P0: là lượng mưa. Vì: u.grad(V) div(V.u) V.div(u) Với chất lỏng không nén được ta có div ( u ) 0 , sử dụng công thức Green- Ostrogradski div(m.u).dS m.u.n.d S  V .dS V.u.n.d P t 0 Từ (2.1) suy ra : S  S (2.2) Vận tốc của dòng chảy trao đổi giữa các ô được tính theo công thức: H 2 / 3 u pente. (2.3) K m Vì lưới sử dụng để tính toán là lưới vuông (DEM) nên thay biểu thức vận tốc vào phương trình tích phân ta thu được:
  16. 16 4 5/ 3 H j t (2.4) H . pente. P . t  K x 0 j 1 m Trong đó: Pente: độ dốc; Km: hệ số nhám Manning x: chiều rộng ô lưới t : Bước thời gian tính j: Hướng chảy của ô lưới (j =1  4) H: Độ sâu mực nước của ô lưới tính. H: Sự thay đổi mực nước của ô lưới tính từ thời điểm t1 đến t2 Đây chính là phương trình tính sự biến thiên mực nước theo thời gian của mỗi ô lưới. Từ sự biến thiên mực nước H của mỗi ô lưới ta tính được tổng lưu lượng trao đổi của mỗi ô (bao gồm lưu lượng nhận từ mưa, lưu lượng chảy vào và lưu lượng chảy ra) tại mỗi bước tính chính bằng sự biến thiên thể tích nước chứa trong ô. Q= H*dx*dx Trong đó: dx là kích thước của lưới tính. Đối với lưu vực kín, lưu vực chỉ có một điểm thoát nước, tại điểm thoát nước của lưu vực ta luôn có lưu lượng ra khỏi lưu vực là: q= Q Đối với lưu vực hở, lưu vực nằm dọc hai bên bờ sông nên có nhiều điểm thoát nước. Với trường hợp này lưu lượng ra khỏi lưu vực là tổng lưu lượng trao đổi của các điểm thoát nước: q=∑Q=∑ H*dx*dx Như vậy kết quả quá trình tính toán của MARINE cho ta lưu lượng ra của các lưu vực. Đây chính là thành phần ra nhập dòng bên q cần trong mô hình thủy lực IMECH-1D. 2.1.2. Mô hình thấm Green Ampt [12]; [16]: Mô hình MARINE tính toán thấm dựa trên lý thuyết thấm Green Ampt từ phương trình liên tục và định luật Darcy. Độ sâu thấm tích lũy tiềm năng được tính bằng phương trình Green - Ampt: F(t) (2.5) F(t)   1 kt   Trong đó: F(t) là độ sâu luỹ tích của nước thấm vào trong đất
  17. 17 : Cột nước mao dẫn của mặt ướt  = -i với  là độ rỗng của đất , i là độ ẩm của đất k: Độ dẫn thuỷ lực Phương trình (2.5) là phương trình phi tuyến, ta có thể giải bằng phương pháp thay thế liên tiếp, hoặc phương pháp lặp Newton. Trong trường hợp độ sâu lớp nước đọng ho không thể bỏ qua ta phải thay thế  bằng giá trị  - ho trước khi giải. Sau khi tìm được độ sâu thấm tích lũy tiềm năng F(t) ta xác định được tốc độ thấm tiềm năng   f k( 1) F (2.6) Theo định luật Darcy ta xác định được lưu lượng thấm của mỗi ô lưới sẽ là: qthấm = f*dx*dx Trước khi sinh nước đọng (t tp), vùng bão hòa trên mặt đất lan dần xuống tầng đất sâu hơn và dòng chảy trên mặt đất bắt đầu xuất hiện từ lượng nước đọng. Thấm được chia 3 thời kỳ: Trước khi xuất hiện nước đọng, toàn bộ lượng mưa đều thấm xuống đất (t tp). Khi t = tp, lượng thấm tích lũy tại thời điểm sinh nước đọng tp được tính bởi công thức: Fp = i*tp nên f = i   Thay vào (2.6) có: i k ( 1 ) it p   dẫn tới t p k (2.7) i(i K ) Sau khi có nước đọng, lượng thấm tích lũy được tính theo công thức:   F F F  ln( ) k(t t ) (2.8) p   F p p
  18. 18 P0 (mưa) Cấ u trúc lưu vực (độ cao địa hình, Thấm - Green Ampt đất, thảm phủ) Dòng chảy mặt Liên kết các ô lưới, tính Lớp nước đọng Q trao đổ̉̉i giữa các ô tại ô lưới Quá trình Q~t tại cửa ra của lưu Dòng chảy sát mặt vực hoặc nút đăng ký tại ô lưới Hình 2.1: Sơ đồ mô tả mô hình MARINE Như vậy, tại mỗi bước thời gian mô đun thấm cho ta lưu lương thấm qthấm của mỗi ô chứa. Với trường hợp tính toán MARINE có sử dụng mô đun thấm tại mỗi bước thời gian tổng lưu lượng trao đổi của mỗi ô bao gồm lưu lượng nhận từ mưa, lưu lương chảy vào, lưu lượng chảy ra và lưu lượng thấm: Q= H*dx*dx - qthấm 2.2. Cấu trúc dữ liệu trong Marine [16]: Trong Marine dữ liệu được đưa vào dưới dạng các lớp dữ liệu chồng lên nhau. Mỗi lớp dữ liệu chứa một loại thông tin khác nhau. Các lớp dữ liệu được để ở dạng raster, sử dụng mã ASCII để lưu trữ. Các lớp thông tin vào của mô hình gồm có: - Lớp Thông tin về địa hình: Đây là lớp thông tin qua trọng nhất vì nó trực tiếp ảnh hưởng đến vận tốc dòng chảy, hướng của dòng chảy. - Lớp thông tin về hiện trạng sử dụng bề mặt lưu vực: Lớp thông tin này phản ánh độ nhám của bề mặt lưu vực, nó trực tiếp quyết định đến hệ số cản của dòng chảy trên bề mặt lưu vực. - Lớp thông tin về phân loại đất trên bề mặt lưu vực: Lớp thông tin này trực tiếp ảnh hưởng đến cường độ thấm của dòng chảy trên bề mặt lưu vực. Với mỗi loại đất sẽ có 3 lớp thông tin cần đưa vào mô hình (theo yêu cầu của phương pháp tính tổn thất Green Ampt). Như vậy có nghĩa là các thông tin về đất được chứa trong ba lớp dữ liệu. - Lớp thông tin về hệ thống sông suối và mặt cắt của sông suối.
  19. 19 - Lớp thông tin về mưa: Đây là lớp thông tin động, vì giá trị của lớp thông tin này luôn thay đổi theo thời gian. Chính vì vậy sẽ có nhiều lớp thông tin về mưa được đưa vào, thời gian thay đổi lớp thông tin mưa chính là ốp thời gian đo mưa trên lưu vực. - Lớp thông tin quan trắc về lưu lượng tại các vị trí kiểm tra và xuất số liệu trong mô hình. 1. Độ cao và độ dốc 1 Bản đồ DEM 50m 2. Mạng lưới sông 3. Phân chia tiểu lưu 1. Phân loại đất 2 Bản đồ Đất dạng số 2. Độ ẩm đất Bản đồ Thảm phủ 1. Lớp phủ thực vật 3 dạng số 2. Lượng trữ nước mặt Phân bố mưa theo 4 Mưa quan trắc không gian Mặt cắt ngang 5 Đo đạc lòng sông 6 Quan trắc lưu lượng vào Dòng chảy trong sông Hình 2.2: Các lớp thông tin yêu cầu của Marine
  20. 20 Hình 2.3: Cấu trúc các lớp dữ liệu trong Marine
  21. 21 Chương 3. PHẦN MỀM THỦY LỰC MỘT CHIỀU IMECH1D IMech1D là phần mềm thủy lực một chiều, được nghiên cứu, xây dựng và phát triển bởi tập thể nghiên cứu về lũ lụt của Viện Cơ học, Viện Khoa học và Công nghệ Việt Nam. IMech1D được xây dựng năm 2000, từ đó đến nay Imech1D luôn được nghiên cứu phát triển ngày một hoàn thiện hơn. IMech1D đã vượt qua được 12 bài toán kiểm định mẫu (Test cases) khắt khe đối với mô hình thủy lực một chiều. Từ năm 2002 đến nay IMech1D đã được chuyển giao cho một số cơ quan để triển khai nghiên cứu và ứng dụng, trong đó có “Cục phòng chống Lụt bão và Quản lý đê điều” của Bộ Nông nghiệp và Phát triển Nông thông, “Trung tâm Dự báo Khí tượng Thủy văn Trung ương” của Bộ Tài nguyên và Môi trường, Đặc biệt từ năm 2004 đến nay IMech1D được sử dụng để dự báo mực nước trên hệ thống sông Hồng và sông Thái Bình trong mùa lũ. IMech1D đã được đánh giá tốt qua nhiều hội đồng khoa học các cấp. 3.1. Các thành phần của hệ thống [5]; [6]; [9]: Để mô phỏng quá trình lan truyền lũ trong hệ thống sông trong IMech1D các thành phần chính được xem xét và mô phỏng lũ gồm có: - Mạng sông. - Ô ruộng. 3.1.1. Mạng sông: Mạng sông được xây dựng từ các đoạn sông và các nút. 3.1.1.1. Nút sông: Trong IMech1D có 2 loại nút sông. + Nút biên: Là vị trí tiếp xúc của hệ thống (ở đây là mạng sông) với các yếu tố bên ngoài của hệ thống. Giả thiết rằng dòng chảy trên mạng sông là dòng chảy êm, do vậy, tại mỗi nút biên sẽ cho 1 điều kiện biên. IMech1D sử dụng 3 loại điều kiện biên. - Cho mực nước (nút biên Z). - Cho lưu lượng (nút biên Q). - Cho quan hệ giữa mực nước và lưu lượng. + Nút trong của mạng (đơn giản gọi là nút): là vị trí tiếp xúc của từ 2 đoạn sông trở lên.
  22. 22 Nút sông chỉ có một đặc trưng duy nhất là cao trình mực nước tại nút đó. 3.1.1.2. Đoạn sông: Đoạn sông là mô hình của đoạn sông thực nằm giữa 2 nút sông. Vì vậy, đoạn sông có các đặc trưng sau: - Nút đầu đoạn. - Nút cuối đoạn. - Mặt cắt đầu đoạn. - Mặt cắt cuối đoạn. - Độ dài của đoạn - Các yếu tố thủy lực của đoạn như: độ dốc, hệ số nhám, lưu lượng bổ xung, Trong mặt cắt của đoạn có các đặc trưng như lòng sông, bãi sông, bối. Một đoạn có thể liên quan 1 hoặc 2 bối. 3.1.2. Ô ruộng (Ô chứa): Trong IMech1D ô ruộng được mô phỏng có vai trò chứa nước, chưa tính đến ảnh hưởng của vận tốc trên các ô ruộng đến các đoạn sông. Vì vậy, một ô ruộng được đặc trưng bởi các thông số: - Thể tích ô theo cao trình mực nước. - Các mối quan hệ giữa ô đang xét với các thành phần khác của hệ thống (thí dụ: trao đổi nước qua đập tràn, chiều cao, chiều rộng, hệ số của đập v.v ) 3.2. Mô hình toán học [5]; [6]; [9]: 3.2.1. Mô hình toán học một đoạn sông: Dòng chảy trong một đoạn sông được mô phỏng bằng hệ phương trình Saint Venant 1 chiều. Trong IMech1D hệ phương trình S.Venant được sử dụng dưới dạng sau: Ac Q q (3.1) t x Q  Q2 Z  gA Sf 0 (3.2) t x A x Ở đây sử dụng các ký hiệu: Q=Q(x,t) – lưu lượng của dòng chảy trong đoạn sông;
  23. 23 Z=Z(x,t) – mực nước trong đoạn sông; q - lưu lượng phụ; Ac - Diện tích mặt cắt (kể cả vùng chứa);  - hệ số điều chỉnh A - Diện tích chảy Sf - Sức cản đáy Sức cản đáy trong IMech1D được tính theo công thức sau: n 2Q Q Sf (3.3) A2R 4 /3 Trong đó R là bán kính thủy lực. 3.2.2. Mô hình toán học của một ô ruộng [5]; [6]; [9]: IMech1D mô phỏng quá trình ngập một ô ruộng (khu chứa) bằng định luật bảo toàn khối lượng nước tại ô đó. dV P Q  k (3.4) dt k Ở đây: V = V(Z) = V(Z(t)) thể tích của ô theo mực nước Z. P - Lượng mưa hoặc bốc hơi tại ô. Qk - Lượng nước trao đổi giữa ô đang xét với các ô liên quan (ô liên quan có thể là một đoạn sông, một nút sông, hoặc một ô ruộng khác). 3.3. Lược đồ sai phân [1]; [5]; [6]: Lược đồ sai phân Preissmann (xem [1]) được sử dụng để sai phân hóa hệ phương trình (3.1), (3.2), (3.4). Cụ thể là: f f n 1 f n f n 1 f n i 1 i 1 i i t 2 t 2 t n 1 n 1 n n f fi 1 fi fi 1 fi  1  , 0  1 x t x  n 1 n 1 1  n n f fi 1 fi fi 1 fi 2 2 Từ phương trình (3.1) ta có:
  24. 24 1 1  A A AT AT   Q Q 1  QT QT  2 t c,i 1 c,i c,i 1 ci x i 1 i i 1 i (3.5) 1  Q q 1  qT qT  2 i 1 i i 1 i Để sai phân hoá phương trình (3.2), ta viết phương trình này thành 4 số hạng a1, a2, a3, a4 như sau: a1 + a2 + a3 + a4 = 0 (3.6) trong đó Q  Q2 a , a  1 t 2 x A Z a gA , a gAS 3 x 4 f Sai phân từng số hạng trong (3.6) ta được 1 a ~a   Q Q QT QT  1 1 2 t i 1 1 i 1 i 2 2  Q 2 Q 2 QT QT a a~   i 1 i 1  i 1 i 2 2 x A A AT AT i 1 i i 1 i ~ g T T 1 T T a3 a3   Ai 1 Ai 1  Ai 1 Ai   Zi 1 Zi 1  Zi 1 Zi  2 x ~  1  T T a4 a4  g  AS f AS f  g  AS f AS f  2 i 1 i 2 i 1 i Do vậy phương trình sai phân của (3.2) có dạng sau: ~ ~ ~ ~ a1 a2 a3 a4 0 (3.7) Sai phân hóa phương trình (3.4) ta có 1 T T T Vi Vi  Pi Qi,k 1  Pi Qi,k (3.8) t k k Trong (3.5) - (3.8), f tính ở thời điểm hiện tại còn fT tính ở thời điểm trước đó. Ta giả thiết rằng nút đầu của đoạn sông có số thứ tự là i còn nút cuối đoạn có số thứ tự là i+1. Khi đó các đẳng thức (3.5) và (3.7) cho ta quan hệ lưu lượng và mực nước tại hai đầu của một đoạn sông. Như vậy, mỗi một đoạn sông có 4 ẩn số cần tìm là: Zd, Zc - mực nước ở đầu và ở cuối đoạn Qd, Qc - Lưu lượng ở đầu và ở cuối đoạn
  25. 25 Mỗi một ô ruộng có một phương trình (3.8). Phương trình này cho phép xác định một ẩn cần tìm của một ô ruộng là cao trình mực nước Zr tại tâm điểm của ô. Còn các hàm Vi, Qi, k trong (3.8) là hàm phi tuyến phụ thuộc vào cao trình mực nước Zr của ô ruộng đang xét và của các ô ruộng liên quan. Bây giờ ta xét tính đóng kín của hệ phương trình (3.5), (3.7) và (3.8) với các ẩn Zd, Zc, Qd, Qc và Zr. Mỗi một ẩn Zd hoặc Zc đều tương ứng với một giá trị cao trình mực nước tại một nút nào đó của mạng sông. Ta đặt: n - số nút của mạng sông (kể cả nút biên) ; nd - số đoạn sông của mạng sông ; nr - số ô ruộng. Tại các nút của mạng sông (nút đơn hoặc nút hợp lưu) giả thiết là chỉ có một cao trình mực nước. Như vậy tại n nút của mạng sông cần phải xác định n cao trình mực nước Zi : ta có n ẩn số. Vì mỗi đoạn sông có lưu lượng vào và lưu lượng ra khỏi đoạn, nên nếu ta có nd đoạn sông trong mạng sông thì ta cần xác định 2nd giá trị lưu lượng : ta có 2nd ẩn số. Nếu ta có nr ô ruộng ta cần xác định nr giá trị cao trình mực nước tại tâm các ô ruộng này : Như vậy để xác định dòng chảy trong mạng sông và ô ruộng ta cần xác định n+2nd+nr ẩn số. Các đẳng thức (3.5), (3.7). (3.8) cho ta 2nd+nr phương trình. Nếu nút là nút biên thì ta có một điều kiện biên. Nếu nút là nút trong thì ta có một phương trình cân bằng lưu lượng tại nút đó: Qvào - Qra = 0 (3.9) Do vậy tại n nút ta có n phương trình (hoặc dạng (3.9) hoặc là điều kiện biên). Cùng với 2nd + nr phương trình đã có, n phương trình này sẽ tạo thành hệ kín để xác định 2nd+nr+n ẩn số. 3.4. Tuyến tính hóa hệ phương trình (3.5), (3.7), (3.8):
  26. 26 Hệ phương trình phi tuyến (3.5), (3.7), (3.8) có thể giải được theo một phương pháp xấp xỉ. Trong IMech1D hệ phương trình phi tuyến (3.5), (3.7), (3.8) được đưa về hệ phương trình tuyến tính theo công thức Newton. Một hàm phi tuyến F(H, Q), tại bước n+1 có thể thay bằng một hàm tuyến tính dạng: F F F Hn 1,Qn 1 F H*,Q* Hn 1 H* Qn 1 Q* H H H* Q Q Q* Ở đây H*, Q* là giá trị lặp, giá trị xuất phát của chúng là Hn, Qn tương ứng. 3.4.1. Tuyến tính hoá các biểu thức đơn giản [10]: Ta xét các biểu thức đơn giản trong (3.5), (3.7), (3.8): 2 A = A(Z), Q /A, AZ, Ad Zc, AcZd. Chỉ số d là giá trị tính tại đầu đoạn, còn chỉ số c là giá trị ở cuối đoạn. Sử dụng công thức Newton để tuyến tính hóa các biểu thức này ta có: A A* + b*(Z - Z*) (3.10) 2 2 Q2 Q* Q* Q* 2 Q Q* b* Z Z* * * * (3.11) A A A A AZ A*Z* Z*b* A* Z Z* (3.12) * * * * * * * AdZc AdZc bdZc Zd Zd Ad Zc Zc (3.13) * * * * * * * AcZd AcZd Ac Zd Zd bcZd Zc Zc (3.14) Trong các công thức (3.10) - (3.14), b* - chiều rộng của sông (tương ứng với mực nước Z*). b* xuất hiện theo đẳng thức A b Z 3.4.2. Tuyến tính hoá biểu thức có lực cản đáy: ~ Trong số hạng a4 của (3.7) có chứa hàm phi tuyến ASf. Để tuyến tính hoá số hạng này, ta viết Sf về dạng sau: Q Q S , K kAR 2 / 3 f K 2 Hệ số K phụ thuộc vào mực nước Z: K = K(Z). Hàm này được tính qua diện tích chảy A = A(Z) và bán kính thủy lực R = R(Z). Vì các mặt cắt thực tế trên lưu vực đồng bằng Bắc bộ nói chung khá phức tạp, cho nên K phải tính theo kiểu lòng dẫn phức hợp. Hiện tại, IMech1D tính K theo từng cấp và giả thiết
  27. 27 rằng K biến đổi tuyến tính trong khoảng giữa hai cấp. Như vậy ASf có thể tuyến tính hoá như sau: Q* Q* Q* * * * * * * * ASf A Sf 2A 2 Q Q 2A 3 Kz Z Z K* K* 3.4.3. Tuyến tính hoá biểu thức trao đổi nước qua đê: Trong phương trình (3.8) có biểu thức phi tuyến Qi,k = F(Zi, Zk). Đây là hàm mô tả lượng nước trao đổi giữa ô i (có mực nước Zi) với ô k (có mực nước Zk). Việc trao đổi này còn phụ thuộc vào công trình giữa ô i và ô k. Giả sử rằng quá trình trao đổi nước giữa ô i với ô k tuân theo qui luật của dòng chảy tràn qua đập. Xét hàm F(Zi, Zk) khi giữa ô i và ô k là một đập tràn với cao độ Zd, độ dài 2 tràn là Ld. Khi đó, nếu Z - Z Z Z (chảy tự do từ i sang k). k d 3 i d 3/ 2 Qi, k F Z i, Z k mL d 2 g Z i Z d (3.15) 2 Nếu ZZZZ (chảy ngập) k d3 i d 1 / 2 Qi,k mLd Zk Zd Zi Zk (3.16) Với các hàm (3.15) Và (3.16) ta có thể sử dụng công thức Newton để tuyến tính hóa: - Trường hợp chảy tự do (3.15) 3 1/ 2 Qi,, k Q i k mL d2 g Z i Z d Z i Z i (3.17) 2 - Trường hợp chảy ngập (1.16) 1/ 21 1/ 2 QQmLZZikik,, dik ZZZZ idik ZZ kk 2 (3.18) 1 1/ 2 mLd Z i Z d Z i Z k Z i Z i 2 Hệ phương trình phi tuyến (3.5), (3.7), (3.8) có thể đưa về hệ phương trình đại số tuyến tính theo mực nước tại các nút sông và các ô ruộng và theo lưu lượng tại hai đầu của mỗi đoạn sông. 3.5. Thuật giải hệ phương trình đại số tuyến tính [1]; [5]; [6]: Sử dụng các hệ thức (3.10) – (3.18) từ hệ phương trình phi tuyến (3.5), (3.7), (3.8). (3.9) và điều kiện biên tại các nút biên ta có thể thu nhận được hệ phương trình tuyến tính để xác định n giá trị cao trình mực nước tại n nút, nr cao
  28. 28 trình mực nước tại r ô ruộng, nd giá trị lưu lượng vào và nd giá trị lưu lượng ra của nd đoạn sông. Ký hiệu Zd, Zc – cao trình mực nước tại đầu đoạn và cuối đoạn của đoạn sông, Qd, Qc – lưu lượng dòng chảy tại đầu đoạn (dòng chảy vào) và tại cuối đoạn (dòng chảy ra). Từ hệ phương trình (3.5) (3.7) có thể thu được hệ phương trình tuyến tính dạng: a1Qd + a2Qc + a3Zd + a4Zc = a5 (3.19) b1Qd + b2Qc + b3Zd + b4Zc = b5 (3.20) Từ (3.19), (3.20) dễ dàng suy ra: Qd = a6Zd + a7Zc + a8 (3.21) Qc = b6Zd + b7Zc + b8 (3.22) Điều kiện cân bằng lưu lượng (3.9) tại các nút trong của mạng sông có thể viết lại dưới dạng: i j Qc Qd 0 (3.23) i j Ở đây ký hiệu i là chỉ số chỉ đoạn sông có lưu lượng chảy vào nút và j là chỉ số chỉ đoạn sông có lưu lượng chảy từ nút ra. Phương trình bảo toàn khối lượng nước trong ô ruộng (3.8) có thể viết dưới dạng T T Zir Zir T T Fi  Pi Qik 1  Pi Qik (3.24) t k k Ta thấy phương trình (3.24) là quan hệ giữa cao trình mực nước tại ô ruộng thứ i với cao trình mực nước tại các nút hoặc ô ruộng liên quan bên cạnh. Nếu tại các nút biên thượng lưu i cho giá trị lưu lượng vào Qi Qdi = Qi (3.25) Tại các nút biên hạ lưu j cho giá trị cao trình mực nước: Zcj = Zj (3.26) Hoặc tại nút biên k cho quan hệ giữa lưu lượng và mực nước Zk = Fk(Qk) (3.27) (Tùy nút biên là nút biên thượng lưu hoặc nút biên hạ lưu mà Zk, Qk có thể chỉ số Zkc, Zkd, Qkc, Qkd)
  29. 29 Sử dụng (3.21) (3.22) từ các phương trình (3.23) – (3.27) ta có thể thu được hệ phương trình tuyến tính n+nr bậc để xác định n cao trình mực nước tại n nút và nr cao trình mực nước tại tâm của r ô ruộng. AZ = B (3.28) ở đây A là ma trận vuông có bậc bằng n+nr. Để tính toán cho 1 cơn lũ, hệ (3.28) sẽ được giải nhiều lần trong một bước thời gian cũng như cho các bước thời gian khác nhau, vì vậy việc chọn phương pháp để giải nhanh hệ trên rất có ý nghĩa trong việc rút ngắn thời gian tính. Để ý rằng, nếu hệ có nhiều nút hợp lưu thì (3.28) là hệ thưa nhưng các hệ số phân bố khá phức tạp,trên một hàng có thể có nhiều phần tử. Nếu chỉ có 1 nhánh sông thì ta được hệ 3 đường chéo và có thể dùng phương pháp khử đuổi (còn gọi là thuật toán Thomas, xem [2]) để giải. Ở trường hợp chung, ý tưởng sẽ là khử đuổi một số nút và khử Gauxơ các nút còn lại. Vấn đề là xác định số nút cần khử đuổi và thứ tự khử chúng. Trước hết ta xét sự thay đổi cấu trúc mạng khi khử đuổi (hình 3.1). Ta quy ước gọi bậc của 1 nút, ký hiệu r, là số đoạn nối với nút đó.Ví dụ ở hình 3.1, nút 2 có bậc là 3, nút 7 có bậc là 5. Dễ thấy các nút biên có bậc là 1,các nút đơn (không phải là nút hợp lưu hoặc biên) có bậc là 2. Các nút hợp lưu có bậc từ 3 trở lên. Khi khử 1 ẩn Zk nào đó, về mặt hình thức, hệ sẽ bớt 1 phương trình và ẩn Zk sẽ không có mặt trong các phương trình liên quan nữa. Đối với mạng sông, điều đó có thể xem như bỏ đi nút k tương ứng và khi đó, các nút kề nút k trước đây sẽ được nối với nhau. Như vậy khi khử, ta bỏ được một số liên kết (đoạn), nhưng cũng phải thêm vào một số liên kết (đoạn) khác. Ta ký hiệu sdb=số đoạn bớt đi; sdt=số đoạn thêm vào; dd=sdt-sdb. Để ý rằng dd có thể âm. Đối với nút đơn cũng như nút biên, dd=-1. Bây giờ ta có thể đánh giá số phép tính khi khử đuổi như sau: Giả sử n-số nút của toàn mạng, r- bậc của nút i. Do phương trình của nút i 2 chứa (r +1) số hạng nên để khử đuổi biến Zi cần (r+1) phép tính (chỉ tính phép nhân hoặc chia) ở bước khử xuôi và r phép tính ở bước khử ngược, tức cỡ r2 phép tính. Nếu khử đuổi l biến, bậc không quá R thì số phép tính sẽ là cỡ lR2 . Để ý rằng khử Gauxơ có số phép tính là n2l (xem [14]). Như vậy nếu R n thì số phép tính theo phương pháp khử đuổi giảm đáng kể.
  30. 30 Một số tác giả khử đuổi các ẩn nằm trong mỗi nhánh sông, kết quả đưa về hệ chứa các ẩn là nút hợp lưu. Tuy nhiên nếu hệ có nhiều nút hợp lưu, ví dụ khi cần xét cả các ô chứa ở ngoài đê, thì hệ sau khi đã khử đuổi vẫn còn khá lớn. Ngoài ra vấn đề trình tự khử đuổi cũng không đơn giản, nếu mạng sông được đánh số tuỳ ý. Thông thường người ta dựa vào cảm nhận của mình để đưa ra trình tự khử cho máy tính. Ở đây chúng tôi đưa ra một thuật toán xác định nút khử đuổi cũng như trình tự khử cho hệ (3.28). Thuật toán này cho phép lập trình để máy tính tự động xác định trình tự khử. Như phân tích ở trên, bậc của nút càng nhỏ thì số phép tính để thực hiện khử nó càng ít, vì vậy cần ưu tiên nút có bậc thấp để khử. Cụ thể thuật toán để chọn nút khử như sau: - Tìm các nút có bậc thấp nhất. - Trong số các nút có bậc thấp nhất, tìm nút có dd (hiệu số giữa số đoạn thêm và số đoạn bớt) bé nhất, đó là nút cần khử tiếp theo. Theo thuật toán này, các nút ở nhánh cụt sẽ được khử trước, bắt đầu từ nút biên (tại thời điểm khử, chúng có bậc là 1). Tiếp theo là các nút đơn nằm giữa 2 nút hợp lưu (tại thời điểm khử, chúng có bậc 2). Tiếp theo là các nút hợp lưu bậc 3, bậc 4 Để minh hoạ, ta xét ví dụ :giả sử sau khi khử các nhánh cụt ta còn mạng như hình 3.1. Áp dụng thuật toán trên ta khử tiếp như sau: - Dễ thấy bậc thấp nhất ở hình 3.1 là 3. Các nút có bậc 3 là nút 2, 6, 10, 9, 13, 17. Nếu khử nút 2 thì sdb=(2-4); (2-6); (2-7) và sdt=(4-6) khi đó dd=-2. Có thể thấy các nút còn lại đều có dd=-1.Vậy ta chọn khử nút 2. Sau khi khử, ta được hình 3.2. - Bậc thấp nhất vẫn là 3. Các nút có bậc 3 là:17,6,10,13,9. Nếu khử nút 17 thì:sdb:(17-12);(17-13);(17-15) sdt: (13-15);dd=-2. Dễ thấy các nút bậc 3 còn lại đều có dd=-1 Vậy ta chọn khử nút 17. Ta được hình 3.3 2 4 11 4 8 4 10 17 4 9 8 6 13 6 9 10 8 12 12 2 12 10 6 13 13 9 17 13 9 4 4 17 Hì Hì Hì Hình 3.1 11 Hình 3.3
  31. 31 Lý luận tương tự, các nút khử tiếp theo như sau: - Khử nút 13 (bậc 3): sdb: (13-9); (13-12); (13-15) sdt: (9-12); (9-15); dd=-1. Ta được hình 3.4 - Khử nút 10 (bậc 3): sdb=(10-6); (10-11); (10-15) sdt=(6-11); (6-15); dd=-1. Ta được hình 3.5 - Khử nút 11 (bậc 4): sdb=(11-6); (11-7); (11-12); (11-15) 15 sdt=(6-12); (7-12); (7-15) 7 khi đó dd=-1.Ta được hình 3.6 15 4 4 4 4 10 4 9 11 12 6 9 6 9 6 4 8 6 7 9 9 11 8 8 7 15 12 Hình 3.5 - Khử nút 6 (bậc 4): sdb=(6-4); (6-7); (6-12); (6-15) sdt=(4-12); (4-15); dd=-2. Ta được hình 3.7 - Khử nút 9 (bậc 4): sdb=(9-4); (9-8); (9-12); (9-15) sdt=(8-15); dd=-3. Ta được hình 3.8 - Khử nút 15 (bậc 4): sdb=(15-4); (15-7); (15-8); (15-12) sdt=0, dd=-4. Ta đư15ợc hình 3.9 12 7 12 4 4 8 4 4 8 9 7 7 15 12 8 8 7 12 8 Hình 3.9 Khử nút 7 (bậc 3): sdb=(7-4); (7-8); (7-12)
  32. 32 sdt=0, dd=-3.Ta được hình 3.10. - Khử nốt nút 4, 12, ta chỉ còn 1 nút là nút 8. Qua ví dụ trên ta thấy, không phải càng khử đuổi, 4 bậc của các nút càng tăng, ví dụ ở hình 3.10, sau khi khử nút 7, bậc của nút 12 giảm 1 đơn vị. Mặc dù thuật toán trên đòi hỏi khá nhiều phép xử 8 lý của máy tính, tuy nhiên chỉ cần xác định 1 lần cho tất 12 cả các lần giải hệ (3.28) mà thôi. Hình 1.10 Hiện nay với số liệu địa hình hiện có, mô hình thủy lực 1 chiều mở rộng gồm mạng sông Hồng – Thái Bình và các ô ruộng có thể mô phỏng dưới dạng mạng có 2093 đoạn sông, 1103 nút và 278 ô ruộng. Thời gian tính cho cơn lũ tháng 8/1971 theo các phương án khử khác nhau cho trong bảng 3.1. Kết quả Số nút Số nút Thời Phương án khử đuổi khử Gauxơ gian tính Khử đuổi nút có bậc 2 404 699 2h50’ Khử đuổi nút có bậc 5 960 143 33’54’’ Khử đuổi nút có bậc 10 1053 50 32’50’’ Bảng 3.1: Thời gian tính của IMech1D Như vậy ta thấy, so với phương án chỉ khử đuổi các nút biên và nút đơn (bậc 2) thì phương án khử đuổi cả các nút hợp lưu làm giảm thời gian tính gần 7 lần. Mặt khác thực tế cho thấy, khi (số nút khử Gauxơ)/(bậc của nút khử đuổi) <10 thì việc tăng tiếp nút khử đuổi là ít có ý nghĩa. Sở dĩ như vậy vì khi khử đuổi, tuy được lợi về số phép tính nhưng do phải lưu lại vị trí của các phần tử nên thời gian xử lý có lâu hơn. 3.6. Các thuật toán phụ trợ sử dụng trong xây dựng bộ chương trình tính toán thủy lực một chiều IMech1D [5]; [6]: yi 3.6.1. Khái toán mặt cắt: Vì địa hình mặt cắt của hệ thống sông thường rất phức tạp zi Hình 3.11
  33. 33 nên việc đưa mặt cắt về dạng hình thang đối xứng sẽ gây ra sai số lớn. Để tính toán các đặc trưng địa hình mặt cắt sông trong IMech1D chia mặt cắt theo các đường song song với trục thẳng đứng z (yi = const) (hình 3.11) Các thông số địa hình của mặt cắt ngang được tính toán theo các công thức sau: Chiều rộng của mặt cắt: b bi i Diện tích mặt cắt: A  Ai i 2 / 3 Môđul lưu lượng: K  K i  ki Ai Ri Ai Ai Bán kính thủy lực: Ri 2 2  i yi zi Trong đó : i – chu vi ướt. Khác với cách xử lý của một số các tác giả trước, trong IMech1D mỗi đoạn sông được đặc trưng bởi hai mặt cắt: mặt cắt đầu và mặt cắt cuối. Trong trường hợp không có hợp lưu thì mặt cắt cuối của đoạn trước sẽ trùng với mặt cắt đầu của đoạn sau. 3.6.2. Tạo giá trị mực nước và lưu lượng làm điều kiện ban đầu: Khi tính thuỷ lực của mạng sông thực, thường không biết trước giá trị đầu về mực nước z và lưu lượng Q tại các nút, chương trình phải tự tạo. Theo lý thuyết, giá trị đầu không cần chính xác vì chỉ ảnh hưởng đến dòng chảy trong 1 thời gian ngắn sau đó. Tuy nhiên, trên thực tế, nếu giá trị đầu chọn không thích hợp, chương trình sẽ không chạy được. Vì thời điểm bắt đầu tính thường là khi chưa có lũ, dòng chảy có thể xem là dòng dừng, vì vậy giá trị đầu chính là nghiệm dừng của hệ phương trình Saint- Venant. Điều kiện biên được lấy bằng giá trị biên tại t=t0. Qbie Zbien Q0 Z0 t t Hình 3.12: Điều kiện đầu trong IMech1D
  34. 34 3.6.3. Vấn đề xác định hệ số nhám và chỉnh kết quả: Vì các hệ số nhám của sông là khó xác định chính xác trước, các chương trình tính toán thuỷ lực phải tìm cách hiệu chỉnh sao cho kết quả tính toán phù hợp với thực đo. Cách làm ở đây là rất khác nhau. Ở một số chương trình, người ta đưa ra các đường cong nhám với các giá trị nhám khác nhau theo 12 cấp (thực ra là giá trị hệ số điều chỉnh nhám). Sau đó mỗi đoạn sông được gán cho 1 đường. Khi kết quả tính toán lệch so với thực đo, người ta chọn lại đường nhám hoặc thay đổi chính bản thân nó. Việc này là khá khó khăn. Cách làm trong IMech1D như sau: - Chia mạng sông thành nhiều vùng, mỗi vùng được gán với 1 đường cong nhám (hiện tại mạng sông Hồng - Thái bình được chia thành 16 vùng). - Mỗi đường cong nhám có dạng (z1,g1); (z2,g2) ; (zm,gm) - Trong đó zi là mực nước của 1 nút chuẩn nằm trong vùng đang xét, z1<z2 zm, gi là giá trị điều chỉnh. Thường lấy nút chuẩn là các trạm thủy văn. - Số cặp giá trị m là tùy ý và khác nhau đối với mỗi đường. - Nếu mực nước z của nút đang xét nằm trong khoảng (zi,zi+1) thì g được nội suy từ gi, gi+1. - Khi cần chỉnh, ta dựa vào độ lệch giữa kết quả tính toán và thực đo của nút chuẩn (được vẽ trên màn hình), tại mực nước đang xét, để chỉnh giá trị gi. Có thể chèn thêm cặp (zk, gk) nếu cần. Các cách tiếp cận trên đã được thử nghiệm trên phần mềm để tính toán các bài toán mẫu cũng như các cơn lũ thực trên hệ thống sông Hồng và sông Thái Bình. Các kết quả là khả quan.
  35. 35 Chương 4. ỨNG DỤNG MÔ HÌNH KẾT NỐI MARINE- IMECH1D CHO LƯU VỰC SÔNG ĐÀ PHẦN TRÊN LÃNH THỔ VIỆT NAM 4.1. Mô hình kết nối Marine và Imech1D [5]: Trong mô hình kết nối Marine và Imech1D, Marine có vai trò là mô hình thu gom nước trên bề mặt lưu vực từ mưa và tập trung lượng nước này tại các vị trí ven hai bên bờ sông trong mạng sông tính toán của mô hình thủy lực một chiều IMech1D. Mô hình thủy lực IMech1D có vai trò tiếp nhận lưu lượng nước vào tại vị trí biên giới như là một biên lưu lượng vào ở vị trí biên trên của mô hình thủy lực đồng thời IMech1D có nhiệm vụ thu gom lượng nước do Marine tập trung dọc hai bên bờ sông như là lưu lượng phụ ra nhập dòng bên trong mô hình thủy lực. MailSông river chính Hình 4.1: Kiểu kết nối nhiều điểm giữa Marine và IMech1D Có hai loại kết nối giữa Marine và IMech1D đó là Kết nối một điểm và kết nối nhiều điểm. Với những lưu vực không co sông (trong mạng sông của mô
  36. 36 hình thủy lực) chảy qua thì sử dụng loại kết nối một điểm, điểm kết nối giữa hai mô hình chính là điểm xuất nước của Marine, đó chính là cửa ra của lưu vực, kiểu kết nối này thường được sử dụng cho lưu vực khép kín hay còn gọi là lưu vực đóng. Với những lưu vực có sông trong mạng sông của mô hình thủy lực chạy qua thì sử dụng loại kết nối nhiều điểm, kiểu kết nối này được sử dụng cho các lưu vực mở, tức là lưu vực không khép kín, trong trường hợp này vị trí kết nối là tất cả các điểm sông trong mạng thủy lực của mô hình thủy lực nằm trong lưu vực. Hai kiểu kết nối được mô tả trên các hình 4.1 và hình 4.2. SôngMail river chính Hình 4.2: Kiểu kết nối một điểm giữa Marine và IMech1D 4.2. Xử lý số liệu cho mô hình kết nối Marine-IMech1D: 4.2.1. Xử lý bản đồ địa hình: Bản đồ địa hình sử dụng làm bản đồ nền cho Marine la bản đồ số dưới dạng các đường đồng mức và điểm độ cao có tỷ lệ là 1:50000. Từ bản đồ này qua quá
  37. 37 trình sử lý trên các phần mềm GIS thu được bản đồ số độ cao DEM của toàn lưu vực sông Đà. Kích thước của DEM là 50x50. Hình 4.3: DEM của toàn lưu bộ lưu vực sông Đà trên lãnh thổ Việt Nam 4.2.1.1. Xác định hướng của dòng chảy và độ tích tụ của dòng chảy trên DEM: Trên nền bản đồ số độ cao (DEM) sử dụng các công cụ GIS để xác định hướng chảy của dòng trên bề mặt lưu vực. Nguyên lý xác định hướng chảy được trình trong hình 4.4 dưới đây: 2 2 4 4 8 1 2 4 8 4 1 1 2 4 8 2 1 4 4 4 1 1 1 2 1 6 Hình 4.4: Nguyên lý xác định hướng của dòng chảy Sau quá trình sử lý trên nền GIS thu được lớp bản đồ xác định hướng chảy của mỗi ô lưới. Mỗi một ô lưới được xác định chỉ có một hướng chảy cố định
  38. 38 theo địa hình. Bản đồ hướng chảy của lưu vực sông Đà được trình bày trên hình 4.5 dưới đây. Hình 4.5: Bản đồ hướng chảy của lưu vực sông Đà Tương tự như xây dựng bản đồ hướng chảy, từ DEM sử dụng các công cụ GIS có thể xác định được hướng tích tụ của dòng chảy. Nguyên lý xác định hướng tích tụ dòng chảy được trình bày trong hình 4.6 và kết quả thu được bản đồ độ tích tụ của dòng chảy trên lưu vực sông Đà trên hình 4.7. Hình 4.6: Nguyên lý xác định độ tích tụ của dòng chảy
  39. 39 Hình 4.7: Bản đồ độ tích tụ của dòng chảy 4.2.1.2. Tạo mạng sông suối từ DEM: Từ hai lớp bản đồ “hướng dòng chảy” và “độ tích tụ dòng chảy” đã xây dựng ở trên xác định được mạng sông suối trên nền DEM.
  40. 40 Hình 4.8: Mạng sông suối được xác định trên nền DEM 4.2.1.3. Phân chia lưu vực trên nền DEM: Cũng từ hai lớp dữ liệu là “hướng dòng chảy” và hướng hội tụ của dòng chảy đã xác định ở trên cho phép xác định được đường phân nước của lưu vực. Hình 4.9: Vị trí của 10 lưu vực bộ phận Qua phân tích địa hình, mạng sông và mạng lưới trạm quan trắc trên lưu vực sông Đà, toàn lưu vực được chia nhỏ ra thành 10 lưu vực bộ phận để tính toán dòng chảy khu giữa từ mưa trong đó có 6 tiểu lưu vực chứa sông chính và 4 tiểu lưu vực chứa sông nhánh gồm: 1. Tiểu lưu vực Mường Tè: 2. Tiểu lưu vực Nậm Pô: chứa sông nhánh Nậm Pô 3. Tiểu lưu vực Nậm Giàng: chứa sông nhánh Nậm Na 4. Tiểu lưu vực Lai Châu 5. Tiểu lưu vực Nậm Mức: chứa sông nhánh Nậm Mức 6. Tiểu lưu vực Quỳnh Nhai 7. Tiểu lưu vực Bản Củng: chứa sông nhánh Nậm Mu 8. Tiểu lưu vực Tạ Bú 9. Tiểu lưu vực Sơn La 10. Tiểu lưu vực Hòa Bình 4.2.2. Xử lý bản đồ phân loại đất:
  41. 41 Bản đồ phân loại đất được dùng trong MARINE cho lưu vực sông Đà có tỷ lệ 1:1.000.000. Ban đầu bản đồ ở dạng Vector, sau đó dùng phần mềm ILWIS 3.1 đưa về dạng Raster. Dạng Raster chính là dạng có cấu trúc giống với cấu trúc của DEM. Toàn bộ bản đồ được đưa về dạng các ô lưới (50mx50m) thay vì các vùng đa giác khép kín có kích thước và hình dáng bất kỳ (dạng Vector). Toàn bộ vùng sông Đà được chia làm 29 loại đất phân bố trên 10 tiểu lưu vực. Do không có đầy đủ các thông số thực nghiệm về các loại đất nên các loại đất được gộp lại với nhau nếu chúng có cùng một số tính chất giống nhau. Như vậy, toàn lưu vực sông Đà nhóm gọn lại làm 4 nhóm đất theo tên của Mỹ gồm: ID (chỉ số loại đất)=3 đất cát mùn; ID=4 đất mùn; ID=5 đất phù sa mùn; ID=12 đất sét. Hình 4.10: Bản đồ phân loại đất của lưu vực sông Đà Sau khi có bản đồ Raster đất toàn lưu vực sông Đà tiến hành cắt tách thành 10 bản đồ đất của từng tiểu lưu vực (xem chi tiết trong Phần Phụ lục). Sử dụng phần mềm ArcView với phần mở rộng XTool và Spatial Analyst xuất bản đồ đất từng tiểu lưu vực ra dạng file ASCII tạo đầu vào cho MARINE.
  42. 42 4.2.3. Xử lý bản đồ hiện trạng sử dụng đất: Bản đồ hiện trạng sử dụng đất của lưu vực sông Đà có tỷ lệ 1:100.000. Ban đầu bản đồ ở dạng Vector cấu tạo bởi các vùng khép kín - polygon, sau đó được đưa về dạng Raster bằng phần mềm ILWIS 3.1. Toàn bộ bản đồ sử dụng đất vùng sông Đà được chia làm 13 loại thảm phủ phân bố trên 10 tiểu lưu vực đã phân chia. Lớp thảm phủ có thông số tính toán trong mô hình riêng tương ứng với mỗi loại. Hình 4.11: Bản đồ hiện trạng sử dụng đất của lưu vực sông Đà Tiếp theo, bản đồ Raster thảm phủ lưu vực sông Đà được cắt tách riêng cho từng tiểu lưu vực nhờ sự hỗ trợ của phần mềm ArcView với phần mở rộng XTool và Spatial Analyst. Khi đã có thảm phủ 10 tiểu lưu vực cắt tách sử dụng tiến hành xuất bản đồ thảm phủ từng tiểu lưu vực ra dạng file ASCII là đầu vào cho MARINE. 4.2.4. Xây dựng bản đồ phân bố mưa trong lưu vực:
  43. 43 Toàn lưu vực sông Đà sử dụng 36 trạm mưa khí tượng và mưa thuỷ văn bảng (4.1). Tích hợp các trạm khí tượng thủy văn lên bản đồ lưu vực sông Đà. Bản đồ phân bố vùng ảnh hưởng mưa được thực hiện theo lý thuyết đa giác Thái Sơn. Nối các trạm đo mưa thành đoạn thẳng chia lưu vực thành lưới các tam giác. Vẽ đường trung trực đi qua các cạnh của các tam giác. Các đường này tạo nên vùng diện tích ảnh hưởng của từng trạm mưa. Kết quả bản đồ phân vùng ảnh hưởng mưa thu được như hình 4.12. ID Tên Trạm Tỉnh Kinh độ Vĩ độ ID Tên Trạm Tỉnh Kinh độ Vĩ độ 1 Tam Đường Lai Châu 103 0 29’ 22 0 25’ 19 Bản Chiềng Hòa Bình 104 0 46’ 20 0 54’ 2 Muờng Tè Lai Châu 102 0 50’ 22 0 22’ 20 Bản Củng Yên Bái 103 0 49’ 21 0 47’ 3 Sin Hồ Lai Châu 103 0 15’ 22 0 21’ 21 Nậm Giàng Lai Châu 103 0 10’ 22 0 15’ 4 Binh Lư Lai Châu 103 0 37’ 22 0 18’ 22 Muờng Tè TV Lai Châu 102 0 37’ 22 0 28’ 5 Lai Châu Lai Châu 103 0 9’ 22 0 3’ 23 Km22 Hòa Bình 105 0 5’ 20 0 45’ 6 Tuần Giáo Lai Châu 103 0 25’ 21 0 35’ 24 Muờng Nhé Lai Châu 102 0 29’ 22 0 11’ 7 Quỳnh Nhai Sơn La 103 0 34’ 21 0 50’ 25 Nậm Pô Lai Châu 102 0 37’ 22 0 7’ 8 Sơn La Sơn La 103 0 54’ 21 0 20’ 26 Vàng Pó Lai Châu 103 0 28’ 22 0 36’ 9 Phù Yên Sơn La 104 0 39’ 21 0 16’ 27 Pa Tần Lai Châu 103 0 12’ 22 0 28’ 10 Cò Nòi Sơn La 104 0 9’ 21 0 8’ 28 Tà Hộc Hòa Bình 104 0 24’ 21 0 11’ 11 Yên Châu Sơn La 104 0 17’ 21 0 3’ 29 Vạn Yên Hòa Bình 104 0 38’ 21 0 1’ 12 Mộc Châu Sơn La 104 0 38’ 20 0 51’ 30 Nậm Mức Lai Châu 103 0 18’ 21 0 49’ 13 Hòa Bình Hòa Bình 105 0 20’ 20 0 49’ 31 Pak Ma Lai Châu 102 0 27’ 22 0 33’ 14 Than Uyên Hòa Bình 103 0 54’ 21 0 56’ 32 Lai Châu TV Lai Châu 103 0 12’ 22 0 4’ 15 Tạ Bú Sơn La 104 0 3’ 21 0 26’ 33 Thuận Châu Sơn La 103 0 42’ 21 0 25’ 16 Muờng Trại Sơn La 103 0 57’ 21 0 36’ 34 Mai sơn Hòa Bình 104 0 2’ 21 0 12’ 17 Muờng Sại Sơn La 103 0 42’ 21 0 33’ 35 Quỳnh Nhai TV Sơn La 103 0 33’ 21 0 51’ 18 Km46 Sơn La 104 0 59’ 20 0 47’ 36 Tà Nàng Hòa Bình 104 0 27’ 20 0 56’ Bảng 4.1: Vị trí các trạm mưa sử dụng trên sông Đà
  44. 44 Hình 4.12: Bản đồphân bố mưa trên lưu vực sông Đà 36 trạm mưa ứng với 36 vùng ảnh hưởng được đánh số theo thứ tự từ 1 đến 36. Bản đồ phân bố mưa được chuyển sang dạng Raster và xuất ra file dạng ASCII tạo đầu vào cho mô hình MARINE. 4.2.5. Tích hợp các mặt cắt sông vào lớp sông suối trên nền DEM: Tổng số mặt cắt trên sông Đà từ Mường Tè-Hoà Bình là 121 (đo năm 1997 do Viện cơ học cấp) được chia thành 2 đoạn: + Đoạn 1: Biên giới-Pa Vinh: 59 mặt cắt + Đoạn 2: Pa Vinh-trước đập Hoà Bình: 62 mặt cắt Từ số liệu trắc dọc của các mặt cắt sử dụng phần mềm Mapinfo 7.0 tích hợp các mặt cắt lên dòng chính sông Đà. Khoảng cách giữa các mặt cắt được đăng ký trong file *.rug chính là điểm kết xuất kết quả tính toán dòng chảy gia nhập khu giữa ở hai mặt cắt tại khoảng cách đó hình 4.13.
  45. 45 SDA_TPV1 SDA_TPV14 SDA_TPV32 SDA_PVHB10 SDA_TPV48 Hình 4.12: Các mặt cắt của sông đã được tích hợp lên nền DEM 4.2.6. Xử lý các số liệu khác: Các số liệu khác dùng cho mô hình kết nối bao gồm các số liệu về lưu lượng thực đo tại các trạm kiểm tra, vị trí các trạm kiểm tra, các điều kiện đầu của Marine và IMech1D cũng đã được xử lý theo yêu cầu của hai mô hình.
  46. 46 Chương 5. NÂNG CAO CHẤT LƯỢNG TÍNH TOÁN CỦA MÔ HÌNH BẰNG KỸ THUẬT LỌC KALMAN Phương pháp lọc Kalman là một trong các công cụ toán học nổi tiếng, được sử dụng để đánh giá ngẫu nhiên từ các kết quả nhiễu đo cảm biến. Cơ sở của phương pháp này là điều chỉnh kết quả tính toán của mô hình dựa trên các kết quả thực đo đã có tại một chuỗi thời gian. Phương pháp lọc Kalman được bắt đầu phát triển vào những năm 1960 bởi nhà thống kê R.E. Kalman, đây là một công cụ được sử dụng khá phổ biến trong thống kê toán học và lý thuyết hiệu chỉnh. Ngày nay phương pháp lọc Kalman đã được nghiên cứu rộng rãi trong bài toán xử lý sau mô hình tại nhiều trung tâm tính toán và dự báo lớn như Trung tâm Dự báo Môi trường Quốc gia NCEP (Cộng hòa Pháp), Trung tâm Dự báo hạn vừa châu Âu ECMWF, Viện Tin học tính toán ứng dụng Grenoble INRIA (Cộng hòa Pháp) [13], . . . Ưu điểm của phương pháp này là các hệ số được cập nhật liên tục theo thời gian dựa trên việc phân tích bản chất sai số dự báo của một vài chu kỳ dự báo trước. Phương pháp lọc Kalman được nghiên cứu và tùy biến sao cho phù hợp với các bài toán thuộc các lĩnh vực khác nhau. Trong luận văn này trình bày trường hợp nghiên cứu ứng dụng phương pháp lọc Kalman cho các mô hình thủy văn và mô hình thủy lực, đặc biệt là sử dụng để hiệu chỉnh kết quả tính đầu ra và kiểm tra ngược lại số liệu đầu vào của các mô hinh. 5.1. Quá trình cần đánh giá (ước lượng) [13]: Lọc Kalman rời rạc là phương pháp nhằm đánh giá trạng thái x n của quá trình điều khiển theo thời gian rời rạc tuân theo phương trình sai phân ngẫu nhiên tuyến tính xk A xk 1 B uk wk 1 (5.1) Với đại lượng đo đạc z m thoả mãn phương trình sai phân sau zk H xk  k (5.2) Ở đây: - xk là vectơ tham số trạng thái n chiều và là đại lượng ngẫu nhiên - zk là vec tơ đại lượng đo đạc m chiều và có giá trị ngẫu nhiên
  47. 47 - wk là vectơ ồn n chiều của quá trình điều khiển, có các tính chất: ồn trắng, phân bố chuẩn, giá trị trung bình thống kê bằng 0, hiệp phương sai Q. Nghĩa là p (w) ~ N (0, Q )  E[w] 0 Q = Cov w = E(w1 1)(w1 1) E(w1 1)(wn n ) Q  ij E(wn n )(w1 1) E(wn n )(wn n ) Vì i 0 suy ra Ew1w1 Ew1wn  Q  Ew w  Q = ij i j (5.3) Ewnw1 Ewnwn  -  k là vectơ ồn (nhiễu) m chiều của quá trình đo đạc, có các tính chất: ồn trắng, phân bố chuẩn, giá trị trung bình thống kê bằng 0, hiệp phương sai R p ( ) ~ N (0, R ) E[v] = 0 E11 E1 m  R = R  E[  ]  (5.4) ij i j E m1 E m m  Thông thường các sai số quan trắc không tương quan với nhau (độc lập thống kê), ta có E11 0 0 0 E 2 2  0 R  E[  ]  ij i j 0 0 E m m 
  48. 48 Hàm phân bố sác xuất có dạng 1 1 p( ) exp ( )T R 1 ( ) (5.5) n n 2 (2 ) 2 R 2 l - uk là vectơ driven function l chiều, đại lượng tiền định, u  - A ma trận [n x n]: quan hệ giữa trạng thái tại thời điểm trước (k-1) với trạng thái hiện thời (k). Đây là mô hình chuyển dịch trạng thái tác động vào trạng thái trước. - B ma trận [n x l]: quan hệ giữa trạng thái hiện thời với driven function cũng thời điểm hiện tại. Đây là mô hình điều khiển-đầu vào (control- input model) tác dụng vào vectơ điều khiển. - H ma trận [m x n]: quan hệ giữa trạng thái hiện thời với đo đạc cũng tại thời điểm hiện tại. Đây là mô hình ánh sạ tham số trạng thái vào quan trắc. 5.2. Các vấn đề tính toán (bản chất tính toán) của lọc Kalman [13]; [17]: 5.2.1. Định nghĩa về các ước lượng tiên nghiệm và hậu nghiệm: f n xk  là ước lượng (đánh giá) tiên nghiệm (priori estimate) của trạng thái tại bước k khi biết trạng thái ở thời điểm trước k (thí dụ như a xk 1 mà những ước lượng này đã sử dụng các giá trị đo đạc tại thời điểm trước k). Còn được gọi là giá trị dự báo. a n xk  là ước lượng (đánh giá) hậu nghiệm (posteriori estimate) của trạng thái tại bước k khi biết zk – giá trị đo đạc tức thời tại điểm k. Còn được gọi là giá trị hiệu chỉnh. 5.2.2. Bước dự báo – cập nhật (ước lượng tiên nghiệm): Sai số ước lượng tiên nghiệm f f ek  xk xk (5.6) Đây là đại lượng ngẫu nhiên có giá trị trung bình thống kê bằng không. Chứng minh: Trước hết theo định nghĩa ta có: f xk  E[xk ]
  49. 49 Do đó f f f E[ek ] E[xk xk ] E[xk ] xk 0 Ma trận hiệp phương sai của sai số ước lượng tiên nghiệm xác định như sau f f f f T Pk Cov ek E ek (ek )  Chứng minh: Trước hết theo định nghĩa và với kỳ vọng bằng không, ma trận hiệp phương sai có dạng: f f f f T f 2 f f T Pk Covek E ek (ek )  (E[ek ]) E ek (ek )  (5.7) Tìm ước lượng tiên mghiệm Lấy kỳ vọng của đại lượng sau xk A xk 1 B uk wk 1 Ta được: f a E[xk ] xk A xk 1 Buk (5.8) Tìm ma trận hiệp phương sai của sai số ước lượng tiên nghiệm Ta có sai số ước lượng sau f f a ek xk xk A(xk 1 xk 1) wk 1 Từ đó tìm được ma trận hiệp phương sai như sau (chú ý tới tính độc lập thống kê giữa các đại lượng ngẫu nhiên x và w ) f a a T T Pk E[A(xk 1 xk 1)(A(xk 1 xk 1)) wk 1 w k 1] a Cov [A(xk 1 xk 1)] Cov[wk 1] a T ACov (ek 1) A Cov (wk 1) Cuối cùng ta thu được f a T Pk A Pk 1 A Q (5.9)
  50. 50 5.2.3. Bước hiệu chỉnh (ước lượng hậu nghiệm): Sai số ước lượng hậu nghiệm a a ek  xk xk Đây là đại lượng ngẫu nhiên có giá trị trung bình thống kê bằng không và ma trận hiệp phương sai xác định như sau (chứng minh tương tự) a a a E[ek ] E[xk xk ] E[xk ] xk 0 a a a a T Pk Cov ek E ek (ek )  (5.10) Tìm ước lượng hậu nghiệm f Trước hết ta xét ý nghĩa của đại lượng H xk . Có thể coi đây là giá trị “đo đạc f dự báo”. Và khi ấy đại lượng zk H xk được gọi là measurement innovation (cập nhật đo đạc) or residual (độ dư đo đạc) (sai số giữa giá trị đo đạc thực tế và đo đạc dự báo). Giả thiết cơ bản của phương pháp lọc Kalman: ước luợng hậu nghiệm là tổ hợp tuyến tính của ước lượng tiên nghiệm và sai số giữa giá trị đo đạc thực tế và đo đạc dự báo đã được trọng số hóa bởi K như sau a f f xk xk K ( zk H xk ) (5.11) K được gọi là gain (blending factor). Đây là một ma trận (n x m) chiều Tìm ma trận hiệp phương sai của sai số ước lượng hậu nghiệm Ta có sai số ước lượng sau a a f f ek xk xk xk (xk K ( zk H xk )) f f xk xk K ( H xk  k H xk ) f (I K H )(xk xk ) K k Chú ý tới tính chất Cov( Kk ) Cov(Kk ) Lúc này ma trận hiệp phương sai có dạng:
  51. 51 a a f Pk Cov (xk xk ) Cov ((I K H ) (xk xk ) K k ) f Cov((I K H ) (xk xk )) Cov(K k ) f T T (I K H )Cov (xk xk ) (I K H ) K Cov( k ) K Cuối cùng, theo các định nghĩa, ta thu được a f T T Pk (I K H ) Pk (I K H ) K R K (5.12) 5.2.4. Tìm Kalman gain (blending factor) K: Mục tiêu xác định K K được tìm từ điều kiện tối thiểu của đại lượng sau a 2 E[ xk xk ] min Đây chính là điều kiện tối thiểu của kỳ vọng của bình phương sai số hậu nghiệm. Có thể chứng minh được rằng giá trị trên chính là trace (vết) của ma trận a hiệp phương sai sai số hậu nghiệm Pk . Thật vậy: a a a tr Pk  Pk [i,i]  Cov ( xk[i] xk [i]) i 1,n i 1,n a a 2  Var ( xk[i] xk [i])  E[ ( xk [i] xk [i]) ] i 1,n i 1,n Tóm lại, ta đã chứng minh được: a a 2 tr Pk E[( xk xk ) ] (5.13) a Biến đổi Pk Trước hết ta nhắc lại công thức ( I A B )T I ( A B )T I BT A T Với công thức trên ta có
  52. 52 a f T T Pk (I K H ) Pk (I K H ) K R K f f T T T (Pk K H Pk ) (I H K ) K R K f f f T T f T T Pk K H Pk Pk H K K (H Pk H R) K Cuối cùng ta được a f f f T T T Pk Pk K H Pk Pk H K K S K Với S có giá trị f T S (H Pk H R) Tìm đạo hàm của tr Pk Trước hết ta có các công thức sau:  tr ( A x B )  tr ( BT x T AT ) AT B T ;  x  x  tr ( A x B x T C ) AT C T x BT C A x B  x Ở đây đ ã sử dụng định nghĩa đạo hàm của một đại lượng vô hướng a theo ma trận x là một ma trận có các thành phần sau  a  a a, x a, x[i, j]  x  x[i, j] Với các công thức trên, lấy đạo hàm của vết của ma trận hiệp phương sai sai số ước lượng hậu nghiệm, ta có:  tr Pa  tr P f  tr (K H P f )  tr (P f H T KT )  tr ( K S KT ) k k k k  K  K  K  K  K Lấy đạo hàm từng thành phần ta được:  tr P f k 0  K  tr (K H P f )  tr ( I K H P f ) k k IT ( H P f )T (P f )T H T  K  K k k
  53. 53  tr (P f H T KT )  tr (P f HT KT I T ) k k I T P f H T P f H T  K  K k k tr ( K S K T ) tr ( I K S K T I ) K ST K S  K  K Tổng hợp lại ta được: a  tr Pk f T T f T T (Pk ) H Pk H K S K S  K f f T T T ( Pk (Pk ) ) H K ( S S ) f Chú ý tới tính đối xứng của các ma trận Pk và S, ta thu nhận được công thức cuối cùng: a  tr Pk a T 2 Pk H 2 K S  K Điều kiện cực tiểu của tr Pk:  tr Pa k 0 K P f H T S 1  K k Như vậy ta đã chứng minh và thu nhận được đầy đủ các bước của lọc Kalman f T f T 1 K Pk H ( H Pk H R ) P f H T k H P f H T R k Một số tính chất của K: 1 limR 0 K H lim f K 0 Pk 0
  54. 54 Hình 5.1: Tóm tắt sơ đồ quá trình lọc Kalman 5.3. Thuật toán lọc Kalman rời rạc: Nhắc lại mô hình quá trình và mô hình đo đạc x A x B u w k k 1 k k 1 z H x  k k k Thuật toán lọc Kalman gồm 2 bước, 5 phép tính: 5.3.1.Cập nhật theo thời gian – dự báo (ước lượng tiên nghiệm) (predict): 5.3.1.1. Phép tính s 1: Ở bước này từ mô hình quá trình, lấy trung bình thống kê, chú ý tới giá trị trung bình thống kê của wk bằng không, ta có: f a x A xk 1 Buk k 5.3.1.2. Phép tính số 2: Ma trận hiệp phương sai của sai số ước lượng tiền nghiệm ở bước một sẽ có dạng P f A Pa AT Q k k 1 Trong đó: a Pk 1– ma trận hiệp phương sai của sai số ước lượng hậu nghiệm tại k-1 Q – ma trận hiệp phương sai của quá trình ồn trắng có dạng Ew1w1 Ew1wl  Q  Ew w  ij i j E w w E w w  l 1  l l  Như vậy, A như một phép chiếu của trạng thái và ma trận hiệp phương sai sai số ước lượng hậu nghiệm tại thời điểm k-1 về thời điểm k.
  55. 55 5.3.2. Cập nhật theo đo đạc – chỉnh sửa (ước lượng hậu nghiệm) (correct): Tại bước này quá trình tính toán được thực hiện như sau: 5.3.2.1. Phép tính số 3: Tính Kalman gain bằng việc cập nhật đo đạc ( qua H và R ) K P f H T ( H P f H T R ) 1 k k k 5.3.2.2. Phép tính số 4: Cập nhật đo đạc tại thời điểm k để xác định trạng thái hậu nghiệm tại k xa x f K ( z H x f ) k k k k 5.3.2.3. Phép tính số 5: Xác định ma trận hiệp phương sai hậu nghiệm của các sai số ước lượng hậu nghiệm Pa ( I K H ) P f k k k Như vậy đây là quá trình truy hồi bắt đầu từ thời điểm ban đầu (k = 0). Nhưng tại thời điểm ban đầu ta không có các ước lượng trước đó, vì vậy ta phải cho trước một ước lượng ban đầu gồm: a a x0  x0 e0 a a a a T P0 Cov e0 E e0 (e0 )  2. Cập nhật theo quan trắc 2.1 Tính f T f T 1 1.Cập nhật theo thời gian Kk Pk H ( H Pk H R ) 1.1 Tính f a 2.2 Tính xk A xk 1 Buk a f f 1.2 Tính xk xk K ( zk H xk ) f a T P A P A Q k k 1 2.3 Tính a f Pk ( I Kk H ) Pk Ước lượng ban đầu a a x0 , P0 Hình 5.2: Sơ đồ tính lặp quá trình lọc Kalman
  56. 56 Chương 6. CÁC KẾT QUẢ ỨNG DỤNG MÔ HÌNH KẾT NỐI MARINE VÀ IMECH1D ĐỂ DỰ BÁO LƯU LƯỢNG VÀO HỒ HÒA BÌNH 6.1. Kết quả kiểm tra bài toán mẫu cho 10 lưu vực bộ phận: Mục đích bài toán kiểm định sự ổn định của bộ thông số mô hình. Tại thời điểm ban đầu, lưu vực hoàn toàn khô, cho mưa là một hằng số đều trên toàn lưu vực trong thời gian đủ dài để lưu lượng ra khỏi lưu vực là một hằng số. Dòng chảy trong lưu vực chính là dòng chảy dừng với biên vào là mưa, biên ra là lưu lượng, vận tốc, mực nước Khi dòng chảy đạt đỉnh và là hằng số theo thời gian, cho lượng mưa trong toàn lưu vực bằng không, tiếp tục tính một thời gian đủ dài để lưu lượng ra khỏi lưu vực đủ nhỏ để có thể bỏ qua. Kết thúc quá trình tính toán, kiểm tra hai thông số đó là: - Sự cân bằng nước, giữa tổng lượng nước vào lưu vực từ mưa và tổng lượng nước ra khỏi lưu vực theo kết quả tính toán. - Sự cân bằng lưu lượng lớn nhất, tức là lưu lượng của dòng chảy dừng - Chỉ tiêu kiểm tra: |(Vvào–Vra)/Vvào*100|<5% với Vvào=R*s*t, Vra=Qi; (R: Lượng mưa; s: Diện tích lưu vực; t: Thời gian tính) Gọi: R – Lượng mưa T – Tổng thời gian tính Tr – Thời gian mưa W – Sai số tổng lượng A – Diện tích của lưu vực Qmax – Sai số Q lớn nhất TT A R Tr T Vi Vo W Qi Qo Qmax LV (km2) (mm/h) (hours) (hours) (Trm3) (Tr m3) (%) (m3/s) (m3/s) (%) 1 3621 10 48 144 1738.08 1730.57 0.43229 10058.3 10032.5 0.25684 2 2441 10 48 144 1171.68 1162.4 0.79161 6780.56 6770.97 0.14137 3 2255 10 36 144 811.8 809.535 0.27897 6263.89 6254.7 0.1467 4 898.8 10 36 144 323.568 322.867 0.21673 2496.67 2493.91 0.11041 5 2211 10 54 144 1193.94 1182.05 0.99588 6141.67 6135.35 0.10285
  57. 57 6 3245 10 54 144 1752.3 1755.7 -0.194 9013.89 9057.11 -0.4795 7 3530 10 60 144 2118 2100.51 0.82561 9805.56 9784.86 0.21106 8 2141 10 60 144 1284.6 1281.07 0.27506 5947.22 5970.91 -0.3983 9 3153 10 48 144 1513.44 1562.01 -3.2094 8758.33 9068.43 -3.5406 10 2889 10 60 144 1733.4 1725.48 0.45667 8025 8005.35 0.24486 Bảng 6.1: Kết quả kiểm tra bài toán mẫu cho 10 lưu vực Các kết quả kiểm tra bài toán mẫu cho thấy cả 10 lưu vực bộ phận đều đã vượt qua được tiêu chí của bài toán mẫu đó là giá trị tuyệt đối của sai số về tổng thể tích nước và sai số của lưu lượng cực đại phải nhỏ hơn 5%. Bảng 6.1 cho thấy trong 10 lưu vực chỉ có lưu vực số 9 có sai số lớn hơn 3% con 9 lưu vực còn lại đều có sai số nhỏ hơn 1%. 6.2. Sử dụng mô hình kết nối Marine-IMech1D để dự báo lại trận lũ năm 2006 và hiệu chỉnh các tham số của mô hình: 6.2.1. Nhận định chung tình hình lũ sông Đà năm 2006: Năm 2006, trên sông Đà có tất cả 10 đợt lũ có Qmax đến hồ 3500 m3/s. Từ 1h/28/VI/2006 đến 19h/28/VII có 3 đợt lũ liên tiếp nhau: - Qmax 1 = 6000 m3/s vào lúc 17h/2/VII - Qmax 2 = 9000 m3/s vào lúc 23h/10/VII - Qmax3 = 14700 m3/s vào lúc 12h/9/VII Đây là lưu lượng đến hồ Hòa Bình lớn nhất năm và lớn thứ 7 trong chuỗi quan trắc từ năm 1902 đến nay (Qmax1 = 22500 m3/s năm 1996; Qmax2 = 17800 m3/s năm 1945; Qmax3 = 17200 m3/s năm 1964; Qmax4 = 16200 m3/s năm 1971; Qmax5 = 15800 m3/s năm 1969 ; Qmax6 = 15200 m3/s năm 2002). Đợt lũ từ 13h/6/X đến 7h/22/X với Qmax đến hồ Hòa Bình = 9200 m3/s 13h/12/X lớn thứ hai trong chuỗi số liệu lưu lượng lớn nhất đến hồ Hòa Bình trong tháng X từ năm 1902 đến nay (Qmax tháng X lớn nhất = 10900 m3/s năm 1932; Qmax tháng X lớn thứ ba = 9000 m3/s năm 1999). 6.2.2. Kết quả tính toán dự báo lại cho trận lũ năm 2006 bằng mô hình kết nối MARINE-IMECH1D: Với số liệu vào hàng ngày là số liệu mưa thực đo của 36 trạm trên lưu vực, số liệu lưu lượng và mực nước của một số trạm kiểm tra, kết hợp với số liệu mưa dự báo số trị. Vận hành bộ chương trình nối hợp MARINE-IMECH1D trên bộ máy tính Intel Pentium 4 3.0GHz mất khoảng 30 phút thu được các kết quả:
  58. 58 - Lưu lượng vào hồ dự báo 48 giờ. - Mực nước và lưu lượng dự báo ở các vị trí trên sông Đà. Dưới đây là các kết quả dự báo lại cho mùa lũ năm 2006 bằng bộ chương trình kết nôi MARINE-IMECH1D: Qmax thực Qmax tính toán (h) đo vào hồ Thời gian dự báo vào hồ Q % sai Số cửa Số cửa (T db-T Hòa Bình xuất hiện Hòa Bình sè xả đáy xả mặt (m3/s) tt) (m3/s) (m3/s) 10h/1/VII/ 5500 5000 -500 -9.1 14 0 0 2006 17h/2/VII/ 6000 5300 -700 -11.7 12 2006 23h/10/VII 9000 8900 -100 -1.1 2 1 0 /2006 12h/19/VII 14700 13808 -892 -6.1 1 5 0 /2006 23h/8/VIII 3600 3672 72 2 -4 0 0 /2006 17h/17/VII 3800 5457 1657 43.6 17 0 0 I/2006 5h/19/VIII 6000 4900 -1100 -18.3 6 0 0 /2006 1h/25/VIII 4000 3238 -762 -19.1 0 0 0 /2006 13h/12/X/ 9200 8900 -300 -3.3 0 4 0 2006 Bảng 6.2: Bảng kết quả dự báo đỉnh lũ đến hồ Hòa Bình của toàn mùa lũ năm 2006
  59. 59 Hình 6.1: Kết quả dự báo lưu lượng vào hồ HB từ 05 đến 30 tháng 7 năm 2006 không sử dụng mô đun lọc Kalman 16000 Thuc DoDo 15000 HieuChinh 14000 13000 12000 11000 10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 0 200 220 240 260 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 Hình 6.2: Kết quả dự báo lưu lượng vào hồ HB tháng 7 năm 2006 có sử dụng mô đun lọc Kalman So sánh kết quả tính toán trên hình 6.1 và hình 6.2 có thể thấy mô đun lọc Kalman đã cho kết quả tính có độ chính xác cao hơn so với kết quả tính nguyên bản của mô hình kết nối.
  60. 60 Ốp thời gian dự kiến Mức sai số Số điểm đúng/Tổng số P% 6 h 18.5 %Q 97/116 82.5 12h 18.5 %Q 103/116 83.6 18h 18.5 %Q 94/116 81 24h 18.5 %Q 89/116 76.7 30h 21.5%Q 88/115 76.5 36h 24.5%Q 81/115 70.4 42h 27.5%Q 83/115 72.1 48h 30.5%Q 83/115 72.1 Bảng 6.3: Kết quả dự báo quá trình lũ đến hồ Hòa Bình (Từ 18/VI đến 23/X/ 2006) Hình 6.3: Đường quá trình dự báo của bản tin dự báo ngày 18/07/2006 6.3. Kết quả sử dụng mô hình kết nối Marine-IMech1D tác nghiệp cho mùa lũ năm 2009: Trong mùa lũ năm 2009, mô hình kết nối Marine và IMech1D đã được sử dụng trong công tác dự báo lũ trên sông Đà tại Viện Cơ học và tại Trung tâm Dự báo Khí tượng Thủy văn Trung ương. Các kết quả tác nghiệp của mô hình là rất
  61. 61 khả quan. Dưới đây là một số kết quả dự báo của mô hình cho trận lũ diễn ra từ ngày 04 tháng 07 đến ngày 10 thang 07 năm 2009. Đây cũng là trận lũ lớn nhất của năm 2009 trên sông Đà. 14000 Tính toá n HieuChinh 13000 ThucDo 12000 11000 10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 0 400 420 440 460 480 500 520 540 560 580 600 620 Hình 6.4: Đường quá trình dự báo lưu lượng vào hồ Hòa Bình ngày 05/07/2009 14000 Tính toá n Hie uChinh 13000 ThucDo 12000 11000 10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 0 400 450 500 550 600 650 700 750 800 Hình 6.5: Đường quá trình dự báo lưu lượng vào hồ Hòa Bình ngày 12/07/2009
  62. 62 9000 TínhHieuCHinh toá n ThucDo 8000 7000 6000 5000 4000 3000 2000 1000 0 610 630 650 670 690 710 730 750 770 790 810 830 Hình 6.6: Đường quá trình dự báo lưu lượng vào hồ Hòa Bình ngày 13/07/2009 9000 HieuCHinh ThucDo 8000 7000 6000 5000 4000 3000 2000 1000 0 610 630 650 670 690 710 730 750 770 790 810 830 Hình 6.7: Đường quá trình dự báo lưu lượng vào hồ Hòa Bình ngày 14/07/2009
  63. 63 KẾT LUẬN Luận văn Thạc sỹ Cơ học này đã giải quyết được bài toán vô cùng khó khăn đó là “xây dựng mô hình phục vụ công tác dự báo lưu lượng vào hồ Hòa Bình” bằng cách sử dụng các phần mềm được xây dựng trên cơ sở các thuật toán tiên tiến nhất mà trên thế giới hiện đang nghiên cứu và phát triển. Qua luận văn này, tác giả đã có được những hiểu biết sâu sắc về nhiều mặt: - Hiểu sâu và rộng về học thuật của các mô hình thủy văn, thủy lực. - Hiểu biết sâu hơn về các phương pháp mô phỏng số, đặc biệt là mô phỏng số trong cơ học chất lỏng. - Có được kiến thức và khả năng sử dụng tốt các công cụ lập trình, các phần mềm GIS. - Đặc biệt qua luận văn này tác giả đã có thêm được nhiều kiến thức thực tế và nhiều kinh nghiệm quý báu trong việc xây dựng, phát triển mô hình thủy văn, thủy lực cũng như triển khai sử dụng các mô hình này cho các bài toán thực tế. Một số kết quả của luận văn đã được công bố trên các tạp chí, hội nghị khoa học cấp quốc gia và quốc tế. Các kết quả này cũng đã được nghiệm thu qua các hội đồng khoa học cấp bộ và cấp nhà nước (của các đề tài cấp bộ, cấp nhà nước mà tác giả đã tham gia). Sản phẩm của luận văn, “Mô hình kết nối thủy văn Marine và thủy lực IMech1D để dự báo lưu lượng vào hồ Hòa Bình” đã được chuyển giao cho Trung tâm Dự báo Khí tượng Thủy văn trung ương và đã được đưa vào khai thác với vai trò là một công cụ chính thức để dự báo lượng lưu lượng vào hồ Hòa Bình trong mùa lũ trước 48 giờ. Luận văn hoàn thành xuất sắc các nội dung đã đặt ra trong bản đề cương. Sản phẩm của luận văn đã trang bị được một công cụ hữu ích cho những người làm công tác phòng chống lũ lụt trên hệ thống sông Hồng sông Thái Bình. Để công cụ này có thể đạt được nững kết quả tốt hơn cũng như có thể sử dụng được cho các lưu vực sông khác cần phải có thêm nhiều sự đầu tư hơn nữa về công sức và trí tuệ.
  64. 64 DANH MỤC CÔNG TRÌNH CỦA TÁC GIẢ 1. Hoang Van Lai, Nguyen Van Diep, Nguyen Tien Cuong, Nguyen Hong Phong et al., “Coupling hydrological-hydraulic models for extreme floods simulating and forecasting in the Central Vietnam.” River Basin Management 2009, The fifth International -Conference, 7-9 September 2009, Malta. Pp. 113-123 2. Nguyen Tien Cuong, Hoang Van Lai, Trinh Thu Phuong, “Model for Forecasting and warning model into Hoa Binh Reservoir”, Proceedings of the National congress on Fluid Mechanics, 2009. 3. Nguyen Tien Cuong, Trinh Thu Phuong, “Forecasting the discharge into Hoa Binh Reservoir by applying the Connecting model MARIN- IMech1D” Vietnam Journal of Mechanics, VAST, Vol.30, No.3 (2008), pp. 149-157 4. Nguyen Tien Cuong, Tran Thu Ha, “On the model forecasting floods in Ferfume river”. Proceedings of the 8th National congress on Mechanics 2007. pp.51-64. 5. Nguyen Tien Cuong, Nguyen Hong Phong, “Application of hydrological and hydraulic model to forecast flood on Da’s river in PC system”, Proceedings of the National congress on Fluid Mechanics, 2007, pp.79-88. 6. Nguyen Tien Cuong, “Simulation flood in Perfume river system basin, Hue city (Vietnam) by MARINE hydrological model”, Proceedings of the National congress on Fluid Mechanics, 2005, pp. 9-19. 7. Nguyen Tien Cuong, “Sensitivity estimation of Parameters in Distributed Hydrological Model”, Proceedings of the National congress on Fluid Mechanics, Phan Thiet city, 2008, pp. 21-29 8. Hoang Van Lai, Nguyen Tien Cuong, Nguyen Hong Phong, “Development and Application forecasting and warning software for flood of Perfume river basin system in the middle of Vietnam”, Proceedings of application for Informatics Technology and Communication, Thai Nguyen city, 2007. 9. Nguyen Tien Cuong, Nguyen Thanh Don, “Application of the connected Marine and Telemac 2D to simulation flash flood in SonLa town (in
  65. 65 Vietnam)”, Proceedings of the National congress on Mechanics 2004, V. 2, pp.8-15. 10. Nguyen Tien Cuong, Marie Madeleine Maubourguet, “The first time application MARINE hydrological model for Da river basin in Vietnam domain”, Proceedings of the National congress on Mechanics 2004, V. 2, pp.16-25. 11. Nguyen Tien Cuong, Nguyen Thanh Don, “Application of the connected Marine and Telemac 2D to simulation flash flood in DongChuan city (in China)”, Proceedings of the National congress on Fluid Mechanics 2004, pp.23-38.
  66. 66 TÀI LIỆU THAM KHẢO Tiếng Việt 1. Phạm Kỳ Anh, Giải tích số, Nhà xuất bản Đại học Quốc gia, Hà Nội. 2. Đỗ Cao Đàm, Trịnh Quang Hòa, Hà Văn Khôi, Đoàn Trung Lưu, Nguyễn Năng Minh, Lê Đình Thành, Dương văn Tiến, (1993), Thủy Văn công trình, Nhà xuất bản Nông nghiệp. 3. Trung tâm QG Dự báo KTTV, Phòng Dự báo Thủy văn (2005), Tổng kết công tác dự báo thủy văn sông Đà phục vụ công trình Hoà Bình năm 1988- 2005. 4. Trịnh Quang Hoà và nnk (7 - 1997), Nghiên cứu xây dựng công nghệ nhận dạng lũ thượng lưu sông Hồng phục vụ việc điều hành Hồ chứa Hoà Bình chống lũ hạ du. Đề tài nghiên cứu khoa học cấp nhà nước, Bộ NN và PTNT, Trường Đại học Thuỷ Lợi, Hà Nội. 5. Báo cáo của dự án quốc tế FLOCODS (2004) phần Các mô hình thủy văn và thủy lực. 6. Báo cáo của đề tài cấp nhà nước (2004), mã số KC.08-13 phần mô hình thủy văn và thủy lực một chiều mở rộng. Tiếng Anh 7. Ven Techow, David R.Maidment, Larry W.Mays. (1968), Applied Hydrology, Mc Graw - Hill, New York. 8. Cung J.A., Holly F.M. Verwey A. (1980), Practical Aspects of Computational River Hydraulics, Pitman Advanced Publishing Program. 9. C.A.J. Fletcher. Computational Techniques for Fluid Dynamics (1987), Spring-Venlag, Vol. 1, pp.183. 10. Bahram Saghafian, Pierre Y. Julien. H. Rajaie (2002), Runoff hydrograph simulation based on time variable isochrone technique, Journal of Hydrology Vol.26, pp.193-203. 11. David R. Maidment, Handbook of Hydrology, Mc Graw-Hill Bookcompany.
  67. 67 12. H. Moradkhani, S. Sorooshian (2006), General Review of Rainfall-Runoff Modelling: Model calibration, Data assimilation and uncertainty analysis. 13. Henrik Madsen, Claus Skotner (2005), Adaptive state updating in real-time river flow forecasting – a combined filtering and error forecasting procedure, J. of Hydrology, Vol.308, pp.302-312. 14. Brankin, R.W., I. Gladwell, and L.F. Shampine, RKSUITE (1991), A Suite of Runge-Kutta Codes for the Initial Value Problem for ODEs, Softreport 91-1, Mathematics Department, Southern Methodist University, Dallas, Texas. Tiếng Pháp 15. Alquier, M., Chorda, J., Dartus, D., Estupina, V., Guennec, B. L. & Maubourguet, M.-M. (2000), Contribution des données d’observation de la terre à la cartographie de l’aléa, Rapport d’étude adema. Institut de Mécanique des Fluides de Toulouse. 16. Roux, H. (2004), Estimation de paramètres en hydraulique fluviale, à partir de donnée caractéristiques de l'imagerie aérienne, Thèse à l'IMFT/HYDRE - 1 Allée du Professeur Camille Soula, 31400 Toulouse, 65-75, 120-121. 17. Kham, L. X., Dartus, D., Maubourguet, M.-M., and Jacques, C. (2006), Sensitivity analysis for Manning coefficient on the Gardons d'Anduze Basin, France, Vietnam - Japan Estuary Workshop, 66-71.
  68. 68 PHỤ LỤC 1. Kết quả kiểm tra bài toán mẫu cho 10 lưu vực bộ phận của lưu vực sông Đà bằng MARINE: Mục đích bài toán kiểm định sự ổn định của bộ thông số mô hình. Tại thời điểm ban đầu, lưu vực hoàn toàn khô, cho mưa là một hằng số đều trên toàn lưu vực trong thời gian đủ dài để lưu lượng ra khỏi lưu vực là một hằng số. Dòng chảy trong lưu vực chính là dòng chảy dừng với biên vào là mưa, biên ra là lưu lượng, vận tốc, mực nước Khi dòng chảy đạt đỉnh và là hằng số theo thời gian, cho lượng mưa trong toàn lưu vực bằng không, tiếp tục tính một thời gian đủ dài để lưu lượng ra khỏi lưu vực đủ nhỏ để có thể bỏ qua. Kết thúc quá trình tính toán, kiểm tra hai thông số đó là: - Sự cân bằng nước, giữa tổng lượng nước vào lưu vực từ mưa và tổng lượng nước ra khỏi lưu vực theo kết quả tính toán. - Sự cân bằng lưu lượng lớn nhất, tức là lưu lượng của dòng chảy dừng. - Chỉ tiêu kiểm tra: |(Vvào–Vra)/Vvào*100|<5% với Vvào=R*s*t, Vra=Qi; (R: Lượng mưa; s: Diện tích lưu vực; t: Thời gian tính) Gọi: R – Lượng mưa T – Tổng thời gian tính Tr – Thời gian mưa W – Sai số tổng lượng A – Diện tích của lưu vực Qmax – Sai số Q lớn nhất Bảng 1: Kết quả bài toán kiểm tra các tiểu lưu vực TT A R Tr T Vi Vo W Qi Qo Qmax LV (km2) (mm/h) (hours) (hours) (Trm3) (Tr m3) (%) (m3/s) (m3/s) (%) 1 3621 10 48 144 1738.08 1730.57 0.43229 10058.3 10032.5 0.25684 2 2441 10 48 144 1171.68 1162.4 0.79161 6780.56 6770.97 0.14137 3 2255 10 36 144 811.8 809.535 0.27897 6263.89 6254.7 0.1467 4 898.8 10 36 144 323.568 322.867 0.21673 2496.67 2493.91 0.11041 5 2211 10 54 144 1193.94 1182.05 0.99588 6141.67 6135.35 0.10285 6 3245 10 54 144 1752.3 1755.7 -0.194 9013.89 9057.11 -0.4795 7 3530 10 60 144 2118 2100.51 0.82561 9805.56 9784.86 0.21106 8 2141 10 60 144 1284.6 1281.07 0.27506 5947.22 5970.91 -0.3983 9 3153 10 48 144 1513.44 1562.01 -3.2094 8758.33 9068.43 -3.5406
  69. 69 10 2889 10 60 144 1733.4 1725.48 0.45667 8025 8005.35 0.24486 Dưới đây là đồ thị kết quả kiểm tra bài toán mẫu cho 10 lưu vực bộ phận: kÕt qu¶ bµi to¸n test tiÓu l­u vùc 2 8000 0 10 6000 M­a (mm) 20 30 4000 L­u l­îng (m3/s) 40 50 2000 60 L­u l­îng(m3/s) 70 0 80 0 10 20 30 40 50 60 70 80 90 10 11 12 13 14 0 0 0 0 0 Thêi gian (h)
  70. 70 kÕt qu¶ bµi to¸n test tiÓu l­u vùc 3 8000 0 6000 M­a (mm) 20 4000 L­u l­îng (m3/s) 40 2000 60 L­u l­îng(m3/s) 0 80 0 10 20 30 40 50 60 70 80 90 10 11 12 13 14 0 0 0 0 0 Thêi gian (h) kÕt qu¶ bµi to¸n test tiÓu l­u vùc 4 4000 0 10 3000 20 M­a (mm) 30 2000 40 L­u l­îng (m3/s) 50 1000 60 L­u l­îng(m3/s) 70 0 80 0 10 20 30 40 50 60 70 80 90 10 11 12 13 14 0 0 0 0 0 Thêi gian (h)
  71. 71 kÕt qu¶ bµi to¸n test tiÓu l­u vùc 5 8000 0 10 6000 M­a (mm) 20 L­u l­îng (m3/s) 30 4000 40 50 2000 60 L­u l­îng(m3/s) 70 0 80 0 10 20 30 40 50 60 70 80 90 10 11 12 13 14 0 0 0 0 0 Thêi gian (h) kÕt qu¶ bµi to¸n test tiÓu l­u vùc 6 12000 0 10000 10 M­a (mm) 20 8000 30 6000 L­u l­îng (m3/s) 40 4000 50 60 L­u l­îng(m3/s) 2000 70 0 80 0 10 20 30 40 50 60 70 80 90 10 11 12 13 14 0 0 0 0 0 Thêi gian (h)
  72. 72 kÕt qu¶ bµi to¸n test tiÓu l­u vùc 7 12000 0 10000 10 M­a (mm) 20 8000 L­u l­îng (m3/s) 30 6000 40 4000 50 60 L­u l­îng(m3/s) 2000 70 0 80 0 10 20 30 40 50 60 70 80 90 10 11 12 13 14 0 0 0 0 0 Thêi gian (h) kÕt qu¶ bµi to¸n test tiÓu l­u vùc 8 8000 0 10 6000 M­a (mm) 20 30 L­u l­îng (m3/s) 4000 40 50 2000 60 L­u l­îng(m3/s) 70 0 80 0 10 20 30 40 50 60 70 80 90 10 11 12 13 14 0 0 0 0 0 Thêi gian (h)
  73. 73 kÕt qu¶ bµi to¸n test tiÓu l­u vùc 9 12000 0 10000 10 M­a (mm) 20 8000 30 6000 L­u l­îng (m3/s) 40 4000 50 60 2000 L­u l­îng(m3/s) 70 0 80 0 10 20 30 40 50 60 70 80 90 10 11 12 13 14 0 0 0 0 0 Thêi gian (h) kÕt qu¶ bµi to¸n test tiÓu l­u vùc 10 12000 0 10000 10 M­a (mm) 20 8000 L­u l­îng (m3/s) 30 6000 40 4000 50 60 L­u l­îng(m3/s) 2000 70 0 80 0 10 20 30 40 50 60 70 80 90 10 11 12 13 14 0 0 0 0 0 Thêi gian (h) Các kết quả kiểm định cho thấy cả 10 lưu vực bộ phận đều đạt tiêu chí của bài toán mẫu đó là giá trị tuyệt đối của phần trăm sai số về tổng lượng và của đỉnh lũ không được vượt quá 5%.
  74. 74 2. Kết quả kiểm định bộ chương trình tính toán thủy lực một chiều IMech1D bằng các bài toán kiểm định mẫu (Test Cases). Để thẩm định các mô hình tính toán thuỷ lực một chiều, các phòng thí nghiệm thuỷ lực của Châu Âu đã thống nhất đưa ra 12 bài toán mẫu (test cases). Các bài toán mẫu này được phân thành 3 nhóm như sau: Nhóm 1: Các bài toán kiểm định mẫu có nghiệm giải tích Bài toán mẫu số 1: Sóng xả trong kênh chữ nhật Bài toán mẫu số 2: Dòng chảy êm, đều trong kênh hình chữ nhật Bài toán mẫu số 3: Dòng chảy đều khi có lưu lượng phụ Bài toán mẫu số 4: Dòng chảy đều khi có công trình. Nhóm 2: Các bài toán mẫu kiểm định độ nhậy của sơ đồ Bài toán mẫu số 5: Sóng động lực học Bài toán mẫu số 6: Sóng khuyếch tán Bài toán mẫu số 7: Sóng động học Bài toán mẫu số 8: Sóng lũ qua hồ chứa Nhóm 3: Các bài toán mẫu kiểm định các thành phần hệ thống của sơ đồ Bài toán mẫu số 9: Nhiễu cục bộ đối với dòng chảy đều Bài toán mẫu số 10: Hình học không đồng nhất trong dòng chảy đều Bài toán mẫu số 11: Dòng chảy không đều trong lòng dẫn phức tạp Bài toán mẫu số 12: Phân lưu Bộ chương trình tính toán thủy lực một chiều IMech1D đã được kiểm định theo các bài toán kiểm định mẫu kể trên. Kết quả kiểm định chỉ ra rằng bộ chương trình tính toán thủy lực 1 chiều IMech1D đã vượt qua tất cả các bài toán mẫu kiểm định do các Phòng thí nghiệm thủy lực của Châu Âu đề nghị. 2.1. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 1: SÓNG XẢ TRONG KÊNH CHỮ NHẬT Mục đích của bài toán kiểm định mẫu: - Thử các thành phần lan truyền trong hệ phương trình Saint - Venant. Mô tả bài toán kiểm định mẫu: Xét dòng chảy trong kênh hình chữ nhật với các đặc trưng sau: - Chiều rộng kênh B = 200 mét - Không có độ dốc, (I = 0)
  75. 75 - Không có ma sát, (k=9999; k-hệ số Strickler) - Độ dài của kênh L = 18000 mét. - Dòng chảy tại thời điểm ban đầu thỏa mãn các điều kiện sau: + Dòng chảy đều đồng nhất + Vận tốc U0 = 2 m/s + Cao trình mặt nước không đổi ở mọi nơi; H0 = 5m - Tại các biên trên (thượng lưu), biên dưới (hạ lưu) dòng chảy thỏa mãn các điều kiện: Tại thượng lưu: quy luật giảm tuyến tính lưu tốc từ Uo=2m/gy xuống 0m/gy sau một thời gian Tmax=1000sec U(x=0, t) = Uo(1 - t/Tmax) Tại hạ lưu: Cao trình mực nước không đổi; H = Ho = 5m. Bài toán này có nghiệm giải tích, sau này ta gọi là nghiệm chính xác. Chúng ta sẽ tìm nghiệm số của bài toán bằng bộ chương trình tính toán thủy lực một chiều IMech1D. So sánh nghiệm giải số và nghiệm chính xác của bài toán ta thấy rằng bộ chương trình tính toán thủy lực IMech1D vượt qua bài toán kiểm định mẫu số. 2.2. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 2: DÒNG CHẢY ÊM, ĐỀU TRONG KÊNH HÌNH CHỮ NHẬT Mục đích của bài toán kiểm định mẫu: - Thử các thành phần gradient mực nước, ma sát trong phương trình Saint- Venant ở trạng thái dừng Mô tả bài toán kiểm định mẫu: Xét dòng chảy trong kênh hình chữ nhật với các đặc trưng sau: - Chiều rộng kênh B = 100 mét - Không có độ dốc I = 0.0005 - Ma sát không đổi theo lòng sông: k=30.59 (k - hệ số Strickler) - Độ dài kênh L = 10000 mét. - Dòng chảy tại thời điểm ban đầu thỏa mãn các điều kiện sau: + Dòng chảy đều đồng nhất + Lưu lượng không đổi; Q0 = 1000 m3/s + Độ sâu cột nước không đổi: H0 = 3 m
  76. 76 - Tại các biên trên (thượng lưu) biên dưới (hạ lưu) dòng chảy thỏa mãn các điều kiện: Tại thượng lưu: (x = 0): cho lưu lượng: Q = 1.000 m3/sec. Tại hạ lưu: cho mực nước: Trường hợp 1: H = 5 m Trường hợp 2: H = 7 m Trường hợp 3: H = 3 m Có thể thu nhận nghiệm giải tích số của bài toán từ việc giải số phương trình vi phân thường mô phỏng dòng chảy dừng bằng phương pháp Runge- Kutta: q 2 q 2 I I dH 2 2 4 / 3 2 10 / 3 Q k R H k H ; q dx q 2 q 2 B 1 1 gH 3 gH 3 Chúng ta sẽ tìm nghiệm số của bài toán bằng bộ chương trình tính toán thủy lực một chiều IMech1D. So sánh nghiệm giải tích số và nghiệm số của bài toán ta thấy rằng bộ chương trình tính toán thủy lực IMech1D vượt qua bài toán mẫu điểm định số 2. 2.3. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 3: DÒNG CHẢY ĐỀU CÓ LƯU LƯỢNG PHỤ Mục đích của bài toán kiểm định mẫu: - Thử các thành phần lưu lượng bổ xung trong phương trình Saint-Venant. Mô tả bài toán kiểm định mẫu: Xét dòng chảy trong kênh có lưu lượng phụ đổ vào - Dòng chính chảy trong kênh hình chữ nhật, độ dốc không đổi I=0.0005, chiều dài 5.000 m - Cứ 1.000 m lại có một lưu lượng phụ bổ xung trên một đoạn 10m. - Nhánh sông mang các lưu lượng phụ có độ dài 1.000 m và hoàn toàn đồng dạng với nhánh chính (trừ chiều rộng). - Các giá trị khác như lưu lượng, chiều rộng của tất cả các nhánh (chính và phụ) được lựa chọn sao cho cột nước là đồng nhất trên mọi nhánh và mọi dH Q 2 điểm, nghĩa là sao cho 0 U 2 I k 2 H 10 / 3 const dx B 2 - Tại đoạn nhập lưu: Không có độ dốc, I=0
  77. 77 Không có ma sát, k=9999 Không có tổn thất áp. - Tại các biên trên (thượng lưu) biên dưới (hạ lưu) dòng chảy thỏa mãn các điều kiện: Tại thượng lưu: Cho lưu lượng vào nhánh chính và nhánh phụ. Tại hạ lưu: Cho chiều cao cột nước tương ứng với trường hợp dòng chảy đều. Đã kiểm định IMech1D cho được trường hợp hệ thống trên là một mạng sông gồm một nhánh chính và 4 nhánh phụ. Cao độ đáy thượng lưu là 10mét. Cao độ đáy hạ lưu là 7.5 mét. Biên lưu lượng ở thượng lưu nhánh chính là 600 m3/giây. Biên lưu lượng ở thượng lưu các nhánh phụ là 100 m3/giây. Biên mực nước ở hạ du là 12.5 mét (độ sâu cột nước ở đây là 5 mét). Theo các số liệu trong phương án kiểm tra thì bài toán sẽ có nghiệm dừng đều với chiều sâu cột nước trên toàn hệ thống là 5 mét. Kết quả tính toán thu được ổn định và gần như trùng với nghiệm chính xác. 2.4. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 4: DÒNG CHẢY ĐỀU CÓ CÔNG TRÌNH Mục đích: - Thử việc đưa các công trình vào hệ phương trình Saint - Venant. Mô tả bài toán kiểm định mẫu: Xét dòng chảy trong kênh với các đặc trưng sau: - Kênh hình chữ nhật B=100 m, độ dốc đều I=0.0005, độ dài 5.000m - Ma sát đều k=30.59. - Công trình đặt ở điểm có toạ độ X = 4.000 m Điều kiện biên: - Thượng lưu: cho lưu lượng - Hạ lưu: cho chiều cao cột nước như trong trường hợp dòng chảy đều. - Công trình: cho quan hệ giữa lưu lượng, mực nước trước và sau công trình. Đã kiểm định IMech1D cho trường hợp lưu lượng chảy qua công trình theo 3/2 quy luật Q = a(H - Hs) với Hs = 3 m; a = 100 Kết quả tính toán bằng bộ chương trình IMech1D ổn định. Kết quả tính toán của bộ chương trình IMech1D phù hợp với nghiệm giải tích.
  78. 78 2.5. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 5: SÓNG ĐỘNG LỰC HỌC Mục đích: - Thử sơ đồ. - Thử độ nhậy của bước thời gian đối với sự lan truyền sóng động lực học. Mô tả bài toán mẫu: Sóng động lực học là loại sóng nhanh. Các thành phần gia tốc của khối nước đóng vai trò chủ chốt trong quá trình lan truyền sóng Xét dòng chảy trong kênh với cá đặc trưng sau: - Kênh hình chữ nhật nằm ngang chiều rộng B=200m, chiều dài L=16000m, - Không có ma sát (k=9999.). - Điều kiện biên: - Thượng lưu: cho lưu lượng thay đổi theo quy luật như hình vẽ sau Q(m3/g y) Qmax Q0 0 T 2T t Q0 t 0 t Q Q Q 0 t T 0 max 0 Q(0,t ) T t Qmax Qmax Q0 1 T t 2T T Q0 t 2T (trong đó T rất nhỏ) - Hạ lưu: cho chiều cao không đổi. - Điều kiện ban đầu: lưu lượng và cột nước không đổi. 3 Đã kiểm định bộ chương trình IMech1D cho trường hợp Q0=1.000m /s, 3 Qmax=2.000 m /s , Hhadu = 2,5 m T = 30s
  79. 79 Kết quả tính toán cao trình mực nước và lưu lượng của dòng chảy tại các vị trị x=0m, 4000m, 8000m và 12000m phù hợp với kết quả của các phòng thí nghiệm thủy lực Châu Âu công bố. 2.6. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 6: SÓNG KHUẾCH TÁN Mục đích: - Thử sơ đồ. - Thử độ nhậy của bước tính. Mô tả bài toán mẫu: Sóng khuyếch tán là loại sóng chậm (như một sóng lũ tự nhiên), trong đó thành phần gia tốc nhỏ so với lực trọng trường, áp suất và ma sát. Xét dòng chảy trong kênh với các đặc trưng sau: - Kênh hình chữ nhật có chiều rộng B=200 m, độ dài L=50.000m, độ dốc I=0,0005, cao độ đáy thượng du Zđ = 100 m. - Hệ số Strickler k = 35. - Điều kiện biên thượng lưu: Q0 t 0 t Q Q Q 0 t T 0 max 0 Lưu lượng Q(0,t ) T t Qmax Qmax Q0 1 T t 2T T Q0 t 2T Q(m3/g y) Qmax Q0 0 T 2T t - Điều kiện biên hạ lưu: quan hệ cột nước - lưu lượng f(HL, QL) - Điều kiện ban đầu: Q0, H0. - Thời gian tính toán cực đại Tmax. Đã kiểm định bộ chương trình IMech1D cho trường hợp 3 3 Q0=750m /s, H0=2,56 m, Qmax= 2.000 m /s.
  80. 80 T = 1.800 s, Tmax = 21.600 s, H(L, t) = H0 = 2,56 m Kết quả tính toán quá trình lưu lượng của dòng chảy bằng bộ chương trình IMech1D tại các điểm x=0, x=12500, x=25000 và x=37500m phù hợp với kết quả trên đồ thị của các phòng thí nghiệm thuỷ lực của Châu Âu đã công bố khi kiểm tra bài toán mẫu số 6. 2.7. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 7: SÓNG ĐỘNG HỌC Mục đích: - Thử sơ đồ - Thử sự lan truyền sóng động học (độ dốc lớn). Mô tả: Sóng động học là sóng xẩy ra trong kênh có độ dốc đủ lớn. Trong trường hợp này thành phần lực trọng trường (gI) cùng với thành phần ma sát đóng vai trò quyết định. Xét dòng chảy trong kênh có các đặc trưng sau: - Kênh hình chữ nhật có chiều rộng B = 100 m, chiều dài L = 100.000m, độ dốc I = 0,005, độc cao đáy thượng du Zđ = 500 m. - Hệ số Strickler k=35. - Điều kiện biên thượng lưu: Cho lưu lượng Q(0, t) biến thiên theo quy luật sau: Q0 t 0 t Q Q Q 0 t T 0 max 0 Q(0,t) T t Qmac Qmax Q0 1 T t 2T T Q0 t 2T Q(m3/g y) Qmax Q0 0 T 2T t - Điều kiện biên hạ lưu: cho quan hệ f(QL, HL).
  81. 81 - Điều kiện ban đầu: Q(x, 0) = Q0 H(x, 0) = H0 3 Đã kiểm định bộ chương trình IMech1D cho trường hợp Q0=1.000m /s; 3 Qmax=2.000 m /s. H0=2,49 m. Thời gian tính toán lớn nhất Tmax = 129.600 s, T = 12 giờ, HL = H0 Kết quả tính toán bằng bộ chương trình IMech1D đường mực nước tại các điểm không gian x=0, x=25000, x=50000 và x=75000mét theo toàn bộ thời gian tính toán phù hợp với kết quả của các phòng thí nghiệm thủy lực Châu Âu đã công bố khi kiểm tra bài toán mẫu số 7. 2.8. Kết quả kiểm định bộ chương trình IMech1D bằng bài toán mẫu số 8: SÓNG LŨ QUA HỒ CHỨA Mục đích: - Thử sơ đồ - Thử kiểm tra mô hình hồ chứa. Mô tả: Xét dòng chảy trong hồ có các đặc trưng sau: - Hồ chứa hình chữ nhật chiều dài L = 4000 m, chiều rộng B = 1000 mét , chiều cao đáy hạ du Zf =25 m. - Đáy phẳng và không có ma sát. - Điều kiện biên: + Thượng du: cho lưu lượng Q = Q(0, t). + Hạ du: cho quan hệ đập tràn. - Điều kiện ban đầu: Dòng chảy đều Q = Q0 H = H0. Đã kiểm định bộ chương trình tính toán IMech1D cho trường hợp: 3/ 2 2 2 2 Biên thượng lưu: Q 0,t a H 0 H1 1 cos t BL H1 sin t T T T ở đây: a là hằng số cho trong quan hệ lưu lượng và mực nước tại hạ lưu. 3/ 2 Q L,t a Z Z s H0, H1 là các hằng số trong biểu thức biến đổi mặt thoáng: