Mô hình số về vận tải bùn cát ba chiều trong sông

46 Trần Thị Mỹ Hồng, Lê Song Giang 
MÔ HÌNH SỐ VỀ VẬN TẢI BÙN CÁT BA CHIỀU TRONG SÔNG 
3D NUMERICAL MODELING OF SEDIMENT TRANSPORT IN RIVERS 
Trần Thị Mỹ Hồng, Lê Song Giang 
Trường Đại học Bách khoa – Đại học Quốc gia Tp.HCM; tranthimyhong@hcmut.edu.vn, lsgiang@yahoo.com 
Tóm tắt - Bài báo trình bày một phương pháp tính toán dòng chảy và 
vận tải bùn cát trong sông, kênh hở trong không gian ba chiều. Các 
phương trình vi phân cơ bản để tính toán dòng chảy và vận tải bùn cát 
được viết trong hệ tọa độ biến đổi “sigma”. Các phương trình này đượ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ới các 
phần tử tứ giác. Phương pháp tính này được tính toán kiểm tra với bài 
toán vận tải bùn cát trong đoạn kênh cong được nghiên cứu thực 
nghiệm bởi Odgaard và Bergs (1988) và áp dụng thử để làm mô hình 
vận tải bùn cát tại khu vực Cù lao Phố trên sông Đồng Nai. Các kết 
quả tính toán cho thấy phương pháp này cho kết quả tính khá chính 
xác và có thể ứng dụng cho các tính toán thực tế. 
Abstract - The article presents a method of calculating three 
dimensional flow and sediment transport in rivers and open 
channels. The governing differential equations for flow and 
sediment transport are written in the "sigma" coordinate system. 
These equations are solved by the finite volume method on the 
unstructured grid with quadrilateral elements. The method is tested 
with the sediment transport in curved channel studied 
experimentally by Odgaard and Bergs (1988). The method is also 
applied to modeling sediment transport at Pho islet in DongNai 
river. Calculation results show that this method gives fairly accurate 
results and can be applied to practical calculations. 
Từ khóa - mô hình số 3D; vân tải bùn cát; lưới phi cấu trúc; sông 
Đồng Nai; phương pháp thể tích hữu hạn 
Key words - 3D numerical model; sediment trannsport; 
unstructured grid; DongNai river; finite volume method 
1. 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 [9] 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 đứng. Lin và Fanconer [4] giải phương trình 
Navier-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 đứng. Wu và ctg. [10] 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 trung. Hiện 
nay cũng có một số phần mềm để tính dòng chảy và vận tải 
bùn cát ba chiều như DELFT3D, TELEMAC 3D tuy 
nhiên thuật toán của các phần mềm này chưa được tối ưu và 
chi phí bản quyền khá cao. 
Trong bào 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 này đảm bảo được các tính chất bảo toàn và 
lưới tính toán có khả năng bám sát theo địa hình, phù hợp 
với các bài toán thực tiễn. 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ệm bởi Odgaard và Bergs [5] và 
bài toán bồi xói khu vực Cù lao Phố trên sông Đồng Nai. 
2. Mô hình toán 
2.1. Phương trình cơ bản 
Dòng chảy được giải từ phương trình dòng chảy ba 
chiều với giả thiết thủy tĩnh. Trong hệ tọa độ biến đổi 
“sigma”, phương trình được viết : 
( ) 0=


++





q
t
D
q
 (1) 
  FUUUqUq DgD
D
A
qDA
t
V
H +−= 


−


+−+





 (2) 
Còn bùn cát lơ lửng và diễn biến đáy được giải từ các 
phương trình: 
  ( ) 00 = 


−−


+−+




C
D
D
CwqCHDC
t
q V
sH
C q
 (3) 
( ) ( ) ( )bbbbbC ED
t
bc
t
z
P −=+


+


− q 1
 (4) 
Trong đó:  – mực nước;  Tyx uu ,=U – hai thành phần 
vận tốc của dòng chảy trên phương ngang; ω – thành phần 
vận tốc trên phương thẳng đứng “sigma”; 
  Uq Dqq Tyx == , và  Dq = ; D – độ sâu; σ – 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; DCqC = ; ws0 – vận tốc lắng; zb – cao độ đáy; C – 
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;  Tbybxb qq ,=q
– 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 đá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 Smagorinsky [6]: 
5,0
222
xy2
1
yx 


+


+ 


+ 


 =
vuvu
yxCAH
 (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) [15, 16]: 
kLCAV  = (6) 
Trong đó C' - hằng số mô hình, được 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à 
ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, VOL. 17, NO. 1.1, 2019 47 
Gerritsen [1] như sau: 
( ) ( ) ( ) ( ) 
−+= huhu
c
k sb 1
1 2
*
2
*
 (7) 
( ) hhDL −= 1 (8) 
Với c=0.09 – hằng số mô hình; =0.4 – hằng số Karman; 
su* - vận tốc ma sát trên mặt nước; 
bu* - dạng sửa chữa của 
vận tốc ma sát đáy bu* , trong đó ( )buu ** ,max=
b
*
u ; 
*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à ( ) Dzh −=  - độ 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ính : 
 (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 Rijn [7 - 
9], qb và cb,e, được tính 
 (10) 
3.0
*
5.1
50
, 015.0
D
T
b
d
c Ceb =
 (11) 
Với: ( )b cr crT   = − (12) 
 (13) 
Trong đó: d50 – đường kính hạt 50%; s – tỷ trọng hạt; τ 
và τcr – ứ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, Eb, và lắng, Db, được 
tính theo Hayter và Mehta [2] và Krone [3]: 
−= 1


 bbE
 (14) 
−=
d
b
sb CwD


10
 (15) 
Trong đó:  - hệ số xói; b - ứng suất tiếp đáy;  - ứng 
suất tiếp đáy ngưỡng xói; d - ứng suất tiếp đáy ngưỡng bồi. 
2.2. Đ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ụng: 
 (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 (= 0), các điều kiện biên sau được sử dụng: 
0 = (17) 
( ) 000 ,, 

yx
V vu
D
A
−= 




 (18) 
00 =


+

C
D
D
Cw Vs
 (19) 
Còn tại đáy (= -1): 
0= (20) 
( ) ( )vuvuCvu
D
A
D
V ,,
5.022 += 




 (21) 
bb
V
s DE
C
D
D
Cw −= 


+−

0
 (22) 
Trong đó CD là hệ số ma sát đáy và được tính : 
( )
2
0ln
=
z
CD

 (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. 
2.3. Phương pháp giải 
Các phương trình (1) - (3) được giải bằng phương pháp 
thể tích hữu hạn theo sơ đồ được cải tiến từ sơ đồ được giới 
thiệu trong tài liệu [13]. 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 qσ đượ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. 
a) Trên mặt bằng b) Trên phương thẳng đứng 
Hình 1. Sơ đồ lưới tính trên mặt bằng (a) và các lớp lưới tính 
theo phương thẳng đứng (b) 
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 η ở 
thời điểm n+1 
• Bước 3: Tính lưu lượng qσ ở 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: 
  + −+ 


kk S
H
S
dSdDAdSd
t 


 UqU
q
( )ebbsbb ccwED ,0 −=−
( )  1.23.0*
5.1
50
5.0
1053.0 TDdgsq Cb
−−= 
( )
3/1
250*
1 
−=

g
sdD
0= nC
48 Trần Thị Mỹ Hồng, Lê Song Giang 
( ) +−= 


−


kk SS
V dSdDgDdSd
D
A
q



 

F
U
U
 (24) 
Sau khi sai phân số hạng đạo hàm theo thời gian, (24) 
đưa tới phương trình: 
( ) ( ) ( ) knknknknkknknk VFlux
t
V
 =−++
− +++
−−+
rJJUqq
2/12/1
1
2/12/12/1
 (25) 
Trong đó 
SV kk .= (26) 
2/1
1
+
+ 


−=
k
V
k
D
A
q
D
S


q
qJ
 (27) 
( ) 


−=
L k
Hnkk dl
n
DAFlux
U
UqU 
 (28) 
( )
kk
gD Fr +−=  (29) 
Các số hạng Fuxk(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 Jk+1 sẽ được 
nội suy và sai phân thành: 
kkPkkTk
AA qqJ
1111 ++++
 + −=
 (30) 
Trong đó 
1kT
A
+
 và
1kP
A
+
 đượ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: 
k
n
kkN
n
kkP
n
kkS
AAA Srcqqq =−+− ++
++
−
2/1
1
2/12/1
1 (31) 
Trong đó: 
( )
kTkP
k
kP
AA
t
V
A + +
=
+1
 (32) 
1+
 =
kTkT
AA
 (33) 
kPkS
AA =
 (34) 
( ) nkknknkkk VFlux
t
V
rUqSrc .2/12/1 +−
= −−
 (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 η, phương trình liên tục (1) được tích 
phân trên phương  trên toàn bộ chiều sâu. Kết quả ta thu được: 
( ) 0
t
D
k
kk = +


q
 (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: 

−=



+
+
k
j
j
2/1n
jknk
2/1n
C Lq
S
1
t
D

 (37) 
Từ (37), độ sâu và mực nước tại nút C được tính: 
t
t
D
DD
2/1n
Cn
C
1n
C 


+=
+
+
 (38) 
1n
Cb
1n
C Dz
++ +=
 (39) 
c. Bước 3 
Để tính lưu lượng đơn vị trên phương , 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: 
0dSd
q
dSddSd
t
D
SSS kkk
=


++









 q
 (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 qσ ở nút 
C trên các lớp ở thời điểm n+1/2: 

+
+
++
+
−


−=
j
j
2/1n
jkn
k
2/1n
C
k
2/1n
k
2/1n
1k
Lq
St
D
qq


 (41) 
Lưu ý là tại đáy 0
2/1
1
=
+n
q 
d. Bước 4 
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
H
S
C
kk
dSdCHDCdSd
t
q



 q
 = 


−


+
S
C
S
V
kk
dSdDSdSd
C
D
D
Cq

 

 (42) 
Sau khi sai phân số hạng đạo hàm theo thời gian (42) 
đưa tới phương trình: 
( ) ( ) ( ) ( )1111
1
1 . ++++
+
+ − =−++−
 n
k
n
k
n
k
n
k
n
k
n
kC
n
kC
k sCrDVJJCFluxqq
t
V
 (43) 
Trong đó: 
S.V kk  = (44) 
1n
2/1k
V1n
1k
C
D
D
CqSJ
+
+
+
+ 


−=


 (45) 
( )   −= +
S
k
n
H
nn
k
n
k dSCHDCCFlux 
2/1
q
 (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 Jk+1 sẽ được nội suy và sai 
phân thành: 
kkPkkTk
CACAJ
1111 ++++
 + −=
 (47) 
Trong đó 
1kT
A
+
 và
1kP
A
+
 đượ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 (47) vào (43), ta được phương trình 
cuối cùng: 
k
n
kkN
n
kkP
n
kkS
SrcCACACA =−+− ++
++
−
1
1
11
1 (48) 
Với 
( )
kTkP
n
kkP
AADVs
t
A + + 
+
= +
+
1
11
 (49) 
1+
 =
kTkT
AA
 (50) 
ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, VOL. 17, NO. 1.1, 2019 49 
kPkS
AA =
 (51) 
( )nknnk
n
kk CFluxDrC
t
D
VSrc −
+
 = +1.
 (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): 
( )
( )
( ) −=+

+


−
S
bb
S
b
S
b
S
b
C dSEDdSdS
t
bc
dS
t
z
P q 1
 (53) 
Sau khi sai phân số hạng đạo hàm theo thời gian, (53) 
đưa tới phương trình: 
( ) ( ) ( )SEDqFluxDC
t
zz
P bbb
n
b
n
b
C −=++
−
−
+1
1 
 (54) 
Trong đó 
( ) ( ) nbnb bcbc
t
S
DC −
=
+1
 (55) 
( ) =
L
nbb
dlqqFlux
 (56) 
Từ (54) ta có cao độ đáy ở nút C ở thời điểm n+1: 
( )
( ) ( ) bbb
C
n
b
n
b qFluxDCSED
P
t
zz −−−
−
+=+
 1
1
 (57) 
3. Tính toán thử nghiệm 
3.1. Vận tải bùn cát trong kênh cong 180o 
Odgaard và Bergs [5] có một nghiên cứu bằng thực 
nghiệm về 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 đã Wu và ctv [10] tính toán bằng 
mô hình 3D. Kênh gồm một đoạn cong 180o, bán kính trung 
bình 13,11m, và 2 đoạn kênh dẫn vào và ra thẳng, dài 20m. 
Mặt cắt ngang kênh hình thang, đáy rộng 1.44m và mặt 
thoáng rộng 2.44m (Hình 2). Đá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 ra theo dòng chảy được bơm ngược 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ới trạng thái ổn định, 
địa hình đáy kênh đã được đo đạc. 
Hình 2. Lưới tính mô hình 3D 
Trong nghiên cứu này mô hình số 3D cho kênh của 
Odgaard và Bergs đã được xây dựng. Trên mặt bằng, đoạn 
cong 180o được chia thành 180×16 phần tử còn mỗi đoạn 
kênh dẫn vào và ra được chia thành 80×16 phần tử (Hình 
2). Mô hình có 11 lớp theo phương thẳng đứng. Địa hình 
đáy đượ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. 
Đầu vào của kênh được áp đặt điều kiện biên lưu lượng 
Q=0,153m3/s với 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. 
a) Vận tốc trên mặt phẳng mặt cắt b) Vận tốc vuông góc với mặt cắt 
Hình 3. Vận tốc tại các mặt cắt ở các góc cong khác nhau 
Bờ trái Bờ phải 
Măt cắt 00 
Bờ trái Bờ phải 
Măt cắt 00 
Bờ trái Bờ phải 
Măt cắt 450 
Bờ trái Bờ phải 
Măt cắt 450 
Bờ trái Bờ phải 
Măt cắt 900 
Bờ trái Bờ phải 
Măt cắt 900 
Bờ phải 
Măt cắt 1350 
Bờ trái Bờ phải 
Măt cắt 1350 
Bờ trái 
Măt cắt 1800 
Bờ trái Bờ phải 
Măt cắt 1800 
Bờ trái Bờ phải 
50 Trần Thị Mỹ Hồng, Lê Song Giang 
Hình 4. Địa hình đáy kênh tính toán (sau 9 giờ - đường liền và 
sau 11 giờ - đường đứt nét) và thí nghiệm (chấm tròn) 
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 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ện một cách rõ ràng. Hình 4 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 hoàn toàn phù hợp với số liệu 
thí nghiệm của Odgaard và Bergs. 
3.2. Vận tải bùn cát khu vực Cù lao Phố trên sông Đồng Nai 
Sông Đồng Nai tại khu vực Cù lao Phố từng có một dự 
án san lấp lấn sông gây tranh cãi. Do địa hình lòng sông 
uốn lượn phức tạp, việc sử dụng mô hình 2D để nghiên cứu 
diễn biến bồi xói, đặc biệt khu vực đầu Cù lao Phố [11, 15] 
là không phù hợp vì mô hình không có khả năng tính toán 
mô phỏng các dòng thứ cấp. Để đánh giá khả năng của mô 
hình 3D phát triển trong bài báo, mô hình đã được áp dụng 
tính toán cho đoạn sông này. 
Trên mặt bằng, đoạn sông được chia thành 7261 phần 
tử tứ giác và 5 lớp theo chiều sâu. Số phần tử theo chiều 
ngang nhánh sông chính là 16 còn nhánh phụ là 8. Các phần 
tử có kích thước khoảng 16 – 25 m. 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 đáy mô hình được tham khảo từ số liệu khảo 
sát bởi Viện Khoa học Thủy lợi miền Nam, thực hiện năm 
2008 [14]. Nghiên cứu [14] cũng cho thấy vật liệu đáy sông 
khu vực 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 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ó đường 
kính trung bình d50=0,35mm được xem xét. 
Để giải quyết vấn đề điều kiện biên thủy lực của mô hình 
3D, mô hình được tích hợp vào mô hình hệ thống sông Đồng 
Nai, kế thừa từ kết quả nghiên cứu [12]. Do mô hình 3D và 
mô hình hệ thống sông Đồng Nai sau tích hợp đã trở thành 
một mô hình thống nhất nên mực nước và vận tốc trên biên 
của mô hình 3D sẽ được tính ngay trong mô hình tích hợp. 
Đ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. Khảo sát của Viện 
Khoa học Thủy lợi Miền Nam [14] cho thấy 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. 
Để đảm bảo điều kiện ổn định bước thời gian tính được 
lấy t = 0,6s. Mặc dù bước thời gian tính khá nhỏ nhưng 
tốc độ tính toán của mô hình khá tốt. Trên máy PC Core i5 
3,2GHz, một giờ máy tính mô phỏng số được khoảng 8 giờ. 
Do mô hình hệ thống sông Đồng Nai đã được hiệu 
chỉnh khá tốt [12] và mô hình 3D có quy mô khá nhỏ nên 
ở quy mô toàn mô hình mô hình tích hợp từ hai mô hình 
này không cần phải hiệu chỉnh lại. Ngược lại, việc tinh 
chỉnh mô hình 3D với vận tốc điểm và hiệu chỉnh mô hình 
bùn cát 3D là cần thiết. Tuy nhiên do điều kiện không cho 
phép nên việc này chưa thực hiện được. Điều này có thể 
chấp nhận vì tính toán chỉ giới hạn ở mục tiêu chứng minh 
khả năng áp dụng mô hình 3D ở quy mô bài toán thực. 
a) Thành phần vận tốc tiếp tuyến với mặt cắt b) Thành phần vận tốc vuông góc với mặt cắt 
Hình 5. Trường vận tốc trên các mặt cắt ngang lúc triều xuống vào thời gian đỉnh lũ 2008
Tính toán mô phỏng đã được thực hiện trong thời gian 
từ ngày 30/8/2008-30/9/2008. Khoảng thời gian này được 
chọn vì trong thời gian đó đã có một trận lũ lớn xảy ra. 
Phân bố vận tốc tại một số mặt cắt ngang lúc triều xuống 
-0,50
-0,40
-0,30
-0,20
-0,10
0,00
-0,500,000,50
h
, m
x/b
-0,50
-0,40
-0,30
-0,20
-0,10
0,00
-0,500,000,50
h
, m
x/b
-0,50
-0,40
-0,30
-0,20
-0,10
0,00
-0,500,000,50
h
, m
x/b
-0,50
-0,40
-0,30
-0,20
-0,10
0,00
-0,500,000,50
h
, m
x/b
-0,50
-0,40
-0,30
-0,20
-0,10
0,00
-0,500,000,50
h
, m
x/b
-0,50
-0,40
-0,30
-0,20
-0,10
0,00
-0,50-0,30-0,100,100,300,50
h
, m
x/b
Mặt cắt 178o 
Mặt cắt 137o 
Mặt cắt 84o 
Mặt cắt 2o 
Mặt cắt 31o 
Mặt cắt 61o 
Măt căt 1-1 
Bờ trái Bờ phải Bờ trái Bờ phải 
Măt căt 1-1 
Măt căt 2-2 
Bờ trái Bờ phải 
Măt căt 2-2 
Bờ trái Bờ phải 
Măt căt 3-3 
Bờ trái Bờ phải 
Măt căt 3-3 
Bờ phải Bờ trái 
Măt căt 4-4 
Bờ phải Bờ trái 
Măt căt 4-4 
Bờ trái Bờ phải 
ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, VOL. 17, NO. 1.1, 2019 51 
vào thời gian lũ cũng được giới thiệu trên Hình 5. Vị trí các 
mặt cắt được giới thiệu trên Hình 6. 
Hình 6. Vị trí các mặt cắt 
Độ bồi xói khu vực đầu Cù lao Phố được giới thiệu trên 
Hình 7. Để 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 8 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 7. 
Hình 7. Trường vận tốc và độ bồi xói tại khu vực đầu cù lao 
Phố (đơn vị thang màu: mét) 
Hình 8. Phân bố vận tốc và biến đối đáy tại các lát cắt 
(đường màu đen: đáy sông ban đầu; đường màu đỏ: đáy sông 
sau một tháng) 
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ăng mô phỏng các dòng thứ cấp này cũng chính là ưu thế 
của mô hình 3D so với mô hình 2D. 
4. Kết luận 
Bài báo đã trình bày mộ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ế. 
Lời cảm ơn 
Nghiên cứu này được tài trợ bởi Trường Đại học Bách 
Khoa TpHCM trong đề tài ‘‘Nghiên cứu tính toán dòng 
chảy và chuyển tải bùn cát bằng mô hình 3D và ứng dụng 
tính toán cho khu vực cù lao Phố, sông Đồng Nai” mã số 
T- KTXD-2017-75 
TÀI LIỆU THAM KHẢO 
[1] Davies, A.M. and Gerritsen H., (1994). An Intercomparison of 
Three-dimensional Tidal Hydrodynamic Models of the Irish Sea, 
Tellus, 46A, pp. 200-221. 
[2] Hayter, E.J. and Mehta, A.J., 1986. Modeling cohesive sediments 
transport in estuarine waters. Applied Mathematical Modelling 51, 
pp. 765-778. 
[3] Krone, R.B., 1962. Flume studies of the transport of sediment in 
estuarine shoaling processes. Final report to San Fransico District 
U.S. Army Corps of Engineers, Washington D.C 
[4] Lin, B.; R.A. Falconer (1997). Three-dimensional layer integrated 
modeling of estuarine flows with flooding and drying. Estuar. Coast. 
Shelf Sci., 44, 737–751. 
[5] Odgaard, A.J. and M.A. Bergs (1988). Flow Processes in a Curved 
Alluvial Channel, Water Resour. Res., 24(1): 45-56. 
[6] Smagorinsky, J. (1963). General Circulation Experiments with the 
Primitive Equations. I: The Basic Experiment. Monthly Weather 
Rev., 91, 99 - 164. 
[7] Van Rijn., J. Hydraulic Eng., ASCE, Vol. 110 (10), pp. 1431-1456. 
(1984a) 
[8] Van Rijn., J. Hydraulic Eng., ASCE, Vol. 110 (11), pp. 1613-1641. 
(1984b) 
[9] Van Rijn, L.C. (1993). Principles of Sediment Transport in Rivers, 
Estuaries and Coastal Seas, Aqua Publications, 690p. 
[10] Wu, W., W. Rodi and T. Wenka (2000). 3D Numerical modeling of 
Flow and Sediment Transport in Open Channel, J. Hydraul. Eng., 
126(1): 4-15. 
[11] Cty cổ phần ĐT-KT-XD Toàn Thịnh Phát (2014). Báo cáo đánh giá 
tác động mô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. 
[12] Lê Song Giang (2017). 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. 
[13] Trần Thị Mỹ Hồng, Lê Song Giang (2017). 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, trg. 
[14] Viện Khoa học Thủy lợi miền Nam (2009). Đánh giá tác động dòng 
chảy sông Đồng Nai đoạn từ cầu Hóa An đến cầu Ghềnh thuộc thành 
phố Biên Hòa. Tp. Hồ Chí Minh, tháng 10 năm 2009. 
[15] Kolmogorov, A. N. (1942), "Equations of Turbulent Motion of an 
Incompressible Fluid," Izvestia Academy of Sciences, USSR; 
Physics, Vol. 6, Nos. 1 and 2, pp. 56-58. 
[16] Prandtl, L. (1945), "Uber ein neues Formelsystem fur die 
ausgebildete Turbulenz," Nacr. Akad. Wiss. Gottingen, Math-Phys. 
Kl. 1945, pp6-19.
(BBT nhận bài: 31/10/2018, hoàn tất thủ tục phản biện: 15/01/2019) 
1 
1 
2 
2 
3 
3 
4 
4 
6 
6 5 
5 
Thượng lưu Hạ lưu 
Lát cắt 5 - 5 
Thượng lưu Hạ lưu 
Lát cắt 6 - 6