Numerical model for simulation of waves in surfzone and nearshore areas based on Boussinesq equations: Results for plane beaches

13 Vietnam Journal of Marine Science and Technology; Vol. 20, No. 1; 2020: 13–24 DOI: https://doi.org/10.15625/1859-3097/20/1/15037 Numerical model for simulation of waves in surfzone and nearshore areas based on Boussinesq equations: results for plane beaches Phung Dang Hieu * , Le Duc Dung, Nguyen Thi Khang Vietnam Institute of Seas and Islands, Hanoi, Vietnam * E-mail: hieupd@visi.ac.vn/phunghieujp@gmail.com Received: 9 April 2019; Accepted: 12 September 2019 ©2020 Vietnam Academy of Science and Technology (VAST) Abstract A numerical model based on the 2D Boussinesq equations has been developed using the Finite Volume Method. The model was verified against experimental data for the case of wave breaking on a sloping beach. Simulated results by the model showed that the model has good capability of simulation of waves in the nearshore area. Numerical simulation was also carried out for the problem of waves on a plane beach with a breakwater and submerged dunes. Simulated results were compared with those computed by MIKE 21. The comparison showed that good agreements were obtained and confirmed the applicability of the Boussinesq model to the simulation of physical phenomena of waves in the nearshore areas, especially, suitable for the simulation of wave-induced current including rip currents. Keywords: Boussinesq model, wave induced current, FVM, nearshore dynamics. Citation: Phung Dang Hieu, Le Duc Dung, Nguyen Thi Khang, 2020. Numerical model for simulation of waves in surfzone and nearshore areas based on Boussinesq equations: results for plane beaches. Vietnam Journal of Marine Science and Technology, 20(1), 13–24. 14 Tạp chí Khoa học và Công nghệ Biển, Tập 20, Số 1; 2020: 13–24 DOI: https://doi.org/10.15625/1859-3097/20/1/15037 Mô hình số mô phỏng sóng ven bờ và trong vùng sóng đổ dựa trên hệ phương trình Boussinesq: một số kết quả thử nghiệm cho bãi biển thoải Phùng Đăng Hiếu*, Lê Đức Dũng, Nguyễn Thị Khang Viện Nghiên cứu Biển và Hải đảo, Hà Nội, Việt Nam * E-mail: hieupd@visi.ac.vn/phunghieujp@gmail.com Nhận bài: 9-4-2019; Chấp nhận đăng: 12-9-2019 Tóm tắt Mô hình số sử dụng phương trình Boussinesq hai chiều được phát triển dựa trên phương pháp thể tích hữu hạn (FVM). Mô hình được kiểm nghiệm bằng việc áp dụng tính toán mô phỏng cho trường hợp sóng lan truyền, biến dạng trên bãi thoải. Kết quả tính toán được so sánh với số liệu thí nghiệm vật lý đã xuất bản, nhằm minh chứng khả năng mô phỏng sóng ven bờ. Mô hình số cũng được áp dụng mô phỏng cho bài toán sóng trên bãi nghiêng có đê chắn sóng nổi và có cồn ngầm. Kết quả so sánh với mô phỏng bằng phần mềm MIKE 21 để có so sánh đánh giá. Kết quả cho thấy có sự phù hợp khá và mô tả được tốt qui luật vật lý của sóng trong khu vực ven bờ, đặc biệt phù hợp cho mô phỏng hệ thống dòng chảy do sóng bao gồm cả dòng rút. Từ khoá: Mô hình Boussinesq, dòng phát sinh do sóng, thể tích hữu hạn, động lực ven bờ. GIỚI THIỆU Xây dựng mô hình tính toán sóng ven bờ từ hệ phương trình Boussinesq cần thiết phải giải quyết được một số vấn đề rất quan trọng và khó đó là: Tính toán được tiêu tán năng lượng do sóng đổ, giải quyết được sóng leo trên bãi biển, sơ đồ số phải bảo toàn, có độ chính xác tốt. Bên cạnh đó phương pháp giải số áp dụng phải ổn định, khả thi đảm bảo tính vật lý của quá trình. Nếu phương pháp số sử dụng không ổn định, sóng tiếp cận bờ do tính chất phi tuyến mạnh, tương tác phức tạp sẽ dẫn đến nhiễu số và phá vỡ mạnh giải số của mô hình làm tràn số trong quá trình tính toán. Các phương pháp tính toán sóng ven bờ bằng sai phân hữu hạn hay phần tử hữu hạn thường mắc phải là không ổn định số đối với các khu vực có địa hình phức tạp, sóng đổ và sóng tràn bãi với tính phi tuyến lớn. Chính vì vậy, cho đến nay rất ít mô hình cho phép mô phỏng được đầy đủ các quá trình sóng ven bờ bao gồm sóng đổ, sóng tràn bãi biển và hệ thống dòng chảy sóng ven bờ mà có thể ứng dụng tốt trên thực tế. Trên thế giới, các nhà khoa học đã quan tâm nghiên cứu phát triển mô hình toán mô phỏng sóng ven bờ dựa trên hệ phương trình Boussinesq trong nhiều thập kỷ qua. Các nghiên cứu phát triển mô hình số dựa trên hệ phương trình Boussinesq tiêu biểu có thể kể ra như Schaffer et al., (1993) [1], Madsen et al., (1997) [2, 3], Kennedy et al., (2000) [4], Kirby et al., (1995) [5] và một số tác giả khác. Thành công từ các nghiên cứu phát triển các mô hình số đó đã đưa ra các mô hình mã nguồn mở cho cộng đồng khoa học biển trên khắp thế giới sử dụng thí dụ như bộ chương trình FUNWAVE do Kirby và cộng sự phát triển, PCOULWAVE của Hoa Kỳ, hay mô hình của Madsen và cộng sự đã được phát triển tiếp để trở thành mô đun BW trong bộ phần mềm thương mại MIKE 21. Các nghiên cứu sử dụng hệ phương trình Boussinesq mở rộng tiếp tục được quan tâm và Numerical model for simulation of waves 15 cải tiến bởi cộng đồng các nhà khoa học về thủy động lực biển ven bờ trên khắp thế giới. Các nghiên cứu chủ yếu tập trung vào cải tiến sơ đồ số để tăng tính ổn định và giải quyết các vấn đề khác như tiêu tán năng lượng sóng đổ tốt hơn,... Ở nước ta, việc nghiên cứu tự xây dựng mô hình số thành chương trình máy tính để ứng dụng cho các nghiên cứu cũng như ứng dụng thực tiễn còn rất ít. Đặc biệt với bài toán sóng biển ven bờ sử dụng hệ phương trình Boussinesq mở rộng thì còn hiếm hơn. Các nhà khoa học động lực biển, ven bờ ở nước ta chủ yếu sử dụng các chương trình máy tính mã nguồn mở hoặc phần mềm MIKE 21 của nước ngoài để mô phỏng, tính toán sóng ven bờ cho các mục tiêu khác nhau. Mặc dù vậy, cũng có một vài tác giả đã bước đầu nghiên cứu phát triển mô hình số dựa trên hệ phương trình Boussinesq cho bài toán sóng dài (sóng thần) hay sóng tàu như Phùng Đăng Hiếu (2008) [6], Nguyễn Bá Thủy và nnk., (2016) [7] hay Vũ Văn Nghi và Lee (2015) [8]. Mục tiêu của nghiên cứu này là phát triển mô hình số để mô phỏng được sóng ven bờ bao gồm các quá trình động lực sóng nêu trên theo phương pháp thể tích hữu hạn kết hợp với thành phần phân tán Boussinesq giải theo sai phân hữu hạn nhằm đảm bảo tính ổn định cao, độ chính xác tốt và có thể áp dụng trên thực tiễn cho mô phỏng sóng ven bờ và động lực phía trong vùng sóng đổ. Trước tiên, hệ phương trình Boussinesq mở rộng có cải tiến tiêu tán năng lượng sóng đổ và cách giải được trình bày, sau đó mô hình số được mô phỏng cho bài toán sóng trên bãi thoải với điều kiện thí nghiệm vật lý nhằm đánh giá khả năng mô phỏng của mô hình số. Các mô phỏng cho bài toán phức tạp hơn với các điều kiện ở tỉ lệ thực được thực hiện và so sánh với kết quả từ mô hình MIKE 21. Cuối cùng là một ứng dụng thử nghiệm cho bài toán mô phỏng dòng chảy phát sinh do sóng giữa hai cồn ngầm trên bãi biển nhằm khẳng định việc mô phỏng được những điều kiện gần với thực tế. MÔ HÌNH TOÁN Hệ phương trình mô tả Xuất phát từ hệ phương trình Boussinesq do Madsen et al., (1997) [2] đề xuất, mô hình tính sóng ven bờ được phát triển cho mô phỏng cả khu vực sóng đổ và phía trong vùng sóng đổ với việc đưa vào các thành phần nhớt rối và tiêu tán sóng do ma sát. Các hệ phương trình được trình bày như sau: Phương trình bảo toàn khối lượng: 0         y Q x Q t yx (1) Phương trình bảo toàn động lượng theo phương x: 3 22 3 2 2 2 2 2 3 3 2 3 2 2 3 2 1 1 3 6 2 x y y yx x x x bx bx Q Q Q QQ Q Q h gd h h t x d y d x t x y y t xt x h h gh gh R x y x yx y x x y                                                                                                     (2) Phương trình bảo toàn động lượng theo phương y: 2 33 2 2 2 2 2 2 3 3 2 3 2 2 2 3 2 1 1 3 6 2 1 6 y x y y yx x Q Q Q Q QQ Qh gd h h t x d y d y t x y x t yt y h h gh gh y x x yx y x y y Qh h y                                                                                                    2 1 3 y yx by by Q R t x t y                  (3) Phung Dang Hieu et al. 16 Các số hạng thêm vào các phương trình nguyên thủy của Madsen et al., (1997) bao gồm: Thành phần mô tả trao đổi động năng nhớt rối do lớp cuộn xoáy sóng đổ gây ra:                                        vh xy uh yy uh xxh R eeebx   2 1 )( 1 (4)                                        vh xx uh yx vh yyh R eeeby   2 1 )( 1 (5) Thành phần mô tả tiêu tán năng lượng sóng do lớp cuộn xoáy gây ra: u d ghu t cBbrx      2 cho phương x (6) v d ghv t cBbry      2 cho phương y (7) Thành phần mô tả tiêu tán động năng do ma sát với đáy: 22 vuuC fx   , 22 vuvC fy   , 3/1 2 d gn C f  (8) Trong đó: η là dao động mặt nước; Qx là thông lượng theo phương x; Qy là thông lượng theo phương y; h là độ sâu nước yên tĩnh; d = (h + η) là độ sâu tổng cộng; g là gia tốc trọng trường. Tham số β được chọn là 1/15. Qx = ud, Qy = vd, với u là vận tốc trung bình độ sâu theo phương x, và v là vận tốc trung bình độ sâu theo phương y; n hệ số Manning được hiệu chỉnh theo tính chất nhám của bề mặt đáy. Hệ số nhớt rối do sóng đổ ve được xác định theo phương pháp của Kennedy et al., (2000) [4] đề xuất như sau: t hBe      )(2 ; δ =0,9–1,5 (9) * ** * * 2 2 0 1 1 tt ttt tt t tB                (10)   ( ) * * ( ) ( ) ( ) *0 0* , ,0 F t t I F I t t t t T t t t t T T                  (11) g h T 5*  ; ( ) 0,65It gh  ; ( ) 0,15Ft gh  (12) Với ηt * được Schaffer et al., (1993) [1] định nghĩa là tham số xác định sóng đổ; T* là khoảng thời gian chuyển đổi sóng đổ; t0 là thời điểm khi sóng đổ xảy ra; t – t0 là tuổi sóng đổ hay khoảng thời gian diễn ra sóng đổ; ηt (I) là giá trị xác định sóng bắt đầu đổ, giá trị của nó nằm trong khoảng hiệu chỉnh từ 0,35 gh đến 0,65 gh ; tham số ηt (F) là giá trị giới hạn cuối cùng của sóng đổ. Các điều kiện biên Biên mở sóng tới: Điều kiện biên nguồn tạo sóng hoặc biên bảng tạo sóng Stokes không phản xạ được thực hiện cho việc đưa dao động sóng tới vào tính toán. Bên cạnh đó, Numerical model for simulation of waves 17 điều kiện phát xạ tự do được cho thỏa mãn tại các biên mở thông thường và tại biên nguồn tạo sóng. Các sóng đều và sóng ngẫu nhiên không đều theo phổ sóng được quan tâm xây dựng cho mô phỏng. Biên cứng tường đứng, kè cứng: Được xác định theo điều kiện biên không thấm. Tức là vận tốc trực giao với biên bị triệt tiêu. Với thành phần tiếp tuyến với biên được áp dụng là điều kiện trượt không nhớt. Biên bãi biển dưới tác động của dâng và rút nước được áp dụng theo phương pháp thông lượng bảo toàn theo phương pháp thể tích hữu hạn sẽ được trình bày chi tiết trong phần rời rạc, giải số hệ phương trình. RỜI RẠC VÀ GIẢI SỐ HỆ PHƯƠNG TRÌNH Xử lý hệ phương trình dưới dạng phương trình nước nông và thành phần Boussinesq Khó khăn nhất trong việc giải quyết hệ phương trình Boussinesq truyền sóng trong vùng ven bờ đó là giải số thành phần phi tuyến tương tự như trong hệ phương trình nước nông truyền thống. Nếu xử lý theo cách thông thường theo phương pháp sai phân hữu hạn đối với thành phần này đòi hỏi phải có phép xấp xỉ dạng ngược dòng bậc nhất để đảm bảo ổn định số. Tuy nhiên, sử dụng phép xấp xỉ ngược dòng bậc một dẫn đến sai số rất lớn làm suy giảm sóng nhanh và khuếch tán số nghiêm trọng. Đối với phép xấp xỉ bậc cao cho phép đảm bảo độ chính xác và giảm khuếch tán số thì lại bị vấn đề không ổn định số và tràn số khi tồn tại sóng đổ và sóng tràn trên bãi. Chính vì vậy, việc xử lý số đối với thành phần này cần có cách phù hợp hơn. Phương pháp thể tích hữu hạn với các hàm giới hạn cho phương trình dạng bảo toàn đã chứng minh cho phép mô phỏng rất tốt thành phần phi tuyến với kết quả có độ chính xác cao và ổn định. Do đó, trong nghiên cứu ở đây sử dụng phương pháp phân tách lai giữa thể tích hữu hạn cho phần các số hạng nước nông truyền thống và sai phân hữu hạn bậc 2 cho thành phần Boussinesq phân tán sóng. Hệ các phương trình Boussinesq trong mô hình toán nêu ở trên có thể được nhìn nhận bao gồm phương trình nước nông truyền thống và thành phần hàm nguồn. Với hàm nguồn bao gồm các thành phần Boussinesq [2], trao đổi nhớt rối, tiêu tán năng lượng và ma sát đáy. Để thuận tiện cho việc rời rạc hóa theo phương pháp thể tích hữu hạn, các phương trình được chia thành hai bước giải số như sau: Bước 1: Giải số hệ phương trình nước nông không có hàm nguồn theo phương pháp thể tích hữu hạn với các thông lượng bảo toàn; Bước 2: Giải số thành phần hàm nguồn theo sai phân hữu hạn. Các phương trình (1), (2) và (3) được viết lại dưới dạng véctơ với hai hàm nguồn riêng biệt để phục vụ cho giải số kết hợp giữa phương pháp thể tích hữu hạn bảo toàn cho thành phần nước nông thuần túy và sai phân hữu hạn cho thành phần Boussinesq và tiêu tán sóng đổ. Các phương trình được viết thành dạng phân tách hai bước và dưới dạng véc tơ bảo toàn như sau: BoussS yxt          S GFU (13) Trong đó: U là véc tơ của các biến bảo toàn; F , G là các véctơ thông lượng tương ứng theo phương x và y ; S thành phần nguồn tương tự theo phương trình nước nông truyền thống. BoussS là thành phần nguồn do sóng đổ và số hạng phân tán của phương trình Boussinesq. 2 21 2 2 21 2 0 0 , , , ,x Bouss ux bx brx uy by bry y dud dv h du du gd duv gd B R x dv B Rduv dv gd h gd y                                                                   U F G S S (14) Phung Dang Hieu et al. 18 3 3 2 2 2 2 2 2 2 2 3 3 2 3 2 2 3 2 1 ( ) ( ) 1 ( ) 1 ( ) 1 ( ) 3 6 3 6 2 ux ud vd h vd h ud vd B h h h t x y y t x x t x t yt x h h gh gh x y x yx y x x y                                                                                              (15) 3 3 2 2 2 2 2 2 2 2 3 3 2 3 2 2 2 3 1 ( ) ( ) 1 ( ) 1 ( ) 1 ( ) 3 6 6 3 2 uy ud vd h ud h ud vd B h h h t x y x t y y t x t yt y h h gh gh y x x yx y x y y                                                                                              (16) Phương pháp giải hai bước Tách phương trình (13) thành hai bước với hai phương trình: Bước 1: S GFU          yxt (13a) Bước 2: BoussS t    *U (13b) Phương pháp thể tích hữu hạn được dựa trên luật bảo toàn vật chất áp đặt cho thể tích hữu hạn. Tích phân phương trình (13a) trên một ô lưới với việc áp dụng định lý Grin (Green’s theorem), cho ta:     ddnnd t yx SGF U )( (17) Trong đó:  là miền ô lưới;  là biên của miền  ; (nx, ny) là véctơ pháp tuyến hướng vào đường biên. Lấy tích phân theo thời gian phương trình (17) trong khoảng thời gian ∆t từ thời điểm t1 đến t2, ta có:        dtyxdtyx 12 ,,,, UU    ddtdnndt t t yx t t SGF 2 1 2 1 )( (18) Trong mô hình hiện tại sử dụng ô lưới đều với bước lưới ∆x, ∆y như thế, phương trình tích phân (18) với bước thời gian ∆t có thể được xấp xỉ với giá trị thời điểm thời gian ở giữa khoảng cho các thông lượng và hàm nguồn nên (18) được xấp xỉ bậc 2 thành.     2/1,2/1 2/1,2/1 2/1,2/1 ,2/12/1 ,2/1,1,         k ji k ji k ji k ji k ji k ji k ji t y t x t SGGFFUU (19) Trong đó: i, j là các chỉ số tại tâm của ô lưới; k ký hiệu bước thời gian hiện tại; các chỉ số một phần hai 2/1i , 2/1i và 2/1j , 2/1j dùng chỉ tại mặt phân cách giữa các ô lưới; và 2/1k chỉ trung bình giữa hai bước thời gian k và k + 1. Chú ý rằng, trong phương trình (19) các biến U và hàm nguồn S là các giá trị tại trung tâm ô lưới. Để giải phương trình (19), ta cần tính toán xác định các thông lượng số 2/1 ,2/1   k jiF , 2/1 ,2/1   k jiF và 2/1 2/1,   k jiG , 2/1 2/1,   k jiG tại các mặt phân cách các ô lưới. Trong nghiên cứu này, sử dụng sơ đồ Godunov-type scheme. Theo sơ đồ Godunov- type scheme, các hàm thông lượng số tại các mặt phân cách các ô lưới được xác định thông qua giải bài toán Riemann địa phương tại các mặt phân cách. Do nghiệm giải trực tiếp đối với bài toán Riemann 2 và 3 chiều chưa có, mô hình toán hiện tại sử dụng phương pháp sơ đồ tách bậc Numerical model for simulation of waves 19 hai của Strang (1968) [9] để giải tách phương trình (19) thành hai bước liên tiếp và được tích phân như sau: k ji tttk ji XYX , 2/2/1 , UU   (20) Với X và Y chỉ toán tử tích phân theo hướng x và y tương ứng. Phương trình theo phương x được tích phân trước với bước thời gian một nửa bước thời gian tích phân và tiếp theo đó là tích phân cả bước thời gian được thực hiện cho phương trình theo hướng y. Diễn tả như sau:   4/1,4/1 ,2/14/1 ,2/1,)2/1(, )( 22 *          k jix k ji k ji k ji k ji t x t SFFUU (21)   2/1,2/1 2/1,2/1 2/1,)2/1(,)1(, )( **          k jiy k ji k ji k ji k ji t y t SGGUU (22) Trong đó: Dấu (*) chỉ rằng các nghiệm phân tách trung gian; Sx, Sy là các nguồn theo hướng x và y. Tích phân theo hướng x trên khoảng nửa bước thời gian được tiếp tục tiến triển cho nửa bước thời gian tiếp theo để thu được nghiệm tại bước thời gian mới.   4/3,4/3 ,2/14/3 ,2/1)1(,1, )( 22 *          k jix k ji k ji k ji k ji t x t SFFUU (23) Các nghiệm từng phần k ji,U , *)2/1( , k jiU và *)1( , k jiU , dùng để cung cấp các thành phần thông lượng trong các phương trình (21), (22) và (23) thông qua giải bài toán Riemann một chiều. Ở đây phép xấp xỉ HLL cho nghiệm bài toán Riemann được sử dụng để xác định các thông lượng số trên mặt phân cách ô lưới. Để giải quyết trường hợp ô lưới chuyển khô ướt, một độ sâu giới hạn nhỏ được áp dụng để chuyển đổi giữa chúng (d = 10–5 m). Đối với phương trình (13b), các hàm nguồn được giải hiện theo các kết quả đã biết tại bước thời gian trước, riêng các thành phần Bux, Buy được sai phân trung tâm ẩn luân hướng thông thường đối với thành phần biến u và v để đảm bảo tăng ổn định cho mô hình số. Giới hạn bước thời gian phụ thuộc chủ yếu vào bước lưới không gian và sóng trọng lực trong miền tính. Điều kiện ổn định tương tự như các phương pháp sai phân thông thường đã sử dụng như sau: max min( , ) 0,5 ( ) x y t g h       (24) CÁC MÔ PHỎNG VÀ KẾT QUẢ Mô phỏng sóng trên bãi nghiêng Điều kiện thí nghiệm của Ting và Kirby (1996) [10] về sóng truyền trên bãi thoải có độ dốc 1/35 được đưa vào để thử nghiệm mô phỏng số và so sánh với kết quả thí nghiệm vật lý về phân bố độ cao sóng trên bãi nghiêng. Trong thí nghiệm này, sóng tới được cho dạng sóng Stokes bậc 2 có độ cao 12,5 m chu kỳ 2 s. Mô phỏng số được thực hiện cho 60 chu kỳ sóng đảm bảo sóng đủ kết hợp giữa sóng tới và sóng phản xạ cũng như tác động của bãi nghiêng lên chuyển động sóng. Độ cao sóng được tính toán là trung bình độ cao của 3 con sóng cuối. Kết quả phân bố độ cao được trình bày trên hình 1. Ta thấy, kết quả tính toán và số liệu thí nghiệm khá phù hợp. Đặc biệt tại điểm sóng đổ, độ cao sóng mô phỏng khá sát với thí nghiệm. Phía trong vùng sóng đổ, độ cao sóng mô phỏng thiên cao thực tế. Nguyên nhân do tiêu tán năng lượng chưa đủ lớn. Mặc dù vậy, gần bờ, độ cao sóng mô phỏng tiếp cận đến số liệu thí nghiệm, điều này cho phép tính toán tốt sóng leo bãi cũng như dòng phát sinh do sóng ở khu vực gần bờ. Phung Dang Hieu et al. 20 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 -5 0 5 10 15 x-xo (m) đ ộ c a o ( m ) Độ cao sóng (mô phỏng) Chân sóng (mô phỏng) Chân sóng (thí nghiệm) Độ cao sóng (thí nghiệm) Hình 1. So sánh độ cao sóng vùng sóng đổ giữa mô phỏng và thí nghiệm vật lý của Ting và Kirby (đường liền: kết quả mô phỏng; chấm tròn: độ cao sóng thí nghiệm, chấm vuông: mực nước chân sóng thí nghiệm) Mô phỏng sóng trên bãi thoải có đê chắn sóng nổi Miền địa hình bãi thoải 1/30 có độ sâu vùng chân bãi là 8 m, trên bãi có một đê chắn sóng dài 120 m đặt cách mép nước đường bờ 150 m. Sóng tới trực diện có độ cao 1,1 m, chu kỳ 6,3 s. Mô phỏng được thực hiện bằng mô hình MIKE 21-SW và mô hình số phát triển ở trên. Các kết quả phân bố độ cao sóng trên toàn miền, dòng chảy phát sinh do sóng toàn miền, phân bố độ cao sóng và dòng chảy sóng trên hai mặt cắt MC1 tại giữa miền tính từ bờ ra khơi và mặt cắt MC2 trên phần bãi thoải không có đê chắn sóng được trình bày trên các hình vẽ để so sánh. Hình 2. Phân bố độ cao sóng trên bãi có đê chắn sóng nổi (hình trái: kết quả MIKE 21-SW; hình phải: kết quả mô hình Boussinesq) Hình 2 cho thấy phân bố độ cao sóng mô phỏng trên miền tính toán. Đối với mô phỏng bằng MIKE 21-SW cho kết quả sóng phân bố khá đơn giản, độ cao có xu thế giảm dần khi độ sâu nông đi. Không thấy sự hiện diện của phản xạ sóng. Sóng phía sau đê chắn sóng giảm mạnh, thấy rõ vùng khuất sóng và sóng khúc xạ, nhiễu xạ. Đối với kết quả từ mô hình Boussinesq, độ cao sóng phân bố rất phức tạp nhìn rõ các vùng sóng kết hợp giữa sóng tới và Numerical model for simulation of waves 21 sóng phản xạ. Sóng bị tăng độ cao tại khu vực bãi thoải do hiệu ứng nước nông. Phía sau đê chắn sóng, tồn tại vùng khuất sóng, có sóng khúc xạ và nhiễu xạ đi vào. Do có trường sóng phức tạp hơn kết quả của MIKE 21-SW nên hệ thống dòng chảy sóng tính từ mô hình Boussinesq cũng sẽ có sự phân bố phức tạp hơn. Trên hình 3 cho thấy bức tranh phân bố của dòng chảy sóng ven bờ. Với hai mô hình, hệ thống dòng chảy dư do sóng chỉ tồn tại phức tạp xung quanh khu vực đê chắn sóng ven bờ. Cả hai mô hình đều cho một hệ thống dòng chảy đi từ bờ ra tới đê chắn sóng. Điều này rất phù hợp với lý thuyết và thực tế là phía sau đê chắn có dòng từ bờ ra mang vật chất nối đê với bờ tạo thành Tombolo. Dòng chảy sóng do MIKE 21 tính tạo ra hoàn lưu rõ ở hai phía đầu đê nhưng có vận tốc khá nhỏ, cực đại cỡ 0,55 m/s. Trong khi đó, dòng chảy sóng do mô hình Boussinesq mô phỏng cho vận tốc lớn hơn, cực đại cỡ 0,9 m/s. Như vậy về mặt vật lý, lý thuyết và thực tế, hai mô hình đều mô phỏng được hiện tượng dòng chảy nối bờ với vật cản phía ngoài để tạo Tombolo. Hình 3. Phân bố dòng chảy phát sinh do sóng trên bãi thoải có đê chắn sóng nổi (hình trái: kết quả MIKE 21; hình phải: kết quả mô hình Boussinesq) Kết quả mô phỏng giữa hai mô hình được xuất ra trên hai mặt cắt MC1 và MC2 để so sánh. Hình 4 so sánh tại mặt cắt MC1. Trên hình cho thấy, phân bố độ lớn dòng chảy dọc theo mặt cắt khá phù hợp với nhau về xu thế giữa hai mô hình. Tuy nhiên, mô hình Boussinesq cho kết quả mô phỏng vận tốc lớn hơn nhiều so với MIKE 21. Điều này có thể giải thích do MIKE 21 sử dụng ứng suất sóng tính theo mô đun SW không chứa đựng các kết hợp phức tạp của sóng ven bờ. Hiệu ứng sóng tăng độ cao do độ sâu giảm không được mô phỏng tốt trong SW. Hơn nữa, độ cao sóng trong SW bị ép giảm theo độ sâu mà không có tính đến hiệu ứng nước nông trước khi sóng đổ, chính điều này làm sóng bị mất năng lượng dẫn đến dòng chảy sóng được ước lượng thiên nhỏ. Điều này cũng được thấy rõ trên hình 4, đối với mô hình Boussinesq thấy rõ độ cao sóng tăng lên khi độ sâu giảm trên bãi thoải do hiệu ứng nước nông và có sự hiện diện của sóng phản xạ. Đối với độ cao mô phỏng bằng MIKE 21-SW không thấy hiện tượng tăng độ cao và sóng phản xạ. Tuy vậy, độ cao sóng nhiễu xạ và khúc xạ sau đê chắn sóng của cả hai mô hình là tương đương nhau. Trên hình 5 so sánh các kết quả tại mặt cắt MC2. Tại mặt cắt này, kết quả giữa hai mô hình có sự tương đồng tốt hơn về độ cao sóng và dòng chảy sóng. Tuy nhiên, phân bố độ Phung Dang Hieu et al. 22 cao sóng lần nữa cho thấy rõ, với MIKE 21- SW không thấy hiệu ứng nước nông và phản xạ sóng. Mô hình Boussinesq cho kết quả độ cao sóng cao hơn SW trước khi sóng đổ và ở vùng sát bờ. Điều này khẳng định với khu vực địa hình phức tạp, có công trình thì việc mô phỏng sóng bằng MIKE 21-SW sẽ gặp nhiều sai sót. Hình 4. Phân bố tốc độ dòng chảy sóng và độ cao sóng tại mặt cắt MC1 (so sánh giữa tính toán bằng mô hình MIKE 21và Boussinesq) Hình 5. Phân bố tốc độ dòng chảy sóng và độ cao sóng tại mặt cắt MC2 (so sánh giữa tính toán bằng mô hình MIKE21 và Boussinesq) Mô phỏng sóng trên bãi thoải có hai cồn ngầm Miền địa hình bãi thoải 1/30 tương tự như phần trên được thiết lập cho mô hình Boussinesq với hai cồn ngầm độ sâu đỉnh 0,5 m có độ dài 120 m cách nhau 60 m cách đều hai bên tạo ra khoảng trống giữa hai cồn. Dạng địa hình này khá hay gặp tại các bãi biển thực tế. Mô phỏng được thực hiện cho sóng tới trực diện có độ cao 1,1 m chu kỳ 6,3 s nhằm xem xét hệ thống dòng chảy sóng xuất hiện thế nào, liệu có dòng Rip nguy hiểm giữa hai cồn ngầm không như các khuyến cáo của các nghiên cứu trước đây. Lưới tính được thiết lập chi tiết với độ phân giải 1 m × 1 m. Kết quả trình bày trên hình 6 cho thấy, sóng bị phản xạ và đổ mạnh tại hai cồn ngầm và tạo ra hai vùng khuất sóng phía sau có độ cao sóng nhỏ. Tuy nhiên, tại khoảng trống giữa hai cồn ngầm, độ cao sóng tăng cao do hội tụ tia sóng và sóng kết hợp với nhau. Nhìn trên hình 6 (hình dưới) ta thấy dòng chảy sóng ven bờ khá mạnh và đặc biệt tạo ra dòng rút tại cửa hở khoảng trống giữa hai cồn ngầm. Chính dòng này đã làm sóng hội tụ và dồn lại làm tăng độ cao sóng. Dòng rút này khá lớn lên đến cỡ hơn 1 m/s rất nguy hiểm cho người tắm trên bãi nếu ra gần phía cửa trống sẽ bị cuốn vào dòng này lôi ra ngoài. Như vậy, bằng mô phỏng số bởi mô hình Boussinesq đã cho thấy được sự xuất hiện của dòng rút giữa khoảng trỗng giữa hai cồn ngầm gần bờ. Dòng rút là một trong những nguyên nhân rất nguy hiểm gây ra các vụ đuối nước khi tắm tại các bãi biển mùa hè. Chính vì vậy, việc nghiên cứu các hiện tượng dòng rút do sóng tại các dạng địa hình khác nhau rất có ý nghĩa cho cảnh báo và giảm tai nạn đuối nước tại các bãi biển. Numerical model for simulation of waves 23 Hình 6. Phân bố độ cao sóng (hình trái) và dòng chảy phát sinh do sóng (hình phải) trên bãi thoải có hai cồn ngầm KẾT LUẬN Bài báo trình bày kết quả nghiên cứu phát triển thành công mô hình số dựa trên hệ phương trình Boussinesq mở rộng và kết hợp cải tiến tiêu tán năng lượng do sóng đổ. Mô hình được xây dựng dựa trên sự kết hợp giải số theo phương pháp thể tích hữu hạn với sơ đồ TVD và sai phân trung tâm cho thành phần phân tán sóng Boussinesq. Các mô phỏng số ban đầu so sánh với số liệu thí nghiệm cũng như với kết quả mô phỏng bằng mô hình thương mại MIKE 21 cho thấy sự phù hợp tốt và có sự tương đồng định tính cao. Mô hình Boussinesq phát triển trong nghiên cứu này có khả năng mô phỏng khá linh hoạt các bài toán sóng vùng ven bờ có sự hiện diện của công trình, có sự phản xạ sóng. Việc nghiên cứu mô phỏng dòng chảy sóng ven bờ khi có công trình và có cồn ngầm cho thấy mô hình cho phép mô phỏng được qui luật của các quá trình vật lý do sóng gây ra như phản xạ, dòng ven, dòng rút. Từ kết quả nghiên cứu cũng cho thấy, mô hình MIKE 21-SW nếu mô phỏng cho khu vực có địa hình phức tạp, có hiện diện công trình, có sự phản xạ, kết hợp sóng phức tạp thì sẽ gặp rất nhiều hạn chế dẫn đến sai sót trong kết quả mô phỏng. Lời cảm ơn: Bài báo được hoàn thành dưới sự hỗ trợ của đề tài TNMT.2016.06.09. Nhóm tác giả xin trân trọng cảm ơn. TÀI LIỆU THAM KHẢO [1] Schäffer, H. A., Madsen, P. A., and Deigaard, R., 1993. A Boussinesq model for waves breaking in shallow water. Coastal engineering, 20(3–4), 185–202. [2] Madsen, P. A., Sørensen, O. R., and Schäffer, H. A., 1997. Surf zone dynamics simulated by a Boussinesq type model. Part I. Model description and cross-shore motion of regular waves. Coastal Engineering, 32(4), 255–287. [3] Madsen, P. A., Sørensen, O. R., and Schäffer, H. A., 1997. Surf zone dynamics simulated by a Boussinesq type model. Part II: Surf beat and swash oscillations for wave groups and irregular waves. Coastal Engineering, 32(4), 289–319. [4] Kennedy, A. B., Chen, Q., Kirby, J. T., and Dalrymple, R. A., 2000. Boussinesq modeling of wave transformation, breaking, and runup. I: 1D. Journal of waterway, port, coastal, and ocean engineering, 126(1), 39–47. [5] Wei, G., Kirby, J. T., Grilli, S. T., and Subramanya, R., 1995. A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves. Journal of Fluid Mechanics, 294, 71–92. Phung Dang Hieu et al. 24 [6] Phung Dang Hieu, 2011. A numerical model for Tsunami propagation and runup: A case study in the Bien Dong sea. Journal of Science, Natural Sciences and Technology, VNU, 27(1S), 96–108. [7] Thuy, N. B., Nandasena, N. A. K., Dang, V. H., Kim, S., Hien, N. X., Hole, L. R., and Thai, T. H., 2017. Effect of river vegetation with timber piling on ship wave attenuation: investigation by field survey and numerical modeling. Ocean Engineering, 129, 37–45. [8] Van Nghi, V. U., and Changhoon, L. E. E., 2015. Solitary wave interaction with porous structures. Procedia Engineering, 116, 834–841. [9] Strang, G., 1968. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5(3), 506–517. [10] Ting, F. C., and Kirby, J. T., 1996. Dynamics of surf-zone turbulence in a spilling breaker. Coastal Engineering, 27(3–4), 131–160.