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

Trở về Mục lục cuốn sách

Giới thiệu chung

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

T(r, t) / ∂t = D ∇2T(r, t), 

trong đó D > 0 là hệ số khuếch tán (đều), mô tả nhiều hiện tượng vật lý thú vị. Chẳng hạn, định luật truyền nhiệt có thể được viết là

q =  − κ ∇T, 

trong đó q là luồng nhiệt, T là nhiệt độ, và κ là hệ số truyền nhiệt. Phương trình trên phát biểu một cách đơn giản rằng nhiệt truyền theo chiều dốc của gra-đi-en nhiệt độ. Nếu không có nguồn phát hoặc thu nhiệt, định luật bảo toàn năng lượng đòi hỏi

 − ∂Q / ∂t = ∫ q ⋅ dS, 

trong đó Q là nhiệt năng chứa trong một thể tích V nào đó được bao bọc bởi một mặt kín S. Phương trình trên phát biểu rằng tốc độ giảm của nhiệt năng chứa trong thể tích V thì bằng luồng nhiệt tức thời truyền qua mặt biên của nó. Tuy nhiên, Q = ∫ cTdV, với c là nhiệt dung của mỗi đơn vị thể tích. Từ các phương trình trên, và định lý Gauss (Ostrogradsky), ta được phương trình khuếch tán nhiệt độ sau:

T / ∂t = D ∇2T, 

trong đó D = κ / c. Ở bài toán truyền nhiệt điển hình, ta được cho trước nhiệt độ T(r, t0) tại thời điểm ban đầu t0, và được yêu cầu tính T(r, t) tại tất cả thời điểm sau đó. Một bài toán như vậy được gọi là bài toán giá trị đầu. Các điều kiện biên không gian có thể thuộc kiểu Dirichlet (T được chỉ định trên biên), kiểu Neumann (T được chỉ định trên biên), hoặc một sự kết hợp giữa hai kiểu trên.

Bài toán 1 chiều với các điều kiện biên hỗn hợp

Xét cách giải phương trình khuếch tán trong một chiều không gian. Giả sử

T(x, t) / ∂t = D ∂2T(x, t) / ∂x2,   (diff)

với xl ≤ x ≤ xh, tuân theo điều kiện biên không gian kiểu hỗn hợp

αl(t) T(x, t) + βl(t) ∂T(x, t) / ∂x = γl(t), 

tại x = xl, và

αh(t) T(x, t) + βh(t) ∂T(x, t) / ∂x = γh(t), 

tại x = xh. Ở đây, αl, βl, v.v., là các hàm theo thời gian đã biết. Dĩ nhiên, T(x, t0) phải được chỉ rõ ở một thời điểm ban đầu t0 nào đó.

Phương trình ([diff]) cần được làm rời rạc theo cả thời gian và không gian. Theo thời gian, ta làm rời rạc trên lưới có các nút cách đều nhau

tn = t0 + nδt, 

trong đó δtbước thời gian. Dùng một lược đồ sai phân bậc nhất đơn giản, ta đưa PT ([diff]) về dạng

\displaystyle{ \frac{T(x,t_{n+1}) - T(x,t_n)}{\delta t} = D \,\frac{\partial^2 T(x,t_n)}{\partial x^2}+ O(\delta t). }         (diff1)

Theo không gian, ta làm rời rạc trên lưới có các điểm nút cách đều nhau như đã đề cập trong PT ([e57]), và xấp xỉ d2 / dx2 bằng sơ đồ sai phân trung tâm bậc hai đã giới thiệu ở PT ([e58]). Các điều kiện biên không gian được làm rời rạc theo cách giống như các PT ([e529]) và ([e530]). Như vậy, PT ([diff1]) cho ta

\displaystyle{ \frac{T_i^{n+1}-T_i^n}{\delta t} = D\,\frac{T_{i-1}^n-2\,T_i^n+T_{i+1}^n}{(\delta x)^2}, }

hay

Tin + 1 = Tin + C (Ti − 1n − 2 Tin + Ti + 1n)  (diff2)

với i = 1, N, trong đó Tin ≡ T(xi, tn), và C = Dδt / (δx)2. Các điều kiện biên được làm rời rạc có dạng

\begin{array}{rcl}  T_0^n &=& \displaystyle{\frac{\gamma_l^n\,\delta x-\beta_l^n\,T_1^n}{\alpha_l^n\,\delta x -\beta_l^n},}\\[0.5ex]  T_{N+1}^n &=& \displaystyle{\frac{\gamma_h^n\,\delta x + \beta_h^n\,T_N^n}{\alpha_h^n\,\delta x +\beta_h^n},}  \end{array}             (diff3)

trong đó γln ≡ γl(tn), v.v. Sơ đồ rời rạc nêu trên được gọi là bậc nhất theo thời gian và bậc hai trong không gian.

Các PT ([diff2])–([diff3]) hợp thành một lược đồ lặp khá dễ thực hiện, mà ta có thể dùng để dẫn ra các giá trị T(x, t) theo thời gian.

Một chương trình ví dụ để giải phương trình khuếch tán 1 chiều

Dưới đây là một ví dụ về đoạn chương trình để giải phương trình khuếch tán 1 chiều trong đó có dùng thư viện Blitz++.

// Diffusion1D.cpp 

// Function to evolve diffusion equation in 1-d: 

//  dT / dt = D d^2 T / dx^2  for  xl <= x <= xh 

//  alpha_l T + beta_l dT/dx = gamma_l  at x=xl 

//  alpha_h T + beta_h dT/dx = gamma_h  at x=xh 

// Array T assumed to be of extent N+2. 

// Now, ith element of array corresponds to 

//  x_i = xl + i * dx    i=0,N+1 

// Here, dx = (xh - xl) / (N+1) is grid spacing. 

// Function evolves T by single time-step. 

// C = D dt / dx^2, where dt is time-step. 

// Uses explicit scheme. 

#include <blitz/array.h> 

using namespace blitz; 

void Diffusion1D (Array<double,1>& T,  
                  double alpha_l, double beta_l, double gamma_l, 
                  double alpha_h, double beta_h, double gamma_h, 
                  double dx, double C) 
{ 
  // Set N. Declare local array. 
  int N = T.extent(0) - 2; 
  Array<double,1> T0(N+2);     

  // Evolve T 
  T0 = T;     
  for (int i = 1; i <= N; i++) 
    T(i) += C * (T0(i-1) - 2. * T0(i) + T0(i+1)); 

  // Set boundary conditions 
  T(0) = (gamma_l * dx - beta_l * T(1)) / (alpha_l * dx - beta_l); 
  T(N+1) = (gamma_h * dx + beta_h * T(N)) / (alpha_h * dx + beta_h); 
}

Ví dụ giải phương trình khuếch tán một chiều

Ta hãy giải phương trình khuếch tán trong không gian một chiều bằng kĩ thuật sai phân hữu hạn như đã thảo luận từ trước. Ta tìm nghiệm của PT ([diff]) trong miền  − x0 ≤ x ≤ x0, tuân theo điều kiện ban đầu

T(x,t_0) = \exp\!\left(\frac{-x^2}{4\,D\,t_0}\right),

trong đó t0 > 0. Các điều kiện biên không gian bao gồm

T(\pm x_0,t) = \sqrt{\frac{t_0}{t}}\, \exp\!\left(\frac{-x_0^2}{4\,D\,t}\right).

Dĩ nhiên là ta có thể giải bài toán này theo cách giải tích để có

T(x,t) =\sqrt{\frac{t_0}{t}}\, \exp\!\left(\frac{-x^2}{4\,D\,t}\right).

Lưu ý rằng phương trình trên biểu thị một xung Gauss vốn sẽ dần giảm chiều cao và tăng bề rộng sao cho diện tích của nó không đổi. Bề rộng của xung biến đổi xấp xỉ bằng

\Delta x \sim \sqrt{D\,t}.

Ngoài ra, xung này còn tiệm cận một hàm δ khi t → 0.

diff1

Hình. Sự phát triển theo khuynh hướng khuếch tán của xung Gauss trong không gian 1 chiều. Tính toán số trị được thực hiện với D = 1; x0 = 5; δt = 4 × 10−3, và N = 100. Thời gian mô phỏng từ t = 0,1 đến t = 1. Đường liền nét biểu thị điều kiện ban đầu lúc t = 0,1, đường đứt nét biểu thị nghiệm số trị lúc t = 1, và đường chấm chấm (bị khuất sau đường đứt nét) là nghiệm giải tích lúc t = 1. [fdiff1]

Hình [fdiff1] so sánh giữa nghiệm giải tích và nghiệm số trị tính toán với D = 1; x0 = 5; t0 = 0,1; δt = 4 × 10−3, và N = 100. Có thể tháy rằng các nghiệm giải tích và số trị rất khớp nhau.

diff2

Hình. Sự phát triển theo khuynh hướng khuếch tán của xung Gauss trong không gian 1 chiều. Tính toán số trị được thực hiện với D = 1; x0 = 5; δt = 4 × 10−3, và N = 125. Thời gian bắt đầu mô phỏng là t = 0,1. Các hình nhỏ ở phía trên bên trái, phía trên bên phải, phía dưới bên trái và phía dưới bên phải lần lượt biểu thị cho nghiệm lúc t = 0,45; t = 0,46; t = 0,47, và t = 0,48.[fdiff2]

Cũng có lý nếu ta mong đợi rằng khi N tăng tại những δt cố định (nghĩa là độ phân giải không gian tăng tại những độ phân giải thời gian cố định) thì nghiệm số trị sẽ càng trở nên chính xác hơn. Điều này đúng—ít nhất là đến tận khi N vượt quá một giá trị tới hạn. Vượt ngoài giá trị này, sẽ có một sự đổ vỡ trầm trọng ở nghiệm số. Sự đổ vỡ này được minh họa trên Hình [fdiff2]. Có thể thấy rằng trong nghiệm đã hình thành những dao động có bước sóng ngắn và nhanh chóng khuếch đại. Thực sự là cuối cùng nghiệm sẽ trở nên vô hạn. Ta hãy khảo sát hiện tượng bất thường và có phần rắc rối này.

Phân tích ổn định von Neumann

Rõ ràng là thuật toán sai phân hữu hạn đơn giản đã dùng để giải phương trình khuếch tán 1 chiều đã bị bất ổn định số trị trong những trường hợp nhất định. Ta hãy thử tìm xem khi nào sự bất ổn định này xảy ra. Xét sự diễn biến theo thời gian của một hình sóng Fourier đơn lẻ có số sóng k:

T(x, t) = \hat{T} (t) e i kx. 

Thay biểu thức trên vào sơ đồ sai phân hữu hạn ([diff2]) ta được

\hat{T}n + 1 e i kxn = \hat{T}n e i kxn [1 + C (e − i kδx − 2 + e + i kδx)], 

hay

\hat{T}n + 1 = A\hat{T}n, 

trong đó

A = 1 − 2 C (1 − coskδx) = 1 − 4 C sin2(kδx / 2).   (stable)

Như vậy, biên độ của hình sóng Fourier này được khuếch đại bởi một hệ số A sau mỗi bước thời gian. Muốn lược đồ sai phân được ổn định, độ lớn của hệ số khuếch đại phải nhỏ hơn đơn vị với tất cả các giá trị mà k có thể nhận. Ở đây, giá trị lớn nhất của sin2(kδx / 2) bằng 1: vì vậy, độ dài sóng ứng với giá trị này chính là hình sóng Fourier không ổn định nhất. Thật ra, hình sóng không ổn định nhất có độ dài bước sóng bằng một nửa khoảng cách nút lưới, tức là λ = δx / 2. Từ PT ([stable]) suy ra hình sóng này sẽ ổn định nếu

C < 1 / 2. 

Sau cùng, từ định nghĩa của C, điều kiện ổn định cần tìm có thể được viết dưới dạng

δt < (δx)2  /  2 D.   (limit)

Lưu ý rằng C = 0. 408 ứng với tính toán ổn định chỉ ra trên Hình [fdiff1], còn C = 0. 635 ứng với tính toán không ổn định ở Hình [fdiff2]. Cần nói thêm rằng dạng phân tích ổn định nói trên được gọi là phân tích ổn định von Neumann. Lưu ý rằng việc bỏ qua các điều kiện biên không gian trong tính toán trên có thể được lượng thứ vì các dạng sóng bất ổn định biến đổi trên những cỡ độ dài rất nhỏ, điển hình là bằng cỡ khoảng cách nút lưới.

Theo PT ([limit]), lược đồ sai phân đang xét để giải phương trình khuếch tán một chiều chỉ ổn định khi bước thời gian nhỏ hơn một giá trị tới hạn nào đó. Lưu ý rằng giá trị tới hạn này tỉ lệ với bình phương của khoảng cách nút lưới. Đây là một tỉ lệ không có lợi, vì nó cho thấy mỗi khi tăng gấp đôi độ phân giải trong không gian thì đồng thời phải giảm bước thời gian đi 4 lần để có thể duy trì sự ổn định số. Dĩ nhiên, ràng buộc trên khiến bước thời gian của ta bị thu hẹp đến mức quá nhỏ trong những tính toán có độ phân giải cao. Liệu có cách nào khắc phục được sự ràng buộc này không?

Lược đồ Crank-Nicolson

Ta hãy xem lại lược đồ sai phân theo thời gian:

\displaystyle{ \frac{T(x,t_{n+1}) - T(x,t_n)}{\delta t}  = D \,\frac{\partial^2 T(x,t_n)}{\partial x^2}+ O(\delta t). }

Lưu ý rằng vế phải được lượng giá hoàn toàn tại điểm đầu của bước thời gian, tức là tại tn. Kiểu lược đồ sai phân thời gian này được gọi là lược đồ hiện. Lược đồ hiện rất dễ tiến hành lập trình, nhưng lại đặc biệt dễ nhiễm phải bất ổn định số. Thật may là ta thường có thể vượt qua những bất ổn định này bằng cách làm cho lược đồ sai phân của ta có tính ẩn. Một lược đồ ẩn là lược đồ mà trong đó vế phải được lượng giá một phần (hay hoàn toàn) ở cuối bước thời gian, tức là tại tn+1. Điều không may là những lược đồ ẩn nói chung khó thực hiện hơn rất nhiều so với lược đồ hiện.

Phương pháp ẩn Crank-Nicolson nổi tiếng để giải phương trình khuếch tán bao gồm việc lấy trung bình vế phải giữa lúc bắt đầu và kết thúc bước thời gian. Nói cách khác,

\displaystyle{ \frac{T(x,t_{n+1}) - T(x,t_n)}{\delta t}  = \frac{D}{2} \,\frac{\partial^2 T(x,t_n)}{\partial x^2}+  \frac{D}{2} \,\frac{\partial^2 T(x,t_{n+1})}{\partial x^2}+ O(\delta t)^2. }

Như được chỉ ra ở số hạng sai số, phương pháp này thật ra là bậc hai theo thời gian.

Dùng lược đồ sai phân không gian ta vẫn thường áp dụng, biểu thức trên sẽ cho

T_i^{n+1} - \frac{C}{2}\,\left(T_{i-1}^{n+1}-2\,T_i^{n+1}+T_{i+1}^{n+1}\right)=  T_i^n + \frac{C}{2}\,\left(T_{i-1}^n-2\,T_i^n+T_{i+1}^n\right).                      (crank)

Khi ta thực hiện phân tích ổn định von Neumann cho lược đồ trên, ta thu được biểu thức sau cho hệ số khuếch đại:

\displaystyle{ A = \frac{1 - 2\,C\,\sin^2(k\,\delta x/2)}{1+2\,C\,\sin^2(k\,\delta x/2)}.}

Lưu ý rằng A∣ < 1 với mọi giá trị của k. Từ đó cho thấy lược đồ Crank-Nicolson ổn định không điều kiện. Không may là, PT ([crank]) chứa một phương trình ma trận ba đường chéo liên hệ các giá trị Tin+1Tin. Do vậy, cái giá mà ta phải trả cho độ chính xác cao và ổn định không điều kiện của lược đồ Crank-Nicolson là phải thực hiện nghịch đảo phương trình ma trận ba đường chéo ở mỗi bước thời gian. Thường thì cái giá phải trả này cũng rất hợp lý.

Một chương trình được cải tiến để giải phương trình khuếch tán 1 chiều

Sau đây là một chương trình đã được cải tiến nhằm giải phương trình khuếch tán 1 chiều trong đó dùng lược đồ Crank-Nicolson, cũng như đoạn chương trình giải ma trận ba đường chéo cùng thư viện Blitz++. Lưu ý sự giống nhau đáng kể về mặt cấu trúc giữa chương trình này và chương trình để giải phương trình Poisson 1 chiều trước đây (xem Mục [poisson1d]).

// CrankNicolson1D.cpp 

// Function to evolve diffusion equation in 1-d: 

//  dT / dt = D d^2 T / dx^2  for  xl <= x <= xh 

//  alpha_l T + beta_l dT/dx = gamma_l  at x=xl 

//  alpha_h T + beta_h dT/dx = gamma_h  at x=xh 

// Array T assumed to be of extent N+2. 

// Now, ith element of array corresponds to 

//  x_i = xl + i * dx    i=0,N+1 

// Here, dx = (xh - xl) / (N+1) is grid spacing. 

// Function evolves T by single time-step. 

// C = D dt / dx^2, where dt is time-step. 

// Uses Crank-Nicolson implicit scheme. 

#include <blitz/array.h> 

using namespace blitz; 

void Tridiagonal (Array<double,1> a, Array<double,1> b, Array<double,1> c,  
                  Array<double,1> w, Array<double,1>& u); 

void CrankNicolson1D (Array<double,1>& T,  
                  double alpha_l, double beta_l, double gamma_l, 
                  double alpha_h, double beta_h, double gamma_h, 
                  double dx, double C) 
{ 
  // Find N. Declare local arrays. 
  int N = T.extent(0) - 2; 
  Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2); 

  // Initialize tridiagonal matrix 
  for (int i = 2; i <= N; i++) a(i) = - 0.5 * C; 
  for (int i = 1; i <= N; i++) b(i) = 1. + C; 
  b(1) += 0.5 * C * beta_l / (alpha_l * dx - beta_l); 
  b(N) -= 0.5 * C * beta_h / (alpha_h * dx + beta_h); 
  for (int i = 1; i <= N-1; i++) c(i) = - 0.5 * C; 

  // Initialize right-hand side vector 
  for (int i = 1; i <= N; i++) 
    w(i) = T(i) + 0.5 * C * (T(i-1) - 2. * T(i) + T(i+1)); 
  w(1) += 0.5 * C * gamma_l * dx / (alpha_l * dx - beta_l); 
  w(N) += 0.5 * C * gamma_h * dx / (alpha_h * dx + beta_h); 

  // Invert tridiagonal matrix equation 
  Tridiagonal (a, b, c, w, T); 

  // Calculate i=0 and i=N+1 values 
  T(0) = (gamma_l * dx - beta_l * T(1)) / 
    (alpha_l * dx - beta_l); 
  T(N+1) = (gamma_h * dx + beta_h * T(N)) / 
    (alpha_h * dx + beta_h); 
}

Một lời giải được cải thiện của phương trình khuếch tán 1 chiều

Bây giờ ta hãy giải phương trình khuếch tán đơn giản đã giới thiệu ở Mục [VD giải PT 1 chiều] với đoạn chương trình Crank-Nicolson nêu trên. Hình [fdiff3] so sánh giữa các nghiệm giải tích và số trị cho bài toán với D = 1, x0 = 5, t0 = 0. 1, δt = 0. 1, và N = 100. Có thể thấy được rằng nghiệm giải tích và số trị rất khớp nhau. Nhưng lưu ý rằng bước thời gian dùng trong tính toán này (δt = 0. 1) lớn hơn rất nhiều so với tính toán trước (δt = 4 × 10 − 3) khi dùng cho sơ đồ sai phân hiện—xem Hình [fdiff1]. Căn cứ vào PT ([limit]), một sơ đồ hiện sẽ bị hạn chế với bước thời gian nhỏ hơn chừng 5 × 10 − 3 với bài toán đang xét. Do vậy, bằng sơ đồ ẩn ta đã vượt được mức này (gấp 20 lần), mà vẫn giữ được tính ổn định số trị. Lưu ý rằng sơ đồ Crank-Nicolson của ta có thể cho phép đạt kết quả chính xác với bước thời gian lên đến 0. 1 vì nó là bậc hai theo thời gian.

diff3

Hình. Sự phát triển khuếch tán từ xung Gauss 1 chiều. Việc tính toán số trị được thực hiện với D = 1 , x0 = 5 , δt = 0,1 , và N = 100. Xung được mô phỏng từ t = 0,1 đến t = 1. Đường nét liền cho thấy điều kiện ban đầu tại t = 0,1, đường nết đứt là nghiệm số trị tại t = 1, còn đường chấm chấm (che khuất bởi đường nét đứt) là nghiệm giải tích tại t = 1. [fdiff3]

Bài toán 2 chiều với điều kiện biên Dirichlet

Ta hãy xét nghiệm của phương trình khuếch tán theo 2 chiều. Giả sử rằng

\displaystyle{ \frac{\partial T(x,y,t)}{\partial t} = D\,\frac{\partial^2 T(x,y,t)}{\partial x^2}+D\,\frac{\partial^2 T(x,y,t)}  {\partial y^2}, \qquad \text{(diff2deq)} }

với xl ≤ x ≤ xh, và 0 ≤ y ≤ L. Coi như T(x, y, t) thỏa mãn điều kiện biên hỗn hợp theo phương x:

αl(t) T(x, y, t) + βl(t) ∂T(x, y, t) / ∂x = γl(y, t), 

tại x = xl, và

αh(t) T(x, y, t) + βh(t) ∂T(x, y, t) / ∂x = γh(y, t), 

tại x = xh. Ở đây, αl, βl, v.v. là các hàm theo t đã biết, còn γl, γh là các hàm theo yt đã biết. Ngoài ra, coi như T(x, y, t) thỏa mãn các điều kiện biên Dirichlet đơn giản sau theo phương y:

T(x, 0, t) = T(x, L, t) = 0. 

Cũng như trước, ta làm rời rạc theo thời gian trên lưới cách đều tn = t0 + nδt, với n = 0, 1, 2, ⋯. Ngoài ra, theo phương x, ta làm rời rạc trên lưới đềy xi = xl + iδx, với i = 0, N + 1, trong đó δx = (xh − xl) / (N + 1). Cuối cùng, theo phương y, ta làm rời rạc trên lưới đều yj = jδy, với j = 0, J, trong đó δy = L / J. Dùng lược đồ sai phân Crank-Nicolson theo thời gian đã xét ở Mục [Crank-Nicolson], và lược đồ sai phân bậc hai theo không gian đã xét ở Mục [poisson1d], PT ([diff2deq]) cho ta

\displaystyle{  \frac{T_{i,j}^{n+1}-T_{i,j}^n}{\delta t} - \frac{D}{2}\,  \frac{T_{i-1,j}^{n+1}-2\,T_{i,j}^{n+1}+T_{i+1,j}^{n+1}}{(\delta x)^2}  - \frac{D}{2}\left(\frac{\partial^2 T}{\partial y^2}\right)_{i,j}^{n+1}= }

\displaystyle{ \frac{D}{2}\,  \frac{T_{i-1,j}^{n}-2\,T_{i,j}^{n}+T_{i+1,j}^{n}}{(\delta x)^2} +  \frac{D}{2}\left(\frac{\partial^2 T}{\partial y^2}\right)_{i,j}^{n}, }            (edisc)

trong đó Ti, jn ≡ T(xi, yj, tn). Điều kiện biên sau khi làm rời rạc có dạng

T_{0,j}^n = \displaystyle { \frac{\gamma_{l\,j}^n\,\delta x-\beta_l^n\,T_{1,j}^n}{\alpha_l^n\,\delta x -\beta_l^n},}
T_{N+1,j}^n = \displaystyle {\frac{\gamma_{h\,j}^n\,\delta x + \beta_h^n\,T_{N,j}^n}{\alpha_h^n\,\delta x +\beta_h^n}, }       (bcxdir1)

Ti, 0n = Ti, jn = 0. $       (bcydir)

Ở đây, αln ≡ αl(tn), v.v., và γljn ≡ γl(yj, tn), v.v.

Áp dụng phương pháp Fourier đã giới thiệu ở Mục [bài toán 2 chiều ĐK biên Dirichlet, chương trước], ta viết các Ti, jn dưới dạng các hàm điều hòa Fourier hình sin:

Ti, jn = ∑ k = 0J\hat{T}i, kn sin(jkπ / J),   (recon)

vốn tự thỏa mãn các điều kiện biên ([bcydir]). Biểu thức trên có thể được nghịch đảo để cho (xem Mục [biến đổi Fourier nhanh, chương trước])

\hat{T}i, jn = (2 / J)∑ k = 0JTi, kn sin(jkπ / J).   (fftthat)

Khi PT ([edisc]) được viết dưới dạng các \hat{T}i, jn, nó rút gọn về

\displaystyle{ -\frac{C}{2}\,\hat{T}_{i-1,j}^{n+1} + \left\{1+C\,(1+j^2\,\kappa^2/2)\right\}\,  \hat{T}_{i,j}^{n+1}-\frac{C}{2}\,\hat{T}_{i+1,j}^{n+1}= }

\displaystyle{ \frac{C}{2}\,\hat{T}_{i-1,j}^{n} + \left\{1-C\,(1+j^2\,\kappa^2/2)\right\}\,  \hat{T}_{i,j}^{n}+\frac{C}{2}\,\hat{T}_{i+1,j}^{n},}      (tri1)

với i = 1, N, và j = 0, J. Ở đây, C = Dδt / (δx)2 , và κ = πδx / L.
Hơn nữa, các điều kiện biên (bcxdir) và (bcxdir1) cho ta

\displaystyle {\hat{T}_{0,j}^n = \frac{\Gamma_{l\,j}^n\,\delta x-\beta_l^n\,\hat{T}_{1,j}^n}{\alpha_l^n\,\delta x -\beta_l^n} }
\displaystyle { \hat{T}_{N+1,j}^n = \frac{\Gamma_{h\,j}^n\,\delta x + \beta_h^n\,\hat{T}_{N,j}^n}{\alpha_h^n\,\delta x +\beta_h^n}, }        (tri3)

trong đó

\hat{\Gamma}_{l\,j}^n = (2/J) \sum_{k=0}^{J} \gamma_{l,k}^n\,\sin(j\,k\,\pi/J),            (fftbc)

v.v. Các PT ([tri1])—([tri3]) tạo thành J + 1 phương trình độc lập nhau dạng ma trận ba đường chéo có các ẩn i, jn + 1; mỗi phương trình có một giá trị của j.

Để có thể tìm nghiệm ở bước thời gian tiếp theo, đầu tiên ta phải biến đổi Fourier các Ti, jn và các điều kiện biên, dựa theo các PT ([fftthat]) và ([fftbc]). Tiếp theo, ta nghịch đảo J + 1 phương trình ba đường chéo ([tri1])—([tri3]) để thu được các i, jn + 1. Cuối cùng, ta tái lập Ti, jn bằng PT ([recon]).

Bài toán 2 chiều với điều kiện biên Neumann

Ta hãy thay các điều kiện biên Dirichlet ([bcydir]) bằng các điều kiện biên Neumann đơn giản sau đây:

T(x, 0, t) / ∂y = ∂T(x, L, t) / ∂y = 0. 

Phương pháp tìm nghiệm đã trình bày ở mục trước vẫn được giữ nguyên, chỉ khác ở chỗ những biến đổi Fourier hình sin được thay bằng biến đổi Fourier hình cosin—xem các Mục [bài toán 2 chiều ĐK biên Neumann] và [biến đổi Fourier nhanh] ở chương trước.

Một chương trình ví dụ để giải phương trình khuếch tán 2 chiều

Dưới đây là một chương trình ví dụ để giải phương trình khuếch tán 2 chiều bằng lược đồ Crank-Nicolson kèm theo công cụ giải ma trận ba đường chéo đã có và thư viện Blitz++. Lưu ý sự giống nhau đáng kể về cấu trúc giữa chương trình này và chương trình giải bài toán Poisson 2 chiều trước đây (xem Mục [poisson2d]).

// CrankNicolson2D.cpp 

// Function to evolve diffusion equation in 2-d: 

//  dT / dt = D d^2 T / dx^2 +  D d^2 T / dy^2  for  xl <= x <= xh   
//                                              and  0 <= y <= L 

//  alphaL T + betaL dT/dx = gammaL(y)  at x=xl 

//  alphaH T + betaH dT/dx = gammaH(y)  at x=xh 

// In y-direction, either simple Dirichlet boundary conditions: 

//  T(x,0) = T(x,L) = 0 

// or simple Neumann boundary conditions: 

//  dT/dy(x,0) = dT/dy(x,L) = 0 

// Matrix T assumed to be of extent N+2, J+1. 
// Arrays gammaL, gammaH assumed to be of extent J+1. 

// Now, (i,j)th elements of matrices correspond to 

//  x_i = xl + i * dx    i=0,N+1 

//  y_j = j * L / J      j=0,J 

// Here, dx = (xh - xl) / (N+1) is grid spacing in x-direction. 

// Now, C =  D dt / dx^2, and kappa = pi * dx / L 

// Finally, Neumann=0/1 selects Dirichlet/Neumann bcs in y-direction. 

// Uses Crank-Nicolson scheme. 

#include <blitz/array.h> 

using namespace blitz; 

void fft_forward_cos (Array<double,1> f, Array<double,1>& F); 
void fft_backward_cos (Array<double,1> F, Array<double,1>& f); 
void fft_forward_sin (Array<double,1> f, Array<double,1>& F); 
void fft_backward_sin (Array<double,1> F, Array<double,1>& f); 
void Tridiagonal (Array<double,1> a, Array<double,1> b, Array<double,1> c,  
                  Array<double,1> w, Array<double,1>& u); 

void CrankNicolson2D (Array<double,2>& T, 
                       double alphaL, double betaL, Array<double,1> gammaL,  
                       double alphaH, double betaH, Array<double,1> gammaH, 
                       double dx, double C, double kappa, int Neumann) 
{ 
  // Find N and J. Declare local arrays. 
  int N = T.extent(0) - 2; 
  int J = T.extent(1) - 1; 
  Array<double,2> TT(N+2, J+1), V(N+2, J+1); 
  Array<double,1> GammaL(J+1), GammaH(J+1); 

  // Fourier transform T   
  for (int i = 0; i <= N+1; i++) 
    { 
      Array<double,1> In(J+1), Out(J+1); 

      for (int j = 0; j <= J; j++) In(j) = T(i, j); 

      if (Neumann) 
        fft_forward_cos (In, Out); 
      else 
        fft_forward_sin (In, Out); 

      for (int j = 0; j <= J; j++) TT(i, j) = Out(j); 
    } 

  // Fourier transform boundary conditions 
  if (Neumann) 
    { 
      fft_forward_cos (gammaL, GammaL); 
      fft_forward_cos (gammaH, GammaH); 
    } 
  else 
    { 
      fft_forward_sin (gammaL, GammaL); 
      fft_forward_sin (gammaH, GammaH); 
    } 

  // Construct source term 
  for (int i = 1; i <= N; i++) 
    for (int j = 0; j <= J; j++)  
      V(i, j) =  
        0.5 * C * TT(i-1, j) +  
        (1. - C * (1. + 0.5 * double (j * j) * kappa * kappa)) * TT(i, j) + 
        0.5 * C * TT(i+1, j); 

  // Solve tridiagonal matrix equations 
  if (Neumann) 
    { 
      for (int j = 0; j <= J; j++) 
       {   
          Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2), u(N+2); 

          // Initialize tridiagonal matrix 
          for (int i = 2; i <= N; i++) a(i) = - 0.5 * C; 
          for (int i = 1; i <= N; i++) 
           b(i) = 1. + C * (1. + 0.5 * double (j * j) * kappa * kappa); 
          b(1) += 0.5 * C * betaL / (alphaL * dx - betaL); 
          b(N) -= 0.5 * C * betaH / (alphaH * dx + betaH); 
          for (int i = 1; i <= N-1; i++) c(i) = - 0.5 * C; 

          // Initialize right-hand side vector 
          for (int i = 1; i <= N; i++) 
           w(i) = V(i, j); 
          w(1) += 0.5 * C * GammaL(j) * dx / (alphaL * dx - betaL); 
          w(N) += 0.5 * C * GammaH(j) * dx / (alphaH * dx + betaH); 

          // Invert tridiagonal matrix equation 
          Tridiagonal (a, b, c, w, u); 
          for (int i = 1; i <= N; i++) TT(i, j) = u(i); 
       } 
    } 
  else 
    { 
      for (int j = 1; j < J; j++) 
       {  
          Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2), u(N+2); 

          // Initialize tridiagonal matrix 
          for (int i = 2; i <= N; i++) a(i) = - 0.5 * C; 
          for (int i = 1; i <= N; i++) 
           b(i) = 1. + C * (1. + 0.5 * double (j * j) * kappa * kappa); 
          b(1) -= betaL / (alphaL * dx - betaL); 
          b(N) += betaH / (alphaH * dx + betaH); 
          for (int i = 1; i <= N-1; i++) c(i) = - 0.5 * C; 

          // Initialize right-hand side vector 
          for (int i = 1; i <= N; i++) 
           w(i) = V(i, j); 
          w(1) += 0.5 * C * GammaL(j) * dx / (alphaL * dx - betaL); 
          w(N) += 0.5 * C * GammaH(j) * dx / (alphaH * dx + betaH); 

          // Invert tridiagonal matrix equation 
          Tridiagonal (a, b, c, w, u); 
          for (int i = 1; i <= N; i++) TT(i, j) = u(i); 
       } 

      for (int i = 1; i <= N  ; i++)  
        { 
          TT(i, 0) = 0.; TT(i, J) = 0.; 
        } 
    } 

  // Reconstruct solution 
  for (int i = 1; i <= N; i++) 
    { 
      Array<double,1> In(J+1), Out(J+1); 

      for (int j = 0; j <= J; j++) In(j) = TT(i, j); 

      if (Neumann) 
        fft_backward_cos (In, Out); 
      else 
        fft_backward_sin (In, Out); 

      for (int j = 0; j <= J; j++) T(i, j) = Out(j); 
    } 

  // Calculate i=0 and i=N+1 values   
  for (int j = 0; j <= J; j++) 
    { 
      T(0, j) = (gammaL(j) * dx - betaL * T(1, j)) / 
        (alphaL * dx - betaL); 
      T(N+1, j) = (gammaH(j) * dx + betaH * T(N, j)) / 
        (alphaH * dx + betaH); 
    } 
}

Ví dụ giải phương trình khuếch tán 2 chiều

Bây giờ ta hãy giải phương trình khuếch tán theo 2 chiều bằng kĩ thuật sai phân hữu hạn đã thảo luận. Ta cần tìm nghiệm của PT ([diff2deq]) trong miền 0 ≤ x ≤ 10 ≤ y ≤ 1, tuân theo điều kiện biên sau lúc t = 0:

T(x,y,0) = 1     với |x – 0,5| < 0,1 và |y – 0,5| < 0,1
T(x,y,0) = 0     trong trường hợp còn lại

Các điều kiện biên đơn giản chỉ là T(0, y) = T(1, y) = T(x, 0) = T(x, 1) = 0.

temp

Hình. Khuếch tán theo 2 chiều. Việc tính toán số trị được thực hiện với D = 1 , δt = 10−4, và N = J = 128. Các đồ thị mật độ của T(x, y, t) tương ứng với các thời điểm: t = 0,000 (trái-trên), t = 0,001 (phải-trên), t = 0,002 (trái-hàng giữa), t = 0,003 (phải-hàng giữa), t = 0,004 (trái-dưới), and t = 0,005 (phải-dưới). [diff2d6]

Hình [diff2d6] cho thấy diễn biến của T(x, y, t) khi tính bằng chương trình để giải phương trình khuếch tán 2 chiều với D = 1 , δt = 10−4 , và N = J = 128.

Bài toán 3 chiều

Kĩ thuật được trình bày ở trên để giải phương trình khuếch tán theo 2 chiều tuân theo một lớp các điều kiện biên nhất định có thể dễ dàng được mở rộng cho trường hợp 3 chiều. Trong trường hợp 3 chiều, ta cần biến đổi Fourier theo hai phương (chẳng hạn yz) để rút gọn bài toán về các phương trình độc lập nhau biểu thị ma trận ba đường chéo. Những phương trình này có thể được nghịch đảo theo cách thông thường, sau đó nghiệm được tái lập bằng phép biến đổi Fourier ngược hai lần.

Advertisements

19 phản hồi

Filed under Vật lý tính toán

19 responses to “Chương 6: Phương trình khuếch tán

  1. Pingback: Chương 1: Giới thiệu chung | Blog của Chiến

  2. Kutelove

    Cảm ơn anh Quang Chiến với bài viết này,
    Hiện tại em cũng đang phải giải bài toán khuếch tán 1D, sử dụng phương pháp Crank – Nicolson. Em lập trình trên nền matlab và đang vấp phải một số khó khăn. Như trong bài viết này của anh thì các hằng số alpha, beta và gamma có ý nghĩa như thế nào?

    • α, β, γ là các biểu thức toán chứ chúng không hẳn có ý nghĩa vật lý riêng. Cho nên bạn lựa chọn để phù hợp với bài toán. Cách chọn là xem bài toán được xét dùng biên Dirichlet hay Neumann.

      Cụ thể, với điều kiện biên Dirichlet ta ấn định một time-series cho T, hay hàm của T theo thời gian. Như vậy T_dirichlet = T(t). Thế thì β = 0 và lúc này T_dirichlet = γ(t) / α(t). Nhiệm vụ của bạn là tự chọn γ và α. Có thể lấy α = 1 luôn cho tiện.

      Thứ hai là điều kiện biên Neumann, ở đó ta ấn định gradient của T. Vậy T_neumann = ∂T(t)/∂x. Khi đó α = 0, T_neumann = γ(t) / β(t) và bạn lại có thể cho β = 1 rồi tự chọn hàm γ.

      Có bài toán vật lý mà không bỏ được thành phần nào và ta có điều kiện biên hỗn hợp. Khi đó bạn cần tự định nghĩa cả α, β, γ. Chúng có thể được khai báo là function theo t.

      Lưu ý là code C++ của tác giả trong bài minh họa cho trường hợp đơn giản, α, β, γ đều là hằng số. Cho nên khi đọc sẽ khó hình dung.

      • Kutelove

        Cảm ơn những chia sẻ của anh,
        Em muốn biết cụ thể bài toán của em như sau:
        – Nồng độ ban đầu ở lớp ngoài cùng C(s), và giá trị này tăng theo thời gian
        – Chiều sao khuếch tán là L, tính trong khoảng thời gian t
        Vậy để tính nồng độ tại thớ thứ i, tại thời điểm t thì lựa chọn điều kiện biên thế nào cho hợp lý?
        EM xin cảm ơn!

      • Biên đặt ở lớp ngoài cùng, x = 0. Điều kiện biên là dạng Dirichlet, C(x=0) = f(t) = 0.01*tanh(t) chẳng hạn, để tiệm cận về 0,01. Ở biên này, α=1, β=0, γ=0.01*tanh(t). Bạn nói “chiều sâu khuếch tán” thì tôi không rõ lắm. Nhưng chắc chắn đây là dữ kiện giúp ta tính được hệ số khuếch tán, D. Còn biên phía đầu kia, không có ràng buộc gì về nồng độ cả, ta cho nó là biên Neumann với gradient bằng 0. Ở biên này, α=0, β=1, γ=0.

  3. Kutelove

    Em xin cảm ơn,
    Biên đặt ở mặt ngoài cùng của em với x=0 là C(x=0) = Cs – đây là một hàm đã biết theo thời gian. Hệ số khuếch tán D coi như không đổi trong suốt chiều dày L của vật liệu. Em cần tính được Nồng độ tại x = L. Vậy những điều kiện như tác giả nói có thể áp dụng được ko ?
    Em cảm ơn!

  4. Nếu như cần tính nồng độ tại x = L thì biên đầu kia phải đặt xa hẳn ra để khi áp dụng điều kiện Neumann mà không ảnh hưởng đến kết quả. Bạn cần thử viết mã lệnh để xem kết quả tính ra có hợp lý không. Kiểm tra xem đúng là đi từ nguồn ra xa thì nồng độ có giảm dần không.

    • Kutelove

      Hi,
      Em thử tạo lập thuật toán theo phương pháp Crank-Nicolson với phương trình như sau:
      A.C(t+1)=B*C(t)
      Với C(t) là nồng độ clorua ở các lớp từ 1-N tại thời điểm t
      C(t+1) là nồng độ clorua ở các lớp từ 1-N tại thời điểm t+1
      Vấn đề của em là nếu tính được C(1), tức là nồng độ clorua ở thời điểm đầu tiên của các lớp thì sẽ giải được bài toán thôi. Tuy nhiên em đang bị vướng điều kiện biên này. Mới chỉ xử lý được khi x = 0 thì C=Cs, còn khi t=0 thì chưa biết áp điều kiện thế nào? Vậy em phải làm thế nào?
      Còn giả sử như em biết tại thời điểm t nào đó, C(L,t) = Cth, vậy thì áp điều kiện biên thế nào cho hợp lý?
      Em cảm ơn!

      • t=0 ứng với điều kiện đầu. Khí đó C=0 cho tất cả các điểm, trừ điểm biên mà tại đó C đã cho. Còn cách kí hiệu của bạn chưa hợp lý. Nên làm như sau: Gọi C(i) là nồng độ tại điểm thứ i, i=0 đến N. Trong đó i=0 là biên mà ta cho trước nồng độ theo điều kiện Dirichlet (tức là lớp ngoài cùng). i=N là biên rất sâu phía trong mà ở đó giả sử là biên Neumann. Dùng một vòng lặp theo biến thời gian t. Ở mỗi thời điểm ta tính được C'(i) theo C(i) đã biết. C’ tức là C tại t+1. Bằng cách này sẽ tính được hết từ C'(1) đến C'(N-1). Còn C'(0) và C'(N) lại áp dụng điều kiện biên.

      • Kutelove

        Chào anh Chiến,
        EM cảm ơn những gợi ý của anh và em đã giải được bài toán 1D. Bây giờ là bài toán 2D và em cũng dự định giải theo phương pháp Crank – Nicholson. anh có thể cho em xin tài liệu chi tiết hơn về bài toán này được ko ạ?
        EM xin cảm ơn,

      • Em xem ở phần code đó. Trong bài này tác giả dùng sai phân ẩn. Dĩ nhiên để giải được theo cách của họ, phải hiểu được biến đổi Fourier. Tại sao 1 chiều ko cần biến đổi kiểu này, mà bây giờ 2 chiều lại cần? Hãy suy nghĩ trả lời câu hỏi này. (Xem lại chương Phương trình Poisson.) Đồng thời vẽ một lược đồ sai phân (discretization scheme, tra google để tìm hiểu thêm).
        Còn nếu em dùng cách sai phân hiện thì có khả năng giải đơn giản hơn. Có những mô hình toán trong kĩ thuật dùng sơ đồ hiện 2 chiều.
        Ngoài ra còn có cách giải 2 chiều theo kiểu “khử luân hướng” (ADI – Alternate Direction Implicit scheme). Hãy tự tìm đọc thêm. Các cách này không cần biến đổi Fourier.

  5. Yên

    Chào anh Chiến ! Em có bài toán này cần sự giúp đỡ của anh
    một vật có nhiệt độ t1 , môi trường có nhiệt độ t2 = hằng số , tìm phương trình truyền nhiệt khi t2 =t1, điều kiện biên, điều kiện đầu ??
    mong sự giúp đỡ của anh , em xin cảm ơn

    • Đầu tiên bạn cần cụ thể hóa bài toán với những giả thiết rồi mới áp dụng phương pháp tính trong bài viết để giải. Ví dụ mình sẽ đặt những giả thiết như sau:
      * Vật dẫn đồng chất, có kích thước bằng L. Nó chiếm không gian 0 <= x <= L.
      * Vật dẫn này có nhiệt độ ban đầu là To. Do nó cô lập và không được cấp nhiệt nên sẽ ấm lên/nguội dần đi để bằng với nhiệt độ môi trường
      * Môi trường có nhiệt độ Tm. Do môi trường lớn vô hạn (chiếm khoảng không gian L < x < ∞) nên nhiệt độ ở xa vật luôn không đổi (bằng Tm ở x = ∞).
      * Định luật Fourier về truyền nhiệt q = k ΔT/Δx (k khác nhau đối với vật dẫn và môi trường).
      * Tại chỗ mép vật dẫn (x = L), sự tỏa nhiện sẽ như thế nào (q = ?) Điều này quan trọng, từ đó bạn mới giải được bài toán.

  6. Ẩn danh

    cảm ơn anh vì bài viết hay

  7. Tuan Vu

    Rất cảm ơn a Chiến vì những chia sẻ , bản thân em đang học master tại Pháp và hiện cũng đang làm thực tập về cái Equation convection diffusion này . Hơi khác một chút là vì lĩnh vực nghiên cứu là tính thấm trong các môi trường lỗ rỗng (milieux porous ) vì vậy trong phương trình có sự xuất hiện thêm của 1 Term để biểu diễn sự lắng đọng lại của các hạt trong chất lỏng ( Deposition of particles ) , gây ra sự thay đổi trong tính thấm . Sử dụng matlab và Comsol để tính toán và mô phỏng chúng .
    Rất vui vì biết thêm 1 người cũng có hứng thú về vấn đề này

  8. Ẩn danh

    Anh có thể giúp đỡ em về cách mô phỏng nghiệm thu được khi giải bài toán truyền thiệt trong thanh hữu hạn và trong thanh có nguồn nhiệt đc không ạ? Em xin cảm ơn!!!

    • Mình không có lời giải sẵn. Bạn có thể tìm đọc cuốn sách của Chapra và Canale “Numerical Methods for Engineers”, trong đó có một chương chuyên về truyền nhiệt và ví dụ tính bằng MatLab.

Trả lời

Mời bạn điền thông tin vào ô dưới đây hoặc kích vào một biểu tượng để đăng nhập:

WordPress.com Logo

Bạn đang bình luận bằng tài khoản WordPress.com Log Out / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Log Out / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Log Out / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Log Out / Thay đổi )

Connecting to %s