MÔ HÌNH 1D CHO LAN TRUYỀN CHẤT Ô NHIỄM HỮU CƠ TRÊN KÊNH
SÔNG (SALBOD) VÀ MỘT SỐ ỨNG DỤNG
Pgs. Ts. NguyÔn tÊt ®¾c
Viện Quy Hoạch Thủy lợi miền Nam
I. Mở đầu
Trong những năm gần đây do sự gia tăng dân số và sự phát triển nhanh chóng của các đô thị và
trung tâm công nghiệp, lượng nước thải chưa được xử lí được thải trực tiếp vào mạng lưới kênh
sông ngày càng gia tăng và làm suy thóai nghiêm trọng chất lượng, ảnh hưởng xấu tới sinh họat và
phát triển kinh tế. Để có biện pháp cải tạo thích hợp thì mô hình tóan là một công cụ trợ giúp đắc
lực và không thể thiếu cho các nhà qui họach và ra quyết định. Mô hình tính tóan dòng chảy đã
được các chuyên gia trong và ngòai nước phát triển và ứng dụng nhiều cho hệ thống sông ngòi Việt
nam, đặc biệt là hệ thống kênh sông của Đồng Bằng Cửu Long và hệ thống sông Đồng Nai-Sài
Gòn. Tuy nhiên bài tóan lan truyền chất chủ yếu được xây dựng và áp dụng cho lan truyền mặn.
Thuật tóan dùng để giải bài tóan lan truyền chất một chiều cơ bản là sơ đồ sai phân. Tuy nhiên do
hiện tượng khuếch tán số của sơ đồ sai phân (nhất là sơ đồ sai phân trung tâm) làm cho kết quả tính
thường không bảo tòan ý nghĩa vật lý của hiện tượng mô phỏng, chẳng hạn nồng độ mặn bị âm,
hoặc không có nguồn mà giá trị tính tóan tại mặt cắt sát biên thường lớn hơn giá trị biên, hoặc buộc
nồng độ tại các điểm sát hợp lưu của các nhánh khác nhau phải bằng nhau. Các phần mềm thương
mại lớn như Mike 11 hoặc ISIS vẫn chưa khắc phục được nhược điểm này.
Nhằm khắc phục các nhược điểm trong thuật tóan tính lan truyền chất, mà chủ yếu là lan truyền
mặn, của các phần mềm hiện có, trong phần mềm SALBOD, việc sử dụng phương pháp phân rã kết
hợp với phương pháp đường đặc trưng cho phương trình tải đã khắc phục được các nhược điểm nêu
trên. Mô hình SALBOD gồm phần tính dòng chảy, tính lan truyền chất (mặn, BOD, DO) và được
thực hiện đồng thời trong một bước thời gian tính, tuy nhiên trong báo cáo này chỉ trình bầy thuật
tóan cho phần tính lan truyền chất. Chi tiết phần tính dòng chảy có thề xem trong [1].
II. Cơ sở học thuật của mô hình SALBOD
Hệ phương trình cơ bản
A. Phần tính tóan dòng chảy:
+ Hệ phương trình cơ bản: Hệ phương trình Saint-Venant một chiều sau đây được sử dụng phổ
biến để tính dòng chảy trong kênh sông:
)1(q

ARC
QQg
x
Z
gA
A
Q
xt
Q
Trong đó : B – là chiều rộng mặt nước: A – là diện tích mặt cắt ngang; Z- là mực nước so với
một cao độ chuẩn; Q là lưu lượng qua mặt cắt ngang; g- là gia tốc trọng trường; C- là hệ số cản
Chézy; R- là bán kính thủy lực; q- là lưu lượng gia nhập trên một đơn vị chiều dài dọc sông (như
bơm, xả nước thải, ); t- là thời gian; x- là tọa độ dọc sông
B. Phần tính toán chất lượng nước:
B1. Hệ phương trình cơ bản
Một số yếu tố của chất lượng nước cũng được tính tóan đồng thời sau mỗi bước tính dòng chảy.
Phương trình cơ bản cho mặn, BOD và DO như liệt kê dưới đây:
+ Mặn với nồng độ S
)4(
2
2
a
A
SQSq
S
A
QQq
x
S
E

KK
x
B
E
x
B
U
t
B
ssrq
rsr
+
+
++
++−
∂
∂
=
∂
∂
+
∂
∂
+ DO với nồng độ D:
)4(.)(
12
2
2
c
A

Trong 3 phương trình (4a,b,c) nêu trên q , Sq, Bq, Dq tương ứng là lưu lượng nước thải gia nhập,
nồng độ mặn, BOD và DO trong dòng gia nhập q; Qsr , Ss, Bs, Ds tương ứng là lưu lượng trao đổi
(trên một đơn vị độ dài) sông-ruộng, độ mặn, BOD và DO trong Qsr. Cũng tương ứng như vậy với
lưu lượng mưa Qr. Các hệ số chuyển hóa K
1
, K
3
của BOD và K
2
là hệ số thấm khí. Trong (4b) nếu
K1=K3=0 ta có dạng (4a). D
sa
là nồng độ bão hòa Oxy trong nước. Cả 3 phương trình trên đều có
thể viết ở dạng chung như sau:
)5(.
2
2
φσ
+−
∂
∂
=
∂
∂
+
∂
∂
C
x
C

x
C
E
t
C
∂
∂
=
∂
∂
B2. Một số sơ đồ sai phân dùng để giải phương trình tải (6)
Để giải (6) thường áp dụng một số sơ đồ sai phân sau đây (để đơn giản ta bỏ qua số hạng nguồn
trong vế phải). Nếu dùng các ký hiệu tóan tử:
x
C
CL
t
C
CL
∂
∂
=
∂
∂
=
21
Khi đó tương ứng với sơ đồ sai phân theo hướng (upwind) ta có:
2
( )
( )

n+1
< 0 (8)
( )
( )
( )
[ ]
n
i
n
i
n
i
n
i
ii
CCCC
xx
CL
1
1
1
1
1
2
1
1
−
+
−
+

))(1()([)(
11
1
1
1
1
1
112
n
i
n
i
n
i
n
iii
CCCCxxCL
−+
+
−
+
+
−
−+
−−+−−=
θθ
(11)
)]())(21()([
1
1

=
ααα
(12)
Trong (12) trọng số α gọi là hệ số phân tách . Bằng cách khai triển Taylor các số hạng trong (11)
và (12) quanh điểm x
i
, t
n
và sử dụng ký hiệu δ = x
i+1
– x
i-1
ta có:

!3
)(
2
)(
2
)(
2
1
+
∆
++
∆
=
∂
∂
−

UCUL
ttxtxxx
θθ
δ
(14)
Trong vế phải của (13) và (14) ta thấy có các số hạng
tt
C
t
2
∆
và
xx
C
U
2
δ
. Số hạng thứ nhất biến
phương trình xuất phát thành phương trình sóng và tạo ra các dao động của nghiệm khi giải số, còn
số hạng thứ hai tạo thêm hệ số khuếch tán do số
2
δ
U
(thường gọi hệ số khuếch tán số). Nếu U và
bước lưới không gian δ lớn thì hệ số khuếch tán số lớn và làm mất ý nghĩa của hiện tượng mô
phỏng. Cũng áp dụng cách khai riển Taylor tương tự với (8),(9),(10) ta cũng có kết luận tương tự.
B3. Sử dụng phương pháp đường đặc trưng trong giải (6)
Với bài tóan truyền tải thuần túy, theo quan điểm Lagrange, có thể hiểu bài tóan như sau: Một
hạt lỏng ở thời điểm t nằm tại điểm A, với vận tốc U hạt lỏng sẽ di chuyển theo một quỹ đạo nào đó
để đạt tới điểm B vào thời điểm t + ∆t (xem hình vẽ dưới). Trong quá trình tải thuần túy hạt lỏng

Trong khỏang tích phân ∆t , xem U chỉ là hàm của x, và được lấy trọng số giũa 2 lớp thời gian n
và n+1 ; giá trị của U được cho tại các điểm x
1
, x
2
, … , x
N
(tính từ mô hình dòng chảy trong cùng
một bước thời gian), giũa các khỏang [ x
i
, x
i+1
] vận tốc U được xem như một hàm tuyến tính theo
x: U = ax + b .
Do U không phụ thuộc t (trong một bước ∆t) cho nên trong mỗi khỏang [ x
k
, x
k+1
] các đường
cong tích phân dx = U(x)dt đều song song:
Cho nên để xác định chân đường đặc trưng chỉ cần tính δt với từng khỏang [ x
k
, x
k+1
] , sau đó
kiểm tra điều kiện:
Nếu điều kiện này thỏa mãn có nghĩa là đặc trưng đã cắt đường thẳng t = t
n
; khi đó tính
Để tính δt, δx (xem hình vẽ) ta có thuật tóan sau

Dễ thấy đường đặc trưng không phải là đường thẳng và có thể cắt nhiều mắt lưới.
+ Nội suy giá trị chân đường đặc trưng A: Vì điểm A có thể trùng với một điểm lưới, tuy nhiên
nói chung A nằm giữa các mắt lưới, vì thế cần tiến hành nội suy giá trị tại A qua giá trị đã biết tại
các điểm lưới. Cách nội suy tuyến tính có thể gây ra khuếch tán số, vì thế trong SALBOD đã dùng
hàm nội suy spline bậc 3 ( phiên bản trước dùng nội suy kết hợp tuyến tính với nội suy Lagrange).
Xem chi tiết trong [1].
B4. Ví dụ giải thích các phương pháp
Để so sánh độ chính xác của các sơ đồ sai phân và phương pháp đường đặc trưng áp dụng cho
phương trình tải ta xét một ví dụ chính xác có nghiệm giải tích sau đây [1]:
Xét phương trình
0),( =
∂
∂
+
∂
∂
x
f
txu
t
f
với x ∈ (0,1) và t > 0
4
const
xU
x
t
=
∆
=

ε
ε
Trong ú ; (Xột trng hp a = ẳ)
iu kin u : f(x,0) = 3sin (4x)
iu kin biờn :










=
ta
ta
tf
sin21
sin
4sin3),0(
;








-2
-1
0
1
2
3
4
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Toaù ủoọ
Gớa trũ
f fL fs fH1
fH2 fT1 fT2
Cỏc ký hiu c s dng trờn biu :
f : l nghim chớnh xỏc; f
L
: l li gii bng phng phỏp c trng s dng phng phỏp ni suy
tuyn tớnh kt hp vi ni suy Lagrange.
f
S
: li gii c trng dựng ni suy Spline bc 3. (f
H1,
f
H2
) : li gii bng phng phỏp sai phõn
theo hng vi giỏ tr khỏc nhau ca trng s . f
T
(f
T1
,f
T2

trong đó C
N
là nồng độ tại các mặt cắt chảy ra khỏi hợp lưu, C
i
, Q
i
v
là nồng độ và lưu lượng tại các
mặt cắt áp sát hợp lưu của các nhánh chảy vào hợp lưu (tất nhiên chảy vào hoặc chảy ra là tùy thuộc
từng thời điểm). Q
j
R
là lưu lượng tại các mặt cắt áp sát hợp lưu của các nhánh chảy ra

Chú ý rằng:
Bằng phương pháp nêu trên vật chất lan truyền tới đâu tính tới đó cho nên tiết kiệm được thời
gian tính tóan trên máy, khác với các phương pháp sai phân là luôn phải giải hệ phương trình đại số
cho tòan mạng, tốn thời gian
III. Một số áp dụng
Mô hình SALBOD khi mới được xây dựng chủ yếu được sử dụng cho các bài tóan lan truyền
mặn, phiên bản lúc đó có tên là SAL. Do yêu cầu đánh giá tác động môi trường của các dự án công
nghiệp có nước thải, mô hình SALBOD được hòan thiện và được sử dụng cho nhiều dự án mà chủ
yếu là các dự án trên sông Thị vải (Ve Dan, Marubeni, Phú Mỹ, ) và các dự án trên hệ thống sông
Đồng Nai-Sàigòn. Do khuôn khổ của một báo cáo, dưới đây chỉ giới thiệu một ví dụ tính tóan ảnh
hưởng của hồ Phước Hòa đến xâm nhập mặn và ô nhiễm vùng hạ du hệ thống sông Đồng Nai-Sài
gòn.
Theo thiết kế, ngòai việc đóng góp vào sơ đồ lưới điện, hồ Phước Hòa còn đóng góp 50m3/s cho
nhiệm vụ cấp nước của hồ Dầu tiếng. Tuy nhiên khi có đập Phước Hòa lượng nước xả xuống sông
Bé trong mùa khô chỉ còn khỏang 15m3/s, nhỏ hơn lưu lượng tự nhiên của sông Bé trong mùa khô
của những năm gần đây, đặc biệt là từ khi có Thác mơ. Việc phát triển các bậc thang thủy điện ở

R
j
i
i
v
i
N
Q
CQ
C
∑∑
=
j
R
j
i
v
i
QQ
Nước hồi quy từ sông Vàm Cỏ Đông và Vàm Cỏ Tây cũng được xét như lưu lượng gia nhập vào
hệ thống từ các sông thuôc Đòng Bằng sông Cử Long.
Sơ đồ, số liệu biên của mô hình:
+ Sơ đồ tính toán (Xem hình vẽ): Gồm toàn bộ các sông Đồng Nai (từ Trị An), sông Sài Gòn (từ
chân đập Dầu tiếng), sông Vàm Cỏ Đông (từ Bến Đá), sông Vàm Cỏ Tây (từ Bình Châu), toàn bộ
mạng kênh rạch thành phố Hồ Chí Minh (Nhiêu Lộc, Thị Nghè, Kênh Đôi, Kênh Tẻ, Tân Hoá, Lò
Gốm, Bến Nghé, Tham Lương, Bến Cát, Vàm Thuật, …), mạng kênh phía nam thành phố và các
sông vùng duyên hải (Bến Lức, Cần Giuộc, Ông Lớn, Cây khô, Thị Vải, Đồng Tranh, Lòng Tầu,
Cái Mép, ). Các biên của sơ đồ tính là lưu lượng Trị An, sông Bé, Dầu Tiếng, Bén Đá, Bình Châu;
mực nước thực đo tại Vũng Tầu có tính tương quan với các cửa Đồng Tranh, Sòai Rạp (Xem bảng
2).

Biên mô hình Sông Loại số liệu Thời gian quan
trắc
Ghi chú
Hạ lưu đập Trị An Đồng Nai Dòng chảy 1/3-30/4/ 2003 Số liệu đo đạc được tại nhà
máy Trị An
Trạm lưu lượng
Phước Hòa
Bé Dòng chảy 1/3-30/4/ 2003 Số liệu đo đạc được tại
trạm Phước Hòa
Hạ lưu đập Dầu
Tiếng
Sài Gòn Dòng chảy 1/3-30/4/ 2003 Số liệu đo đạc được tại đập
Dầu Tiếng
Bình Châu Vàm Cỏ Tây Dòng chảy 1/3-30/4/ 2003 Tính từ mô hình ĐBSCL
Bến Đá Vàm Cỏ Đông Dòng chảy 1/3-30/4/ 2003 Số liệu tính tóan
Thị Tính Nhánh của sông
Sài Gòn
Dòng chảy 1/3-30/4/ 2003 Số liệu tính tóan
Cửa Đồng Tranh
(ltương quan Vũng
Tàu )
Đồng Tranh Mực nước giờ 1/3-30/4/ 2003 Biên hạ lưu
Cửa Soài Rạp
(Tương quan Vũng
tầu)
Đồng Nai Mực nước giờ 1/3-30/4/ 2003 Biên hạ lưu
Nhà Bè Đồng Nai Mực nước giờ
và độ mặn
1/3-30/4/ 2003 Vị trí kiểm tra mô hình
Biên Hòa Đồng Nai Mực nước giờ 1/3-30/4/ 2003 Vị trí kiểm tra mô hình

Đồng Nai.
Từ kết quả tính toán này một khuyến nghị quan trọng là quá trình xây dựng đập Phước Hoà
(ngăn dần dòng sông Bé) phải kết hợp với việc phát triển các bậc thang thuỷ điện trên sông Đồng
Nai để sao cho tổng lượng nước ở ngã ba Đồng Nai-Sông Bé bình quân không dưới 290m3/s.
Về ô nhiễm: Để tính toán ta xem nồng độ BOD nền là 2mg/L (có thể giá trị này hơi cao so với
một số giá trị đo đạc, tuy nhiên các giá trị đo dưới 1mg/L khó có độ tin cậy khi so với các giá trị DO
tương ứng tại cùng vị trí và thời điểm đo đạc, mặt khác mục tiêu là đánh giá xu thế biến đổi nên ta
dùng giá trị BOD là 2mg/L làm nền).
9
Bảng 4. BOD (mg/L) bình quân trong hai tháng 3 và 4 tại các vị trí quan trọng
trên sông Đồng Nai, Sài Gòn, tương ứng với các phương án nêu trên

Trong bảng 4 có thêm phương án ký hiệu là PD. Trong phương án này có đập Phước Hoà nhưng
chưa có nguồn nước bổ xung (như phương án PB), tất cả các nguồn thải thuộc khu vực Biên Hoà
buộc xử lí để BOD trong nguồn chỉ còn 50mg/L, các nguồn thải thuộc các khu vực khác được xử lí
để BOD trong nguồn đạt 100mg/L.
Từ bảng 4 có thể thấy do các nguồn ô nhiễm tập trung ở khu vực Biên Hoà và nội thành thành
phố Hồ Chí Minh, hơn nữa khi có đập Phước Hoà thì về cơ bản sông Sài Gòn được thêm nước từ
các nguồn hồi qui của hồ Dầu Tiếng nên mức độ BOD bình quân trên sông Sài Gòn hầu như không
thay đổi. Còn đối với sồng Đồng Nai, đặc biệt từ khu vực Hoá An xuống tới Long Đại, thì tình
trạng ô nhiễm có gia tăng. Khi gia tăng lượng nước lấy cấp cho các nhà máy nước nhưng chưa ngăn
đập Phước Hoà (phương án PA) thì BOD bình quân tại Hoá An tăng từ 3.52mg/L lên 3.59mg/L,
hoặc tại Cầu Đồng Nai tăng từ 4.68mg/L lên 4.77mg/L. Sự gia tăng được xem như không đáng kể,
nhưng khi có dập Phước Hoà (phương án PB) thì sự gia tăng tương ứng là 3.52 lên 4.48 và 4.68 lên
5.74 sẽ là đáng kể.
Trong phương án PC khi có thêm 80m3/s bổ xung từ thượng lưu thì bức tranh ô nhiễm gần như
trở lại bình thường như phương án PA.
Với phương án PD khi các nguồn thải được xử lí thì nước sông không còn bị ô nhiễm mặc dù có
gia tăng nước lấy cấp cho dân sinh và công nông nghiệp.
Một số nhận xét:

Summary
1D SALBOD MODEL FOR WASTE WATER SPREADING IN THE CANAL SYSTEM
AND ITS APPLICATION
Ass. Prof. Dr. NGUYEN TAT DAC
Southern Institute for Water resources planning
It is presented in this study the theoretical background of 1D SALBOD Model for waste water
spreading in the canal system. The analysis of numerical schemes applied to 1D transport equation
shows that numerical diffusion encountered by finite difference methode leads to, physically,
sometimes unacceptable solutions. To overcome this shortcoming it is suggested to use, in SALBOD
model, splitting approach with characteristic method applied to pure transport. For illustration the
results of application of the developed SALBOD model to the estimation of impact caused by Phuoc
Hoa dam is also briefly presented.
Người phản biện: GS. TSKH Nguyễn Ân Niên
11