Ả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Ệ 
 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 107 
Ả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 ĐỀ * 
Sử dụng mô hình toán mô phỏng các bài toán 
thủy lực như: sự lan truyền sóng lũ trên địa 
hình có độ dốc phức tạp; sóng gián đoạn đã 
được thực hiện rộng rãi trong thủy lợi. Với 
những bài toán thủy lực hai chiều hay ba 
chiều, kích thước ô lưới của miền tính toán 
luôn là yếu tố có ảnh hưởng lớn tới độ chính 
xác của kết quả tính bằng mô hình toán. Đặc 
biệt, trên những địa hình phức tạp, cao độ đáy 
biến đổi nhiều. 
Valiani và nnc (2002) đã dùng hai loại ô lưới 
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 
nghiên cứu, [1]. Wang (2011) chia miền tính 
toán cũng của hồ này thành lưới vuông với các 
kích cỡ khác nhau để tính sự lan truyền lũ và 
kết luận rằng, lưới 40m cho kết quả tốt nhất so 
Ngày nhận bài: 26/7/2017 
Ngày thông qua phản biện: 06/9/2017 
Ngày duyệt đăng: 26/9/2017 
với thực đo chứ không phải lưới có kích thước 
nhỏ hơn 30m, [2]. Lê Thanh Hùng (2017) lại 
khảo sát sự ảnh hưởng của 4 loại ô lưới tới kết 
quả tính sự lan truyền sóng lũ do tình huống 
vỡ đập Nậm Chiến (Sơn La), [3]. Về lý thuyết, 
chia lưới có kích thước càng nhỏ thì kết quả 
càng chính xác, tuy nhiên, khi số ô tính toán 
càng lớn thì thời gian tính sẽ càng nhiều, thậm 
chí các máy tính cá nhân thông dụng không xử 
lý được. Vì vậy, việc đánh giá ảnh hưởng của 
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 
thước ô lưới hợp lý là rất cần thiết. 
Trong nội dung bài báo này, các tác giả dùng 
một số ví dụ để đánh giá ảnh hưởng của kích 
thước ô lưới tới kết quả tính thủy lực như mực 
nước, lưu lượng dòng chảy lũ bằng chương 
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. 
2. PHƯƠNG PHÁP NGHIÊN CỨU 
Sử dụng chương trình tính 2D-FV dựa trên 
KHOA HỌC CÔNG NGHỆ 
TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 108 
phương pháp thể tích hữu hạn (FVM) loại 
Godunov giải hệ phương trình nước nông phi 
tuyến (NSWE) trên lưới Catersian đã được giới 
thiệu trong [3], [4], [5]. Điều kiện ổn định của 
phương pháp số được lấy theo Courant-
Fredrichs-Lewy (CFL). Với lưới Cartesian, bước 
thời gian t bị ràng buộc bởi phương trình (1): 
1~~~~
max
Δy
hgv
Δx
hgu
CrΔt , (1) 
trong đó: 
Cr là số Courant có giá trị trong khoảng 0 < Cr 
 1, các kết quả tính toán trình bày ở mục 3, hệ 
số Cr được lấy bằng 0,9. 
hvu ~,~,~ là các giá trị lưu tốc và độ sâu trung 
bình tính theo Roe, (1981), [6]. 
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 
thời gian chạy máy tính sẽ càng lâu. 
3. KẾT QUẢ TÍNH TOÁN 
3.1. Sóng dao động trong chảo parabol có 
ma sát 
Với miền tính toán có kích thước 10000m 
10000m, địa hình đáy mô tả bằng phương 
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 
quả của phương pháp số được đưa ra là mực 
nước tại các thời điểm t = 1000s và t = 6000s 
được so sánh với kết quả chính xác của 
nghiệm giải tích. 
 20 x/ahxzb , (2) 
trong đó các thông số h0 và a là các hằng số. 
Nghiệm giải tích được đưa ra lần đầu bởi 
Sampson và nnc [7] phụ thuộc vào hệ số ma 
sát đáy , trong trường hợp này số hạng ma sát 
uhτS f  . 
Nghiệm giải tích về mực nước dao động của 
bài toán này tuân theo biểu thức (3): 
 xstτBstBs
g
e
g
eBsts/τstsτAhx,tη
τt/τt
sin
2
cos
4
2cos42sin.
22
22
0 , (3) 
trong đó thông số B là hằng số và 222 /τps ; 208 /aghp , 
0
2
22
8 hg
eBaA
τt 
T=6000s - x=100m
0
10
20
30
-5000 -3000 -1000 1000 3000 5000x(m)
ca
o 
độ
(m
) bed
analytical
numerical
T=6000s - x = 50m
0
10
20
30
-5000 -3000 -1000 1000 3000 5000
x(m)
ca
o 
độ
(m
)
bed profile
analytical
numerical
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 để 
đánh giá độ chính xác của kết quả tính mực 
nước theo chương trình 2D-FV so với kết 
quả chính xác tương ứng với các kích thước 
lưới khác nhau: 
KHOA HỌC CÔNG NGHỆ 
 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 109 
 

n
i analyticalianalytical
n
i ielianalytical
XX
XX
E
1
2
,
1
2
,mod,1 , (4) 
trong đó Xanalytical,i là giá trị chính xác tính theo 
(3), còn Xmodel,i là nghiệm tính theo phương 
pháp số ở các vị trí khác nhau i. 
Từ kết quả ở bảng 1 cho thấy: Với 4 kích 
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 
= 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 
nước tương đối tốt. Vì vậy, trong những bài 
toán có địa hình phức tạp, nếu chỉ cần kết quả 
liên quan đến mực nước có thể dùng lưới thô 
để tính. 
Bảng 1: Chỉ số Nash tương ứng với các 
kích thước lưới khác nhau 
 x Số ô lưới E (%) (t=1000s) 
E (%) 
(t=6000s) 
25 
50 
100 
200 
400 400 
200 200 
100 100 
50 50 
99,47 
99,43 
98,64 
98,12 
98,65 
98,29 
98,04 
97,79 
3.2. Sóng gián đoạn trên bề mặt có mỏm núi 
Để đánh giá ảnh hưởng của kích cỡ lưới tính 
toán tới hai loại kết quả là mực nước, lưu lượng, 
ví dụ dưới đây do Castro và nnc (2009) giới 
thiệu được sử dụng, [8]. Trong miền tính toán 
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): 
   
  
222
222
5,015,1nêu 0
5,015,1nêu 12cos1)5,0(2cos
8
1
myx
myxyx
zb
 (5) 
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). 
222
222
1,0125,1nêu6,0
1,0125,1nêu1,1
),(
yx
yx
yx (6) 
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 -
0,2m3/s.m trong khi ảnh hưởng của hai kích 
thước lưới này tới mực nước tại cùng vị trí 
này không đáng kể. 
Hình 2: Mực nước ban đầu và địa hình đáy 
KHOA HỌC CÔNG NGHỆ 
TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 110 
t = 0.15s
0.3
0.4
0.5
0.6
0.7
0.8
0.4 0.6 0.8 1 1.2 1.4 1.6 1.8x(m)
el
ev
at
io
n(
m
)
bed
dx=0.005m
dx=0.0025m
dx=0.002m
t = 0.15s
-0.4
-0.2
0
0.2
0.4
0.4 0.6 0.8 1 1.2 1.4 1.6 1.8x(m)
di
sc
ha
rg
e(
m
3 /s
/m
)
dx=0.005m
dx=0.0025m
dx=0.002m
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 
Ví dụ này dùng để kiểm tra khả năng của mô 
hình trong việc xử lý các vấn đề khô ướt do 
địa hình biến đổi phức tạp. Địa hình được định 
nghĩa bằng hệ các phương trình sau để chia 
miền tính toán thành 3 lòng hồ. 
 ],,,,,,min[, 321 yxzyxzyxzyxz bbbb (7) 
trong đó: 
 10
2251225
250,
;
225
50
225
,
;
4001600
250,
22
3
22
2
22
1
yxyxz
yxyxz
yxyxz
b
b
b
Hai vật cản được mô phỏng bằng biểu thức: 
   
322
2
22
1
101020010,
5050
25080,,max,
yxneuyxz
yxyxzyxz
bB
bbB (8) 
Với điều kiện ban đầu, mực nước trong hồ là 
35m ở vị trí x -100m, phần còn lại đáy khô 
(hình 4). Tất cả các biên là biên đóng với hệ số 
nhám Manning bằng 0,033. Thời gian tính 
toán là 500s. Kết quả tính sự lan truyền sóng 
vỡ đập trên miền tính toán có kích thước 
1000m 400m được tính với 3 kích thước lưới 
khác nhau 2,5m; 5m và 10m. Hình 4 chỉ ra sự 
lan truyền sóng vỡ đập tại các thời điểm khác 
nhau. Sau 10s dòng chảy do vỡ đập lan truyền 
tới hồ chứa 2 và sau 40s dòng nước đã xuống 
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 
kiện địa hình rất phức tạp, mô hình toán vẫn 
có khả năng cho kết quả hợp lý khi không sinh 
ra nhiễu động. Bảng 2 lại chỉ ra các kết quả 
tính tương ứng với 3 kích thước lưới này, bao 
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 
lan truyền sóng lũ tới điểm nghiên cứu P(x = 
250m; y = 0m) (hình 4). Với kích thước lưới 
2,5m, thời gian chạy máy lên tới 40 phút so 
với 12 phút của lưới lớn gấp đôi 5,0m. Lưới 
8m và 10m thời gian này chỉ là 3 phút và 2 
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 
trăm sai số được xác định theo biểu thức (9): 
%100%
min
min 
x
xxi
err X
XX , (9) 
trong đó: 
X là đặc trưng thủy lực; X xmin là đặc trưng 
thủy lực tính với kích thước lưới nhỏ nhất. 
Bảng 2 cho thấy, sai số của lưu lượng lớn 
nhất tại đập ứng với các kích thước lưới 5,0m 
so với lưới mịn nhất 2,5m nhỏ hơn nhiều 
(-2,33%) so với lưới 8,0m và 10,0m (-6,92% 
và +11,07%). Vì vậy, trong trường hợp này 
lưới 5,0m có thể được dùng để tính lan truyền 
sóng lũ. 
KHOA HỌC CÔNG NGHỆ 
 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 111 
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 
 x (m) 
Thời gian 
chạy máy 
(phút) 
Qmax 
(m3/s) 
Sai số 
(%) 
Hmax 
(m) 
Sai số 
(%) 
Thời gian 
lũ đến 
(s) 
Sai số 
(%) 
(400 160) 
(200 80) 
(125 50) 
(100 40) 
40 
12 
3 
2 
9799,29 
9570,68 
9121,20 
10883,96 
-2,33 
-6,92 
+11,07 
20,33 
20,05 
19,79 
18,09 
-1,38 
-2,65 
-11,02 
18,82 
20,34 
19,95 
12,82 
+8,08 
+6,00 
-31,88 
3.4. Kịch bản vỡ đập vòm Nậm Chiến 
Đập vòm Nậm Chiến trên suối Chiến, huyện 
Mường La, Sơn La là đập vòm duy nhất ở Việt 
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 
945m, hạ lưu khô và coi đập vòm vỡ tức thời, 
hoàn toàn. Kết quả tính lưu lượng đỉnh lũ lớn 
nhất tại đập được thể hiện trên bảng 3. 
Hình 5: Bản đồ DEM lòng hồ Nậm Chiến 
P 
KHOA HỌC CÔNG NGHỆ 
TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 112 
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 
 x Qmax.103 (m3/s) (n=0,06) % sai số Qmax.103 (m3/s) (n=0,04) % sai số 
10 
20 
25 
30 
40 
50 
90 
114,41 
123,11 
121,78 
128,65 
119,60 
132,26 
106,52 
+7,6 
+6,4 
+12,4 
+4,5 
+15,6 
-6,9 
115,55 
124,58 
123,23 
128,73 
120,00 
133,83 
122,83 
+7,8 
+6,6 
+11,4 
+3,8 
+15,8 
+6,3 
Rõ ràng, lưới càng thô, số ô tính toán càng ít 
thì thời gian tính toán sẽ càng giảm. Tuy 
nhiên, từ những kết quả trên cho thấy, không 
phải lưới có kích thước càng nhỏ sẽ cho kết 
quả càng chính xác. Lưới 20m và 30m cho sai 
số là +7,6% và +12,4% tương ứng với độ 
nhám n = 0,06 và kết quả này là +7,8% và 
+11,44% khi n = 0,04. Tuy nhiên, khi x = 
40m, giá trị sai số này chỉ là +4,5% và +3,8% 
tương ứng với 2 độ nhám trên. Lưới 50m cho 
kết quả khác biệt nhất so với lưới 10m (hơn 
15% ứng với cả hai độ nhám), trong khi lưới 
thô nhất 90m lại cho kết quả thiên nhỏ. Vì vậy, 
ta có thể tính toán sự lan truyền sóng lũ cho cả 
miền tính toán lớn bao gồm cả hạ lưu hồ Nậm 
Chiến khi x = 40m. Hình 6 là kết quả tính sự 
lan truyền sóng lũ do vỡ đập vòm Nậm Chiến, 
tính với độ nhám n = 0,06 tại hai thời điểm t = 
500s và t = 1500s tương ứng với kích thước 
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 
Độ chính xác của kết quả tính theo phương 
pháp số luôn là mục tiêu hàng đầu khi mô 
phỏng các hiện tượng thủy lực theo phương 
pháp này. Ảnh hưởng của kích thước ô lưới tới 
kết quả tính mực nước chỉ ra rằng có sự khác 
biệt không nhiều giữa kích thước lưới mịn và 
lưới thô (mục 3.1, 3.2). Tuy nhiên, khi xét ảnh 
hưởng này tới lưu lượng thì có sự khác biệt 
đáng kể, đặc biệt khi địa hình thay đổi phức 
KHOA HỌC CÔNG NGHỆ 
 TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ THỦY LỢI SỐ 40 - 2017 113 
tạp như ví dụ trong mục 3.2 và 3.3. Áp dụng 
chương trình tính 2D-FV cho một trường hợp 
thực tế là vỡ đập hồ Nậm Chiến, Sơn La tương 
ứng với 7 kích thước lưới khác nhau chỉ ra 
rằng: Sự ảnh hưởng của kích thước lưới tới lưu 
lượng lớn nhất tại đập là rất rõ ràng. Lưới 40m 
cho kết quả tốt hơn cả lưới nhỏ hơn là 30m và 
20m và cả lưới lớn hơn 50m, 90m khi so sánh 
với kết quả tính theo kích thước lưới mịn nhất 
10m. Việc tìm ra kích thước lưới hợp lý sẽ đưa 
ra được kêt quả có độ sai số cho phép, đồng 
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.