Xây dựng một sơ đồ tính truyền chất 1 D trên cơ sở hàm spline bậc 2

1 XÂY DỰNG MỘT SƠ ĐỒ TÍNH TRUYỀN CHẤT 1 D TRÊN CƠ SỞ HÀM SPLINE BẬC 2 GS. TSKH NGUYỄN ÂN NIÊN (1), THS. NCS BÙI VIỆT HƯNG (2) (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 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 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 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 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 đó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 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 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ơ đồ hiện dùng. I . MỞ ĐẦU 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ó dạng.   01                 ccc A q x c AD xAx c v t c q  (1) Trong đó : v – vận tốc trung bình mặt cắt A – diện tích mặt cắt ướt phụ thuộc mực nước Z tại mặt cắt đó. D – hệ số khuếch tán vật lý ( phân tử + rối) q – lưu lượng bổ sung ngang trên một đơn vị chiếu dài dọc dòng chảy cq – nồng độ chất dòng bổ sung ngang  - cường độ ( thể tích ) tăng / giảm nồng độ chất  - 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 các điểm mặt cắt.   A dAcu vAc 1  (2) Để 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 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 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 bằng các sơ đồ quen biết như VRSAP, TLUC, KOD01, MIKE 11 ,v,v. [2] Để giải phương trình truyền chất (1) ta phân rã thành phương trình tải :   0      ccc A q x c v t c q  (3) và phương trình khuếch tán 0 1              x c DA xAt c (4) 2 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 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 nghiệm cuối cùng. 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 trung vào việc giải phương trình tải (3). Phương trình này có đặc trưng : C =  v (5) 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 nguồn tập trung đổ vào đoạn dòng chảy) ta có : 0      c dt dc c x c v t c C   (6) Trong đó, đạo hàm toàn phần là trên đường đặc trưng v dt dx C   (7) Nghiệm của phương trình (6) trên đường đặc trưng là :    t Ot C dt oecc   (8) Trong trường hợp  = 0 thì nghiệm (8) là : constc C   (9) Để 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à 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 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 dùng hàm Largrange, tóm tắt như sau : (xem hình 1 ) [3]. L Lj 1 t+t N Lj-1 Lj+1 C t 0 x j-1 M x j j+1 j-1 j j+1 Hình 1 Để tìm c*j = cN ta vẽ đường đặc trưng đi tới điểm N và theo (9), ta có (khi =0) :       1111   jjjjjjj cxLcxLcxLc  (10) 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 : + 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 tính chất đặc trưng của phương trình tải. + 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 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]. Trong trường hợp  ≠ 0 ta có : 3   tjjjjjjtMj ecLcLcLecc    1111 (11) Trong đó  là trị số  trung bình trên đoạn MN 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, để khắc phục các nhược điểm trên. Đó là lý do đề xuất một sơ đồ giải mới. II. SƠ ĐỒ GIẢI PHƯƠNG TRÌNH TẢI SỬ DỤNG HÀM NỘI SUY SPLINE BẬC 2 Để đả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 đặ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) : E N x G H S t+t x’j-1 x’j Sj Sj-1 t j-1 xj-1 j xj j+1 j-1 j Hình 2 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ý hiệu :    vvvtvNGx tvvxEGx j jjjj    2 1 ; 111  v’ : lưu tốc mặt cắt tại lớp t+t Ta có các hàm Spline sau (x’ là khoảng cách đến điểm G) và : 0 ≤ x’ ≤ x’j-1 1 11 1 1 1 1 1                        jj jj j j j SS x x x x k x x S  (12) 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 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] :       j j jj jGHjEG EGj x x xx xccxcc cck            1 1 1 1 (13) 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 ). Ngược lại, ta lấy kj-1 như sau :        1******1   jGHjEGEGj xccxccccsignk (14) Như vậy ta luôn có : -1 ≤ k ≤ 1 (15) và điều kiện này đảm bảo trị số S nằm giữa 0 và 1. K=-1 K=0 K=1 K=1 K=0 K=-1 4 Công thức (13) có thể biến đổi như sau : Nếu ta ký hiệu * Hc là ngoại suy tuyến tính cG * và cE *, ta có :   1 ****     j j EGGH x x cccc và ta ký hiệu : *** HHH ccc  thì để tính kj-1 có thể dùng biểu thức :    jjj j HEGj xxx x ccck       1 2 1* 1 (16) 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 cG * - cE * ), và trường hợp ngược lại, thì ta lấy:   ***1 EGHj cccsignk  (17) Sau khi có trị số kj-1 , ta tính cj * = cN *     **1* GjEjj cxScxSc    (18) Trong (12) với tvx j , nếu ta nhân trên dưới tỷ số 1   jx x với A’(diện tích mặt cắt ướt ) ta sẽ có thể thay tỷ số đó bằng : 111             j j j j j w tQ x tv A A x x (19) trong đó :   tQQww jjjj   111 (20) 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 gian t. 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 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 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* = c te  với    2 1 . III. TÍNH TOÁN NỒNG ĐỘ CHẤT TẠI NÚT SÔNG. 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 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ể (xem hình 3) 5 Hình 3 – nút sông Các mặt cắt đều lấy sát nút sông và ta có :     raravaovaoravao QccQQQ ; (21) 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 phẩy để chỉ đặc trưng ở lớp thời gian tính toán t+t) :       n K jjjjj n K jjjjj KKKKKKKKK QdQdccQdQd 11 (22) Trong đó, djK là thông số chỉ hướng chảy vào nút. djK = 1 nếu chiều dương QjK trùng hướng vào nút djK = -1 trong trường hợp ngược lại 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 nồng độ chất (chẳng hạn độ mặn) của các nhánh ra không đều nhau. 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 ( 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 1 Kj k , 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. Trước tiên ta phân biệt nhánh chảy vào và ra theo dấu của KK jj Qd và sắp xếp lại như sau : - Nhánh vào đánh số ji với i  (1,2,.,.,n) có n1 nhánh - Nhánh ra đánh số jo với o  (1,2,.,.,n) có n – n1 nhánh Ta ký hiệu i Oi Oi j jj jj Q Q    , , (23) Đó là tỷ lệ lưu lượng nhánh ji đổ vào nhánh ra jO. Ta ký hiệu          n K jjjj jjjj jj jj j KKKK iiii ii ii i QdQd QdQd Qd Qd 1  (24) Đó 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ệ số ji ≠ 0 . Bây giờ, ta có thể tính cjo * như sau :  i iOiiO j jjjjj cc * , *  (25) Việc tìm giá trị Oi jj ,  như thế nào chúng tôi sẽ có một bài viết riêng. Trường hợp sơ đồ 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ố ,) với việc có nội suy bằng hàm spline trong (25) – xem thêm [5] ) IV. KẾT LUẬN 6 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 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 vẫn dùng và có thể phát triển ứng dụng rộng rãi. TÀI LIỆU THAM KHẢO [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 KH và Cong nghệ – Viện KH Thủy Lợi Miền Nam – NXB Nông Nghiệp – 1999. [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 Cửu Long – Tuyển tập KH và CN – Viện KH Thủy Lợi Miền Nam – NXB Nông Nghiệ – 2001. [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 mặn cho hệ thống sông kênh – Luận án Tiến sỹ Cơ học – 1987. [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. 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 Nông nghiệp – 2003. [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. 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 Nông Nghiệp – 2003. [6] Nguyễn An Niên, Tăng Đức Thắng – Computation of mass transmission by a forced – mixed model ( One-dimensional problem). Proc. of International colloquium in Mechanics of Solids. Fluids, Structures and Interactions – Nha Trang – 2000. SUMMARY Creation of an one-dimensional mass transmission mathematical scheme base on second degree Spline Prof. Dr. of Sci. Nguyen An Nien - Southern Institute for Water Resources Research MEng. Bui Viet Hung - Hydraulic Construction Consultant Company 2 When solving one-dimensional mass transmission problem, the finite difference method has more priority. After solving one – dimensional hydraulic problem for receiving hydraulic characteristics (discharge, water level ,.,.,) we resolute mass transmission equation by separating into transport equation and diffussion one. Solutions of transport problem, as well known, has dominant role and solutions of diffusion one based on solution of transport one only adjust these solutions. We present a scheme for solving transport equation based on secondary degree spline, that can overcome defects of currently used computational schema. Người phản biện: PGS.TS. Phó Đức Anh