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

NGUYỄN BÁ HƯNG

ỨNG DỤNG PHƯƠNG PHÁP
TỐI ƯU BIẾN PHÂN ĐỂ HIỆU CHỈNH THAM SỐ
CHO BÀI TOÁN LAN TRUYỀN Ô NHIỄM HAI CHIỀU
LUẬN VĂN THẠC SĨ


Ngành : Cơ học kỹ thuật
Chuyên ngành: Cơ học kỹ thuật
Mã số : 60 52 02
LUN VN THC S
Ngời hớng dẫn khoa học: TS. Trần Thu Hà.
H NI 2011
1

MỤC LỤC
DANH MỤC HÌNH VẼ 2
DANH MỤC BẢNG BIỂU 3
MỞ ĐẦU 4
CHƯƠNG I. MÔ HÌNH LAN TRUYỀN CHẤT 2 CHIỀU 6
I.1 Mô hình lan truyền chất hai chiều 6
I.2 Thuật toán tính lan truyền chất 2 chiều 7
I.2.1 Thuật toán giải hệ phương trình dòng chảy 2 chiều 7
I.2.2 Thuật toán giải phương trình truyền tải khuyếch tán 2 chiều 9
I.3 Phát triển mô hình truyền tải đa chất 10

Hình I.5: Bài toán mẫu thủy lực 3. 17
Hình I.6: Kết quả bài toán mẫu thủy lực 3 18
Hình I.7: Điều kiện đầu của bài toán mẫu 19
Hình I.8: So sánh nghiệm tính toán và nghiệm chính xác của bài toán mẫu. 20
Hình I.9: Bản đồ địa hình khu vực hồ Thanh Nhàn 22
Hình I.10: Số liệu địa hình 23
Hình I.11: Kết quả chia lưới 23
Hình I.12: Kết quả tính trường vận tốc 24
Hình I.13: Vị trí điểm đo nồng độ ô nhiễm 24
Hình I.14: Kết quả tính toán mô phỏng BOD
5
. 25
Hình I.15: Kết quả tính toán mô phỏng NH
3
25
Hình I.16: So sánh nồng độ BOD
5
tính toán với số liệu thực đo. 26
Hình I.17: So sánh nồng độ NH
3
tính toán với số liệu thực đo. 26
Hình II.1: BOD
5
trong mô hình tham khảo bài toán mẫu. 43
Hình II.2: BOD
5
trong mô hình bài toán mẫu có hiệu chỉnh tham số 44
Hình II.3: BOD
5
trong mô hình bài toán mẫu không hiệu chỉnh tham số 44

Hình II.14: Kết quả so sánh nồng độ NH
3
với thực đo trong trường hợp 50
có hiệu chỉnh hoặc không hiệu chỉnh tham số tại một số điểm đo trên mặt hồ 50
3
DANH MỤC BẢNG BIỂU
Bảng I.1: Bảng số liệu đo đạc tại 3 điểm trên hồ Thanh Nhàn 21
Bảng II.1: Bảng các hệ số D trong mô hình. 40
Bảng II.2: Bảng các hệ số K trong mô hình. 42
Bảng II.3: Chất lượng nước hồ Thanh Nhàn tháng 11/2001 51
nhiễm bởi nước thải, khí thải và chất thải rắn.
Ở thành phố Hà Nội, tổng lượng nước thải của thành phố lên tới 300.000 -
400.000 m
3
/ngày; hiện mới chỉ có 5/31 bệnh viện có hệ thống xử lý nước thải, chiếm
25% lượng nước thải bệnh viện; 36/400 cơ sở sản xuất có xử lý nước thải; lượng rác
thải sinh hoại chưa được thu gom khoảng 1.200m
3
/ngày đang xả vào các khu đất ven
các hồ, kênh, mương trong nội thành; chỉ số BOD, oxy hoà tan, các chất NH
4
, NO
2
,
NO
3
ở các sông, hồ, mương nội thành đều vượt quá quy định cho phép. Vì vậy việc
giữ gìn và bảo vệ môi trường đă trở nên vấn đề bức xúc
Hiện nay phương pháp thực nghiệm trực tiếp với hệ sinh thái tự nhiên thường
gặp khó khăn và trong nhiều trường hợp là không thể. Cho nên khả năng mô phỏng bài
toán môi trường trong phòng thí nghiệm thường rất hạn chế. Các mô hình toán học là
công cụ cơ bản cho tính toán định lượng cũng như áp dụng vào thực tế khi nghiên cứu
các hệ sinh thái nước. Trên thực tế nhiều kênh sông có chiều dài lớn, thẳng và nước
nông cho nên một số thành phần trong phương trình mô tả có thể được rút gọn. Hiện
nay bài toán giải số các mô hình môi trường đang được đặt ra cấp thiết. Trong luận
văn này trình bày một ứng dụng phương pháp đồng hóa số liệu để hiệu chỉnh bộ tham
số cho bài toán ô nhiễm. Phương pháp này đã được nghiên cứu ứng dụng cho mô hình
tính toán ô nhiễm nước của nhóm nghiên cứu lũ lụt Viện Cơ học trong chương trình
nghiên cứu đề tài cấp Viện Khoa học Việt nam 2009-2010.
Việc tìm hiệu chỉnh ra bộ thông số tối ưu để nâng cao chất lượng tính toán của

6

CHƯƠNG I. MÔ HÌNH LAN TRUYỀN CHẤT 2 CHIỀU

I.1 Mô hình lan truyền chất hai chiều
Mô hình lan truyền chất gây ô nhiễm 2 chiều bao gồm mô hình thủy lực 2 chiều
và mô hình truyền tải khuyếch tán vật chất 2 chiều. Hệ phương trình mô tả dòng chảy
2 chiều được viết trong ([1]-[2]), phương trình truyền tải khuyếch tán vật chất được
viết trong ([3],[4]).
Phương trình truyền tải khuyếch tán vật chất:
)D(f
2
2
2
2
i iiii
i
C
y
C
x
C
y
vC
x
u
t
C







3/42
2/122
)(
hK
vugu
x
z
g
y
u
v
x
u
u
t
u
x
















(I.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ận tốc (trung bình) theo hướng x.
v - vận tốc (trung bình) theo hướng y.
g- gia tốc trọng trường.
K
x
- hệ số Strickler trong lực cản đáy theo hướng x.
K
y
- hệ số Strickler trong lực cản đáy theo hướng y.
C
i
– nồng độ chất gây ô nhiễm.
f
i
– thành phần nguồn.
I.2 Thuật toán tính lan truyền chất 2 chiều
I.2.1 Thuật toán giải hệ phương trình dòng chảy 2 chiều
Chúng tôi đã sử dụng phương pháp khối hữu hạn (FVM - Finite Volume Method)
với lưới tính toán là lưới tam giác. 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 được
chia nhỏ thành các khối có hình học đơn giản thí dụ như các 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ột lưới không cấu trúc. Tính không cấu trúc của lưới được thể hiện qua việc đánh số
các nút lưới và đánh số các khối trong lưới không cần theo một quy luật nào.







v
u
z
V ,












 gz
uv
u
uh
A
2
2
1


























x
u
v
hK
vu

Để xấp xỉ các tích phân trên ta giả thiết rằng trong một khối nhỏ S có thể xấp xỉ
các hàm cần tìm z, u, v và hàm F bằng hằng số. Như vậy, đạo hàm riêng trong tích
phân thứ nhất có thể thay bằng đạo hàm thường theo t; tích phân thứ 2 được tính qua
tích phân đường theo công thức Green.
FSSnBAS
dt
dV


).,( (I.6)
Trong đẳng thức (I.6), S ở số hạng thứ nhất là diện tích của ô lưới. Ở số hạng thứ
2 của (I.6) ∂S là cạnh của ô lưới, n là véctơ pháp tuyến ngoài của ∂S. 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.
I.2.2 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 (I.1) có thể viết lại như sau:
f
y
C
D
yx
C
D
x
CU
t
C


C



(I.7)

),()0,,(
0
yxCyxC 

10

)(),,(
),(
tCtyxC
yx


, 0



n
C

Trong các phương trình trên:
u - vận tốc (trung bình) theo hướng x.
v - vận tốc (trung bình) theo hướng y.
U=(u,v).
D - Hệ số khuếch tán.

x
C
DCUdSC
ht
zz
dS
t
CC
ty
t
tyx
t
txt
tt












.)(.)(
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











 )()()1(

I.3 Phát triển mô hình truyền tải đa chất
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 ví dụ như sau:
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ỉ
11

tiêu và chất cần tính. Thành phần f
i
trong phương trình (I.1) được gọi là thành phần
chuyển đổi (Conversion term) và đã được diễn giải và tính như trong công thức sau:
F
D
(DO) = k
2
.(DO
S
- DO) - k

3
.BOD
5
- Y
4
.(P
o
-R
o
)+ k
6
PO
4

F
D
(NH
3
) = Y
3
.k
3
.BOD
5
- k
4
.NH
3
- Y
5

- hệ số phân huỷ của BOD
5
.
k
4
- tốc độ nitrat hoá.
k
5
- tốc độ phản nitrat hoá.
k
6
hệ số phân hủy PO
4
.
Y
1
- hiệu suất dùng ôxy trong quá trình nitrat.
Y
2
- hàm lượng phốt pho trong hợp chất hữu cơ.
Y
3
- hàm lượng nitơ trong hợp chất hữu cơ.
Y
4
- hệ số hấp thụ phốt phát bởi thực vật.
Y
5
- hệ số hấp thụ amoniac bởi thực vật.
Y

; k
4
= k
40

4
(T-20)
; k
5
= k
50

5
(T-20)

12

Phương trình lan truyền chất (I.1) có thể viết lại như sau:
i
ii
i
i
f
y
C
D
yx




.

Hay là
iii
i
fCDCU
t
C



(I.8)

),()0,,(
0,
yxCyxC
ii


)(),,(
),(
tCtyxC
iyxi


, 0
2

trên và C
i
được tính theo phương pháp sau:
Tích phân 2 vế của phương trình (I.8) theo dx và dy chúng ta có:
 







S
j
m
j
ji
S
ii
S
i
S
i
dxdyCkdxdyCDdxdyCUdxdy
t
C
1
,

Sử dụng công thức Green cho phương trình trên ta có:

Trong tính toán tích phân phương trình trên từng phần tử được viết lại như sau:
dSCkSn
y
C
DvCSn
x
C
DuCdSC
ht
zz
dS
t
CC
n
j
tjjiy
ti
itix
ti
itii
t
tii










13

Do vậy giá trị ô nhiễm được tính như sau:
, ,
,1 , ,2 2, , , , , ,
( ) 1 ( ) ( )
y
i t i t
t x
i i i t i t i n n t i t i t i i t i
n S t
C C
z z n S t
C k C k C k C t C uC D vC D
h x dS y dS
 
 
  
 
          
 
 
 

Ở đây trong phần tử tam giác thứ j với chỉ số tại đỉnh j
1
, j
2
, j

được tính theo công thức sau:
))(())((
))(()(*)(
1231312
12
1
,
3
,13
1
,
2
,,
jjijjjjjj
ii
j
ti
j
tiii
j
ti
j
ti
j
ti
yyxxyyxx
yyCCyyCC
x
C


yyxxyyxx
xxCCxxCC
y
C






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


15

Bài toán mẫu thủy lực 2 (xem [16]): Cho kênh hình chữ nhật với chiều dài 150 m,
rộng 10 m. Zđáy của kênh được cho như sau:
Hình I.3: Bài toán mẫu thủy lực 2.

Đ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 dx = dy = 2,5 m. 16

Kết quả tính mực nước ở trạng thái dừng so với nghiệm chính xác.
0.5
1
1.5
2
2.5
0 50 100 150 200 250 300
ZdayHình I.5: Bài toán mẫu thủy lực 3.

Đầu vào của kênh cho lưu lượng Q=20m3/s,
Đầu ra cho mực nước Z=0.71 m.
Lưới tính được chia với dx = dy = 2 m.
Kết quả tính mực nước ở trạng thái dừng so với nghiệm chính xác.

18

0
0.5
1

x
C
AD
x
A
x
UC
t
C











)(
1
(I.10)
Trong phương trình (I.9):
C: Nồng độ trung bình (trên mặt cắt đang xét) của chất.
U: Vận tốc trung bình của dòng chảy.
A: Diện tích mặt cắt.
D: Hệ số phân tán (Dispersion Coefficient)
Ra: Tốc độ sản sinh ra vật chất đang xét.
Rs: Tốc độ mất đi của chất đang xét.


Hình I.7: Điều kiện đầu của bài toán mẫu.
20

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 (trong phương trình (I .10) ta có D = 0, Ra = 0 và Rs
=0). Trong bài toán mẫu này phương trình (I.10) có dạng phương trình tải thuần túy:
0
dC dC
U
dt dx
 

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à
2
0
( )
0
( )
x Ut x
k
C x e
  


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.
Trong mô hình 2D không sử dụng khái niệm lưu lượng Q, nhưng các đại lượng qx =
uh và qy =vh được sử dụng trong hệ phương trình (I.5) (tương đương với hệ phương
trình (I.2) – (I .4)).

/ngày-đêm.
1 Đặc tính thủy lý:
- Nhiệt độ nước hồ dao động từ 23,4
o
C đến 23,8
o
C trong mùa khô.
- Độ pH dao động ở mức kiềm ( 7,88 - 8, 57 ) và không chênh lệch nhau nhiều.
- Độ đục dao động từ 9 - 34 mg/l .
- Độ dẫn không thay đổi tại các trạm thu mẫu 0,04 S/m .
- Độ muối ( NaCL) không thay đổi tại các điểm thu mẫu 0.01%
- Hàm lượng ô xy hoà tan dao động từ 7,9 mg/l - 11,6 mg/l .
2 Đặc tính thủy hóa :
Các số liệu đo đạc được thể hiện trên bảng sau:
Điểm thu mẫu
P1
Cổng chảy vào ở
phía Bắc hồ
P2
Khu vực giữa hồ
P3
Cổng chảy ra
phía Bệnh viện
Thanh Nhàn
BOD
5
20
o
C (mg/l) 24.0 15.5 17.0
NH

5
, NH
3
để đưa vào nghiên cứu khả năng mô phỏng của mô
hình. Quá trình hóa học để phân hủy chất gây ô nhiễm được mô tả như sau:
f
1
(BOD
5
) = - k
1
.BOD
5

f
2
(NH
3
) = Y
3
.k
1
.BOD
5
- k
3
.NH
3
- Y
4

=Y
3
k
1
-Y
4
k
1

Trong đó:
- BOD
5
- Nhu cầu ôxy sinh học.
- NH
3
- Hàm lượng Nitơ trong nước dưới dạng hợp chất.
Dựa trên các số liệu thu thập, tác giả thiết lập các số liệu đầu vào cho mô hình. Các
số liệu địa hình được phân tách thành các vùng biên (bờ) và địa hình lòng hồ sau đó
tiến hành chia lưới thành các lưới tam giác không cấu trúc bao gồm: 1058 nút và 1964
phần tử. Hình I.10: Số liệu địa hình Hình I.11: Kết quả chia lưới

Kịch bản tính toán được xây dựng như sau:
Lưu lượng vào ra khỏi hồ là 2100 m3/ngày đêm (bằng với số liệu khảo sát tại cửa
vào).
Mực nước hồ lấy trung bình là 4m.
Chỉ tiêu BOD
5