Nghiên cứu xây dựng mô hình thủy động lực ba chiều tính toán trường dòng chảy xung quanh các công trình thủy lực phức tạp

                                                                                                                                                                                                 '''''' '''''' '''''' '''''' '''''' '''''' '''''' '''''' '''''' 1 ww zz wv yz wu xz vw zy vv yy vu xy uw zx uv yx uu xx ww zz wv yz wu xz vw zy vv yy vu xy uw zx uv yx uu xx ww zz wv yz wu xz vw zy vv yy vu xy uw zx uv yx uu xx J zzyzxz yzyyxy zxyxxx zzyzxz yzyyxy zxyxxx zzyzxz yzyyxy zxyxxx                                                                                  (7) Để khép kín hệ phương trình nói trên, cần phải thêm các phương trình mô tả ứng suất rối Reynold. Nhằm mục đích dự báo được dòng chảy bất ổn định ba chiều, chú trọng đến xoáy Karman và xoáy quy mô lớn hình thành trong kênh hỗn hợp (có bãi) vốn có quy mô thời gian lớn hơn là các thành phần dao động rối, một mô hình k- phi tuyến chứa các hàm thực nghiệm đã được sử dụng và cho thấy được khả năng mô phỏng tương đối tốt (Hosoda và cộng sự 1999, Kimura và Hosoda 2003, Kimura và nnk 2004). Mô hình rối phi tuyến này, bao gồm thành phần bậc hai, đã được chứng tỏ là tương đương với mô hình nhớt rối hiện phi tuyến trong mô hình ứng suất Reynold đại số (Pope 1975, Gastki và Speziale 1993). Năng lượng rối động học k và hệ số khuếch tán rối  nhận được từ các phương trình vận chuyển sau đây trong hệ tọa độ khớp biên di động:                                       m lm k t ll i j ljij G j j k J gD J u xJ uu J kUU J k t  '' (8) 23                                                   m lm k t l l i j lji l j G j j k J gD kJ C u xJ uu k C J UU Jt           2 2 '' (9) trong đó: 00,1k , 30,1 , 44,1lC và 92,12 C . Tensor ứng suất Reynold có thể được biểu diễn như sau:                    3 1 3 1 3 2 ''    ijijtijm j i m m i j m t ji SSCD k k u x u x Duu (10) trong đó: r j r i ij u x u S     1                 r i j r r j i r ij x u x u x u x u S 2 1 2 j r i r ij x u x u S     3 Hệ số nhớt rối Dt và các thông số Cµ, C1, C2, C3 có thể được biểu diễn theo Nagata và nnk (2005) như sau:  2k CDt          209,01 3,0 ;09,0min M C 21 01,01 4,0 M C   02 C 23 01,01 13,0 M C    Thông số M trong các phương trình ở trên được xác định thông qua tham số sức căng S và số xoáy Ω: 24   ;max SM 2 2 1                   i j j i x u x uk S  2 2 1                   i j j i x u x uk  2.2. Thuật toán và phương pháp giải 2.2.1. Mô hình nguyên gốc của Hosoda a) Rời rạc hóa hệ phương trình – sơ đồ sai phân. Phương trình động lượng và phương trình chuyển động của k và ε được rời rạc hóa bằng phương pháp thể tích hữu hạn trên hệ thống lưới so le. Để thỏa mãn tính liên tục cục bộ, nghiên cứu này sử dụng phương pháp Mark-and-Cell đơn giản hóa (HSMAC) (Hirt and Cook 1972) vì cả vận tốc và áp suất có thể tính đồng thời trong khi trong phương pháp MAC nguyên bản (Harlow and Welch, 1965) thì phương trình Poission của biến áp suất phải được giải trước. Yếu tố đối lưu được rời rạc hóa bằng phép nội suy ngược bậc 2 (quadratic upstream interpolation) sơ đồ động học đối lưu QUICK (Leonard, 1979) và yếu tố khuếch tán được rời rạc hóa bằng sơ đồ sai phân trung tâm. Trong khi xác định vận tốc cho các bước thời gian tiếp theo khi đã biết vận tốc và áp suất cục bộ, phương trình liên tục được sử dụng để tính áp suất thông qua phép lặp số (SOLA), đồng thời giá trị mực nước cũng được xác định theo điều kiện động lực ở mặt thoáng. Khi giá trị mực nước này được xác định, thì sơ đồ lưới trong mặt phẳng thẳng đứng được xây dựng bằng cách biến đổi khoảng cách lưới. Lưu ý không có sự thay đổi lưới trong mặt phẳng (x1−x2) trong nghiên cứu này. Các điều kiện biên của động năng dòng chảy rối kb, và tốc độ tiêu hao động năng dòng chảy rối εb được xác định ở biên cố định sử dụng hàm biên (wall function) sau (Rodi 1993): 𝑘𝑏 = 𝑢∗ 2 𝐶𝜇 ; 𝜀𝑏 = 𝑢∗ 3 𝜅Δ𝑙 (11) trong đó Δl là khoảng cách từ biên. Có một vài phương pháp (ví dụ, Naot và 25 Rodi 1982; Gibson và Rodi 1989; Naot và cộng sự. 1993) được đưa ra để tính tác động tắt dần của hiện tượng chảy rối do sự xuất hiện của các mặt thoáng. Vận tốc trượt xác định từ profin vận tốc loga được sử dụng để xác định kb và εb. Trong nghiên cứu hiện nay, một số phương pháp tính toán được Kimura và Hosoda sử dụng: (1) độ nhớt rối được biến đổi bằng hàm tắt dần Isawa và Hosoda đưa ra (1990); và (2) tốc độ tiêu tán ở mặt thoáng được xác định bằng công thức Siguyama và cộng sự đề xuất (1997), tương tự với quan niệm của Naot và Rodi (1982). Rõ ràng cường độ rối theo hướng thẳng đứng giảm theo hướng mặt thoáng. Để tính toán tác động này, nhân độ nhớt rối với hàm tắt dần fd được Iwasa và Hosoda đưa ra (1990): 𝑓𝑑 = 1− 𝑒𝑥𝑝 −𝐵 𝑍𝑠−𝜒 3 𝜀𝑠 𝑘𝑠 3/2 (12) trong đó zs là mực nước; B là hằng số (=10); ks là động năng của dòng rối ở bề mặt; εs là tốc độ tiêu tán động năng dòng rối ở mặt thoáng; và 𝜒 3 là tọa độ thẳng đứng. Nhằm chứng tỏ tính ổn định và hợp lý của mô hình số trị, một số các mô phỏng đã thực hiện và so sánh với các giá trị đo đạc trong phòng thí nghiệm đã được công bố, cụ thể là mô hình nguyên bản của Hosoda được kiểm định với thí nghiệm của Munita và Shimizu (1994) về dòng chảy xung quanh công trình kè mỏ hàn vuông góc. Mực nước được xây dựng như biên hạ lưu trong kiểm tra thực nghiệm, được giữ không đổi. Điều kiện biên trên, giá trị lưu lượng tương ứng với thực nghiệm được đưa ra và phân bố như vậy để có profin vận tốc loga theo phương thẳng đứng. Kiểm định này đã được công bố trong công trình của Nagata và nnk. (2005). Ở đây xin được nhắc lại một số kết quả cơ bản nhằm chứng tỏ độ tin cậy của mô hình này. b) Kết quả kiểm định với thí nghiệm của Munita và Shimizu (1994) Nghiên cứu của Munita và Shimizu được thực hiện trong phòng thí nghiệm nhằm mục đích đo đạc chi tiết trường dòng chảy xung quanh kè mỏ hàn sử dụng các máy đo vận tốc Laser Doppler và chú trọng đến các trường phân bố vận tốc tại các 26 mặt cắt ngang ngay phía sau khu vực kè mỏ hàn. Kênh nghiên cứu là kênh hình chữ nhật thẳng, có kích thước là rộng 0.4m, dài 5m, trên đó thiết đặt một mỏ hàn vuông góc phía bờ trái theo hướng dòng chảy có chiều dài kè là 0.2m và chiều rộng kè 0.04 m (hình 1). Kè chảy hoàn toàn trong trạng thái không ngập. Kênh thí nghiệm là kênh xi măng đáy cố định. Các mặt cắt đo đạc số liệu vận tốc dòng chảy là các mặt cắt dọc (A-A) và các mặt cắt ngang (B-B), (C-C) và (D-D). Các điều kiện về dòng chảy như sau: - Dạng lòng dẫn: Lòng cứng - Lưu lượng: 0.00187 m3/s - Bề rộng kênh: 0.4 m - Độ dốc đáy: 1/1.000 - Độ sâu dòng chảy: 0.07 m - Vận tốc trung bình đầu kênh: 0.067 m/s -Số Froude: 0.08 Hình 1. Mô tả thí nghiệm của Munita và Shimizu (1994) trên mặt bằng (hình có tính minh họa không chính xác về tỷ lệ) 27 Trong mô phỏng bằng mô hình, lưới tính trên mặt bằng là lưới vuông góc thay đổi dần, lưới có độ phân giải lớn 1cm x 1cm ở xung quanh khu vực kè và giảm dần về hai phía đầu và cuối kênh (1cm x 10cm) để tiết kiệm thời gian tính toán (hình 2). y ( m ) x (m) Hình 2. Minh họa lưới tính toán sử dụng trong mô hình Miền tính toán sử dụng có kích cỡ là 75 x 45 x 13 nút lưới, trong đó trục thẳng đứng được chia gồm 13 lớp từ đáy đến mặt nước. Một số kết quả tính toán được biểu diễn trên các hình 3 đến 5 và được so sánh với các số liệu đo đạc (hình 6). Kết quả tính toán cho thấy rằng mô hình đã mô phỏng được trường vận tốc có tính chất ba chiều rõ rệt xung quanh khu vực có kè mỏ hàn vuông góc, dòng chảy hướng từ trên xuống xuất hiện ở phía mặt trước của kè (hình 3) và đã hình thành nên một xoáy ngay phía sát mặt đáy. Mặt khác, xét mặt cắt ngang ngay phía sau kè, dòng chảy hướng từ phía trái sang phía phải (hình 4) ở phía gần đáy do ảnh hưởng của dòng chảy đi xuống phía trước mặt kè và đồng thời do tại mũi kè thì dòng chảy đã bị thu hẹp đột ngột. Phía xa hơn ở hạ lưu (mặt cắt C- C) một xoáy nhỏ theo chiều kim đồng hồ cũng đã hình thành và đã được chính Munita và Shimizu (1994) xem là một tính chất nổi bật của dòng chảy dọc theo đường hoàn lưu từ sau kè, trong khi đó một xoáy khác ngược chiều kim đồng hồ lại xuất hiện ở gần sát mặt nước. Các so sánh giữa số liệu tính toán và đo đạc trên các hình (6 a-c) cũng khẳng định lại các tính chất cơ bản nói trên. Sự phù hợp giữa số liệu thực đo và tính toán là có thể chấp nhận được. 28 Z (m ) x(m) Mặt trước kè mỏ hàn Hình 3. Trường vận tốc theo mặt cắt dọc A-A thượng lưu kè mỏ hàn Z (m ) y(m) Hình 4. Trường vận tốc trên mặt cắt ngang B-B Z (m ) y(m) Hình 5. Trường vận tốc trên mặt cắt ngang C-C 29 Z (m ) Bờ trái x(m) Bờ phải  o Số liệu đo đạc  Số liệu tính toán Z (m ) Bờ trái x(m) Bờ phải Z (m ) Bờ trái y(m) Bờ phải Hình 6. So sánh lưu tốc tính toán và thực đo trên các mặt cắt ngang a- Mặt B-B b- Mặt C-C c- Mặt D-D 30 2.2.2. Chỉnh sửa mô hình để có thể chạy với trường hợp kè hoàn lưu và kè chảy ngập Mô hình nguyên gốc của GS Hosoda và các cộng sự mặc dầu đã tỏ rõ ưu thế trong việc mô tả trường dòng chảy ba chiều xung quanh công trình chỉnh trị như kè mỏ hàn và cũng đã được các tác giả chứng minh tính ứng dụng với các hố xói quanh các trụ cầu (Nagata và nnk, 2005). Tuy nhiên, so với thực tiễn của Việt Nam, mô hình vẫn còn một số các hạn chế và cần thiết phải có các thay đổi đặc biệt là trong trường hợp ứng dụng cho công trình kè hoàn lưu (tấm hướng dòng) khi mà thuật toán nhận diện công trình phụ thuộc vào mực nước và vào từng lớp tọa độ theo phương thẳng đứng. Hạn chế của mô hình là chỉ áp dụng được cho các công trình liền khối và chảy ở trạng thái không ngập, nhưng ở hầu hết các sông ngòi Việt Nam, do hạn chế về kinh phí và mục tiêu ứng dụng (chỉ phục vụ nâng cao mực nước kiệt) nên về mùa lũ mực nước thường cao hơn nhiều so với đỉnh kè nên kè mỏ hàn (và các công trình tương tự) chuyển sang chế độ chảy ngập hoặc chuyển tiếp giữa các chế độ chảy trong cùng một trận lũ. Với mục tiêu, chỉnh sửa mô hình nhằm mô phỏng trường dòng chảy ba chiều xung quanh công trình kè hoàn lưu, nghiên cứu này đã phân tích các đặc trưng hình thái đoạn sông, kích thước và quy mô các công trình và đi đến giả thiết như sau: So với các công trình kè mỏ hàn khác, kè hoàn lưu thường có dạng bản mỏng, độ rộng (chiều dày khoảng 0,5 m) không đáng kể so với chiều dài kè (khoảng 50-80m) cũng như so với quy mô đoạn sông nghiên cứu (thông thường bán kính cong 700-800 m, chiều rộng mặt cắt ngang ~ 500 m). Do vậy, có thể xem rằng dòng chảy dưới đáy kè là dòng chảy không áp trong trường hợp mực nước trong sông cao hơn ngưỡng dưới của kè hoàn lưu. Xuất phát từ giả thiết đó, sau khi xác định vị trí của các công trình (cụ thể là kè hoàn lưu) trong không gian 3 chiều (i, j, k), các điều kiện tương tự như điều kiện biên đối với tường cứng đã được áp dụng để mô tả ảnh hưởng của công trình. Các điều kiện biên này được đưa trực tiếp vào các thủ tục giải phương trình động lượng. 31 2.2.2. Sơ đồ khối mô hình Toàn bộ các thủ tục tính toán trong mô hình được trình bày trong sơ đồ khối mô tả như trong hình 7 sau: Hình 7. Sơ đồ khối mô hình thủy động lực ba chiều 32 2.2.3. Chương trình tính toán Giải thích cụ thể về chương trình tính toán như sau: a) Chương trình chính Nội dung mã nguồn của chương trình chính (main program) như sau: implicit real*8(a-h,o-z) parameter (mx=325,ny=36,lz=17) common/par1/ n,nend,nmod,isr,ibr,idown,iconv,kemdl,nst,nst2,ibo & ,iost,ioend,jost,joend,kbed,inqb,nd,np,mp & ,ist,mpls,mend,m common/par2/ dgu,dyi,dtu,roh,g,dt,zdown,q,slp,bl,al,an,anu & ,akap,as,sigk,sige,cmu,cmun,ce1,ce2,aks,ar,alp,omg & ,atvd,btvd,dm,sigm,amus,amuk,ramd,hsd,akl,ak2 & ,thr,w0,ak0,ae0,dtd common/par3/ pai,cd0,sha,asd2,asd3,cm,one,ustc,tstac,sdm & ,aaa,znen,ztai,asnen0,astai0 common /grd/ x(mx,ny,lz),y(mx,ny,lz),z(mx,ny,lz) common/hyd1/quc(mx,ny,lz),qvc(mx,ny,lz),qwc(mx,ny,lz), & quc1(mx,ny,lz),qvc1(mx,ny,lz),qwc1(mx,ny,lz), & uc(mx,ny,lz),vc(mx,ny,lz),wc(mx,ny,lz), & uc1(mx,ny,lz),vc1(mx,ny,lz),wc1(mx,ny,lz), & u(mx,ny,lz),v(mx,ny,lz),w(mx,ny,lz), & p(mx,ny,lz),pt(mx,ny,lz),pb(mx,ny) common/zsb/ zs(mx,ny),zb(mx,ny),wl(mx,ny),zbed(mx,ny) common/obj/ mob(mx,ny),mobg(mx,ny),nmob(mx,ny) dimension tarray(4) 33 call openf keizok=0 call prmt n=0 call anhgrd write(6,*) x(iost,jost,1),y(iost,jost,1) write(6,*) x(ioend,jost,1),y(ioend,jost,1) write(6,*) x(iost,joend,1),y(iost,joend,1) write(6,*) x(ioend,joend,1),y(ioend,joend,1) call metrix call metri2 call intcnd call regrd do n=1,nend do i=1,mx do j=1,ny do k=1,lz if(u(i,j,k).gt.1..or.v(i,j,k).gt.1..or.w(i,j,k).gt.1.) then write(73,*) '---error---',n,i,j,k,u(i,j,k),v(i,j,k),w(i,j,k) end if end do end do end do call metrix call bound if(n.ge.100) then call calpro 34 call calak call calae end if call fri call visco call calquc call calqvc call calqwc call hsmac call metri2 call conv call surf call regrd if(mod(n,20000).eq.0) call out if(mod(n,20000).eq.0) call out2 end do stop end b) Các chương trình con Chương trình gồm có những chương trình con cụ thể như sau: Openf: mở các file đầu vào và đầu ra mà chương trình làm việc (trong đó có file fdz13.dat là file tỷ lệ phân chia lưới theo phương thẳng đứng và file out2.dat là file kết quả suất ra.) Prmt: xác định các thông số Anhgrd: đọc vào file số liệu địa hình và tạo lưới tính cho mô hình Readc: bắt đầu tính toán từ số liệu đầu ra 35 Metrix: tính toán ma trận Metri2: tính toán ma trận Intcnd: điều kiện ban đầu Regrd: sinh lại lưới Conv: tính toán véc-tơ trong hệ tọa độ Cartesian tại mỗi điểm lưới (i, j, k), (i+1/2, j, k), (i, j+1/2, k), (i, j, k+1/2), (i, j+1/2, k+1/2), (i+1/2, j, k+1/2), (i+1/2, j+1/2, k). Bound: điều kiện biên Calpro: số hạng j i ji x u uu    Calak: phương trình k Calae: phương trình ε Fri: tính toán ứng suất tiếp đáy và vận tốc ma sát Visco: tính toán các số hạng ứng suất tiếp trong phương trình mô men. Calquc: thành phần phản biến của véc-tơ vận tốc U1 Calqvc: thành phần phản biến của véc-tơ vận tốc U2 Calqwc: thành phần phản biến của véc-tơ vận tốc U3 Hsmac: tính toán áp suất bằng thủ tục lặp Surf: tính cao trình mặt nước Out, out2: xuất kết quả ra file đầu ra. c) Các thông số dgu,dyi,dtu: ,, dt: t (bước thời gian) n: bước 36 nend: bước cuối cùngfinal step nmod: số liệu được ghi lại tại mỗi nmod bước g: gia tốc trọng trường slp: độ dốc q: lưu lượng zdown: cao trình mặt nước tại hạ lưu iost: i, j, k là số hiệu lưới theo chiều dọc, ngang sông và theo phương thẳng đứng và iost là số hiệu bắt đầu của công trình theo chiều dọc sông ioend: điểm kết thúc của công trình theo chiều dọc sông jost, joend: số hiệu lưới bắt đầu và kết thúc theo hướng ngang sông. an: hằng số Manning akap: hằng số Karman anu: nhớt rối động học phân tử ak0, ae0: giá trị ban đầu của k và  roh: mật độ nước as: As theo luật loga (đáy nhẵn), ar : Ar theo luật loga (đáy thô) sigk: k (=1.0), sige :  (=1.3) cmu: c cmun: giá trị nhỏ nhất của c (=0.09), và c được tính bằng hàm của S và  ce1: 1, ce2: 2 ibo: điều kiện biên thượng lưu (=0, lưu lượng tại i=1 được cho trước) idown: mặt nước tại hạ lưu là một hằng số (=0), gradient =0 (=1) iconv: lựa chọn sơ đồ tính các số hạng đối lưu 37 nst: k,  được gộp vào sau khi n>nst nst2: lựa chọn turbulence model isr : lựa chọn về ma sát thành bên (=0, trơn), (=1, nhám) isr : lựa chọn về ma sát đáy (=0, trơn), (=1, nhám) aks: hệ số nhám tương đương ks omg : hằng số được sử dụng trong tính toán áp suất kemdl: mô hình chuẩn (=0), mô hình k-ε phi tuyến (=1) atvd,btvd: hằng số trong sơ đồ TVD amus,amuk : s, k ramd:  hsd: kd = 0.8 akl,ak2 = 0.85, 0.7 d) Các biến sử dụng trong chương trình x(mx,ny,lz),y(mx,ny,lz),z(mx,ny,lz): Tọa độ Cartesian x0(mx,ny,lz),y0(mx,ny,lz),z0(mx,ny,lz): tọa độ Cartesian trong bước trước xgint(mx,ny,lz),ygint(mx,ny,lz): x và y ban đầu tại điểm lưới xcnt(mx,ny,lz),ycnt(mx,ny,lz),zcnt(mx,ny,lz) : x, y và z tại tâm phần tử xcu(mx,ny,lz),ycu(mx,ny,lz),zcu(mx,ny,lz) : x, y và z tại vị trí của U xcv(mx,ny,lz),ycv(mx,ny,lz),zcv(mx,ny,lz) : x, y và z tại vị trí của V xcw(mx,ny,lz),ycw(mx,ny,lz),zcw(mx,ny,lz) : x, y và z tại vị trí của W xcnt1(mx,ny,lz),ycnt1(mx,ny,lz),zcnt1(mx,ny,lz),xcw1(mx,ny,lz),ycw1(mx,n y,lz),zcw1(mx,ny,lz): bước thời gian tiếp theo quc(mx,ny,lz),qvc(mx,ny,lz),qwc(mx,ny,lz) : U i /J 38 quc1(mx,ny,lz),qvc1(mx,ny,lz),qwc1(mx,ny,lz) : U i /J tại bước tiếp theo uc(mx,ny,lz),vc(mx,ny,lz),wc(mx,ny,lz) : U i (thành phần phản biến) uc1(mx,ny,lz),vc1(mx,ny,lz),wc1(mx,ny,lz) : U i tại bước tiếp theo u(mx,ny,lz),v(mx,ny,lz),w(mx,ny,lz) : ui (thành phần Cartesian) ui(mx,ny,lz),vi(mx,ny,lz),wi(mx,ny,lz), uj(mx,ny,lz),vj(mx,ny,lz),wj(mx,ny,lz), uk(mx,ny,lz),vk(mx,ny,lz),wk(mx,ny,lz): ui tại các điểm khác nhau umeng(mx,ny,lz),vmeng(mx,ny,lz),wmeng(mx,ny,lz), umeny(mx,ny,lz),vmeny(mx,ny,lz),wmeny(mx,ny,lz), ument(mx,ny,lz),vment(mx,ny,lz),wment(mx,ny,lz): ui tại các điểm khác nhau của tâm ô lưới p(mx,ny,lz),pt(mx,ny,lz),pb(mx,ny): áp suất dp(mx,ny,lz),dpt(mx,ny,lz): thay đổi theo thời gian của áp suất tại mỗi bước zs(mx,ny): cao trình mặt nước tại điểm lưới zb(mx,ny): cao trình đáy tại điểm lưới wl(mx,ny),zbed(mx,ny): mặt nước và cao trình đáy tại tâm ô lưới fdz(mx,ny,lz),fdn(mx,ny,lz),fdzk(lz): phân bố của ô lưới theo phương thẳng đứng zs0(mx,ny),zb0(mx,ny),p0(mx,ny,lz) : mặt nước, cao trình đáy và áp suất ban đâu zbint(mx,ny),zgint(mx,ny,lz): cao trình đáy ban đầu tại tâm ô lưới và điểm lưới mob(mx,ny): ô có công trình (=1) 39 mobg(mx,ny) : xác định công trình (4 ô lân cận) (i-1, j-1), (i-1, j), (i, j-1), (i, j) nmob(mx,ny) : xác định công trình (9 ô lân cận) guxyac(mx,ny,lz),guyyac(mx,ny,lz),guzyac(mx,ny,lz) : x/J, y/J, z/J yixyac(mx,ny,lz),yiyyac(mx,ny,lz),yizyac(mx,ny,lz) : x/J, y/J, z/J tuxyac(mx,ny,lz),tuyyac(mx,ny,lz),tuzyac(mx,ny,lz) : x/J, y/J, z/J yac(mx,ny,lz) : J yacuc(mx,ny,lz),yacvc(mx,ny,lz),yacwc(mx,ny,lz): Jacobian tại vị trí U, V, W gux(mx,ny,lz),guy(mx,ny,lz),guz(mx,ny,lz): x, y, z (điểm xác định là giống với U) yix(mx,ny,lz),yiy(mx,ny,lz),yiz(mx,ny,lz) : x, y, z (điểm xác định là giống với V) tux(mx,ny,lz),tuy(mx,ny,lz),tuz(mx,ny,lz) : x, y, z (điểm xác định là giống với W) g11(mx,ny,lz), g12(mx,ny,lz), g13(mx,ny,lz), g21(mx,ny,lz), g22(mx,ny,lz), g23(mx,ny,lz), g31(mx,ny,lz), g32(mx,ny,lz), g33(mx,ny,lz) : g 11 -g 33 xgu(mx,ny,lz),xyi(mx,ny,lz),xtu(mx,ny,lz) : x, x, x ygu(mx,ny,lz),yyi(mx,ny,lz),ytu(mx,ny,lz) : y, y, y zgu(mx,ny,lz),zyi(mx,ny,lz),ztu(mx,ny,lz) : z, z, z guxvc(mx,ny,lz),guyvc(mx,ny,lz),guzvc(mx,ny,lz) :x, y, z (điểm xác định là giống với V) guxwc(mx,ny,lz),guywc(mx,ny,lz),guzwc(mx,ny,lz) :x, y, z (điểm xác định là giống với W) 40 yixuc(mx,ny,lz),yiyuc(mx,ny,lz),yizuc(mx,ny,lz) : x, y, z (điểm xác định là giống với U) yixwc(mx,ny,lz),yiywc(mx,ny,lz),yizwc(mx,ny,lz) :x, y, z (điểm xác định là giống với W) tuxuc(mx,ny,lz),tuyuc(mx,ny,lz),tuzuc(mx,ny,lz) x, y, z (điểm xác định là giống với U) tuxvc(mx,ny,lz),tuyvc(mx,ny,lz),tuzvc(mx,ny,lz) : x, y, z (điểm xác định là giống với V) guxyac1(mx,ny,lz),guyyac1(mx,ny,lz),guzyac1(mx,ny,lz), yixyac1(mx,ny,lz),yiyyac1(mx,ny,lz),yizyac1(mx,ny,lz), tuxyac1(mx,ny,lz),tuyyac1(mx,ny,lz),tuzyac1(mx,ny,lz) yac1(mx,ny,lz), yacuc1(mx,ny,lz),yacvc1(mx,ny,lz),yacwc1(mx,ny,lz), gux1(mx,ny,lz),guy1(mx,ny,lz),guz1(mx,ny,lz), yix1(mx,ny,lz),yiy1(mx,ny,lz),yiz1(mx,ny,lz), tux1(mx,ny,lz),tuy1(mx,ny,lz),tuz1(mx,ny,lz) : bước tiếp theo qug(mx,ny,lz),qvg(mx,ny,lz),qwg(mx,ny,lz) : U i G/J thành phần phản biến của vận tốc lưới ugrid(mx,ny,lz),vgrid(mx,ny,lz),wgrid(mx,ny,lz) : thành phần Cartesian của vận tốc lưới. qwgi(mx,ny,lz),qwgj(mx,ny,lz),qugk(mx,ny,lz),qvgk(mx,ny,lz): U i G/J tại các điểm khác nhau buu(mx,ny,lz),bvv(mx,ny,lz),bww(mx,ny,lz), buv(mx,ny,lz),buw(mx,ny,lz),bvw(mx,ny,lz): ứng suất Reynolds visgu(mx,ny,lz),visyi(mx,ny,lz),vistu(mx,ny,lz): các số hạng nhớt rối theo các phương , ,  41 ak(mx,ny,lz) : k ae(mx,ny,lz) :  ad(mx,ny,lz) : = +Dt/k d(mx,ny,lz) : = Dt (=C*k 2 /) ak1(mx,ny,lz),ae1(mx,ny,lz): bước tiếp theo akn0(mx,ny,lz),aen0(mx,ny,lz): bước trước pro(mx,ny,lz) : các số hạng nhân taurx(mx,lz),taury(mx,lz),taurz(mx,lz): ứng suất tiếp từ thành bên phải taulx(mx,lz),tauly(mx,lz),taulz(mx,lz): ứng suất tiếp từ thành bên trái taubx(mx,ny),tauby(mx,ny),taubz(mx,ny): ứng suất tiếp từ đáy taurox(mx,lz),tauroy(mx,lz),tauroz(mx,lz): ứng suất tiếp trong trường hợp mà công trình được đặt ở mặt phải taulox(mx,lz),tauloy(mx,lz),tauloz(mx,lz): ứng suất tiếp trong trường hợp mà công trình được đặt ở mặt trái taufox(ny,lz),taufoy(ny,lz),taufoz(ny,lz): ứng suất tiếp trong trường hợp công trình được đặt ở mặt trước taubox(ny,lz),tauboy(ny,lz),tauboz(ny,lz): ứng suất tiếp trong trường hợp công trình được đặt ở mặt sau ustro(mx,lz),ustlo(mx,lz),ustfo(ny,lz),ustbo(ny,lz): vận tốc ma sát trong trường hợp mà công trình được đặt ở mặt phải, trái, trước, sau tương ứng ustb(mx,ny): vận tốc ma sát đáy umapb(mx,ny),vmapb(mx,ny),wmapb(mx,ny): vận tốc theo từng phương song song với đáy ubx(mx,ny),uby(mx,ny),ubz(mx,ny): véc-tơ vận tốc sát đáy 42 bnx(mx,ny), bny(mx,ny),bnz(mx,ny): véc-tơ chuẩn đơn vị tại điểm lưới bnxc(mx,ny),bnyc(mx,ny), bnzc(mx,ny): véc-tơ chuẩn đơn vị tại tâm ô lưới bpxx(mx,ny),bpxz(mx,ny) : =                       bpxz bpxx p b b b 0 sin 0 cos 1 1 1   bpyy(mx,ny),bpyz(mx,ny): =                       bpyz bpyyp b bb 0 sin cos 0 2 22   ustbg(mx,ny),ustsd(mx,ny,msd): vận tốc ma sát tại điểm lưới, vận tốc ma sát trung bình ustbgo(mx,ny,msd): vận tốc ma sát rmdkjn(mx,ny,msd): độ dài bước gs1(mx,ny,msd),gs2(mx,ny,msd),gs3(mx,ny,msd), gs4(mx,ny,msd): tỷ lệ phân bố cho điểm lưới mhos(mx,ny,msd) : đánh dấu cho xác định hướng (tức là từ i hay i+1) 43 Chương 3 – MỘT SỐ KẾT QUẢ THỬ NGHIỆM MÔ HÌNH VỚI CÁC CÔNG TRÌNH THỦY LỰC PHỨC TẠP 3.1. Thử nghiệm mô hình với thí nghiệm đoạn sông thẳng có công trình Để chứng tỏ khả năng mô phỏng trường dòng chảy ba chiều xung quanh kè hoàn lưu, đặc biệt chú trọng đến các xoáy và dòng hoàn lưu thượng và hạ lưu kè cũng như các dòng chảy vòng xung quanh khu vực công trình, nghiên cứu này đã thử nghiệm mô hình khi có công trình kè hoàn lưu vuông góc với bờ. 3.1.1. Thử nghiệm với thí nghiệm số Thực nghiệm số được tiến hành với trường hợp kênh thẳng hình chữ nhật, dài 10m rộng 3m (hình 8), trên đó thiết lập kè mỏ hàn có chiều dày nhỏ (có thể bỏ qua) và xét các trường hợp: - Kè không ngập (giống như các thực nghiệm số trước đây) - Kè mỏ hàn chảy ngập - Kè hoàn lưu Hình 8. Mô tả kênh và công trình thực nghiệm số a) Trường hợp kè không ngập: Thực nghiệm số được thực hiện nhằm chứng tỏ những thủ thuật thay đổi đã áp dụng (đưa các điều kiện biên tường cứng vào vị trí có công trình) không làm thay đổi bức tranh dòng chảy so với mô hình nguyên gốc vốn đã được kiểm định trước 3m 1.5m 44 đây. Các kết quả mô phỏng biểu diễn trên các hình 9 và 10 cho thấy, khi dòng chảy gặp phải kè không ngập, tại mặt kè cả phía thượng và hạ lưu sẽ xuất hiện dòng đi xuống và do đó làm phát sinh xoáy, lõi xoáy có vận tốc rất bé (giống như đã quan sát thấy trong thí nghiệm của Munita và Shimizu (1994). Mặt khác, trên mặt phẳng nằm ngang cũng quan sát thấy các xoáy ở chân kè phía hạ lưu, nguyên nhân dẫn đến hiện tượng bồi tụ chân kè. Hình 9. Véc tơ vận tốc trên mặt cắt dọc tại mũi kè (a) và tại thân kè (b) Hình 10. Véc tơ vận tốc trên mặt ngang tại độ sâu (a) giữa thân kè và (b) trên mặt nước a) b) a) b) 45 b) Trường hợp kè ngập: Trước khi mô phỏng kè chảy ngập, nghiên cứu này đã tăng lưu lượng và mực nước để trường hợp a) trên đây trở thành kè ngập. Hình 11. Véc tơ vận tốc trên mặt cắt dọc tại (a) mũi kè và (b) giữa thân kè Hình 12. Véc tơ vận tốc trên mặt ngang tại độ sâu (a) giữa thân kè và (b) trên mặt ngang sát đỉnh kè c) Trường hợp kè hoàn lưu: Kè mỏ hàn dạng hoàn lưu được nghiên cứu ở đây mới thuộc dạng đơn giản, chỉ có hướng vuông góc với bờ trong kênh thẳng (trong khi thực tế là kè hoàn lưu hướng dòng chữ L, và gốc kè cũng không vuông góc với bờ). a) b) a) b) 46 Hình 13. Véc tơ vận tốc trên mặt cắt dọc tại mũi kè (a) và giữa thân kè (b) Hình 14. Véc tơ vận tốc trên mặt ngang tại độ sâu (a) giữa thân kè và (b) trên mặt ngang sát đỉnh kè a) b) a) b) 47 3.1.2. Thử nghiệm với thí nghiệm vật lý của Tominaga a) Thí nghiệm vật lý của Tominaga và cộng sự. Hình 15. Các thiết đặt của mô hình thí nghiệm. Mô hình vật lý dùng để kiểm chứng mô hình 3D ở đây là của Tominaga và cộng sự (2000) xây dựng trong phòng thí nghiệm. Trong mô hình vật lý này (Hình 15) có thiết đặt hai kè ở bờ trái của một kênh hở hình chữ nhật dài 8 m, rộng 0.3 m và độ dốc đáy (i) bằng 1/2000. Độ dài phương ngang (l) của mỗi kè là 0.05 m; độ dày của mỗi kè (b) là 0.02 m; khoảng cách giữa chúng (S) là 0.1 m. Độ cao của các kè (d) là 0.02, 0.03, 0.04, 0.05, 0.06 m (trong các trường hợp chảy ngập) và 0.08 m (trong trường hợp chảy không ngập) tương ứng với các trường hợp SD2, SD3, SD4, SD5, SD6 và SD8. Lưu lượng đầu vào phía thượng lưu (Q) và độ sâu cột nước phía hạ lưu được giữ không đổi bằng 4.1 l/s và 0.08 m. b) Kết quả thử nghiệm với mô hình vật lý của Tominaga và cộng sự. Tất cả các điều kiện thí nghiệm ở trên đã được mô phỏng sử dụng mô hình 3D đã xây dựng với một lưới gồm 321 x 61 x 17 nút lưới tương ứng theo các phương dọc, ngang và thẳng đứng (Hình 16) có độ phân giải cao hơn ở các ô lưới gần sát với các công trình. 48 Hình 16. Lưới tính toán: a) Nhìn theo mặt bằng (x-y) b) Nhìn theo mặt bên (x-z) Sau khi ứng dụng mô hình 3D đã xây dựng với lưới tính toán như trên cho các trường hợp thí nghiệm khác nhau của Tominaga, thu được một số kết quả cụ thể như sau: Trong hình 17 biểu thị hình ảnh các véc tơ vận tốc dòng chảy so sánh giữa thực đo và thí nghiệm tại khu vực giữa hai kè. Qua đó ta có thể thấy các xoáy quy mô lớn theo phương thẳng đứng giữa các kè đã được tái tạo lại bằng các mô phỏng rất tốt, kể cả những tác động của việc tăng độ cao của các kè, làm tăng kích thước của các xoáy, tăng vận tốc và chuyển dịch của các lõi xoáy về phía hạ lưu. Các mô phỏng cũng đã tái hiện lại động thái dòng chảy ở phía thượng lưu của kè đầu tiên, trong đó dòng chảy được tách thành hai phần: một phần hướng lên trên từ đỉnh kè và một phần hướng xuống dưới, hình thành nên xoáy hình móng ngựa và là nguyên nhân tiềm năng gây xói cục bộ. Những xoáy hướng lên này được gia tăng khi tăng chiều cao của kè nhưng gần như biến mất trong trường hợp chảy không ngập (SD8). Hiện tượng này không được quan trắc thấy ở phía trước của kè thứ hai. 49 Thí nghiệm Mô phỏng Hình 17. Trường véc tơ vận tốc trên mặt cắt Y. Trái: các kết quả đo đạc, Phải: các kết quả mô phỏng Các xoáy theo phương ngang cũng đã được tái tạo lại rất tốt bằng mô hình như ta có thể thấy trong Hình 18 biểu thị trường vận tốc trên mặt bằng nằm ngang tại độ sâu bằng 1/2 chiều cao kè. Các xoáy theo phương ngang xuất hiện trong tất cả 50 các trường hợp mô phỏng và mạnh hơn trong các trường hợp SD6 và SD8. Đây chính là nguyên nhân tiềm năng gây ra xói cục bộ nghiêm trọng hơn ở phía mũi kè. Thí nghiệm Mô phỏng Hình 18. Trường véc tơ vận tốc trên mặt Z. Trái: các kết quả đo đạc; Phải: các kết quả mô phỏng 51 Trường vận tốc trong hình 19 cho thấy rằng dòng nghịch được sinh ra gần bờ kênh với dòng chảy tràn và các xoáy giữa hai kè khá mạnh (Hình 19 – trên). Có thể đây là một nguyên nhân gây nên xói lở bờ gần chỗ tiếp giáp kè. Nhưng điều này không được quan sát thấy ở khu vực gần mũi kè (Hình 19 – dưới). Tại cùng một thời điểm, một dòng hướng xuống được thể hiện ở phía trước của mũi kè thứ hai (Hình 19). Điều này cũng có thể là nguyên nhân gây nên xói cục bộ tại mũi của kè thứ hai. Những đặc trưng này đã được khẳng định bằng việc xuất hiện một xoáy lớn gần tầng đáy như được chỉ ra trong Hình 20. Hơn nữa những sự tương đồng về phân bố vận tốc giữa các kết quả đo đạc và mô phỏng cũng được thẻ hiện trong Hình 21 và Hình 22. Thí nghiệm Mô phỏng Hình 19. Trường véc tơ vận tốc trên mặt Y. Trái: các kết quả đo đạc; Phải: các kết quả mô phỏng 52 Thí nghiệm Mô phỏng Hình 20. Trường véc tơ vận tốc trên mặt Z. Trái: các kết quả đo đạc; Phải: các kết quả mô phỏng 53 Hình 21. Phân bố vận tốc Ux. Trái: phân bố thẳng đứng tại x=0.06m, y=0.025m; Phải: phân bố hướng ngang tại x=0.06, z= d/2 0 1 2 3 4 5 6 7 8 9 10 -10 -5 0 5 10 15 20 25 Z (c m ) Ux (cm/s) SD2: Y2 (X = 6 cm) Cal Obs 0 1 2 3 4 5 6 7 8 9 10 -10 -5 0 5 10 15 20 25 Y ( cm ) Ux (cm/s) SD2: Z2 (X = 6 cm) Cal Obs 0 1 2 3 4 5 6 7 8 9 10 -10 -5 0 5 10 15 20 25 Z (c m ) Ux (cm/s) SD4: Y2 (X = 6 cm) Cal Obs 0 1 2 3 4 5 6 7 8 9 10 -10 -5 0 5 10 15 20 25 Y ( cm ) Ux (cm/s) SD4: Z2 (X = 6 cm) Cal Obs 0 1 2 3 4 5 6 7 8 9 10 -10 -5 0 5 10 15 20 25 Z (c m ) Ux (cm/s) SD6: Y2 (X = 6 cm) Cal Obs 0 1 2 3 4 5 6 7 8 9 10 -10 -5 0 5 10 15 20 25 Y ( cm ) Ux (cm/s) SD6: Z2 (X = 6 cm) Cal Obs 54 Hình 22. Phân bố theo phương thẳng đứng của vận tốc Uz x=-0.01, y=0.025 (trái) và phân bố theo phương ngang của vận tốc Uy x=-0.01, z=d/2 (phải) 0 1 2 3 4 5 6 7 8 9 10 -4 -2 0 2 4 6 8 Z (c m ) Uz (cm/s) SD2: Y2 (X = -1 cm) Cal Obs 0 1 2 3 4 5 6 7 8 9 10 -4 -2 0 2 4 6 8 Y ( cm ) Uy (cm/s) SD2: Z2 (X = -1 cm) Cal Obs 0 1 2 3 4 5 6 7 8 9 10 -4 -2 0 2 4 6 8 Z (c m ) Uz (cm/s) SD4: Y2 (X = -1 cm) Cal Obs 0 1 2 3 4 5 6 7 8 9 10 -4 -2 0 2 4 6 8 Y ( cm ) Uy (cm/s) SD4: Z2 (X = -1 cm) Cal Obs 0 1 2 3 4 5 6 7 8 9 10 -4 -2 0 2 4 6 8 10 12 14 Y ( cm ) Uy (cm/s) SD6: Z2(X = -1 cm) Cal Obs 55 Tất cả những đặc trưng trên của cấu trúc dòng chảy đã được mô phỏng khá phù hợp với những kết quả đo đạc của Tominaga và nnk (2000). Điều này khẳng định khả năng của mô hình số đã xây dựng trong mô phỏng đặc tính ba chiều của dòng chảy xung quanh các kè chảy ngập. 3.2. Thử nghiệm mô hình với thí nghiệm vật lý đoạn sông cong có công trình 3.2.1. Thí nghiệm vật lý Trong nghiên cứu này tiến hành tính toán thử nghiệm mô hình động lực học dòng chảy ba chiều với kết quả thí nghiệm vật lý đã được tiến hành trong khuôn khổ nội dung đề tài KC.08.14/06-10 do GS.TS Lương Phương Hậu chủ trì thực hiện [1]. Máng thí nghiệm được thiết kế với chiều rộng B = 120 cm, chiều cao 50 cm. Mô hình được thiết kế thí nghiệm liên hoàn với 3 khúc cong nối tiếp nhau, đoạn quá độ dài bằng 5B = 600 cm, các khúc cong đều là những cung tròn với góc tâm 600. Hình 23. Mặt bằng mô hình thí nghiệm vật lý + Đoạn cong thứ nhất: R = 5B = 600 cm; + Đoạn cong thứ hai: R = 4B = 480 cm; + Đoạn cong thứ ba: R = 3B = 360 cm. Như vậy, mô hình dài tổng cộng 38 m. Hình 23 thể hiện bản vẽ mặt bằng của mô hình. Đây là mô hình thủy lực máng cong đầu tiên được thực hiện ở Việt Nam. Trong nội dung thực hiện của đề tài KC.08.09-06/10, đã tiến hành rất nhiều 56 các thí nghiệm mô hình với nhiều trường hợp khác nhau. Trong khuôn khổ của nghiên cứu này, chỉ xét ở khúc sông cong đầu tiên và chỉ xét hai trường hợp (một trường hợp không có tấm hướng dòng và một trường hợp có tấm hướng dòng) để kiểm nghiệm với mô hình toán thủy lực ba chiều đã xây dựng. Hai trường hợp này đều có lưu lượng thí nghiệm là Q = 20 l/s. Trong trường hợp có tấm hướng dòng: - Các tấm hướng dòng được bố trí tại các vị trí 2T, 3T, 4T, 5T, 6T như trong hình 24. Hình 24. Mặt bằng mô hình tại khúc sông cong thứ nhất - Mỗi tấm hướng dòng có cấu tạo như trong hình 25, đặt nghiêng so với bờ (trong trường hợp này α = 1350, θ = 600). Bề dầy tấm hướng dòng H = 20 cm, có một khoảng trống kể từ cạnh dưới của tấm hướng dòng và đáy mô hình (độ mở tấm hướng dòng) a = 0.2 H = 4 cm. Hình 25. Cấu tạo tấm hướng dòng θ 57 Hình 26 biểu thị mặt bằng lưới tính toán bằng mô hình thủy động lực ba chiều và vị trí các tấm hướng dòng. Kích thước của lưới tính theo chiều dòng chảy, chiều ngang và chiều thẳng đứng tương ứng là 86, 60, 12. Hình 26. Mặt bằng lưới tính và vị trí các công trình. 3.2.2. Kết quả kiểm nghiệm Hình 27 biểu thị so sánh các giá trị mực nước giữa kết quả thí nghiệm vật lý và kết quả mô phỏng bằng mô hình thủy động lực ba chiều trong trường hợp có bố trí công trình. Ta có thể thấy kết quả so sánh về mực nước này là khá phù hợp với nhau. Hình 28 biểu thị so sánh phân bố vận tốc (Uave = u2 + v2) mô phỏng bằng mô hình thủy động lực ba chiều trong trường hợp có công trình và trường hợp không có công trình. Qua đó ta có thể thấy, trong trường hợp không có công trình vận tốc tại bờ ngoài cao hơn dọc theo đoạn sông cong, điều này phù hợp với các nghiên cứu về đoạn sông cong. Trong trường hợp có công trình, vận tốc ở bờ trong cao hơn so với bờ ngoài. Điều này là do hướng dòng chảy bị chuyển sang bờ trong do tác động của các công trình. 58 a) Mặt cắt 2 (i = 31) b) Mặt cắt 3 (i = 41) c) Mặt cắt 4 (i = 51) d) Mặt cắt 5 (i = 61) e) Mặt cắt 6 (i = 71) f) Mặt cắt 7(i = 86) Hình 27. So sánh mực nước giữa thí nghiệm và mô phỏng trong trường hợp có công trình tại các mặt cắt 59 a) Mặt cắt 2T (i=36) b) Mặt cắt 3T (i=46) c) Mặt cắt 5T (i=66) d) Mặt cắt 6T (i=76) Open = Không có công trình, Cal. = Có công trình, (B) = đáy, (M) = giữa, (S) = mặt Hình 28. So sánh phân bố vận tốc (Uave = u2 + v2) trong trường hợp có công trình và không có công trình Hình 29 thể hiện sự so sánh về sự phân bố vận tốc (Uave = u2 + v2) giữa thí nghiệm vật lý và kết quả mô phỏng bằng mô hình thủy động lực ba chiều trong trường hợp có công trình. Qua đó có thể thấy rằng: - Vận tốc ở bờ trong là cao hơn do tác động hướng dòng của công trình - Vận tốc tại tầng đáy của bờ ngoài (đường màu xanh nước biển) cao hơn so với tầng giữa và tầng mặt của phía bờ ngoài (đường màu nâu và đường màu xanh lá cây). - Sự phù hợp giữa kết quả số và thực nghiệm không thực sự được thỏa mãn. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.5 0.7 0.9 1.1 1.3 Ua ve (m /s ) distance (m) Open (B) Open (M) Open (S) Cal. (B) Cal. (M) Cal. (S) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.7 0.9 1.1 1.3 Ua ve (m /s ) distance (m) Open (B) Open (M) Open (S) Cal. (B) Cal. (M) Cal. (S) 0 0.05 0.1 0.15 0.2 0.25 0.3 .3 0.4 0.45 0.5 0.5 0.7 0.9 1.1 1.3 Ua ve (m /s ) distance (m) Open (B) Open (M) Open (S) Cal. (B) Cal. (M) Cal. (S) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.5 0.7 0.9 1.1 1.3 Ua ve (m /s ) distance (m) Open (B) Open (M) Open (S) Cal. (B) Cal. (M) Cal. (S) 60 a) Mặt cắt 2 (i=31) b) Mặt cắt 2T (i=36) c) Mặt cắt 3 (i=41) d) Mặt cắt 3T (i=46) e) Mặt cắt 4T (i=56) f) Mặt cắt 6 (i=71) Hình 29. So sánh về sự phân bố vận tốc tốc (Uave = u2 + v2) giữa thí nghiệm vật lý và mô phỏng bằng mô hình trong trường hợp có công trình Qua Hình 30 thể hiện mặt bằng các véc-tơ vận tốc ở trên mặt và dưới đáy (được mô phỏng bằng mô hình) ta thấy rằng nước gần đáy chảy dọc theo hướng 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.5 0.7 0.9 1.1 1.3 U av e (m /s ) distance (m) Exp. (B) Exp. (M) Exp. (S) Cal. (B) Cal. (M) Cal. (S) 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.5 0.7 0.9 1.1 1.3 Ua ve (m /s ) distance (m) Exp. (B) Exp. (M) Exp. (S) Cal. (B) Cal. (M) Cal. (S) 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 40 0.5 0.7 0.9 1.1 1.3 U av e (m /s ) distance (m) Exp. (B) Exp. (M) Exp. (S) Cal. (B) Cal. (M) Cal. (S) 0.00 0.05 0.10 0.15 0.20 25 0.30 0.35 0.40 0.45 0.5 0.7 0.9 1.1 1.3 U av e (m /s ) distance (m) Exp. (B) Exp. (M) Exp. (S) Cal. (B) Cal. (M) Cal. (S) 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 4 0.5 0.7 0.9 1.1 1.3 U av e (m /s ) distance (m) Exp. (B) Exp. (M) Exp. (S) Cal. (B) Cal. (M) Cal. (S) 0.00 0.05 0.10 0.15 0.20 0.25 .3 0.35 0.40 0.45 0.50 0.5 0.7 0.9 1.1 1.3 U av e (m /s ) distance (m) Exp. (B) Exp. (M) Exp. (S) Cal. (B) Cal. (M) Cal. (S) 61 lòng dẫn, còn hướng dòng chảy tại lớp trên mặt thì thay đổi, hướng tới bờ trong do tác động của các công trình. Như vậy ảnh hưởng của các công trình đã được tái hiện lại bằng mô hình thủy động lực ba chiều. Các véc-tơ vận tốc ở trên mặt Các véc-tơ vận tốc ở dưới đáy Hình 30. Mặt bằng các véc-tơ vận tốc ở trên mặt và dưới đáy Mặc dầu mới chỉ là những bước thử nghiệm ban đầu, song những kết quả trên đây cho thấy khả năng phát triển và ứng dụng mô hình 3 chiều cho mô phỏng trường dòng chảy xung quanh kè mỏ hàn và ở đoạn sông cong. 62 KẾT LUẬN Những nội dung cơ bản đã thực hiện được trong nghiên cứu này gồm có: - Tìm hiểu và tổng quan lại tình hình nghiên cứu trong và ngoài nước về mô hình thủy động lực ba chiều mô phỏng trường dòng chảy xung quanh các công trình thủy lực phức tạp và trong đoạn sông cong. Từ đó cho thấy hiện vẫn chưa có mô hình thủy động lực 3 chiều thể hiện được tính phức tạp của dòng chảy xung quanh công trình và có tính ứng dụng đối với các công trình dạng kè mỏ hàn ngập và kè hoàn lưu, đặc biệt là khi các công trình này được đặt ở vị trí các đoạn sông cong. - Nghiên cứu tìm hiểu chi tiết về mô hình thủy thạch động lực ba chiều của Hosoda và nnk. Trên cơ sở đó chỉnh sửa lại mã nguồn mô hình để có thể ứng dụng cho trường hợp kè chảy ngập và kè hoàn lưu, có khả năng phù hợp hơn đối với điều kiện thực tế ở Việt Nam. - Kiểm nghiệm mô hình đã xây dựng với thí nghiệm trong đoạn sông thẳng có công trình: (1) Thí nghiệm số công trình kè hoàn lưu vuông góc với bờ; (2) Thí nghiệm vật lý của Tominaga và nnk (2000) đối với hai kè (chảy ngập hoặc không ngập) bố trí liên tiếp nhau. Kết quả cho thấy mô hình đã mô phỏng lại một cách khá phù hợp các đặc tính ba chiều của dòng chảy khi có tác động của các công trình. - Kiểm nghiệm mô hình đã xây dựng với thí nghiệm vật lý trong đoạn sông cong: (1) Đoạn sông cong không có công trình, (2) Đoạn sông cong có bố trí các công trình hướng dòng. Kết quả kiểm nghiệm cho thấy mô hình xây dựng về mặt định tính đã tái hiện lại được các hiện tượng dòng chảy tuy nhiên sự phù hợp về mặt định lượng giữa thí nghiệm và mô phỏng chưa thực sự được thỏa mãn. Cần có những nghiên cứu, cải thiện hơn nữa mô hình. 63 TÀI LIỆU THAM KHẢO Tiếng Việt 1. Lương Phương Hậu (2010), Nghiên cứu giải pháp khoa học, công nghệ cho hệ thống công trình chỉnh trị sông trên các đoạn trọng điểm vùng đồng bằng Bắc Bộ và Nam Bộ, Báo cáo tổng kết đề tài KC08.14/06-10 do GS. Lương Phương Hậu chủ trì, Hà Nội. 2. Nguyễn Thọ Sáo (2008), Động lực học chất lỏng tính toán, Giáo trình biên dịch, trường Đại học Khoa học Tự nhiên, Hà Nội. 3. Phạm Văn Tiến (2012), Ứng dụng mô hình (VNU/MDEC) tính toán chế độ thủy động lực và vận chuyển trầm tích vùng cửa sông ven biển Hải Phòng, Luận văn thạc sĩ khoa học ngành Thủy văn học, Trường Đại học Khoa học Tự nhiên, Hà Nội. Tiếng Anh 4. Ahmed, F., and Rajaratnam, N. (1998). “Flow around bridge piers”, J. Hydraul. Eng., 124(3), 288–300. 5. Bosch, G., and Rodi, W. (1998). “Simulation of vortex shedding past a square cylinder with different turbulence models”, Int. J. Numer. Methods Fluids, 28, 601–616. 6. Breusers, H. N. C., Nicollet, G., and Shen, H. W. (1977). “Local scour around cylindrical piers.” J. Hydraul. Res., 15(3), pp. 211–252. 7. Brookes, A., Knight, S. S., and Shields, F. D., Jr. (1996). “Habitat enhancement”. River channel restoration, A. Brookes and F. D. Shields, Jr., eds., Wiley, Chichester, U.K., pp. 103–126. 8. Chrisohoides, A., Sotiropoulos, F., and Sturm, T. W. (2003). “Coherent structures in flat-bed abutment flow: Computational fluid dynamics simulations and experiments.” J. Hydraul. Eng., 129(3), pp. 177–186. 9. Dargahi, B. (1990). “Controlling mechanism of local scouring.” J. Hydraul. Eng., 116(10), pp. 1197–1214. 10. Dey, S. (1997). “Local scour at piers, Part I: A review of developments of research.” Int. J. Sediment Res., 12(2), pp. 23–46. 11. Dey, S., Bose, S. K., and Sastry, G. L. N. (1995). “Clear water scour at circular piers: A model.” J. Hydraul. Eng., 121(12), pp. 869–876. 12. Elliott, K. R., and Baker, C. J. (1985). “Effect of pier spacing on scour around bridge piers.” J. Hydraul. Eng., 111(7), pp. 1105–1109. 13. Ettema, R., Mostafa, E. A., Melville, B. W., and Yassin, A. A. (1998). “Local scour at skewed piers.” J. Hydraul. Eng., 124(7), pp. 756–759. 64 14. Franke, R., and Rodi, W. (1993). “Calculation of vortex shedding past a square cylinder with various turbulence models.” Selected Papers from the 8th Int. Symp. on Turbulent Shear Flows, Munich, Germany, September 9–11, 1991, F. Durst, et al., eds., Springer-Verlag, Berlin, pp. 189–204. 15. Garde, R. J., Subramanya, K., and Nambudripad, K. D. (1961). “Study of scour around spur-dikes.” J. Hydraul. Div., Am. Soc. Civ. Eng., 87 (6), pp. 23–37. 16. Gatski, T. B., and Speziale, C. G. (1993). “On explicit algebraic stress models for complex turbulent flows.” J. Fluid Mech., 254, pp. 59–78. 17. Ge, L., and Sotiropoulos, F. (2005). “3D unsteady RANS modeling of complex hydraulic engineering flows. I: Numerical model.” J. Hydraul. Eng., 131(9), pp. 800–808. 18. Ge, L., Lee, S. O., Sotiropoulos, F., and Sturm, T. (2005). “3D unsteady Reynolds-averaged Navier-Stokes modeling of complex hydraulic engineering flows. II: Model validation and flow physics.” J. Hydraul. Eng., 131(9), pp. 809–820. 19. Gill, M. A. (1972). “Erosion of sand beds around spur dikes.” J. Hydraul. Div., Am. Soc. Civ. Eng., 98(9), 1587–1602. 20. Hosoda, T., Sakurai, T., Kimura, I., and Muramoto, Y. (1999). “3-D computations of compound open channel flows with horizontal vortices and secondary currents by means of non-linear k-epsilon model.” J. Hydrosci. Hydr. Eng., 17(2), 87–96. 21. Jia, Y., and Wang, S. S. Y. (1993). “3D numerical simulation of flow near a spur dike.” Proc., 1st Int. Conf. on Hydro-Sci. and -Engineering, Washington, D.C., 2150–2156. 22. Jia, Y., and Wang, S. S. Y. (1996). “A modeling approach to predict local scour around spur dike-like structures.” Proc., 6th Federal Interagency Sedimentation Conf., Las Vegas, Nev., II-90–97. 23. Jia, Y., and Wang, S. S. Y. (1999). “Numerical model for channel flow and morphological change studies.” J. Hydraul. Eng., 125(9), 924–933. 24. Jain, S. C. (1981). “Maximum clear-water scour around circular piers.” J. Hydraul. Div., Am. Soc. Civ. Eng., 107(5), 611–626. 25. Kimura, I., and Hosoda, T. (2003). “A non-linear k- model with realizability for prediction of flows around bluff bodies.” Int. J. Numer. Methods Fluids, 42, 813–837. 26. Kimura, I., Hosoda, T., Onda, S., and Tominaga, A. (2004). “3D numerical analysis of unsteady flow structures around inclined spur dikes by means of a non-linear k-_ model.” Shallow Flows, Jirka and Uijttewaal, eds., Selected Papers of the International Symposium on Shallow Flows, 16–18 June 2003, Delft, The Netherlands, Taylor & Francis Group, London, 651–660. 65 27. Kothyari, U. C., Garde, R. J., and Raju, K. G. R. (1992). “Temporal variation of scour around circular bridge piers.” J. Hydraul. Eng., 118(8), 1091–1106. 28. Kuhnle, R. A., Alonso, C. V., and Shields, F. D., Jr. (1999). “Geometry of scour holes associated with 90° spur dikes.” J. Hydraul. Eng., 125(9), 972–978. 29. Kwan, T. F., and Melville, B. W. (1994). “Local scour and flow measurements at bridge abutments.” J. Hydraul. Res., 32(5), 661–673. 30. Lai, Y. G., Weber, L. J., and Patel, V. C. (2003a). “Nonhydrostatic threedimensional model for hydraulic flow simulation. I: Formulation and verification.” J. Hydraul. Eng., 129(3), 196–205. 31. Lai, Y. G., Weber, L. J., and Patel, V. C. (2003b). “Nonhydrostatic threedimensional model for hydraulic flow simulation. II: Validation and application.” J. Hydraul. Eng., 129(3), 206–214. 32. Laursen, E. M. (1963). “An analysis of relief bridge scour.” J. Hydraul. Div., Am. Soc. Civ. Eng., 89(3), 93–118. 33. Lim, S. Y. (1997). “Equilibrium clear-water scour around an abutment.” J. Hydraul. Eng., 123(3), 237–243. 34. Mayerle, R., Toro, F. M., and Wang, S. S. Y. (1995). “Verification of a three- dimensional numerical model simulation of the flow in the vicinity of spur dikes.” J. Hydraul. Res., 33(2), 243–256. 35. Melville, B. W. (1975). “Local scour at bridge site.” Rep. No. 117, School of Engineering, The Univ. of Auckland, New Zealand. 36. Melville, B. W. (1992). “Local scour at bridge abutments.” J. Hydraul. Eng., 118(4), 615–631. 37. Melville, B. W. (1997). “Pier and abutment scour: Integrated approach.” J. Hydraul. Eng., 123(2), 125–136. 38. Melville, B. W., and Chiew, Y. M. _1999_. “Time scale for local scour at bridge piers.” J. Hydraul. Eng., 125_1_, 59–65. 39. Melville, B. W., and Raudkivi, A. J. (1977). “Flow characteristics in local scour at bridge piers.” J. Hydraul. Res., 15(4), 373–380. 40. Melville, B. W., and Raudkivi, A. J. (1996). “Effects of foundation geometry on bridge pier scour.” J. Hydraul. Eng., 122(4), 203–209. 41. Melville, B. W., and Sutherland, A. J. (1988). “Design method for local scour at bridge piers.” J. Hydraul. Eng., 114(10), 1210–1226. 42. Michiue, M., and Hinokidani, O. (1992). “Calculation of 2-dimensional bed evolution around spur-dike.” Ann. J. Hydraul. Eng., 36, 61–66 (in Japanese). 66 43. Nagata N, Hosoda T, Nakato T và Muramoto Y (2005). “Three-dimensional numerical model for flow and bed deformation around river hydraulic structures”. J. of Hydraulic Engineering, 131(12), 1074-1087. 44. Olsen, N. R. B. (2003). “Three-dimensional CFD modeling of selfforming meandering channel.” J. Hydraul. Eng., 129(5), 366–372. 45. Olsen, N. R. B., and Kjellesvig, H. M. K. (1998). “Three-dimensional numerical flow modeling for estimation of maximum local scour depth.” J. Hydraul. Res., 36(4), 579–590. 46. Olsen, N. R. B., and Melaaen, M. C. (1993). “Three-dimensional calculation of scour around cylinders.” J. Hydraul. Eng., 119(9), 1048–1054. 47. Ouillon, S., and Dartus, D. (1997). “Three-dimensional computation of flow around groyne.” J. Hydraul. Eng., 123(11), 962–970. 48. Pope, S. B. (1975). “A more general effective viscosity hypothesis.” J. Fluid Mech., 72, 331–340. 49. Rahman, M. M., Nagata, N., Muramoto, Y., and Murata, H. (1998). “Effect of side slope on flow and scouring around spur-dike-like structures.” Proc., 7th Int. Symp. on River Sedimentation, Hong Kong, China, 165–171. 50. Rajaratnam, N., and Nwachukwu, B. A. (1983a). “Flow near groin-like structures.” J. Hydraul. Eng., 109(3), 463–480. 51. Rajaratnam, N., and Nwachukwu, B. A. (1983b). “Erosion near groyne-like structures.” J. Hydraul. Res., 21(4), 277–287. 52. Raudkivi, A. J. (1986). “Functional trends of scour at bridge piers.” J. Hydraul. Eng., 112(1), 1–13. 53. Raudkivi, A. J., and Ettema, R. (1977). “Effect of sediment gradation on clear water scour.” J. Hydraul. Div., Am. Soc. Civ. Eng., 103(10), 1209–1213. 54. Raudkivi, A. J., and Ettema, R. (1985). “Scour at cylindrical bridge piers in armored beds.” J. Hydraul. Eng., 111(4), 713–731. 55. Richardson, J. E., and Panchang, V. G. (1998). “Three-dimensional simulation of scour-inducing flow at bridge piers.” J. Hydraul. Eng., 124(5), 530–540. 56. Roulund, A., Sumer, B. M., Fredsoe, J., and Michelsen, J. (1998). “3D mathematical modelling of scour around a circular pile in current.” Proc., 7th Int. Symp. on River Sedimentation, Hong Kong, China, 131–137. 57. Sinha, S. K., Sotiropoulos, F., and Odgaard, A. J. (1998). “Threedimensional numerical model for flow through natural rivers.” J. Hydraul. Eng., 124(1), 13–24.