Chương này kế thừa các giả định mô hình đã nêu ở Chương 2 và mô tả các phương trình chi phối, các cách tiếp cận số tiêu chuẩn, cùng những giả định đơn giản hóa gắn với các cách tiếp cận đó. Chương này không nhằm trình bày chi tiết các phương pháp hay lý thuyết số.
3.1. Nguyên lý vật lý
Việc hiểu các nguyên lý vật lý được thể hiện trong một mô hình điển hình là rất quan trọng. Đa số mô hình bao gồm cả trạng thái chất lỏng đứng yên và chuyển động.
- Các lực gắn với chất lỏng đứng yên gồm lực nổi và tác dụng của trọng lượng riêng/mật độ tạo áp lực lên đập, bờ đắp, cống và các công trình khác.
- Trái lại, các lực gắn với chất lỏng chuyển động gồm cản do ma sát và nhiễu loạn (chảy rối) do dòng chảy sinh ra.
Các đại lượng và tham số vật lý chủ chốt dùng để mô tả và định lượng chuyển động của chất lỏng gồm: gia tốc, vận tốc, khối lượng, động lượng và năng lượng:
- Gia tốc là tốc độ thay đổi của vận tốc theo thời gian.
- Vận tốc là tốc độ của một vật theo một hướng cho trước. Vận tốc là một đại lượng véc-tơ; cả hướng và độ lớn đều quan trọng.
- Theo định luật II Newton, một lực tổng hợp tác dụng lên một khối lượng sẽ gây gia tốc cho khối lượng đó theo hướng của lực, với độ lớn a=F/m. Khi lực tổng hợp không đổi (và khối lượng không đổi), gia tốc có độ lớn không đổi.
- Khối lượng đo mức “kháng lại gia tốc” của một vật.
- Trọng lượng là tích của khối lượng vật và giá trị địa phương của gia tốc trọng trường.
- Động lượng là tích của khối lượng vật và vận tốc của nó.
- Năng lượng cơ học là tổng của thế năng và động năng.
Các đại lượng vật lý này liên hệ với nhau như sau:
$$F = m\,a \tag{3.1}$$
$$M = m\,v \tag{3.2}$$
$$E \;=\; \frac{m\,v^{2}}{2} \;+\; m\,g\,d \tag{3.3}$$
trong đó:
- F = lực
- m = khối lượng
- a = gia tốc
- M = động lượng
- v = độ lớn vận tốc
- g = gia tốc trọng trường (hằng số)
- d = khoảng cách/độ cao (chiều sâu) so với mặt chuẩn
- E = năng lượng cơ học
Dựa trên các mối liên hệ này và các nguyên lý cơ bản về bảo toàn khối lượng và bảo toàn động lượng, các phương trình Navier–Stokes trung bình Reynolds (RANS) (trình bày ở mục sau) được suy ra từ các phương trình chuyển động của Newton. Hệ RANS là cơ sở cho phân tích thủy động lực học số. Các phương trình này có các ẩn số là mực nước và vận tốc. Vận tốc là đại lượng véc-tơ nên trong bài toán hai chiều sẽ có hai ẩn số vận tốc; vì vậy, biểu diễn số 2D của một miền dòng chảy có ba ẩn số tại mỗi vị trí tính toán. Các điểm tính toán có thể là nút lưới, tâm phần tử, hoặc vị trí đặc biệt bên trong mỗi phần tử tùy theo dạng công thức. Việc cài đặt số tập hợp ba ẩn số tại mỗi vị trí thành một hệ phương trình và giải đồng thời.
Khi áp dụng các phương pháp số này, cần cung cấp thêm các hằng (hoặc điều kiện) biểu diễn trường dòng để xác định trường dòng; trong lý thuyết, chúng thường được gọi là các hằng số tích phân hoặc trong số học là điều kiện biên. Trong thủy lực, người mô hình chỉ định giá trị dòng chảy (mực nước và vận tốc) tại mọi biên của miền mô hình. Việc xác định biên như vậy cho phép bộ giải số tính các giá trị tại tất cả các vị trí tính toán.
Để đơn giản phần phải khai báo, đa số mô hình giả định không có dòng qua các biên (boundary flow = 0) trừ khi được chỉ định khác đi. Do đó người dùng thường chỉ định điều kiện biên tại các vị trí dòng chảy đi vào hoặc đi ra khỏi mô hình. Với hệ thống sông suối, đó thường là đầu thượng lưu và đầu hạ lưu của mô phỏng. Các điều kiện biên này gồm mực nước và lưu lượng. Tùy cách cài đặt số, lưu lượng được chỉ định sẽ được chuyển đổi thành vận tốc tại các phần tử biên, dựa trên số liệu đầu vào của người dùng hoặc các giả định về hướng và phân bố dòng.
(nd: boundary flow)
“Không có dòng qua các biên” nghĩa là trên các biên kín của miền tính, ta đặt điều kiện không thông lượng:
\(\mathbf{u}\cdot\mathbf{n} = 0\) (không có lưu lượng pháp tuyến xuyên qua biên).
Thường áp cho bờ sông, đê, rìa miền lưới… — tức các biên bên nếu bạn không muốn nước ra/vào ở đó.
Dòng bên (lateral inflow) theo nghĩa thủy văn (nhập lưu từ suối nhánh, xả điểm, chảy tràn từ floodplain, mưa–thấm phân bố…) không tự động đi qua các biên kín. Muốn có “dòng bên”, bạn phải:
mở biên hở riêng tại vị trí suối nhánh và gán thủy đồ/mực nước, hoặc đưa vào như nguồn khối lượng nội miền (mưa, xả điểm theo đường/điểm, v.v.), hoặc mở kênh/ô lưới nối và mô phỏng công trình (weir, culvert…) để trao đổi lưu lượng.
Cơ chế của hệ số có thể hình dung như bài toán cân bằng khối lượng toàn cục. Nếu đã biết lưu lượng vào/ra của hệ, có thể tính thể tích nước tổng trong hệ tại bất kỳ thời điểm nào bằng cách cộng tổng tích lũy inflow và trừ tổng tích lũy outflow. Mô hình thực hiện cân bằng khối lượng này cho từng phần tử của lưới. Tùy phương pháp số được dùng (xem Mục 3.2), bảo toàn khối lượng có thể được cưỡng bức theo từng phần tử (local), xấp xỉ theo từng phần tử, hoặc cưỡng bức theo toàn miền (global). Một số mô hình cũng cho phép khối lượng đi vào hoặc đi ra qua nền đáy, qua mặt nước, hoặc qua các biên nội bộ. Khối lượng qua nền đáy biểu diễn thấm vào (infiltration) hoặc thấm ra (exfiltration) với tầng chứa nước ngầm. Khối lượng qua mặt nước biểu diễn mưa và bốc hơi.
Người mô hình cũng chỉ định các tham số đại diện cho các lực tác dụng lên khối nước trong từng phần tử. Các lực chính gồm ma sát đáy, ma sát bề mặt nước, và áp suất. Ma sát đáy được lượng hóa thông qua giá trị độ nhám quy định (thường là Manning n). Ma sát bề mặt được lượng hóa bằng vận tốc gió quy định. Áp suất được suy ra lặp từ mực nước tính được. Các lực này được biểu diễn dưới dạng các hạng tử vi phân riêng trong phương trình bảo toàn động lượng.
3.2. Phương trình chi phối cho mô hình 2D trung bình theo chiều sâu
Các phương trình Navier–Stokes, mang tên Claude-Louis Navier và George Gabriel Stokes, sử dụng các giả định đơn giản hóa để mô tả chuyển động của chất lỏng. Phần lớn mô hình 2D dùng một biến thể của hệ phương trình nước nông (SWE) trung bình theo chiều sâu được suy ra từ phương trình Navier–Stokes. Mỗi mô hình số có cách suy dẫn và biểu diễn riêng của các phương trình này để làm cơ sở cho việc phân tích số. Nguyên lý vận tốc trung bình theo chiều sâu đã nêu trước đó được minh họa ở Hình 3.1. Trong cột nước, mô hình trung bình theo chiều sâu thường giả định phân bố vận tốc theo phương đứng có dạng logarit.
$$U \;=\; \frac{1}{H}\int_{Z}^{Z+H} u\,dz \tag{3.4}$$
trong đó:
- U, V = vận tốc trung bình theo chiều sâu (U theo phương x, V theo phương y)
- H = độ sâu nước
- Z = cao trình đáy kênh
- u = vận tốc nước tại độ sâu z tương ứng
- z = độ sâu nước tương ứng với vận tốc uu (biến tích phân theo phương đứng)
Depth-averaged velocity principle.
Nội dung hình 3.1:
Vận tốc trong một cột nước đang chảy được biểu diễn bằng 0 tại điểm tiếp giáp với đáy và tăng theo dạng parabol đến giá trị cực đại ở mặt nước. Nếu kẻ một đường thẳng đứng từ nơi vận tốc bằng 0 ở đáy lên đến mặt nước, đồ thị trường vận tốc tạo thành nửa dương của một hình có dạng parabol.
Trên cùng đường thẳng đứng đó, vận tốc trung bình theo chiều sâu được biểu diễn bằng một giá trị vận tốc không đổi từ đáy đến mặt nước sao cho diện tích hình chữ nhật bằng với diện tích hình parabol vừa nêu.
Cách tiếp cận 2D có thể được truyền đạt trực quan bằng cách biểu diễn các thuộc tính của một cột nước hữu hạn, như minh họa ở Hình 3.2.
Nội dung hình 3.2:
Sơ đồ gồm hai phần.
+ Hình bên trái cho thấy một kênh, trong đó bề mặt được chia (trong ví dụ này) thành các phần tử hình chữ nhật. Chọn một cột hình chữ nhật (một thể tích) ngay dưới một phần tử.
+ Hình bên phải lý tưởng hóa cột nước đã chọn, với các mặt định hướng theo mặt phẳng ngang x–y. Cột được gán kích thước theo các phương x, y, z: Δx, Δy và H (chiều cao theo z). Tâm hình học (centroid) của cột tại tọa độ x, y.
Phương trình nước nông (SWE) xét đến nguyên tắc liên tục và động lượng (Hình 3.3).
Nguyên tắc liên tục (bảo toàn khối lượng) tính đến dòng chảy vào và ra khỏi tất cả sáu mặt của cột nước hữu hạn, và yêu cầu rằng mức biến thiên của khối lượng bên trong cột nước phải bằng thông lượng ròng của khối lượng đi vào cột.
Tổng lưu lượng đi vào cột nước bằng tổng lưu lượng đi ra cộng với mọi thay đổi trữ lượng
Continuity and momentum principles
Nội dung hình 3.3: Nguyên lý Liên tục và Động lượng
Sơ đồ gồm hai phần. Cả hai phần đều biểu diễn thể tích hữu hạn lý tưởng hóa như một khối lập phương, với các mặt định hướng theo mặt phẳng ngang x–y.
+ Khối bên trái minh họa liên tục tổng quát với đại lượng mật độ nước × lưu lượng đi qua khối.
+ Khối bên phải minh họa động lượng tổng quát đi qua khối.
Trong cả hai trường hợp, các đại lượng được mô tả chảy vuông góc qua các mặt đứng của khối lập phương thể tích hữu hạn lý tưởng hóa.
Nguyên tắc động lượng yêu cầu rằng đạo hàm theo thời gian của động lượng trong cột nước hữu hạn, trừ thông lượng động lượng đi vào, cộng thông lượng động lượng đi ra, bằng tổng các lực tác dụng lên cột. Các hạng tử quán tính có thể được tính đến. Sự thay đổi ròng của động lượng thường được viết dưới dạng \(\rho\,\dfrac{du}{dt}\) và (\(\rho\,\mathbf{v})\,\mathbf{v}\).
Các lực tác dụng lên cột nước gồm lực thể tích (trọng lực và các lực liên quan đến sự quay của Trái Đất – lực Coriolis) và lực bề mặt, bao gồm thủy tĩnh-hydrostatic, ma sát trượt đáy-bed shear, ứng suất gió-wind shear và chảy rối-turbulence như ở Hình 3.4. Việc biểu diễn chảy rối có thể bao gồm các hạng tử phân tán để xét sự hòa trộn của nước. Các phần khác nhau của chất lỏng chảy với tốc độ khác nhau, đưa các khối nước từ các vùng khác nhau và có động lượng khác nhau lại gần nhau, phản ánh chuyển động không đồng nhất của các hạt trong cột nước.
Việc bảo toàn khối lượng và động lượng trên toàn miền cuối cùng tạo ra một hệ phương trình đạo hàm riêng (PDE) của SWE cần được giải trong miền để tính thủy lực 2D. SWE xét đến nhiều sắc thái của dòng chảy như minh họa ở Hình 3.5.
Nội dung hình 3.4:
Sơ đồ gồm bốn phần.
Mỗi hình đều biểu diễn thể tích hữu hạn lý tưởng hóa như một khối lập phương, các mặt định hướng theo mặt phẳng ngang x–y.
+ Lực thủy tĩnh là lực áp do trọng lượng nước, tăng tuyến tính theo chiều sâu. Áp suất bằng khối lượng riêng của nước × gia tốc trọng trường × chiều cao cột nước (\(p=\rho g h\)).
+ Lực ma sát trượt đáy (bed shear) xuất hiện tại mặt tiếp xúc giữa cột nước và bề mặt nền; chúng ngược hướng chuyển động.
+ Ứng suất gió (wind shear) xuất hiện tại mặt tiếp xúc giữa cột nước và không khí; chúng cùng hướng với chuyển động của không khí và độ lớn/hướng có thể phân giải thành các thành phần vectơ theo x và y.
+ Lực do chảy rối (turbulence) tác dụng bên trong cột nước. Trong mô hình hai chiều, chảy rối chỉ được xét theo các phương x và y.
Example 2D continuity and momentum equations, and their contributing terms
Như đã nêu, các mô hình khác nhau dùng các biến thể khác nhau của SWE. Hệ dưới đây là một ví dụ điển hình.
Phương trình liên tục (The continuity equation):
$$\frac{\partial h}{\partial t}+\frac{\partial (HU)}{\partial x}+\frac{\partial (HV)}{\partial y}=e \tag{3.5}$$
Phương trình động lượng (The momentum equation) theo các phương x và y:
$$\frac{\partial (HU)}{\partial t} +\frac{\partial (HU^{2})}{\partial x} +\frac{\partial (HUV)}{\partial y} = \frac{\partial (H T_{xx})}{\partial x} +\frac{\partial (H T_{xy})}{\partial y} -\,gH\frac{\partial z}{\partial x} -\frac{\tau_{bx}-\tau_{wx}}{\rho} + D_{xx}+D_{xy} -\Omega HV \tag{3.6}$$
$$\frac{\partial (HV)}{\partial t} +\frac{\partial (HV^{2})}{\partial y} +\frac{\partial (HVU)}{\partial x} = \frac{\partial (H T_{yy})}{\partial y} +\frac{\partial (H T_{xy})}{\partial x} -\,gH\frac{\partial z}{\partial y} -\frac{\tau_{by}-\tau_{wy}}{\rho} + D_{yy}+D_{yx} +\Omega HU \tag{3.7}$$
trong đó:
- x,y: tọa độ Đề-các.
- t: thời gian.
- H: độ sâu nước.
- U,V: vận tốc trung bình theo chiều sâu theo các phương x,y.
- e: net mass input (rainfall + exfiltration – infiltration).
- z: cao trình mặt nước (WSE) == cao trình đáy + h.
- g: gia tốc trọng trường.
- \(T_{xx},T_{xy},T_{yy}\): các ứng suất cắt rối trung bình theo chiều sâu.
- \(D_{xx},D_{xy},D_{yy},D_{yx}\): hạng tử phân tán/khuyếch tán (trộn).
- \(\tau_{bx},\tau_{by}\): ứng suất cắt đáy.
- \(\tau_{wx},\tau_{wy}\): ứng suất gió.
- \(\rho\): khối lượng riêng của nước.
- \(\Omega\): tham số Coriolis.
(nd: e – net mass input)
e không phải là một lưu lượng tổng (m³/s), mà là cường độ nguồn phân bố theo mặt phẳng.
Trong phương trình liên tục
$$\frac{\partial h}{\partial t}+\frac{\partial(HU)}{\partial x}+\frac{\partial(HV)}{\partial y}=e$$
- \(\partial h/\partial t\) có đơn vị m/s,
- \(\partial(HU)/\partial x\) cũng m/s,
⇒ để cân bằng đơn vị, e phải là m/s — tức thể tích trên giây trên đơn vị diện tích (volumetric source per unit area).
Hiểu thực dụng:
- Với một ô lưới diện tích A, lưu lượng do nguồn trong ô là Q=e×A (m³/s).
- Mưa, bốc hơi, thấm… đều tự nhiên có dạng m/s (ví dụ mm/h)
Ví dụ: mưa 15 mm/h ⇒ e=0.015/3600=\(0{.}015/3600=4{.}17\times10^{-6} m/s\).
Ô 200×200 m ⇒ A=40,000 m² ⇒ \(Q = eA \approx 0{.}167 m³/s\).
Có thể ghi:
e – cường độ nguồn ròng (m/s), và ghi chú: nhân với diện tích ô để ra m³/s.
(nd: + exfiltration – infiltration)
Từ cột nước → xuống nền đáy: infiltration → thấm vào (ra khỏi cột nước) ⇒ góp dấu âm vào e.
Từ nền đáy → vào cột nước: exfiltration → thấm ra (vào cột nước) ⇒ góp dấu dương vào e.
Tóm tắt:
| Hướng trao đổi | Tiếng Anh | Thuật ngữ VN | Dấu trong e |
|---|---|---|---|
| Cột nước → nền | Infiltration | Thấm vào | − |
| Nền → cột nước | Exfiltration | Thấm ra | + |
Ghi chú: đôi khi để đỡ nhầm, người ta gọi exfiltration là “nước ngầm cấp lên kênh” hoặc “thấm lên”.
3.3. Cách tiếp cận mô hình
Nhiều cách tiếp cận đã được phát triển để xây dựng một phương pháp số có hệ thống giải các PDE. Nhìn chung, các cách này thuộc ba nhóm: sai phân hữu hạn (finite difference), phần tử hữu hạn (finite element) và thể tích hữu hạn (finite volume). Các cách tiếp cận này rời rạc hóa miền mô hình thành các hình học đơn giản kích thước nhỏ, gọi là phần tử (hoặc ô). Tập hợp các phần tử tạo thành lưới (grid/mesh) đại diện cho miền mô hình.
Nếu các phần tử đồng nhất về hình dạng và kích thước (ví dụ hình vuông hay chữ nhật), cách gọi thông dụng là lưới có cấu trúc. Lưới không cấu trúc dùng khi bộ phần tử có kích thước và hình dạng biến đổi. (Thuật ngữ mesh cũng thường được dùng thay cho lưới không cấu trúc.)
3.3.1. Sai phân hữu hạn (Finite Difference)
Phương pháp sai phân hữu hạn là cách trực tiếp nhất để ánh xạ các phương trình vi phân lên từng phần tử trong miền. Cách này thay thế biểu diễn liên tục của phương trình bằng một hệ phương trình dạng sai phân tại một dãy điểm trong miền. Phương pháp thường làm việc trên lưới có khoảng cách điểm đều đặn để tăng hiệu quả tính toán. Tuy nhiên, sự phụ thuộc vào lưới đều khiến nó kém phù hợp với các miền có hình dạng bất quy tắc (ví dụ lòng sông). Các điểm chính của phương pháp sai phân hữu hạn:
- Vị trí tính toán thường đặt tại tâm của mỗi ô.
- Khối lượng được bảo toàn tại từng ô, vì vậy phương pháp bảo toàn khối lượng toàn cục.
3.3.2. Phần tử hữu hạn (Finite Element)
- Phương pháp phần tử hữu hạn lấy các phương trình trường liên tục mô tả vật lý của thủy lực (như đã nêu) và thiết lập phương trình cho từng phần tử. Cách làm là xấp xỉ trường bên trong mỗi phần tử bằng một hàm đơn giản (chẳng hạn đa thức bậc nhất) với số bậc tự do hữu hạn. Cách tiếp cận này cung cấp một mô tả cục bộ xấp xỉ của cơ học dòng chảy trong phần tử dưới dạng một hệ phương trình tuyến tính. Khi lắp ráp đóng góp của tất cả phần tử, FEM tạo ra một ma trận lớn với một hàng cho mỗi ẩn trong hệ (mỗi điểm tính có 3 ẩn). Tuy nhiên, vì các ẩn của một phần tử chỉ bị ảnh hưởng bởi các phần tử lân cận trực tiếp, đa số hệ số trong mỗi phương trình bằng 0 → ta thu được hệ phương trình ma trận thưa, có thể giải bằng các bộ giải ma trận thông dụng.
- FEM là phương pháp rất tổng quát. Nó đã được áp dụng thành công cho nhiều bài toán vật lý và dễ dàng mở rộng sang phần tử bậc cao (bậc hai, bậc ba, …), thay vì chỉ giới hạn ở dạng tuyến tính. Phần tử bậc cao giúp tăng độ chính xác khi biểu diễn các ẩn trong một phần tử. Những điểm chính của FEM:
- Vị trí tính là các vị trí đặc biệt bên trong phần tử, thường gọi là điểm Gauss.
- Phương pháp lặp cục bộ trên toàn bộ phần tử để bảo toàn khối lượng ở mức toàn cục; còn ở cấp phần tử, bảo toàn khối lượng xấp xỉ.
3.3.3. Thể tích hữu hạn (Finite Volume)
Phương pháp thể tích hữu hạn (FVM) giống phần tử hữu hạn (FEM) ở chỗ đều dùng lưới; tuy nhiên nền tảng số của FVM khác hẳn FEM. Trái với FEM, FVM tận dụng việc nhiều định luật vật lý có dạng bảo toàn khối lượng và bảo toàn động lượng. Dựa trên ý tưởng đó, FVM xây dựng công thức gồm các phương trình bảo toàn thông lượng được xác định theo nghĩa trung bình trên phần tử.
Ưu điểm của FVM là chỉ cần tính thông lượng qua các biên phần tử. Ưu điểm này vẫn đúng cho bài toán phi tuyến, khiến FVM đặc biệt hữu ích cho các luật bảo toàn phi tuyến thường gặp trong các bài toán truyền tải (ví dụ vận chuyển bùn cát, chất lượng nước). Tinh chỉnh lưới có thể dùng để bắt các biến thiên thủy lực quy mô nhỏ tương tự như FEM. Tuy nhiên, vì FVM chỉ dựa vào bảo toàn qua cạnh/biên, nó không hưởng lợi từ việc tăng bậc của phần tử (như bậc hai, bậc ba).
Các điểm chính của FVM:
- Vị trí tính toán thường là tâm phần tử hoặc tâm mặt chung giữa các phần tử.
- Phương pháp cưỡng bức bảo toàn khối lượng đúng trong từng phần tử; do đó bảo toàn khối lượng cả cục bộ lẫn toàn cục.
3.4. Các giả định của mô hình
Nhà thống kê George E. P. Box nổi tiếng với câu: “Tất cả các mô hình đều sai; một số mô hình thì hữu ích.” Ý này nhấn mạnh rằng mô hình chỉ là xấp xỉ của thực tế nên không thể khớp hoàn toàn.
Các nhà khoa học và kỹ sư đã xây dựng mô hình 2D với những giả định giản lược để có thể giải được bằng toán. Hiểu các giả định này giúp người mô hình xác định tính phù hợp của mô hình đối với một kịch bản hay một bộ điều kiện cụ thể. Mục 2.3 và 3.2 đã thảo luận các giả định số thông dụng và trình bày giả định trung bình theo chiều sâu. Khi các điều kiện vật lý cần biểu diễn không phù hợp với những giả định trong công thức mô hình, thì mô hình sẽ “sai” nhiều hơn và kém hữu ích hơn.
3.4.1. Giả định về bề mặt được xấp xỉ
Mọi mô hình số đều xấp xỉ địa hình nền của miền quan tâm bằng các nút tính toán (xem Mục 5.3).
Ví dụ: mô hình 1D chỉ tính điều kiện thủy lực tại một điểm đại diện bởi mặt cắt do người dùng định nghĩa; còn mô hình 2D tính điều kiện thủy lực tại các phần tử như đã mô tả.
Một số mô hình có thể biểu diễn biến thiên giữa các nút hoặc bên trong phần tử theo quan hệ phi tuyến để mô tả chi tiết điều kiện thủy lực hơn. Chẳng hạn, nếu một cạnh được biểu diễn bằng một cao độ duy nhất, thì diện tích mặt cắt phía trên cạnh đó đơn giản bằng WSE trừ cao độ cạnh, rồi nhân với chiều dài cạnh. Tuy nhiên, nếu diện tích mặt cắt được biểu diễn bằng đường cong phi tuyến WSE–diện tích, thì diện tích dẫn dòng (conveyance area) theo đường cong đó có thể bắt đầu ở cao trình thấp hơn. Đến một mức nào đó, diện tích mặt cắt qua cạnh theo cả hai cách sẽ hội tụ về cùng một giá trị.
Năng lực này được gọi là “độ phân giải dưới lưới (sub-grid scale resolution)”.
Không có biểu diễn dưới lưới, một phần tử chỉ có thể ướt hoặc khô. Có độ phân giải dưới lưới, mô hình có thể biểu diễn khả năng dẫn dòng một phần (conveyance) đi qua phần tử hoặc qua một cạnh đang ngập một phần. Các mô hình sub-grid tính quan hệ phi tuyến này bằng cách đưa thêm thông tin địa hình chi tiết dọc theo cạnh hoặc trên chính phần tử.
Cần thận trọng khi dùng độ phân giải dưới lưới. Nếu đặt cạnh phần tử không hợp lý, sub-grid có thể thổi phồng khả năng dẫn dòng hoặc tạo dẫn dòng ở nơi lẽ ra không có. Một mô hình số chỉ tính một nghiệm thủy lực cho mỗi phần tử; vì vậy, sub-grid không thể thay thế cho độ phân giải lưới đầy đủ do kích thước, cấu hình, và hướng phần tử quyết định để có nghiệm hợp lý. Cả hai dạng công thức mô hình — có hoặc không có tính năng sub-grid — đều có thể cho kết quả thủy lực chi tiết nếu lưới mô hình biểu diễn đầy đủ các yếu tố khống chế thủy lực quan trọng.
Đối với các mô hình số dùng cho phân tích/thiết kế thủy lực công trình giao thông, lưới cần có độ chi tiết đủ để bắt được các đặc trưng địa hình có ý nghĩa thủy lực và để biểu diễn biến thiên điều kiện dòng chảy (độ phân giải trường dòng). Chương 5 nêu các thực hành khuyến nghị để xây dựng lưới mô hình 2D hiệu quả và chất lượng. Chương 6 trình bày phương pháp rà soát mô hình nhằm bảo đảm mô phỏng phù hợp với mục đích phân tích hoặc thiết kế.
3.4.2. Giả định về vận tốc theo phương đứng
Như đã nêu ở Mục 2.3.1, đa số mô hình 2D thuộc nhóm trung bình theo chiều sâu, nghĩa là mô hình chỉ tính vận tốc trong mặt phẳng ngang và coi đó là vận tốc đã được lấy trung bình theo chiều sâu tại tọa độ (x,y). Các mô hình này hàm chứa giả định rằng phân bố vận tốc theo phương đứng là không đáng kể.
Trong thực tế sông ngòi có thể xuất hiện những điều kiện trái với giả định này, ví dụ quanh công trình (như trụ cầu), tại các đoạn cong gắt, hoặc xuống các dốc tràn. Khi thủy lực thực địa không phù hợp với giả định về vận tốc theo phương đứng, kết quả mô hình 2D cần được xem xét cẩn trọng, vì theo nghĩa chặt chẽ thì chúng “sai” (không biểu diễn hết mọi biến thủy lực trong tình huống phức tạp đó). Tuy vậy, nhờ phương pháp lấy trung bình, các mô hình này vẫn hữu ích.
3.5. Quy trình lời giải
Đầu ra chính của một mô hình 2D gồm ba giá trị tại mỗi vị trí tính toán: mực nước, vận tốc theo phương x, và vận tốc theo phương y. Để tính được các giá trị này, mô hình lắp ráp một hệ PDE biểu diễn bảo toàn khối lượng và bảo toàn động lượng cho từng vị trí tính toán. Giải hệ phương trình này cho ta xấp xỉ của các ẩn tại các vị trí tính. Cách thực hiện phụ thuộc vào cách tiếp cận số của mô hình (xem Mục 3.2).
Không hiếm khi nghiệm tại phần tử được nội suy hoặc lấy trung bình sang vị trí khác để trực quan hóa. Với phần tử hữu hạn, các vị trí tính là điểm Gauss bên trong phần tử chứ không phải tâm/cạnh phần tử; vì vậy cần nội suy để hiển thị. Thể tích hữu hạn và sai phân hữu hạn cũng thường lấy trung bình giá trị từ tâm phần tử hoặc từ cạnh về nút. Việc nội suy/lấy trung bình này tạo ra sai số làm tròn và sai số trung bình tiềm ẩn. Hầu hết người dùng chấp nhận mức sai số này để có thể hiển thị nghiệm như một trường biến thiên liên tục.
3.5.1. Implicit và Explicit
Trong thực hành, các sơ đồ lời giải theo thời gian thường được gọi là explicit hoặc implicit. Một cách phân biệt là xem phương pháp tiến thời gian:
- Sơ đồ explicit tính trực tiếp các ẩn ở bước thời gian kế tiếp từ các đại lượng đã biết ở bước hiện tại.
- Sơ đồ implicit thiết lập một hệ phương trình liên hợp chứa các giá trị ở cả hai bước thời gian (n và n+1) và giải hệ đó bằng các bộ giải ma trận hay lặp.
Các phương trình chi phối trong cơ học chất lưu tính toán thường phi tuyến và có nhiều ẩn, nên phương pháp implicit hầu như luôn cần kỹ thuật lặp.
Ví dụ sau minh họa cách hoạt động của hai phương pháp. Cho biến V tại thời điểm \(t = n\,\Delta t\) (với \(\Delta t\) là bước thời gian, n là chỉ số bước). Mục tiêu là tính V ở bước kế tiếp \(t=(n+1)\Delta t\). Có thể diễn tả đạo hàm theo thời gian tại thời điểm t là \(\Delta = dV/dt\) và viết:
$$ V_{n+1}=V_n+\Delta t\,\Delta .$$
- Với explicit, \(\Delta\) được đánh giá bằng các đại lượng tại bước n.
- Với implicit, một phần của \(\Delta\) được đánh giá bằng các ẩn ở bước n+1, nên phải giải hệ phương trình (ma trận hoặc lặp) để tìm các giá trị mới.
Sơ đồ explicit dùng ít thời gian CPU hơn để có \(V_{n+1}\) nhưng kém ổn định hơn. Có thể thấy sự khác biệt về ổn định khi chia phương trình cho \(\Delta t\):
$$\frac{V_{n+1}-V_n}{\Delta t}= \Delta .$$
Khi \(\Delta t \to \infty\) thì vế trái \(\to 0\), ta được \(0=\Delta\).
- Với explicit, \(\Delta\) chỉ phụ thuộc các giá trị ở bước n, nên không tồn tại nghiệm thỏa điều kiện này → luôn có giới hạn thực tế cho kích thước \(\Delta t\).
- Với implicit, \(\Delta\) chứa các biến của cả hai bước, nên luôn có thể tạo được một nghiệm xấp xỉ (không bảo đảm tốt, nhưng tồn tại), thể hiện tính ổn định của sơ đồ implicit.
Về chi phí tính toán: explicit thường rẻ hơn, dễ lập trình, dễ song song hóa và mỗi bước giải tốn ít công, nhưng thường cần nhiều bước thời gian hơn (bước nhỏ). Implicit có thể cho ít bước hơn theo bậc độ lớn, nhưng mỗi bước nặng do phải giải hệ ma trận/lặp.
Do implicit kết hợp thông tin lân cận theo thời gian (từ bước n+1) vào nghiệm, nó có xu hướng làm trơn/giảm dao động quá độ và làm dịu các gradient thủy lực cục bộ. Nếu các hiệu ứng không ổn định theo thời gian là quan trọng cho thiết kế, một mô hình dùng implicit có thể không phản ánh đủ những hiệu ứng này. Tuy nhiên, với bài toán sông ngòi điển hình, thiết kế lấn chiếm đường sá, hoặc trạng thái ổn định, cả hai loại đều có thể phù hợp. Người thực hành cần đánh giá kết quả thủy lực để bảo đảm mô phỏng đáp ứng nhu cầu thiết kế/phân tích, thường bằng phân tích nhạy cảm hoặc thử các độ phân giải lưới khác nhau (xem Mục 5.3.3).
3.5.2. Điều kiện Courant – giới hạn CFL
Điều kiện Courant–Friedrichs–Lewy (CFL) (Courant, Friedrichs & Lewy, 1928), thường gọi tắt là điều kiện Courant, là điều kiện cần để hội tụ khi giải một số PDE bằng sai phân hữu hạn với sơ đồ thời gian tường minh (explicit). Với các sơ đồ này, bước thời gian phải nhỏ hơn một ngưỡng nêu dưới đây; nếu không, nghiệm có thể mất ổn định hoặc không chính xác/phi vật lý.
Số Courant (đại lượng không thứ nguyên) so sánh vận tốc truyền sóng v (tốc độ truyền thông tin trong miền nghiệm) với chiều dài đặc trưng \(\Delta x\) (kích thước phần tử, thường lấy căn bậc hai diện tích phần tử) và bước thời gian \(\Delta t\): \(C \;=\; v\,\Delta t/\Delta x\).
Các phương pháp explicit đòi hỏi chọn \(\Delta t\) sao cho dòng không tiến xa quá một phần tử trong mỗi bước thời gian. Ràng buộc này liên quan đến độ chính xác vì đa số phương trình sai phân chỉ dùng các đại lượng ở phần tử lân cận. Nếu trong một bước thời gian, “hạt” nước đi xa hơn một phần tử, nó đã sang vùng không ảnh hưởng tới biến đang xét—phi vật lý và gây mất ổn định số. Do đó, để ổn định, với explicit cần C<1 trên toàn miền.
Ngược lại, các phương pháp implicit liên kết tất cả phần tử qua lời giải lặp, cho phép tín hiệu truyền qua lưới. Đổi lại, sự “giao tiếp” giữa các phần tử xa nhau này thường gây giảm dao động/làm trơn biến do các bước thư giãn khi giải hệ liên hợp.
Đa số mô hình số được lập trình ở một trong hai dạng (implicit hoặc explicit) tùy mục đích và lựa chọn của nhà phát triển; một số cho phép người dùng chọn sơ đồ—chi tiết xem trong tài liệu hướng dẫn của mô hình.
Khi nghiệm biến thiên nhanh theo thời gian và chính sự biến thiên này là đối tượng nghiên cứu, dùng explicit thường chính xác hơn (đổi lại phải dùng bước thời gian nhỏ). Với bài toán trạng thái ổn định hoặc nơi biến thiên từ tốn, implicit có thể hiệu quả hơn. Nói chung, yêu cầu chủ chốt là người dùng chọn \(\Delta t\) đủ nhỏ để đạt hội tụ số với mô hình đang dùng.
3.6. Đơn giản hóa mô hình
Do chi phí tính toán tăng cao khi dùng các công thức vật lý phức tạp trong mô hình 2D, nên không hiếm gặp các công thức có giả định giản lược. Người mô hình cần biết rõ công thức mà mô hình cụ thể đang dùng và bảo đảm rằng dạng giản lược đó phù hợp với điều kiện sẽ mô phỏng. Mục này mô tả một số giản lược phổ biến trong các mô hình 2D hiện có.
Dù một giả định giản lược sẽ giới hạn phạm vi áp dụng của công thức, điều đó không có nghĩa là ta không nên dùng nó. Trong một số điều kiện, các hạng tử bị lược bỏ thực sự nhỏ không đáng kể. Ngược lại, nếu các hạng tử ấy đáng kể, không nên dùng công thức giản lược. Nói chung, các ràng buộc về thời gian và tài nguyên tính toán không đủ để biện minh cho việc chỉ vì hiệu quả mà đơn giản hóa công thức. Xem Mục 2.4 để biết thêm các cân nhắc khi chọn mô hình.
3.6.1. Dạng “dynamic wave” và “diffusive wave”
Các cụm “dynamic wave equation” hay “shallow water equation (SWE)” thường được dùng để chỉ hoặc là hệ Navier–Stokes 2D đầy đủ, hoặc là hệ Saint-Venant 1D. Các hệ này khó giải về mặt số, nhưng áp dụng được cho mọi kịch bản dòng chảy trong kênh (tương ứng 2D hay 1D).
Diffusive wave và kinematic wave là dạng giản lược của dynamic wave và đã được dùng rộng rãi trong thực hành kỹ thuật:
- Diffusive wave giả định các hạng tử quán tính (gia tốc) là không đáng kể so với áp lực, ma sát và trọng lực → loại bỏ các hạng tử quán tính khỏi phương trình.
- Kinematic wave còn giản lược thêm: coi áp lực cũng không đáng kể so với ma sát và trọng lực.
Kinematic wave là giản lược nhất, vẫn dùng tốt cho dòng tràn mặt (overland flow), sóng lũ dâng chậm, hoặc kênh có dốc rất lớn không có hiệu ứng dềnh nước (backwater). Dạng này phù hợp khi trọng lực và ma sát cân bằng nhau.
Diffusive wave áp dụng cho mọi điều kiện mà kinematic áp dụng được, và mô tả tốt hơn sự suy giảm/tắt dần của sóng lũ. Nó cũng dùng được cho phổ dốc đáy và chu kỳ sóng rộng hơn kinematic. Diffusive thích hợp khi áp lực có vai trò bên cạnh trọng lực và ma sát ⇒ các lực động lượng (quán tính) vẫn chưa là nhân tố chi phối.
Dynamic wave cần dùng khi cả quán tính (động lượng) và áp lực đều quan trọng, và hiệu ứng dềnh nước không thể bỏ qua (ví dụ kênh dốc thoải có điều kiện khống chế hạ lưu). Các lực này quan trọng trong đa số bài toán thủy lực liên quan đến công trình giao thông; vì vậy cần cân nhắc kỹ khi chọn mô hình 2D và dạng SWE tương ứng.
Trong mô phỏng dòng sông, năng lượng hệ càng cao thì càng nên dùng dynamic wave. Với khu vực cầu và các dòng chảy vận tốc lớn trong lòng kênh, luôn nên dùng dynamic wave equation.
3.6.2. Nền đáy cứng (rigid bed) và nền đáy di động (mobile bed)
Một giản lược khác thường dùng khi tính toán thủy lực là tách thủy lực khỏi vận chuyển bùn cát.
Mô hình nền đáy cứng giả định địa hình dùng cho phân tích thủy lực không di động, tức không thay đổi trong suốt quá trình mô phỏng. Nhờ vậy, các phương trình bảo toàn động lượng được đơn giản hóa đáng kể.
Ngược lại, mô phỏng nền đáy di động làm tăng số biến và số hạng cần giải trong các phương trình bảo toàn khối lượng (cho nước, bùn cát…). Toàn bộ bài toán phức tạp hơn ít nhất một bậc độ lớn.
Mô hình nền đáy di động cho phép địa hình thay đổi dưới tác dụng của lực thủy lực trên bề mặt lòng dẫn hoặc đồng bằng ngập. Cách làm phổ biến là tách tính toán vận chuyển bùn cát thành một bước xử lý xen giữa các bước thời gian thủy lực.
Hiện nay, trong thực hành thủy lực công trình giao thông, người ta chủ yếu dùng phân tích nền đáy cứng, và chỉ áp dụng nền đáy di động cho các trường hợp phức tạp. Tuy nhiên, khi năng lực tính toán tiếp tục tăng, phân tích nền đáy di động sẽ ngày càng phổ biến (xem Mục 8.4).
3.7. Điều kiện biên
Những thử nghiệm đơn giản với mô hình số—chẳng hạn thay đổi độ sâu khống chế tại biên thoát—cho thấy hiệu ứng biên có thể lan truyền vào trong miền mô hình (thậm chí xa bằng nhiều lần bề rộng floodplain). Vì vậy, người mô hình nên đặt các điều kiện biên sao cho càng đơn giản càng tốt và thực hiện phân tích độ nhạy (xem Mục 6.4) để bảo đảm rằng một khoảng bất định hợp lý của điều kiện biên không ảnh hưởng đáng kể đến các điều kiện thủy lực dự đoán trong vùng quan tâm. Do ảnh hưởng này, có thể cần mở rộng ranh miền ra xa hơn khỏi khu vực nghiên cứu.
Tốt nhất, hãy đặt điều kiện biên tại vị trí mà dòng chảy có thể xấp xỉ dạng 1D đơn giản. Nếu không thể, nên kéo dài miền đến một vùng đơn giản nhân tạo để thuận lợi cho việc gán biên (Alemseged, 2007).
Để biết thêm về yêu cầu dữ liệu và cách gán điều kiện biên, xem các Mục 4.3 và 5.4.
Hỗ trợ duy trì trang:
Tôi xây dựng trang này để chia sẻ các tài liệu kỹ thuật cốt lõi trong thiết kế hạ tầng giao thông.
Nếu bạn thấy nội dung hữu ích và muốn góp phần duy trì trang hoạt động bền vững, tôi rất trân trọng mọi sự ủng hộ.