Luận văn Thạc sĩ: 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

ĐẠI HỌC QUỐC GIA HÀ NỘI

VIỆN KHOA HỌC VÀ CÔNG NGHỆ VIỆT NAM

TRƢỜNG ĐẠI HỌC CÔNG NGHỆ

VIỆN CƠ HỌC

NGUYỄN TIẾN CƢỜNG

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

LUẬN VĂN THẠC SỸ

Hà Nội - 2010

MỤC LỤC

MỤC LỤC ......................................................................................................... 1

MỞ ĐẦU ........................................................................................................... 7

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 ....................................................................................... 9

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]: ..................................................... 9

1.1.1. Đặc điểm mƣa gây lũ [4]: ..................................................... 9

1.1.2 Đặc điểm dòng chảy lũ sông Đà [4]: .................................... 10

1.1.3. Vai trò của hồ Hòa Bình [4]:............................................... 11

1.2. Tổng quan về mô hình thủy văn [3]-[5]: ....................................... 12

1.3. Tổng quan về mô hình thủy lực [3]; [6]: ....................................... 13

Chƣơng 2. PHẦN MỀM THỦY VĂN THAM SỐ PHÂN BỐ MARINE ......... 16

2.1. Cơ sở khoa học của phần mềm thủy văn tham số phân bố Marine: .. 16

2.1.1. Mô hình dòng chảy trên bề mặt lƣu vực [3]; [5]; [6]; [8]: ... 16

2.1.2. Mô hình thấm Green Ampt [12]; [16]: ................................ 17

2.2. Cấu trúc dữ liệu trong Marine [16]: ................................................. 19

Chƣơng 3. PHẦN MỀM THỦY LỰC MỘT CHIỀU IMECH1D ..................... 22

3.1. Các thành phần của hệ thống [5]; [6]; [9]: ....................................... 22

3.1.1. Mạng sông: ......................................................................... 22

3.1.1.1. Nút sông: ............................................................... 22

3.1.1.2. Đoạn sông: ............................................................ 23

3.1.2. Ô ruộng (Ô chứa): .............................................................. 23

3.2. Mô hình toán học [5]; [6]; [9]: ......................................................... 23

3.2.1. Mô hình toán học một đoạn sông: ....................................... 23

3.2.2. Mô hình toán học của một ô ruộng [5]; [6]; [9]: .................. 24

3.3. Lƣợc đồ sai phân [1]; [5]; [6]: ......................................................... 24

3.4. Tuyến tính hóa hệ phƣơng trình (3.5), (3.7), (3.8): ......................... 26

3.4.1. Tuyến tính hoá các biểu thức đơn giản [10]: ....................... 27

3.4.2. Tuyến tính hoá biểu thức có lực cản đáy: ............................ 27

3.4.3. Tuyến tính hoá biểu thức trao đổi nƣớc qua đê: .................. 28

3.5. Thuật giải hệ phƣơng trình đại số tuyến tính [1]; [5]; [6]: ................ 28

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]: ............................................. 33

3.6.1. Khái toán mặt cắt: ............................................................... 33

3.6.2. Tạo giá trị mực nƣớc và lƣu lƣợng làm điều kiện ban đầu: . 34

3.6.3. Vấn đề xác định hệ số nhám và chỉnh kết quả:.................... 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 ............................... 36

4.1. Mô hình kết nối Marine và Imech1D [5]: ........................................ 36

4.2. Xử lý số liệu cho mô hình kết nối Marine-IMech1D: ...................... 37

4.2.1. Xử lý bản đồ địa hình: ........................................................ 37

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: ......................................................... 38

4.2.1.2. Tạo mạng sông suối từ DEM: ................................ 40

4.2.1.3. Phân chia lƣu vực trên nền DEM: .......................... 41

4.2.2. Xử lý bản đồ phân loại đất: ................................................. 41

4.2.3. Xử lý bản đồ hiện trạng sử dụng đất: .................................. 43

4.2.4. Xây dựng bản đồ phân bố mƣa trong lƣu vực: .................... 43

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: ..................................................................................................... 45

4.2.6. Xử lý các số liệu khác:........................................................ 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 .......................................................................... 47

5.1. Quá trình cần đánh giá (ƣớc lƣợng) [13]: ......................................... 47

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]: ............................................................................................................... 49

5.2.1. Định nghĩa về các ƣớc lƣợng tiên nghiệm và hậu nghiệm: .. 49

5.2.2. Bƣớc dự báo – cập nhật (ƣớc lƣợng tiên nghiệm): .............. 49

5.2.3. Bƣớc hiệu chỉnh (ƣớc lƣợng hậu nghiệm): .......................... 51

5.2.4. Tìm Kalman gain (blending factor) K: ................................ 52

5.3. Thuật toán lọc Kalman rời rạc: ........................................................ 55

5.3.1.Cập nhật theo thời gian – dự báo (ƣớc lƣợng tiên nghiệm) (predict): ...................................................................................... 55

5.3.1.1. Phép tính s 1: ......................................................... 55

5.3.1.2. Phép tính số 2: ....................................................... 55

5.3.2. Cập nhật theo đo đạc – chỉnh sửa (ƣớc lƣợng hậu nghiệm) (correct): ...................................................................................... 56

5.3.2.1. Phép tính số 3: ....................................................... 56

5.3.2.2. Phép tính số 4: ....................................................... 56

5.3.2.3. Phép tính số 5: ....................................................... 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 ...................... 57

6.1. Kết quả kiểm tra bài toán mẫu cho 10 lƣu vực bộ phận: .................. 57

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: ......................................... 58

6.2.1. Nhận định chung tình hình lũ sông Đà năm 2006: .............. 58

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: .............................................. 58

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: ................................................................................... 61

KẾT LUẬN ...................................................................................................... 64

DANH MỤC CÔNG TRÌNH CỦA TÁC GIẢ ................................................. 65

TÀI LIỆU THAM KHẢO ................................................................................ 67

Tiếng Việt .............................................................................................. 67

Tiếng Anh .............................................................................................. 67

Tiếng Pháp ............................................................................................. 68

PHỤ LỤC ........................................................................................................ 69

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 1. sông Đà bằng MARINE: ........................................................................ 69

Kết quả kiểm định bộ chƣơng trình tính toán thủy lực một chiều 2. IMech1D bằng các bài toán kiểm định mẫu (Test Cases). ...................... 75

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 .................... 75

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 .......................................................................................... 76

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Ụ.............. 77

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 ..................... 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 ............................................ 79

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 ................................................ 80

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 ..................................................... 81

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 ........................................ 82

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 83

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 ......................................................................................... 84

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 ............................................................ 86

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 .................................................................... 89

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

để 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.

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-

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

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ũ

- 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ỹ)

- 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 đã

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ế

- 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.

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:

(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ì:

Với chất lỏng không nén đƣợc ta có , sử dụng công thức Green-

Ostrogradski

Từ (2.1) suy ra : (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:

(2.3)

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:

(2.4)

Trong đó: Pente: độ dốc; Km: hệ số nhám Manning

x: chiều rô ̣ng ô lƣớ i t : Bƣớ c thờ i gian tính j: Hƣớ ng chảy củ 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 t 1 đế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ô hình MARINE tính toán thấm dƣ̣a trên lý 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:

(2.5)

Trong đó: F(t) là độ sâu luỹ tích của nƣớc thấm vào trong đất

: 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

(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) cƣờng độ mƣa nhỏ hơn tốc độ thấm tiềm năng và mặt đất chƣa bão hòa. Quá trình nƣớc đọng trên mặt đất bắt đầu xảy ra khi cƣờng độ mƣa vƣợt quá cƣờng độ thấm tiềm năng (t = tp) lúc đó đất ở trạng thái bão hòa. Khi quá trình mƣa tiếp tục (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).Tốc độ thấm tiềm năng là một hàm của lƣợng thấm tích lũy F (t = tp).

Nƣớc đọng xuất hiện khi tốc độ thấm tiềm năng nhỏ hơn hoặc bằng cƣờng

độ mƣa i (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ó:

dẫn tới (2.7)

Sau khi có nƣớc đọng, lƣợng thấm tích lũy đƣợc tính theo công thức:

(2.8)

P0 (mƣa)

Thấm - Green Ampt

Cấu trú c lƣu vƣ̣c (đô ̣ cao đi ̣a hình, đất, thảm phủ)

Dòng chảy mặt

Lớ p nƣớc đọng

tại ô lƣới

Liên kết các ô lƣớ i, tính ̉i giƣ̃a các ô Q trao đổ

Dòng chảy sát mặt

Quá trình Q~t ta ̣i cƣ̉ a ra củ a lƣu 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.

- 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.

Bản đồ DEM 50m

1 1. Độ cao và độ dốc 2. Mạng lƣới sông 3. Phân chia tiểu lƣu

vực

Bản đồ Đất dạng số

2 1. Phân loại đất 2. Độ ẩm đất

3 1. Lớp phủ thực vật 2. Lƣợng trữ nƣớc mặt

Bản đồ Thảm phủ dạng số

Mƣa quan trắc

4 Phân bố mƣa theo không gian

Mặt cắt ngang

Đo đạc lòng sông

5 Quan trắc lƣu lƣợng vào

Dòng chảy trong sông

6

Hình 2.2: Các lớp thông tin yêu cầu của Marine

Hình 2.3: Cấu trúc các lớp dữ liệu trong Marine

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.

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:

(3.1)

(3.2)

Ở đâ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;

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:

(3.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 ô đó.

(3.4)

Ở đâ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à:

Từ phƣơng trình (3.1) ta có:

(3.5)

Để 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:

(3.6) a1 + a2 + a3 + a4 = 0

trong đó

Sai phân từng số hạng trong (3.6) ta đƣợc

Do vậy phƣơng trình sai phân của (3.2) có dạng sau:

(3.7)

Sai phân hóa phƣơng trình (3.4) ta có

(3.8)

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

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

đó:

(3.9) Qvào - Qra = 0

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):

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:

Ở đâ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): A = A(Z), Q2/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)

(3.11)

(3.12)

(3.13)

(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

3.4.2. Tuyến tính hoá biểu thức có lực cản đáy:

Trong số hạng 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:

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

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:

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

(chảy tự do từ i sang k). tràn là Ld. Khi đó, nếu

(3.15)

Nếu (chảy ngập)

(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.17)

- Trƣờng hợp chảy ngập (1.16)

(3.18)

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

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:

(3.19) a1Qd + a2Qc + a3Zd + a4Zc = a5

(3.20) b1Qd + b2Qc + b3Zd + b4Zc = b5

Từ (3.19), (3.20) dễ dàng suy ra:

(3.21) Qd = a6Zd + a7Zc + a8

(3.22) Qc = b6Zd + b7Zc + b8

Đ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:

(3.23)

Ở đâ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

(3.24)

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

(3.25) Qdi = Qi

Tại các nút biên hạ lƣu j cho giá trị cao trình mực nƣớc:

(3.26) Zcj = Zj

Hoặc tại nút biên k cho quan hệ giữa lƣu lƣợng và mực nƣớc

(3.27) Zk = Fk(Qk)

(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)

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 chứa (r +1) số hạng nên để khử đuổi biến Zi cần (r+1)2 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ỡ

. Để ý rằng khử Gauxơ có số phép tính là n2l (xem [14]). Nhƣ vậy nếu  n thì số phép tính theo phƣơng pháp khử đuổi giảm đáng kể.

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

4

Hì nh Hình 3.1 3.1 Hì 11 nh 3.2 Hì Hình 3.3 nh 3.3

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)

sdt=(6-12); (7-12); (7-15) khi đó dd=-1.Ta đƣợc hình 3.6

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 đƣợc hình 3.9

Hình 3.9

Khử nút 7 (bậc 3): sdb=(7-4); (7-8); (7-12)

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, 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ử lý của máy tính, tuy nhiên chỉ cần xác định 1 lần cho tất 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.

Phương án

Kết quả

Số nút khử đuổi Số nút khử Gauxơ Thời gian tính

404 699 2h50’ Khử đuổi nút có bậc  2

960 143 33’54’’ Khử đuổi nút có bậc  5

1053 50 32’50’’ Khử đuổi nút có bậc  10

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]:

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

Hình 3.11

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

Chiều rộng của mặt cắt:

Diện tích mặt cắt:

Môđul lƣu lƣợng:

Bán kính thủy lực:

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.

Zbien

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 n

Hình 3.12: Điều kiện đầu trong IMech1D

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

- 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.

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.

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ô

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.

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á

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:

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

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

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.

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:

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.

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:

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

ID Tên Trạm

Tỉnh Kinh độ Vĩ độ

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’

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’

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’

Lai Châu Lai Châu 103 0 9’ 22 0 3’

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’

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’

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’

Cò Nòi

Sơn La 104 0 9’ 21 0 8’

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’

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’

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’

Km46

Sơn La 104 0 59’ 20 0 47’ 36

Tà Nàng

Hòa Bình 104 0 27’ 20 0 56’

ảnh hƣởng mƣa thu đƣợc nhƣ hình 4.12.

Bảng 4.1: Vị trí các trạm mưa sử dụng trên sông Đà

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.

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.

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 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

(5.1)

Với đại lƣợng đo đạc thoả mãn phƣơng trình sai phân sau

(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

- 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 )

Q = Cov w =

Vì suy ra

Q = (5.3)

- 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

R = (5.4)

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ó

Hàm phân bố sác xuất có dạng

(5.5)

- uk là vectơ driven function l chiều, đại lƣợng tiền định,

- 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:

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ƣ

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.

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

(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ó:

Do đó

 Ma trận hiệp phƣơng sai của sai số ƣớc lƣợng tiên nghiệm xác định nhƣ

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:

(5.7)

 Tìm ƣớc lƣợng tiên mghiệm

Lấy kỳ vọng của đại lƣợng sau

Ta đƣợc:

(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

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 )

Cuối cùng ta thu đƣợc

(5.9)

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

Đâ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ự)

(5.10)

 Tìm ƣớc lƣợng hậu nghiệm

Trƣớc hết ta xét ý nghĩa của đại lƣợng . Có thể coi đây là giá trị “đo đạc

dự báo”. Và khi ấy đại lƣợng đƣợ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

(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

Chú ý tới tính chất

Lúc này ma trận hiệp phƣơng sai có dạng:

Cuối cùng, theo các định nghĩa, ta thu đƣợc

(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

Đâ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

hiệp phƣơng sai sai số hậu nghiệm .

Thật vậy:

Tóm lại, ta đã chứng minh đƣợc:

(5.13)

 Biến đổi

Trƣớc hết ta nhắc lại công thức

Với công thức trên ta có

Cuối cùng ta đƣợc

Với S có giá trị

 Tìm đạo hàm của tr Pk

Trƣớc hết ta có các công thức sau:

Ở đâ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

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ó:

Lấy đạo hàm từng thành phần ta đƣợc:

Tổng hợp lại ta đƣợc:

Chú ý tới tính đối xứng của các ma trận và S, ta thu nhận đƣợc công

thức cuối cùng:

Điều kiện cực tiểu của tr Pk:

Nhƣ vậy ta đã chứng minh và thu nhận đƣợc đầy đủ các bƣớc của lọc

Kalman

Một số tính chất của K:

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

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ó:

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

Trong đó:

– 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

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.

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 )

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

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

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:

2. Cập nhật theo quan trắc 2.1 Tính

1.Cập nhật theo thời gian 1.1 Tính

1.2 Tính

2.2 Tính

Ƣớc lƣợng ban đầu

2.3 Tính

Hình 5.2: Sơ đồ tính lặp quá trình lọc Kalman

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 (km2)

R (mm/h)

Tr (hours)

T (hours)

Vi (Trm3)

Vo (Tr m3)

Qi (m3/s)

Qo (m3/s)

 W (%)

Qmax (%)

LV

3621

144

1738.08 1730.57 0.43229 10058.3 10032.5 0.25684

1

2441

144

1171.68 1162.4 0.79161 6780.56 6770.97 0.14137

2

2255

144

811.8

809.535 0.27897 6263.89 6254.7

0.1467

3

898.8

144

323.568 322.867 0.21673 2496.67 2493.91 0.11041

4

2211

144

1193.94 1182.05 0.99588 6141.67 6135.35 0.10285

5

3245

144

1752.3

1755.7

-0.194 9013.89 9057.11

-0.4795

6

3530

144

2118

2100.51 0.82561 9805.56 9784.86 0.21106

7

2141

144

1284.6 1281.07 0.27506 5947.22 5970.91

-0.3983

8

3153

144

1513.44 1562.01

-3.2094 8758.33 9068.43

-3.5406

9

2889

144

1733.4 1725.48 0.45667

8025

8005.35 0.24486

10

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ả:

- 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:

Q (m3/s)

(h) (T db-T tt)

% sai sè Số cửa xả đáy Số cửa xả mặt Thời gian xuất hiện

Qmax thực đo vào hồ Hòa Bình (m3/s) Qmax tính toán dự báo vào hồ Hòa Bình (m3/s)

5500 5000 -500 -9.1 14 0 0 10h/1/VII/ 2006

6000 5300 -700 -11.7 12 17h/2/VII/ 2006

9000 8900 -100 -1.1 2 1 0 23h/10/VII /2006

1 5 0 14700 13808 -892 -6.1 12h/19/VII /2006

3600 3672 72 2 -4 0 0 23h/8/VIII /2006

3800 5457 1657 43.6 17 0 0 17h/17/VII I/2006

6000 6 4900 -1100 -18.3 0 0 5h/19/VIII /2006

4000 0 3238 -762 -19.1 0 0 1h/25/VIII /2006

9200 0 8900 -300 -3.3 4 0 13h/12/X/ 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

Thuc Do

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

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.

Ốp thời gian dự kiến Mức sai số Số điểm đúng/Tổng số P%

18.5 %Q 97/116 82.5 6 h

18.5 %Q 103/116 83.6 12h

18.5 %Q 94/116 81 18h

18.5 %Q 89/116 76.7 24h

21.5%Q 88/115 76.5 30h

24.5%Q 81/115 70.4 36h

27.5%Q 83/115 72.1 42h

30.5%Q 83/115 72.1 48h

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

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 Đà.

Tính toá n

Hình 6.4: Đường quá trình dự báo lưu lượng vào hồ Hòa Bình ngày

Tính toá n

05/07/2009

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

Tính toá n

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

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

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ệ.

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

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.

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.

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.

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

Qo (m3/s)

Qi (m3/s)

Vo (Tr m3)

 W (%)

Qmax (%)

A (km2) 3621 2441 2255 898.8 2211 3245 3530 2141 3153

R (mm/h) 10 10 10 10 10 10 10 10 10

Tr (hours) 48 48 36 36 54 54 60 60 48

T (hours) 144 144 144 144 144 144 144 144 144

Vi (Trm3) 1738.08 1730.57 0.43229 10058.3 10032.5 0.25684 1171.68 1162.4 0.79161 6780.56 6770.97 0.14137 0.1467 809.535 0.27897 6263.89 6254.7 811.8 323.568 322.867 0.21673 2496.67 2493.91 0.11041 1193.94 1182.05 0.99588 6141.67 6135.35 0.10285 -0.4795 1752.3 -0.194 9013.89 9057.11 1755.7 2100.51 0.82561 9805.56 9784.86 0.21106 2118 -0.3983 1284.6 1281.07 0.27506 5947.22 5970.91 -3.5406 -3.2094 8758.33 9068.43 1513.44 1562.01

LV 1 2 3 4 5 6 7 8 9

2889

144

1733.4 1725.48 0.45667

8025

8005.35 0.24486

10

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:

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%.

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)

- 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

- 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:

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

điểm, nghĩa là sao cho

Tại đoạn nhập lƣu: Không có độ dốc, I=0 -

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

quy luật Q = a(H - Hs)3/2 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.

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,

Q(m3/g y)

Qmax

- 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

(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. Đã kiểm định bộ chƣơng trình IMech1D cho trƣờng hợp Q0=1.000m3/s,

Qmax=2.000 m3/s , Hhadu = 2,5 m T = 30s

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:

Q(m3/g y)

Qmax

Lƣu lƣợng

Q0, H0.

- Đ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: - 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 Q0=750m3/s, H0=2,56 m, Qmax= 2.000 m3/s.

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:

Q(m3/g y)

Cho lƣu lƣợng Q(0, t) biến thiên theo quy luật sau:

Qmax Q0

- Điều kiện biên hạ lƣu: cho quan hệ f(QL, HL).

- Điều kiện ban đầu:

Q(x, 0) = Q0

H(x, 0) = H0

Đã kiểm định bộ chƣơng trình IMech1D cho trƣờng hợp Q0=1.000m3/s; Qmax=2.000 m3/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:

Biên thƣợng lƣu:

ở đây: a là hằng số cho trong quan hệ lƣu lƣợng và mực nƣớc tại hạ lƣu.

H0, H1 là các hằng số trong biểu thức biến đổi mặt thoáng:

Kết quả tính toán bằng bộ chƣơng trình IMech1D quá trình mực nƣớc hồ và quá trình lƣu lƣợng tại địa điểm x=0m, 2000m, phù hợp với kết quả của nghiệm giải tích.

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

Trường hợp 1

Mục đích:

- Thử sơ đồ. - Thử định tính và định lƣợng đƣờng mực nƣớc dừng xung quanh vị trí nhiễu cục bộ.

Trường hợp 2

Mô tả bài toán mẫu kiểm định :

Xét dòng chảy trong kênh

Trường hợp 3

chữ nhật có các đặc trƣng sau:

- Kênh hình chữ nhật, độ dài L, chiều rộng B, Zđ- cao độ đáy thƣợng lƣu. - Chiều rộng có những điểm thay đổi (rộng ra hoặc hẹp vào) nhƣ trên hình vẽ.

- Độ dốc I đều. - Ma sát đều (K = const) - Điều kiện biên:

+ Thƣợng du: cho lƣu lƣợng

+ Hạ du: cho chiều cao cột nƣớc nhƣ trong dòng đều H0.

Giải tích số và đánh giá định tính:

Phƣơng trình Saint - Venant trong trƣờng hợp dòng chảy dừng:

(vẫn coi ; hoặc không có ma sát trên thành thẳng đứng)

Khi B biến đổi theo X ta có:

Nếu biết B(x) có thể giải phƣuơng trình này bằng phƣơng pháp Runge-

Kutta và xác định đƣợc cao trình mực nƣớc.

Đã kiểm định chƣơng trình tính toán IMech1D cho các trƣờng hợp sau:

L = 5.000 m

B = 100 m

I = 0,0005 Q = 1.000 m3/s

H0 = 5 m

Zđ = 10 m

X2 = 4.000 m K = 30,59 (M1/3/s)

L = X2 - X1 = X3 - X2

Kết quả kiểm định cho thấy IMech1D đã vƣợt qua đƣợc bài toán mẫu số 9.

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

Mục đích:

- Thử sơ đồ.

- Thử định tính đƣờng mực nƣớc trong dòng chảy dừng vời hình học

không đều.

Mô tả bài toán mẫu:

Hình học:

- Kênh hình thang đều, đáy có chiều rộng BD, mặt có chiều rộng BM tại

chiều cao HM. - Độ dốc đều, I. - Hình học không đều:

+ Tiết diện mở rộng (hoặc thu hẹp) một cách tuyến tính sao cho tiết diện đạt giá trị lớn nhất (hoặc nhỏ nhất) tại điểm giữa của đoạn sông.

+ Từ điểm giữa ta thu hẹp (hoặc mở rộng) tiết diện sao cho trở

lại giá trị ban đầu tại điểm cuối của đoạn.

+ Số Stricker thay đổi sao cho đạt giá trị lớn nhất (hoặc nhỏ

nhất) tại điểm giữa đoạn sông.

- Điều kiện biên:

+ Thƣợng du: cho lƣu lƣợng Q.

+ Hạ du: cho mực nƣớc (nhƣ trong trƣờng hợp đều) H0.

Đánh giá định tính:

- Việc mở rộng diện tích mặt cắt ở giữa đoạn sông và việc giảm hệ số nhám (tăng hệ số Stricker) có cùng hiệu ứng giống nhau về mặt định tính đối với đƣờng mực nƣớc: cột nƣớc giảm về phía thƣợng du.

L/2 2

BM1

BM2

- Ngƣợc lại, việc thu hẹp diện tích mặt cắt sẽ tƣơng đƣơng với việc tăng hệ số nhám (giảm hệ số Stricker): cột nƣớc tăng về phía thƣợng du.

Kết quả kiểm định:

Đã kiểm định bộ chƣơng trình IMech1D cho trƣờng hợp sau:

L = 8.000 m

I = 0,00025

Zf = 60 m Q = 2.000 m3/s

H0 = 9,352

BM = 200 m

BD = 50 m

HM = 10 m

BM1 = 400 m

BM2 = 100 m

K = 35

K1 = 55

K2 = 15

Các kết quả tính toán đƣờng mực nƣớc cho trƣờng hợp mở rộng (thu hẹp) lòng kênh và tăng (giảm) hệ số Strickler ở giữa đoạn sông phù hợp với kết quả mẫu của tác giả bài toán mẫu.

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

Mục đích:

- Thử sơ đồ. - Thử định tính sự truyền sóng lũ trong kênh có lòng dẫn phức hợp(có

lòng và bãi),kiểm tra hiệu ứng vòng Q/H.

Mô tả:

Hình học:

- Lòng dẫn phức hợp (có lòng và bãi), hình lăng trụ đều (hình 2.38) - Có hệ số ma sát trên lòng và bãi khác nhau. - Độ dốc đều, I.

Điều kiện biên:

- Thƣợng du: cho lƣu lƣợng Q(0,t).

Qmax

- Hạ du: cho quan hệ lƣu lƣợng mực nƣớc trong trƣờng hợp dòng đều.

Bãi

Lòng

Đánh giá định tính:

- Tốc độ lan

truyền sóng lũ - Độ giảm đỉnh lũ - Hiệu ứng vòng

Q/H

Số liệu:

L = 400 km

I = 0,0001

Zf = 57 m Q0 = 300 m3/s Qmax = 1000 m3/s Chiều rộng lòng B1=200 m

Chiều rộng bãi B2=1000 m

Hệ số Strickler của lòng k1 = 35

Hệ số Strickler của bãi k2 = 15

T=4h

Dt=300s

Dx=500m

Modul lƣu lƣợng đƣợc tính theo công thức nhƣ sau:

ở đây Ai, Ri, Xi-lần lƣợt là diện tích ƣớt,bán kính thuỷ lực và chu vi ƣớt của

lòng dẫn thứ i,H-độ sâu mực nƣớc; .

Nếu đặt , trong đó k, A,  là hệ số Strickler, diện tích ƣớt và

chu vi ƣớt của lòng dẫn phức hợp, thì suy ra:

, nghĩa là k không còn là hằng số mà là một hàm của H. Ta

thu đƣợc hàm k và K phụ thuộc vào H dƣới dạng bảng.

Để tính vận tốc của từng lòng dẫn, ta làm nhƣ sau:

Ký hiệu -Độ dốc thuỷ lực,ta có:

-lƣu lƣợng của lòng dẫn thứ i

Do . ,với Ui-vận tốc lòng dẫn thứ i. Suy ra:

Để tìm quan hệ Q(H) của dòng đều,ta xuất phát từ hệ thức:

I=Sf hay

Nếu biết Q, từ công thức trên suy ra K và tra từ bảng sẽ tìm đƣợc H.

Trên hình 2.39 là kết quả tính quá trình lƣu lƣợng bằng bộ chƣơng trình

IMech1D tại x=0, x=5km, x=10km.

Trên hình 2.40 là kết quả tính quá trình mực nƣớc bằng bộ chƣơng trình

IMech1D tại x=0, x=5km, x=10km.

Trên hình 2.41 là đƣờng quan hệ giữa mực nƣớc và lƣu lƣợng thu đƣợc

bằng bộ chƣơng trình IMech1D tại x=5km.

Trên hình 2.42 là kết quả tính vận tốc của dòng chảy bàng bộ chƣơng trình

IMech1D vùng lòng sông, bãi sông tại x=5km.

Từ các kết quả tính của IMech1D có thể nhận xét nhƣ sau:

- Vì biên dƣới là quan hệ Q(H) của dòng đều, nếu ở biên trên Q=const thì ta đƣợc dòng đều với độ sâu mực nƣớc H tƣơng ứng với Q ở biên. Tuy nhiên do tính không dừng của dòng chảy nên khi lũ lên, với cùng H thì Q lớn hơn so với dòng đều, khi lũ xuống thì ngƣợc lại.

- Các thời điểm Vmax, Qmax, Hmax không xẩy ra đồng thời mà kế tiếp nhau: V đạt max khi lũ đƣơng lên, tiếp đến là lƣu lƣợng Q, sau cùng là mực nƣớc H. Nhƣ vậy, khác với dòng đều, mực nƣớc max không tƣơng ứng với Q max.

Kết quả này khá phù hợp với kết quả do các phòng thí nghiệm thủy lực của châu Âu đƣa ra cho bài toán về dòng chảy không dừng có lòng dẫn phức hợp với biên Q_H.

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

Mục đích:

- Thử sơ đồ. - Thử kiểm tra sự phân bố lƣu lƣợng tại phân lƣu.

Mô tả:

Hình học:

- Thƣợng du: có 1 nhánh sông hình

chữ nhật đều nhau.

Hình 2.43

- Hạ du: Có 2 nhánh với hình học

thay đổi.

- Trong trƣờng hợp cơ sở: hình học cơ sở của hệ đều nhau và sao cho nếu gộp cả 2 nhánh ở hạ lƣu lại thì hoàn toàn giống vơí nhánh ở thƣợng du (cùng chiều rộng, cùng độ dốc).

- Cho ma sát đáy.

Điều kiện biên:

- Thƣợng du: cho lƣu lƣợng. - Hạ du: cho mực nƣớc chuẩn trong trƣờng hợp chảy đều.

Đánh giá định tính:

- Trường hợp cơ sở: Do dòng chảy đều, phân bố lƣu lƣợng sẽ tỉ lệ thuận

với chiều rộng đoạn sông.

- Trường hợp thay đổi hình học: 2 nhánh hạ du sẽ làm thay đổi phân bố

lƣu lƣợng có thể định tính trƣớc.

Số liệu:

Trường hợp cơ sở:

B2 = 150 m L1 = 2.000 m

B3 = 50 m

K = 45 L2 = L3 = 3.000 m I = 10-4

H01 = 1 m

H02 = 5,41 m

Dt = 60 s

Dx = 100 m Zf = 10 m Q1 = 90 m3/s Q2 = 1.500 m3/s B1 = 200 m

Trường hợp thay đổi:

- Phương án 1: Giảm tuyến tính chiều rộng của nhánh 3 từ 50m xuống 25m tại điểm x=500m, sau đó lại tăng tuyến tính trở về 50m ở cuối nhánh.

- Phương án 2: Tại nhánh 3 cho độ dốc I=0, chiều cao đáy là 9,8m. - Phương án 3: Giảm chiều rộng của nhánh 2 từ 150m xuống 50m tại hạ

du.

Kết quả kiểm định :

Kết quả tính toán các phƣơng án bằng bộ chƣơng trình tính toán IMech1D khi Q1=90m3/s, H01=1 m và khi Q2=1.500m3/s, H02=5,41m thể hiện trên các bảng (dƣới đây)

Lưu Lưu Cột

lượng nhánh 2 lượng nhánh 3 nước tại hợp lưu Cột nước tại thượng du

Trƣờng hợp cơ sở 67.5 1.003 1.001 22.5

Phƣơng án 1 73.504 1.004 1.033 16.496

Phƣơng án 2 70.394 1.017 1.007 19.606

Phƣơng án 3 57.042 1.082 1.154 32.958

Lưu Lưu Cột

lượng nhánh 2 lượng nhánh 3 nước tại hợp lưu Cột nước tại thượng du

Trƣờng hợp cơ sở 1125 375 5.416 5.413

Phƣơng án 1 1219.195 280.805 5.456 5.464

Phƣơng án 2 1142.578 357.422 5.425 5.421

Phƣơng án 3 829.352 670.648 5.916 5.971

Qua kết quả trên có thể nhận xét nhƣ sau:

- ở trƣờng hợp cơ sở, nghiệm chính xác là dòng đều. Ta thấy kết quả tính

bằng IMech1D là phù hợp cả về định tính lẫn định lƣợng.

- ở các phƣơng án 1-3, kết quả tính bằng IMech1D phù hợp về định tính: chiều cao cột nƣớc tăng lên ở đoạn có mặt cắt mở rộng và giảm ở đoạn có mặt cắt thu hẹp.