Đề tài Ứng dụng mô hình toán phục vụ quy hoạch lưu vực sông Trà Khúc

của dòng chảy. Ph−ơng pháp số thứ hai là ph−ơng pháp phần tử hữu hạn (Wang và Anderson, 1982) [31]. H−ớng nghiên cứu cơ bản của ph−ơng pháp này hiện nay đ−ợc ứng dụng trong cơ học xây dựng, nơi mà các phần tử là các phần có thật của một cấu trúc, nó quyết định về các rầm và các cột trong khung s−ờn của một toà nhà, hoặc về l−ới của các rầm trong sàn của một cây cầu. Sự biến dạng của các yếu tố đ−ợc biểu diễn trong các số hạng của lực tác động lên hai đầu. Điều này cho phép biểu thị sự thay thế của mỗi điểm nút bằng các nút lân cận, và sự biến dạng của các phần tử liên quan. Hệ thống ph−ơng trình cuối cùng nhận đ−ợc từ điều kiện cân bằng tại mỗi nút. Trong bài toán về dòng chảy s−ờn dốc cũng có thể hình dung rằng một vùng đ−ợc phân chia thành các phần tử nhỏ với mỗi đặc tính vật lý riêng, bằng cách đó đối với mỗi một phần tử dòng chảy đ−ợc mô tả trong đặc tính của các điểm giao. Sử dụng hệ Saint - Venant vào mỗi phần tử với hệ các ph−ơng trình đại số nhận đ−ợc từ điều kiện mà dòng chảy phải liên tục tại mỗi nút. 30 Cách th−ờng dùng để mô tả ph−ơng pháp phần tử hữu hạn không dùng nh− là một lập luận mang tính vật lý. Thay vì sử dụng đối số toán học thì sử dụng hàm trọng số nào đó, trong đó hệ thống các ph−ơng trình nhận đ−ợc do yêu cầu ph−ơng trình sai phân thoả mãn "ở sát trung bình". Hệ thống các ph−ơng trình nhận đ−ợc trong ph−ơng pháp phần tử hữu hạn có cấu trúc giống nh− trong ph−ơng pháp sai phân hữu hạn. Trên thực tế, hai ph−ơng pháp rất giống nhau và đối với một bài toán nào đó thì chúng có thể đ−ợc xem xét nh− là hai quá trình biểu diễn của một mô hình toán đơn. Tuy nhiên, cách thức xuất phát và phát triển th−ờng biểu thị một sự khác nhau nào đó. Thí dụ chẳng hạn, dạng tự nhiên và đơn giản nhất của phần tử là dạng hình tam giác, làm cho sự miêu tả tr−ờng một cách linh hoạt hơn, trong khi đó các mắt l−ới tự nhiên và đơn giản nhất trong ph−ơng pháp sai phân hữu hạn là mạng vuông hoặc hình chữ nhật, nó kém linh động hơn. Thuận lợi khác của ph−ơng pháp phần tử hữu hạn là công thức chuyển của nó có tính chất trung gian mà mỗi một phần tử có thể có các giá trị riêng cho các tham số vật lý nh− là các tham số về dẫn truyền và tích trữ. Để xấp xỉ l−u vực sông bằng các phần tử hữu hạn, lòng dẫn đ−ợc chia thành các phần tử lòng dẫn và s−ờn dốc đ−ợc chia thành các dải t−ơng ứng với mỗi phần tử lòng dẫn sao cho: trong mỗi dải dòng chảy xảy ra độc lập với dải khác và có h−ớng vuông góc với dòng chảy trong phần tử lòng dẫn. Trong mỗi dải lại chia ra thành các phần tử s−ờn dốc sao cho độ dốc s−ờn dốc trong mỗi phần tử t−ơng đối đồng nhất. Việc mô phỏng l−u vực bằng các phần tử hữu hạn nh− vậy cho phép chuyển bài toán hai chiều (2D) trên s−ờn dốc thành bài toán một chiều (1D) trên s−ờn dốc và trong sông. Vì vậy, theo lý thuyết Bephanhi A. N. [26] cho phép áp dụng mô hình sóng động học một chiều cho từng dải s−ờn dốc.Mô hình phần tử hữu hạn sóng động học đánh giá tác động của việc sử dụng đất trên l−u vực đến dòng chảy đ−ợc xây dựng dựa trên hai ph−ơng pháp: ph−ơng pháp phần tử hữu hạn để mô tả quá trình lan truyền vật chất trên s−ờn dốc và trong lòng dẫn và ph−ơng pháp SCS để mô tả quá trình tổn thất trên bề mặt l−u vực [21]. 2.1. Ph−ơng pháp SCS Cơ quan bảo vệ thổ nh−ỡng Hoa Kỳ (1972) đã phát triển một ph−ơng pháp để tính tổn thất dòng chảy từ m−a rào (gọi là ph−ơng pháp SCS) [28]. Ta đã thấy, trong một trận m−a rào, độ sâu m−a hiệu dụng hay độ sâu dòng chảy trực tiếp Pe không bao giờ v−ợt quá độ sâu m−a P. T−ơng tự nh− vậy, sau khi quá trình dòng chảy bắt đầu, độ sâu n−ớc bị cầm giữ có thực trong l−u vực, Fa bao giờ cũng nhỏ hơn hoặc bằng một độ sâu n−ớc cầm giữ có thực trong l−u vực, mặt khác Fa bao giờ cũng nhỏ hơn hoặc bằng một 31 độ sâu n−ớc cầm giữ tiềm năng tối đa nào đó S (hình 1.4). Đồng thời còn có một l−ợng Ia bị tổn thất ban đầu nên không sinh dòng chảy, đó là l−ợng tổn thất ban đầu tr−ớc thời điểm sinh n−ớc đọng trên bề mặt l−u vực. Do đó, ta có l−ợng dòng chảy tiềm năng là P - Ia. Trong ph−ơng pháp SCS, ng−ời ta giả thiết rằng tỉ số giữa hai đại l−ợng có thực Pe và Fa thì bằng với tỉ số giữa hai đại l−ợng tiềm năng P - Ia và S. Vậy ta có: a ea IP P S F −= (2.1) Từ nguyên lí liên tục, ta có: aae FIPP ++= (2.2) Kết hợp (2.1) và (2.2) để giải Pe ( ) SIP IP P a a e +− −= 2 (2.3) Đó là ph−ơng trình cơ bản của ph−ơng pháp SCS để tính độ sâu m−a hiệu dụng hay dòng chảy trực tiếp từ một trận m−a rào [28]. Hình 2.1: Các biến số có tổn thất dòng chảy trong ph−ơng pháp SCS Ia - độ sâu tổn thất ban đầu, Pe - độ sâu m−a hiệu dụng, Fa - độ sâu thấm liên tục, P - tổng độ sâu m−a. Qua nghiên cứu các kết quả thực nghiệm trên nhiều l−u vực nhỏ, ng−ời ta đã xây dựng đ−ợc quan hệ kinh nghiệm : Ia = 0,2S Trên cơ sở này, ta có : ( ) SP SP Pe 8.0 2.0 2 + −= (2.4) Lập đồ thị quan hệ giữa P và Pe bằng các số liệu của nhiều l−u vực, ng−ời ta đã tìm ra đ−ợc họ các đ−ờng cong. Để tiêu chuẩn hoá các đ−ờng cong này, ng−ời ta sử dụng số hiệu của đ−ờng cong, CN làm thông số. Đó là một số không thứ nguyên, lấy 32 giá trị trong khoảng 1000 ≤≤ CN . Đối với các mặt không thấm hoặc mặt n−ớc, CN = 100 ; đối với các mặt tự nhiên, CN < 100. Số hiệu của đ−ờng cong và S liên hệ với nhau qua ph−ơng trình : 10 1000 −= CN S inch) hay ⎟⎠ ⎞⎜⎝ ⎛ −= 1010004.25 CN S (mm) (2.5) Các số hiệu của đ−ờng cong CN đã đ−ợc cơ quan bảo vệ thổ nh−ỡng Hoa Kỳ lập thành bảng tính sẵn [28] dựa trên phân loại đất và tình hình sử dụng đất. 2.2. Ph−ơng pháp phần tử hữu hạn Một ví dụ của mô hình thuỷ động lực học là mô hình của Ross B.B và nnk., Đại học Quốc gia Blacksburg, Mỹ [34]. Mô hình dùng để dự báo ảnh h−ởng của việc sử dụng đất đến quá trình lũ. M−a v−ợt thấm là đầu vào của mô hình. Ph−ơng pháp phần tử hữu hạn số kết hợp với ph−ơng pháp số d− của Galerkin đ−ợc sử dụng để giải hệ ph−ơng trình sóng động học của dòng chảy một chiều. Việc áp dụng lý thuyết phần tử hữu hạn để tính toán dòng chảy đ−ợc Zienkiewicz và Cheung (1965) [1] khởi x−ớng. Các tác giả đã sử dụng ph−ơng pháp này để phân tích vấn đề dòng chảy thấm. Nhiều nhà nghiên cứu khác cũng đã áp dụng ph−ơng pháp phần tử hữu hạn để giải quyết các vấn đề của dòng chảy Oden và Somogyi (1969), Tong (1971) [9, 13, 28, 30, 34]. Judah (1973) [9, 21] đã tiến hành việc phân tích dòng chảy mặt bằng ph−ơng pháp phần tử hữu hạn. Tác giả đã sử dụng ph−ơng pháp số d− của Galerkin trong việc xây dựng mô hình diễn toán lũ và đã thu đ−ợc kết quả thoả mãn khi mô hình đ−ợc áp dụng cho l−u vực sông tự nhiên. Tác giả cho rằng mô hình phần tử hữu hạn dạng này gặp ít khó khăn khi l−u vực có hình học phức tạp, sử dụng đất đa dạng và phân bố m−a thay đổi. Ph−ơng pháp phần tử hữu hạn kết hợp với ph−ơng pháp Galerkin còn đ−ợc Al- Mashidani và Taylor (1974) áp dụng để giải hệ ph−ơng trình dòng chảy mặt ở dạng vô h−ớng[35]. So với các ph−ơng pháp số khác, ph−ơng pháp phần tử hữu hạn đ−ợc coi là ổn định hơn, hội tụ nhanh hơn và đòi hỏi ít thời gian chạy hơn. Cooley và Moin (1976) [30] cũng áp dụng ph−ơng pháp Galerkin khi giải bằng ph−ơng pháp phần tử hữu hạn cho dòng chảy trong kênh hở và thu đ−ợc kết quả tốt. ảnh h−ởng kỹ thuật tổng hợp thời gian khác nhau cũng đ−ợc đánh giá. Ph−ơng pháp phần tử hữu hạn đặc biệt đ−ợc ứng dụng vào việc đánh giá ảnh h−ởng của những thay đổi trong sử dụng đất đến dòng chảy lũ vì l−u vực có thể đ−ợc chia thành một số hữu hạn các l−u vực con hay các phần tử. Những đặc tính thuỷ văn của một hoặc tất cả các phần tử có thể đ−ợc thay đổi để tính toán các tác động đến phản ứng thủy văn của toàn bộ hệ thống l−u vực. 33 2.2.1. Xây dựng mô hình Desai và Abel (1972) [32] đã kể ra những b−ớc cơ bản trong ph−ơng pháp phần tử hữu hạn nh− sau: 1. Rời rạc hoá khối liên tục. 2. Lựa chọn các mô hình biến số của tr−ờng. 3. Tìm các ph−ơng trình phần tử hữu hạn. 4. Tập hợp các ph−ơng trình đại số cho toàn bộ khối liên tục đã rời rạc hoá. 5. Giải cho vector của các biến của tr−ờng tại nút. 6. Tính toán các kết quả của từng phần tử từ biên độ các biến của tr−ờng tại nút. Những b−ớc này sẽ đ−ợc sử dụng trong việc phát triển mô hình dòng chảy mặt và dòng chảy trong sông. (A) (B) Hình 2.2. L−u vực và l−ới phần tử hữu hạn t−ơng ứng. 1. Rời rạc hoá khối liên tục: Khối liên tục, tức là hệ thống vật lý đang nghiên cứu đ−ợc chia thành một hệ thống t−ơng đ−ơng gồm những phần tử hữu hạn. Việc rời rạc hoá thực sự là một quá trình cân nhắc vì số l−ợng, kích th−ớc và cách sắp xếp của các phần tử hữu hạn đều có liên quan đến chúng. Dù vậy cần xác định một phần tử sao cho bảo toàn đ−ợc tính chất đồng nhất thủy văn trong mỗi phần tử. Tính chất đồng nhất thuỷ lực cũng là một mục tiêu cần xem xét tiếp theo khi tạo ra l−ới. Có thể sử dụng một số l−ợng lớn các phần tử, nh−ng số l−ợng các phần tử th−ờng hạn chế do những IA1 IA2 IB1 IA3 IIA1 IIIA1 IIIA2 IIIA3 IIIC1 IIID1 IIIB1 IIID1 IB1 IC2 IIC1 IID1 34 điều kiện ràng buộc thời gian và kinh phí. Một l−u vực giả thuyết đ−ợc sử dụng để minh hoạ cho quá trình này. L−u vực bao gồm một dòng chính và một nhánh lớn. Cả hai nhánh này đều đ−ợc đ−a vào sơ đồ dòng chảy. Ba l−u vực con hay bãi dòng chảy trên mặt đ−ợc xác định. Ngoài ra, ba kênh có thể đ−ợc xác định. Dù vậy, bất kỳ số l−ợng bãi dòng chảy bề mặt hay kênh có thể xác định nếu nh− có số liệu mặt cắt ngang của kênh. Trên hình 2.2B, những đ−ờng đậm là ranh giới gần đúng của l−u vực và các bãi dòng chảy mặt. B−ớc tiến hành tiếp theo là xác định các thành phần của kênh. Cách thức đơn giản nhất là chia mỗi một trong 3 kênh thành một số l−ợng các đoạn bằng nhau thích hợp. Từ những nút của các phần tử kênh này kẻ các các đ−ờng ra phía ngoài làm ranh giới của các l−u vực con thành một phần tử kênh. Trong tr−ờng hợp có một l−u vực thực tế thì các bản đồ địa hình của khu vực sẽ cung cấp cơ sở cho việc vạch ra các ranh giới này. Các đ−ờng này xác định các dải trong đó dòng chảy mặt diễn ra một cách độc lập với các dải khác và theo h−ớng vuông góc với dòng chảy trong các phần tử kênh. Khái niệm này cho phép có thể sử dụng việc phân tích một chiều. Các phần tử bổ sung đ−ợc hình thành bằng cách vẽ các đ−ờng song song với các phần tử kênh, bằng cách đó chia mỗi một dải thành một hệ thống các phần tử. Xét bãi dòng chảy mặt thứ nhất, quá trình giải là quá trình phân tích phần tử hữu hạn cho từng dải với m−a v−ợt thấm là đầu vào để tìm ra dòng chảy mặt chảy vào kênh dẫn. Sau đó phân tích phần tử hữu hạn cho kênh dẫn đ−ợc thực hiện t−ơng tự nh− với một dải dòng chảy mặt riêng lẻ để tìm ra l−u l−ợng trong kênh dẫn tại vị trí các nút phần tử kênh. Quá trình này đ−ợc lặp lại cho các bãi dòng chảy còn lại để tìm đ−ợc quá trình l−u l−ợng tại nút hạ l−u của toàn bộ l−u vực. Việc đánh số đúng các phần tử bãi dòng chảy sẽ chỉ ra đ−ợc chính xác từng phần tử, dải và bãi dòng chảy. Theo thí dụ trên hình 2.2B, các số La Mã biểu thị các bãi dòng chảy, các chữ in hoa biểu thị các dải và các chữ sô th−ờng biểu thị các phần tử trong dải. 2.Lựa chọn mô hình biến số của tr−ờng: B−ớc này bao gồm việc lựa chọn các mẫu giả định về các biến của tr−ờng trong từng phần tử và gán các nút cho từng phần tử. Các hàm số mô phỏng xấp xỉ sự phân bố của các biến của tr−ờng trong từng phần tử hữu hạn là các ph−ơng trình thủy động học liên tục và động l−ợng. Hệ ph−ơng trình này đã đ−ợc chứng tỏ có thể áp dụng đ−ợc cho cả dòng chảy trên mặt và dòng chảy trong kênh. 35 Ph−ơng trình liên tục: 0=−+ q t A x Q ∂ ∂ ∂ ∂ (2.6) Ph−ơng trình động l−ợng x y gASSgA A Q xt Q f ∂ ∂ ∂ ∂ ∂ ∂ −−=⎟⎟⎠ ⎞⎜⎜⎝ ⎛+ )( 2 (2.7) trong đó: Q - L−u l−ợng trên bãi dòng chảy trên mặt hoặc trong kênh; q - dòng chảy bổ sung ngang trên một đơn vị chiều dài của bãi dòng chảy (m−a v−ợt thấm đối với bãi dòng chảy trên mặt và và đầu ra của dòng chảy trên mặt đối với kênh dẫn); A- Diện tích dòng chảy trong bãi dòng chảy trên mặt hoặc trong kênh dẫn; x- khoảng cách theo h−ớng dòng chảy; t thời gian; g gia tốc trọng tr−ờng; S độ dốc đáy của bãi dòng chảy; Sf độ dốc ma sát; y độ sâu dòng chảy. Việc xấp xỉ sóng động học đ−ợc áp dụng đối với ph−ơng trình động l−ợng. Đó là sự lựa chọn tốt nhất vì các điều kiện biên và điều kiện ban đầu chỉ cần áp dụng đối với ph−ơng trình liên tục. Tính đúng đắn của quá trình này đã đ−ợc nói đến trong nhiều tài liệu (Lighthill và Witham, 1955; Woolhiser và Liggett, 1967) [1,13, 32]. Việc xấp xỉ động học đòi hỏi sự cân bằng giữa các lực trọng tr−ờng và quán tính trong ph−ơng trình động l−ợng và dòng chảy là hàm số chỉ phụ thộc vào độ sâu. Do đó ph−ơng trình động l−ợng có thể rút gọn về dạng: S = Sf (2.8) Ph−ơng trình (2.1) có thể biểu diễn d−ới dạng ph−ơng trình dòng chảy đều nh− ph−ơng trình Chezy hoặc Manning. Ph−ơng trình Manning đ−ợc chọn cho việc giải này: ASRQ 2/13/2 1 n = (2.9) trong đó: R - bán kính thuỷ lực (diện tích/chu vi −ớt); n- hệ số nhám Manning. Sau khi xấp xỉ sóng động học sẽ còn lại hai biến của tr−ờng cần xác định là A và Q. Cả hai đều là những đại l−ợng có h−ớng, do vậy có thể áp dụng sơ đồ một chiều. Khi đ−ợc biểu diễn trong dạng ẩn tại các điểm nút, A và Q có thể đ−ợc coi là phân bố trong từng phần tử theo x nh− sau: A(x,t) ≈ A* (x,t) = [ ]{ }ANtAxNn i ii∑ = = 1 )()( (2.10) 36 Q(x,t) ≈ Q*(x,t) = [ ]{ }QNtQxNn i ii∑ = = 1 )()( (2.11) trong đó: Ai(t) - diện tích, là hàm số chỉ phụ thuộc vào thời gian; Qi(t) - l−u l−ợng, hàm số chỉ phụ thuộc vào thời gian; Ni(x) - hàm số nội suy; n - số l−ợng nút trong một phần tử. Đối với một phần tử đ−ờng một chiều, n = 2 và: A* (x,t) = Ni(x) Ai(t) + Ni+1(x)Ai+1(t) (2.12) Q* (x,t) = Ni(x)Qi(t) + Ni+1(x)Qi+1(t) (2.13) trong đó: i i i x xx xN Δ −= +1)( và i i i x xx xN Δ −=+ )(1 với x ∈ (xi , xi+1) Các hàm nội suy th−ờng đ−ợc coi là các hàm toạ độ vì chúng xác định mối quan hệ giữa các toạ độ tổng thể và địa ph−ơng hay tự nhiên. Các hàm nội suy đối với các phần tử đ−ờng đã đ−ợc bàn luận t−ơng đối kỹ trong nhiều bài viết về phần tử hữu hạn (Desai và Abel, 1972; Huebner, 1975)[1, 21] . 3.Tìm hệ ph−ơng trình phần tử hữu hạn: Việc tìm các ph−ơng trình phần tử hữu hạn bao gồm việc xây dựng hệ ph−ơng trình đại số từ tập hợp các ph−ơng trình vi phân cơ bản. Có bốn quy trình th−ờng đ−ợc sử dụng nhất là ph−ơng pháp trực tiếp, ph−ơng pháp cân bằng năng l−ợng, ph−ơng pháp biến thiên và ph−ơng pháp số d− có trọng số. Ph−ơng pháp số d− có trọng số của Galerkin đ−ợc dùng để thiết lập các ph−ơng trình vì nó đã chứng tỏ là một ph−ơng pháp tốt đối với các bài toán về dòng chảy mặt (Judah, 1973; Taylor và nnk, 1974)[35]. Ph−ơng pháp Galerkin cho rằng tích phân: ∫ = D iRdDN 0 (2.14) D: khối chứa các phần tử. R: số d− đ−ợc gán trọng số trong hàm nội suy Ni Do ph−ơng trình (2.14) đ−ợc viết cho toàn bộ không gian nghiệm nên nó có thể đ−ợc áp dụng cho từng phần tử nh− d−ới đây, ở đó hàm thử nghiệm sẽ đ−ợc thay thế vào ph−ơng trình (2.14) và lấy tích phân theo từng phần tử của không gian : 0 1 = ⎭⎬ ⎫ ⎩⎨ ⎧ ⎥⎦ ⎤⎢⎣ ⎡ −+∑ ∫ = NE i D eie dDqA x Q N &∂ ∂ (2.15) trong đó: NE : số phần tử trong phạm vi tính toán. A& : đạo hàm theo thời gian của A. De : phạm vi của một phần tử. Xét riêng một phần tử, ph−ơng trình (2.15) trở thành: 37 { } { }N Nx Q N N A N q dDi j i j iD ee ∂∂ + −⎡⎣⎢ ⎤ ⎦⎥ =∫ & 0 (2.16) Đối với 1 phần tử là đoạn thẳng, ph−ơng trình này có thể viết nh− sau { } { }N Nx Q N N A N q dxi j i j i ix x ∂ ∂ + − ⎡ ⎣⎢ ⎤ ⎦⎥ =∫ &1 2 0 (2.17) Lấy tích phân của từng số hạng trong (2.17): { } { }N N x dx Q N N x N N x N N x N N x dx Qi j x x x x∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎛ ⎝⎜ ⎞ ⎠⎟ = ⎡ ⎣ ⎢⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥⎥ ∫ ∫ 1 2 1 2 1 1 1 2 2 1 2 2 N N x dx x x x x x x x x x dx x x x x dx x x x x x x x x 1 1 1 2 2 1 2 2 1 1 2 1 2 2 1 2 1 2 ∂ ∂ ∂ ∂∫ ∫ ∫= − − − − ⎛ ⎝⎜ ⎞ ⎠⎟ = − − − = −( ) T−ơng tự, lấy tích phân của tất cả các số hạng khác, cuối cùng nhận đ−ợc: { } { }N N x dx Q Qi j x x ∂ ∂ ⎛ ⎝⎜ ⎞ ⎠⎟ = − − ⎡ ⎣ ⎢⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥⎥ ∫ 1 2 1 2 1 2 1 2 1 2 =[FQ]{Q} ( ) { } { }N N dx A x Ai j x x 1 2 1 3 1 6 1 6 1 3 ∫ = ⎡ ⎣ ⎢⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥⎥ & &Δ = [FA]{A*} N dxq xqi x x 1 2 1 2 1 2 ∫ = ⎧ ⎨ ⎪⎪ ⎩ ⎪⎪ ⎫ ⎬ ⎪⎪ ⎭ ⎪⎪ Δ = q {Fq} Kết hợp ba số hạng cho ph−ơng trình đối với một phần tử hữu hạn tuyến tính: [FA] { &A }+[FQ]{Q} - q{Fq} = 0 (2.18) Nếu đạo hàm của diện tích theo thời gian đ−ợc lấy xấp xỉ ở dạng: &A (t) = [A(t+Δt) - A(t)]/Δt ph−ơng trình (2.18) trở thành: 1 Δt [FA] {A}t+Δt - 1 Δt [FA] {A}t +[FQ]{Q}t - q{Fqt}t+Δt = 0 (2.19) 38 Hệ ph−ơng trình thiết lập cho l−ới phần tử hữu hạn gồm n phần tử đ−ợc thiết lập sao cho có thể bao hàm đ−ợc toàn bộ số phần tử. ở đây, do các dải đ−ợc diễn toán một cách độc lập nên ph−ơng trình tổng hợp cần phải viết cho từng dải và từng kênh dẫn. 4. Giải hệ ph−ơng trình cho véc tơ các biến của tr−ờng tại các nút. Hệ ph−ơng trình phần tử hữu hạn (2.19) với các ẩn số là các biến tại các nút có thể đ−ợc giải bằng ph−ơng pháp khử Gauss. Hệ ph−ơng trình phi tuyến cần phải giải thông qua các b−ớc lặp. Các điều kiện ban đầu có thể làm hệ ph−ơng trình trở nên đơn giản hơn. Ví dụ đối với một dải chứa n phần tử tuyến tính và n+1 nút, trên các bãi dòng chảy s−ờn dốc của kênh tại thời điểm t = 0, có một vài số hạng sẽ bằng 0. Ph−ơng trình phần tử hữu hạn trở thành: 1 Δt [FA] {A}t+Δt = {fq} (2.20) Sau khi giải đồng thời hệ ph−ơng trình này tìm các ẩn {A}, ph−ơng trình Manning đ−ợc sử dụng để tìm các ẩn {Q}. Điều kiện biên tiếp theo có thể làm đơn giản hoá việc giải hệ ph−ơng trình là l−u l−ợng bằng 0 ở mọi thời điểm tại các biên trên hoặc tại các nút của các dải và kênh dẫn. Có một ngoại lệ là tr−ờng hợp t−ơng tự nh− đối với 3 bãi dòng chảy s−ờn dốc và 3 kênh dẫn khi l−u l−ợng ở mọi thời điểm t tại nút trên cùng của kênh thứ 3 là tổng của các l−u l−ợng tại các nút d−ới của 2 kênh khác. Các giá trị A và Q tìm đ−ợc tại một b−ớc thời gian sẽ đ−ợc đ−a vào ph−ơng trình phần tử hữu hạn để tìm các giá trị A, Q ở b−ớc thời gian tiếp theo. Các giá trị {A}t+Δt, {Q}t+Δt tại một b−ớc thời gian tính toán sẽ trở thành các giá trị {A}t và {Q}t trong b−ớc thời gian tính toán tiếp theo. Quá trình này đ−ợc thực hiện cho đến khi tìm đ−ợc kết quả cần thiết. 5. Tổng hợp hệ ph−ơng trình đại số cho toàn bộ miền tính toán: Hệ ph−ơng trình thiết lập cho l−ới phần tử hữu hạn gồm n phần tử đ−ợc thiết lập sao cho có thể bao hàm đ−ợc toàn bộ số phần tử. ở đây, do các dải đ−ợc diễn toán một cách độc lập nên ph−ơng trình tổng hợp cần phải viết cho từng dải và từng kênh dẫn. Quá trình tổng hợp hệ ph−ơng trình cho n phần tử tuyến tính với (n+1) nút đ−ợc thực hiện nh− sau: Viết ph−ơng trình (2.19) cho một phần tử: 0 2 1 2 1 2 1 2 1 2 1 2 1 3 1 6 1 6 1 3 1 Δ 3 1 6 1 6 1 3 1 Δ 112 1 2 11 Δ2 11 = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ − ⎭⎬ ⎫ ⎩⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − − + ⎭⎬ ⎫ ⎩⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − ⎭⎬ ⎫ ⎩⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ + ql Q Q A A t l A A t l ttt Triển khai ph−ơng trình cho 3 phần tử, 4 nút: 39 0 22 1 2 1 6 1 3 1 Δ6 1 3 1 Δ 11 2121 1 Δ 21 1 =−⎟⎠ ⎞⎜⎝ ⎛ +−+⎟⎠ ⎞⎜⎝ ⎛ +−⎟⎠ ⎞⎜⎝ ⎛ + + ql QQAA t l AA t l tttt 0 22 1 2 1 3 1 6 1 Δ3 1 6 1 Δ 11 2121 1 Δ 21 1 =−⎟⎠ ⎞⎜⎝ ⎛ +−+⎟⎠ ⎞⎜⎝ ⎛ +−⎟⎠ ⎞⎜⎝ ⎛ + + ql QQAA t l AA t l tttt 0 22 1 2 1 6 1 3 1 Δ6 1 3 1 Δ 22 3232 1 Δ 32 2 =−⎟⎠ ⎞⎜⎝ ⎛ +−+⎟⎠ ⎞⎜⎝ ⎛ +−⎟⎠ ⎞⎜⎝ ⎛ + + ql QQAA t l AA t l tttt 0 22 1 2 1 3 1 6 1 Δ3 1 6 1 Δ 22 3232 1 Δ 32 2 =−⎟⎠ ⎞⎜⎝ ⎛ +−+⎟⎠ ⎞⎜⎝ ⎛ +−⎟⎠ ⎞⎜⎝ ⎛ + + ql QQAA t l AA t l tttt 0 22 1 2 1 6 1 3 1 Δ6 1 3 1 Δ 33 4343 1 Δ 43 3 =−⎟⎠ ⎞⎜⎝ ⎛ +−+⎟⎠ ⎞⎜⎝ ⎛ +−⎟⎠ ⎞⎜⎝ ⎛ + + ql QQAA t l AA t l tttt 0 22 1 2 1 3 1 6 1 Δ3 1 6 1 Δ 33 4343 1 Δ 43 3 =−⎟⎠ ⎞⎜⎝ ⎛ +−+⎟⎠ ⎞⎜⎝ ⎛ +−⎟⎠ ⎞⎜⎝ ⎛ + + ql QQAA t l AA t l tttt Hệ ph−ơng trình trên đ−ợc viết d−ới dạng ma trận nh− sau: [ ] [ ] [ ] { } 0 Δ 1 Δ 1 4 3 2 1 4 3 2 1 Δ4 3 2 1 =− ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ + ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ + qQA tt A F Q Q Q Q F A A A A F t A A A A F t trong đó: [ ] ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ + = 36 00 636 0 0 6336 00 63 33 322 2211 11 ll lll llll ll FA ; [ ] ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − − = 2 1 2 100 2 10 2 10 0 2 10 2 1 00 2 1 2 1 QF ; { } ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ + + = 2 2 2 2 33 3322 2211 11 ql qlql qlql ql fq Tổng quát cho n phần tử tuyến tính ta có ph−ơng trình dạng: 1 Δt [FA] {A}t+Δt - 1 Δt [FA] {A}t +[FQ]{Q} - q{Fq} = 0 (2.21) trong đó: [ ] . 36 00..0000 6336 0..0000 6336 00..000 .......... 0.. 6336 0000 0..0 6336 000 0..00 6336 00 00..00 6336 0 00..000 6336 00...000 63 11 1122 6655 5544 4433 3322 22 11 ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ + + + + + + + = −− −−−− nn nnnn nnnn ll A ll llll llll llll llll llll llll llll ll F 40 ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − − − − − − − − = 2 1 2 1000...00 2 10 2 100...00 0 2 10 2 100..00 .......... .......... 0..0 2 10 2 1000 0.000 2 10 2 100 00..00 2 10 2 10 00...00 2 10 2 1 00...000 2 1 2 1 QF { } ⎪⎪ ⎪⎪ ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎪⎪ ⎪⎪ ⎨ ⎧ + + + + = −− 2 22 . . 22 22 22 2 11 4433 3322 2211 11 nn nnnn q ql qlql qlql qlql qlql ql f trong đó các chỉ số của A và Q là số thứ tự của nút, các chỉ số của l và q là các chỉ số của phần tử. Một cách tiệm cận khác để giải quyết bài toán khi số liệu địa hình lòng dẫn thiếu. Khi đó cần thiết tiến hành một số thủ thuật để thay biến A bằng Q. Từ (2.9), ph−ơng trình Manning có thể viết lại là: 2/13/5 3/2 2/1 3/2 1.1 SA nP SA P A n Q =⎟⎠ ⎞⎜⎝ ⎛= (2.22) 6.0 6.0 2/1 3/2 2/1 3/2 3/5 Q S nP AQ S nP A ⎟⎟⎠ ⎞⎜⎜⎝ ⎛=⇒= Viết d−ới dạng tổng quát β β Q S nP A ⎟⎟⎠ ⎞⎜⎜⎝ ⎛= 2/1 3/2 41 Đặt β β 2/1 3/2 QA S nP αα =⇒=⎟⎟⎠ ⎞⎜⎜⎝ ⎛ Trong ph−ơng trình Manning β = 0.6. Khi đó (2.6) có thể viết lại là: q t Q Q x Q =∂ ∂+∂ ∂ −1βαβ (2.23) Đặt μ1 =−βαβQ , ph−ơng trình (2.23) trở thành: q t Q x Q =∂ ∂+∂ ∂ μ (2.24) Dễ dàng nhận thấy, áp dụng ph−ơng pháp phần tử hứu hạn (2,10), (2.11) vào (2.24) nhận đ−ợc: [ ]{ } [ ]{ } { }qAQ fQFQF =+ &μ [ ]{ } { } [ ]{ }QFfQF QqA −=&μ (2.25) Đặt )( 2 1 tttt QQQ += +Δ , ph−ơng trình (2.25) trở thành: [ ]{ } { } [ ]{ } [ ]{ }QFQFftQF QtAttqttA −+= ++ μμ ΔΔ Δ { } [ ]{ } [ ]{ } [ ]{ }ttQttQtAttq QFtQFtQFft ΔμΔ ΔΔ 5.0Δ5.0 −−+= ++ [ ] [ ][ ]{ } { } { } [ ][ ]{ }tQAttqttQA QFtFftQFtF ΔμΔΔ 5.0ΔΔ5.0μ −+=+⇒ ++ [ ]{ } [ ]{ } [ ] ttqttt FQCQB ΔΔ ++ +=⇒ Biểu thức cuối cùng sẽ là: { } [ ] [ ]{ } [ ] { } ttqttt FBQCBQ ΔΔ + −− + += 11 (2.26) Ph−ơng trình (2.26 ) có thể giải đ−ợc chỉ phụ thuộc vào l−u l−ợng. 2.3. Ch−ơng trình diễn toán lũ Trong ch−ơng trình đ−a vào các đặc tr−ng thuỷ văn nh− độ dốc, hệ số Manning, m−a v−ợt thấm trong từng phần tử. Các công trình chậm lũ hoặc hồ chứa cũng có thể đ−ợc mô hình hoá. Đầu vào của quá trình diễn toán lũ là l−ợng m−a v−ợt thấm đ−ợc tính theo ph−ơng pháp SCS. Hệ số Manning của từng phần tử cũng đ−ợc xác định theo cách lấy trung bình có trọng số. Độ dốc của từng phần tử có thể xác định theo bản đồ địa hình của khu vực. Độ dốc của các phần tử lòng dẫn có thể tìm đ−ợc theo cách t−ơng tự [21]. 42 2.4. Kiểm tra mô hình Số liệu đo đạc dòng chảy từ các bãi dòng chảy s−ờn dốc của Crawford và Linsley (1966)[30] đã đ−ợc sử dụng để kiểm tra tính đúng đắn của ch−ơng trình diễn toán lũ đối với dòng chảy s−ờn dốc. Ph−ơng pháp xấp xỉ bằng phần tử hữu hạn cho kết quả có thể thoả mãn mặc dù việc lấy hệ số Manning biến đổi theo độ sâu có thể còn cho kết quả tốt hơn nữa. Mô hình này còn có thể áp dụng cho cả l−u vực lớn trong tự nhiên (Ross, 1975). Các phép kiểm tra sự hội tụ, tính ổn định và ảnh h−ởng của của việc phân bố các l−ới ô khác nhau đến dòng chảy lũ cũng đ−ợc xét đến (Ross, 1975)[34]. 2.5. Nhận xét về khả năng sử dụng mô hình Với giả thiết của mô hình phần tử hữu hạn sóng động học có thể chia l−u vực ra thành các phần tử rất chi tiết, khi đó có thể tính toán mô phỏng dòng chảy sinh ra từ m−a ứng với từng phần tử của l−u vực, thông qua việc áp dụng mô hình sóng động học một chiều. M−a hiệu quả trên l−u vực đ−ợc tính thông qua ph−ơng pháp SCS, ph−ơng pháp này có tính đến cả tổn thất ban đầu c−ờn