ĐẠI HỌC QUÓC GIA HÀ NỘI
TRƯỜNG ĐẠI HỌC KHOA HỌC Tự NHIÊN
—OỈ-CD-SO—
TÊN ĐẺ TẢI
MÔ HÌNH SỐ TRỊ BA CHIÈU ĐÓNG KÍN RỐI
k-8

TRONG TÍNH TOÁN LAN TRUYỀN Ô NHIẺM
MÃ SÓ: QT-07-02
CHỦ TRÍ ĐÈ TÀI: TS TRẦN VÁN cúc
Những ngừô tham gia chính:
NCS:LƯU QUANG HƯNG
T s TRÀN VÃN TRẢN
.M , HOC Q U Ô C GIA HÀ NÕI
lí ƯNG TÀÍV THÔNG TIN THƯ VIỀN
1 ) 1 / 1 5 4
HÀ NỘI, 2007
1
TÓM TẮT ĐÈ TÀI
1- Tên đề tài: Mô hình số trị ba chiểu đóng kín rối k-S trong tính toán ỉan
truyên ô nhiêm.
2- Mã số: QT07- 02
3- Chủ trì đè tài: TS.Trần Văn Cúc
4- Những người tham gia chính: NCS Lưu Quang Hưng ,TS. Trần Văn Tràn.
5- Mục tiêu và nội dung nghiên cứu:
- Xây đựng hệ phương trình đầy đủ cho bài toán lan truyền chất tan ở biên
theo phương pháp khép kín rối k- £ trong hệ toạ độ sig m a
- Chọn thuật toán để giải số bài toán.
- Mô phỏng kểt quả trong một số trường hợp đơn giản.
6- Các kết quả đã đạt được:
- Đã thiết lập được hệ phương trình đầy đủ cho bài toán lan truyền chất tan

- Finding the algorimth to solve numericaly for these mathematics.
- Simulating the result in some simple cases.
6- Obtained result:
- Establing fulfill simutaneous equation for problem on dispersion of
solutes in coastal areas according to the k-C turbullent closure method in
Sigma coordinate system.
- Fiding copatible algorimth and making numerical solution of copatible
problem
- The result is to be appeared in the Journal of Vietnam national University,
Hanoi: "Caculating the 3D dispersion of pollutant solutes in coastal areas
using k-G turbullent closure in Sigma coordinate system ”
- A scientific mater thesis.
7-Financials:
a) Receiving: 20.000.000d
b) Spending:
-For research work: 12.000.000d
-For conferences and seminars: 3.000.000d
-For other works: 5.000.000d
Total : 20.000.000d
MỤC LỤC
Mở đầu 7
Chương 1 8
Mô hình tính lan truyền chất tan ven biển, ba chiều, khép kin rối
k -£
trong hệ toạ dộ sigma
8
1.1 Mò hình lan truyền chất tan ven biển, ba chiều, khép kín rối k -£ trong hệ toạ
độ Descartes 8
1.2 Hệ toạ độ sigm a


2.4 Hệ phương trinh khép kín roi 56
2.5 Hệ phương trình trạng thái, nhiệt độ, độ m ặn 60
2.6 Phương trình khuếch tán chất tan 62
2.7 Hệ phương trình tích phàn 63
2 8 Tiêu chuẩn ổn định 66
Chương 3 68
Mô phỏng 68
3.1 Đặt bài toán 68
3.2 Kết quả tính toán 69
Kết luận 80
Tài liệii tham khảo 83
Phụ lục 89
5
MỞ ĐẦU
Vùng ven biển nước ta, đặc biệt là những nơi có cửa sông, vịnh, hàng ngày
có một số lượng lớn tạp chất đổ ra biển. Bài toán tìm quy luật phân bố tạp chât cho
những miền đó rất có ý nghĩa thực tiễn,vì thế đã có nhiều công trình nghiên cứu bài
toán này .Trong các hướng tiếp cận để giải bài toán, hướng nghiên cứu bằng phương
pháp mô hình hoá số trị có nhiều ưu thế hơn, do dựa trên cơ sở toán học chặt chẽ và
kết quả nghiên cứu có khả năng triển khai ứng dụng. Cho đến nay,một số tác giả đã
sử dụng mô hình sổ trị sẵn có của nước ngoài[4 ], hoặc tự xây dựng nhưng sử dụng
phương pháp sai phân hiện đơn giản [18 ] trong việc giải bài toán tìm phân bố tạp
chất trong chất lỏng nói chung và trong vùng ven biển nói riêng.
Khi nghiên cứu bài toán trên, nhiều tác giả đã sừ dụng hệ phương trình mô
tả bài toán trong hệ toạ độ Đềcác.Nhưng trong nhiều trường hợp hệ toạ độ đó không
mô tả được thuận lợi, nhất là vùng có địa hình đáy phức tạp.
Trong đề tài này ,để khắc phục nhược điểm nói trên chủng tôi sử dụng hệ
toạ độ sigma cho bài ba chiều (3D) đầy đủ, gồm một hệ phương trình động lực học
có tính đến ảnh hường cùa gió, ma sát đáy và nhớt rối, kết hợp đồng thời với việc
giải phương trình khuếch tán tạp chất,phương trình trạng thái.phương trình mô tả

trong đó p là khối lượng riêng cùa nước biển, ữ; là khối lượng rièna trune bình. 9 là gia
tốc trọng trường, p là áp suất. / là tham số Coriolis. được tính từ phượng trình chuyển
động quay của Trái Đât
7
/ = 2D cos <p
(a.5)
p = PflfW - 9ữữn - 9 pci:
(a.6)
với là vận tốc quay, 9 là góc vĩ độ, là áp suất khí quyển trên mặt thoáng tại điểm
tính. Ngoài ra, áp suất tại độ sâu - có thề thu được bằng cách tích phàn phương trinh
chuyển động của hạt lỏng theo phương thẳng đứng (a.4) từ độ sâu điểm đang xét tới mặt
thoáng như biểu diễn (a.6). Hệ số khuếch tán rối D:-i được xác định qua hệ số nhớt tuyệt đối
/■' và hệ sổ nhớt rối ílz
Trong hệ phương trình trên, bên cạnh V-1 ^ ' p là các giá trị cần tìm cùa mô hình, thì còn
các đại lượng chưa biết là khối lượng riêng p, phương trình dao động mặt thoáng rK và hệ
số nhớt rối f*t.
b) Hệ phưong trình khép kín rối
Hệ phương trình khép kín rối hai phương trình f bao gồm phương trình cho
động năng rối K và tốc độ hao tán năng lượng rối f . Hệ phương khép kín rối ’* — £ được
xây dụng trên cơ sờ tham số hoá các đại lượng mạch động chưa biết bàng biểu thức của
động năng rối k và tốc độ hao tán năng lượng rối £ rồi đưa về các phương trình năng lượng
rối. Phương trình cho động năng rối bằng cách lấy trung bình Reynolds cho toán tử Navier-
Stokes, sau đó tham số hoá các phi tuyến bậc cao của mạch động. Việc khép kín mô hình
rối đòi hỏi phải đưa vào phương trình liên quan đến hao tán năng lượng, trong mô hình này
là tốc độ hao tán năng lượng rối Chi tiết cách thiết lập mô hình khép kín rối - < có thể
xem trong Wilcox D.c. (1994). Các đậc trung và tính chất của mô hình k — £ được phân
tích trong Mohammadi B., Pironneau o. (1994).
Các già thiết cơ bản để thiết lập hệ phương trình trong mô hình bao gồm
i) Giả thiết đẳng hướng địa phương cùa chuyển động rối: nhớt rối theo phương
ngang cân bàne với nhớt xoáy theo phương thẳng đứne;

t e l - M S 7 ) - ( 5 7 - 5 7 ) - I f c ) - (* !
(b.3)
(b.4)
với ^ nhiệt độ thế vị, là nhiệt độ mà khối nước có được nếu nó thay đổi đoạn nhiệt (hoàn
toàn cách nhiệt với môi trường) từ mực khởi điểm đến mực áp suất chuẩn. Các đại lượng
D-* liên quan đến sự khuếch tán cùa động năng rối và tổc độ hao tán nãng lượng rối
(b.5)
(b.6)
Trong các công thức trên, Prĩ là số Prandtl rối, C1«- là các hàng số bán thực nghiệm
liên quan đến hao tán năng luợng
cu = 1.44
c,. =
C = 1.92
11.0 neu P-
[0 2 r.eu P;
> 0
< 0
(b.7)
(b.8)
(b.9)
9
Pr, - 0.35
(b.lO)
Nhớt tuyệt đối hoặc là giá trị cho trước ứng với từng môi trường nước, hoặc là hàm cùa
nhiệt độ được tính qua công thức Seeton. Trong khi đó, nhớt rối ự Ị được định nghĩa mô
hỉnh k — £ thông qua động năng rối k và tốc độ hao tán năng lượng rối f
Các đại lượng °k-ơt là các số Prandlt rối liên quan đến việc tham số hoá khép kín rối
k ~ f , còn S. trong biểu thức trên là các hàng số cùa mô hỉnh
Như vậy ta đã đưa thêm vào mô hình động năng rối X và tốc độ hao tán năng lượng rối f .
Mối quan hệ giữa hệ phương trình chuyển động và hệ phương trình khép kín rối A' - £■

ỏy' oz\ CZ •
ỞS ds ds ỞS d ị ds*. d ị 05, d ị ôs,
(c.l)
(C.2)
trong đó Đị ,Ds là hệ số khuếch tán của nhiệt độ và độ mặn. Hệ số khuếch tán độ mặn
được xét trong phần khuếch tán chất tan. Còn mối quan hệ giữa hệ số khuếch tán nhiệt Dỹ
với hệ số khuếch tán rối ữ .'.ỉ được thể hiện qua số Prandtl Pr, một đại lượng không thứ
nguyên. Trong các bài toán truyền nhiệt, số Prantdl cho biết độ dày tương đối giữa lớp biên
cùa vận tốc và lớp biên nhiệt. Tương tự như vậy trong các bài bài toán truyền chất, số
không thứ nguyên Schmidt mô tả quan hệ giữa các hệ khuếch tán vật chất với
khuếch tán rối &.'.Ị
° . - % M
D - = ■“ (c.4)
s Sc
Với mô hình này, chúng tôi sử dụng các hệ số Prandtl cho nhiệt độ và Schmidt cho độ mặn
của nước biển như sau
Pr = 0.72 (C.5)
Sc = 1.0 (c.6)
Phương trình trạng thái cho biết khối lượng riêng được tính qua nhiệt độ. độ mặn và
áp suất. Cụ thể, khối lượng riêng cùa nước biền được tính qua công thức Jackett-
McDougallz, là biểu thức chứa 26 số hạng. Trong tính toán số, chúng tôi sử dụng phương
trình trạng thái xấp xi của của B n ’don D., Sun s., Bỉeck R. (1999)
Pi9 5 pì = Cji.p) - C-(PÌ<? - C3(P)5 - c ÍP)6: - CẠP)Sê
' . (C.7)
- CẠP)Ôi - C-(P)S$-
với c. (p) (.' = • -7 ) là các hàm số bậc hai cùa áp suất
11
c. (p) = ct. — &.P — y.p~ (i = 12 7) (C.8)
Các hằng số C(. fi Yi (■ = 1-2. • -7) đối với các khoảng biển thiên nhiệt độ thế vị độ mặn
5 và áp suất ^

3 .2 6 36 3 X 1CT11
12
d) Phương trình khuếch tán chất tan
Phương trình khuếch tán chất tan c trong nước biền sau khi lấy trung bình
trong đó Dc là hệ sổ khuếch tán cùa chất tan, được tính qua hệ sô khuêch tán rối qua
công thức
Dc = (d.2)
c '■ Sc
Qr = /.c -ỹ.1C-ỳ.:C: (d.3)
với /. là tham số khuếch tán ứng với từng loại chất tan xác định. Qc là đại lượng đặc trưng
cho sự hình thành và phân huỹ chất tan, cũng như các tương tác hoá-lý-sinh khác,
/ ữ(.v.y r;f )1/ 1(.r y t), À:[ V y.z. t) là các hàm số liên tục liên quan đến sự phân huỷ,
tạo thành hoặc tương tác của chất tan.
g) Các điều kiện biên
Điều kiện biên bao gồm các điều kiện cho vận tốc - l phương trình mặt
thoáng 'ì. nhiệt độ độ mặn 5, chất tan c, động năng rối * và tốc độ hao tán năng lượng
rối f .
i)Trẽn mặt thoáng
Gió trẽn mặt thoáng và ma sát tác động làm thav đổi vận tốc dòng chảy. Gọi là
vận tốc gió trèn mặt biển với hai thành phần theo hai phương ngang tương ứne là "V W; ,
còn ~i là úng suất do gió bề mặt
= ữ-Yt W.: (g.l)
W. = (8-2)
Khi đó, ứng suất cùa dòng chày tại bề mặt tỷ lệ thuận với úng suất do gió thỏna qua hệ số
tỷ lệ •}. còn thành phần vận tốc theo phương thẳng đứng tại bề mặt được tính qua
13
à v I
O- ' r/sfi
đ r
PoD'ỈZ

T:
= 0
(g.n:
0.07 K- D
(g.12)
Trên mặt thoáng diễn ra quá trình trao đổi nhiệt độ cùa nước với lớp khône khí trên biên,
và xem không có sự thay đôi nông độ chât tan do bôc hơi tại bề mặt nước biên
dỡ
0=D-dz
Hị,
= _ Ị t? - e !
c_ -
(g.l3j
14
ds
P'ds t :
= 0
(g.!4)
Với ^5 là hệ số truyền nhiệt bề mặt, Cp là nhiệt dung riêng cùa nước, 8ỉ là nhiệt độ mặt
nước, là nhiệt độ cân bằng. Chất tan cũng không thay đồi nổng độ trên mặt thoáng do
bốc hơi hay biến đổi
ỖC
= 0
(g-15)
ii)Tạĩ đáy biển
Tại đáy biển, vận tốc dòng chảy bị suy giảm do ma sát đáy. Gọi 11'b ỉà vận tốc tại
lớp sát đáy biển với hai thành phần theo hai phương ngang tương ứng là , còn rt là
ứng suất ma sát đáy, chúng được tính qua
r; = PoYb u ;
\yt « :(u ;‘ ): - (u ;v j;

(g.25)
1
K o:
(8-26)
Nếu có dòng chày đồ vào miền tính trên đường bờ, thì có thể xem điều kiện biên cho và s
liên quan đến điều kiện đầu ra cùa bài toán dòng chảy trong kênh
= 0 004 V: (g.27)
_ c. k'-1 ,
0.09 ả K
(8-28)
trong đó U* là vận tốc dòng, a là chiều rộng dòng chảy đổ vào miền tính. Đáy biển coi như
cách nhiệt và không có quá trình xâm nhập mặn vào đất
p.Dg
de
dz
= c
(g.29)
0:DÌ
<?5|
= 0
(g.30)
Cũng như chất tan không thấm vào đất
= 0
(g 31)
16
Như vậy với phương pháp tuyến n bất kỳ, ta coi không có dòng nhiệt độ, độ mặn và nồng
độ chất tan trên biên rắn
de
dn
= 0

- OAi HO C Q UỘ C G ia h a n ộ i
Ị ftuNG T A M THQ N G Tin Thư viền
17
ở\’ ỖV
— —
c —— —

0
àt dy
(g-41)
(g.42)
(g.43)
Ta cũng có thể tính vận tốc qua dao động của mặt thoáng. Phương trình dao động
cùa mặt thoáng mặt thoáng hoặc cho trước bời các giá trị quan trắc, hoặc tính toán thông
qua tổng hợp của các dao động thuỷ động lực thành phần. Ta cũng có thể coi không có thay
đổi gradient phương trình mặt thoáng, hay sử dụng điều kiện phát xạ. Như vậy các điều
kiện biên cho phương trình mặt thoáng đối với dòng chày vào hoặc dòng chày ra bao gồm
với /: là hệ số suy giảm, là hệ số điều hoà biên độ, Q. là tần số góc, 0 ’c " u); là hệ số
điều hoà pha, 9. là góc lệch pha của dao động, c = ỹK là vận tốc truyền sóng trong
vùng nước nông. Ở đây °, đỏng vai trò v hoặc y đối với biên vuông góc với trục y hoặc trục
X tương ứng. Sau khi tìm được phương trinh mặt thoáng, chúng ta thay vào phương trinh
chuvền động với một số già thiết đơn giàn để xác định giá trị vận tốc.
Các điều kiện biên lòng cho động năng rối phức tạp hơn về mặt vật lý. Có nhiều
cách thiết lập các điều kiện biên này trong tính toán mô hình. Đối với dòng chảy vào động
năniĩ rối 'x và hao tán năng lưọng rối f chúng được cho trước trên biên
7
(g-44)

àĩ àq
(g-54)
ds ds
dt dạ
(g-55)
CI = ct
(g-56)
Trong trường hợp dòng chày ra khỏi miên tính
ỏc dc
dĩ c dq
. °
(g.57)
19
1.2 Hê toa đô sigma
• • 9 o
Những mô hình thuỷ nhiệt động lực học biển ba chiều dưới dạng Descartes không
thật sự tốt cho các tính toán ứng dụng cụ thể, đặc biệt là đối với các lớp bài toán hoàn lưu
vùng ven biển. Nguyên nhân chù yếu ở chỗ: độ sâu của đáy biển biến thay đổi từ khoảng
vài mét ở vùng nước nông gần bờ cho đến hàng kilômét đối với những vùng xa bờ. dẫn tới
việc tính toán đồng thời các đặc trưng động lực học của cột nước tại các miền có độ sâu
khác nhau gặp nhiều khó khăn. Do đó, xử lý điều kiện biên tại từng toạ độ trở nên hết sức
phức tạp.
Hệ toạ độ sigma đã hạn chế được những nhược điểm kể trên. Trong hệ sigma. chiều
cao cột nước là như nhau đối với toàn miền tính toán, với giá trị ơ nằm trong khoảng
[0; -1 ], Điều này tạo nên sự đồng nhất các điểm trên biên mặt nước và biên đáy với vị trí
tính toán. Việc giải số với các phương trình trong hệ sigma sẽ dễ dàng hơn do các điểm
biên này sẽ nằm ngay trên lưới tính.
20
Hình 4. Mô hình trong hệ toạ độ Descartes
Xét phép biến đổi toạ độ từ hệ toạ độ (v,y ,r, f) sang hệ toạ độ ụ y \ <7 f ) với

Như vậy đạo hàm toàn phần của - là hàm
dĩ Ở-: ỗĩ dĩ
T - - V r r - v - r - w - r
dt dx ỞY d:ơx ơy
_
1

(d- L
" D ' T r
ờ ị 3'ì \ ỏ ị ỏ - 1 ỏ Ị ỏ - \
1 2
dĩ' D d-rưD d:Y'D dT
(1.5)
ởx ■
ớy'
ỞC7
đ:
dx •
g
dr'
dí'
- 1
<d; ở-y
dơ'\ D dơ
(1.6)
. Công thức (1.5) là kết quả của biến đổi chính xác, còn công thức (1.6) là kết quả của biến
đổi xấp xỉ. Trường hợp riêng của công thức (1.6) là xấp xỉ cùa Mellor với được thav
bằng H. Ở đây chúng tôi biến đổi dưới dạng tổng quát và giữ nguyên trong tính toán. Vận
tốc theo phương sigma trong hệ toạ độ mới được định nghĩa
/ ỞD d'ì \ / dD ỡ'?> Ị ỖD dì] .

ở: d.x ởy dơ
. Trong hệ toạ độ sigma, phương trình của áp suất thu được bằng cách thay đạo hàm theo r
bởi đạo hàm theo ơ, sau đó lay tích phân từ mặt thoáng đến vị trí đang xét
p = - ỹữ: n - gD ị pảơ (1.11)
Áp suất cùa không khí trên mặt thoáng thu được từ việc phân tích và đồng bộ hoá dữ
liệu thực đo. Đối với mô hình hoàn lưu kết hợp cùa cà đại dương và khí quyền, áp suất
không khí có trên mặt thoáng là ẩn số và có thể tính được. Trong trường hợp không có dừ
liệu thực đo, ta có thề xem đại lượng này là hằng số. Đe đi tìm phương trình chuyền động
theo phương ngang trong hệ toạ độ sigma, ta cân phải đi tìm các đạo hàm riêng của áp suất
ỞP ở>5 d ị'c ÕD ởp
£ * JỊ*>
dP
d\
3 dif d ị'c dD ị'z dp
- = p.g -1- gD—
Ị
páa — g —— G ~r~ d
dy ơy I ởvJ. dơ
’ c .
da (1.12b)
. Sử dụng biến đồi (3.8) sang (3.9) với - lần lượt là vận tốc ~ ^ . kết hợp với việc đưa vào
các đại lượng lực Coriolis và đạo hàm riêng cùa áp suât. ta thu được phương trình chuvển
động trong hệ toạ độ sigma
dUD ỞU:D dUYD ỞUt-o di} qD~ d gDỞD f c dữ
— = -gD— — — peer

— a -—ào

1.4 Hệ phương trình khép kín rối
Hệ phương trình khép kín rối k - E trong hệ toạ độ sigma bao gồm phương trình
cho động năng rối k và tốc độ hao tán năng lượng rối Như đã trình bày trong chương
trước, hệ phương khép kín rối 'K - £ được xây dựng trên cơ sở tham số hoá các đại lượng
mạch động chưa biết bằng biểu thức của động năng rối X và tốc độ hao tán năng lượng rối
f rồi đưa về các phương trình năng lượng rối. Các giả thiết cơ bản để thiết lập hệ phương
trình trong mô hình bao gồm: giả thiết về tính đẳng hướng địa phương của chuyển động
rối: nhớt rối theo phương ngang cân bằng với nhớt xoáy theo phương thẳng đứng; giả thiết
Reynolds và giả thiết Boussinessq. Sử dụng biến đoi (1.8) sang (1.9) với - lần lượt là vận
tốc K ĩ, kết hợp với việc đưa vào thành phần nổi của động năng rối và thành phần úng suất
của động năng rối, ta thu được hai phương trình khép kín rối trong hệ toạ độ sigma
dkD dkUD dkVD dkoj
= (ppr - DPS) - Ds
dĩ ởx dy dơ
ỏ . dk i d ị _
- — \DD,— \DL: , ,
dx dx • dv ở d o D dơ
dỉD dsUD ởs\’D ỞSíj _ _ ĩ ĩ :
-JT - = Cu (DPr - CZiDPs) - - cĩtD -
ơt ox ơ) Ơơ K K
(1.16)
(1.17)
25