Ảnh hưởng của kích thước ô lưới tới kết quả tính toán thủy lực dòng chảy lũ

 KHOA HỌC CÔNG NGHỆ 
 ẢNH HƯỞNG CỦA KÍCH THƯỚC Ô LƯỚI TỚI 
 KẾT QUẢ TÍNH TOÁN THỦY LỰC DÒNG CHẢY LŨ 
 Lê Thị Thu Hiền, Lê Thanh Hùng 
 Trường đại học Thủy Lợi 
Tóm tắt: Kích thước ô lưới là một trong những yếu tố quan trọng ảnh hưởng tới tính chính xác 
của kết quả tính thủy lực bằng mô hình toán. Sử dụng phần mềm tính toán thủy lực 2D-FV do 
các tác giả tự xây dựng dựa trên phương pháp thể tích hữu hạn để giải hệ phương trình nước 
nông phi tuyến hai chiều (2D-NSWE) trên lưới có cấu trúc. Bài báo phân tích, đánh giá ảnh 
hưởng của kích thước ô lưới tới kết quả tính toán thủy lực như: độ sâu, lưu lượng của dòng chảy 
lũ trên địa hình phức tạp. 
Từ khóa: Kích thước lưới, mô hình toán, chương trình 2D-FV, đặc trưng thủy lực 
Summary: Grid size is one of the most important factors affect to the accuracy of hydraulic 
results obtained by mathematical models. Using the 2D-FV hydraulic computation software built 
by the authors based on the finite volume method to solve 2D nonlinear shallow water equations 
(2D-NSWE) on a structured grid. The paper estimates influencing of grid size on hydraulic 
results such as water depth, dischage of flood flow on complex terrain. 
Keywords: Grid size, numerical model, 2D-FV model, hydraulic character 
1. ĐẶT VẤN ĐỀ * với thực đo chứ không phải lưới có kích thước 
Sử dụng mô hình toán mô phỏng các bài toán nhỏ hơn 30m, [2]. Lê Thanh Hùng (2017) lại 
thủy lực như: sự lan truyền sóng lũ trên địa khảo sát sự ảnh hưởng của 4 loại ô lưới tới kết 
hình có độ dốc phức tạp; sóng gián đoạn đã quả tính sự lan truyền sóng lũ do tình huống 
được thực hiện rộng rãi trong thủy lợi. Với vỡ đập Nậm Chiến (Sơn La), [3]. Về lý thuyết, 
những bài toán thủy lực hai chiều hay ba chia lưới có kích thước càng nhỏ thì kết quả 
chiều, kích thước ô lưới của miền tính toán càng chính xác, tuy nhiên, khi số ô tính toán 
luôn là yếu tố có ảnh hưởng lớn tới độ chính càng lớn thì thời gian tính sẽ càng nhiều, thậm 
xác của kết quả tính bằng mô hình toán. Đặc chí các máy tính cá nhân thông dụng không xử 
biệt, trên những địa hình phức tạp, cao độ đáy lý được. Vì vậy, việc đánh giá ảnh hưởng của 
biến đổi nhiều. các kích thước ô lưới tới các dạng kết quả thủy 
 lực như: mực nước, lưu lượng để tìm ra kích 
Valiani và nnc (2002) đã dùng hai loại ô lưới thước ô lưới hợp lý là rất cần thiết. 
thô và mịn chia lưu vực hồ Malpasset (Pháp) 
để tính quá trình mực nước tại các điểm Trong nội dung bài báo này, các tác giả dùng 
nghiên cứu, [1]. Wang (2011) chia miền tính một số ví dụ để đánh giá ảnh hưởng của kích 
toán cũng của hồ này thành lưới vuông với các thước ô lưới tới kết quả tính thủy lực như mực 
kích cỡ khác nhau để tính sự lan truyền lũ và nước, lưu lượng dòng chảy lũ bằng chương 
kết luận rằng, lưới 40m cho kết quả tốt nhất so trình tính toán 2D-FV do chính các tác giả xây 
 dựng bằng ngôn ngữ lập trình Fortran. 
Ngày nhận bài: 26/7/2017 2. PHƯƠNG PHÁP NGHIÊN CỨU 
Ngày thông qua phản biện: 06/9/2017 Sử dụng chương trình tính 2D-FV dựa trên 
Ngày duyệt đăng: 26/9/2017 
 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 107 
 KHOA HỌC CÔNG NGHỆ 
phương pháp thể tích hữu hạn (FVM) loại 3. KẾT QUẢ TÍNH TOÁN 
Godunov giải hệ phương trình nước nông phi 3.1. Sóng dao động trong chảo parabol có 
tuyến (NSWE) trên lưới Catersian đã được giới ma sát 
thiệu trong [3], [4], [5]. Điều kiện ổn định của 
phương pháp số được lấy theo Courant- Với miền tính toán có kích thước 10000m 
Fredrichs-Lewy (CFL). Với lưới Cartesian, bước 10000m, địa hình đáy mô tả bằng phương 
thời gian t bị ràng buộc bởi phương trình (1): trình (2) được chia làm 4 loại ô lưới có kích 
 thước lần lượt 25m; 50m; 100m và 200m. Kết 
 1
 ~ ~ ~ ~ quả của phương pháp số được đưa ra là mực 
 u gh v gh 
 Δt Cr max , (1) nước tại các thời điểm t = 1000s và t = 6000s 
 Δx Δy 
 được so sánh với kết quả chính xác của 
trong đó: nghiệm giải tích. 
 2
Cr là số Courant có giá trị trong khoảng 0 < Cr zb x h0 x/a , (2) 
 1, các kết quả tính toán trình bày ở mục 3, hệ 
 trong đó các thông số h0 và a là các hằng số. 
số Cr được lấy bằng 0,9. Nghiệm giải tích được đưa ra lần đầu bởi 
 ~
u~, v~, h là các giá trị lưu tốc và độ sâu trung Sampson và nnc [7] phụ thuộc vào hệ số ma 
bình tính theo Roe, (1981), [6]. sát đáy , trong trường hợp này số hạng ma sát 
 S f τ  h u . 
Công thức (1) chỉ ra rằng khi x và y càng 
nhỏ thì bước thời gian tính t sẽ càng nhỏ, nên Nghiệm giải tích về mực nước dao động của 
thời gian chạy máy tính sẽ càng lâu. bài toán này tuân theo biểu thức (3): 
 2 τt τt/ 2
 2 2 B e e τB 
η x,t h0 A. sτ sin 2st τ / 4 s cos 2st Bs cos st sin st x , (3) 
 4g g 2 
 2 2 τt
 2 2 2 a B e
trong đó thông số B là hằng số và s p τ / 2 ; p 8gh0 /a , A 2 
 8g h0
 30
 30
 T=6000s - x=100m
 T=6000s - x = 50m
 bed profile
 20 bed 20 analytical
 analytical numerical
 numerical
 cao độ(m)
 cao độ(m) 10 10
 0
 0
 -5000 -3000 -1000 1000 3000 5000
 -5000 -3000 -1000 x(m) 1000 3000 5000 x(m) 
 Hình 1: Sự dao động mực nước trong chảo parabol tính với hai loại lưới 100m và 50m 
Chỉ số Nash-Sutcliffe (E) được sử dụng để quả chính xác tương ứng với các kích thước 
đánh giá độ chính xác của kết quả tính mực lưới khác nhau: 
nước theo chương trình 2D-FV so với kết 
108 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 
 KHOA HỌC CÔNG NGHỆ 
 n
 X X 2 Bảng 1: Chỉ số Nash tương ứng với các 
 E 1 i 1 analytical,i model,i kích thước lưới khác nhau 
 n 2 , (4) 
 X X
 i 1 analytical,i analytical E (%) E (%) 
 x Số ô lưới 
 (t=1000s) (t=6000s) 
trong đó Xanalytical,i là giá trị chính xác tính theo 
 25 400 400 99,47 98,65 
(3), còn Xmodel,i là nghiệm tính theo phương 
pháp số ở các vị trí khác nhau i. 50 200 200 99,43 98,29 
 100 100 100 98,64 98,04 
Từ kết quả ở bảng 1 cho thấy: Với 4 kích 200 50 50 98,12 97,79 
thước lưới từ mịn đến thô: 25m; 50m; 100m và 
200m chỉ số Nash tại thời điểm t = 1000s và t 3.2. Sóng gián đoạn trên bề mặt có mỏm núi 
= 6000s không có sự khác biệt lớn. Nghĩa là 
lưới thô vẫn có thể cho kết quả mô phỏng mực Để đánh giá ảnh hưởng của kích cỡ lưới tính 
nước tương đối tốt. Vì vậy, trong những bài toán tới hai loại kết quả là mực nước, lưu lượng, 
toán có địa hình phức tạp, nếu chỉ cần kết quả ví dụ dưới đây do Castro và nnc (2009) giới 
liên quan đến mực nước có thể dùng lưới thô thiệu được sử dụng, [8]. Trong miền tính toán 
để tính. hình chữ nhật có kích cỡ [0; 2]m × [0; 2]m, cao 
 độ đáy được cho bởi phương trình (5): 
 1 2 2 2
 cos 2 (x 0,5) 1  cos 2 y 1  nêu x 1,5 y 1 0,5m 
 zb 8 (5) 
 2 2 2
 0 nêu x 1,5 y 1 0,5m 
Vận tốc ban đầu bằng 0, biên đóng tại vị trí x = 0m và x = 2,0m; hai biên còn lại là biên mở. 
Mực nước ban đầu được cho bởi biểu thức (6). 
 2 2 2 3
 1,1 nêu x 1,25 y 1 0,1 0,2m /s.m trong khi ảnh hưởng của hai kích 
(x, y) (6) thước lưới này tới mực nước tại cùng vị trí 
 0,6 nêu x 1,25 2 y 1 2 0,1 2
 này không đáng kể. 
Cắt dọc hình dạng cột nước ban đầu và địa 
hình đáy được thể hiện trên hình 2. 
Trên hình 3 cho thấy: Kết quả tính bằng 
phương pháp số theo chương trình 2D-FV do 
các tác giả xây dựng trùng khớp với kết quả 
của Castro và nnc (2009), [8]. Mặt khác, 
trong kết quả tính toán quá trình mực nước 
với 3 kích thước lưới: 0,0050m; 0,0025m; 
0,0020m cho thấy không có sự sai khác nhiều 
tại vị trí đáy có vật cản. 
Tuy nhiên, kết quả tính lưu lượng đơn vị lại 
có sự khác biệt rõ rệt tại vị trí này. Tại x = 
1,35m, chênh lệch kết quả tính lưu lượng đơn 
vị của lưới thô nhất và mịn nhất lên tới - Hình 2: Mực nước ban đầu và địa hình đáy 
 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 109 
 KHOA HỌC CÔNG NGHỆ 
 0.8
 t = 0.15s 0.4
 t = 0.15s
 0.7
 0.2
 0.6 /s/m)
 3
 0
 0.5 bed
 elevation(m) dx=0.005m
 dx=0.0025m -0.2discharge(m dx=0.005m
 0.4 dx=0.0025m
 dx=0.002m
 dx=0.002m
 0.3 -0.4
 0.4 0.6 0.8 1 x(m) 1.2 1.4 1.6 1.8 0.4 0.6 0.8 1 x(m) 1.2 1.4 1.6 1.8 
 Hình 3:Quá trình mực nước và lưu lượng úng với x = 0,005m; 0,0025m và 0,002m 
3.3. Sóng vỡ đập trên địa hình phức tạp nhau. Sau 10s dòng chảy do vỡ đập lan truyền 
Ví dụ này dùng để kiểm tra khả năng của mô tới hồ chứa 2 và sau 40s dòng nước đã xuống 
hình trong việc xử lý các vấn đề khô ướt do hồ 3 và đến 500s, mực nước trong 3 hồ đạt vị 
 trí ổn định. Kết quả này chỉ ra rằng trong điều 
địa hình biến đổi phức tạp. Địa hình được định 
 kiện địa hình rất phức tạp, mô hình toán vẫn 
nghĩa bằng hệ các phương trình sau để chia 
 có khả năng cho kết quả hợp lý khi không sinh 
miền tính toán thành 3 lòng hồ. 
 ra nhiễu động. Bảng 2 lại chỉ ra các kết quả 
 zb x, y min[zb1 x, y ,zb2 x, y , zb3 x, y ], (7) tính tương ứng với 3 kích thước lưới này, bao 
trong đó: gồm: thời gian chạy máy tính; lưu lượng lớn 
 nhất và mực nước lớn nhất tại đập; thời gian 
 2 2
 x 250 y lan truyền sóng lũ tới điểm nghiên cứu P(x = 
 zb1 x, y ;
 1600 400 250m; y = 0m) (hình 4). Với kích thước lưới 
 x2 y 50 2 2,5m, thời gian chạy máy lên tới 40 phút so 
 z x, y ; 
 b2 225 225 với 12 phút của lưới lớn gấp đôi 5,0m. Lưới 
 x 250 2 y 2 8m và 10m thời gian này chỉ là 3 phút và 2 
 z x, y 10
 b3 1225 225 phút. Bên cạnh đó, so sánh các kết quả thủy 
 lực với kích thước nhỏ nhất x = 2,5m, phần 
Hai vật cản được mô phỏng bằng biểu thức: 
 trăm sai số được xác định theo biểu thức (9): 
 2 2 
 x 250 y X xi X x min
 zbB1 x, y max zb x, y , 80 % 100% , (9) 
 50 50 err
 (8) X x min
 2 2 3
 zbB2 x, y 10neu x 200 y 10  10 trong đó: 
Với điều kiện ban đầu, mực nước trong hồ là X là đặc trưng thủy lực; X xmin là đặc trưng 
35m ở vị trí x -100m, phần còn lại đáy khô thủy lực tính với kích thước lưới nhỏ nhất. 
(hình 4). Tất cả các biên là biên đóng với hệ số Bảng 2 cho thấy, sai số của lưu lượng lớn 
nhám Manning bằng 0,033. Thời gian tính nhất tại đập ứng với các kích thước lưới 5,0m 
toán là 500s. Kết quả tính sự lan truyền sóng so với lưới mịn nhất 2,5m nhỏ hơn nhiều 
vỡ đập trên miền tính toán có kích thước (-2,33%) so với lưới 8,0m và 10,0m (-6,92% 
1000m 400m được tính với 3 kích thước lưới và +11,07%). Vì vậy, trong trường hợp này 
khác nhau 2,5m; 5m và 10m. Hình 4 chỉ ra sự lưới 5,0m có thể được dùng để tính lan truyền 
lan truyền sóng vỡ đập tại các thời điểm khác sóng lũ. 
110 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 
 KHOA HỌC CÔNG NGHỆ 
 P 
 Hình 4: Sự lan truyền sóng gián đoạn trên kênh có địa hình phức tạp tại 
 t = 0s; 10s; 100s và 500s ứng với kích thước lưới x = y = 2,5m 
 Bảng 2: Các đặc trưng thủy lực ứng với các kích thước lưới khác nhau 
 Thời gian Thời gian 
 Q Sai số H Sai số Sai số 
 x (m) chạy máy max max lũ đến 
 (m3/s) (%) (m) (%) (%) 
 (phút) (s) 
 (400 160) 40 9799,29 20,33 18,82 
 (200 80) 12 9570,68 -2,33 20,05 -1,38 20,34 +8,08 
 (125 50) 3 9121,20 -6,92 19,79 -2,65 19,95 +6,00 
 (100 40) 2 10883,96 +11,07 18,09 -11,02 12,82 -31,88 
3.4. Kịch bản vỡ đập vòm Nậm Chiến 945m, hạ lưu khô và coi đập vòm vỡ tức thời, 
Đập vòm Nậm Chiến trên suối Chiến, huyện hoàn toàn. Kết quả tính lưu lượng đỉnh lũ lớn 
Mường La, Sơn La là đập vòm duy nhất ở Việt nhất tại đập được thể hiện trên bảng 3. 
Nam tính đến nay. Lê Thanh Hùng (2017) đã 
dùng phương pháp số tính sự lan truyền sóng 
lũ theo kịch bản vỡ đập hoàn toàn trên miền 
tính toán có kích thước lên đến 14000m 
12000m, [3]. Vì vậy, việc tìm ra một kích 
thước lưới phù hợp không tốn nhiều thời gian 
chạy máy tính mà vẫn cho được kết quả hợp lý 
là cần thiết. 
Với bản đồ DEM 90m 90m khu vực lòng hồ 
Nậm Chiến như hình 5. Nội suy chia miền tính 
toán này ra làm 7 kích thước lưới khác nhau: 
 x = 10m; 20m; 25m; 30m; 40m; 50m và 
90m. Mực nước ban đầu trong hồ ở cao trình Hình 5: Bản đồ DEM lòng hồ Nậm Chiến 
 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 111 
 KHOA HỌC CÔNG NGHỆ 
 Bảng 3: Lưu lượng đỉnh lũ lớn nhất tại đập với các kích thước lưới khác nhau 
 3 3 3 3
 x Qmax.10 (m /s) (n=0,06) % sai số Qmax.10 (m /s) (n=0,04) % sai số 
 10 114,41 115,55 
 20 123,11 +7,6 124,58 +7,8 
 25 121,78 +6,4 123,23 +6,6 
 30 128,65 +12,4 128,73 +11,4 
 40 119,60 +4,5 120,00 +3,8 
 50 132,26 +15,6 133,83 +15,8 
 90 106,52 -6,9 122,83 +6,3 
Rõ ràng, lưới càng thô, số ô tính toán càng ít kết quả khác biệt nhất so với lưới 10m (hơn 
thì thời gian tính toán sẽ càng giảm. Tuy 15% ứng với cả hai độ nhám), trong khi lưới 
nhiên, từ những kết quả trên cho thấy, không thô nhất 90m lại cho kết quả thiên nhỏ. Vì vậy, 
phải lưới có kích thước càng nhỏ sẽ cho kết ta có thể tính toán sự lan truyền sóng lũ cho cả 
quả càng chính xác. Lưới 20m và 30m cho sai miền tính toán lớn bao gồm cả hạ lưu hồ Nậm 
số là +7,6% và +12,4% tương ứng với độ Chiến khi x = 40m. Hình 6 là kết quả tính sự 
nhám n = 0,06 và kết quả này là +7,8% và lan truyền sóng lũ do vỡ đập vòm Nậm Chiến, 
+11,44% khi n = 0,04. Tuy nhiên, khi x = tính với độ nhám n = 0,06 tại hai thời điểm t = 
40m, giá trị sai số này chỉ là +4,5% và +3,8% 500s và t = 1500s tương ứng với kích thước 
tương ứng với 2 độ nhám trên. Lưới 50m cho lưới x = y = 40m. 
 1
 1
 Hình 6: Bản đồ ngập lụt theo tình huống vỡ đập Nậm Chiến lúc t = 500s và 1500s 
4. KẾT LUẬN kết quả tính mực nước chỉ ra rằng có sự khác 
Độ chính xác của kết quả tính theo phương biệt không nhiều giữa kích thước lưới mịn và 
pháp số luôn là mục tiêu hàng đầu khi mô lưới thô (mục 3.1, 3.2). Tuy nhiên, khi xét ảnh 
phỏng các hiện tượng thủy lực theo phương hưởng này tới lưu lượng thì có sự khác biệt 
pháp này. Ảnh hưởng của kích thước ô lưới tới đáng kể, đặc biệt khi địa hình thay đổi phức 
112 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 
 KHOA HỌC CÔNG NGHỆ 
tạp như ví dụ trong mục 3.2 và 3.3. Áp dụng cho kết quả tốt hơn cả lưới nhỏ hơn là 30m và 
chương trình tính 2D-FV cho một trường hợp 20m và cả lưới lớn hơn 50m, 90m khi so sánh 
thực tế là vỡ đập hồ Nậm Chiến, Sơn La tương với kết quả tính theo kích thước lưới mịn nhất 
ứng với 7 kích thước lưới khác nhau chỉ ra 10m. Việc tìm ra kích thước lưới hợp lý sẽ đưa 
rằng: Sự ảnh hưởng của kích thước lưới tới lưu ra được kêt quả có độ sai số cho phép, đồng 
lượng lớn nhất tại đập là rất rõ ràng. Lưới 40m thời tiết kiệm thời gian chạy máy tính. 
TÀI LIỆU THAM KHẢO 
[1] A. Valiani, V. Caleffi and A. Zanni (2002). Case study: “Malpasset Dam-break Simulation 
 using a two dimensional finite volume method”. J. Hydraulic Engineering. (5) 128, 460- 472. 
[2] Wang. Y (2011). “Numerical Improvements for Large-Scale Flood Simulation”. Thesis of 
 Doctor Philosophy of Newcastle University. 
[3] Lê Thanh Hùng (2017). “Nghiên cứu sự lan truyền sóng lũ tới hạ lưu công trình trong tình 
 huống vỡ đập vòm Nậm Chiến bằng mô hình toán”. Tạp chí Khoa học và Công nghệ Thủy 
 lợi, 38, 87-94. 
[4] Le T.T.H (2014). “2D Numerical modeling of dam break flows with application to case 
 studies in Vietnam”, Ph.D thesis, University of Brescia, Italia. 
[5] Lê Thị Thu Hiền (2015). “Ứng dụng phương pháp số giải bài toán sóng gián đoạn trong 
 tính toán thủy lực khi đập bê tông vỡ”. Tạp chí khoa học kỹ thuật thủy lợi và môi trường, 
 50, 88-94. 
[6] P.L. Roe. (1981). “Approximate Riemann Solvers, parameter vectors and difference 
 schemes”. Journal of Computational Physics, 43, 357-372. 
[7] J. Sampson; A. Easton; M. Singh (2006). “Moving boundary shallow water flow above 
 bottom topography”. ANZIAM (EMAC2005), 47, C373-C387. 
[8] M.J. Castro, E.D. Fernandez Nieto, A.M. Ferreio, J.A. Garcia Rodriguez, C. Pares (2009). 
 “High order Extensions of Roe schemes for two dimensional Non Conservative Hyperbolic 
 Systems”. J. Sci. Comput, 39, 67 – 114. 
[9] Lê Thị Thu Hiền, Lê Thanh Hùng (2017). “Dự báo quá trình lưu lượng do vỡ đập bằng mô 
 hình toán nước nông hai chiều: Áp dụng cho hồ A Vương – Quảng Nam”. Tạp chí nông 
 nghiệp và Phát triển nông thôn, 7, 71-76. 
 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 113