Chương 1 Tổng quan tài liệu Trong các nghiên cứu ứng dụng về môi trường, việc giải quyết các bài toán về mô hình thủy lực luôn là một yêu cầu rất cấp thiết. Mô hình thủy lực là sự mô phỏng các quá trình, hiện tượng thủy lực - sự vận động rất phức tạp của nước trong tự nhiên dưới dạng các phương trình toán học, logic và giải chúng trên các máy tính. Mô hình thủy lực có khả năng xem xét những diễn biến của hiện tượng thủy lực từ vi mô tới vĩ mô. Đối với các bài toán mô hình hóa môi trường, để đạt độ chính xác cao nhất cần phải có những mô hình toán đáng tin cậy để đảm bảo kết quả tính toán tương ứng với kết quả đo đạc thực tế.
Trong những năm gần đây, các thuật toán để giải quyết phương trình dòng chảy cho bài toán nước nông (chiều sâu của chất lỏng là nhỏ so với chiều ngang của quá trình lan truyền) trên lưới phi cấu trúc đã đạt được thành tựu và kết quả ấn tượng. Bằng cách khai thác sự linh hoạt cao của lưới phi cấu trúc và áp dụng phương pháp thể tích hữu hạn, hệ phương trình thủy lực có thể được giải quyết để tin học hóa. Tính linh hoạt của phương pháp thể tích hữu hạn thực hiện trên lưới phi cấu trúc là thích hợp với những lời giải kĩ thuật dựa trên lưới cấu trúc hoặc phương pháp phần tử hữu hạn. Trong [6] Zhao và các cộng sự xem xét phương trình nước nông 2D có tính tới nhớt và sử dụng lược đồ Osher trên lưới phi cấu trúc.
Trong công trình [7], trong khuôn khổ sơ đồ thể tích hữu hạn trung tâm dạng Godunov, đã xử lí mô hình nước nông hai chiều trên lưới phi cấu trúc. Trong nghiên cứu này, sơ đồ Albada được sử dụng. Xấp xỉ Riemann, van Leer, Harten, Lax cũng được xem xét để tính toán và xấp xỉ thành phần thông lượng và moment. Công trình này được công bố trên tạp chí Advances in Water Resorces 12 [7] được coi là tạp chí uy tín bậc nhất trong giới chuyên môn liên quan.
Để nhận được hệ phương trình mô tả chuyển động thủy triều (Jingming Hou et al, 2013) xem xét phương trình dòng chảy và phương trình liên tục của chất lỏng không nén được. Dưới dạng vector, hệ này có dạng ∂q ∂f ∂g + + =S (1.1) ∂t ∂x ∂y Trong đó: t là thời gian (s) x và y là chiều trong hệ tọa độ đề các (m) q, f và g lần lượt là thông lượng, luồng theo phương x và y (m2 /s) S là vector nguồn, gồm 2 thành phần nguồn góc nghiêng Sb và ma sát Sf " # " # " # h qx qy q= qx ,f = uqx + gh2 /2 ,g = vqx , qy uqy uqy + gh2 /2 (1.2) √0 " # " # 0 S = Sb + Sf = −gh∂zb /∂x + −Cf u√u2 + v 2 −gh∂zb /∂y −Cf v u2 + v 2 Trong đó: qx = uh là thông lượng theo phương x qy = vh là thông lượng theo phương y u, v lần lượt là thành phần của vector vận tốc dòng chảy trung bình theo độ sâu trong hệ tọa độ vuông góc Oxy (m/s) h là độ sâu tính từ đáy đến mặt nước (m) zb là cao trình đáy (m) xem hình bên dưới Hình 1.1: Mô hình dòng chảy gn Cf = là hệ số ma sát. h1/3 13 η = h + zb là cao trình (m). Các tác giả xây dựng xấp xỉ thể tích hữu hạn trên miền ΓhN với thể tích kiểm tra Ti , i = 1, 2, .2: Chiếu miền 3D xuống mặt phẳng 2D Với ô bất kì, dạng tích phân phương trình (1.3) ∂t ∂x ∂y Ω Ω Ω Sử dụng định lí Divergence, phương trình (1.3) được viết lại như sau: Z I Z ∂q dΩ + F (q) · ndΓ = Sb + Sf dΩ ∂t Ω Γ Ω Viết lại dạng 2 chiều: Z I Z ∂q dTi + F (q) · ni dΓ = Sb + Sf dTi (1.4) ∂t Ti Γ Ti Trong đó: Ω là thể tích kiểm tra thử, trong hai chiều thì nó là diện tích m2 của tam giác Ti Γ là biên của thể tích kiểm tra, trong hai chiều nó là cạnh của tam giác Ti ni là vector pháp tuyến đơn vị ứng với cạnh Γ của tam giác Ti.
Và gồm hai thành phần (nx , ny ) Và F (q).ni là vector pháp tuyến thông lượng ứng với biên đang xét. qx nx +qy ny F (q) · ni = (f nx + gny ) = uqx + gh2 /2 nx + vqx ny (1.5) 2 uqx qy + uqy + gh /2 ny 14 → Hình 1.3: Định nghĩa vector pháp tuyến n của tam giác Ti Sau mỗi bước thời gian kế tiếp, q của ô i được tính lại theo công thức sau Z I ∆t qin+1 = qin + Sb + Sf dTi − F (q) · ni dΓ (1.6) |Ti | T Γ Và tích phân của F(q) · ni tại ô i được xấp xỉ bởi phương pháp Euler tường minh như sau: I X3 F (q) · ni dΓ = Fk (q n ) · nk lk Γ k=1 Trong đó: k=1 là chỉ số và chiều dài của cạnh tam giác Ti Hình 1.4: Những tam giác chung cạnh với tam giác Ti Ta có ma trận Jacobian của phương trình (1.1) ∂f1 ∂f1 ∂f1 ∂h ∂qx ∂qy " # ∂f ∂f ∂f2 ∂f2 0 1 0 2 = = c2 − u2 2u 0 = A ∂q ∂h ∂qx ∂qy −uv v u ∂f3 ∂f3 ∂f3 ∂h ∂qx ∂qy 15 ∂g1 ∂g1 ∂g1 ∂h ∂qx ∂qy " # ∂g ∂g2 ∂g2 ∂g2 0 0 1 = = −uv v u =B ∂q ∂h ∂qx ∂qy c 2 − u2 0 2v ∂g3 ∂g3 ∂g3 ∂h ∂qx ∂qy Trong đó: √ c = gh là tốc độ truyền sóng và ∂F(q) · ni J= = Anx + Bny ∂q 0 2 − u2 n − uvn nx ny J= c x y 2unx + vny un y 2 2 −uvnx + c − v ny vnx unx + 2vny Trị riêng của ma trận J: λ1 = unx + vny − c, λ2 = unx + vny , λ3 = unx + vny + c Vector riêng: T T r1 = [ 1 u − cnx v − cny ] , r2 = [ 0 −nx ny ] , T r3 = [ 1 u + cnx v + cny ] Các biên lưu lượng nằm tại tam giác và tập hợp các tam giác có chung cạnh với tam giác Ti , i = 1, 2,. , N lk = Ti ∩ Tk , i = 1, 2,. , N là cạnh của tam giác Ti và k = 1, 2, 3 Hình 1.5: Định nghĩa tam giác kiểm tra Từ phương trình (1.6), ta viết lại như sau: 3 Z Z ∆t ∆t X ∼ ∼ qin+1 = qin − f nx +g ny dlk + Sb + Sf dTi (1.7) |Ti | |Ti | k=1 l T k 16 Trong đó: |Ti | : diện tích tam giác (m2 ) ∼ ∼ ∼ ni = nx , ny : vector pháp tuyến đơn vị (tương ứng của Ti tại cạnh ∂Ti ∩ ∂Tk ) Đặt: Z ∼ ∼ Φk = f nx + g ny dlk (1.8) lk Ta viết lại phương trình (1.7) như sau: 3 Z ∆t X ∆t qn+1 i = qni − Φk + (Sb + Sf )dTi (1.9) |Ti | |Ti | k=1 T Sử dụng bài toán Riemann của Roe[1] để xấp xỉ trao đổi lưu lượng như sau: Φk = Z(qLi , nk ) + J˜LR (qLk − qLi ) (1.10) Trong đó: J˜LR ="(P̃ Λ̃− P̃ 1 )LR , Λ̃− = diag{λ̃− − i }, λ̃i = min(λ̃i , 0), i =# 1, 2, 3 ũnx + ṽny − c̃ knk k 0 0 Λ̃− = 0 ũnx + ṽny 0 0 0 ũnx + ṽny + c̃ knk k " # 1 0 0 P̃ = ũ − c̃nx −ny ũ + c̃nx ṽ − c̃ny nx ṽ + c̃ny ũn + ṽn 1 nx ny x y + − − 2c̃ 2 2c̃ 2c̃ (ũny − ṽnx ) ny nx − P̃ = − knk k k nk k k nk k ũny + ṽnx 1 n ny x − + knk k 2 2c̃ 2c̃ Kí hiệu của các ma trận được tính toán bằng cách sử dụng các giá trị Roe trung bình của các biến nguyên thủy như sau: √ √ √ hL uL + hR uR h̃ = hL · hR , ũ = √ √ , hL +r hR √ √ hL v L + hR v R hL + hR ṽ = √ √ , c̃ = g hL + hR 2 Và qxL nx + qyL ny ) L L 1 2 Z(qLi , nk ) = fL nx + gL ny = qx (u nx + v L ny ) + g L (hL ) nx 2 1 L L 2 L L L qy (u nx + v ny ) + g (h ) nqy 2 17 Mặt khác, qLi = qM R M i và qk = qk với M là trung điểm của cạnh ∂Ti ∩ ∂Tk Để tính các giá trị lưu lượng tại điểm M chúng ta phải tính thông qua điểm D.
* Tính qM i và qk M |rLD | qiD = qi + ψ (∇q)upw cent L · rLR , (∇ q ) · r LR |rLR | (1.11) |rRD | qkD = qk + ψ (∇q)upw cent R · rLR , (∇ q ) · rLR |rLR | Hình 1.6: Xác định các giá trị lưu lượng tại vị trí D Trong đó: [~rRL : vector từ tâm của tam giác trái (tam giác Ti ) tới tâm của tam giác phải(tam giác Tk ). Các đạo hàm được biểu diễn như sau: (∇q)cent · ~rLR = qi − qk = qR − qL (∇q)upw L = (∇q)upw i = 2(∇q)L − (∇q)cent (1.12) (∇q)upw R = (∇q)upw k = 2(∇q)R − (∇q)cent Thêm vào đó, ψ(a, b) làm hàm giới hạn với 2 biến số a và b. Giới hạn Van Albada cho phép đạt được độ chính xác bậc 2, và được viết như sau: (a2 + e)b + (b2 + e)a ( , nếu ab > 0 ψ(a, b) = 2 2 a + b + 2e 0 nếu ab 6 0 Cuối cùng ta tính toán các giá trị lưu lượng tại vị trí M. qiM = qD i + rDM ∇qL (1.13) qM M k = qk + rDM ∇qR L , η̄ R , h̄L , h̄R chiều sâu của mực nước tới đáy sông tại điểm M của Mực nước η̄M M M M 18 Hình 1.