$\mathrm{TE} = \left( \frac{\partial f}{\partial x} \right)_j - \frac{f_{j+1} - f_j} {\Delta x} = -\frac{\Delta x}{2} \left( \frac{\partial^2 f}{\partial x^2} \right)_j -\frac{(\Delta x)^2}{3!} \left( \frac{\partial^3 f}{\partial x^3} \right)_j - \cdots$

Bài tập. Viết biểu thức sai số cắt cụt cho phương pháp sai phân lùi. Cho thấy rằng mặc dù hai phương pháp này có cùng độ chính xác cấp 1 nhưng các sai số cắt cụt là khác nhau.

Thống nhất (consistency). Một sơ đồ sai phân được gọi là thống nhất nếu nó đảm bảo sai số cắt cụt $\to 0$ khi bước tính toán $\to 0$. Tính thống nhất đảm bảo cho sơ đồ sai phân thực sự biểu diễn được PVR.

Ổn định. Một sơ đồ được gọi là ổn định nếu giá trị sai số của đại lượng không bị khuếch đại trong mọi giai đoạn tính toán. Một ví dụ tính toán độ ổn định sẽ được minh hoạ với phương pháp Euler trong mục 2.3.

Hội tụ. Sơ đồ sai phân hội thụ nếu phương trình sai phân đủ sát với nghiệm thực của PVR. Hội tụ là điều kiện mạnh hơn thống nhất (hội tụ = thống nhất + ổn định), tuy nhiên trong thực tế khi thiết lập các sơ đồ sai phân, nói chung hội tụ và thống nhất là đồng nghĩa nhau; khi càng chia nhỏ bước thời gian và bước không gian, nghiệm số trị sẽ hội tụ về nghiệm thực với một sai số cắt cụt hữu hạn.

Ta có thể thiết lập công thức sai phân bằng cách “trực giác” như đạo hàm bậc nhất bằng hiệu giá trị chia khoảng cách, và trong nhiều trường hợp, cách này vẫn đảm bảo tính thống nhất. Chẳng hạn, Chapra [@Chapra1997] đã dùng cách đơn giản để xây dựng công thức đạo hàm bậc hai:
$\displaystyle{ f''_j \approx \frac{f'_{j+} - f'_{j-}} {\Delta x} = \frac{\frac{f_{j+1} - f_j}{\Delta x} - \frac{f_j - f_{j-1}}{\Delta x} } {\Delta x} = \frac{f_{j-1} - 2f_j + f_{j+1}}{\Delta x^2} }$; kết quả thu được giống như thực hiện khai triển Taylor.

Tuy nhiên, không thể lạm dụng cách này để sai phân hoá những phương trình phức tạp. Trong mọi trường hợp, khai triển chuỗi Taylor là cần thiết để kiểm tra xem sai số tiệm cận 0 khi $\Delta x \to 0$. Xét ví dụ sau đây [@Ramshaw2011]:

Ví dụ. Xét trường hợp miền tính gồm các điểm $\{ x_1, x_2, ..., x_j, ..., x_N \}$. Tại phía trái, biên vật lý nằm ở vị trí $x_{+}$, trung điểm giữa $x_1, x_2$. Cần xác định giá trị đạo hàm $f''$ tại điểm $x_{2}$. Nếu biên vật lý đặt tại $x_1$ thì ta đã có thể viết $f''_2 \approx (f_1 - 2f_2 + f_3)/ \Delta x^2$, với độ chính xác cấp hai. Nhưng bây giờ phải biến đổi công thức này để tận dụng giá trị $f_{+}$. Dường như có thể dễ dàng sử dụng thông tin: $x_+$ là trung điểm $x_1, x_2$ vì vậy $f_1 = 2f_+ - f_2$. Từ đó: $f''_2 \approx (2f_+ - 3f_2 + f_3)/ \Delta x^2 \; (*)$. Tuy nhiên, cách này dẫn đến sự không hội tụ! Thật vậy, khai triển chuỗi Taylor quanh điểm $x_2$ cho ta: $\begin{gathered} f_3 = f_2 + f'_2 \Delta x + \frac{1}{2} f''_2 \Delta x^2 + ... \\ f_+ = f_2 - f'_2 (\frac{1}{2} \Delta x) + \frac{1}{2} f''_2 (\frac{1}{4} \Delta x^2) + ...\end{gathered}$ Nhân đôi hai vế PT dưới sau đó cộng PT trên để triệt tiêu số hạng đạo hàm bậc nhất ta được: $2f_+ + f_3 = 3f_2 - \frac{1}{4}f''_2 \Delta x^2$. Vậy nếu sử dụng cách khai triển như (*) thì sẽ xấp xỉ được đạo hàm bậc hai, mà thật ra: $(2f_+ - 3f_2 + f_3)/ \Delta x^2 = \frac{3}{4}f''_2 + O(\Delta x)$, tức là chưa xấp xỉ được giá trị đúng của đạo hàm. Tuy vậy, vẫn có thể xấp xỉ được đạo hàm với độ chính xác cấp 1. (Hãy nhân PT trên gần nhất với 4/3, bạn đọc tự thực hiện như một bài tập.)

Bài toán điều kiện đầu trong PVT

Bài toán này có ý nghĩa dự báo. Cho phương trình vi phân $du/dt = f(u,t)$ với điều kiện đầu $u_\mathrm{đầu}$ cho trước. Một công thức xấp xỉ sai phân đơn giản là: $(u^{n+1} - u^n)/\Delta t \approx f(u^{n}, t^n)$, từ đó cho phép ta dự đoán $u^{n+1} = u^n + \Delta t f(u^n, t^n)$. Đây là phương pháp Euler tiến.

Lập chương trình tính toán giá trị hàm $u$ theo biến $t$; cho điều kiện đầu $u_{(t=0)} = u^0$ và hàm $f$. Quan sát động thái của kết quả dự tính khi $\Delta t$ thay đổi. Khi giảm $\Delta t$, đến một lúc nào đó nghiệm sẽ đột ngột nhảy vọt đến giá trị cực lớn/nhỏ. Ghi lại giá trị $\Delta t$ này. Ta nói phép sai phân đã mất ổn định (xem đoạn giải thích tiếp theo).

Bây giờ hãy xét sơ đồ sai phân sai phân Euler “lùi”: $u^{n+1} = u^n + \Delta t f(u^{n+1}, t^{n+1})$. Cả hai phương pháp Euler tiến và lùi đều có độ chính xác cấp 1, song chúng khác nhau về độ ổn định. Để đơn giản, hãy xét bài toán: $du/dt = au$ trong đó $a$ là hằng số. Phương trình này có nghiệm giải tích $u = u^0 \exp(at)$ trong đó $u^0 = u_{(t=0)}$ là điều kiện đầu.

Để tính toán ổn định, ta dựa theo hệ số khuếch đại, G = $u^{n+1} / u^n$. Với phương pháp Euler tiến, G = $1 + a \Delta t$. Để đảm bảo điều kiện ổn định, phải có $\vert G \vert < 1 \Rightarrow \Delta t < - 2/a$.

Bài tập. Hãy xác định điều kiện ổn định cho phương pháp Euler lùi.

Để cải thiện phương pháp Euler, trước hết ta thấy rằng giá trị đạo hàm tại điểm đầu ($t^n, u^n$) hay điểm cuối ($t^{n+1}, u^{n+1}$) đều có thể không biểu diễn tốt độ dốc của đồ thị hàm số $u$. Từ đó, có hai nhóm kĩ thuật cải thiện sau:

• Dự đoán - hiệu chỉnh. Một phiên bản đơn giản của kĩ thuật này là phương pháp “hình thang”, với độ dốc lấy bằng trung bình hai độ dốc đầu và cuối đoạn. Tại bước 1, tính được giá trị dự đoán $u^\ast = u^n + \Delta t f(u^{n}, t^{n})$, giống như Euler tiến, nhưng giá trị này được dùng tiếp. Bước 2, chỉnh lại giá trị vừa dự đoán $u^{n+1} = u^n + \Delta t f(u^\ast, t^{n+1})$.
• Trung điểm. $u^{n+1/2} = u^n + \Delta t f(u^{n+1/2}, t^{n+1/2})$.

Hãy lập chương trình tính theo phương pháp Euler tiến, và hai kĩ thuật cải thiện nêu trên. Lưu ý rằng, trong những cách này $f(u^{n+?}, t^{n+?})$ sẽ cần có giá trị $u^{n+?}$ tính trước, khác với phương pháp Euler tiến.

Bài toán điều kiện biên trong PVT

Xét PVT biểu diễn bài toán điều kiện biên sau: $-d^2 u/dx^2 + qu = f(x)$ trên miền $a \le x \le b$ với các điều kiện biên $u = \alpha$ tại $x = a$ (ĐK biên dạng Dirichlet) và $du/dx + \beta u = \gamma$ tại $x = b$ (ĐK biên dạng Robin). Với những PT dạng đơn giản, đạo hàm phân li được, có thể giải bằng phương pháp bắn (shooting): Bắt đầu với ĐK biên phía $x=a$, giải PT theo chiều $x$ tăng dần (tương tự như ta đã làm với $t$). Giải đến kết quả tại $x=b$ thì so sánh với điều kiện biên tại đó và điều chỉnh lại nghiệm sao cho giá trị thu được tại $b$ đúng bằng ĐK biên [@Chapra2005].

Tuy vậy, các PT vi phân thường phức tạp và nên được giải bằng cách sai phân. Theo đó có thể viết (công thức sai phân từ bài [fwd-bak-ctr]):

$\frac {-u_{j-1} + u_j - u_{j+1} } {(\Delta x)^2} + q u_j - f_j = 0$

Từ đó, có thể viết lại dưới dạng phương trình đại số tuyến tính: $u_{j - 1} -(2 + q (\Delta x)^2) u_j + u_{j+1} = -f_j (\Delta x)^2$ Và hình thành nên ma trận gồm $N$ ẩn ($j = 1, \dots, N$) với $N-2$ phương trình. Để có thể giải được, cần bổ sung hai điều kiện biên đã nêu ở đề bài. Do vậy, PT đầu trong hệ trở thành $u_1 = \alpha$ còn PT cuối trở thành $-u_{N-1} + (1 + \beta) u_N = \gamma \Delta x$.

Đặt $p = -(2 + q (\Delta x)^2)$ và $r_j = -f_j (\Delta x)^2$ để tiện biểu diễn các PT “giữa”, ta có hệ phương trình ba đường chéo. Ở đây các hệ số bằng 0 thì được bỏ trống. $\begin{bmatrix} 1 & 0 & & & & \\ 1 & p & 1 & & & \\ & 1 & p & 1 & & \\ & & . & . & . & \\ & & & 1 & p & 1 \\ & & & & -1 & 1+\beta \\ \end{bmatrix} \begin{Bmatrix} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{N-1} \\ u_N \end{Bmatrix} = \begin{Bmatrix} \alpha \\ r_2 \\ r_3 \\ \vdots \\ r_{N-1} \\ \gamma \Delta x \end{Bmatrix}$

Mặc dù đã có nhiều thuật toán tốt để giải hệ phương trình đại số tuyến tính nhưng với các ma trận thưa như hệ 3 đường chéo thì có kĩ thuật giải riêng, như thuật toán Thomas.

Bài tập. Lập hàm trong Julia, nhận các tham số là các vec-tơ (đều có độ dài $N$) là A, B, C, R và giải ra vec-tơ X. Thuật toán Thomas được liệt kê dưới dạng mã lệnh dưới đây. Sau đó, lập chương trình để giải bài toán vi phân trên với $q = 10^4, \alpha=-3, \beta=2, \gamma=100$.

for i = 2:N
    A[i] /= B[i]
    B[i] -= A[i]*C[i-1]
end
for i = 2:N
    R[i] -= A[i]*R[i-1]
end
X[N] = R[N] / B[n]
for i = N-1:-1:1
    X[i] = (R[i] - C[i] * X[i+1]) / B[i]
end

Vấn đề sai phân tại biên. Ta đã thấy tại các vị trí biên không thể sai phân trung tâm để đạt độ chính xác cấp hai được. Thay vào đó, cần phải xây dựng công thức mới biểu diễn sai phân tại biên mà vẫn giữ được độ chính xác cấp 2. Chẳng hạn, tại biên trái, ta có thể chọn ba điểm 1, 2, 3 để lập sai phân tại biên [@Anderson1995]. Có thể coi một đường cong $f(x) = ax^2 + bx + c$ đi qua ba điểm này với vị trí tương ứng $x_1 = 0, x_2 = \Delta x, x_3 = 2\Delta x$. Có thể rút ra $f_1 = c$ cùng 2 PT: $\{ f_2 = a(\Delta x)^2 + b\Delta x + c; f_3 = a(2\Delta x)^2 + b.2\Delta x + c \}$. Hai PT sau này được biến đổi để rút ra: $4f_2 - f_3 = 2b\Delta x + 3c \Rightarrow 4f_2 - f_3 = 2b\Delta x + 3f_1$. Mặt khác, ta cần tính $f' = 2ax + b$; $f'_1 = b$. Vậy $f'_1 = (-3f_1 + 4f_2 - f_3)/(2\Delta x)$.

Bài tập. Khai triển Taylor để thấy đạo hàm trên được tính với độ chính xác cấp 2.

Lưới điểm không đều

Với mục đích cụ thể, người ta có thể dùng lưới với khoảng chia $\Delta x$ không đều. Khi đó, kí hiệu $\Delta x_{j+} = x_{j+1} - x_j$. Với lưới điểm không đều, sai phân trung tâm chỉ có độ chính xác cấp 1 như sai phân tiến hoặc lùi.

Bài tập. Làm thế nào để điều chỉnh công thức sai phân chỉ sử dụng giá trị hai điểm $j-1$ và $j$ nhưng vẫn đạt được độ chính xác cấp 2?

Đạo hàm trên lưới điểm không đều có thể được tính theo cách đề xuất bởi Chapra và Canale [@Chapra2005]: $\frac{df}{dx} \approx f_{j-1} \frac{2x - x_j - x_{j+1}} {\Delta x_{-} (\Delta x_{-} + \Delta x_{+} )} - f_{j} \frac{2x - x_{j-1} - x_{j+1}} {\Delta x_{-} \Delta x_{+}} + f_{j+1} \frac{2x - x_{j-1} - x_{j}} {\Delta x_{+} (\Delta x_{-} + \Delta x_{+} )}$

Ngoài sai phân hữu hạn còn gì nữa?

Ta đã giải PVT bằng cách số trị—tức là xấp xỉ nghiệm của PV—qua việc tính tại nhiều thời điểm cách đều nhau. Trong việc giải PVR có hai vấn đề: biểu diễn nghiệm (hàm) số liên tục dưới dạng một tập hợp các điểm số liệu (véc-tơ), và việc tính đạo hàm tại các điểm đó. Có hai phương châm cơ bản là: nhóm phương pháp điểm lưới và nhóm phương pháp khai triển chuỗi. Ở lời giải Euler, ta sử dụng phương pháp sai phân hữu hạn (FDM) thuộc nhóm phương pháp điểm lưới. Trong sai phân hữu hạn, chỉ có thông tin tại các điểm lưới chứ không có giả thiết gì về giá trị xấp xỉ giữa các nút lưới này.

Một dạng phương pháp quan trọng là thể tích hữu hạn (FVM), trong đó có giả thiết về cấu trúc của nghiệm xấp xỉ giữa các điểm nút lưới này. Trong phương pháp FVM, giá trị tại nút lưới $f_i$ biểu diễn trị trung bình của hàm $f(x)$ trên đoạn, hay ô lưới, $[(i - \frac{1}{2})\Delta x, (i + \frac{1}{2})\Delta x]$. FVM rất hợp sử dụng để xấp xỉ các hàm có điểm gián đoạn. Nếu hàm trơn thì hai phương pháp FDM và FVM cho kết quả như nhau [@Durran2010].

Hình. Xấp xỉ FVM cho một hàm tuần hoàn, dùng đến (a) các hàm bậc thang cho từng đoạn, (b) các hàm tuyến tính cho từng đoạn [@Durran1999].

Về việc tính đạo hàm tại nút lưới, nếu như FDM sử dụng một trong các cách sai phân tiến / lùi / trung tâm, thì FVM quy định giá trị đạo hàm qua cấu trúc được giả sử cho nghiệm xấp xỉ trong các ô tính. Nói riêng, trong FVM, thường tính thông lượng (flux) xuyên qua cạnh các ô tính thay vì đạo hàm; nhưng để tính được các thông lượng này thì phải giả sử cấu trúc nghiệm trong từng ô lưới. Giả sử này không thể chỉ là nội suy tuyến tính giữa các điểm nút lưới; bởi nếu vậy thì giá trị tại nút lưới sẽ không thể bằng trung bình giá trị trong từng ô lưới bao lấy nó được. Phương pháp FVM sẽ được xét đến cùng với các PVR.

Trong nhóm phương pháp khai triển chuỗi, hàm chưa biết được xấp xỉ bởi một tổ hợp tuyến tính từ tập hợp các hàm khai triển liên tục, còn bộ số liệu mô tả hàm xấp xỉ là tập hợp các hệ số khai triển. Các đạo hàm được tính giải tích qua việc vi phân các hàm khai triển. Khi các hàm khai triển này là tập trực giao thì phương pháp được gọi là phương pháp phổ. Một ví dụ là chuỗi Fourier, chẳng hạn xấp xỉ 5 điểm bởi $a_1 + a_2 \cos x + a_3 \sin x + a_4 \cos 2x + a_5 \sin 2x$. Nếu các hàm khai triển chỉ khác 0 trong một khu nhỏ của cả miền tính, thì kĩ thuật khai triển chuỗi được gọi là phương pháp phần tử hữu hạn (FEM). Các hàm khai triển khi đó có thể được chọn là hàm tuyến tính từng đoạn. Mỗi hàm nhận giá trị = 1 tại một điểm nút nhất định và = 0 với tất cả điểm nút khác. (Giữa các nút, giá trị của hàm được xác định bằng nội suy.) Với dạng này, FEM thường cho kết quả giống với phương pháp điểm lưới.

Phương trình vi phân riêng

Trong các bài toán thuỷ động lực, ta thường gặp các PVR tuyến tính bậc 2 với ít nhất hai biến độc lập. Chẳng hạn:

• PT sóng phi tuyến mô tả quá trình diễn biến cuả mặt front shock do tính phi tuyến của sóng, áp dụng trong khí động lực học. $\frac{\partial f}{\partial t} + f \frac{\partial f}{\partial x} = 0$

• PT Burgers mô tả quá trình chuyển tải phi tuyến và khuếch tán, được áp dụng cho các sóng âm, sóng nước. Với $\nu = 0$ có thể áp dụng cho chất lỏng không nhớt, và dạng PT Burgers không nhớt chính là dạng PT sóng phi tuyến nêu trên. $\frac{\partial f}{\partial t} + f \frac{\partial f}{\partial x} = \nu \frac{\partial^2 f}{\partial x^2}$

• Hệ phương trình St Venant cho dòng chảy 1 chiều có mặt thoáng. Có hai dạng trình bày phương trình này: dạng bảo toàn trong đó các biến lượng là diện tích mặt cắt, $A$, và lưu lượng, $Q$. Dạng không bảo toàn được đề cập dưới đây có các biến lượng độ sâu nước, $h$, và vận tốc dòng chảy, $u$; dạng này có hình thức đơn giản hơn nhưng nhược điểm là không áp dụng được khi dòng chảy có đột biến, ví dụ nước nhảy. Hệ PT St-Venant còn được gọi là hệ phương trình nước nông, hệ PT sóng dài, và (bằng cách mở rộng ra 2 chiều) được áp dụng cho thuỷ động lực ven bờ như động thái dòng triều và dao động nước trong bể cảng. $\frac{\partial h}{\partial t} + \frac{\partial (uh)}{\partial x} = q$ $\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + g \frac{\partial h}{\partial x} - g S_b + Bh \frac{\tau_b}{\rho h} = q u_q$

• Hệ phương trình Navier-Stokes cho dòng chảy 3 chiều của chất lỏng nhớt nén được; có thể coi là phương trình tổng quát chi phối động thái của những chất lưu Newton: $\frac{\partial (\rho f)}{\partial t} + \frac{\partial (\rho u_i f)}{\partial x_i} = \frac{\partial}{\partial x_i} \left( D \frac{\partial f}{\partial x_i} \right) + q_f$ trong đó các chỉ số ${}_i$ thể hiện tổng các số hạng theo các chiều không gian $\{x_i\} = x, y, z$ còn $f$ là đại lượng được bảo toàn. Tuỳ vào các điều kiện cụ thể mà hệ PT Navier-Stokes được giản hoá để mô tả dòng chảy không nén được, dòng không nhớt, dòng chảy thế, dòng chảy rất chậm (creeping flow, có độ nhớt lớn); và dùng các phép xấp xỉ như Boussinesq và lớp biên.

Phương trình phi thứ nguyên

Một kĩ thuật đôi khi được dùng trong mô phỏng thuỷ động lực là làm phi thứ nguyên phương trình; nhằm mục đích biểu thị độ lớn tương đối giữa các số hạng thành phần trong phương trình, từ đó có thể giúp ta định hướng giải đúng đắn [@Wood1993]. Chẳng hạn, từ phương tình Navier-Stokes 1 chiều cho dòng không nén được: $\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + g \frac{\partial h}{\partial x} + \frac{1}{\rho} \frac{\partial p}{\partial x} = \nu \frac{\partial^2 u}{\partial x^2}$ với $p$ là áp suất và $\nu$ là độ nhớt động học. Phương trình có thể được phi thứ nguyên hoá bằng đưa vào các đại lượng tham chiếu về chiều dài ($L$), vận tốc ($U$); khi đó ta thiết lập được các biến phi thứ nguyên: $x^* = x/L \quad h^* = h/L \quad u^* = u/L, \quad t^* = t/(L/U) \quad p^* = p/(\rho U^2)$

Thay các biến phi thứ nguyên vào PT N-S 1 chiều, sau đó nhân hai vế với $L/U^2$ ta thu được: $\frac{\partial u^*}{\partial t^*} + u^* \frac{\partial u^*}{\partial x^*} + \frac{1}{\mathrm{Fr}^2} \frac{\partial h^*}{\partial x^*} + \frac{\partial p^*}{\partial x^*} = \frac{1}{\mathrm{Re}} \frac{\partial^2 u^*}{\partial x^{*2}}$ Ở đây đã đưa vào số Froude Fr = $U / \sqrt{gL}$ và số Reynolds Re $= UL/\nu$. Trong dòng chảy tự nhiên, số Re thường lớn (dòng chảy rối) và số hạng đạo hàm bậc hai ($\frac{1}{\mathrm{Re}} (\partial^2 u^* / \partial x^{*2})$) thường được bỏ qua. Còn độ lớn của số Fr cho biết trạng thái của dòng chảy là êm, xiết, hay phân giới.

Việc chọn các đại lượng tham chiếu trong các dòng chảy đơn giản thì rất hiển nhiên: $U$ được chọn là vận tốc dòng đều, $L$ là chiều sâu dòng đều hoặc chiều dài thuỷ vực. Sau khi đã biến đổi phương trình thì hoàn toàn có thể áp dụng sai phân cho phương trình phi thứ nguyên này. Mặc dù vậy, với các dạng dòng chảy phức tạp, có nhiều quá trình vật lý diễn ra dẫn đến số biến phi thứ nguyên tăng lên khiến cho tính toán khó khăn hơn. Bởi vậy, trong các chương trình phần mềm không phải lúc nào cũng áp dụng cách làm này và phương pháp chỉ được nêu ra đây mà thôi.

Phân loại PVR

Với hai biến $x, t$, và đại lượng $u$ là hàm (phụ thuộc $x, t$) cần xác định, dạng tổng quát của PVR như sau: $A \frac{\partial^2 u}{\partial x^2} + B \frac{\partial^2 u}{\partial x \partial t} + C \frac{\partial^2 u}{\partial t^2} = D \frac{\partial u}{\partial x} + E \frac{\partial u}{\partial t} + F u + G$

trong đó $A, B, \dots, G$ có thể là các hàm của $x, t$, hoặc cả hai, hoặc chỉ là hằng số. Định thức của PT vi phân bậc hai nêu trên là: $\Delta = B^2 - 4AC$.

Tuỳ theo dấu của $\Delta$, PT vi phân (tổng quát) được chia thành:

1. PVR elip, nếu $\Delta < 0$ với nghiệm trơn, dễ giải;
2. PVR hyperbol ($\Delta > 0$), với nghiệm có thể chứa điểm gián đoạn, thường khó giải;
3. PVR parabol ($\Delta = 0$), có thể mang đặc tính của cả 2 loại PVR trên.

Mỗi loại PVR nêu trên mô tả một loại hiện tượng vật lý trong thuỷ lực. Dạng elip biểu thị diễn biến ổn định, không phụ thuộc thời gian (như dòng nước ngầm ổn định, các đường backwater trong thuỷ lực kênh hở). Dạng hyperbol biểu thị các quá trình vật lý phuộc thuộc thời gian, có tính bảo toàn (như quá trình chuyển tải / đối lưu), không ổn định và không có xu hướng tiến tới trạng thái ổn định. Cuối cùng, dạng parabol biểu diễn các quá trình mang tính tiêu tán, phụ thuộc vào thời gian (chẳng hạn khuếch tán) vốn hướng đến một trạng thái ổn định.

Khác với PVT, các PVR thường dùng lược đồ bậc thấp hơn. Có hai lý do: thứ nhất sai số do vi phân thời gian không phải là nguồn duy nhất gây sai số đến nghiệm. Thứ hai là để giảm công thực hiện tính toán.

Nghiệm của PVR có thể giải được nếu xác lập đầy đủ điều kiện ban đầu và điều kiện biên. Có ba loại điều kiện biên: Dirichlet (ấn định giá trị tại biên), Neumann (ấn định gradient tại biên) và Robin (loại hỗn hợp).

Sự ảnh hưởng của điều kiện đầu, điều kiện biên đến nghiệm của bài toán được thể hiện như trên Hình characteristics. Trong bài toán PVR hyperbol trong 1 chiều không gian, có 2 đường đặc tính; với bài toán parabol thì 2 đường đặc tính này suy biến thành một; còn bài toán elip thì không có đường đặc tính nào (toàn bộ miền tính đều ảnh hưởng đến nghiệm).

Hình … Các miền ảnh hưởng đến nghiệm trong PT vi phân [@Versteeg2007]: (a) hyperbol, (b) parabol, và ( c ) elip

Một vấn đề quan trọng trong việc giải PVR là sai số. PVR dạng hyperbol dễ mắc sai số nhất vì các điều kiện biên sẽ gây nên sai số truyền đi vào bên trong miền tính.

Phương trình khuếch tán

Phương trình khuếch tán trong không gian 1 chiều có dạng: $\frac{\partial f}{\partial t} - D\frac{\partial^2 f}{\partial x^2} = 0$ Đây là phương trình dạng parabol. Hằng số $D>0$ là hệ số khuếch tán. Có thể sai phân hoá theo kiểu FTCS như sau: $\frac{f_j^{n+1} - f_j^n}{\Delta t} - D\frac{f_{j-1}^n - 2f_j^n + f_{j+1}^n}{\Delta x^2} = 0$ Từ đó rút ra công thức dự tính giá trị tại thời điểm tương lai: $f_j^{n+1} = f_j^n + \frac{D \Delta t}{\Delta x^2} (f_{j-1}^n - 2f_j^n + f_{j+1}^n)$

Bài tập. Bài toán “khuếch tán đường bờ”: Giả sử miền tính toán đuọc chọn là bờ biển dài 1 km; người ta thực hiện nuôi bãi nhân tạo trong phạm vi dọc theo 200 m ở khoảng chính giữa đường bờ trên. Bề rộng của vùng đổ cát là 5 m. Tác động của sóng khiến cho bùn cát bị khuếch tán và sẽ dàn sang hai bên. Cho hệ số khuếch tán đường bờ $D = 4.2 \times 10^6$ m$^2$/năm, hãy biểu diễn vị trí của đường bờ.

Có thể nhận thấy rằng, trong bài tập trên, $\Delta t$ không thể chọn nhỏ tuỳ ý vì nghiệm sẽ mất ổn định. Biểu thức ràng buộc giá trị của $\Delta t$ sẽ được làm sáng tỏ thông qua phép phân tích ổn định. Cũng như trong mục 2.3, chỉ số để đánh giá ổn định là hệ số khuếch đại $G = f^{n+1} / f^n$, song ở đây có khác vì có mặt những giá trị ở các điểm lân cận, $f_{j-1}$ và $f_{j+1}$. Để giải quyết điều này, phương pháp Von Neumann được dùng. Coi rằng sai số gắn với $f_j^n$ là $A^n e^{-ikj\Delta x}$ trong đó $i = \sqrt{-1}$ là số ảo còn $k$ là một số sóng của nhiễu động sai số. Khi đó, từ PT [rời rạc hoá] ta có: $\begin{aligned} A^{n+1} e^{-ikj\Delta x} = A^n e^{-ikj\Delta x} + \frac{D \Delta t}{\Delta x^2} \left(A^n e^{-ik(j-1)\Delta x} - 2A^n e^{-ikj\Delta x} + \right. \\ \left. + A^n e^{-ik(j+1)\Delta x} \right)\end{aligned}$ với lưu ý $e^{-ik(j+1)\Delta x} = e^{-ikj\Delta x} e^{-ik\Delta x}$, có thể rút gọn số hạng chung $e^{-ikj\Delta x}$, và

$A^{n+1} / A^n = 1 + \frac{D \Delta t}{\Delta x^2} \left( e^{ik\Delta x} - 2 + e^{-ik\Delta x} \right)$

với $\cos y = (e^{iy} + e^{-iy})/2$, biểu thức được viết thành: $G = A^{n+1} / A^n = 1 + \frac{2D \Delta t}{\Delta x^2} \left( \cos {k\Delta x} - 1 \right)$

Phép sai phân này chỉ ổn định khi $|G| < 1$ hay $-1 < G < 1$. Vì cos(…) < 1 nên $G<1$ luôn đúng. Chỉ cần xét $G > -1$ hay $-\frac{2D \Delta t}{\Delta x^2} \left( \cos {k\Delta x} - 1 \right) < 2 \Leftrightarrow \frac{2D \Delta t}{\Delta x^2} \left( 2\sin^2 \frac{k\Delta x}{2} \right) < 2$. Rõ ràng vì $\sin^2 (...) \leqslant 1$ nên $D\Delta t / \Delta x^2 \leqslant \frac{1}{2}$.

Phương pháp Von Neumann đã dùng ở trên áp dụng được với các bài toán điều kiện đầu dạng tuần hoàn và hệ số khuếch tán $D$ là hằng số. Mặc dù vậy, trường hợp $D$ biến đổi trên miền tính thì vẫn có thể dùng phương pháp này để phân tích ổn định cục bộ.

Cũng từ ví dụ về sơ đồ FTCS cho bài toán khuếch tán, có thể thấy rõ rằng một trị số $f_{j}^{n+1}$ chỉ phụ thuộc vào ba trị số ở bước thời gian trước đó, là $f_{j-1}^n, f_j^n, f_{j+1}^n$ và do vậy phụ thuộc vào 5 trị số ở bước trước nữa: $f_{j-2}^{n-1}, f_{j-1}^{n-1}, f_j^{n-1}, f_{j+1}^{n-1}, f_{j+2}^{n-1}$. Cứ như vậy, miền ảnh hưởng nghiệm số trị biểu diễn trên mặt phẳng $x~t$ là một tam giác cân (một cách biểu diễn tương tự, nhưng cho bài toán khác, có thể thấy trên Hình [delta-t]). Một điều phi lý của sơ đồ sai phân là lẽ ra mọi điểm $f_{*}^n$ cùng lớp thời gian phải có ảnh hưởng đến $f_j^n$ mới đúng bản chất hiện tượng khuếch tán trong thực tế. Điều này được khắc phục bằng một sờ đồ sai phân ẩn (mục 3.6).

Quá trình giải kiểu “duyệt binh”. Không chỉ có FTCS mà mọi phương pháp khác để giải PT vi phân theo thời gian đều thực hiện kiểu “duyệt binh” (marching). Khi biểu diễn miền tính toán là lưới có trục hoành $x$ và trục tung $t$, mỗi đường lưới nằm ngang là tập hợp nghiệm số ở từng lớp thời gian; chúng phải được giải xong (dù giải đồng thời hay lần lượt) trước khi chuyển sang lớp thời gian tiếp theo. Quá trình chuyển động theo hàng này gợi đến hình ảnh duyệt binh, và tương phản với quá trình “căn chỉnh” khi giải PT elip 2 chiều (mục 3.8), khi đó trị số mỗi điểm nút trong hệ phải được liên tục điều chỉnh theo các trị số nút kề bên, sao cho giá trị tổng thể của toàn miền tuân theo PT cơ bản (kiểu Poisson hoặc Laplace).

Phương trình chuyển tải 1 chiều

Phương trình chuyển tải một chiều (1-D convection) rất thông dụng trong các bài toán thuỷ động lực học như một phép kiểm tra chất lượng của các phương pháp sai phân. PT có dạng:

$\frac{\partial f}{\partial t} + C\frac{\partial f}{\partial x} = 0$ trong đó $C$ = const là tốc độ chuyển tải; tuỳ theo dấu của $C$ mà hiện tượng chuyển tải theo hoặc ngược chiều dương trục $x$.¹

Sơ đồ sai phân FTCS cho PT [1D-conv] là: $\frac{f_j^{n+1} - f_j^n}{\Delta t} + C\frac{f_{j+1}^n - f_{j-1}^n}{2\Delta x} = 0$

Khác với trường hợp PT khuếch tán, dạng sai phân FTCS của PT chuyển tải luôn bất ổn định vô điều kiện, nghĩa là dù thu ngắn bước thời gian đi bao nhiêu nữa, nghiệm số trị vẫn không thể hội tụ.

Bài tập. Bằng phương pháp phân tích ổn định Von Neumann, hãy cho thấy rằng sơ đồ FTCS cho bài toán chuyển tải 1 chiều là bất ổn định không điều kiện (nghĩa là $|G| > 1$ với mọi $\Delta x, \Delta t$). (Gợi ý: một số biến đổi toán học hữu ích: $e^{iy} - e^{-iy} = 2i \sin y$ và $|a + ib| = \sqrt{a^2 + b^2}$.)

Để khắc phục nhược điểm của FTCS, người ta đã xây dựng những sơ đồ sai phân tốt hơn để giải PT chuyển tải một chiều: $\frac{f_j^{n+1} - \frac{1}{2}(f_{j-1}^n + f_{j+1}^n)}{\Delta t} + C\frac{f_{j+1}^n - f_{j-1}^n}{2\Delta x} = 0 \\ \tag{sơ đồ Lax}$ $\frac{f_j^{n+1} - f_j^n}{\Delta t} + C\frac{f_{j+1}^n - f_{j-1}^n}{2\Delta x} = C^2 \frac{\Delta t}{2\Delta x^2} (f_{j-1}^n - 2f_j^n + f_{j+1}^n) \\ \tag{sơ đồ Lax-Wendroff}$
$\frac{f_j^{n+1} - f_j^n}{\Delta t} + \begin{cases} C\frac{f_{j}^n - f_{j-1}^n}{\Delta x} = 0 & C > 0 \\ C\frac{f_{j+1}^n - f_{j}^n}{\Delta x} = 0 & C < 0 \\ \end{cases} \\ \tag{sơ đồ ngược dòng}$

Sơ đồ Lax không sử dụng giá trị “tại chỗ” để tính sai phân tiến theo thời gian, mà dùng trung bình hai giá trị lân cận, $\frac{1}{2}(f_{j-1}^n + f_{j+1}^n)$. Trong sơ đồ Lax-Wendroff, số hạng được bổ sung thêm ở vế phải có dáng vẻ của thành phần “khuếch tán”. Việc đưa vào số hạng này để bù cho số hạng bậc hai – số hạng đầu trong sai số cắt cụt với sai phân thời gian ở vế trái. Thật vậy, $\frac{\Delta t}{2} \frac{\partial^2 f}{\partial t^2} = \frac{\Delta t}{2} \frac{\partial}{\partial t} \left( -C \frac{\partial f}{\partial x} \right) = \frac{\Delta t}{2} \frac{\partial}{\partial x} \left( -C \frac{\partial f}{\partial t} \right) = \frac{C^2 \Delta t}{2} \frac{\partial^2 f}{\partial x^2}$

Sơ đồ ngược dòng sử dụng sai phân tiến/lùi trong không gian và chỉ có độ chính xác bậc nhất, nhưng một chi tiết quan trọng của kĩ thuật ngược dòng là một điểm nào đó trong đoạn được sai phân sẽ mang thông tin với tốc độ truyền sóng đến điểm cần dự tính trong tương lai.

Bài tập. Sự lan truyền hình dạng sóng. Lập chương trình thực hiện sai phân FTCS với bài toán chuyển tải 1 chiều có $C = 0.1$ trên miền dài $N\Delta x = 10$. Điều kiện đầu là hàm “mũ phù thủy”: $\begin{aligned} f_j = \begin{cases} 0.1(j-20) & 20 \leqslant j \leqslant 25 \\ 0.1(30-j) & 25 < j \leqslant 30 \\ 0 & \forall j \ne \end{cases}\end{aligned}$ Giảm nhỏ $\Delta t$ để thấy rằng không thể triệt tiêu được những nhiễu động. Sửa đổi sơ đồ FTCS bằng một trong các phương pháp cải thiện (Lax, Lax-Wendroff, ngược dòng). Sau đó, hãy đổi dấu của $C$ và quan sát động thái kết quả.

Các hàm gián đoạn kiểu “bậc” thường được sử dụng để kiểm chứng chất lượng của lược đồ sai phân cho bài toán chuyển tải 1 chiều.

Hình. Nghiệm số trị cho PT chuyển tải ($C$ = 0.5) tính theo các phương pháp khác nhau: FTCS (phía trái trên), ngược dòng bậc nhất (phải trên), Lax–Wendroff (trái dưới) và Crank–Nicolson (phải dưới). [@Watanabe2013]

Hình. Ảnh hưởng của độ lớn bước thời gian tới tương quan giữa miền phụ thuộc của sơ đồ ngược dòng (các điểm $\circ$) và của PT chuyển tải 1 chiều (đường đứt nét). (a) $\Delta t$ ổn định, (b) không ổn định.[@Durran2010]

Bài tập. Bằng phương pháp Von Neumann, phân tích ổn định với sơ đồ ngược dòng. Biểu thức của điều kiện ổn định có thể được đưa về dạng $u \Delta t / \Delta x \leqslant 1$. Đại lượng không thứ nguyên vế trái là số Courant.

Bất đẳng thức trong bài trên là điều kiện Courant; với hàm ý tốc độ dòng chảy, $U$, phải nhỏ hơn tốc độ truyền thông tin trên lưới tính, $\Delta x / \Delta t$. Cụ thể, như trong hình [delta-t], với cùng tốc độ $U$ thì miền ảnh hưởng đến nghiệm số PT sẽ có dạng tam giác tạo bởi đường đứt nét và hình chiếu của nó. Với $\Delta t$ nhỏ, các điểm nút mà giá trị được dùng cho sai phân che kín hết miền ảnh hưởng.

Miền ảnh hưởng đến nghiệm số trị được biểu diễn dưới dạng một vùng các điểm $\circ$ kề nhau trên lưới tính. Dĩ nhiên lược đồ sai phân hiện sẽ tạo nên vùng ảnh hưởng tam giác có đỉnh nhọn hướng lên. Với lược đồ cho bài toán khuếch tán, miền tam giác cân sẽ trải đều về hai phía trục $x$.

Phương trình chuyển tải-khuếch tán

Đây là dạng kết hợp của hai PVR vừa đề cập đến. Một ví dụ ứng dụng là sự lan truyền chất ô nhiễm dọc theo đoạn sông. Nồng độ chất, $c$, biến đổi theo cả phương dọc sông, $x$, và thời gian, $t$: $\frac{\partial c}{\partial t} + U\frac{\partial c}{\partial x} + D\frac{\partial^2 c}{\partial x^2} + Kc = 0$

Dòng chảy mang chất ô nhiễm đi có tốc độ chuyển tải $U$ = const (cỡ 1 m/s). Sự khuếch tán rối làm lan rộng chất ô nhiễm với hệ số khuếch tán $D$ = const (cỡ 1 m$^2$/s). Tải lượng chất ô nhiễm là $Kc$ với $K$ = const phân bố đều trên đoạn sông. Số hạng cuối (không chứa đạo hàm của $c$), sau khi chuyển vế, được gọi là số hạng nguồn. Cách sai phân số hạng nguồn nằm chính ở trong nội dung của bài toán điều kiện đầu, mục 2.3.

Bài tập. Bài toán dầu thải [@Chapra1997]: Một dòng sông có đặc trưng sau: bề rộng: 60 m, chiều sâu: 1 m, vận tốc dòng chảy trung bình: 2400 m/giờ, hệ số khuếch tán 150000 m$^2$/giờ. Giả sử có lượng dầu là $m =$ 5 kg thải ra tại vị trí 500 m hạ lưu khu dân cư. Hãy viết phương trình sai phân và lập chương trình tính ra phân bố của nồng độ dầu dọc theo sông, cụ thể tại các thời điểm 0.2 giờ, 0.4 giờ, 0.6 giờ sau khi xả thải. Chọn các bước điểm nút lưới và bước thời gian thích hợp. Nghiệm giải tích có dạng $\displaystyle{ c(x,t) = \frac{m/(Bh)} {2\sqrt{\pi D t } } \exp\left( \frac{(x-Ut)^2}{4Dt} \right) }.$

Gợi ý: Dùng điều kiện Pe < 2 (xem thêm mục dưới) để xác định chiều dài đoạn sông phù hợp. Sau đó dùng điều kiện ổn định khuếch tán, $D\Delta t / \Delta x^2 < \frac{1}{2}$ để xác định bước thời gian phù hợp. Để xác lập được điều kiện đầu, cần ấn định nồng độ cho 1 đoạn sông đúng chỗ xả thải, $c_0 = m/Bh\Delta x$ (nồng độ ban đầu của các đoạn khác thì bằng 0). $\blacksquare$

Quan điểm trên về việc coi $\Delta x$ là chiều dài đoạn thay vì khoảng cách giữa các điểm nút là điều thường gặp trong lĩnh vực mô hình hoá chất lượng nước; một dòng sông được chia thành nhiều “thể tích kiểm soát”, như phương pháp thể tích hữu hạn (FVM) ở mục 4. Ưu điểm của cách làm này là trực quan, có thể căn cứ vào bề rộng sông để phân đoạn sao cho mỗi đoạn có dạng hình lăng trụ ($Bh \sim$ const), dù phải chấp nhận việc tâm của các thể tích kiểm soát này sẽ không cách đều nhau.

Giải thích điều kiện về số Péclet. Số Péclet trong bài toán chuyển tải-khuếch tán được định nghĩa là Pe = $C \Delta x / D$. Nó thể hiện độ mạnh tương đối giữa quá trình chuyển tải và khuếch tán.² Khi sử dụng phương pháp ngược dòng cho số hạng chuyển tải cùng với sai phân trung tâm cho số hạng khuếch tán, phép ngược dòng chỉ đảm bảo độ chính xác cấp 1 với số hạng lớn nhất trong sai số cắt cụt là (xem lại PT [TE]). Điều đó dẫn đến PT: $\displaystyle{ \frac{\partial f}{\partial t} + C\frac{\partial f}{\partial x} = D \left( 1 + \frac{\mathrm{Pe}}{2} \right) \frac{\partial^2 f}{\partial x^2} }$, nghĩa là hệ số khuếch tán đã tăng thêm một lượng bằng Pe $\times D/2$, đây gọi là thành phần khuếch tán giả (false diffusion) hay khuếch tán số trị. Đây cũng là nguyên nhân của biểu hiện nghiệm số trị tính theo phương pháp ngược dòng sẽ bị bẹt đi theo thời gian.

Bài tập. Biến đổi để thu được công thức trên. Muốn giữ cho lượng khuếch tán số trị nhỏ hơn khuếch tán thực thì Pe phải thỏa mãn điều kiện gì?

Sơ đồ hiện và ẩn

Các sơ đồ thiết lập đến giờ đều là sơ đồ hiện (explicit scheme), theo nghĩa trong mỗi phương trình, đại lượng duy nhất ở bước thời gian kế tiếp cần dự tính chỉ là $f_j^n$. Với đặc điểm này, sơ đồ hiện có ưu điểm là ta có thể giải ngay mỗi phương trình để tìm ra nghiệm. Ngược lại, trong sơ đồ sai phân ẩn (implicit scheme), cần kết hợp tất cả các phương trình ở nút lưới mới có thể tìm được các giá trị. Điều đó dẫn đến giải hệ PT tuyến tính cho các $f_j^{n+1}$ với ma trận hệ số thường là ma trận thưa với các phần tử khác 0 nằm gần đường chéo chính.

Chẳng hạn với phương trình khuếch tán, một sơ đồ ẩn có dạng: $\frac{f_j^{n+1} - f_j^n}{\Delta t} - D\frac{f_{j-1}^{n+1} - 2f_j^{n+1} + f_{j+1}^{n+1}}{\Delta x^2} = 0$ Riêng dạng này còn được gọi là sơ đồ toàn ẩn (fully implicit scheme), với tất cả các đạo hàm theo không gian đều được sai phân tại lớp thời gian kế tiếp ($n+1$). Sơ đồ này gọi là BTCS, đối lập với sơ đồ FTCS trong sơ đồ hiện.

Dạng tổng quát hơn với sai phân của PT khuếch tán là: $\frac{f_j^{n+1} - f_j^n}{\Delta t} - \frac{D}{\Delta x^2} \left[ (1 - \theta)(f_{j-1}^{n} - 2f_j^{n} + f_{j+1}^{n}) + \theta(f_{j-1}^{n+1} - 2f_j^{n+1} + f_{j+1}^{n+1}) \right]= 0$ Hệ số $\theta$ ($0 \leqslant \theta \leqslant 1$) thể hiện mức độ “ẩn” của sơ đồ. $\theta=0$ ứng với sơ đồ hiện, $0<\theta<1$: sơ đồ ẩn, và riêng $\theta=1$: toàn ẩn. Đặc biệt, với $\theta=0.5$, trọng số của các giá trị ở lớp thời gian $n$ và $n+1$ là như nhau. Sơ đồ như vậy có tên là Crank-Nicolson, còn gọi là CTCS.

Phương pháp tách thời gian MacCormack

Phương pháp đơn giản này cho độ chính xác cao (bậc 2 cả trong không gian lẫn thời gian) [@Lin2008], nhưng cũng như các lược đồ sai phân trung tâm trong không gian (Lax-Wendroff và Crank-Nicolson), một nhược điểm là có thể gây xuất hiện nhiễu động giả (spurious).

Phương pháp được thực hiện qua hai bước (“dự đoán - hiệu chỉnh”). Ở bước “đoán”, một giá trị sơ bộ, $\tilde{f}$, được tính ra. Giá trị sơ bộ này được chữa lại ở bước hiệu chỉnh, khi đó đạo hàm không gian được lấy trung bình cộng của các sai phân tiến ở lớp thời gian $t^n$ và sai phân lùi ở $t^{n+1}$: $\begin{aligned} \tilde{f}_j = f_j^n - C\frac{\Delta t}{\Delta x} (f_{j+1}^n - f_j^n) \\ f_j^{n+1} = f_j^n - C\frac{\Delta t}{2\Delta x} (f_{j+1}^n - f_j^n + \tilde{f}_j - \tilde{f}_{j-1}) \end{aligned}$

Một dạng của phương pháp tách này là ADI, sẽ được trình bày ở mục 3.9.

Sai phân theo không gian 2 chiều

Từ các sơ đồ sai phân đã xây dựng trong các mục trên có thể dễ dàng mở rộng sang không gian 2 chiều.

Dạng đơn giản nhất là phương trình elip: $K_x \frac{\partial^2 \phi}{\partial x^2} + K_y \frac{\partial^2 \phi}{\partial y^2} = q(x,y)$

Đây là phương trình Poisson, nó có thể dùng để mô tả phân bố nhiệt ($\phi=T$) hoặc cột nước trong tầng nước ngầm có áp ($\phi=H$). Trong trường hợp thứ hai, ý nghĩa của $K_x$ và $K_y$ là độ dẫn thuỷ lực (= const, đơn vị m/s) còn $q$ (m/s) là nguồn hoặc điểm rút.

• Khi $q=0$, ta có phương trình Laplace, $\nabla^2 \phi = 0$.
• Khi $K_x, K_y \neq$ const, mà biến đổi theo $(x,y)$, thì cần viết lại PT: $\frac{\partial}{\partial x} \left( K_x \frac{\partial \phi}{\partial x} \right) + \frac{\partial}{\partial y} \left( K_y \frac{\partial \phi}{\partial y} \right) = q(x,y).$

Một khó khăn khi sai phân trong bài toán 2 chiều là trong trường hợp biên có dạng hình học phức tạp, ví dụ như trường hợp biên dạng Dirichlet (cho số), xem cuốn sách của Chapra và Canale [@Chapra2005], [công thức (29.23-24)]; tuy nhiên các sai phân như vậy cũng chỉ đạt được độ chính xác cấp 1. Trường hợp biên dạng Neumann (cho đạo hàm) thì còn khó hơn nữa.

Một ví dụ khác về PVR hai chiều; xét phương trình sóng dài và biên độ nhỏ,³ ở dạng đơn giản (bỏ qua ma sát và số hạng chuyển tải): $\begin{gathered} \frac{\partial \eta}{\partial t} + h \frac{\partial u}{\partial x} + h \frac{\partial v}{\partial y}= 0 \\ \frac{\partial u}{\partial t} = -g \frac{\partial \eta}{\partial x} \\ \frac{\partial v}{\partial t} = -g \frac{\partial \eta}{\partial y} \end{gathered}$

trong đó $\eta$ là mực nước, $h$ là độ sâu nước tĩnh, $u,v$ là các thành phần vận tốc dòng chảy theo phương $x,y$. PT này có thể được đưa về: $\displaystyle{\frac{\partial^2 \eta}{\partial t^2} = C^2 \left( \frac{\partial^2 \eta}{\partial x^2} + \frac{\partial^2 \eta}{\partial y^2} \right) }$, với $C=\sqrt{gh}$ là tốc độ truyền sóng nước sâu. (Bạn hãy tự biến đổi!) Đó là PT dạng gì? elip, hyperbol hay parabol?

Đối với PT Poisson, sau khi sai phân ta sẽ được hệ phương trình, như trường hợp bài toán điều kiện biên (mục 2.4). Tuy nhiên, cái khác ở đây là không gian 2 chiều sẽ khiến cho mỗi phương trình có 5 ẩn (một nút tại chỗ và 4 nút lân cận). Điều đó dẫn đến một hệ PT 5 đường chéo với $N = N_x \times N_y$ ẩn. Ngoài ba đường chéo trung tâm thì hai đường chéo ngoài sẽ cách giữa khoảng $N_x$ hoặc $N_y$ phần tử, tuỳ theo cách lưu các biến ma trận theo các hàng ngang những giá trị $x$ hay theo các cột dọc những giá trị $y$. Không có thuật toán như Thomas để giải nhanh hệ này, và do vậy để tránh giải theo cách chính xác (nhưng lâu), người ta thường dùng cách xấp xỉ Gauss–Seidel (Liebmann) (mục 29.2, Chapra & Canale (2005) [@Chapra2005]). Để đẩy nhanh quá trình lặp, người ta có thể dùng cách kéo dãn quá SOR (successive over-relaxation). Để minh hoạ cho cách làm này, xét trường hợp đơn giản, $K_x = K_y = 1$ và $\Delta x = \Delta y$. Khi đó, phương trình Laplace được sai phân thành: $\begin{aligned} \frac{\phi_{j-1,k} - 2\phi_{j,k} + \phi_{j+1,k}}{\Delta x^2} + \frac{\phi_{j,k-1} - 2\phi_{j,k} + \phi_{j,k+1}}{\Delta y^2} = 0 \\ \Rightarrow \phi_{j-1,k} + \phi_{j+1,k} + \phi_{j,k-1} + \phi_{j,k+1} - 4\phi_{j,k} = 0 \\ \Rightarrow \phi_{j,k} = \frac{1}{4} ( \phi_{j-1,k} + \phi_{j+1,k} + \phi_{j,k-1} + \phi_{j,k+1} )\end{aligned}$

Việc giải lặp được thực hiện từ $k=1\to N_y$ và $j=1\to N_x$. Cuối mỗi vòng lặp, có thể cập nhật giá trị $\phi_{j,k}$ bằng cách dãn quá: $\phi_{j,k} = \lambda \phi_{j,k} + (1-\lambda) \phi_{j,k}^{old}$ trong đó $\phi_{j,k}^{old}$ là giá trị cũ từ lần lặp trước còn $\lambda \in [1,2]$ là một hệ số dãn.

Sơ đồ ẩn luân hướng (ADI) cho bài toán điều kiện đầu

Trong bài toán điều kiện đầu, việc xây dựng sơ đồ sai phân ẩn cũng dẫn đến một ma trận tương tự. Tuy nhiên, trong trường hợp này người ta thường dùng một phương pháp có tên Alternating-Direction Implicit (ADI, ẩn luân hướng). Theo cách này, mỗi bước thời gian sẽ được chia ra thực hiện bằng hai khâu. Ví dụ như với PT khuếch tán 2 chiều, $\partial \phi / \partial t = K_x(\partial^2 \phi / \partial x^2) + K_y (\partial^2 \phi / \partial y^2)$, ở khâu thứ nhất: $\frac{\phi_{j,k}^{n+} - \phi_{j,k}^{n}}{\Delta t / 2} = K_x \left( \frac{\phi_{j-1,k}^n - 2\phi_{j,k}^n + \phi_{j+1,k}^n}{\Delta x^2} \right) + K_y \left( \frac{\phi_{j,k-1}^{n+} - 2\phi_{j,k}^{n+} + \phi_{j,k+1}^{n+}}{\Delta y^2} \right)$

Như vậy, số hạng $\partial^2 \phi / \partial x^2$ được sai phân hiện còn $\partial^2 \phi / \partial y^2$ thì sai phân ẩn. Ta nhận được một hệ PT 3 đường chéo để giải ra các $\phi^{n+}$ ở “nửa bước thời gian”. Những trị số này được sử dụng trong khâu thứ hai, ở đó số hạng $\partial^2 \phi / \partial x^2$ được sai phân ẩn: $\frac{\phi_{j,k}^{n+1} - \phi_{j,k}^{n+}}{\Delta t / 2} = K_x \left( \frac{\phi_{j-1,k}^{n+1} - 2\phi_{j,k}^{n+1} + \phi_{j+1,k}^{n+1}}{\Delta x^2} \right) + K_y \left( \frac{\phi_{j,k-1}^{n+} - 2\phi_{j,k}^{n+} + \phi_{j,k+1}^{n+}}{\Delta y^2} \right)$

và sau khi giải hệ ta sẽ thu được các $\phi^{n+1}$. Phương pháp ADI thuộc nhóm các phương pháp “tách” (splitting method).

Phân tích ổn định đối với sơ đồ hiện giải PT khuếch tán parabol 2 chiều; so sánh điều kiện ổn định với kết quả phân tích bài toán 1 chiều.

Phương pháp thể tích hữu hạn

Phương pháp thể tích hữu hạn (FVM) có thể áp dụng được cho các PT vi phân biểu diễn các dạng định luật bảo toàn. Ưu điểm của nó là có thể dùng được cho miền có dạng hình học bất kì và hỗ trợ cả lưới tính toán cấu trúc lẫn phi cấu trúc. Trong phương pháp FVM:

• Miền tính toán được chia thành các phần tử riêng biệt gọi là các thể tích kiểm soát (control volume)
• PT vi phân sẽ được tích phân để thành dạng một phương trình cân bằng đối với mỗi thể tích kiểm soát; sau đó các tích phân này được xấp xỉ bằng phương pháp số trị.

Đặc điểm quan trọng của FVM là thông lượng được bảo toàn từ ô tính toán này sang ô liền kề. Nghĩa là có sự bảo toàn tại chỗ đối với thông lượng. Điều đó rất quan trọng vì thông lượng rất cần thiết khi tính toán thuỷ động lực học.[@Popescu2014]

Một sơ đồ thể tích hữu hạn có thể coi là biểu thức tương đương với sai phân hữu hạn nhưng viết dưới dạng tích phân (chứ không phải dạng vi phân) của các phương trình cơ bản.

Hãy ví dụ bằng bài toán điều kiện biên (BVP) 1 chiều gần giống như mục 2.4. Xét PT vi phân: $-\frac{d}{dx}\left[p(x) \frac{du}{dx}\right] + q(x) u = f(x)$

Bây giờ, ngoài các điểm nút lưới $\{x_j\}$, ta cần hình dung các thể tích kiểm soát bao quanh mỗi điểm nút (Hình trái của [fig:FVM1]{reference-type=“ref” reference=“fig:FVM1”}). Trong không gian 1 chiều thì thể tích này có dạng một đoạn; chẳng hạn thể tích bao quanh điểm $x_j$ sẽ là đoạn $[x_{j-}, x_{j+}]$. Khi đó ta lấy tích phân của PT trên đoạn này:

$-\left[p(x)\frac{du}{dx}\right]_{x_{j-}}^{x_{j+}} + \int_{x_{j-}}^{x_{j+}} q(x) u dx = \int_{x_{j-}}^{x_{j+}} f(x) dx$

Dùng hai giá trị trung gian cho tham biến $p(x)$ và sử dụng quy tắc “điểm giữa” để tính hai tích phân ta được: $-p_{j+} \frac{u_{j+1} - u_j}{\Delta x} + p_{j-}\frac{u_j - u_{j-1}}{\Delta x} + \Delta x q_j u_j = \Delta x f_j$

Từ đó đưa về phương trình: $-a_j u_{j-1} + b_j u_j - c_j u_{j+1} = \Delta x^2 f_j$

Có thể thấy rằng đây là hệ 3 đường chéo đối xứng với $a_{j+1} = p_{j+} = c_j$.

Bài tập. Cho thấy rằng PT sai phân [discr-FVM-BVP] cũng đạt độ chính xác cấp 2.

Mở rộng bài toán cho trường hợp 2 chiều, xét bài toán điều kiện biên (với $h$ là biến lượng): $-\frac{\partial }{\partial x}\left[p \frac{\partial h}{\partial x}\right] -\frac{\partial }{\partial y}\left[q \frac{\partial h}{\partial y}\right] + gh = f(x,y)$

Hình. Phân định thể tích kiểm soát bao quanh mỗi điểm nút [@Popescu2014]. Bên trái: lưới 1 chiều. Bên phải: lưới 2 chiều.

Để thích hợp với bài mới, ta xét trường hợp 2 chiều với lưới điểm nút (Hình phải của [FVM1]). Mỗi điểm $(x_j, y_k)$ sẽ được bao quanh bởi “thể tích kiểm soát” (CV) là hình chữ nhật $[x_{j-}, x_{j+}] \times [y_{k-}, y_{k+}]$. Diện tích của hình chữ nhật này bằng $(\Delta x_- + \Delta x_+) \times (\Delta y_- + \Delta y_+)$ (dĩ nhiên bằng $\Delta x \times \Delta y$ trong trường hợp lưới đều). Lấy tích phân trên miền chữ nhật này theo cách tương tự như PT một chiều [BVP-1D-int]. Nhưng lần này tích phân của các biểu thức ngoặc vuông không thể được thay cận một cách dễ dàng nữa. Ta cần dùng một công cụ: định lý Gauss (Ostrogradsky). Với vec-tơ biến lượng: $\mathbf{u} = \left( p\frac{\partial h}{\partial x}, q\frac{\partial h}{\partial y} \right)^T$

thì hai số hạng đầu trong PT [eq:BVP-2d] sẽ là⁴: $-\mathrm{div} \mathbf{u}$. Định lý Gauss cho phép chuyển từ tích phân miền sang tích phân biên: $\int_{V_{j,k}} \mathrm{div} \mathbf{u} \; dV = \int_{S_{j,k}} \mathbf{u \cdot n} \; dS$

Hình. Vị trí các nút lân cận, các điểm giữa của 4 mặt một thể tích kiểm soát trên lưới

Tích phân vế phải tiện cho tính toán vì biên $S$ chỉ gồm 4 đoạn thẳng là 4 cạnh của miền chữ nhật $V$. Chẳng hạn với mặt “bên phải” hay “phía đông (e)” (Hình [grid]), số hạng của tích phân là: $(\mathbf{u \cdot n})_e \times \Delta S_e$ = (hình chiếu của vec-tơ $\mathbf{u}$ trên trục $x) \times \Delta y = \displaystyle{ \left( p\frac{\partial h}{\partial x} \right)_e } \times \Delta y = p_{j+, k} \displaystyle{\frac{h_{j+1,k} - h_{j,k}} {\Delta x}} \Delta y$. Như vậy, thừa số thứ nhất đã được tính tại điểm e (điểm $x_{j+}, y_k$) với sai phân trung tâm có độ chính xác cấp hai.

Tương tự, có thể tính các số hạng còn lại cho các cạnh s, w, n. Hai số hạng còn lại, $\displaystyle{\int_{V_{j,k}} gh dV \approx g h_{j,k} \Delta x \Delta y }$ và $\displaystyle{\int_{V_{j,k}} f dV \approx f_{j,k} \Delta x \Delta y }$.

Ở đây, cần giả thiết rằng $f$ phân bố đều trong ô còn mỗi thông lượng phân bố đều trên từng mặt bên.

Ví dụ trên cho thấy FVM hoạt động tốt trên lưới chữ nhật, cũng như phương pháp sai phân hữu hạn. Tuy nhiên lợi thế của FVM thấy rõ khi các ô tính bị biến dạng (đặc biệt là tại biên) hoặc một phần cạnh ô bị chặn thông lượng qua. FVM còn có ưu điểm là cho phép lập các lưới phi cấu trúc, giúp tính toán trên những vùng có dạng đường biên phức tạp và những bài toán yêu cầu làm mịn lưới một cách linh hoạt. Trong mọi trường hợp, định lý Gauss đều cho phép xây dựng công thức mà không cần phải quy về các phương $x,y$.

Biểu diễn thông lượng. PT cân bằng thông lượng như một tổ hợp tuyến tính các biến lượng. Từ PT khuếch tán, có thể dùng FVM đẻ đưa về phương trình tuyến tính, ví dụ như PT [FVM-1D], và có thể dễ dàng mở rộng sang trường hợp 2 chiều. Trong trường hợp có số hạng nguồn $f$, người ta thường tuyến tính hóa số hạng này dưới dạng $f_P = f_0 + k_P \phi_P$. Khi đó phương trình cân bằng cho biến lượng $\phi$ đối với thể tích kiểm soát sẽ được đưa về dạng: $a_P \cdot \phi_P = \sum a_{nb} \cdot \phi_{nb} + f_0$, trong đó chỉ số nb là để chỉ các ô lân cận (N, S, E, W).

Hãy viết PT tuyến tính liên hệ biến lượng $h$ tại các CV áp dụng cho trường hợp PT [BVP-2d]. Gợi ý: số hạng nguồn là ($f_P - g \cdot h_P$). Sau đó cho thấy rằng các hệ số $a_{nb} = \{a_W, a_E, a_S, a_N\}$ chỉ phụ thuộc vào $p, q, \Delta x, \Delta y$; còn $a_P = a_W + a_E + a_S + a_N + g$.

Hệ PT Navier-Stokes cho chất lỏng không nén

Xét hệ PT Navier-Stokes 2D cho chất lỏng không nén được ($\rho=$ const), thường dùng trong thuỷ khí động lực công nghiệp: $\begin{aligned} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} &= 0 \\ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + \frac{1}{\rho} \frac{\partial p}{\partial x} &= \nu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \\ \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + \frac{1}{\rho} \frac{\partial p}{\partial y} &= \nu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right) \end{aligned}$