Nghiên cứu mô hình toán mô phỏng dòng chảy hở một chiều có kể đến vận tốc theo chiều đứng tại đáy

KHOA HỌC KỸ THUẬT THỦY LỢI VÀ MÔI TRƯỜNG - SỐ 61 (6/2018) 91 BÀI BÁO KHOA HỌC NGHIÊN CỨU MÔ HÌNH TOÁN MÔ PHỎNG DÒNG CHẢY HỞ MỘT CHIỀU CÓ KỂ ĐẾN VẬN TỐC THEO CHIỀU ĐỨNG TẠI ĐÁY Huỳnh Phúc Hậu1, Nguyễn Thế Hùng2, Trần Thục3, Lê Thị Thu Hiền4 Tóm tắt: Trong bài báo này, phương pháp phần tử hữu hạn Taylor-Galerkin được áp dụng để rời rạc hóa hệ phương trình Saint-Venant có kể vận tốc chiều đứng ở đáy lòng dẫn, với độ chính xác bậc ba theo thời gian và không gian. Mô hình toán được kiểm định bởi hai ví dụ: Dòng chảy ổn định trên kênh có vật cản và dòng chảy vỡ đập trên kênh dốc. Kết quả cho thấy tính hiệu quả và chính xác của mô hình toán. Mô hình vật lý được xây dựng nhằm tạo ra vận tốc chiều đứng ở đáy kênh để kiểm chứng tính đúng đắn của mô hình. Kết quả đo đạc về sự biến đổi mực nước dọc máng thí nghiệm được thực hiên với các cấp lưu lượng khác nhau. Kết quả này được so sánh với kết quả tính toán theo mô hình toán cho thấy sự phù hợp tốt khi chỉ số Nash trong các trường hợp lên tới gần 90%. Từ khóa: Saint-Venant, Taylor-Galerkin, thí nghiệm, xáo trộn đáy lòng dẫn. 1. ĐẬT VẤN ĐỀ 1 Hệ phương trình vi phân phi tuyến Saint- Venant (hay còn được xem là hệ phương trình nước nông một chiều) đã và đang được sử dụng rộng rãi trong việc mô phỏng dòng chảy không ổn định một chiều trên lòng dẫn hở. Trong những năm gần đây, đã có nhiều nghiên cứu về việc giải hệ phương trình này khi xét tới dòng chảy chịu ảnh hưởng của trọng lực hay lực Coriolit (Lai và nnk, 2012; Pilotti và nnk, 2011). Tuy nhiên, ảnh hưởng của sự xáo trộn ở đáy lòng dẫn do có dòng chảy bổ sung ở đáy thì chưa được xem xét. Vì vậy, các tác giả đã xét tới thành phần này và bổ sung vào số hạng nguồn của hệ phương trình Saint -Venant. Mặt khác, việc lựa chọn phương pháp số phù hợp để giải hệ phương trình này cũng là vấn đề được nhiều nhà khoa học quan tâm nghiên cứu. Lai và nnk(2012) dùng phương pháp phần tử hữu hạn discontinuous Galerkin để giải; Pilotti và nnk(2011) lại dùng phương pháp sai phân hữu hạn Mac-Cormack để có được nghiệm chính xác bậc hai theo không gian và thời gian. Tuy nhiên, số hạng nguồn mới chỉ xét tới ảnh hưởng của độ 1Trường Cao đẳng Giao thông Vận tải Trung ương V 2Trường Đại học Bách Khoa, Đại học Đà Nẵng 3Viện khoa học khí tượng thủy văn và biến đổi khí hậu 4Bộ môn Thủy lực, Trường Đại học Thủy lợi. dốc đáy và ma sát. Vì vậy, trong nội dung bài báo này, các tác giả đã dùng phương pháp phần tử hữu hạn Taylor-Galerkin để rời rạc hóa hệ phương trình Saint-Venant có kể sự xáo trộn ở đáy lòng dẫn, với độ chính xác bậc ba theo thời gian và không gian. Sau đó dùng ngôn ngữ lập trình Fortran để xây dựng chương trình tính. Tính chính xác, tính ổn định và hiệu quả của sơ đồ số được kiểm định bằng một số ví dụ có nghiệm giải tích hay thực đo cũng được chỉ ra trong bài báo này. Bên cạnh đó, để đánh giá khả năng của mô hình toán trong việc mô phỏng ảnh hưởng của dòng chảy bổ sung theo chiều đứng, mô hình vật lý được thiết lập và đo đạc tại Phòng Thí nghiệm trọng điểm Quốc gia. Kết quả về đường mặt nước giữa tính toán và thực đo khá phù hợp khi chỉ số Nash trong các trường hợp thí nghiệm lên tới 90%. Hệ phương trình vi phân đạo hàm riêng của dòng chảy một chiều khi có kể đến xáo trộn ở đáy lòng dẫn được giải số bằng phương pháp phần tử hữu hạn Taylor–Galerkin và lập trình bằng ngôn ngữ Fortran (Huỳnh Phúc Hậu, Nguyễn Thế Hùng, 2017). Để kiểm chứng tính chính xác của mô hình toán, thí nghiệm trên mô hình vật lý đã được thực hiện và trình bày trong bài báo này. KHOA HỌC KỸ THUẬT THỦY LỢI VÀ MÔI TRƯỜNG - SỐ 61 (6/2018) 92 2. MÔ HÌNH TOÁN 2.1. Hệ phương trình Sant Venant có kể đến xáo trộn ở đáy lòng dẫn qw h A x Q t A           qvR A QQ gngAi x hAag x AQ t Q           3/42 2 / (1) Trong đó: h: độ sâu dòng chảy (m); Q: lưu lượng dòng chảy (m3/s); q: lưu lượng bổ sung dọc sông (m2/s); g: gia tốc trọng trường (m/s2); i: độ dốc đáy lòng dẫn; n: hệ số nhám lòng dẫn. A: diện tích mặt cắt ướt (m2): mhhbA 20 5.0 ; b0: bề rộng đáy; m: tổng 2 hệ số mái dốc; R: bán kính thủy lực (m). Dòng chảy bổ sung tại đáy lòng dẫn gây xáo trộn, có vận tốc w và gia tốc t wa    Viết lại hệ phương trình Saint Venant theo cặp biến (h, Q), ta được:                qw h A hAx Q hAt h / 1 / 1     qvR A QQ gngAi x hAag x AQ t Q           3/42 2 / (2) Viết thành dạng vector: )()( pS x pf t p       (3) Hay: )()( pS x ppD t p       (4) Trong đó vec-tơ ẩn p=(h,Q)T ; f là thông lượng: Ma trận Jacobian D(p) được tính bằng biểu thức (5)                    A QAag h A A Q hApD p pf 2 / 10 )()( 2 2 (5) Số hạng nguồn trong phương trình (3) được xác định bằng: T A QqR A QQ gngAiqw h A hA pS                 3/42, / 1)( (6) 2.2. Rời rạc theo thời gian Thực hiện khai triển véc tơ ẩn pn+1 bằng chuỗi Taylor theo t lân cận bên phải điểm thời gian t=tn; đến bậc ba, nhận được:       3 32 1 62 tOptptptpp tttn tt n t nnn            ttnttntnnn ptptptpp 2 1 2 1 23 1 6 1 2                (7) Trong đó:  là trọng số ẩn, tnp là đạo hàm bậc nhất theo thời gian của p đánh giá tại t= t n. Và tương tự như vậy, ttnp  là đạo hàm bậc hai:              )()()()( pS x pfpS x pf t p (8) Vậy:   t ppB t ppD xt p p pS t p p pf x pS tx pf tt p                                  )()()()()()(2 2 KHOA HỌC KỸ THUẬT THỦY LỢI VÀ MÔI TRƯỜNG - SỐ 61 (6/2018) 93                                                   )()()()()()( )()()()( 2 2 2 2 pS x pfpBpS x pfpD xt p t p p pS t p p pf x pS tx pf tt p (9) Thay thế (8) và (9) vào phương trình (7):     nn n n n pS x ppDpBpS x ppDpD x tpS x ppD tppS x ppDpBpS x ppDpD x tp                                                                                )()()()()()( 23 1)()( )()()()()()( 6 1 2 2 1 2 1   (10) 2.3. Rời rạc theo không gian Gọi chiều dài phần tử một chiều bậc 2 là 2L, có 3 nút 1,2,3. Chọn gốc tọa độ địa phương tại nút đầu 1, hướng x dương từ nút đầu 1 đến nút cuối 3. Chọn hàm nội suy bậc 2, ta có:          21 2 2 200 2 L LxLx LL LxLx                222 22 20 20 L xLx L Lxx LLL Lxx                 22 2022 0 L Lxx LLL xLx      Áp dụng tích phân trọng số cho phương trình (10) ở trên, áp dụng tích phân từng phần cho đạo hàm bậc 2 ta được hệ phương trình đại số tuyến tính để xác định phương trình hệ ma trận phần tử, sau khi ghép nối được hệ phương trình tổng thể, gán điều kiện biên để giải ra vec tơ ẩn số ở từng bước thời gian. Các tác giả đã sử dụng ngôn ngữ lập trình Fortran90 xây dựng chương trình tính dựa trên mô hình toán đã chọn. Phương pháp số đã được kiểm định tính bảo toàn khối lượng, không xuất hiện nhiễu động, tính chính xác của kết quả phương pháp số v.v Một số các ví dụ nhằm kiểm định tính đúng đắn của mô hình được chỉ ra trong mục 2.5. 2.4. Kiểm định mô hình toán a. Dòng chảy ổn định trên kênh có vật cản Ví dụ này nhằm mô phỏng dòng chảy ổn định trên kênh có vật cản (Hou và nnk, 2013). Kênh dẫn mặt cắt chữ nhật dài 25m, độ nhám coi như bằng 0. Cao độ đáy kênh được định dạng bằng biểu thức:            128xk0 128k1005.02.0 2 xhi mxmhixzb  (11) Trường hợp 1: Dòng chảy trên kênh chuyển là chuyển tiếp, không có sóng gián đoạn. Độ sâu hạ lưu là 0,66m, lưu lượng đơn vị phía thượng lưu là q = 1,53m3/s.m. Trường hợp 2: Dòng chảy trên kênh chuyển là chuyển tiếp, có sóng gián đoạn. Độ sâu mực nước hạ lưu 0,33m, lưu lượng đơn vị phía thượng lưu là q = 0,18m3/s.m. Kết quả quá trình mực nước và lưu lượng đơn vị tính theo phương pháp số được so sánh với kết quả giải tích cho thấy có sự phù hợp cao. Vì vậy, mô hình toán do các tác giả lựa chọn có khả năng mô phỏng dòng chảy ổn định trên kênh có địa hình phức tạp. KHOA HỌC KỸ THUẬT THỦY LỢI VÀ MÔI TRƯỜNG - SỐ 61 (6/2018) 94 Hình 1. Quá trình mực nước và lưu lượng đơn vị trong hai trường hợp 1 và 2 b. Dòng chảy vỡ đập trên kênh dốc Thí nghiệm này được thực hiện tại phòng thí nghiệm US Army Engineer Waterway Experiment Station (Bellos và nnk, 1987) nhằm kiểm tra khả năng của mô hình trong việc mô phỏng dòng chảy do vỡ đập trên kênh dốc. Kênh lăng trụ mặt cắt chữ nhật dài 122m, rộng 1,22m có độ dốc đáy So =0,005, hệ số nhám Manning lấy bằng 0,009. Đô sâu mực nước trước đập là h1 = 0,305m, kênh hạ lưu là khô. Đường quá trình độ sâu nước tại các vị trí x=70,1m và 85,1m được chỉ ra trên hình 2. Kết quả giữa mô hình toán và thực đo chỉ ra rằng mô hình toán đã chọn cho kết quả hoàn toàn phù hợp với thực đo với chỉ số Nash tương ứng là 87,25% và 89,1%. Hình 2. Quá trình mực nước tại các vị trí x=70,1m và x=85,4m. Những ví dụ trên cho thấy, phương pháp số các tác giả lựa chọn hoàn toàn phù hợp. Để đánh giá sự ảnh hưởng của nhiễu động sinh ra do có dòng bổ sung theo chiều đứng tại đáy kênh. Các tác giả đã xây dựng mô hình vật lý. Kết quả đo đạc mực nước được so sánh với kết quả tính toán theo mô hình toán được trình bày trong mục 3. 3. MÔ HÌNH VẬT LÝ 3.1. Mô tả thí nghiệm Thí nghiệm kiểm chứng mô hình toán về dòng chảy hở một chiều có sự xáo trộn ở đáy lòng dẫn được thực hiện tại Phòng thí nghiệm trọng điểm quốc gia về động lực học sông, biển. Mô hình thí nghiệm: Máng kính mặt cắt ngang chữ nhật rộng 50 cm, cao 1m, dài 15m. Để tạo điều kiện biên là vận tốc chiều đứng tại đáy dòng chảy, Máng kính được chia thành 2 phần: phần dòng chảy trên và dưới được ngăn cách bởi lớp bê tông dày 5cm và lớp vữa xi măng dày 25cm xoa phẳng. Phần dưới gọi là đường hầm có bề rộng 0,44m, chiều cao 0,15m. Thiết bị đo lưu lượng sử dụng trong thí nghiệm là đập lường thành mỏng tiết diện chữ nhật có bề rộng b=0,6m; chiều cao đập lường P=0,75m. Công thức đo lưu lượng: gHmbHQ 2 với hệ số lưu lượng m = 0,402+0,054.H/P, trong đó H: chiều sâu nước trên đỉnh đập lường (m). KHOA HỌC KỸ THUẬT THỦY LỢI VÀ MÔI TRƯỜNG - SỐ 61 (6/2018) 95 Hình 3. Máng thí nghiệm M¸ng l­êng h×nh thang ®o l­u l­îng TÊm lÆng sãng §æ c¸t x©y tr¸t mÆt §æ c¸t x©y tr¸t mÆt Cöa ra khe ®¸yi=1%bê tông ®­êng hÇm m¸ng kÝnh cã s½n Hình 4. Thông số kỹ thuật máng kính thí nghiệm 3.2. Tiến hành thí nghiệm Mặt cắt số 1 (MC1) cách tâm khe đáy 350cm về thượng lưu. MC2 cách tâm khe 300cm về thượng lưu. MC3 cách tâm khe 200cm về thượng lưu. MC 4 cách tâm khe 100cm về thượng lưu. MC5 tại tâm khe đáy. MC6 cách tâm khe 100cm về hạ lưu. MC 7 cách tâm khe 200cm về hạ lưu. MC8 cách tâm khe 300cm về hạ lưu. MC 9 cách tâm khe 400cm về hạ lưu. MC10 cách tâm khe 450cm về hạ lưu. Giữa MC 4 và MC6 chia nhỏ thành các mặt cắt cách nhau 10cm do giữa hai mặt cắt này mực nước biến đổi nhiều. Các cấp lưu lượng tổng Q: 70; 75; 80; 90; 95; 100; 105 (l/s). Các cấp lưu lượng dòng chính phía trên Q1: 45; 50; 60; 65; 70; 75(l/s). Lưu lượng bổ sung Q2 = Q-Q1 Chiều sâu được đo bằng thước thép, máy thủy bình và mia. Mỗi mặt cắt ngang đo 3 thủy trực để lấy trị số trung bình. 4. KẾT QUẢ THÍ NGHIỆM VÀ THẢO LUẬN 4.1. Kết quả đo độ sâu mực nước Khe đáy tạo vận tốc chiều đứng Tấm giảm sóng Khe đáy máng hình thang đo lưu lượng Hạ lưu KHOA HỌC KỸ THUẬT THỦY LỢI VÀ MÔI TRƯỜNG - SỐ 61 (6/2018) 96 Bảng 1. Độ sâu mực nước khi Q=75÷105(l/s); Q2 = 30 (l/s) Độ sâu mực nước (cm) tại cấp lưu lượng tổng Q (l/s) STT Tên mặt cắt 75 80 90 95 100 105 Ghi chú 1 MC1 22,64 22,84 23,74 23,97 24,31 24,87 Cách tâm khe 350cm về thượng lưu 2 MC2 23,49 23,67 24,09 25,17 25,54 26,04 Cách tâm khe 300cm về thượng lưu 3 MC3 23,54 23,82 24,39 25,24 25,89 26,49 Cách tâm khe 200cm về thượng lưu 4 MC4 22,64 23,49 24,49 24,84 25,49 26,04 Cách tâm khe 100cm về thượng lưu 5 MC5 20,99 21,99 23,09 23,59 24,14 25,29 Tâm khe 6 MC6 10,34 11,24 11,79 12,64 12,89 13,49 Cách tâm khe 100cm về hạ lưu 7 MC7 9,99 10,79 11,54 12,57 12,76 13,67 Cách tâm khe 200cm về hạ lưu 8 MC8 9,64 10,44 11,24 12,46 12,62 13,81 Cách tâm khe 300cm về hạ lưu 9 MC9 9,77 9,89 11,16 11,97 12,34 13,27 Cách tâm khe 400cm về hạ lưu 10 MC10 9,74 9,84 11,11 11,87 12,29 13,24 Cách tâm khe 450cm về hạ lưu Bảng 2. Độ sâu mực nước chi tiết giữa mặt cắt 4 và 6 MC Q=75-45-30 Q=80-50-30 Q=90-60-30 Q=95-65-30 Q=100-70-30 Q=105-75-30 1-4 22,65 23,50 24,50 24,85 25,50 26,05 2 22,65 23,50 24,50 24,80 25,50 26,05 3 22,65 23,50 24,50 24,75 25,45 26,05 4 22,70 23,45 24,50 24,75 25,50 25,95 5 22,65 23,45 24,50 24,75 25,45 25,95 6 22,80 23,40 24,40 24,75 25,55 25,90 7 23,00 23,35 24,40 24,65 25,50 25,95 8 22,70 23,25 24,30 24,70 25,40 25,90 9 22,35 22,85 24,10 24,50 25,30 25,80 10 22,20 22,65 24,00 24,30 25,00 25,75 11-5 21,00 22,00 23,10 23,60 24,15 25,25 12 19,50 19,50 20,55 21,50 22,50 23,50 13 17,00 17,65 18,75 19,20 20,05 20,50 14 14,15 14,95 16,25 16,90 17,90 18,50 15 12,00 13,00 14,50 15,30 15,80 16,50 16 11,25 12,10 13,50 14,35 14,70 15,50 17 10,85 11,60 12,90 13,50 14,20 14,70 18 10,70 11,45 12,50 13,10 13,80 14,25 19 10,55 11,30 12,10 12,80 13,50 13,95 20 10,50 11,20 11,90 12,65 13,25 13,65 21-6 10,35 11,25 11,80 12,65 12,90 13,50 4.2. So sánh kết quả thí nghiệm và kết quả giải số trên mô hình toán Thí nghiệm nhằm kiểm chứng thuật toán và chương trình tính đã thiết lập (Huỳnh Phúc Hậu, 2016, 2017); qua kết quả so sánh giữa thí nghiệm và tính toán ở các hình 5 đến 10 (sai số tương đối max là 5,5%) cho thấy tính đúng đắn của thuật toán và chương trình tính. KHOA HỌC KỸ THUẬT THỦY LỢI VÀ MÔI TRƯỜNG - SỐ 61 (6/2018) 97 Hình 5. Chiều sâu nước với lưu lượng dòng chính phía trên Q=45(l/s), Hình 6. Chiều sâu nước với lưu lượng dòng chính phía trên Q=50(l/s), Hình 7. Chiều sâu nước với lưu lượng dòng chính phía trên Q=60 (l/s) Hình 8. Chiều sâu nước với lưu lượng dòng chính phía trên Q=65(l/s), Hình 9. Chiều sâu nước với lưu lượng dòng chính phía trên Q=70(l/s) Hình 10. hiều sâu nước với lưu lượng dòng chính phía trên Q=75(l/s), KHOA HỌC KỸ THUẬT THỦY LỢI VÀ MÔI TRƯỜNG - SỐ 61 (6/2018) 98 5. KẾT LUẬN Bài báo này đã giải quyết mô hình toán dòng chảy hở một chiều Saint-Venant dưới ảnh hưởng bởi vận tốc theo chiều đứng tại đáy. Với phương số là phần tử hữu hạn Taylor-Galerkin, kết quả thu được có độ chính xác bậc ba theo thời gian và không gian. Chương trình tính được kiểm nghiệm tính đúng đắn và hiệu quả thông qua 2 ví dụ có nghiệm giải tích và thực đo. Sau đó, một thí nghiệm được thực hiện ở phòng Thí nghiệm Trọng điểm Quốc gia và Động lực Sông biển để kiểm chứng thuật toán và chương trình tính xây dựng cho mô hình toán dòng chảy hở một chiều dưới ảnh hưởng bởi vận tốc theo chiều đứng tại đáy. Kết quả thí nghiệm được dùng để so sánh với kết quả giải số trên mô hình toán, cho thấy có sự phù hợp tốt, điều này cho thấy thuật toán và chương trình tính có độ tin cậy cao. TÀI LIỆU THAM KHẢO Huỳnh Phúc Hậu, Nguyễn Thế Hùng (2017), "Áp dụng phương pháp phần tử hữu hạn Taylor - Galerkin giải bài toán dòng chảy hở một chiều không ổn định có sự xáo trộn ở đáy lòng dẫn", Hội nghị cơ học toàn quốc 2017, Hà Nội. Huỳnh Phúc Hậu, Nguyễn Thế Hùng (2017), "Áp dụng phương pháp Runghe-Kutta bậc bốn giải bài toán dòng chảy hở một chiều ổn định có sự xáo trộn ở đáy lòng dẫn", Hội nghị cơ học toàn quốc 2017, Hà Nội. W. Lai, A.A, Khan (2014), "Discontinuous Galerkin Method for 1D shallow water flow in Natural Rivers”, J, Engineering Application of Computational Fluid Mechanics, 6, 74-86. M. Pilotti, A. Maranzoni, M. Tomirotti và G. Valerio (2011), “1923-Gleno Dam Break: Case Study and Numerical Modelling”, J, Hydraulic Engineering, 137( 4) (ASCE), 480-492. J. Hou, F. Simons, M. Mahgoub, R. Hinkelmann (2013), “A robust well balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography”, Comput. Methods Appl. Mech. Engrg, 257, 126-149. C. Bellos, G. J. Sakkas (1987), “1D Dam break flood wave propagation on dry bed”, Journal of Hydraulic Engineering, 113(12), 1510-1524 ASCE. Abstract: STUDY A NUMERICAL MODEL FOR SOLVING THE ONE DIMENSIONAL FLOW ACOUNTING FOR VERITCAL VELOCITY AT THE BOTTOM OF CHANNEL This paper is investigated the Taylor–Galerkin finite element method to solve Saint-Venant equations accouting for additional discharge at the bottom of channel. The numerical solution with third order accuracy in space and time is validated by some reference test cases. A physical model is implemented at Key Labolatory of River and Coastal Engineering to verify the capacity of the proposed numerical model in terms of capturing correctly water hydrographs with different cases of discharges. The very good agreement between numerical rerults and experimental ones of can be observered. Keywords: Taylor-Galerkin, Saint Venant, Expreriment, bed disturbance. Ngày nhận bài: 18/5/2018 Ngày chấp nhận đăng: 11/6/2018