Mô hình tính toán dòng chảy và vận tải bùn cát ba chiều trong sông và kênh hở

Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Bài Nghiên cứu
Trường Đại học Bách khoa,
ĐHQG-HCM
Liên hệ
Lê Song Giang, Trường Đại học Bách khoa,
ĐHQG-HCM
Email: lsgiang@yahoo.com
Lịch sử
 Ngày nhận: 06-11-2018
 Ngày chấp nhận: 30-5-2019
 Ngày đăng: 22-6-2019
DOI :
https://doi.org/10.32508/stdjsee.v3i1.508
Bản quyền
© ĐHQG Tp.HCM. Đây là bài báo công bố
mở được phát hành theo các điều khoản của
the Creative Commons Attribution 4.0
International license.
Mô hình tính toán dòng chảy và vận tải bùn cát ba chiều trong
sông và kênh hở
Lê Song Giang*, Trần Thị Mỹ Hồng
TÓM TẮT
Mô hình toán là một công cụ hữu ích trong nghiên cứu dòng chảy và vận tải bùn cát, bồi xói lòng
dẫn và được xây dựng dựa trên việc giải các phương trình vi phân cơ bản. Mô hình toán có nhiều
cấp độ khác nhau và mô hình 3 chiều là cấp độ cao nhất, cho phép mô phỏng chi tiết dòng chảy
và quá trình vận tải bùn cát trong không gian 3 chiều. Bài báo này trình bày một sơ đồ tính toán
dòng chảy và vận tải bùn cát ba chiều trong kênh hở. Mực nước và vận tốc dòng chảy được giải
từ phương trình chuyển động ba chiều với giả thuyết thủy tĩnh. Nồng độ bùn cát lơ lửng, bùn cát
đáy và diễn biến đáy được giải từ các phương trình vận tải. Các phương trình vi phân cơ bản được
viết trong hệ tọa độ biến đổi ``sigma'' và được giải theo phương pháp thể tích hữu hạn trên lưới
phi cấu trúc phần tử tứ giác. Đối với dòng chảy điều kiện biên mực nước hoặc lưu lượng sẽ được
áp đặt trên các biên hở. Đối với bùn cát lơ lửng nồng độ bùn cát trên biên được áp đặt trong pha
chảy vào và điều kiện thoát tự do được sử dụng trong pha chảy ra. Sơ đồ tính toán được kiểm tra
với bài toán vận tải bùn cát trong kênh cong được nghiên cứu bằng thực nghiệm bởi Odgaard và
Bergs. Kết quả tính toán khá phù hợp với số liệu đo. Nhằm kiểm tra khả năng ứng dụng thực tế, sơ
đồ tính toán mới cũng được kiểm tra với bài toán vận tải bùn cát trên sông Đồng Nai tại khu vực
Cù lao Phố. Để giải quyết vấn đề điều kiện biên thủy lực của bài toán này, mô hình khu vực Cù Lao
Phố được tích hợp vào mô hình hệ thống sông Sài Gòn – Đồng Nai. Kết quả tính diễn biến lòng
dẫn khu vực cù lao Phố trên sông Đồng Nai cũng cho thấy phương pháp tính toán này cho kết quả
phù hợp với quy luật và hoàn toàn có thể sử dụng trong nghiên cứu thực tế.
Từ khoá: mô hình 3D, vận tải bùn cát, lưới phi cấu trúc, tọa độ ``sigma''
GIỚI THIỆU
Tính toán dòng chảy và vận tải bùn cát là một trong
các bài toán phổ biến trong thủy lực. Đối với các
trường hợp đơn giản, người ta có thể dùng mô hình
1D hoặc 2D. Tuy nhiên đối với các bài toán phức tạp,
nơi dòng chảy có cấu trúc 3 chiều với các dòng thứ
cấp xuất hiện rõ rệt thì chỉ có mô hình 3D mới đủ
độ tin cậy. Gần đây, một số mô hình 3D đã được
công bố. Van Rijn đã giải phương trình vận tải bùn
cát lơ lửng 3D trong khi dòng chảy 3D nhận được
từ lời giải 2D cùng với giả thiết phân bố theo quy
luật logarithm trên phương thẳng đứng1. Lin và Fan-
coner giải phương trìnhNavier-Stokes trung bình hóa
Reynolds và phương trình vận tải bùn cát lơ lửng bằng
phương pháp sai phân hữu hạn trên lưới cố định theo
phương thẳng đứng2. Wu và các tác giả đã xây dựng
một phương pháp tính dựa trên việc giải bằng phương
pháp thể tích hữu hạn phương trình Navier-Stokes
bình hóa Reynolds và phương trình vận tải bùn cát
lơ lửng trung3. Tuy nhiên lưới tính của các thuật giải
này chưa đủ mềm dẻo để mô phỏng tốt các bài toán
có hình học phức tạp. Ngoài ra thuật giải cũng không
hoàn toàn thích hợp để xây dựng mô hình tích hợp.
Trong bài báo này, một phương pháp tính toán dòng
chảy và vận tải bùn cát 3D khác được giới thiệu. Dòng
chảy được giải từ phương trình 3 chiều với giả thiết
thủy tĩnh. Các phương trình được giải theo phương
pháp thể tích hữu hạn trên lưới phi cấu trúc và trục
thẳng đứng biến đổi sigma. Phương pháp đã được áp
dụng tính toán thử nghiệm với bài toán kênh cong
180o được nghiên cứu bằng thực nghiệmbởiOdgaard
và Bergs4 và bài toán bồi xói khu vực Cù lao Phố trên
sông Đồng Nai.
MÔHÌNH TOÁN
Phương trình cơ bản
Mực nước và vận tốc của dòng chảy được giải từ
phương trình chuyển động ba chiều với giả thiết thủy
tĩnh5. Trong hệ tọa độ biến đổi “sigma”, phương trình
được viết :
¶D
¶ t
+Ñs (q)+
¶qs
¶s
= 0 (1)
¶q
¶ t
+Ñs [qUDAHÑsU]
+
¶
¶s
[
qsU AVD
¶U
¶s
]
=gDÑsh+DF
(2)
Trích dẫn bài báo này: Giang L S, Mỹ Hồng T T. Mô hình tính toán dòng chảy và vận tải bùn cát ba
chiều trong sông và kênh hở. Sci. Tech. Dev. J. - Sci. Earth Environ.; 3(1):23-36.
23
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Nồngđộbùn cát lơ lửngđược giải từ phương trình vận
tải1,3. Cũng trong hệ tọa độ biến đổi “sigma”, phương
trình được viết:
¶qC
¶ t
+Ñs [qCHDHÑsC]
+
¶
¶s
[
(qs ws0)C DVD
¶C
¶s
]
= 0
(3)
Diễn biến đáy được giải từ phương trình 3 :
(1P)rC ¶ zb¶ t +
¶ (bcb)
¶ t
+Ñqb = (DbEb) (4)
Trong các phương trình tên: h – mực nước; U =[
ux;uy
]T– hai thành phần vận tốc của dòng chảy trên
phương ngang; w – thành phần vận tốc trên phương
thẳng đứng “sigma”; q=
[
qx;qy
]T
=DUvà qs =Dw
D – độ sâu; Ñs – toán tử vi phân trên mặt phẳng
“sigma”; F– vector ngoại lực; AH và AV – độ nhớt rối;
C – nồng độ bùn cát lơ lửng trong cột nước; qC =DC
ws0 – vận tốc lắng; zb – cao độ đáy; rC – khối lượng
riêng hạt bùn cát; P – hệ số rỗng; cb – nồng độ bùn cát
trong tầng đáy; b – bề dày tầng đáy; qb =
[
qbx;qby
]T
vector lưu lượng đơn vị của bùn cát đáy; Db và Eb –
suất lắng và xói của bùn cát tại tầng đáy (ở cách đáy
một khoảng là b); DH và DV – hệ số khuếch tán rối.
Độ nhớt rối ngang, AH , được tính theo Smagorin-
sky6:
AH =C∆x∆y
[(
¶u
¶x
)2
+
(
¶v
¶y
)2
+
1
2
(
¶u
¶y
+
¶v
¶x
)2]0;5 (5)
Hệ số C nằm trong khoảng 0,01 – 0,5. Độ nhớt
rối theo phương đứng, AV , được tính theo mô hình
Prandtl-Kolmogorov (1942) 7,8:
AV =C
′
mL
p
k (6)
Trong đó C′m - hằng số mô hình, xác định trong quá
trình hiệu chỉnh mô hình; L – chiều dài xáo trộn; k –
động năng rối. Các thông số này được lấy theo Davies
và Gerritsen 9 như sau:
k =
1pcm
[(
ub
)2
(h)+(us)2 (1h)
]
(7)
L= kD(1h)ph (8)
Với cm = 0;09 – hằng số mô hình; k = 0;4– hằng
số Karman; us- vận tốc ma sát trên mặt nước; ub-
dạng sửa chữa của vận tốc ma sát đáy ub trong đó
ub = max(u;ub)u- giá trị trung bình của các vận
tốc ma sát tính từ vận tốc tại các mắt lưới trên đường
thủy trực nếu biên dạng vận tốc là logarit; và h =
(h z)=D- độ sâu tương đối.
Hệ số khuếch tán rối được lấy bằng với độ nhớt rối.
Đối với bùn cát rời, suất lắng và xói bùn cát ở độ cao
b, Db và Eb, được tính3 :
DbEb = ws0
(
cb cb;e
)
(9)
Trong đó cb và cb;e – nồng độ bùn cát và nồng độ
bùn cát bão hoà ở mặt phân cách lớp đáy. Theo Van
Rijn1,10,11, lưu lượng bùn cát đáy, q b, và nồng độ bùn
cát bão hoà, c b;e, được tính:
qb = 0:053rC [(s1)g]0:5 d1:550 D0:3 T 2:1 (10)
cb;e = 0.015rC
d50
b
T 1.5
D0.3
(11)
Với:
T = (tb tcr)=tcr (12)
D = d50
[
(s1) g
v2
]1=3
(13)
Trongđó: d50 –đường kínhhạt 50%; s – tỷ trọnghạt; t
và tcr – ứng suất tiếp đáy và ứng suất tiếp đáy ngưỡng
(xác định từ đồ thị Shield).
Đối với bùn cát kết dính, suất xói, E b, và lắng, Db,
được tính theo Hayter và Mehta12 và Krone13:
Eb = e
(
tb
te
1
)
(14)
Db = ws0C
(
1 tb
td
)
(15)
Trong đó: e - hệ số xói; tb - ứng suất tiếp đáy; te -
ứng suất tiếp đáy ngưỡng xói; td- ứng suất tiếp đáy
ngưỡng bồi.
Điều kiện biên
Trên biên hở diễn biến lưu lượng hoặc mực nước sẽ
được áp đặt. Đối với nồng độ bùn cát lơ lửng, trong
pha chảy vào, nồng độ bùn cát được áp đặt. Còn trong
pha chảy ra, điều kiện thoát tự do được sử dụng1:
¶C=¶n= 0 (16)
Trên biên kín, điều kiện không thẩm thấu được sử
dụng.
Trên mặt thoáng (s = 0), các điều kiện biên sau được
sử dụng1,5:
w = 0 (17)
AV
D
[
¶u
¶s
;
¶v
¶s
]
=(t0x;t0y)=r0 (18)
ws0C+
DV
D
¶C
¶s
= 0 (19)
24
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Còn tại đáy (s =1)1,5:
w = 0 (20)
AV
D
[
¶u
¶s
;
¶v
¶s
]
=CD
(
u2+ v2
)0.5
(u;v) (21)

(
ws0C+
DV
D
¶C
¶s
)
= EbDb (22)
Trong đó CD là hệ số ma sát đáy và được tính 5:
CD =
[
k
ln(∆=z0)
]2
(23)
Với ∆ – khoảng cách từ mắt lưới đầu tiên tới đáy; và
z0 – thông số nhám.
Phương pháp giải
Các phương trình (1) - ( 3 ) được chúng tôi giải bằng
phương pháp thể tích hữuhạn theo sơ đồ được cải tiến
từ sơ đồ được giới thiệu trong tài liệu14. Trong sơ đồ
mới này, các thành phần của lưu lượng q và nồng độ
bùn cát lơ lửng C được giải ẩn theo phương”sigma” để
gia tăng tính ổn định của lời giải.
Hình 1 giới thiệu lưới tính của mô hình. Mực nước
được tính tại các nút lưới. Lưu lượng q được tính tại
các cạnh của các phần tử, lưu lượng qs được tính tại
các nút lưới ở khoảng giữa các lớp còn nồng độ bùn
cát lơ lửng được tính tại các nút lưới ở ngay từng lớp.
Thứ tự giải các phương trình như sau:
- Bước 1: Giải phương trình (2) tính lưu lượng q ở
thời điểm n+1/2.
- Bước 2: Giải phương trình (1) tính mực nước h ở
thời điểm n+1
- Bước 3: Tính lưu lượng qs ở thời điểm n+1/2.
- Bước 4: Giải phương trình (3) tính nồng độ bùn cát
lơ lửng C ở thời điểm n+1
- Bước 5: Giải phương trình (4) tính cao độ đáy zb ở
thời điểm n+1
a) Bước 1
Phương trình (2) được tích phân trên thể tích kiểm
soát của cạnh ict (Hình 1a) ở lớp k:∫
dsk
∫
S
¶q
¶ t
dSds +
∫
sk
∫
S
Ñs [qUDAHÑsU]
dSds +
∫
dsk
∫
S
¶
¶s
[
qsU AVD
¶U
¶s
]
dSds =
∫
dsk
∫
S
(gDÑsh+DF)dSds
(24)
Sau khi sai phân số hạng đạo hàm theo thời gian, (24)
đưa tới phương trình:(
qn+1=2k q
n1=2
k
)
∆Vk
∆t
+Fluxk
(
Un1=2
)
+
(
Jn+1=2k+1  J
n+1=2
k
)
= rnk∆Vk
(25)
Trong đó
∆V k = dsk .S (26)
Jk+1 =
S
D
[
qsq AVD
¶q
¶s
]
k+1=2
(27)
Fluxk (U) = ds k
H
L
[
qnUDAH
¶U
¶n
]
k
dl (28)
rk = D(gÑsh+F)k (29)
Các số hạng Fux k(U) và rk sẽ được ước tính bằng các
phép tích phân số và sai phân. Riêng số hạng J k+1 sẽ
được nội suy và sai phân thành:
Jk+1 =A′′Tk+1qk+1+A
′′
Pk+1qk (30)
Trong đóA′′Tk+1àA
′′
Pk+1được xác định tùy theo sơ đồ nội
suy. Riêng tại đáy và mặt thoáng, điều kiện biên sẽ
được sử dụng để tính J. Thay (30) vào (25), ta được
phương trình cuối cùng:
ASkqn+1=2k1 +APkq
n+1=2
k ANkq
n+1=2
k+1 = Srck (31)
Trong đó:
APk =
∆V k
∆t
+
(
A
′′
Pk+1 +A
′′
Tk
)
(32)
ATk = A
′′
Tk+1 (33)
ASk = A
′′
Pk (34)
Srck =
∆Vk
∆t
qn1=2k Fluxk
(
Un1=2
)
+∆Vk .rnk (35)
Trên mỗi thủy trực ta sẽ thiết lập được một hệ nz
phương trình (31) chứa nz ẩn số mà sau khi giải ta
sẽ có nz nghiệm qk (k=1,nz) tại cạnh ict trên các lớp.
b) Bước 2
Để tính mực nước h , phương trình liên tục (1) được
tích phân trên phương s trên toàn bộ chiều sâu. Kết
quả ta thu được:
¶D
¶ t
+å
k
(ds kÑsqk) = 0 (36)
(36) sau đó được tích phân trên diện tích kiểm soát
của nút C (Hình 1a), từ đó ta được:
¶Dn+1=2C
¶ t
=1
Såk
(
dskå
j
qn+1=2nk j L j
)
(37)
Từ (37), độ sâu và mực nước tại nút C được tính:
Dn+1C = D
n
C+
¶Dn+1=2C
¶ t
∆t (38)
hn+1C = zb+D
n+1
C (39)
c) Bước 3
25
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 1: Lưới tính trên mặt bằng và các lớp lưới tính theo phương thẳng đứng. a. Trên mặt bằng; b. Trên
phương thẳng đứng
Hình 2: Lưới tínhmô hình 3D.
Để tính lưu lượng đơn vị trên phương s , phương
trình liên tục (1) được tích phân trên thể tích kiểm
soát ở nút C tại lớp k:
∫
S
∫
ds k
¶D
¶ t
dsdS+
∫
S
∫
ds k
ÑsqdsdS+
∫
S
∫
ds k
¶qs
¶s
dsdS= 0
(40)
Các tích phân được ước tính bằng các tích phân số và
sau đó sai phân theo thời gian ta được lời giải cho qs
ở nút C trên các lớp ở thời điểm n+1/2:
qn+1=2sk+1 = q
n+1=2
sk dsk
¶Dn+1=2C
¶ t
 dsk
S
å j q
n+1=2
nk j L j
(41)
Lưu ý là tại đáy qs n+1=21 = 0
d) Bước 4
26
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 3: Trường vận tốc trênmặt thoáng.
Phương trình (3) được tích phân trên thể tích kiểm
soát tại nút C ở lớp k:∫
S
∫
dsk
¶qc
¶ t
dsdS+
∫
S
∫
dsk
Ñs [qCHDHÑsC]
dsds++
∫
S
∫
dsk
¶
¶s
[
qsC DVD
¶C
¶s
]
dsd
S=
∫
S
∫
sk
DSCdsdS
(42)
Sau khi sai phân số hạng đạo hàm theo thời gian (42)
đưa tới phương trình:
∆Vk
∆t
(
qCn+1k qCnk
)
+Fluxk (Cn)
+
(
Jn+1k+1  Jn+1k
)
= ∆Vk .Dn+1
(
r sCn+1k
) (43)
Trong đó:
∆V k = ds k .S (44)
Jn+1k+1 = S
[
qsC DVD
¶C
¶s
]n+1
k+1=2
(45)
Fluxk (Cn) = ds k∫
S
Ñs
[
qn+1=2CnHDHÑsCn
]
k
dS (46)
Số hạng Fuxk(Cn) sẽ được ước tính bằng các phép tích
phân số và sai phân còn số hạng J k+1 sẽ được nội suy
và sai phân thành:
Jk+1 =A′′Tk+1Ck+1+A
′′
Pk+1Ck (47)
Trong đóA′′Tk+1àA
′′
Pk+1được xác định tùy theo sơ đồ nội
suy. Tại đáy và mặt thoáng, điều kiện biên sẽ được sử
dụng để tính J. Thay (50) vào (46), ta được phương
trình cuối cùng:
ASkCn+1k1 +APkCn+1k ANkCn+1k+1 = Srck (48)
Với
APk =
(
1
∆t
+ s
)
∆V kDn+1+
(
A
′′
Pk+1 +A
′′
Tk
)
(49)
Ark = A
′′
Tk+1 (50)
Ask = A
′′
pk (51)
Srck = ∆V k
(
Dn
∆t
Cnk + r.Dn+1
)
Fluxk (Cn) (52)
Trên mỗi thủy trực ta sẽ thiết lập được một hệ nz
phương trình (48) chứa nz ẩn số mà sau khi giải ta
sẽ có nz nghiệm qC (k=1,nz) ở nút C trên các lớp ở
thời điểm n+1.
e) Bước 5
Để tính biến đổi đáy, phương trình (4) được tích phân
trên diện tích kiểm soát của nút C (Hình 1a):
(1P)rC
∫
S
¶ zb
¶ t
dS+
∫
S
¶ (bcb)
¶ t
dS
+
∫
S
ÑqbdS=
∫
S
(DbEb)dS
(53)
27
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 4: Vận tốc trênmặt phẳngmặt cắt ở các góc cong khác nhau.
Sau khi sai phân số hạng đạo hàm theo thời gian, (53)
đưa tới phương trình:
(1P)rC
zn+1b  znb
∆t
+DC
+Flux(qb) = (DbEb)S
(54)
Trong đó
DC= S
∆t
[
(bcb)n+1 (bcb)n
]
(55)
Flux(qb) =
H
L qbndl (56)
Từ (54) ta có cao độ đáy ở nút C ở thời điểm n+1:
zn+1b = z
n
b+
∆t
(1P)rC
[(DbEb)SDCFlux(qb)]
(57)
TÍNH TOÁN THỬNGHIỆM
Vận tải bùn cát trong kênh cong 180o
Odgaard và Bergs4 có một nghiên cứu bằng thực
nghiệm dòng chảy và vận tải bùn cát trong một kênh
cong. Bài toán này cũng đã được Wu và các tác giả
tính toán bằng mô hình toán 3D 3. Kênh gồm một
đoạn cong 180º, bán kính trung bình 13,11m và 2
đoạn kênh thẳng dài 20mdẫn vào và ra (Hình2 ). Mặt
cắt ngang kênh hình thang, đáy rộng 1,44m và mặt
thoáng rộng 2,44m. Đáy kênh được phủ một lớp cát
dày 23cm với đường kính trung bình 0,3 mm. Nước
và cát trôi theo dòng chảy được bơm đưa trở lại đầu
kênh trong một hệ thống tuần hoàn khép kín với lưu
lượng 0,153m3 /s. Sau khi đáy kênh biến đổi đạt tới
trạng thái ổn định, địa hình đáy kênh đã được đo đạc.
Để đánh giá phương pháp phát triển trong bài báomô
hình toán 3D cho kênh củaOdgaard và Bergs đã được
28
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 5: Vận tốc vuông góc với mặt cắt ở các góc cong khác nhau.
xây dựng. Trên mặt bằng, đoạn cong 180º được chia
thành 18016 phần tử còn mỗi đoạn kênh dẫn vào
và ra được chia thành 8016 phần tử (Hình 2 ). Mô
hình có 11 lớp theo phương thẳng đứng. Địa hình
đáy của mô hình được lấy theo hình dạng của kênh
cứng hình thang và điều kiện ban đầu là kênh đã bị
bồi 23cm tương ứng với lớp cát phủ ban đầu trong thí
nghiệm của Odgaard và Bergs.
Đầu vào của kênh được áp đặt điều kiện biên lưu
lượngQ=0,153m3 /s và nồng độ bùn cát lơ lửng tương
ứng với lưu lượng bùn cát ra ở mặt cắt cuối kênh
3,7g/m/min như đã được ghi nhận trong thí nghiệm
của Odgaard và Bergs. Điều kiện biên mực nước ở
cuối kênh được áp đặt sao cho độ sâu trong kênh
khoảng 15cm như trong thí nghiệm.
Tính toán cho thấy sau khoảng 9 giờ kể từ thời điểm
ban đầu địa hình đáy kênh gần như không còn thay
đổi. Hình 3 giới thiệu trường vận tốc trên mặt nước
sau khi đáy kênh đã ổn định còn Hình 4 và Hình 5
giới thiệu trường vận tốc tại một số mặt cắt kênh ở
các vị trí khác nhau. Các hình cho thấy xoáy thứ cấp
trong kênh đã được tái hiệnmột cách rõ ràng. Hình 6
giới thiệu kết quả tính diễn biến đáy kênh sau khi đã
ổn định. Đường liền là đáy kênh sau 9 giờ còn đường
đứt đoạn là đáy kênh sau 11 giờ. Kết quả tính toán
khá phù hợp với số liệu thí nghiệm của Odgaard và
Bergs.
Vận tải bùn cát khu vực Cù Lao Phố trên
sông Đồng Nai
Sông Đồng Nai khu vực Cù lao Phố từng có một dự
án san lấp lấn sông. Mô hình 2D đã được sử dụng
để nghiên cứu dự báo khả năng diễn biến bồi xói khu
vực15,16. Tuy nhiên trong bài toán này mô hình 2D là
không đủ tin cậy do không có khả năng tính toán mô
phỏng các dòng thứ cấp phát sinh do địa hình lòng
sông phức tạp, đặc biệt là ở khu vực đầu Cù lao Phố.
29
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 6: Địa hình đáy kênh tính toán (đường liền và đứt nét) và thí nghiệm (chấm tròn).
Đây chính là bài toán dành cho mô hình 3D.
Để đánh giá phương pháp tính 3Dphát triển trong bài
báo trong bài toán thực tế, một tính toán thử nghiệm
dòng chảy và quá trình bùn cát cho đoạn Cù lao Phố
đã được thực hiện.
Hình 7 giới thiệu lưới tính của mô hình. Trên mặt
bằng, đoạn sông được chia thành 7261 phần tử. Số
phần tử theo chiều ngang sông trên nhánh chính là
16, còn trên nhánh phụ là 8. Kích thước của các phần
tử khoảng 16 - 25 m. Mô hình có 5 lớp. Lưới tính
này là khámịn, đủ đểmô phỏng chi tiết cấu trúc dòng
chảy khu vực nghiên cứu. Địa hình đáymô hình được
thamkhảo từ số liệu củaViệnKhoa họcThủy lợimiền
Nam khảo sát năm 2008 16. Địa hình đáy củamô hình
cũng được giới thiệu trênHình 8.
Để giải quyết vấn đề điều kiện biên thủy lực, mô hình
được tích hợp vào mô hình hệ thống sông Đồng Nai
xây dựng bằng phần mềm F2817, kế thừa từ kết quả
nghiên cứu 18. Mô hình hệ thống sông Đồng Nai này
đã được hiệu chỉnh kỹ với khoảng 10 bộ số liệu thực
đo.
Điều kiện biên bùn cát cho mô hình 3D được áp đặt
một cách độc lập ngay tại các mặt cắt biên. Theo
kết quả khảo sát của Viện Khoa học Thủy lợi Miền
Nam16, nồng độ bùn cát lơ lửng tại khu vực này vào
khoảng 30 - 40 mg/l. Giá trị này đã được dùng làm
điều kiện biên.
Vật liệu đáy sông khu vực làm mô hình không đồng
nhất. Giữa dòng là cát mịn tới thô trong khi ở gần bờ
là bùn. Tuy nhiên do hạn chế của phương pháp, mô
hình chỉ có khả năng tính với một loại vật liệu đáy. Vì
vậy trong nghiên cứu này vật liệu đáy được xem là cát
mịn với d50 =0,35mm tương ứng với đường kính hạt
trung bình ở khu vực.
Bước thời gian tính được lấy khá nhỏ, ∆t = 0,6s, để
đảm bảo điều kiện ổn định. Mặc dù bước thời gian
tính nhỏ nhưng tốc độ tính toán vẫn khá tốt. Trên
máy PC Core i5-3,2GHz, một giờ máy tính tính mô
30
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 7: Lưới tínhmô hình 3D khu vực Cù Lao Phố.
Hình 8: Địa hình đoạn sông Đồng Nai khu vực Cù lao Phố (đơn vị thangmàu: mét).
phỏng được khoảng 8 giờ.
Thông thường, mô hình toán trước khi sử dụng sẽ
phải được hiệu chỉnh. Đặc biệt nếu nghiên cứu là
nhằmđưa ra bức tranh chi tiết và tin cậy về vận tải bùn
cát ở Cù lao Phố thì hiệu chỉnh mô hình sẽ là điều bắt
buộc. Tuy nhiênmôhình 3DCù lao Phố này sẽ không
được hiệu chỉnh mà các thông số của mô hình được
xác định theo kinh nghiệm. Trong tính toán thông số
nhám được lấy z0 = 0,01mm, đường kính hạt trung
bình d50 =0,35mm, tỷ trọng hạt d = 0;65, độ thô thủy
lực ws0 =5,19cm/s và độ rỗng P=0,2. Điều này là có
thể chấp nhận vì các khó khăn về số liệu đo đạc và vì
mục tiêu của bài báo là phát triển một phương pháp
tính toán vận tải bùn cát trong điều kiện lòng dẫn
phức tạp và mô hình Cù lao Phố chỉ là minh chứng
về khả năng ứng dụng vào thực tế của phương pháp.
Tính toán mô phỏng đã được thực hiện cho 1 tháng
mùa lũ năm 2008. Hình 9 giới thiệu kết quả tính
31
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 9: Trường vận tốc trênmặt nước khu vực Cù lao Phố.
trường vận tốc trên mặt nước đặc trưng tại khu vực
cù lao Phố lúc triều lên và lúc triều xuống. Phân bố
vận tốc tại một số mặt cắt ngang lúc triều xuống vào
thời gian lũ cũng được giới thiệu trênHình 10 và 11.
Vị trí các mặt cắt được giới thiệu trênHình 12.
Độ bồi xói khu vực đầu Cù lao Phố được giới thiệu
trên Hình 13. Để thấy rõ hơn bức tranh bồi xói tại
khu vực, ta có thể xem một vài lát cắt. Hình 14 giới
thiệu phân bố vận tốc và biến đổi cao độ đáy tại các
lát cắt 5-5 và 6-6. Vị trí các lát cắt này được chỉ ra trên
Hình 13.
Lát cắt 5 – 5 cho thấy rõ các mô ngầm bị dòng chảy
bào mòn và di chuyển dần về hạ lưu. Trong khi đó lát
cắt 6 - 6 lại cho thấy cấu trúc 3D của dòng chảy ở đầu
cù lao. Có một xoáy ngang xuất hiện ở đầu cù lao và
xoáy này mang bùn cát từ chân cù lao ra ngoài. Xoáy
này là một trong các nguyên nhân chính làm cho đầu
cù lao bị xói mạnh. Khả năngmô phỏng các dòng thứ
cấp này cũng chính là ưu thế của mô hình 3D trước
mô hình 2D.
KẾT LUẬN
Bài báo đã trình bàymột phương pháp tính toán dòng
chảy và vận tải bùn cát, biến hình lòng dẫn 3D. Các
phương trình vi phân cơ bản được giải bằng phương
pháp thể tích hữu hạn trên lưới tính phi cấu trúc phần
tử tứ giác. Kết quả tính toán kiểm tra với số liệu thực
nghiệm cho thấy phương pháp cho lời giải hợp lý và
khá chính xác. Bên cạnh đó việc tính toán trường hợp
diễn biến lòng dẫn khu vực Cù lao Phố cũng cho thấy
phương pháp cho kết quả phù hợp quy luật và hoàn
toàn có thể sử dụng trong nghiên cứu thực tế.
32
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 10: Vận tốc trênmặt phẳngmặt cắt vàomùa lũ tại các vị trí khác nhau.
Hình 11: Vận tốc vuông góc với mặt cắt vàomùa lũ tại các vị trí khác nhau.
33
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 12: Vị trí các mặt cắt.
Hình 13: Trường vận tốc và độ bồi xói tại khu vực đầu cù lao Phố (đơn vị thangmàu: mét).
34
Tạp chí Phát triển Khoa học và Công nghệ – Khoa học Trái Đất và Môi trường, 3(1):23- 36
Hình 14: Phân bố vận tốc và biến đối đáy tại các lát cắt (đường màu đen nhạt: đáy sông ban đầu; đường
màu đen đậm: đáy sông saumột tháng).
XUNGĐỘT LỢI ÍCH
Nhóm tác giả cam đoan rằng không có xung đột lợi
ích trong công bố bài báo “Mô hình tính toán dòng
chảy và vận tải bùn cát ba chiều trong sông và kênh
hở”.
ĐÓNGGÓP CỦA TÁC GIẢ
Tác giả Lê Song Giang xây dựng thuật toán của mô
hình và viết chương trình tính toán.
Tác giả TrầnThịMỹHồng tính toán các bài toánmẫu.
TÀI LIỆU THAMKHẢO
1. Rijn LCV. Principles of sediment transport in rivers, estuaries
and coastal seas. Aqua Publications. 1993;p. 690–690.
2. Lin B, Falconer RA. Falconer Three-dimensional layer inte-
gratedmodelling of estuarine flowswith flooding and drying.
Estuar Coast Shelf Sci. 1997;44:737–751.
3. Wu W, Rodi W, Wenka T. 3D numerical modeling of flow
and sediment transport in open channel. J Hydraul Eng.
2000;126(1):4–15. Available from: 10.1061/(ASCE)0733-
429(2000)126:1(4).
4. Odgaard AJ, Bergs MA. Flow processes in a curved alluvial
channel. Water Resour. 1988;Res., 24(1):45–56. Res. Available
from: 10.1029/WR024i001p00045.
5. Ahsan AKMQ, Blumberg AF. Three-dimensional hydrother-
mal model of Onondaga Lake, New York. J Hydralic Eng.
1999;125(9):912–923. Available from: 10.1061/(ASCE)0733-
9429(1999)125:9(912).
6. Smagorinsky J. General circulation experiments with the
primitive equations. I: TheBasic Experiment. MonthlyWeather
Rev. 1963;91:99–164.
7. Kolmogorov AN. Equations of turbulent motion of an incom-
pressible fluid. Equations of turbulent motion of an incom-
pressible fluid Izvestia Academy of Sciences, USSR. 1942;6(2
and 2):56–58. Physics.
8. Prandtl L. Uber ein neues formelsystem fur die ausgebildete
turbulenz. Nacr Akad Wiss Gottingen, Math-Phys. 1945;p. 6–
19.
9. Davies AM, Gerritsen H. An intercomparision of three-
dimensional tidal hydrodynamic models of the Irish sea. Tel-
lus. 1994;46:200–221.
10. Rijn JV. Hydralic Eng., ASCE, Vol. 110 (11) , p. 1613-1641; 1984b.
11. Rijn JV, Eng. Hydralic Eng., ASCE, Vol. 110 (10), p. 1431-1456.
ASCE; 1984a.
12. Hayter EJ, Mehta AJ. Modelling cohesive sediments trans-
port in estuarine waters. Applied Mathematical Modelling.
1986;51:765–778.
13. Krone RB. Flume studies of the transport of sediment in es-
tuarine shoaling processes. In: Final report to San Fransico
District U.S. Army Corps of Engineers. Washington D.C; 1962.
14. Hồng TTM. Lê Song Giang. Một sơ đồ tính toán dòng chảy
sông, biển 3 chiều trên lưới tính phi cấu trúc. Tuyển tập Công
trình Hội nghị khoa học Cơ học Thủy khí toàn quốc lần thứ 20.
2017; tr. 335-343; 2017.
15. Cty cổ phần ĐT-KT-XD Toàn Thịnh Phát. Báo cáo đánh giá tác
độngmôi trường dự án cải tạo cảnh quan và phát triển đô thị
ven sông Đồng Nai, quy mô 8,4ha tại phường Quyết Thắng,
Tp. Biên Hòa, tỉnh Đồng Nai.; 2014.
16. Viện Khoa học Thủy lợi miền Nam. Đánh giá tác động dòng
chảy sôngĐồngNai đoạn từ cầuHóaAnđến cầuGhềnh thuộc
thành phố Biên Hòa. Tp. Hồ Chí Minh, tháng 10 năm. ; 2009.
17. Giang LS. Development of an integrated software for calcu-
lation of urban flood flow. Report B2007-20-13TĐ. VNU-HCM;
2011.
18. Giang LS. Nghiên cứu đề xuất lựa chọn chiến lược quản lý
ngập lụt thích hợp trên cơ sở các dự án đã, đang và dự kiến
triển khai tại Tp.HCM. Báo cáo khoa học tổng kết đề tài.; 2017.
35
Science & Technology Development Journal – Science of The Earth & Environment, 3(1):23- 36
Original Research
Ho Chi Minh City Univesity of
Technology, VNU-HCM
Correspondence
Le Song Giang, Ho Chi Minh City
Univesity of Technology, VNU-HCM
Email: lsgiang@yahoo.com
History
 Received: 06-11-2018
 Accepted: 30-5-2019
 Published: 22-6-2019
DOI :
https://doi.org/10.32508/stdjsee.v3i1.508
Copyright
© VNU-HCM Press. This is an open-
access article distributed under the
terms of the Creative Commons
Attribution 4.0 International license.
3D numerical modeling of flow and sediment transport in rivers
and open channels
Le Song Giang*, Tran Thi My Hong
ABSTRACT
Numerical model is a useful tool in studying the flow and sediment transport, change in river bed
and is built on solving governing differential equations. Numerical model has many different levels
and three-dimensional model is the highest level, allowing detailed simulation of flow and sedi-
ment transport process in 3D space. The paper presents a method calculating three - dimensional
flow and sediment transport in the open channel. Water level and flow velocity are solved from
three-dimensional equations with hydrostatic hypothesis. Concentration of suspended sediment,
bottom sediment and bottom evolution is solved from transport equations. The governing differ-
ential equations in the "sigma" transform coordinate system are solved by finite volume method
on unstructured grid of quadrilateral elements. Boundary condition of water level or flow will be
imposed on open boundary. For suspended sediment concentrations in the injected phase, sus-
pended sediment concentrations are applied and the outflow phase applies free drainage condi-
tions. This method of calculation was tested with the problem of curved channel sediment trans-
port which was studied experimentally by Odgaard and Bergs. Calculation results are quite consis-
tent with the measured data. In order to test the practical applicability, this method is also tested
with the problem of sediment transport in Cu lao Pho islet on Dong Nai river. To solve the matter
of hydraulic boundary condition of this problem, the model of Cu lao Pho islet is integrated into
the Sai Gon - Dong Nai river system model. Results of the calculation of the river bed evolution
of the Cu lao Pho islet on the Dong Nai river also show that this calculation method gives results
consistent with the rule and can be used in practical research.
Key words: 3D model, sediment trannsport, unstructured grid, the "sigma" coordinate
Cite this article : Giang L S, Hong T T M. 3D numerical modeling of flow and sediment transport in
rivers and open channels. Sci. Tech. Dev. J. - Sci. Earth Environ.; 3(1):23-36.
36