Tin học ứng dụng MIKE21-SW

Bài giảng dùng cho Khóa 56B

Năm học 2017-2018

Bộ môn Quản lý biển và đới bờ
Tầng 3 nhà C1
Trường Đại học Thủy lợi
Số 175 Tây Sơn, Đống Đa, Hà Nội

Biên soạn tài liệu:
GV Nguyễn Quang Chiến

Ấn bản LibreOffice 11/2017. Chuyển bản điện tử dạng Markdown trên StackEdit.

1. Giới thiệu chung

Máy tính ngày càng giúp ích nhiều hơn trong mọi lĩnh vực kĩ thuật ở ba công đoạn: xử lí số liệu, tính toán và biểu thị kết quả. Trong chuyên ngành kĩ thuật biển, tin học được ứng dụng trong các bài toán hải dương học, để xác định dòng thủy triều, phân bố mặn và nhiệt độ, các dự án kĩ thuật bờ biển (truyền sóng để xác định tác động của sóng lên công trình đê biển, đập chắn sóng, hải cảng), xử lý số liệu tính toán sóng cực trị, tính toán ổn định nền công trình, thiết kế các bản vẽ công trình, v.v. Ở trường Đại học thủy lợi, các phần mềm máy tính như Delft3D, MIKE, Plaxis, SWAN, CRESS, Wadibe, GENESIS, v.v. dần được đưa vào giảng dạy, giúp cho sinh viên nâng cao kĩ năng chuyên ngành đồng thời củng cố kiến thức lý thuyết được học.

Trong môn học 2 tín chỉ này, sinh viên được giới thiệu mô hình tính sóng ven bờ. Mô hình sóng có nhiều ứng dụng quan trọng như sau:

• dự báo sóng cho hoạt động hàng hải;
• thiết kế công trình;
• đánh giá các quá trình biến đổi địa mạo, môi trường biển ven bờ.

Chương trình phần mềm được sử dụng là bộ MIKE do công ti DHI Software (Đan Mạch) phát hành. Đây là phần mềm được nhiều công ti tư vấn trong và ngoài nước dùng phổ biến và đánh giá cao. Độ chính xác của MIKE luôn được cải thiện và tiệm cận những kĩ thuật tính toán tiên tiến mới được nghiên cứu. Về tốc độ tính toán, MIKE chạy trên nền Windows chưa nhanh bằng các chương trình khác trên các dòng hệ UNIX/Linux nhưng gần đây (bản MIKE 2012) cũng đã được thiết kế để chạy song song trên nhiều nhân CPU. Giao diện là một thế mạnh của MIKE (cả về giao diện đồ họa lẫn cách mô đun hóa chương trình và nhập xuất nhiều định dạng dữ liệu khác nhau, hỗ trợ tiền và hậu xử lý).

MIKE có giá thành tương đối đắt (~ 1 tỉ đồng cho một bộ gồm các module thông dụng như sóng, dòng chảy ven biển). Hiện tại ở trung tâm tin học của trường, sinh viên thực hành trên bản demo với tính năng tạo lưới địa hình, nhập các thông số tính toán và hiển thị kết quả. Với bản demo, khả năng tính toán bị hạn chế với một miền rất đơn giản. Sinh viên cũng có thể cài bản demo mới nhất (MIKE 2020), lên máy tính cá nhân của mình. (Chọn phiên bản 32/64-bit thích hợp từ https://support.dhigroup.com/download/). Hiện tại ở Trung tâm Tin học ĐH Thủy lợi cài đặt sẵn phiên bản MIKE 2007, và tập tài liệu hướng dẫn những đặc điểm cơ bản của phiên bản này.

Tác giả thấy có hai thiên hướng: hoặc diễn giải quá sâu sắc lý thuyết mô hình sóng, hoặc ghi nhớ máy móc những thao tác sử dụng phần mềm. Cả hai đều không phải mục đích chính của môn học, và cần phải giảm bớt. Thay vào đó, tác giả muốn nhấn mạnh ý tưởng của chương trình phần mềm như một công cụ giúp người dùng thao tác trên một lượng lớn các số liệu hiện có.

2. Dự án

Để tiện cho việc thực hành, sinh viên thực hiện một ví dụ tính toán với huyện Hải Hậu – Nam Định. Vùng bờ Hải Hậu thẳng và không có cửa sông lớn, tác động của sóng rất đáng kể và có lẽ là nguyên nhân chính gây ra xói lở bờ biển mạnh như đã quan sát thấy trong nửa cuối thế kỉ qua.

Để tiện cho việc thực hành, sinh viên thực hiện một ví dụ tính toán với huyện Hải Hậu – Nam Định. Vùng bờ Hải Hậu thẳng và không có cửa sông lớn, tác động của sóng rất đáng kể và có lẽ là nguyên nhân chính gây ra xói lở bờ biển mạnh như đã quan sát thấy trong nửa cuối thế kỉ qua.

Trong nội dung 2 tín chỉ của môn này, sinh viên cần thực hành độc lập, khoanh vùng và tạo lưới tính toán cho bờ biển địa phương, lấy các số liệu địa hình, các yếu tố khí tượng, thủy hải văn như gió, sóng, thủy triều (gọi chung là các “lực đẩy” hay forcing) để xem sóng truyền như thế nào (động thái của hệ thống như thế nào). Sẽ rất thú vị nếu mô hình cho ta thấy được các hiện tượng truyền sóng đã được học qua bài giảng, nhưng trước hết cần phải đề cập đến một số vấn đề mô hình hóa sóng trong §3.

Bài tập
Dùng Google Earth để xác định trên bản đồ vị trí của đường bờ Hải Hậu: kinh-vĩ độ__________________ hay theo hệ UTM là ____________________ m. Xác định độ dài: _________, hướng đường bờ: _________.

Bạn có nghĩ rằng bờ biển thẳng thì các đường đồng mức đáy biển cũng song song với nhau không? Hãy thử di chuột để khảo sát độ sâu và vẽ các đường đồng mức –5 m, –10 m. Ước tính độ dốc đáy biển: _________.

Một mục đích quan trọng của môn học là sinh viên không chỉ thực hành “như một cái máy” mà cần có tư duy giải quyết vấn đề. Cuốn tài liệu này trình bày xen kẽ thao tác máy tính với bản chất của bài toán. Xét cho cùng, phần mềm máy tính rồi sẽ lỗi thời nhưng sự hiểu biết về bản chất vấn đề sẽ có ích để giải quyết những bài toán sau này, có thể bằng phần mềm khác.

3. Mô hình toán sóng ven bờ: sơ lược

Bạn chưa phân biệt rõ mô hình toán, mô hình số trị, và phần mềm? Mô hình nói nôm na là thể thức do người tạo ra để phỏng theo sự vật hiện tượng trên thực tế (mô phỏng về hình dạng, đặc tính ứng xử, sự vận động), nhưng đơn giản hơn cái trong thực tế. Mô hình toán dùng các biểu thức toán học, thể hiện đặc tính ứng xử của vật qua các biến số. Để mô phỏng hiệu quả được các hệ thống phức tạp gồm nhiều thành phần, phương pháp số được áp dụng trong mô hình số trị. Còn phần mềm (chương trình máy tính) bắt nguồn từ việc viết mã lệnh thực hiện cụ thể mô hình số trị đó rồi có thể là đóng gói thành một sản phẩm như MIKE chẳng hạn. Lấy ví dụ sự truyền sóng trên mặt cắt ngang bờ (theo trục x), từ phương trình vi phân cơ bản người ta sai phân hóa trên các đoạn được đánh thứ tự 1, 2, …, i-1, i, … n; mỗi đoạn dài Δx. Mã lệnh được viết bằng ngôn ngữ lập trình, dĩ nhiên là “xấu” hơn biểu thức toán học nhưng lại rất hiệu quả khi thực hiện trên máy tính.

• PT sóng: $\frac{\partial}{\partial x} E C_g \cos \theta = -D$
• Sai phân: (E~i~ C_gi cosθ_i − E_i-1C_g,i-1 cosθ_i-1) / Δx = − ½(D_i + D_i-1)
  → E_i = …
• Mã lệnh: E(i)=(E(i-1)\*Cg(i-1)\*cos(dir(i-1))-0.5\*(Dw(i-1)+Dw(i))\*dx)/(Cg(i)\*cos(dir(i)));

Có nhiều mô hình sóng, được phân loại thành:

• dạng công thức hay đồ thị Hs ~ Tp ~ đà gió ~ thời gian (0-D),
• mô hình truyền sóng trên mặt cắt ngang (1-D, như Wadibe),
• mô hình truyền sóng vùng / mặt bằng (2-D)

Trong mô hình 2 chiều có 2 loại chính:

• 2-D phân giải pha (mô phỏng từng con sóng nhấp nhô): mô hình họ Boussinesq (như MIKE21 BW), phi thủy tĩnh (như SWASH): áp dụng cho sóng ở cảng, có nhiễu xạ mạnh;
  

• 2-D phổ (chỉ mô phỏng đặc trưng thống kê của trường sóng, như SWAN, MIKE21 SW): áp dụng cho các vùng bờ biển, vịnh lớn.
  

Với môn học này, sinh viên sử dụng MIKE21 SW. Hãy sẵn sàng tư duy trong không gian 2 chiều: các đặc trưng sóng như H_s và T_p là chưa đủ, còn phải biết phân bố của chúng thế nào! Thường người ta quy định tọa độ x ngang bờ, y dọc bờ, và θ là góc giữa hướng truyền sóng và pháp tuyến đường bờ (θ = 0, sóng vỗ trực diện vào bờ).

Lưu ý x và y trên rất khác tọa độ địa lý X và Y.

Bài tập

1. Loại mô hình sóng nào sau đây sẽ cho kết quả là bản đồ màu chỉ thị chiều cao sóng ý nghĩa?
   a) mô hình 1-D
   b) 2-D Boussinesq
   c) 2-D phổ

2. Mô hình nào yêu cầu tính toán với độ phân giải cao hơn?
   a) 2-D Boussinesq
   b) 2-D phổ

3. Mô hình Boussinesq thường được dùng để dự báo sóng.
   a) Đúng
   b) Sai

4. Mô hình sóng phổ

Cơ chế động lực của sóng: truyền năng lượng sóng (hiểu là E/ρg), các tham biến: tốc độ truyền sóng C_g, hướng sóng θ. Năng lượng sóng được biểu diễn qua phổ sóng.

Phổ sóng 2 chiều: phức tạp nhưng quan trọng, ví dụ khi có sóng lừng và sóng gió thì phổ 1 chiều sẽ không phân biệt được hướng truyền của 2 thành phần sóng này.

• Nguồn bổ sung năng lượng sóng: gió.
• Tiêu tán năng lượng sóng: “bạc đầu”, ma sát đáy, sóng vỡ.
• Thay đổi năng lượng sóng do tương tác sóng.

Phương trình cân bằng năng lượng sóng:

• cho 1 đơn vị diện tích mặt biển:

mức độ thay đổi NLG = (thông lượng NLG vào − thông lượng NLG ra) + (phát sinh NLG − tiêu tán NLG tại chỗ)

• cho nước sâu:
  $\frac{\partial E{({f,\theta})}}{\partial t} + c_{g,x}\frac{\partial E{({f,\theta})}}{\partial x} + c_{g,y}\frac{\partial E{({f,\theta})}}{\partial y} = S{({f,\theta})}$

• cho nước ven bờ:
  $\frac{\partial E{({f,\theta})}}{\partial t} + \frac{c_{g,x}\partial E{({f,\theta})}}{\partial x} + \frac{c_{g,y}\partial E{({f,\theta})}}{\partial y} + \frac{c_{\theta}\partial E{({f,\theta})}}{\partial\theta} = S{({f,\theta})}$

• cho trường hợp có dòng chảy: sử dụng “hoạt độ sóng” N = E/σ:

$\frac{\partial N{({f,\theta})}}{\partial t} + \frac{c_{g,x}\partial N{({f,\theta})}}{\partial x} + \frac{c_{g,y}\partial N{({f,\theta})}}{\partial y} + \frac{c_{\theta}\partial N{({f,\theta})}}{\partial\theta} + \frac{c_{\sigma}\partial N{({f,\theta})}}{\partial\sigma} = \frac{S{({f,\theta})}}{\sigma}$

Qua kết quả mô phỏng có thể nhận diện các hiện tượng đã được học: khúc xạ, nước nông → những trường hợp đặc biệt của phương trình tổng quát!

5. Tổng quát vận hành chương trình

Bộ phần mềm MIKE được tổ chức thành một hệ thống các mô đun.

• MikeZero là tập hợp các mô đun xử lý số liệu (địa hình, chuỗi thời gian). Trong đó ta chủ yếu sẽ dùng mô đun Mesh Generator để tạo lưới địa hình dạng tam giác.
• Mike21 là tập hợp các mô hình tính toán dòng chảy, sóng và vận chuyển cát. Bên cạnh mô đun SW, mô đun dòng chảy (Flow FM, xem §31) rất quan trọng để tính dòng triều ở các cửa sông, lạch triều và có thể gộp với mô đun SW (tương tác dòng chảy-sóng). Ngoài ra, Mike21 Toolbox (§32) cũng cung cấp tính năng để tạo ra và điều chỉnh số liệu.

Vận hành chương trình qua các bước sau:

• Khoanh vùng tính toán, chọn thời gian tính toán
• Nhập số liệu
  • địa hình: (x y z): bắt buộc để mô hình chạy được
  • thủy động (điều kiện biên): để mô hình cho kết quả
  • các thông số khác (độ nhám v.v.): để hiệu chỉnh, kiểm định KQ
• Chỉ định yêu cầu lấy kết quả (trên toàn vùng hay tại vị trí nào?)
• Tất cả các số liệu địa hình, thủy động lực đều được ghi ra những file chữ ASCII (dạng text), kể cả số liệu mạng lưới vốn có yếu tố hình học.
• Chạy chương trình tính (chạy thử, chạy hiệu chỉnh, chạy kiểm định, chạy các kịch bản). Chương trình (.exe) sẽ đọc các file số liệu nêu trên và tính toán, ghi ra các file kết quả. Lẽ ra các file này đều là dạng text ASCII, nhưng do khối lượng số liệu (2 chiều) rất lớn nên các file thường được nén và mã hóa theo format riêng của MIKE.
• Biểu thị kết quả: các file mã hóa trên được phần mềm đọc và lập thành đồ thị khi chương trình tính kết thúc.

Với lần đầu tiên chạy chương trình, cơ chế hoạt động của chương trình nói chung như sau:

• Ban đầu hệ thống “tĩnh tại”: mặt biển phẳng lặng.
• Tại bước thời gian tiếp theo, các lực đẩy (forcing) tác động lên mô hình làm hệ thống vận động (sóng, dòng chảy).
• Mỗi bước thời gian, tập hợp các đặc tính thủy động của hệ (như trường sóng, trường dòng chảy) được lưu vào các ma trận tồn tại trong bộ nhớ RAM.
• Bước thời gian mới, trạng thái thủy động thay đổi, các ma trận được cập nhật. Bởi vậy muốn không bị mất kết quả tính quá trình thủy động thì phải lưu các ma trận này vào đĩa cứng (chỉ định yêu cầu output).
• Khi dùng thạo, ta có thể dùng kết quả trạng thái thủy động một thời điểm nhất định đã lưu làm mốc để tính toán cho kịch bản mới (hot-start).
• TUY NHIÊN riêng với tính toán sóng người ta cũng có thể quan tâm đến kết quả cuối cùng của trường sóng khi hệ thống đã ổn định (Stationary condition). Khi đó chỉ cần xuất kết quả tại 1 thời điểm (thường là cuối thời gian mô phỏng).

6. Xác định vùng tính toán, lấy số liệu địa hình

Dữ liệu có ở http://coastal-study.uk/Thuc-hanh-MIKE.zip. Giải nén vào thư mục và chọn thư mục này để bố trí chạy MIKE.

Trên thực tế ta phải định rõ quy mô dự án quy hoạch, từ đó mới đi lấy các số liệu địa hình đủ rộng để bao phủ khu vực đó. Tuy nhiên bài thực hành này có tính chất minh họa nên ta làm ngược lại và đã có số liệu địa hình sẵn rồi.

Chuẩn bị các phần mềm phụ trợ như Excel, và cài vào máy một text editor như SciTE hoặc Notepad++, v.v. (tốt hơn Notepad rất nhiều).

Bài tập
Mở file KML thể hiện các điểm đo sâu của khu vực bằng Google Earth. Quan sát sự phân bố điểm đo trong Google Earth và phạm vi của số liệu: __________________________. Độ sâu đạt được đến bao nhiêu ______, so sánh với độ sâu quan sát khi di chuột trong Google Earth.
Mở file địa hình xyz bằng SciTE (trong đó có các điểm độ sâu đúng như thể hiện trong file KML. Cột thứ ba là độ sâu, hai cột đầu là kinh-vĩ hoặc tọa độ X-Y theo lưới chiếu UTM. Có bao nhiêu điểm dữ liệu?

Lưu ý

Cần kiểm tra số liệu xem tọa độ ghi theo kinh vĩ hay hệ UTM?

Trong hệ UTM, bờ biển Việt Nam chủ yếu ở múi 48, chỉ có một phần nhỏ miền Trung (Phú Yên?) ở múi 49. Tọa độ X dao động khoảng > 500.000 m và tọa độ Y tăng dần từ Nam ra Bắc, giá trị khoảng 1-2 triệu m.

Có thể nhập số liệu này vào Excel để kiểm tra (lưu ý tách cột cẩn thận: Data - Text to Column - Delimiter – Space (hoặc đôi khi là comma).

Lưu ý rằng khác phần mềm Delft3D, MIKE đọc z là cao độ (dưới nước số âm) chứ không phải độ sâu dương.

Một file xyz chỉ có 3 cột. Số liệu bắt đầu ngay từ dòng 1. Chấp nhận cách kí hiệu như 5.854E+05 = 585400.

Một số trường hợp cần chỉnh sửa địa hình trong file XYZ: chuyển độ sâu sang cao độ, chuyển gốc cao độ, chuyển từ một tọa độ bình đồ địa phương sang tọa độ X-Y chuẩn của lưới chiếu UTM. Hãy chỉnh sửa số liệu trong Excel rồi copy paste vào SciTE. Lưu file với tên mới có đuôi xyz trong thư mục đã chọn. (Nếu sử dụng Notepad chứ không dùng SciTE, trong khi đặt tên file phải kèm trong dấu nháy kép, chẳng hạn “diahinh.xyz”.

Cũng cần nói thêm rằng có những nguồn số liệu không cho dưới dạng xyz mà dưới dạng ma trận ô độ sâu (chẳng hạn format *.asc). Một file dạng asc thì ngoài ma trận độ sâu còn có những thông tin khác như kích thước của từng ô và tọa độ của một ô cố định (như ô góc Tây Nam của ma trận có tọa độ bao nhiêu?). Có nguồn cung cấp số liệu dạng này (như GEBCO www.gebco.net); định dạng asc có thể chuyển sang xyz được.

7. Nhập số liệu địa hình

Từ giao diện chính của MIKE chọn biểu tượng New (giấy trắng). Trong hộp thoại mới, phần bên trái chọn MikeZero, bên phải tìm chọn Mesh Generator. Chọn Long/Lat hoặc UTM phù hợp (48). Ấn OK; một cửa sổ mới xuất hiện có lưới chiếu và tọa độ. Đây là module Mesh Generator.

Chọn menu Data → Import Scatter Data… Chọn file xyz đã lưu. Chọn Long/Lat hoặc UTM 48, sẽ thấy file đó hiện ra trong bảng. Ấn Apply và quan sát thông tin thống kê được từ file số liệu.

Bài tập
Tổng số điểm đo bằng bao nhiêu; so sánh với số dòng trong file xyz.
Các phạm vi mặt bằng Xmin, Xmax, Ymin, Ymax của số liệu thế nào? Phạm vi cao độ thế nào (Zmin - Zmax)? Cao độ đã được chỉnh đúng chưa?

Đóng cửa sổ và quan sát những điểm đo mới tạo được. Sự phân bố các điểm này có giống trong Google Earth không? Về màu sắc, các điểm xanh ứng với nước sâu và điểm đỏ ứng với nước nông hoặc “trên bờ”.

Chỉnh sửa nếu thấy có chỗ điểm đo quá dày: Data → Scatter Set Reduction … Limit scatter data set based on scatter point proximity với Distance threshold là khoảng cách tối thiểu 2 điểm gần nhau, chẳng hạn 100 (m). Đặt tên file địa hình mới (.xyz) sau đó ấn Reduce. Ngược lại, vùng điểm quá thưa có thể “bổ sung điểm” (Add Scatter Data, trên thanh công cụ).

* Quản lý lớp bản đồ

Trước đây khi chưa có các công cụ như Google Earth thì tính năng GIS (“bản đồ”) của MIKE giúp ích rất nhiều cho việc mô hình hóa vì qua đó ta xác định được vị trí đường bờ, tuyến công trình v.v. Ở ví dụ thực hành này vì đường bờ rất đơn giản và xuất hiện rõ trong số liệu xyz nên chức năng này chưa cần thiết ngay; nhưng với các vùng nước trong đất liền như sông, hồ… thì vị trí đường bờ rất quan trọng và cần tận dụng tính năng này của MIKE.

Nếu đã có ảnh bản đồ giấy/ảnh vệ tinh của khu vực với dạng file JPG, PNG, GIF v.v. thì tiến hành các bước sau:

• Từ trong Mesh Generator chọn Options → Import Graphic Layers. Tạo một lớp ảnh mới rồi chọn File … Trong mục Edit cần chỉnh các thông tin bao gồm tọa độ góc dưới phía trái của tấm ảnh (tọa độ Long/Lat hoặc XY của UTM). Khi lấy ảnh ta cần lưu ý chính xác tọa độ này) và bề rộng, chiều cao của ảnh (tính theo độ hoặc mét). Với Google Earth, việc xác định tọa độ góc ảnh có thể thực hiện bằng cách bật lưới kinh vĩ (Show Grid) hoặc đánh dấu các điểm trên bản đồ và ghi lại tọa độ các điểm đó.
• Chuyển sang thẻ [Overlay Manager] xuống dưới cùng danh sách, chọn lấy lớp ảnh mới lập nên rồi đẩy lên trên cùng ().
• Khi đóng hộp thoại, nếu chưa nhìn thấy hình ảnh nền thì cần chọn Options → Workspace … rồi nhập các tọa độ X,Y thích hợp vào.

8. Khoanh vùng tính toán

• Nếu để thiết kế công trình thì vùng tính toán phải rộng hơn hẳn công trình (cả về phía biển lẫn dọc ven bờ), để tránh tình trạng công trình gây nhiễu động thủy lực truyền ra biên.
• Vùng tính toán vẫn nằm trong phạm vi số liệu địa hình. (Thực ra là ngược lại, với một vùng tính toán định trước thì phải tìm cách thu thập số liệu địa hình thật nhiều, bao phủ được vùng tính toán.)
• Với mô hình truyền sóng bờ biển thẳng đơn giản thế này, nên khoanh vùng chữ nhật có “cạnh” song song với bờ biển, 1 cạnh sát bờ biển cũng được. Biên gần sát bờ biển này đôi khi người ta tỉ mỉ vẽ đường khúc khuỷu như trên bản đồ, nếu là bến cảng thì điều này cần thiết nhưng với đường bờ tự nhiên thì không cần. Khi đó đường bờ sẽ được MIKE tự nội suy qua các điểm độ sâu địa hình gần nhau có độ sâu chuyển từ âm lên dương.

Các cạnh của vùng khoanh có thể là đường cong. Lưu ý tùy theo điều kiện địa hình người ta không nhất thiết khoanh 4 cạnh hình chữ nhật. Chẳng hạn trong 1 vịnh người ta có thể khoanh 2 “cạnh”: một cạnh cong phía biển và 1 cạnh zích zắc theo đường bờ.

Sử dụng tính năng trên thanh công cụ: Chọn điểm, cạnh, vùng. Vẽ mới điểm, cạnh, vùng. Di chuyển điểm, xóa điểm. Chọn, tạo đường Break.

Bắt đầu khoanh vùng bằng việc chọn vẽ mới từng cạnh. Thao tác đúng: click trái để bắt đầu và kết thúc bằng 1 click phải → Insert end node. Nếu ở điểm cuối bạn click trái rồi mới click phải sẽ bị lỗi trùng điểm. Một cạnh có thể là 1 đoạn thẳng hoặc 1 đoạn gấp khúc.

Kết thúc một cạnh rồi thì bắt đầu cạnh tiếp theo có thể ở điểm cuối cạnh trước (nhưng không cắt qua cạnh trước).

Kiểm tra lại mã màu: Một cạnh có 2 điểm đầu vuông, màu tím than; các điểm trung gian của cạnh là chấm tròn màu đỏ. Nếu khoanh vùng 4 cạnh thì phải thấy có 4 điểm vuông tím than ở 4 góc.

Lưu ý
Lỗi phổ biến:

1) Thông báo “Triangulation failed due to duplicated point”. Do các điểm trên cạnh trùng nhau, thường là trùng 2 điểm cuối trên 1 cạnh. Phái xóa bớt điểm.
2) Không thấy xuất hiện lưới tam giác. Có lẽ do vùng khoanh bị “hở”. Phóng to (Zoom) vào để kiểm tra và chỉnh lại. Lưu ý các đoạn góc của vùng khoanh phải là đỉnh vuông màu tím than.

9. Đánh số biên

Lưu ý
Rất quan trọng; thiếu bước này thì dù lưới đã lập xong nhưng không thể đặt điều kiện biên để chạy chương trình được.

Vùng có bao nhiêu cạnh thì có bấy nhiêu số hiệu biên. Số 0 và 1 dành cho biên trên đất liền. Số từ 2 trở lên đánh cho các biên dưới nước/ngoài biển.

Để đánh số, chọn 1 cạnh trên vùng. Click phải → Properties. Điền vào số hiệu. Lưu ý rằng điển xong một lượt các cạnh, khi quay trở lại thì hai số hiệu hai bên có thể khác đi như số hiệu ở giữa (Arc attribute) chuẩn là được.

10. Tạo lưới

Ta sẽ tạo lưới tam giác trong phạm vi vùng vừa được khoanh. Chọn menu Mesh → Triangulate. Hộp thoại mở ra có 3 thông số quan trọng:

• Maximum element area _____ m² : diện tích tối đa của một ô lưới. Nếu để tự động, MIKE sẽ “lười” và tạo ô lưới rất thô. Cần phải hạn chế diện tích lớn nhất này. Cạnh ô lưới trung bình có thể xấp xỉ $1,5 \times \sqrt{d.\mathit{tích}}$. Riêng với hệ tọa độ kinh vĩ (LONG/LAT) thì bạn cần nhập vào diện tích tính theo độ vuông (deg²). Tùy theo vĩ độ của vùng dự án, 1 độ vuông có diện tích khác nhau, từ 1,23×10¹⁰ m² (lớn nhất, ở xích đạo) đến 1,18×10¹⁰ m² (ở vĩ độ 16, miền Trung Việt Nam), đến 0,88×10¹⁰ m² (ở vĩ độ 45, các nước châu Âu).
• Smallest allowable angle _____ (độ): góc nhọn nhất của tam giác ô lưới, khuyến cáo không nên < 26, tam giác càng “nhọn” thì càng dở; hoàn hảo nhất là các tam giác đều (60 độ).
• Maximum number of nodes _______ : số điểm tối đa trong lưới. Lần này là để MIKE không “chăm” tạo quá nhiều điểm, sẽ tính toán rất lâu. Cố gắng hạn chế không quá 50000 phần tử.

Ấn Generate và phía dưới thông báo số ô (elements) và số điểm (nodes). Số nào nhiều hơn, tại sao?

Bài tập Để xây dựng lưới cho vùng ven bờ có công trình đập mỏ hàn, người ta muốn có lưới mịn (cạnh ô lưới cỡ 10 m) thì cần nhập vào diện tích tối đa ô lưới bằng bao nhiêu? (*Lưu ý rằng độ dài bao giờ cũng dễ hình dung hơn diện tích, bởi vậy trong các báo cáo kĩ thuật hãy cố gắng đề cập độ dài cạnh ô lưới.)

Đóng hộp thoại quan sát lưới tạo được, thường thì sát bờ do có độ khúc khuỷu thì ô lưới mịn, ra xa ô lưới sẽ thô hơn. Tuy vậy, để kiểm soát tốt hơn độ phân giải ta cần thực hiện kĩ thuật phân vùng như ở §11.

Làm trơn lưới: Chọn Mesh → Smooth Mesh. Làm trơn 10 lần chẳng hạn. Kết quả sẽ được lưới có các tam giác “đẹp” hơn. Một ví dụ cho lưới chưa làm trơn (hình trái) và sau khi làm trơn (hình phải) cho thấy các tam giác phía dưới hình vẽ đã được điều chỉnh cân đối hơn.

Cuối cùng, nếu làm hỏng lưới, hãy chọn Mesh → Delete Mesh, rồi vẽ lại lưới.

11. Phân vùng lưới tính

Một quy tắc tính toán là kích thước ô lưới phải tương đồng với độ sâu. Vùng nước sâu → ô lưới thô, nước nông → lưới mịn. Thường cạnh ô lưới > độ sâu.

Để thỏa mãn điều đó, độ phân giải lưới phải có sự chuyển tiếp từ mịn → thô từ bờ đi ra. Ví dụ trên hình vẽ có ít nhất 3 vùng với độ phân giải khác nhau.

Để xác lập các vùng này trước hết phải kẻ thêm các đường phân chia miền tính thành mấy phần khép kín, kề sát nhau. Để kẻ đường phân chia này (hai đầu là nút vuông màu tím than, cần chuyển từ điểm màu đỏ. Kích chuột chọn điểm đỏ, kích phải, chọn Vertex to Node. Mỗi miền được đánh dấu một đa giác (chấm một điểm xanh vào trong miền đó, bằng công cụ ). Sau đó đặt thuộc tính độ phân giải cho từng miền tương tự như ở §10. Cụ thể, ấn chọn dấu đa giác, click phải chuột → Properties, bỏ tick “Exclude from Triangulation”, chọn độ phân giải theo m² hoặc theo độ² tùy theo hệ tọa độ dùng cho file xyz.

Một quy tắc thực dụng là các miền lân cận nhau thì diện tích ô lưới chênh lệch không quá 10 lần. Như vậy hai tam giác lớn và nhỏ của hai vùng cạnh nhau có độ dài cạnh gấp nhau 3 lần là được. Các phiên bản mới của MIKE cho phép tự động hóa làm mịn lưới căn cứ vào thông tin độ sâu qua menu Data → Refine Mesh; tuy nhiên điều này chỉ thực hiện được sau bước nội suy triển điểm độ sâu như ở §12.

Cách phân vùng này cũng được áp dụng trường hợp có đảo. Khi đó đảo là 1 đa giác không được chia lưới (Properties → Exclude from Mesh).

Lưu ý Lưu lại file lưới (*.mdf). Sau này nếu đã có địa hình, nếu muốn làm mịn lưới theo yêu cầu mô phỏng có thể quay lại dùng file này. Khi thử nghiệm các phương án tạo lưới khác nhau nên save as ra file khác…

Tính năng Go to (biểu tượng ống nhòm) trong module Mesh Generator cho ta xem vị trí những ô lưới lớn nhất, nhỏ nhất, hay một khu vực nhỏ có tọa độ cụ thể, thuận tiện cho việc điều chỉnh lưới địa hình.

12. Triển địa hình lên lưới (Nội suy)

Để lưới dùng được cho việc tính toán thì trên lưới phải có thông tin về độ sâu. Cụ thể, MIKE sẽ bắn độ sâu từ các điểm đo lên các nút lưới (đỉnh tam giác).

Chọn menu Mesh → Interpolate… Phương pháp nội suy thường được sử dụng là Natural Neighborhood (nội suy từ trung bình độ sâu của các điểm lân cận, với trọng số cao hơn cho những điểm gần hơn). Có thể lựa chọn Extrapolate (ngoại suy ra ngoài phạm vi cụm điểm số liệu).

Việc nội suy thành công sẽ tạo nên bản đồ màu. Hãy kiểm tra lại độ sâu. Bạn có thể lấy được con số giá trị cụ thể của độ sâu bằng cách click biểu tượng “vùng” rồi khoanh một đa giác trên hình.

Đây là lưới độ sâu *.mesh (khác lưới tam giác *.mdf). Lưu lại file này bằng cách Mesh → Export … File mesh cũng là một file text, hãy dùng SciTE để xem xét cấu trúc của nó. Lưu ý đừng chỉnh sửa tùy ý và save lại, sau này MIKE có thể sẽ không mở được.

Trong MIKE21, file mesh cũng là một text file xem được bằng SciTE và nếu cần, có thể chỉnh sửa được nhưng phải rất thận trọng. Sau khi chỉnh sửa nên lưu 1 bản copy rồi mở thử trong MIKE xem được không.

13. Chỉnh địa hình

Chỉnh địa hình không phải là “ăn gian” số liệu. Có 2 trường hợp:

• Khi có nạo vét đáy biển hoặc đổ cát, xây dựng một số công trình ngầm.
• Khi số liệu đo có những điểm biến đổi đột ngột, không phù hợp với địa hình, địa mạo khu vực
• (Riêng trong mô hình flow mô phỏng ngập lụt vùng nước nông, vị trí biên mô hình rất dễ gây lỗi nếu có những bãi cạn; phải cắt bỏ những bãi cạn này để dòng chảy xuôi thuận.)

Cách chỉnh địa hình:

• Chỉnh thủ công: click phải → Select Area, khoanh vùng, sửa số và kiểm tra lại bản đồ
• Chỉnh theo công thức toán: Sau khi khoanh vùng, chọn menu Edit → Calculator. Trong biểu thức này, CurrItem chính là độ sâu.
• Chỉnh xong click phải → clear selection. Kiểm tra xem địa hình mới xuất hiện có ổn không?

Một số biểu thức ví dụ (phần xám có sẵn rồi không cần gõ vào):

• CurrItem = -5 → san/cào bằng đến cao trình -5 m.
• CurrItem = CurrItem + 2 → tôn cao đồng thời thêm 2 m.
• CurrItem = MIN(-15, CurrItem) → đảm bảo cao trình từ -15 m trở xuống (chỗ nào nông hơn 15 m thì nạo vét xuống sâu đến 15 m).

Bài tập Một cảng được xây dựng ven bờ. Phạm vi 2 km dọc bờ × 1 km vươn ra xa bờ. Đập phá sóng có cao trình đỉnh +3 m. Bể cảng được nạo sâu -15 m. Luồng tàu đảm bảo độ sâu 15 m ứng với mực nước = 0. Hãy chỉnh sửa địa hình rồi Export ra một file mesh mới.

Lưới tính toán (lưới không gian) là yếu tố quan trọng, đặc thù của các mô hình số trị. Lưới tính toán phải đảm bảo độ phân giải cần thiết (khi chạy mô hình, các đặc trưng sóng sẽ được tính toán trên lưới, mỗi ô lưới sẽ lưu giữ một giá trị cho đặc trưng sóng nhất định như năng lượng sóng E). Khi tính toán trên lưới tam giác, những giá trị nào ở các ô lưới cạnh nhau sẽ liên hệ qua phương trình vi phân; bởi vậy các ô lưới càng “đều đặn” thì càng tốt (xem thêm §17). Trong nhiều nghiên cứu CFD (tính toán thủy động), độ phân giải của lưới chia cũng là một thông số quan trọng ảnh hưởng đến sai số của mô phỏng.

14. Mô đun sóng phổ SW

Trong giao diện MIKE chọn biểu tượng New. Ở hộp thoại mới, bên trái chọn MIKE21, bên phải chọn Spectral Wave FM. Mô đun mới được mở ra có một danh sách các dữ liệu như trên hình. Các mục có dấu tick xanh lá cây là “được”, mục có dấu gạch chéo màu đỏ là “bị lỗi”. Cần sửa hết lỗi mới chạy được.

• Domain: miền tính toán, nhập file mesh. Kiểm tra minimum cut-off depth (những cao độ nào lớn hơn số này sẽ được coi là vách dựng đứng, nước không tràn qua). Kiểm tra số hiệu các biên (Boundary Names) và nếu cần đổi tên dễ hiểu.
• Time: thời điểm bắt đầu - tùy theo số liệu thực đo để kiểm định như thế nào, số bước tính toán v.v. Các số cần nhập vào là ngày giờ bắt đầu mô phỏng, số bước thời gian và độ dài bước thời gian tính toán. Với sóng gió cần thời gian mô phỏng để sóng phát triển hoàn toàn. Ví dụ chọn 10-15 bước thời gian, mỗi bước Δt = 1 giờ.

Chi tiết các yếu tố trong module sẽ được trình bày trong các mục sau.

15. Phương trình cơ bản

Các phương trình cơ bản (Basic equations) được xét theo hai khía cạnh:

• Theo phổ:
  • Phương pháp Directionally Decoupled Parametric, về cơ bản là tách hướng, tham số hóa theo miền tần số f và các biến được chọn là mô men bậc 0 và bậc nhất của phổ N(f).
  • Phương pháp Fully Spectral tính toán theo phổ đầy đủ N(f, θ).
• Theo thời gian
  • Phương pháp Quasi-stationary không xét biến thời gian, mỗi bước tính đều được thực hiện theo phương trình ổn định.
  • Phương pháp Instationary tính toàn bộ diễn biến theo thời gian.

Tùy theo yêu cầu tính toán khác nhau mà lựa chọn phù hợp.

• Nếu chỉ tính sóng lừng, chọn Directionally Decoupled.
• Nếu chỉ sóng gió, Fully Spectral.
• Nếu chọn mô phỏng 1 cơn sóng lớn (trường hợp cực hạn) hoặc nếu các con sóng riêng lẻ coi là độc lập, chọn Quasi Stationary.
• Các trường hợp khác:
  • Miền nhỏ có thể chọn Quasi Stationary, kết quả chấp nhận được.
  • Miền lớn: chọn Instationary (tính lâu hơn nhiều).

16. Rời rạc về phổ

Ta biết rằng phổ sóng tại một địa điểm, một thời điểm (x, y, t) là hàm của 2 biến (xem hình vẽ mặt cong ở §4). Khi tính toán cần rời rạc phổ sóng này.

• Frequency discretation: nên chọn logarithmic. Công thức là f_n = f₀ cⁿ với f₀ là tần số nhỏ nhất, n số đoạn chia, c là thừa số (factor). Cần xem sóng ở phạm vi chu kì nào để điều chỉnh các giá trị f₀, c và n. Giá trị mặc định cho phép mô phỏng ở quãng 1,8 s → 20 s. So sánh với sóng ngoài khơi điển hình có chu kì từ 4 s → 25 s và sóng ở vùng nước khép kín 2 → 3 s. Nếu sóng có phổ tần số hẹp thì có thể giảm n và c. Một hình vẽ minh họa cho các đoạn chia có ở bài tập tiếp sau.
• Directional discretation: Lựa chọn full (cả vòng tròn, 360°) cho phép mô phỏng sóng từ nhiều hướng có sóng gió, sóng lừng kết hợp. Tuy nhiên nếu chỉ có sóng lừng thì chọn một hình quạt bao lấy hướng sóng chính. Độ phân giải nên lấy Δθ = 5°. Với sóng gió, độ phân giải cho phép kém hơn Δθ = 10°. Theo nghiên cứu được Goda (2000) tổng hợp, với một hướng sóng chính trung tâm MWD là θ₀, luôn có sự phân tán năng lượng sóng hai bên hướng θ₀. Với sóng gió, có 85% năng lượng tập trung trong hình quạt từ θ₀ − 30° đến θ₀ + 30° (phân bố năng lượng như đường cong ở hình vẽ bài tập tiếp theo). (Lưu ý rằng góc sóng ở đây đều là góc so với hướng Bắc. Nêu rõ cách chọn Δθ_min, θ_max, n_θ.)

Như vậy qua phép rời rạc hóa, miền phổ sẽ được chia thành các phần tử năng lượng hình “cung vành khăn” (liên hệ với hình vẽ ở §4). Lưu ý rằng lưới chia này là trên không gian phổ tại một điểm, khác với trên không gian x-y, lưới chia dạng tam giác.

Bài tập
Với điều kiện sóng gió có chu kì sóng biến đổi từ 5 s đến 20 s và hướng sóng chính là đông bắc (NE), các phổ theo tần số và theo hướng như hình vẽ, hãy nhập các thông số thích hợp cho Discretization.

• Phân biệt sóng lừng và sóng gió: có thể tính gộp hoặc phân theo “ngưỡng tần số cố định” hoặc “ngưỡng tần số động”.

17. Kĩ thuật giải

Phương pháp thể tích hữu hạn (finite volume method) được dùng trong MIKE21 với các ô lưới lưu giữ biến ở vị trí trung tâm. Thông lượng được tính theo sơ đồ sai phân “đầu nguồn” (upwind) với độ chính xác bậc nhất.

Với tính sóng Instationary: tính với sơ đồ sai phân hiện (phương pháp Euler) để tìm giá trị của các biến lượng ở bước thời gian mới. Bước thời gian Δt phải thỏa mãn điều kiện CFL:

$\mathit{CFL} = \left| {c_{x}\frac{\Delta t}{\Delta x}} \right| + \left| {c_{y}\frac{\Delta t}{\Delta y}} \right| + \left| {c_{\sigma}\frac{\Delta t}{\Delta\sigma}} \right| + \left| {c_{\theta}\frac{\Delta t}{\Delta\theta}} \right| \leq 1$

Thường chọn giá trị CFL = 0,8. Giá trị Δt sẽ được chương trình tự động tính và do đó sẽ liên tục thay đổi khi mô hình đang chạy.

Với tính sóng quasi-stationary, có thể chọn phương pháp Modified Newton-Raphson. Nghiệm ổn định nhận được và quá trình lặp kết thức khi đạt tới số lần lặp tối đa định trước, hoặc khi sai số giữa hai lần lặp kế tiếp đã giảm xuống dưới mức định trước.

18. Điều kiện mực nước, dòng chảy

• Nếu tính toán trong điều kiện thiết kế, hãy chọn Specify water level variation và format là Constant. Nhập vào giá trị mực nước thiết kế.
• Trong trường hợp mực nước có biến đổi theo thời gian thì chọn Varying in time, constant in domain. Khi đó, số liệu quá trình mực nước này nên được trích xuất từ một mô hình thủy triều từ trước.
• Trong điều kiện bãi biển thẳng, dòng triều không đáng kể thì ở Current conditions chọn No current variation.
• Tùy chọn Soft Start Interval có ý nghĩa như sau: để đặc tính dòng chảy biến đổi từ 0 đến giá trị chỉ định, cần một thời gian để “khởi động”.

19. Gió

Gió đóng vai trò lực đẩy (forcing) trong mô hình, tuy không phải điều kiện biên nhưng tác động nhiều đến điều kiện thủy động lực. Trường gió có thể chọn không đổi, chỉ thay đổi theo thời gian, và thay đổi theo cả thời gian lẫn trong không gian. Với trường hợp sau này cần chỉ định file trường gió 2 chiều (trên lưới đều, dfs2 hoặc trên lưới tam giác, dfsu).

• Với mô hình sóng phổ toàn phần, chọn cơ chế tương tác khí quyển - đại dương hoặc là coupled (động lượng truyền từ gió → sóng không những phụ thuộc vào gió mà vào cả sóng “độ gồ ghề của mặt biển”, hay uncoupled (động lượng truyền từ gió → sóng chỉ phụ thuộc vào tốc độ gió).
• Với mô hình thông số phổ, cơ chế tạo sóng chỉ dựa theo hệ thức kinh nghiệm (SPM 73, SPM 84, Kahma & Calkoen, Jonswap).

Với bài toán gió bão, trường gió nên được tạo bằng công cụ MIKE 21 Toolbox → Generating Wind Fields, xem thêm §32.

20. Chuyển năng lượng sóng

Khi sóng truyền đi, phổ sóng có thể biến dạng do tương tác giữa từng con sóng.

• Ở nước sâu: quadruplet interaction (tương tác 4 sóng): tính toán rất lâu, chỉ nêu dùng cho vùng nhỏ
• Ở nước nông: triad interaction (tương tác 3 sóng), hệ số tương tác mặc định là 0,25.

21. Tiêu tán năng lượng sóng

• Ma sát đáy (Bottom Friction). Có thể khai báo theo một trong các cách: hệ số ma sát (mặc định C_fw = 0,0075 m/s), hệ số ma sát không thứ nguyên (mặc định f_w = 0,0212)¹, độ nhám Nikuradse (mặc định k_n = 0,04 m, tuy nhiên có thể giảm đi khi tính cho vùng gần bờ); kích thước hạt (mặc định d₅₀ = 0,00025 m). Số liệu có thể đồng đều hoặc biến đổi trên miền không gian. Ma sát đáy tham gia như một số hạng nguồn S trong phương trình cân bằng năng lượng.
• Sóng vỡ (Wave Breaking). Hệ số Alpha kiểm soát tốc độ tiêu tán năng lượng (mặc định = 1); hệ số gamma kiểm soát độ dốc sóng (mặc định = 1). Hệ số gamma cũng có thể được tính theo công thức thực nghiệm. Theo Ruessink và nnk. (2003), γ = 0,76kh + 0,29. Theo Nelson (1994), γ = F_c / (22 + 1,82F_c), trong đó F_c = g^1,25 H^0,5 T^2,5 h^−1,75 là hệ số phi tuyến.
• Sóng bạc đầu (Whitecapping). Hệ số Cdis chi phối tốc độ tiêu tán năng lượng, DELTAdis chi phối trọng số tiêu tán năng lượng trong phổ sóng.

22. Điều kiện ban đầu

• Zero Spectra: ban đầu năng lượng sóng = 0 trên toàn miền
• Tính theo công thức kinh nghiệm: có thể theo JONSWAP, SPM 1973 nước sâu hoặc nước nông.
• Khai báo từ một file trường dfsu sẵn có: thường dùng cho hotstart, chạy từ kịch bản chạy dở trước đó.

23. Điều kiện biên

Có một số loại điều kiện biên trong MIKE:

• Biên đóng kín (closed boundary): dùng cho đất liền
• Biên sóng: biên ngoài khơi. Có thể …
  • Khai báo theo tham số: H_s , T_p, MWD (hướng sóng trung tâm), và đặc trưng cho độ phân tán hướng sóng là chỉ số n (Version 1) hoặc độ lệch chuẩn về hướng, DSD (Version 2).
  • Khai báo theo phổ sóng: phổ năng lượng E(f,θ) hoặc phổ hoạt độ N(σ,θ). Với phổ hoạt độ có thể lấy toàn bộ hình dạng hoặc các thông số momen phổ là m₀ và m₁.

Dạng biên sóng có thể không đổi hoặc biến đổi (theo thời gian, theo dọc biên). Đơn giản nhất là dạng không đổi. Phức tạp hơn một chút là biến đổi theo thời gian, khi đó phải nhập file số liệu dfs0. Phức tạp nhất là biến đổi theo dọc biên và cả theo thời gian (file số liệu dfs1). Cả hai file dfs0 và dfs1 cần được biên tập (New > Times Series (dfs0) hoặc Profile Series (dfs1) > Blank…); lưu ý phải khai báo các đặc trưng sóng theo tên gọi quy định: wave height, wave period, wave direction, spreading index, spreading factor (DSD). Chú ý chọn đơn vị phù hợp. Nội dung số liệu được nhập vào bảng (có thể copy paste từ Excel). Có công cụ giúp nội suy số liệu (Tools – Interpolation) cho phép điền giá trị trống ví dụ khi giá trị ở biên bị khuyết. Riêng với biên tập dfs1, có thể chuyển giữa các đại lượng khác nhau bằng cách chọn Edit – Items > …

• Biên ngang (lateral): thường dùng để đặt 2 biên chạy theo mặt cắt ngang bờ. Lưu ý rằng biên loại lateral này không phát ra năng lượng sóng nên sẽ hình thành hai dải khuất sóng. Chính vì vậy phạm vi mô hình phải rộng để vùng khuất này không ảnh hưởng vào vùng dự án).

Bài tập
Thực hiện mô phỏng một kịch bản (lưu với tên file KB1.sw). Các thông số điều kiện sóng ngoài khơi là: Hs = 1 m, Tp = 8 s, hướng sóng SE (bao nhiêu độ?). Hãy nhập ĐK biên phù hợp. Đồng thời, hãy chọn dạng mô phỏng (§15) cùng các thông số (§16-22).

Thử nhập các biên số liệu biến đổi theo hướng dẫn của GV. (Lưu ý bản demo chỉ cho phép nhập ít hơn 50 dòng số liệu thời gian.)

24. Chỉ định xuất kết quả

Do số lượng dữ liệu tạo ra trong quá trình chạy MIKE 21 thường lớn, nên máy sẽ mặc định không lưu kết quả này. Muốn có số liệu kết quả, ta phải chỉ định qua mục Output. Có 3 dạng kết quả chính:

• Dạng điểm: để lấy số liệu diễn biến một yếu tố theo thời gian. Cần chỉ định tọa độ của một (hoặc một vài) điểm trích dữ liệu.
• Dạng tuyến: để lấy sự biến đổi một yếu tố theo một mặt cắt, hoặc lưu lượng qua một mặt cắt. Chỉ định tọa độ đầu và cuối tuyến cũng như nêu rõ bao nhiêu điểm trung gian trên tuyến đó.
• Dạng vùng: để lấy trường yếu tố theo thời gian. Chọn các yếu tố đặc trưng sóng để xuất kết quả. Quan trọng nhất có lẽ là H_s và vectơ vận tốc truyền sóng. Riêng đại lượng Radiation Stress (ứng suất phát xạ) rất cần thiết để làm đầu vào cho mô hình kết hợp sóng - thủy động lực.

Các kiểu file kết quả ứng với 3 dạng trên lần lượt là *.dfs0, *.dfs1 và *.dfsu. Ấn nút […] và chọn thư mục lưu, đặt tên file (không gõ chữ tiếng Việt có dấu).

25. Chạy chương trình

Chương trình chỉ chạy được khi các số liệu vào không bị lỗi (toàn bộ dấu tick xanh). Tuy nhiên ngay cả khi số liệu ok thì vẫn kiểm tra xem bạn đã chỉ định việc ghi kết quả (Output) chưa.

Ngoài những yếu tố như dung lượng đĩa phải đủ để ghi kết quả tính toán v.v. thì quan trọng nhất là thời gian thực hiện tính toán. Với một hệ thống máy nhất định tốc độ tính toán phụ thuộc nhiều về độ phân giải lưới không gian và độ phân giải lưới phổ.

Khi thực hiện chạy (chọn menu Run – Start simulation …), phía dưới cửa sổ số liệu sẽ xuất hiện thanh hiển thị số % tiến trình chạy, trên thanh trạng thái (Status Bar) sẽ hiện thời gian ước lượng còn lại đến khi chạy xong, số bước tính toán đã thực hiện và thời điểm tính toán trong mô phỏng kịch bản.

Trong bản Demo có hạn chế về số lượng ô lưới của miền tính, thường là vài trăm. Nếu vượt quá số giới hạn này thì sẽ báo lỗi “CheckLicense: No. of nodes exceeded for demo-version” và “Abnormal completion”. Còn nếu quá trình chạy bình thường thì thông báo dưới cùng sẽ là “Normal run completion”. Trong cả hai trường hợp, có thể kiểm tra thêm thông tin trong file *.log cùng thư mục với file mô phỏng.

Bài tập Chọn file xuất kết quả như ở §24. Chạy kịch bản KB1.sw. Sau khi chạy xong, thử lại kĩ thuật hotstart (§22) như sau: khai báo chính file dfsu kết quả vừa chạy ra. Sau đó đổi tên file kịch bản thành KB11.sw, đổi tên các file kết quả mới. Sau đó chạy lại rồi so sánh kết quả kịch bản KB11 với KB1.

26. Quan sát và biểu thị kết quả

Quan sát kết quả giúp ta nắm bắt được quá trình tổng thể và có thể giúp phát hiện các lỗi khi mô phỏng, nếu có. Chọn mục Output - View để xem kết quả dạng trường (có thể phải tua đến thời gian cuối). Những đại lượng vô hướng như H_s sẽ cho kết quả bản đồ màu. Đại lượng vec tơ sẽ cho biểu đồ vec tơ mũi tên. Khi đó chọn Data – Options → Vectors. Biểu thức giá trị vector thường được chọn sẵn. Với cách bố trí vector,

• Vectors at element mesh: mỗi ô có 1 vector, vùng gần bờ lưới dày các vecto sẽ sít lại.
• Vectors interpolated to structured mesh: vector dàn đều ra, khi đó chọn được số vector theo trục ngang và theo trục dọc. Hình minh họa có khoang màu và vector dàn đều như vậy.

Hãy di chuyển và phóng to hình, nếu cần, để xem sự phân bố các yếu tố sóng trên hình, hướng truyền sóng (đặc biệt các hiện tượng khúc xạ, nhiễu xạ), kiểm tra điều kiện bên chỉ định đúng chưa v.v. Nếu cần, có thể trích xuất đường quá trình thời gian tại một điểm cụ thể chọn trên bản đồ (bằng cách click biểu tượng sau đó click vào vị trí cần trích. Chọn Data - Time series by coordinates … (dĩ nhiên ta có thể bổ sung thêm điểm). Quan sát chiều cao sóng dần đạt đến ổn định, nếu có.

Nếu có kết quả trích xuất dạng tuyến, chọn Output - View và nhận xét phân bố sóng dọc theo mặt cắt và quá trình biến đổi sóng mặt cắt theo thời gian.

Biểu diễn kết quả là công đoạn rất quan trọng sau khi chạy xong mô hình. Nó cho phép diễn giải trực quan khối lượng lớn các con số được mã hóa trong file kết quả. Biểu diễn có thể trên màn hình hoặc in. Với cách thứ nhất, có thể phải chỉnh lại thang màu (Chọn Colors - Edit Current Palette …). Nên vẽ viền đen giữa các khoang màu bằng cách View → Isolines. Sau khi căn chỉnh cân đối và rõ ràng thì chụp màn hình lưu file (ví dụ kiểu file *.PNG). Xem ví dụ như hình trường sóng dưới đây.

Với hình thức hai, chọn các nét rõ ràng và sử dụng dạng file đường nét (file vector); ví dụ như hình các đường đồng mức mực nước ở cuối trang. Một số kiểu file là *.EPS, *.EMF. Một số công cụ phần mềm cho phép tinh chỉnh các file vector này, như Inkscape, LibreOffice Draw hay Adobe^® Illustrator.

Riêng với mô hình sóng thì kết quả từng điểm cũng quan trọng. Module MIKE 21 SW cho phép xuất kết quả phổ sóng; biểu đồ tròn được vẽ bằng cách dùng Plot Composer. Đây là công cụ tạo nhiều dạng đồ thị. Chọn New → Mike Zero → Plot Composer. Chọn Insert Plot Object. Các dạng Basic Graphics cho phép các đồ thị cơ bản đã xem được trực tiếp từ kết quả. Ngoài ra, Advanced Graphics cho ta vẽ các đồ thị kết nối các yếu tố khác nhau; chẳng hạn biểu đồ hoa sóng, biểu đồ phổ, biểu đồ tương quan. Chất lượng bản vẽ tạo bởi Plot Composer cũng tốt hơn, cho phép tạo hình nét phục vụ in ấn.

27. Gỡ lỗi

Nếu kết quả chạy ra dường như rất lạ, hãy kiểm tra lại:

• Lưới địa hình đã được chưa?
• Kiểm tra các file kết quả xem mô hình có hội tụ không và có chỗ nào có vấn đề phải chia lại lưới không?
• Tăng số bước lặp, Number of Iterations, để hội tụ được (khi có vẻ gần hội tụ)
• Thay đổi dung sai để cho phép hội tụ, nếu kết quả chạy có vẻ ổn định.
• Tắt bỏ tính toán nhiễu xạ để giảm hệ số làm trơn (trong khoảng từ 0-1, thử 0.9) và tăng số bước
• Thay đổi điều kiện ban đầu thành ‘Zero Spectra’.

28. Hiệu chỉnh và kiểm định mô hình

Trong khâu hiệu chỉnh, cần thay đổi giá trị các thông số để mô hình cho kết quả tính toán phù hợp với kết quả thực đo. Sự phù hợp (hay cố gắng để giảm thiểu sai số) được thể hiện qua nhiều tiêu chí. Chẳng hạn trong tài liệu SWAN (2000) của Delft Hydraulics có đề cập những tiêu chí sau (với kí hiệu x_i là các số liệu thực đo, y_i là các số liệu tính toán, N là tổng số điểm số liệu, x là giá trị trung bình):

• Sai số dự tính:

  • Độ lệch: BIAS = y – x
  • Độ phân tán: Phương sai = $\frac{1}{N-1}{\sum{({x_{i}–y_{i}})}^{2}}$

• Sai số trung bình:

  • Sai số tuyệt đối trung bình MAE = $\frac{1}{N}{\sum{|{y_{i} - x_{i}}|}}$
  • Sai số căn quân phương RMSE = $\sqrt{\frac{1}{N}{\sum{({x_{i}–y_{i}})}^{2}}}$

• Sai số tương đối

  • Chỉ số phân tán (scatter index): SCI = RMSE/|x|
  • Chỉ số hiệu năng mô hình: MPI = 1 − RMSE/RMS_change trong đó RMS_change giống như RMSE chỉ khác là các giá trị quan sát ban đầu x_i0 được dùng thay cho các giá trị tính toán y_i. MMPI có ưu điểm hơn SCI ở chỗ dùng các thay đổi quan trắc được của sóng để đánh giá chất lượng của những thay đổi tính được.
  • Chỉ số hiệu năng vận hành OPI = RMSE / x_i0

• Chỉ số phù hợp

  • Wilmott (1981, 1984, 1985) đề xuất các chỉ số phù hợp có xét đến tỉ số giữa sai số dự tính thực và sai số dự tính tiềm năng:
    $d_{1} = 1 - \frac{\sum{|{y_{i} - x_{i}}|}}{\sum{\lbrack{{|{y_{i} - \overline{x}}|} + {|{y_{i} - x_{i}}|}}\rbrack}}$
    và
    $d_{2} = 1 - \frac{\sum{|{y_{i} - x_{i}}|^2}}{\sum{\lbrack{{|{y_{i} - \overline{x}}|} + {|{y_{i} - x_{i}}|}}\rbrack^2}}$

• Hệ số xác định (coefficient of determination), R², còn gọi là chỉ số hiệu năng Nash-Sutcliffe (NSE). Chỉ số này được dùng phổ biến trong các mô hình thủy văn, mô hình dòng triều:
  $\mathit{R}^2 = 1 - \frac{\sum{({y_{i} - x_{i}})}^{2}}{\sum{({\overline{x} - x_{i}})}^{2}}$

• Hệ số tương quan (Pearson) có thể tính theo công thức:
  $r = \frac{{\sum{({x_{i} - \overline{x}})}}{({y_{i} - \overline{y}})}}{\sqrt{\sum{({x_{i} - \overline{x}})}^{2}}\sqrt{\sum{({y_{i} - \overline{y}})}^{2}}}$
  Hệ số tương quan tiến gần 1 thể hiện quan hệ chặt chẽ giữa kết quả đo và kết quả tính, và khi biểu diễn trên đồ thị điểm chấm x~y, sẽ có dạng tập trung sát một đường thẳng chéo góc 45°. Khi mô hình hóa chế độ sóng (wave climate) thường sử dụng cách kiểm định theo tương quan này, ví dụ Cavaleri & Sclavo (2006).
  

• Kiểm định bằng quan sát: Cách này thường có giá trị về khía cạnh khoa học, để kiểm tra tính chặt chẽ của các mô hình toán truyền sóng (ngang bờ), ví dụ một mô hình cân bằng năng lượng sóng của Ruessink và nnk. (2003).

Với mô hình sóng SW có thể chỉnh các thông số vật lý sau (ngoài các thông số toán học như độ phân giải lưới, độ phân giải phổ):

• Độ nhám đáy biển: tăng độ nhám ở vùng nước nông thường dẫn đến tăng độ tiêu hao năng lượng sóng, làm giảm chiều cao sóng.
• Hệ số sóng vỡ do nước nông (γ) làm giảm chiều cao sóng trong vùng sóng vỡ, còn α ảnh hưởng đến độ tiêu hao năng lượng. Ngoài ra, tham số ngưỡng độ dốc sóng vỡ cũng có thể được điều chỉnh (riêng đối với mô hình tham số tách hướng).
• Các tham số chi phối mức tiêu tán năng lượng do sóng bạc đầu, C_ds và δ, có thể được điều chỉnh (riêng đối với mô hình phổ toàn phần). Nếu vùng nước sâu nhận thấy sóng tính toán thiên nhỏ thì có thể giảm C_ds để thu được kết quả tốt hơn.

Trong giai đoạn kiểm định mô hình, ta lấy bộ thông số được chọn qua khâu kiểm định, sau đó chạy mô hình cho một thời đoạn khác và so sánh với số liệu thực đo; đồng thời đánh giá độ phù hợp. Không được phép chỉnh thông số mô hình trong giai đoạn này, song sự sai khác giữa tính toán và thực đo có thể cho phép rút ra những nhận xét khả năng vận dụng mô hình vào bài toán cụ thể.

Bài tập Dùng bảng tính như Google Sheet (https://docs.google.com/spreadsheets/u/0/) hoặc Excel để thực hiện tính sai số RMSE, chỉ số Nash, hệ số tương quan r cũng như các chỉ số khác. Thực hành với số liệu do GV cung cấp.

29. Chia kịch bản

Để chạy nhiều kịch bản ta cần xem xét: những thông số nào giữ nguyên (địa hình?), những thông số nào thay đổi? Nếu cần, copy số liệu gốc thành nhiều thư mục, mỗi thư mục ứng với một kịch bản và chạy từng kịch bản. Kết quả cuối cùng được tập hợp lại và so sánh với nhau.

• Số liệu chuỗi thời gian dễ dàng so sánh bằng cách copy cột số liệu kết quả từng kịch bản vào Excel và vẽ.
• Số liệu dạng trường có thể chồng xếp lên nhau, tính ra hiệu số và biểu diễn hiệu số này, chẳng hạn bản đồ hiệu số H_s giữa hai kịch bản. Thực hiện xem kết quả như Mục 26. Sau đó bật menu Edit – Calculator. Mở thêm file dfs của kịch bản còn lại để có biển, chẳng hạn tên là Var2. Trong mục biểu thức (Expression), gõ một biểu thức ví dụ CurrItem – Var2. Các biểu thức có thể có dạng thứ nguyên khác với đại lượng CurrItem và do vậy cần dùng hàm void để khử thứ nguyên một cách thích hợp, chẳng hạn (CurrItem – Var2)/void(Var2). Ngoài ra, để tránh chồng lũy tích kết quả, cần chọn: Only Selected Timestep.

Bài tập Kịch bản 1 đã tính ở §25. Hãy bổ sung kịch bản 2 (từ file KB1.sw, Save as file mới KB2.sw), với các thông số: H_s = 2 m, T_p = 10 s, hướng sóng E. Sửa tên file kết quả KQ2.dfsu. Chạy mô hình và so sánh kết quả với kịch bản trước.

30. Kiểm tra độ nhạy

Các kịch bản cũng có thể được dùng để kiểm tra độ nhạy của từng thông số. Nhiều số liệu đầu vào, do độ bất định lớn, cũng cần được xem xét và điều chỉnh. Có ý kiến cho rằng đó cũng là những tham số để hiệu chính, nhưng quan điểm của người viết là chúng nên được dùng để đánh giá độ nhạy mà thôi.

Kinh nghiệm trong chuyên ngành cho thấy một số đại lượng có độ bất định như sau (Soulsby 1998):

• độ sâu nước h ± 5%
• chiều cao sóng H_s ± 10%
• chu kì sóng T ± 10%
• hướng sóng MWD ± 15°
• vận tốc dòng chảy U ± 10%
• đường kính hạt D ± 20%

Bài tập Biểu đồ của độ nhạy một thông số có dạng như hình vẽ, trong đó giả sử tại kịch bản gốc, với hai giá trị thông số X₁* và X₂* (ví dụ như chiều cao sóng và chu kì tại biên) ta tính được kết quả Y* (ví dụ như chiều cao sóng tại điểm công trình ở sát bờ). Thực hiện chỉnh hai tham số X₁ ± 15% và X₂ ± 8% ta tính được các kết quả mới. Nhận xét về các độ nhạy kết quả theo sự thay đổi của từng thông số 1 và 2 ?
Bây giờ, hãy làm với số liệu cụ thể của KB2. Chọn thay đổi X₁* = H_s và X₂* = T_p đã đưa vào điều kiện biên (§29). Chọn Y là điểm cố định ở gần bờ và tính ra các độ nhạy.

31. Module thủy động Flow-FM

Việc tính toán mô phỏng dòng chảy rất quan trọng đối với những trường hợp vùng biển có vũng vịnh, lạch triều, vùng cửa sông. Dòng chảy đóng vai trò then chốt trong các quá trình biến đổi địa mạo, truyền chất ô nhiễm, và có tương tác với sự lan truyền sóng. Sinh viên có thể tự tìm hiểu và thực hành; ở đây chỉ nêu những điểm khác biệt so với mô hình sóng SW.

• Tidal potential (thế thủy triều): nhập các thông tin đặc trưng cho dao động mực nước trong miền tính toán

• Bed resistance: với các công thức của Chezy/Manning

• Eddy viscosity: độ nhớt xoáy. Có thể chọn giá trị hằng số (cỡ một vài m²/s), hoặc tốt hơn là dùng công thức tính như Smagorinsky, $\nu_{t} = C_{s}^{2} \Delta x \Delta y\sqrt{\left( \frac{\partial u}{\partial x} \right)^{2} + \left( \frac{\partial v}{\partial y} \right)^{2} + \frac{1}{2}\left( {\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}} \right)^{2}}$

  trong đó C_s là một hệ số còn u, v là vận tốc (trung bình độ sâu) theo các hướng x và y.

• Wave radiation: ứng suất bức xạ sóng, có thể lấy kết quả từ file dfsu lưu giữ biến Radiation Stress khi chạy MIKE SW. Sóng giữ vai trò như “lực” tham gia phương trình động lượng dòng chảy.

Các thông số vật lý thường dùng để hiệu chỉnh module Flow-FM là:

• Độ nhám ở đáy (Chezy hoặc Manning). Số Chezy nói chung cỡ khoảng 50 m^0,5s^-1 với đa số các vùng biển và thay đổi không nhiều theo độ sâu. Số Manning² M = 1/n trong đó n là con số độ nhám tra trong các bảng thủy lực kênh hở thông dụng. Đơn vị của M là m^1/3s^-1.
• Hệ số nhớt xoáy (công thức Smagorinsky). Giá trị mặc định là C_s = 0,28. Sự thay đổi hệ số này có thể làm thay đổi động thái của thủy vực ở chỗ có gradient vận tốc lớn như sát bờ vách đá hay nơi có dòng chảy mạnh. Một gợi ý khác của DHI [MIKE21-HR] là trong trường hợp ngập lụt do nước dâng trong bão, thay vì chọn Smagorinsky thì dùng hằng số xoáy dạng Flux-based, ước lượng bởi công thức: 0.02Δx Δy/Δt (m²/s).

Trong quá trình chạy module Flow-FM, có thể xảy ra lỗi (“blow up”) khi trên miền xuất hiện vận tốc hoặc độ sâu quá lớn. Khi đó, cần kiểm tra trường dòng chảy/mực nước tại bước thời gian cuối cùng để nhận diện vị trí có giá trị dị thường và có thể cần chỉnh lại để địa hình không biến đổi quá gấp và đường bờ không quá khúc khuỷu.

32. Công cụ phụ trợ

Như đã thấy muốn có số liệu “chuẩn” phục vụ chạy mô hình, cũng như tổ chức trình bày kết quả đầu ra, thì cần thực hiện rất nhiều công đoạn xử lý. Chẳng hạn, trên cơ sở biết định dạng (format) các file xyz mà MIKE sử dụng, ta đã dùng Excel để chỉnh lý số liệu. Cách làm này có ưu điểm là tận dụng được kiến thức sẵn có của Excel mà người dùng có thể hiểu rõ hơn về số liệu. Nhưng cũng có cách làm khác là dùng các công cụ phụ trợ, có dạng biểu tượng “hình chóp da cam” như Data Extraction, Plot Composer, Mike 21 Toolbox.³ Đây là những tập hợp các chương trình nhỏ rất tiện dụng để xử lý số liệu của MIKE.

Data Extraction. Đây là công cụ trích xuất số liệu, nghĩa là tạo một “tập con” số liệu từ kết quả sau khi chạy một kịch bản. Chẳng hạn, từ một file kết quả dfsu, ta muốn: (i) Cắt ngắn kết quả này, chỉ lấy một số thời đoạn quan tâm, (ii) chỉ lấy một vùng nhỏ cần quan tâm trên toàn miền tính, (iii) chỉ định một đường và một số điểm trích số liệu mà ta chưa khai báo khi chạy; tất cả đều có thể thực hiện bằng Data Extraction. Về thao tác, ban đầu ta chỉ định file dfsu cần trích. Sau đó chỉ định cách rút gọn số liệu; việc này khá giống với chỉ định xuất kết quả ở §24.

Mike 21 Toolbox. Một số chương trình hữu ích trong công cụ này bao gồm: tạo sóng đều, sóng ngẫu nhiên (Generation of regular/random waves), tạo phổ sóng (Generate wave energy spectrum), tạo trường sóng trong bão (Generate cyclone wind field). Dưới đây chỉ chọn giới thiệu công cụ tạo trường bão. Để làm điều này, trước hết cần chọn file lưới tên đó sẽ tính trường áp, trường gió. Tiếp theo cần điền vào một bảng gồm 7 cột số liệu: thời gian (giờ) từ khi bắt đầu tính bão, tọa độ x của tâm bão, tọa độ y của tâm bão, bán kính gió lớn nhất (km), vận tốc gió lớn nhất (m/s), áp suất tại tâm bão (hPa hay mbar), áp suất tự nhiên (hPa).

Một vài thông số khác có thể cần được nhập vào khi chọn các mô hình bão khác nhau. Về công thức tính bão, có thể xem tài liệu của Nghiêm Tiến Lam (2008).

Tài liệu tham khảo

Cavaleri, L., Sclavo, M. (2006) The calibration of wind and wave model data in the Mediterranean Sea, Coastal Engineering 53, 613–627.

DHI (2012) MIKE 21 Mesh Generator Step-by-Step Training Guide.

DHI (2015) MIKE 21 Wave Modelling, Spectral Wave FM, Short Description.

DHI (2007) MIKE 21 SW User guide.

DHI (2007) MIKE Zero Preprocessing & Postprocessing User Guide.

DHI UK (2012) MIKE 21 Quick Start Guide.

DHI MIKE 21 FlowFM Manual.

DHI [MIKE21-HR] Hints and recommendations in applications with significant flooding and drying.

Goda (2000) Random Waves and Design of Maritime Structures, 2nd ed., World Scientific.

Hearn, C. (2008) Dynamics of Coastal Models. Cambridge.

Holthuijsen, L. (2005) Waves in Coastal and Oceanic Waters. Cambridge.

Lam, N. T. (2008) Tính toán nước dâng do bão. Hướng dẫn thực hành kĩ thuật bờ biển, ĐH Thủy lợi. link.

Nelson, R.C. (1994) Depth limited design wave heights in very flat regions. Coastal Engineering 23, 43–59.

Roelvink, D., Reniers, A. (2011) A Guide to Modeling Coastal Morphology, World Scientific.

Ruessink, B.G., Walstra, D.J.R., Southgate, H.N. (2003) Calibration and verific­ation of a parametric wave model on barred beaches, *Coastal Engineering *48, 139–149.

Soulsby (1997) Dynamics of Marine Sands. Thomas Telford.

UCAR (2013) Nearshore wave modeling. Online course, https://www.meted.ucar.edu/oceans/nearshore

---

1. Hệ số ma sát sóng có thể được tính dựa theo độ nhám và biên độ dịch chuyển phần tử nước ở đáy, a₀; chẳng hạn f_w = exp(5,5(k_n/a₀)^0,2 − 6,3). ↩︎

2. Thật ra là hệ số Strickler, https://en.wikipedia.org/wiki/Manning_formula. ↩︎

3. Mesh Generation (§7) cũng có biểu tượng tương tự nhưng đó là công cụ quan trọng và được trình bày riêng. Plot Composer được trình bày ở §26. ↩︎