ĐẠI HỌC QUỐC GIA HÀ NỘI
TRƯỜNG ĐẠI HỌC CÔNG NGHỆ

NGUYỄN HỒNG PHONG

MỘT SỐ VẤN ĐỀ KHOA HỌC TRONG
NGHIÊN CỨU MÔ PHỎNG VÀ TÍNH TOÁN
CHẤT LƯỢNG NƯỚC
Ngành: Cơ học kỹ thuật
Chuyên ngành: Cơ học kỹ thuật
Mã số: 60 52 02

TÓM TẮT LUẬN VĂN THẠC SĨ NGƯỜI HƯỚNG DẪN KHOA HỌC: TS. Trần Thu Hà

Hà Nội - 2012

DANH MỤC BẢNG BIỂU VÀ HÌNH VẼ
Bảng 3.1: Số liệu đo đạc các chỉ tiêu ô nhiễm được sử dụng trong mô hình.
Bảng 3.2: So sánh kết quả tính toán tại vị trí có điểm đo đạc.
Hình 1.1: Lưới không cấu trúc dạng tam giác.
Hình 1.2: Quan hệ giữa các phần tử trong miền.
Hình 2.1: Cấu trúc lưới tính được sử dụng trong các bài toán mẫu.
Hình 2.2 : So sánh mực nước tính toán và nghiệm chính xác theo chiều dài kênh
Hình 2.3: Sai số tuyệt đối của mực nước giữa nghiệm chính xác và tính toán theo
chiều dài kênh.
Hình 2.4: Cao độ đáy của kênh dẫn trong bài toán mẫu thủy lực 2
Hình 2.5: So sánh mực nước tính toán và nghiệm chính xác theo chiều dài kênh
Hình 2.6: Sai số tuyệt đối của mực nước giữa nghiệm chính xác và tính toán theo
chiều dài kênh.
Hình 2.7: Cao độ đáy kênh của bài toán mẫu thủy lực 3.
Hình 2.8: So sánh mực nước tính toán và nghiệm chính xác theo chiều dài kênh.
Hình 2.9: Sai số tuyệt đối của mực nước giữa nghiệm chính xác và tính toán theo
chiều dài kênh.
Hình 2.10: Điều kiện đầu của bài toán mẫu 1.
Hình 2.11: So sánh nghiệm tính toán và nghiệm chính xác bài toán mẫu 1.
Hình 2.12: Phân bố của chất ô nhiễm dọc kênh tại thời điểm 1000s.
Hình 2.13: So sánh kết quả tính toán và nghiệm chính xác của bài toán mẫu 2.
Hình 2.14: Kết quả tính toán so sánh với nghiệm chính xác của bài toán mẫu 3
tương ứng với 3 hệ số k khác nhau.
Hình 2.15: Sai số tương đối (%) của bài toán mẫu tương ứng với 3 giá trị khác
nhau của hệ số k.
Hình 3.1: Bản đồ địa hình khu vực hồ Thanh Nhàn.
Hình 3.2: Xử lý số liệu địa hình trên lưới không cấu trúc.
Hình 3.3a: Kết quả tính trường vận tốc
Hình 3.3b: Vị trí điểm đo nồng độ ô nhiễm
Hình 3.4: Kết quả tính BOD

22
22
23

24

26
28
29
29
30
31
32
33
34
35
2
MỞ ĐẦU
Sự phát triển nhanh chóng của nền kinh tế sản xuất và tiêu dùng hiện nay của chúng ta
gây sức ép trực tiếp và gián tiếp lên mức độ ô nhiễm của môi trường sống. Đặc biệt là
trong các thành phố lớn, nơi tập trung các khu sản xuất và tiêu dùng các sản phẩm
công nghiệp, nông nghiệp, dịch vụ phục vụ cuộc sống của dân cư. Tình trạng ô nhiễm
môi trường nói chung và ô nhiễm môi trường nước từ trước đến nay đã xuất hiện và
nhiều nơi diễn ra khá nghiêm trọng. Bên cạnh việc ô nhiễm sông ngòi thì ô nhiễm các
nguồn nước và hồ chứa trong đô thị do nguồn nước thải của sinh hoạt và sản xuất cũng
rất đáng lo ngại vì chúng ảnh hưởng trực tiếp đến sức khỏe của bộ phận không nhỏ dân
cư sinh sống quanh hồ và ảnh hưởng cảnh quan đô thị. Do đặc điểm địa hình của lòng

nhiễm. Mô hình WASP cũng có thể liên kết với các mô hình thủy động lực và vận
chuyển trầm tích để tính toán ra được trường dòng chảy, nhiệt độ, độ muối và các
thông lượng trầm tích.
- Mô hình QUAL2K (hay Q2K) (River and Stream Water Quality Model) được
phát triển từ mô hình QUAL2E (hay Q2E (Brown và Barnwell 7 1987)). Đây là mô
hình mô phỏng chất lượng nước suối và sông 1 chiều. Mô hình tính toán chu trình Nitơ
và thông qua các chu trình chuyển hóa nitơ để biểu diễn các hợp chất cacbon, các loại
cacbon hữu cơ không sống
- Bộ phần mềm MIKE (DHI) được phát triển trên môi trường đồ họa của hệ
điều hành Windows bao gồm nhiều modul tính toán trong đó có cả tính toán chất
lượng nước. MIKE cũng tích hợp tốt với hệ thống thông tin bản đồ địa lý (GIS) và
được nhiều đề tài và nghiên cứu trong nước chọn để sử dụng như công cụ hỗ trợ tính
toán.
1.1.2. Các nghiên cứu ở trong nước.
Ở trong nước gần đây cũng có nhiều nghiên cứu quan tâm đến việc sử dụng mô
hình để mô phỏng và tính toán chất lượng nước. Các đề tài và nghiên cứu này chia làm
hai hướng:

4
- Hướng nghiên cứu thứ nhất sử dụng mô hình và phần mềm nước ngoài để mô
phỏng và tính toán một số yếu tố của chất lượng nước trong một khu vực cụ thể.
Hướng này thiên về nghiên cứu ứng dụng của mô hình và tác động của chất ô nhiễm
đến khu vực nghiên cứu để phục vụ nhu cầu nghiên cứu, đánh giá tổn thất do ô nhiễm
hoặc hoạch định chính sách phát triển của địa phương. Ví dụ chương trình hợp tác với
Cơ quan hợp tác Quốc tế 9 Nhật Bản - JICA của Viện Tài nguyên và Môi trường biển
– Viện Khoa học và Công nghệ Việt Nam sử dụng phần mềm chuyên dụng
CABARET của LOICZ (Mỹ) để đánh giá mức độ tích tụ và khuếch tán vật chất tại
một số điểm thuộc vịnh Hạ Long. Một số nghiên cứu khác sử dụng mô hình eco Lab -
MIKE
1

theo chiều sâu và điều kiện áp suất tuân theo quy luật thuỷ tĩnh. Nếu bỏ qua thành
phần nhớt và rối, ảnh hưởng của gió và ảnh hưởng của lực Coriolis, hệ phương trình S.
Venant 2D có thể viết dưới dạng sau([1]-[2]).
Phương trình liên tục:

   
0








vh
y
uh
xt
z
(1)
Phương trình động lượng theo hướng x:

x
x
gS
hK
vugu
x
h

y
y
gS
hK
vugv
y
h
g
y
v
v
x
v
u
t
v
,0
3/42
2/122
)(




















(4)
Trong các phương trình trên:
z – mực nước, h là độ sâu dòng chảy.
u,v – vận tốc.
g – gia tốc trọng trường.
K
x
, K
y
– hệ số Strickler trong lực cản đáy.
C
i
– nồng độ chỉ tiêu ô nhiễm thứ i.
f
i
– thành phần nguồn.
S
o,x
,S
o,y
– độ dốc đáy.

Biên có dòng chảy ra khỏi miền (outflow): u
n
 0, dòng chảy theo chiều pháp tuyến
ngoài.
- Trạng thái chảy êm: u
n
< c: cần cho 1 điều kiện biên
- Trạng thái chảy xiết: u
n
 c : không cần cho điều kiện biên.
Đối với phương trình truyền tải khuyếch tán vật chất cần cho một điều kiện biên vào
và một điều kiện biên ra. Ví dụ đối với bài toán thủy lực có dòng chảy vào êm, cần cho
một điều kiện biên vào và một điều kiện biên ra, để bài toán ổn định và giảm sai số thì
điều kiện biên vào thường được chọn là vận tốc và điều kiện biên ra được chọn là cao
trình mực nước. Với bài toán truyền tải khuyếch tán thì điều kiện biên vào thường
chọn là nồng độ chỉ tiêu ô nhiễm và điều kiện biên ra là đạo hàm của nồng độ chỉ tiêu
ô nhiễm.

1.3. Thuật toán giải hệ phương trình dòng chảy 2 chiều.

Hệ phương trình Saint – Venant 2D (1) - (3) được giải số theo phương pháp
khối hữu hạn (FVM - Finite Volume Method) trên lưới không cấu trúc dạng tam giác
cho một đoạn sông hoặc hồ chứa. Lưới tam giác mô tả tốt dòng chảy trên miền có địa
hình phức tạp (cả biên và đáy) . Trong phương pháp khối hữu hạn, miền tính toán D
được chia nhỏ thành các khối có hình học đơn giản S
i
dạng tam giác và các khối này
tạo thành một lưới phủ kín miền tính toán. Các tham số để mô tả một lưới là nút lưới
và liên kết giữa các nút lưới để tạo thành các khối. Trên hình dưới đây là một phần của


M
i,j
, nếu S
i
nằm cạnh S
j
, M
i,j
là cạnh chung của 2 phần tử này.

8

Hình 1.2: Quan hệ giữa các phần tử trong miền.
Ở đây ta dùng các ký hiệu của lý thuyết tập hợp như  là hợp,  là giao,  là tập
rỗng.
i
S
là phủ của S
i
, như vậy một đỉnh của phần tử chỉ có thể nằm ở một trong hai
đầu của cạnh, không có đỉnh của phần tử này nằm trên cạnh của phần tử khác.
Tính không cấu trúc của lưới được thể hiện qua việc đánh số nút k và đánh số phân tử
của các phần tử không theo một quy luật nào cả. Những lưới được đánh số nút (hoặc
phần tử) theo 2 chỉ số i, j, trong đó i theo chiều x, j theo chiều y thường phải liền kề
nhau và có trật tự thì được gọi là lưới có cấu trúc. Lưới không cấu trúc mỗi nút, mỗi
phần tử được đánh số bằng một chỉ số và không có một điều kiện ràng buộc nào cho
chỉ số này.
Trong khi chia lưới miền D, phần tử không bắt buộc có cùng một cấu trúc hình học là
tam giác hoặc tứ giác. Thí dụ miền D có thể chia bằng một lưới gồm các tam giác và
các tứ giác. Tính chất hình học của một phần tử cũng không tuân theo một nguyên tắc

b
x



,0
,
y
z
S
b
y



,0
,
Hệ phương trình (1)-(3) có thể viết lại dưới dạng sau:

F
y
B
x
A
t
V













 gz
uv
u
q
A
x
2
2
1
,

























x
u
v
hK
vu
gv
y
v
u
hK
vu
guF
y
y
3/42
22
3/42
22

(6)
Trong đẳng thức (6), S ở số hạng thứ nhất là diện tích của phần tử. Ở số hạng thứ 2
của (6) M là cạnh của phần tử, n là véctơ pháp tuyến ngoài của M. Trong các bài toán
có quá trình biến đổi chậm, để tính gần đúng các tích phân ta có thể sử dụng giá trị
trung bình tại khối S cho các hàm. Trong tích phân này còn có các đạo hàm riêng bậc
nhất. Đạo hàm này được xấp xỉ bằng giá trị của độ dốc của mặt phẳng đi qua các đỉnh
của tam giác. Sơ đồ cụ thể của thuật toán được trình bày chi tiết ở trong báo cáo tổng
kết đề tài [1].
1.4. Thuật toán giải phương trình truyền tải khuyếch tán 2 chiều.
Phương trình lan truyền chất (4) cho một chất có thể viết lại như sau:

f
y
C
D
yx
C
D
x
UC
t
C









Điều kiện đầu :
),()0,,(
0
yxCyxC 

Điều kiện biên :
)(),,(
),(
tCtyxC
vyx


và
0


r
n
C

Trong các phương trình trên:
U=(u,v) – Vận tốc.
D – Hệ số khuếch tán.
C – Nồng độ chất ô nhiễm.
f = kC – Thành phần nguồn ô nhiễm.
k – Hệ số hao tán.

v
, 
r

.
.
Khai triển tích phân phương trình trên từng phần tử được viết lại như sau

dSkCdMn
y
C
DdMnCv
dMn
x
C
DdMnCudS
t
CC
ty
iNj
t
y
iNj
t
x
iNj
t
x
iNj
t
t



































tt
)()(
)()(
)1(

1.5. Phát triển mô hình truyền tải đa chất.
Mô hình lan truyền chất ban đầu chỉ tính toán với một chất ô nhiễm cụ thể, khi mô
phỏng các chất gây ô nhiễm trong hồ chứa, chúng ta phải xét đến một số chất và chỉ
tiêu đánh giá chất lượng nước. Do đó cần phải mở rộng thêm khả năng tính toán nhiều
chất đồng thời và có tương tác qua lại lẫn nhau. Tùy thuộc vào số lượng các chất tham
gia trong mô phỏng chúng ta cần xây dựng một kịch bản các phản ứng và quá trình
biến đổi giữa các chất trong tập hợp các chỉ tiêu và chất cần tính.
Phương trình lan truyền chất (7) có thể viết lại như sau:
i
ii
i
i
f
y
C
D
yx
C
D
x
UC
t
C



fDgradCUC
t
C



.).(
(9)

12
Điều kiện đầu :
),()0,,(
0,
yxCyxC
ii


Điều kiện biên :
)(),,(
),(
tCtyxC
ivyxi


và
0


r
i






S
j
m
j
ji
S
ii
S
i
S
i
dxdyCkdxdygradCDdxdyUCdxdy
t
C
1
,
.

Sử dụng công thức Green cho phương trình trên ta có:







CC
n
j
tjjiy
ti
itix
ti
iti
tii











1
,,
,
,
,
,
,
.)(.)(

Tại đây các giá trị hàm được tính tại điểm trung tâm của phần tử, dS là diện tích

ti
iti







 )()(
,
,
,
,

Ở đây trong phần tử tam giác thứ j với chỉ số tại đỉnh j
1
, j
2
, j
3
và các tọa độ (x
j1
,y
j1
),
(x
j2
,y
j2

,
2
,,
jjijjjjjj
ii
j
ti
j
tiii
j
ti
j
ti
j
ti
yyxxyyxx
yyCCyyCC
x
C




13
))(())((
))(()(*)(
12131312



Với
3
,
3
,
1
,
,,
j
ti
j
ti
j
ti
CCC
là nồng độ C
i,t
tại các điểm j
1
, j
2
, j
3

1.6. Hiệu chỉnh một số tính toán của thuật toán.
Trong khi thực hiện các tính toán mô phỏng lan truyền chất ô nhiễm với nhiều
chất thì sự bổ sung thành phần fi dẫn đến thuật giải kém ổn định và khuyếch tán số
nhất là với những bài toán kiểm định có thành phần vận tốc u và hệ số D lớn. Chính vì























)(
2/1
)(
2/1
)(
2/1
)(
2/1

C
x
C
tt
5.0
2/1
thì (11) có thể phân tích
thành công thức lặp.
ty
iNj
n
t
y
iNj
t
n
x
iNj
n
t
x
iNj
t
n
t
n
CCdMn
y
C
y
































1
)(
1
;5.0
2
5.0
2
)1(

Phép tính sẽ dừng lại khi giá trị


1nn
CC
. Ngoài ra một số hàm tính toán bổ trợ
ví dụ như xác định giá trị vector normal của phần tử, giá trị diện tích của phần tử được
hiệu chỉnh lại để tính nhanh và chính xác hơn. Chi tiết phần lập trình và chú giải các
phần tính toán này được trình bày cụ thể trong phần phụ lục của luận văn.
KẾT LUẬN CHƯƠNG
Phần 1-4 của chương này đã trình bày các phương trình chính và thuật toán giải
của mô hình IMECH_2DUSZ, mô hình lan truyền chất ô nhiễm đã thực hiện trong

14
chương trình nghiên cứu cơ bản. Dựa trên những kết quả nghiên cứu này, học viên đã
tiếp tục cải tiến và mở rộng mô hình sang hướng tính toán được nhiều chất hơn và hiệu
chỉnh lại một số tính toán trong phần tính lan truyền chất ô nhiễm của mô hình được
trình bày trong phần 5-6 của chương.
CHƯƠNG 2. CÁC TÍNH TOÁN KIỂM ĐỊNH CHẤT LƯỢNG MÔ HÌNH
Đối với bất cứ mô hình tính toán nào thì việc kiểm định mô hình là công việc
cần thiết để xác định độ chính xác và tin cậy của mô hình thông qua các bài toán mẫu

Điều kiện biên của kênh được cho như sau:
* Đầu vào của kênh cho lưu lượng Q=20m3/s
* Đầu ra cho mực nước Z=0.8 m
Lưới tính được chia với tam giác có chiều dài cạnh góc vuông = 2,5m.
Kết quả tính mực nước ở trạng thái dừng so với nghiệm chính xác.

Hình 2.5: So sánh mực nước tính toán và nghiệm chính xác theo chiều dài kênh.

17 Hình 2.6: Sai số tuyệt đối của mực nước giữa nghiệm chính xác và tính toán theo chiều
dài kênh.
Bài toán mẫu thủy lực 3 [8]:
Cho kênh hình chữ nhật, chiều dài 300 m, chiều rộng 10 m. Zđáy của kênh được cho
như sau:

Hình 2.7: Cao độ đáy kênh của bài toán mẫu thủy lực 3.
Điều kiện biên của kênh được cho như sau:
Zday
0
0.5
1
1.5
2
2.5
0
50
100
150

C
AD
xAx
UC
t
C











)(
1

Trong phương trình trên:
C: Nồng độ trung bình (theo tiết diện ngang) của chất.
U: Vận tốc.
A: Diện tích mặt cắt.
D: Hệ số khuyếch tán.
f = kC – thành phần nguồn ô nhiễm.

Cho một kênh có chiều dài L = 10000 mét. Dòng chảy trên kênh là dòng dừng đều với
vận tốc dòng chảy U = 0.5 mét/giây. Tại thời điểm bắt đầu tính toán (t
0

DK_Dau

Hình 2.10: Điều kiện đầu của bài toán mẫu 1.
Giả thiết tiếp rằng chất C không bị khuếch tán, không được sản sinh ra và
không bị mất đi trong quá trình tính toán ta có D = 0, f =0). Trong bài toán mẫu này có
dạng phương trình tải thuần túy:
0
dx
dC
U
dt
dC

Có thể kiểm tra ngay được rằng nghiệm chính xác của bài toán mẫu này là

k
xUtx
exC
2
0
)(
0
)(




Bài toán mẫu trình bày ở trên có thể sử dụng để kiểm định mô hình 2D bằng
cách giả thiết rằng độ sâu cột nước và vận tốc dòng chảy không thay đổi theo chiều y.
Dưới đây là kết quả kiểm định tại thời điểm t = 3600 giây. Tại thời điểm này,














trong đó erfc (Complementary error function -

dtexerfc
x
t




2
2
)(


Với t=1000s,
C
1

2.5
-5000 -3000 -1000 1000 3000 5000
C00

Hình 2.12: Phân bố của chất ô nhiễm dọc kênh tại thời điểm 1000s.
Kết quả tính toán kiểm định tại thời điểm t=3000s so sánh với nghiệm chính xác.

Hình 2.13: So sánh kết quả tính toán và nghiệm chính xác của bài toán mẫu 2.
Bài toán mẫu 3[10]:
Cho một kênh có chiều dài L = 10000 mét (0<x<10000).
0
0.5
1
1.5
2
2.5
-5000
-2500
0
2500
5000
X
nongdo
NghiemCX
Nghiem_TT

23
Dòng chảy trên kênh là dòng dừng đều với vận tốc dòng chảy U = 0.03 mét/giây. Tại
thời điểm bắt đầu tính toán cho một chất hòa tan có nồng độ C=30 bắt đầu đổ vào kênh
ở biên vào. Tại biên ra chọn điều kiện biên như công thức (10). Theo [10], nghiệm của

Trong đó D là hệ số khuyếch tán, k là hệ số hao tán của chất ô nhiễm.
Nếu chọn D= 0.0000017 (m2/s) và k
1
=0 ; k
2
=1 (ngày
-1
) ~ 1.16x10
-5
(s
-1
); k
3
=2 (ngày
-1
)
~ 2.31x10
-5
(s
-1
) ta có kết quả tính toán với mô hình tại thời điểm t=555000 giây (150
giờ).

Hình 2.14: Kết quả tính toán so sánh với nghiệm chính xác của bài toán mẫu 3
tương ứng với 3 hệ số k khác nhau.

0
5
10
15