Phương Pháp Thể Tích Hữu Hạn (Finite Volume Method - FVM) trong Vật lý

1. Giới thiệu

Phương pháp Thể Tích Hữu Hạn (FVM) là một phương pháp số được sử dụng rộng rãi để giải gần đúng các phương trình vi phân đạo hàm riêng (Partial Differential Equations - PDEs), đặc biệt là các phương trình bảo toàn (Conservation Equations) trong vật lý và kỹ thuật. FVM nổi bật nhờ tính bảo toàn (conservativeness) cục bộ, tức là nó đảm bảo các định luật bảo toàn vật lý (ví dụ: bảo toàn khối lượng, năng lượng, động lượng) được tuân thủ trong quá trình tính toán. Tài liệu này trình bày chi tiết về FVM, các bước thực hiện, ưu nhược điểm và ứng dụng của nó trong vật lý.

2. Nguyên lý cơ bản

FVM dựa trên việc tích phân các phương trình bảo toàn trên các thể tích kiểm soát (Control Volumes - CV) rời rạc, không chồng lấn, bao phủ toàn bộ miền tính toán. Thay vì tìm nghiệm điểm như phương pháp sai phân hữu hạn (Finite Difference Method - FDM), FVM tính các giá trị trung bình của biến số trên mỗi CV.

2.1. Phương trình bảo toàn tổng quát

Một phương trình bảo toàn tổng quát có thể được viết dưới dạng:

$\frac{\partial (\rho \phi)}{\partial t} + \nabla \cdot (\rho \phi \vec{v}) = \nabla \cdot (\Gamma \nabla \phi) + S_{\phi}$

Trong đó:

• $\phi$: Biến số được bảo toàn (ví dụ: nhiệt độ, nồng độ, vận tốc)
• $\rho$: Mật độ
• $\vec{v}$: Vận tốc
• $\Gamma$: Hệ số khuếch tán
• $S_{\phi}$: Nguồn (Source term)

2.2. Tích phân trên thể tích kiểm soát

Bước quan trọng nhất của FVM là tích phân phương trình bảo toàn tổng quát trên một CV bất kỳ $\Omega$:

$\int_{\Omega} \frac{\partial (\rho \phi)}{\partial t} dV + \int_{\Omega} \nabla \cdot (\rho \phi \vec{v}) dV = \int_{\Omega} \nabla \cdot (\Gamma \nabla \phi) dV + \int_{\Omega} S_{\phi} dV$

2.3. Định lý Gauss - Ostrogradsky

Để chuyển các tích phân thể tích của các số hạng đạo hàm thành tích phân mặt, ta sử dụng định lý Gauss - Ostrogradsky (định lý phân kỳ):

$\int_{\Omega} \nabla \cdot \vec{F} dV = \oint_{\partial \Omega} \vec{F} \cdot \vec{n} dA$

Trong đó:

• $\vec{F}$ là một trường vector
• $\partial \Omega$ là biên của CV $\Omega$
• $\vec{n}$ là vector pháp tuyến đơn vị hướng ra ngoài biên

Áp dụng định lý Gauss, phương trình trên trở thành:

$\int_{\Omega} \frac{\partial (\rho \phi)}{\partial t} dV + \oint_{\partial \Omega} (\rho \phi \vec{v}) \cdot \vec{n} dA = \oint_{\partial \Omega} (\Gamma \nabla \phi) \cdot \vec{n} dA + \int_{\Omega} S_{\phi} dV$

2.4. Rời rạc hóa thời gian

Số hạng đạo hàm theo thời gian có thể được rời rạc hóa bằng các lược đồ sai phân hữu hạn. Ví dụ, sử dụng lược đồ sai phân lùi (Backward Difference):

$\int_{\Omega} \frac{\partial (\rho \phi)}{\partial t} dV \approx \frac{(\rho \phi)^{n+1} - (\rho \phi)^{n}}{\Delta t} V$

Trong đó:

• $n$ là chỉ số thời gian hiện tại
• $n+1$ là chỉ số thời gian tiếp theo
• $\Delta t$ là bước thời gian
• $V$ là thể tích của CV

2.5. Rời rạc hóa không gian

Các tích phân mặt được rời rạc hóa bằng cách tính toán các thông lượng (fluxes) qua các mặt của CV. Thông lượng bao gồm hai thành phần chính:

• Thông lượng đối lưu (Convection Flux): $\int (\rho \phi \vec{v}) \cdot \vec{n} dA$
• Thông lượng khuếch tán (Diffusion Flux): $\int (\Gamma \nabla \phi) \cdot \vec{n} dA$

Việc tính toán chính xác các thông lượng này rất quan trọng để đảm bảo độ chính xác và ổn định của phương pháp. Có nhiều lược đồ khác nhau để ước tính các thông lượng, ví dụ:

• Lược đồ Upwind: Ước tính giá trị của $\phi$ tại mặt bằng giá trị của $\phi$ tại nút ngược dòng.
• Lược đồ Central Difference: Ước tính giá trị của $\phi$ tại mặt bằng trung bình của các giá trị của $\phi$ tại hai nút lân cận.
• Lược đồ QUICK (Quadratic Upstream Interpolation for Convective Kinematics): Sử dụng phép nội suy bậc hai để ước tính giá trị của $\phi$ tại mặt.

3. Các bước thực hiện FVM

1. Chia miền tính toán thành các CV rời rạc: Lựa chọn lưới (mesh) phù hợp với hình dạng miền tính toán và độ chính xác mong muốn.
2. Tích phân phương trình bảo toàn trên mỗi CV: Áp dụng tích phân trên CV cho phương trình bảo toàn.
3. Áp dụng định lý Gauss - Ostrogradsky: Chuyển các tích phân thể tích thành tích phân mặt.
4. Rời rạc hóa các tích phân mặt: Sử dụng các lược đồ khác nhau (Upwind, Central Difference, QUICK,...) để tính toán thông lượng.
5. Rời rạc hóa số hạng nguồn: Ước tính giá trị trung bình của số hạng nguồn trong CV.
6. Thiết lập hệ phương trình đại số: Sau khi rời rạc hóa, ta thu được một hệ phương trình đại số tuyến tính hoặc phi tuyến cho các giá trị $\phi$ tại các nút.
7. Giải hệ phương trình: Sử dụng các phương pháp số như phương pháp lặp Gauss-Seidel, phương pháp SOR (Successive Over-Relaxation), hoặc các phương pháp ma trận thưa để giải hệ phương trình.

4. Ưu điểm và nhược điểm của FVM

4.1. Ưu điểm

• Tính bảo toàn: FVM đảm bảo tính bảo toàn cục bộ, là một ưu điểm lớn khi giải các bài toán bảo toàn.
• Tính linh hoạt: FVM có thể được áp dụng cho các hình dạng miền tính toán phức tạp và các loại lưới khác nhau (lưới cấu trúc, lưới phi cấu trúc).
• Dễ dàng xử lý các điều kiện biên phức tạp: Các điều kiện biên (Boundary Conditions - BCs) có thể được áp dụng trực tiếp vào các tích phân mặt.

4.2. Nhược điểm

• Phức tạp hơn FDM: FVM đòi hỏi nhiều bước rời rạc hóa hơn FDM, đặc biệt là việc tính toán thông lượng.
• Độ chính xác: Độ chính xác của FVM phụ thuộc vào lược đồ rời rạc hóa và chất lượng lưới.
• Tính toán: Việc giải hệ phương trình đại số có thể tốn nhiều thời gian tính toán đối với các bài toán lớn.

5. Ứng dụng của FVM trong Vật lý

FVM được sử dụng rộng rãi trong nhiều lĩnh vực của vật lý, bao gồm:

• Cơ học chất lưu (Computational Fluid Dynamics - CFD): Mô phỏng dòng chảy chất lỏng và chất khí, truyền nhiệt, chuyển khối.
• Truyền nhiệt: Giải các phương trình dẫn nhiệt, đối lưu nhiệt, bức xạ nhiệt.
• Điện từ học (Electromagnetics): Giải các phương trình Maxwell.
• Vật lý plasma (Plasma Physics): Mô phỏng plasma.
• Vật lý chất rắn (Solid State Physics): Tính toán phân bố nhiệt trong vật rắn, mô phỏng quá trình khuếch tán.

6. Ví dụ minh họa (1D)

Xét phương trình dẫn nhiệt 1D ổn định (steady-state) không nguồn:

$\frac{d}{dx} \left( k \frac{dT}{dx} \right) = 0$

Trong đó:

• $T$ là nhiệt độ
• $k$ là độ dẫn nhiệt

1. Chia miền 1D thành các CV: Chia miền từ $x_W$ đến $x_E$ thành các CV, mỗi CV có nút trung tâm $P$ và hai mặt $w$ và $e$.

2. Tích phân trên CV: Tích phân phương trình trên CV từ $x_w$ đến $x_e$:

$\int_{x_w}^{x_e} \frac{d}{dx} \left( k \frac{dT}{dx} \right) dx = 0$

3. Áp dụng định lý Gauss:

$\left( k \frac{dT}{dx} \right) \Big|_{x_e} - \left( k \frac{dT}{dx} \right) \Big|_{x_w} = 0$

4. Rời rạc hóa:

$k_e \frac{T_E - T_P}{\Delta x} - k_w \frac{T_P - T_W}{\Delta x} = 0$

Trong đó:

• $T_P$, $T_W$, $T_E$ là nhiệt độ tại các nút $P$, $W$, $E$
• $k_e$, $k_w$ là độ dẫn nhiệt tại các mặt $e$ và $w$
• $\Delta x$ là kích thước CV

5. Sắp xếp lại:

$a_P T_P = a_W T_W + a_E T_E$

Với:

• $a_W = \frac{k_w}{\Delta x}$
• $a_E = \frac{k_e}{\Delta x}$
• $a_P = a_W + a_E$

Đây là phương trình rời rạc hóa cho một CV trong bài toán dẫn nhiệt 1D. Các phương trình tương tự được thiết lập cho tất cả các CV, tạo thành một hệ phương trình đại số cần giải.

7. Kết luận

Phương pháp Thể Tích Hữu Hạn (FVM) là một công cụ mạnh mẽ và linh hoạt để giải các phương trình bảo toàn trong vật lý. Tính bảo toàn, khả năng xử lý các hình dạng miền phức tạp và điều kiện biên đa dạng làm cho FVM trở thành một lựa chọn phổ biến trong mô phỏng số. Việc hiểu rõ nguyên lý, các bước thực hiện và các lược đồ rời rạc hóa là rất quan trọng để áp dụng FVM hiệu quả trong giải quyết các bài toán vật lý.