Mô hình số về vận tải bùn cát ba chiều trong sông - Trần Thị Mỹ Hồng

46 Trần Thị Mỹ Hồng, Lê Song Giang MÔ HÌNH SỐ VỀ VẬN TẢI BÙN CÁT BA CHIỀU TRONG SÔNG 3D NUMERICAL MODELING OF SEDIMENT TRANSPORT IN RIVERS Trần Thị Mỹ Hồng, Lê Song Giang Trường Đại học Bách khoa – Đại học Quốc gia Tp.HCM; tranthimyhong@hcmut.edu.vn, lsgiang@yahoo.com Tóm tắt - Bài báo trình bày một phương pháp tính toán dòng chảy và vận tải bùn cát trong sông, kênh hở trong không gian ba chiều. Các phương trình vi phân cơ bản để tính toán dòng chảy và vận tải bùn cát được viết trong hệ tọa độ biến đổi “sigma”. Các phương trình này được giải theo phương pháp thể tích hữu hạn trên lưới phi cấu trúc với các phần tử tứ giác. Phương pháp tính này được tính toán kiểm tra với bài toán vận tải bùn cát trong đoạn kênh cong được nghiên cứu thực nghiệm bởi Odgaard và Bergs (1988) và áp dụng thử để làm mô hình vận tải bùn cát tại khu vực Cù lao Phố trên sông Đồng Nai. Các kết quả tính toán cho thấy phương pháp này cho kết quả tính khá chính xác và có thể ứng dụng cho các tính toán thực tế. Abstract - The article presents a method of calculating three dimensional flow and sediment transport in rivers and open channels. The governing differential equations for flow and sediment transport are written in the "sigma" coordinate system. These equations are solved by the finite volume method on the unstructured grid with quadrilateral elements. The method is tested with the sediment transport in curved channel studied experimentally by Odgaard and Bergs (1988). The method is also applied to modeling sediment transport at Pho islet in DongNai river. Calculation results show that this method gives fairly accurate results and can be applied to practical calculations. Từ khóa - mô hình số 3D; vân tải bùn cát; lưới phi cấu trúc; sông Đồng Nai; phương pháp thể tích hữu hạn Key words - 3D numerical model; sediment trannsport; unstructured grid; DongNai river; finite volume method 1. Giới thiệu Tính toán dòng chảy và vận tải bùn cát là một trong các bài toán phổ biến trong thủy lực. Đối với các trường hợp đơn giản, người ta có thể dùng mô hình 1D hoặc 2D. Tuy nhiên đối với các bài toán phức tạp, nơi dòng chảy có cấu trúc 3 chiều với các dòng thứ cấp xuất hiện rõ rệt thì chỉ có mô hình 3D mới đủ độ tin cậy. Gần đây, một số mô hình 3D đã được công bố. Van Rijn [9] giải phương trình vận tải bùn cát lơ lửng 3D trong khi dòng chảy 3D nhận được từ lời giải 2D cùng với giả thiết phân bố theo quy luật logarithm trên phương thẳng đứng. Lin và Fanconer [4] giải phương trình Navier-Stokes trung bình hóa Reynolds và phương trình vận tải bùn cát lơ lửng bằng phương pháp sai phân hữu hạn trên lưới cố định theo phương thẳng đứng. Wu và ctg. [10] xây dựng một phương pháp tính dựa trên việc giải bằng phương pháp thể tích hữu hạn phương trình Navier-Stokes bình hóa Reynolds và phương trình vận tải bùn cát lơ lửng trung. Hiện nay cũng có một số phần mềm để tính dòng chảy và vận tải bùn cát ba chiều như DELFT3D, TELEMAC 3D tuy nhiên thuật toán của các phần mềm này chưa được tối ưu và chi phí bản quyền khá cao. Trong bào báo này, một phương pháp tính toán dòng chảy và vận tải bùn cát 3D khác được giới thiệu. Dòng chảy được giải từ phương trình 3 chiều với giả thiết thủy tĩnh. Các phương trình được giải theo phương pháp thể tích hữu hạn trên lưới phi cấu trúc và trục thẳng đứng biến đổi sigma. Phương pháp này đảm bảo được các tính chất bảo toàn và lưới tính toán có khả năng bám sát theo địa hình, phù hợp với các bài toán thực tiễn. Phương pháp đã được áp dụng tính toán thử nghiệm với bài toán kênh cong 180o được nghiên cứu bằng thực nghiệm bởi Odgaard và Bergs [5] và bài toán bồi xói khu vực Cù lao Phố trên sông Đồng Nai. 2. Mô hình toán 2.1. Phương trình cơ bản Dòng chảy được giải từ phương trình dòng chảy ba chiều với giả thiết thủy tĩnh. Trong hệ tọa độ biến đổi “sigma”, phương trình được viết : ( ) 0=   ++      q t D q (1)   FUUUqUq DgD D A qDA t V H +−=        −   +−+      (2) Còn bùn cát lơ lửng và diễn biến đáy được giải từ các phương trình:   ( ) 00 =        −−   +−+     C D D CwqCHDC t q V sH C q (3) ( ) ( ) ( )bbbbbC ED t bc t z P −=+   +   − q1 (4) Trong đó:  – mực nước;  Tyx uu ,=U – hai thành phần vận tốc của dòng chảy trên phương ngang; ω – thành phần vận tốc trên phương thẳng đứng “sigma”;   Uq Dqq Tyx == , và  Dq = ; D – độ sâu; σ – toán tử vi phân trên mặt phẳng “sigma”; F – vector ngoại lực; AH và AV – độ nhớt rối; C – nồng độ bùn cát lơ lửng trong cột nước; DCqC = ; ws0 – vận tốc lắng; zb – cao độ đáy; C – khối lượng riêng hạt bùn cát; P – hệ số rỗng; cb – nồng độ bùn cát trong tầng đáy; b – bề dày tầng đáy;  Tbybxb qq ,=q – vector lưu lượng đơn vị của bùn cát đáy; Db và Eb – suất lắng và xói của bùn cát tại đáy (ở cách đáy một khoảng là b); DH và DV – hệ số khuếch tán rối. Độ nhớt rối ngang, AH, được tính theo Smagorinsky [6]: 5,0 222 xy2 1 yx                 +   +        +        = vuvu yxCAH (5) Hệ số C nằm trong khoảng 0.01 - 0.5. Độ nhớt rối theo phương đứng, AV, được tính theo mô hình Prandtl- Kolmogorov (1942) [15, 16]: kLCAV = (6) Trong đó C' - hằng số mô hình, được xác định trong quá trình hiệu chỉnh mô hình; L – chiều dài xáo trộn; k – động năng rối. Các thông số này được lấy theo Davies và ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, VOL. 17, NO. 1.1, 2019 47 Gerritsen [1] như sau: ( ) ( ) ( ) ( )      −+= huhu c k sb 1 1 2 * 2 *  (7) ( ) hhDL −= 1 (8) Với c=0.09 – hằng số mô hình; =0.4 – hằng số Karman; su* - vận tốc ma sát trên mặt nước; bu* - dạng sửa chữa của vận tốc ma sát đáy bu* , trong đó ( )buu ** ,max= b * u ; *u - giá trị trung bình của các vận tốc ma sát tính từ vận tốc tại các mắt lưới trên đường thủy trực nếu biên dạng vận tốc là logarit; và ( ) Dzh −=  - độ sâu tương đối. Hệ số khuếch tán rối được lấy bằng với độ nhớt rối. Đối với bùn cát rời, suất lắng và xói bùn cát ở độ cao b, Db và Eb, được tính : (9) Trong đó cb và cb,e – nồng độ bùn cát và nồng độ bùn cát bão hoà ở mặt phân cách lớp đáy. Theo Van Rijn [7 - 9], qb và cb,e, được tính (10) 3.0 * 5.1 50 , 015.0 D T b d c Ceb = (11) Với: ( )b cr crT   = − (12) (13) Trong đó: d50 – đường kính hạt 50%; s – tỷ trọng hạt; τ và τcr – ứng suất tiếp đáy và ứng suất tiếp đáy ngưỡng (xác định từ đồ thị Shield). Đối với bùn cát kết dính, suất xói, Eb, và lắng, Db, được tính theo Hayter và Mehta [2] và Krone [3]:         −= 1    bbE (14)         −= d b sb CwD   10 (15) Trong đó:  - hệ số xói; b - ứng suất tiếp đáy;  - ứng suất tiếp đáy ngưỡng xói; d - ứng suất tiếp đáy ngưỡng bồi. 2.2. Điều kiện biên Trên biên hở diễn biến lưu lượng hoặc mực nước sẽ được áp đặt. Đối với nồng độ bùn cát lơ lửng, trong pha chảy vào, nồng độ bùn cát được áp đặt. Còn trong pha chảy ra, điều kiện thoát tự do được sử dụng: (16) Trên biên kín, điều kiện không thẩm thấu được sử dụng. Trên mặt thoáng (= 0), các điều kiện biên sau được sử dụng: 0 = (17) ( ) 000 ,,   yx V vu D A −=          (18) 00 =   +  C D D Cw Vs (19) Còn tại đáy (= -1): 0= (20) ( ) ( )vuvuCvu D A D V ,, 5.022 +=           (21) bb V s DE C D D Cw −=        +−  0 (22) Trong đó CD là hệ số ma sát đáy và được tính : ( ) 2 0ln        = z CD  (23) Với Δ – khoảng cách từ mắt lưới đầu tiên tới đáy; và z0 – thông số nhám. 2.3. Phương pháp giải Các phương trình (1) - (3) được giải bằng phương pháp thể tích hữu hạn theo sơ đồ được cải tiến từ sơ đồ được giới thiệu trong tài liệu [13]. Trong sơ đồ mới này, các thành phần của lưu lượng q và nồng độ bùn cát lơ lửng C được giải ẩn theo phương”sigma” để gia tăng tính ổn định của lời giải. Hình 1 giới thiệu lưới tính của mô hình. Mực nước được tính tại các nút lưới. Lưu lượng q được tính tại các cạnh của các phần tử, lưu lượng qσ được tính tại các nút lưới ở khoảng giữa các lớp còn nồng độ bùn cát lơ lửng được tính tại các nút lưới ở ngay từng lớp. a) Trên mặt bằng b) Trên phương thẳng đứng Hình 1. Sơ đồ lưới tính trên mặt bằng (a) và các lớp lưới tính theo phương thẳng đứng (b) Thứ tự giải các phương trình như sau: • Bước 1: Giải phương trình (2) tính lưu lượng q ở thời điểm n+1/2. • Bước 2: Giải phương trình (1) tính mực nước η ở thời điểm n+1 • Bước 3: Tính lưu lượng qσ ở thời điểm n+1/2. • Bước 4: Giải phương trình (3) tính nồng độ bùn cát lơ lửng C ở thời điểm n+1 • Bước 5: Giải phương trình (4) tính cao độ đáy zb ở thời điểm n+1 a. Bước 1 Phương trình (2) được tích phân trên thể tích kiểm soát của cạnh ict (Hình 1a) ở lớp k:   +  −+    kk S H S dSdDAdSd t     UqU q ( )ebbsbb ccwED ,0 −=− ( )  1.23.0* 5.1 50 5.0 1053.0 TDdgsq Cb −−=  ( ) 3/1 250* 1       −=  g sdD 0= nC 48 Trần Thị Mỹ Hồng, Lê Song Giang ( )  +−=          −   kk SS V dSdDgDdSd D A q       F U U (24) Sau khi sai phân số hạng đạo hàm theo thời gian, (24) đưa tới phương trình: ( ) ( ) ( ) knknknknkknknk VFlux t V =−++   − +++ −−+ rJJUqq 2/12/1 1 2/12/12/1 (25) Trong đó SV kk .= (26) 2/1 1 + +         −= k V k D A q D S   q qJ (27) ( )          −= L k Hnkk dl n DAFlux U UqU  (28) ( ) kk gD Fr +−=  (29) Các số hạng Fuxk(U) và rk sẽ được ước tính bằng các phép tích phân số và sai phân. Riêng số hạng Jk+1 sẽ được nội suy và sai phân thành: kkPkkTk AA qqJ 1111 ++++ +−= (30) Trong đó 1kT A +  và 1kP A +  được xác định tùy theo sơ đồ nội suy. Riêng tại đáy và mặt thoáng, điều kiện biên sẽ được sử dụng để tính J. Thay (30) vào (25), ta được phương trình cuối cùng: k n kkN n kkP n kkS AAA Srcqqq =−+− ++ ++ − 2/1 1 2/12/1 1 (31) Trong đó: ( ) kTkP k kP AA t V A ++   = +1 (32) 1+ = kTkT AA (33) kPkS AA = (34) ( ) nkknknkkk VFlux t V rUqSrc .2/12/1 +−   = −− (35) Trên mỗi thủy trực ta sẽ thiết lập được một hệ nz phương trình (31) chứa nz ẩn số mà sau khi giải ta sẽ có nz nghiệm qk (k=1,nz) tại cạnh ict trên các lớp. b. Bước 2 Để tính mực nước η, phương trình liên tục (1) được tích phân trên phương  trên toàn bộ chiều sâu. Kết quả ta thu được: ( ) 0 t D k kk = +   q (36) (36) sau đó được tích phân trên diện tích kiểm soát của nút C (Hình 1a), từ đó ta được:          −=    + + k j j 2/1n jknk 2/1n C Lq S 1 t D  (37) Từ (37), độ sâu và mực nước tại nút C được tính: t t D DD 2/1n Cn C 1n C    += + + (38) 1n Cb 1n C Dz ++ += (39) c. Bước 3 Để tính lưu lượng đơn vị trên phương , phương trình liên tục (1) được tích phân trên thể tích kiểm soát ở nút C tại lớp k: 0dSd q dSddSd t D SSS kkk =   ++               q (40) Các tích phân được ước tính bằng các tích phân số và sau đó sai phân theo thời gian ta được lời giải cho qσ ở nút C trên các lớp ở thời điểm n+1/2:  + + ++ + −   −= j j 2/1n jkn k 2/1n C k 2/1n k 2/1n 1k Lq St D qq   (41) Lưu ý là tại đáy 0 2/1 1 = +n q d. Bước 4 Phương trình (3) được tích phân trên thể tích kiểm soát tại nút C ở lớp k:   +−+      S H S C kk dSdCHDCdSd t q     q    =        −   + S C S V kk dSdDSdSd C D D Cq     (42) Sau khi sai phân số hạng đạo hàm theo thời gian (42) đưa tới phương trình: ( ) ( ) ( ) ( )1111 1 1 . ++++ + + −=−++−   n k n k n k n k n k n kC n kC k sCrDVJJCFluxqq t V (43) Trong đó: S.V kk  = (44) 1n 2/1k V1n 1k C D D CqSJ + + + +         −=   (45) ( )   −= + S k n H nn k n k dSCHDCCFlux  2/1 q (46) Số hạng Fuxk(Cn) sẽ được ước tính bằng các phép tích phân số và sai phân còn số hạng Jk+1 sẽ được nội suy và sai phân thành: kkPkkTk CACAJ 1111 ++++ +−= (47) Trong đó 1kT A +  và 1kP A +  được xác định tùy theo sơ đồ nội suy. Tại đáy và mặt thoáng, điều kiện biên sẽ được sử dụng để tính J. Thay (47) vào (43), ta được phương trình cuối cùng: k n kkN n kkP n kkS SrcCACACA =−+− ++ ++ − 1 1 11 1 (48) Với ( ) kTkP n kkP AADVs t A ++      +  = + + 1 11 (49) 1+ = kTkT AA (50) ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, VOL. 17, NO. 1.1, 2019 49 kPkS AA = (51) ( )nknnk n kk CFluxDrC t D VSrc −         +  = +1. (52) Trên mỗi thủy trực ta sẽ thiết lập được một hệ nz phương trình (48) chứa nz ẩn số mà sau khi giải ta sẽ có nz nghiệm qC (k=1,nz) ở nút C trên các lớp ở thời điểm n+1. e. Bước 5 Để tính biến đổi đáy, phương trình (4) được tích phân trên diện tích kiểm soát của nút C (Hình 1a): ( ) ( ) ( ) −=+  +   − S bb S b S b S b C dSEDdSdS t bc dS t z P q1 (53) Sau khi sai phân số hạng đạo hàm theo thời gian, (53) đưa tới phương trình: ( ) ( ) ( )SEDqFluxDC t zz P bbb n b n b C −=++  − − +1 1  (54) Trong đó ( ) ( ) nbnb bcbc t S DC −  = +1 (55) ( ) = L nbb dlqqFlux (56) Từ (54) ta có cao độ đáy ở nút C ở thời điểm n+1: ( ) ( ) ( ) bbb C n b n b qFluxDCSED P t zz −−− −  +=+ 1 1 (57) 3. Tính toán thử nghiệm 3.1. Vận tải bùn cát trong kênh cong 180o Odgaard và Bergs [5] có một nghiên cứu bằng thực nghiệm về dòng chảy và vận tải bùn cát trong một kênh cong. Bài toán này cũng đã Wu và ctv [10] tính toán bằng mô hình 3D. Kênh gồm một đoạn cong 180o, bán kính trung bình 13,11m, và 2 đoạn kênh dẫn vào và ra thẳng, dài 20m. Mặt cắt ngang kênh hình thang, đáy rộng 1.44m và mặt thoáng rộng 2.44m (Hình 2). Đáy kênh được phủ một lớp cát dày 23cm với đường kính trung bình 0.3 mm. Nước và cát trôi ra theo dòng chảy được bơm ngược trở lại đầu kênh trong một hệ thống tuần hoàn khép kín với lưu lượng 0,153m3/s. Sau khi đáy kênh biến đổi tới trạng thái ổn định, địa hình đáy kênh đã được đo đạc. Hình 2. Lưới tính mô hình 3D Trong nghiên cứu này mô hình số 3D cho kênh của Odgaard và Bergs đã được xây dựng. Trên mặt bằng, đoạn cong 180o được chia thành 180×16 phần tử còn mỗi đoạn kênh dẫn vào và ra được chia thành 80×16 phần tử (Hình 2). Mô hình có 11 lớp theo phương thẳng đứng. Địa hình đáy được lấy theo hình dạng của kênh cứng hình thang và điều kiện ban đầu là kênh đã bị bồi 23cm tương ứng với lớp cát phủ ban đầu trong thí nghiệm. Đầu vào của kênh được áp đặt điều kiện biên lưu lượng Q=0,153m3/s với nồng độ bùn cát lơ lửng tương ứng với lưu lượng bùn cát ra ở mặt cắt cuối kênh 3,7g/m/min như đã được ghi nhận trong thí nghiệm của Odgaard và Bergs. Điều kiện biên mực nước ở cuối kênh được áp đặt sao cho độ sâu trong kênh khoảng 15cm như trong thí nghiệm. a) Vận tốc trên mặt phẳng mặt cắt b) Vận tốc vuông góc với mặt cắt Hình 3. Vận tốc tại các mặt cắt ở các góc cong khác nhau Bờ trái Bờ phải Măt cắt 00 Bờ trái Bờ phải Măt cắt 00 Bờ trái Bờ phải Măt cắt 450 Bờ trái Bờ phải Măt cắt 450 Bờ trái Bờ phải Măt cắt 900 Bờ trái Bờ phải Măt cắt 900 Bờ phải Măt cắt 1350 Bờ trái Bờ phải Măt cắt 1350 Bờ trái Măt cắt 1800 Bờ trái Bờ phải Măt cắt 1800 Bờ trái Bờ phải 50 Trần Thị Mỹ Hồng, Lê Song Giang Hình 4. Địa hình đáy kênh tính toán (sau 9 giờ - đường liền và sau 11 giờ - đường đứt nét) và thí nghiệm (chấm tròn) Tính toán cho thấy sau khoảng 9 giờ kể từ thời điểm ban đầu địa hình đáy kênh gần như không còn thay đổi. Hình 3 giới thiệu trường vận tốc tại một số mặt cắt kênh ở các vị trí khác nhau. Các hình cho thấy xoáy thứ cấp trong kênh đã được tái hiện một cách rõ ràng. Hình 4 giới thiệu kết quả tính diễn biến đáy kênh sau khi đã ổn định. Đường liền là đáy kênh sau 9 giờ còn đường đứt đoạn là đáy kênh sau 11 giờ. Kết quả tính toán hoàn toàn phù hợp với số liệu thí nghiệm của Odgaard và Bergs. 3.2. Vận tải bùn cát khu vực Cù lao Phố trên sông Đồng Nai Sông Đồng Nai tại khu vực Cù lao Phố từng có một dự án san lấp lấn sông gây tranh cãi. Do địa hình lòng sông uốn lượn phức tạp, việc sử dụng mô hình 2D để nghiên cứu diễn biến bồi xói, đặc biệt khu vực đầu Cù lao Phố [11, 15] là không phù hợp vì mô hình không có khả năng tính toán mô phỏng các dòng thứ cấp. Để đánh giá khả năng của mô hình 3D phát triển trong bài báo, mô hình đã được áp dụng tính toán cho đoạn sông này. Trên mặt bằng, đoạn sông được chia thành 7261 phần tử tứ giác và 5 lớp theo chiều sâu. Số phần tử theo chiều ngang nhánh sông chính là 16 còn nhánh phụ là 8. Các phần tử có kích thước khoảng 16 – 25 m. Lưới tính này là khá mịn, đủ để mô phỏng chi tiết cấu trúc dòng chảy khu vực nghiên cứu. Địa hình đáy mô hình được tham khảo từ số liệu khảo sát bởi Viện Khoa học Thủy lợi miền Nam, thực hiện năm 2008 [14]. Nghiên cứu [14] cũng cho thấy vật liệu đáy sông khu vực không đồng nhất. Giữa dòng là cát mịn tới thô trong khi ở gần bờ là bùn. Tuy nhiên do hạn chế của phương pháp, mô hình chỉ có khả năng tính một loại vật liệu đáy. Vì vậy trong nghiên cứu này vật liệu đáy có đường kính trung bình d50=0,35mm được xem xét. Để giải quyết vấn đề điều kiện biên thủy lực của mô hình 3D, mô hình được tích hợp vào mô hình hệ thống sông Đồng Nai, kế thừa từ kết quả nghiên cứu [12]. Do mô hình 3D và mô hình hệ thống sông Đồng Nai sau tích hợp đã trở thành một mô hình thống nhất nên mực nước và vận tốc trên biên của mô hình 3D sẽ được tính ngay trong mô hình tích hợp. Điều kiện biên bùn cát cho mô hình 3D được áp đặt một cách độc lập ngay tại các mặt cắt biên. Khảo sát của Viện Khoa học Thủy lợi Miền Nam [14] cho thấy nồng độ bùn cát lơ lửng tại khu vực này vào khoảng 30 – 40 mg/l. Giá trị này đã được dùng làm điều kiện biên. Để đảm bảo điều kiện ổn định bước thời gian tính được lấy t = 0,6s. Mặc dù bước thời gian tính khá nhỏ nhưng tốc độ tính toán của mô hình khá tốt. Trên máy PC Core i5 3,2GHz, một giờ máy tính mô phỏng số được khoảng 8 giờ. Do mô hình hệ thống sông Đồng Nai đã được hiệu chỉnh khá tốt [12] và mô hình 3D có quy mô khá nhỏ nên ở quy mô toàn mô hình mô hình tích hợp từ hai mô hình này không cần phải hiệu chỉnh lại. Ngược lại, việc tinh chỉnh mô hình 3D với vận tốc điểm và hiệu chỉnh mô hình bùn cát 3D là cần thiết. Tuy nhiên do điều kiện không cho phép nên việc này chưa thực hiện được. Điều này có thể chấp nhận vì tính toán chỉ giới hạn ở mục tiêu chứng minh khả năng áp dụng mô hình 3D ở quy mô bài toán thực. a) Thành phần vận tốc tiếp tuyến với mặt cắt b) Thành phần vận tốc vuông góc với mặt cắt Hình 5. Trường vận tốc trên các mặt cắt ngang lúc triều xuống vào thời gian đỉnh lũ 2008 Tính toán mô phỏng đã được thực hiện trong thời gian từ ngày 30/8/2008-30/9/2008. Khoảng thời gian này được chọn vì trong thời gian đó đã có một trận lũ lớn xảy ra. Phân bố vận tốc tại một số mặt cắt ngang lúc triều xuống -0,50 -0,40 -0,30 -0,20 -0,10 0,00 -0,500,000,50 h , m x/b -0,50 -0,40 -0,30 -0,20 -0,10 0,00 -0,500,000,50 h , m x/b -0,50 -0,40 -0,30 -0,20 -0,10 0,00 -0,500,000,50 h , m x/b -0,50 -0,40 -0,30 -0,20 -0,10 0,00 -0,500,000,50 h , m x/b -0,50 -0,40 -0,30 -0,20 -0,10 0,00 -0,500,000,50 h , m x/b -0,50 -0,40 -0,30 -0,20 -0,10 0,00 -0,50-0,30-0,100,100,300,50 h , m x/b Mặt cắt 178o Mặt cắt 137o Mặt cắt 84o Mặt cắt 2o Mặt cắt 31o Mặt cắt 61o Măt căt 1-1 Bờ trái Bờ phải Bờ trái Bờ phải Măt căt 1-1 Măt căt 2-2 Bờ trái Bờ phải Măt căt 2-2 Bờ trái Bờ phải Măt căt 3-3 Bờ trái Bờ phải Măt căt 3-3 Bờ phải Bờ trái Măt căt 4-4 Bờ phải Bờ trái Măt căt 4-4 Bờ trái Bờ phải ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, VOL. 17, NO. 1.1, 2019 51 vào thời gian lũ cũng được giới thiệu trên Hình 5. Vị trí các mặt cắt được giới thiệu trên Hình 6. Hình 6. Vị trí các mặt cắt Độ bồi xói khu vực đầu Cù lao Phố được giới thiệu trên Hình 7. Để thấy rõ hơn bức tranh bồi xói tại khu vực, ta có thể xem một vài lát cắt. Hình 8 giới thiệu phân bố vận tốc và biến đổi cao độ đáy tại các lát cắt 5-5 và 6-6. Vị trí các lát cắt này được chỉ ra trên Hình 7. Hình 7. Trường vận tốc và độ bồi xói tại khu vực đầu cù lao Phố (đơn vị thang màu: mét) Hình 8. Phân bố vận tốc và biến đối đáy tại các lát cắt (đường màu đen: đáy sông ban đầu; đường màu đỏ: đáy sông sau một tháng) Lát cắt 5 – 5 cho thấy rõ các mô ngầm bị dòng chảy bào mòn và di chuyển dần về hạ lưu. Trong khi đó lát cắt 6 - 6 lại cho thấy cấu trúc 3D của dòng chảy ở đầu cù lao. Có một xoáy ngang xuất hiện ở đầu cù lao và xoáy này mang bùn cát từ chân cù lao ra ngoài. Xoáy này là một trong các nguyên nhân chính làm cho đầu cù lao bị xói mạnh. Khả năng mô phỏng các dòng thứ cấp này cũng chính là ưu thế của mô hình 3D so với mô hình 2D. 4. Kết luận Bài báo đã trình bày một phương pháp tính toán dòng chảy và vận tải bùn cát, biến hình lòng dẫn 3D. Các phương trình vi phân cơ bản được giải bằng phương pháp thể tích hữu hạn trên lưới tính phi cấu trúc phần tử tứ giác. Kết quả tính toán kiểm tra với số liệu thực nghiệm cho thấy phương pháp cho lời giải hợp lý và khá chính xác. Bên cạnh đó việc tính toán trường hợp diễn biến lòng dẫn khu vực Cù lao Phố cũng cho thấy phương pháp cho kết quả phù hợp quy luật và hoàn toàn có thể sử dụng trong nghiên cứu thực tế. Lời cảm ơn Nghiên cứu này được tài trợ bởi Trường Đại học Bách Khoa TpHCM trong đề tài ‘‘Nghiên cứu tính toán dòng chảy và chuyển tải bùn cát bằng mô hình 3D và ứng dụng tính toán cho khu vực cù lao Phố, sông Đồng Nai” mã số T- KTXD-2017-75 TÀI LIỆU THAM KHẢO [1] Davies, A.M. and Gerritsen H., (1994). An Intercomparison of Three-dimensional Tidal Hydrodynamic Models of the Irish Sea, Tellus, 46A, pp. 200-221. [2] Hayter, E.J. and Mehta, A.J., 1986. Modeling cohesive sediments transport in estuarine waters. Applied Mathematical Modelling 51, pp. 765-778. [3] Krone, R.B., 1962. Flume studies of the transport of sediment in estuarine shoaling processes. Final report to San Fransico District U.S. Army Corps of Engineers, Washington D.C [4] Lin, B.; R.A. Falconer (1997). Three-dimensional layer integrated modeling of estuarine flows with flooding and drying. Estuar. Coast. Shelf Sci., 44, 737–751. [5] Odgaard, A.J. and M.A. Bergs (1988). Flow Processes in a Curved Alluvial Channel, Water Resour. Res., 24(1): 45-56. [6] Smagorinsky, J. (1963). General Circulation Experiments with the Primitive Equations. I: The Basic Experiment. Monthly Weather Rev., 91, 99 - 164. [7] Van Rijn., J. Hydraulic Eng., ASCE, Vol. 110 (10), pp. 1431-1456. (1984a) [8] Van Rijn., J. Hydraulic Eng., ASCE, Vol. 110 (11), pp. 1613-1641. (1984b) [9] Van Rijn, L.C. (1993). Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas, Aqua Publications, 690p. [10] Wu, W., W. Rodi and T. Wenka (2000). 3D Numerical modeling of Flow and Sediment Transport in Open Channel, J. Hydraul. Eng., 126(1): 4-15. [11] Cty cổ phần ĐT-KT-XD Toàn Thịnh Phát (2014). Báo cáo đánh giá tác động môi trường Dự án cải tạo cảnh quan và phát triển đô thị ven sông Đồng Nai, quy mô 8,4ha tại phường Quyết Thắng, Tp. Biên Hòa, tỉnh Đồng Nai. [12] Lê Song Giang (2017). Nghiên cứu đề xuất lựa chọn chiến lược quản lý ngập lụt thích hợp trên cơ sở các dự án đã, đang và dự kiến triển khai tại Tp.HCM. Báo cáo khoa học tổng kết đề tài. [13] Trần Thị Mỹ Hồng, Lê Song Giang (2017). Một sơ đồ tính toán dòng chảy sông, biển 3 chiều trên lưới tính phi cấu trúc. Tuyển tập Công trình Hội nghị khoa học Cơ học Thủy khí toàn quốc lần thứ 20, trg. [14] Viện Khoa học Thủy lợi miền Nam (2009). Đánh giá tác động dòng chảy sông Đồng Nai đoạn từ cầu Hóa An đến cầu Ghềnh thuộc thành phố Biên Hòa. Tp. Hồ Chí Minh, tháng 10 năm 2009. [15] Kolmogorov, A. N. (1942), "Equations of Turbulent Motion of an Incompressible Fluid," Izvestia Academy of Sciences, USSR; Physics, Vol. 6, Nos. 1 and 2, pp. 56-58. [16] Prandtl, L. (1945), "Uber ein neues Formelsystem fur die ausgebildete Turbulenz," Nacr. Akad. Wiss. Gottingen, Math-Phys. Kl. 1945, pp6-19. (BBT nhận bài: 31/10/2018, hoàn tất thủ tục phản biện: 15/01/2019) 1 1 2 2 3 3 4 4 6 6 5 5 Thượng lưu Hạ lưu Lát cắt 5 - 5 Thượng lưu Hạ lưu Lát cắt 6 - 6