XÂY DỰNG MỘT SƠ ĐỒ TÍNH TRUYỀN CHẤT 1 D<br />
TRÊN CƠ SỞ HÀM SPLINE BẬC 2<br />
<br />
GS. TSKH NGUYỄN ÂN NIÊN (1), THS. NCS BÙI VIỆT HƯNG (2)<br />
(1) Viện Khoa học Thuỷ lợi Miền Nam (2) Công ty Tư vấn xây dựng Thủy lợi 2<br />
<br />
Tóm tắt : Trong việc giải bài toán truyền chất một chiều, phương pháp sai phân hữu<br />
hạn có nhiều lợi thế. Sau việc giải bài toán thủy lực độc lập để xác định đặc trưng<br />
dòng chảy ( lưu lượng, mực nước,.,) ta giải phương trình truyền chất bằng cách phân<br />
rã thành phương trình tải và phương trình khuếch tán. Nghiệm của phương trình tải<br />
đóng vai trò chính và phương trình khuếch tán được giải trên nền nghiệm này, nói<br />
chung chỉ vi chỉnh lại nghiệm đó để cho lời giải cuối cùng. Chúng tôi đề xuất việc giải<br />
phương trình tải dùng hàm spline bậc 2 để khắc phục một số nhược điểm của sơ đồ<br />
hiện dùng.<br />
<br />
I . MỞ ĐẦU<br />
<br />
Phương trình truyền chất một chiều với nồng độ c trên một đơn vị thể tích chất lỏng có<br />
dạng.<br />
c c 1 c q<br />
v AD c q c c 0 (1)<br />
t x A x x A<br />
Trong đó :<br />
v – vận tốc trung bình mặt cắt<br />
A – diện tích mặt cắt ướt phụ thuộc mực nước Z tại mặt cắt đó.<br />
D – hệ số khuếch tán vật lý ( phân tử + rối)<br />
q – lưu lượng bổ sung ngang trên một đơn vị chiếu dài dọc dòng chảy<br />
cq – nồng độ chất dòng bổ sung ngang<br />
- cường độ ( thể tích ) tăng / giảm nồng độ chất<br />
- hệ số mô men do sự phân bố không đều của nồng độ chất c’ và lưu tốc u’ tại<br />
các điểm mặt cắt.<br />
1<br />
uc dA (2)<br />
vAc A<br />
Để có tốc độ v, diện tích mặt cắt ướt A(x,z) và cả hệ số khuếch tán D (phụ thuộc vào<br />
các yếu tố thủy lực, chẳng hạn trong công thức của Fisher [1]) ta cần giải bài toán thủy lực<br />
trước để nhận được các đặc trưng mực nước z, lưu lượng Q tại các mặt cắt ở mọi lớp thời gian<br />
bằng các sơ đồ quen biết như VRSAP, TLUC, KOD01, MIKE 11 ,v,v. [2]<br />
<br />
Để giải phương trình truyền chất (1) ta phân rã thành phương trình tải :<br />
c c q<br />
v c q c c 0 (3)<br />
t x A<br />
và phương trình khuếch tán<br />
c 1 c <br />
DA 0 (4)<br />
t A x x <br />
<br />
<br />
<br />
<br />
1<br />
Nghiệm của phương trình tải đóng vai trò chính trong nghiệm của bài toán và giải<br />
phương trình khuếch tán (4) trên nền nghiệm phương trình tải chỉ tu chỉnh nghiệm này để cho<br />
nghiệm cuối cùng.<br />
<br />
Vì việc giải phương trình khuếch tán không có gì phải bàn thêm nên, chúng tôi chỉ tập<br />
trung vào việc giải phương trình tải (3). Phương trình này có đặc trưng :<br />
C = v (5)<br />
Không làm giảm tích chất tổng quát ta xem q=0 (tập trung lưu lượng bên thành một<br />
nguồn tập trung đổ vào đoạn dòng chảy) ta có :<br />
c c dc<br />
v c c0 (6)<br />
t x dt C<br />
Trong đó, đạo hàm toàn phần là trên đường đặc trưng<br />
dx<br />
C v (7)<br />
dt<br />
Nghiệm của phương trình (6) trên đường đặc trưng là :<br />
t<br />
<br />
<br />
dt<br />
tO<br />
c co e (8)<br />
C<br />
<br />
Trong trường hợp = 0 thì nghiệm (8) là :<br />
c const (9)<br />
C<br />
<br />
Để giải bằng sai phân ta chia dòng chảy thành từng đoạn bởi các mặt cắt có toạ độ xj và<br />
cách nhau đoạn xj . Để giải ra nghiệm c*j của phương trình tải khử được khuếch tán số của<br />
phép nội suy tuyến tính, một phương pháp vẫn hay dùng được gọi là phương pháp nội suy<br />
dùng hàm Largrange, tóm tắt như sau : (xem hình 1 ) [3].<br />
L Lj<br />
1<br />
t+t N Lj-1 Lj+1<br />
C<br />
<br />
<br />
t 0 x<br />
j-1 M x j j+1 j-1 j j+1<br />
Hình 1<br />
<br />
Để tìm c*j = cN ta vẽ đường đặc trưng đi tới điểm N và theo (9), ta có (khi =0) :<br />
<br />
c j L j 1 x c j 1 L j x c j L j 1 x c j 1 (10)<br />
Biến thiên của hàm nội suy L như trên hình vẽ. Cách nội suy đó có hai nhược điểm :<br />
<br />
+ Huy động trị số cj+1 đứng phía hạ lưu ( trong hình vẽ xem v>0) là không đúng với<br />
tính chất đặc trưng của phương trình tải.<br />
+ Các hàm L có thể có trị số âm (Lj-1, Lj+1) và lớn hơn 1 (Lj), nên có thể gặp nhiều tình<br />
huống trị số c*j = cM vượt ra ngoài giới hạn cj-1, cj, cj+1 và thậm chí có trị số âm [1].<br />
<br />
Trong trường hợp ≠ 0 ta có :<br />
<br />
<br />
2<br />
c j c M e t L j 1 c j 1 L j c j L j 1 c j 1 e t (11)<br />
Trong đó là trị số trung bình trên đoạn MN<br />
<br />
Như vậy, rõ ràng cần phải tìm phương pháp giải phương trình tải một cách khác, để<br />
khắc phục các nhược điểm trên. Đó là lý do đề xuất một sơ đồ giải mới.<br />
<br />
II. SƠ ĐỒ GIẢI PHƯƠNG TRÌNH TẢI SỬ DỤNG HÀM NỘI SUY SPLINE BẬC 2<br />
<br />
Để đảm bảo chỉ nội suy trị số nồng độ giữa mặt cắt thượng lưu (ảnh hưởng theo đường<br />
đặc trưng truyền xuống) và mặt cắt đang xét, ta sử dụng sơ đồ động sau (hình 2) :<br />
E N x G H S<br />
t+t x’j-1 x’j K=-1<br />
K=0<br />
K=1 Sj<br />
Sj-1<br />
K=1<br />
t K=-1<br />
K=0<br />
j-1 xj-1 j xj j+1 j-1 j<br />
Hình 2<br />
<br />
Ta dựng các đường đặc trưng qua các mặt cắt j-1, j và cắt lớp t+t tại E và G . Ta ký<br />
hiệu :<br />
<br />
xj 1 EG x j 1 v j v j 1 t <br />
1<br />
x NG v j t ; v v v<br />
2<br />
v’ : lưu tốc mặt cắt tại lớp t+t<br />
<br />
Ta có các hàm Spline sau (x’ là khoảng cách đến điểm G) và :<br />
0 ≤ x’ ≤ x’j-1<br />
x x x <br />
S j 1 k j 1 1<br />
x j 1 x j 1 x j 1 <br />
S j 1 S j 1<br />
(12)<br />
với kj-1 hệ số liên quan đến đạo hàm bậc 2 của hàm c, có thể huy động một điểm phía thượng<br />
lưu hoặc hạ lưu, ở đây, ta lấy điểm H ở hình 2. Ta có biểu thức tính kj-1 như sau [4] :<br />
<br />
<br />
k j 1 cG c E <br />
x j 1 x j<br />
<br />
x j<br />
<br />
cG c E x j c H c G x j 1 x j 1<br />
(13)<br />
<br />
Nếu | c*G – c*E | > | vế phải (13) | thì tính kj-1 thông thường ( chia vế phải cho c*G – c*E ).<br />
Ngược lại, ta lấy kj-1 như sau :<br />
<br />
k j 1 sign cG* c E* c G* c *E x j c H* c G* x j 1 <br />
(14)<br />
Như vậy ta luôn có : -1 ≤ k ≤ 1 (15)<br />
và điều kiện này đảm bảo trị số S nằm giữa 0 và 1.<br />
<br />
<br />
<br />
3<br />
*<br />
Công thức (13) có thể biến đổi như sau : Nếu ta ký hiệu c H là ngoại suy tuyến tính cG*<br />
và cE*, ta có :<br />
* x j<br />
c H cG* cG* c E* <br />
x j 1<br />
và ta ký hiệu : c H* c *H c H* thì để tính kj-1 có thể dùng biểu thức :<br />
x j 21<br />
<br />
*<br />
k j 1 cG c E c H<br />
x j x j 1 x j <br />
(16)<br />
<br />
tương tự, nếu | cG* - cE* | > | vế phải (16) | thì khi tính kj-1 như thông thường (chia vế phải cho<br />
cG* - cE* ), và trường hợp ngược lại, thì ta lấy:<br />
*<br />
<br />
k j 1 sign c H c G* c E* (17)<br />
Sau khi có trị số kj-1 , ta tính cj* = cN*<br />
<br />
c *j S j 1 x c *E S j x c G* (18)<br />
<br />
x <br />
Trong (12) với x v j t , nếu ta nhân trên dưới tỷ số với A’(diện tích mặt cắt<br />
x j 1<br />
ướt ) ta sẽ có thể thay tỷ số đó bằng :<br />
x A v j t Q j t<br />
(19)<br />
x j 1 A x j 1 wj 1<br />
trong đó :<br />
<br />
wj 1 w j 1 Q j Q j 1 t (20)<br />
là thể tích chứa nước trong khối EG còn wj-1 là thể tích nước giữa 2 mặt cắt j-1 và j ở lớp thời<br />
gian t.<br />
<br />
Cần lưu ý là khi giới hạn kj-1 = ± 1 là trường hợp | cj – cj-1 | khá nhỏ nên sai số phạm<br />
phải trong tính hệ số này sẽ có bậc của sai số tính toán, hơn nửa hệ số kj-1 chỉ nằm trong phần<br />
chỉnh bậc 2 của phép nội suy tuyến tính. Các trị số c* tại E, G, H có thể ứng với ≠ 0, tức là c*<br />
= c e t với 1 .<br />
2<br />
<br />
<br />
III. TÍNH TOÁN NỒNG ĐỘ CHẤT TẠI NÚT SÔNG.<br />
<br />
Nút sông là chỗ các nhánh sông giao nhau. Trong các sơ đồ tính truyền thống đều xem<br />
có sự xáo trộn hoàn toàn ở nút sông, do đó, nồng độ của các nhánh ra là bằng nhau. Cụ thể<br />
(xem hình 3)<br />
<br />
<br />
<br />
<br />
4<br />
Hình 3 – nút sông<br />
<br />
Các mặt cắt đều lấy sát nút sông và ta có :<br />
Qvao Qra ; Qvao c vao Qra<br />
c ra (21)<br />
<br />
Nếu ta ký hiệu các mặt cắt liên quan đến nút j là jK ( K=1,2,3,.,.,n) thì ta có thể viết (dấu<br />
phẩy để chỉ đặc trưng ở lớp thời gian tính toán t+t) :<br />
n n<br />
<br />
d<br />
K 1<br />
<br />
K 1<br />
<br />
jK Q jK d jK Q j K c j K c j d j K Q jK d jK Q jK<br />
(22)<br />
<br />
Trong đó, djK là thông số chỉ hướng chảy vào nút.<br />
djK = 1 nếu chiều dương QjK trùng hướng vào nút<br />
djK = -1 trong trường hợp ngược lại<br />
<br />
Giả thiết xáo trộn hoàn toàn tại nút sông chưa phù hợp với thực tế và số liệu đo đạc các<br />
nồng độ chất (chẳng hạn độ mặn) của các nhánh ra không đều nhau.<br />
<br />
Cũng như các sơ đồ truyền thống ta giải bài toán tải tìm cjK* cho các nhánh chảy đến (<br />
chỉ khác là huy động nồng độ cD* của các mút đường đặc trưng đứng trước nữa để tính k jK 1 ,<br />
xem [4]) còn tuỳ thuộc vào việc phân chia nước của từng nhánh vào đối với nhánh ra.<br />
<br />
Trước tiên ta phân biệt nhánh chảy vào và ra theo dấu của d jK Q jK và sắp xếp lại như<br />
sau :<br />
- Nhánh vào đánh số ji với i (1,2,.,.,n) có n1 nhánh<br />
- Nhánh ra đánh số jo với o (1,2,.,.,n) có n – n1 nhánh<br />
Q ji , jO<br />
Ta ký hiệu j ,j (23)<br />
i O<br />
Q ji<br />
Đó là tỷ lệ lưu lượng nhánh ji đổ vào nhánh ra jO.<br />
d ji Q ji d ji Q ji d ji Q ji<br />
Ta ký hiệu ji n (24)<br />
d ji Q ji d jK Q jK d jK Q jK<br />
<br />
K 1<br />
Đó là tỷ lệ lưu lượng đổ vào của nhánh ji , và tất nhiên, chỉ những nhánh vào mới có hệ<br />
số ji ≠ 0 .<br />
<br />
Bây giờ, ta có thể tính cjo* như sau :<br />
c *jO ji ji , jO c *ji (25)<br />
ji<br />
<br />
Việc tìm giá trị ji , jO như thế nào chúng tôi sẽ có một bài viết riêng. Trường hợp sơ đồ<br />
tính có mặt cắt không sát các nút ( ví dụ sơ đồ KODWQ) thì xử lý cũng gần giống ( các hệ số<br />
,) với việc có nội suy bằng hàm spline trong (25) – xem thêm [5] )<br />
<br />
IV. KẾT LUẬN<br />
<br />
<br />
<br />
5<br />
Sơ đồ giải bài toán tải bằng hàm Spline đưa ra một hướng mới, hoàn chỉnh sơ đồ xáo<br />
trộn cưỡng bức mà trước kia chúng tôi đã đề ra [6] nhằm khắc phục nhược điểm của sơ đồ tính<br />
vẫn dùng và có thể phát triển ứng dụng rộng rãi.<br />
<br />
TÀI LIỆU THAM KHẢO<br />
<br />
[1] Nguyễn An Niên – Một số cải tiến sơ đồ tính truyền chất ( bài toán một chiếu) – Tuyển tập<br />
KH và Cong nghệ – Viện KH Thủy Lợi Miền Nam – NXB Nông Nghiệp – 1999.<br />
<br />
[2] Nguyễn An Niên – Phân tích các mô hình tính toán thủy lực sử dụng cho Đồng bằng sông<br />
Cửu Long – Tuyển tập KH và CN – Viện KH Thủy Lợi Miền Nam – NXB Nông Nghiệ –<br />
2001.<br />
<br />
[3] Nguyễn Tất Đắc – Mô hình toán học không dừng một chiều cho truyền triều và xâm nhập<br />
mặn cho hệ thống sông kênh – Luận án Tiến sỹ Cơ học – 1987.<br />
<br />
[4] Nguyễn An Niên, Bùi Việt Hưng – Cải tiến tiếp theo của sơ đồ tính truyền chất một chiều.<br />
Phần I – Sơ đồ tính nhánh sông. Tuyển tập KH và CN – Viện KH Thủy lợi Miền Nam – NXB<br />
Nông nghiệp – 2003.<br />
<br />
[5] Nguyễn An Niên, Bùi Việt Hưng – Cải tiến tiếp theo của sơ đồ tính truyền chất một chiều.<br />
Phần II – Sơ đồ tính mạng sông – Tuyển tập KH và CN – Viện KH Thủy lợi Miền Nam – NXB<br />
Nông Nghiệp – 2003.<br />
<br />
[6] Nguyễn An Niên, Tăng Đức Thắng – Computation of mass transmission by a forced –<br />
mixed model ( One-dimensional problem). Proc. of International colloquium in Mechanics of<br />
Solids. Fluids, Structures and Interactions – Nha Trang – 2000.<br />
<br />
SUMMARY<br />
<br />
Creation of an one-dimensional mass transmission mathematical<br />
scheme base on second degree Spline<br />
<br />
Prof. Dr. of Sci. Nguyen An Nien - Southern Institute for Water Resources Research<br />
MEng. Bui Viet Hung - Hydraulic Construction Consultant Company 2<br />
<br />
When solving one-dimensional mass transmission problem, the finite difference method<br />
has more priority. After solving one – dimensional hydraulic problem for receiving hydraulic<br />
characteristics (discharge, water level ,.,.,) we resolute mass transmission equation by<br />
separating into transport equation and diffussion one. Solutions of transport problem, as well<br />
known, has dominant role and solutions of diffusion one based on solution of transport one<br />
only adjust these solutions. We present a scheme for solving transport equation based on<br />
secondary degree spline, that can overcome defects of currently used computational schema.<br />
<br />
Người phản biện: PGS.TS. Phó Đức Anh<br />
<br />
<br />
<br />
6<br />