Chương 8: Chương trình tính dùng phương pháp hạt-trong-ô (particle-in-cell)

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

Giới thiệu chung

Ta hãy xét một dòng plasma một chiều đều, phi từ tính, chứa N electron và N ion tích điện bằng đơn vị. Như vậy, ở những khoảng thời gian ngắn, có thể coi các ion như là nền trung tính đứng yên, và chỉ xét chuyển động của các electron. Đặt ri là tọa độ phương x của electron thứ i. Các phương trình chuyển động của electron thứ i này là:

\begin{array}{rcl}  \frac{dr_i}{dt} &=& v_i,\\[0.5ex]  \frac{dv_i}{dt} &=& - \frac{e\,E(r_i)}{m_e},  \end{array}

trong đó e > 0 là độ lớn của điện tích electron, me là khối lượng electron, và E(x) là thành phần theo phương x của cường độ điện trường tại vị trí x. Đến đây, có thể biểu diễn cường độ điện trường dưới dạng một điện thế:

E(x) =  − dφ(x) / dx. 

Ngoài ra, từ phương trình Poisson-Maxwell, ta có

\displaystyle{  \frac{d^2\phi(x)}{dx^2} = -\frac{e}{\epsilon_0}\,\{n_0- n(x)\}, }

trong đó ε0 là độ điện thẩm của môi trường, n(x) là mật độ electron (n(x) dx là số electron có trong khoảng từ x đến x + dx), và n0 là mật độ đồng nhất của ion. Dĩ nhiên, giá trị trung bình của n(x) đúng bằng n0, vì số ion và electron phải bằng nhau.

Ta hãy xét một hàm phân bố electron ban đầu bao gồm hai luồng Maxwell truyền ngược hướng nhau với vận tốc trung bình vb và độ lan tỏa nhiệt vth,

\displaystyle{  f(x,v) = \frac{n_0}{2}\left\{\frac{1}{\sqrt{2\,\pi}\,v_{th}}  \,{\rm e}^{-(v-v_b)^2/2\,v_{th}^{\,2}}+ \frac{1}{\sqrt{2\,\pi}\,v_{th}}\,{\rm e}^{-(v+v_b)^2/2\,v_{th}^{\,2}}  \right\}. }

Ở đây, f(x, v) dxdv là số electron trong khoảng giữa xx + dx có vận tốc trong khoảng từ v đến v + dv. Dĩ nhiên, n(x) = ∫−∞f(x, v) dv. Nhiệt độ của luồng, T, liên hệ với vận tốc nhiệt qua hệ thức v_{th} = \sqrt{k_B\,T/m_e}, trong đó kB là hằng số Boltzmann. Người ta đã phát hiện rằng nếu vb lớn hơn nhiều so với vth thì dạng phân bố trên sẽ bất ổn định theo tình trạng bất ổn plasma có tên bất ổn định hai dòng. 1 Ta hãy nghiên cứu sự bất ổn định này về mặt số trị.

Lược đồ chuẩn hóa

Sẽ tiện hơn khi ta chuẩn hóa thời gian theo ωp−1, trong đó

\displaystyle{ \omega_p^2 = \frac{n_0\,e^2}{\epsilon_0\,m_e} }

được gọi là tần số plasma, tức là tần số điển hình trong các dao động tĩnh điện của các electron. Tương tự, sẽ tiện hơn khi ta chuẩn hóa chiều dài theo đại lượng được gọi là chiều dài Debye:
λD = vth / ωp, 
vốn là cỡ độ dài mà vượt quá nó, các electron sẽ thể hiện hiệu ứng tụ tập, tức là hiệu ứng kiểu plasma, thay vì chỉ là các hạt riêng rẽ.

Các phương trình đã chuẩn hóa có dạng:

\begin{array}{rcl}  \frac{dx_i}{dt} &=&v_i,\label{eqm1}\\[0.5ex]  \frac{d v_i}{dt} &=& - E(x_i),\label{eqm2}\\[0.5ex]  E(x) &=& - \frac{d\phi(x)}{dx},\label{eqm4}\\[0.5ex]  \frac{d^2\phi(x)}{dx^2} &=& \frac{n(x)}{n_0}-1.\label{eqm3}  \end{array}

còn hàm phân bố ban đầu đang xét có dạng

f(x,v) = \frac{n_0}{2}\left\{\frac{1}{\sqrt{2\,\pi}}  \,{\rm e}^{-(v-v_b)^2/2}+ \frac{1}{\sqrt{2\,\pi}}\,{\rm e}^{-(v+v_b)^2/2}  \right\}.

Lưu ý rằng vth = 1 trong hệ đơn vị chuẩn hóa.

Ta hãy giải hệ phương trình trên trong miền 0 ≤ x ≤ L. Ngoài ra, để đơn giản, ta chọn các điều kiện biên tuần hoàn, cụ thể là các biên trái và phải của miền nghiệm. Điều này có nghĩa là n(0) = n(L), φ(0) = φ(L), và E(0) = E(L). Ngoài ra, mỗi electron khi vượt khỏi biên phải của miền nghiệm cần xuất hiện trở lại ở phía biên trái với cùng vận tốc, và ngược lại.

Nghiệm của hệ phương trình chuyển động cho electron

Ta có thể giải các phương trình chuyển động của electron, ([eqm1]) và ([eqm2]), dưới dạng một hệ gồm 2N PVT bậc nhất bằng các phương pháp RK4 đã đề cập đến ở Mục [ode]. Tuy nhiên, để tính được vế phải của các phương trình này ta cần biết điện trường E(x) ở mỗi bước thời gian. Để làm điều này ta có thể giải phương trình Poisson, ([eqm3]), ở mỗi bước thời gian. Tuy nhiên, muốn tính được số hạng nguồn cho phương trình này, ta cần tính mật độ electron, n(x), vốn dĩ là một hàm theo vị trí tức thời của các electron.

Tính mật độ electron

Để tính được mật độ electron n(x) từ các tọa độ của electron ri, ta áp dụng một phương pháp có tên là hạt-trong-ô (hay còn có thể gọi là phương pháp “chất điểm-trong-ô”, particle-in-cell method, PIC). Ta hãy định nghĩa một tập hợp J điểm nút lưới cách đều nhau ở các tọa độ

xj = jδx, 

với j = 0, J − 1, trong đó δx = L / J. Chọn nj ≡ n(xj). Giả sử rằng electron thứ i nằm giữa các nút lưới thứ j(j + 1), hay xj < ri < xj + 1. Ta đặt

\displaystyle{  n_j \rightarrow n_j+\left.\left(\frac{x_{j+1}-r_i}{x_{j+1}-x_j}\right)\right/\delta x, }

\displaystyle{  n_{j+1} \rightarrow n_{j+1}+ \left.\left(\frac{r_i-x_j}{x_{j+1}-x_j}\right)\right/\delta x. }

Như vậy là, njδx tăng thêm 1 nếu như electron ở nút lưới thứ j, nj + 1δx tăng thêm 1 nếu như electron ở nút lưới thứ (j + 1), và cả njδx lẫn nj + 1δx đều tăng thêm 1/2 nếu electron ở giữa hai nút lưới này, v.v. Bằng cách thực hiện một phép gán như vậy lần lượt cho mỗi electron, ta có thể lập nên các gía trị nj từ các tọa độ của electron (coi rằng tất cả nj ban đầu đều được cho bằng không).

Nghiệm của phương trình Poisson

Xét nghiệm của phương trình Poisson:

d2φ(x) / dx2 = ρ(x), 

trong đó ρ(x) = n(x)/n0 − 1. Lưu ý rằng n0 = N / L theo đơn vị chuẩn hóa. Đặt φj ≡ φ(xj)ρj ≡ ρ(xj). Ta có thể viết

\begin{array}{rcl}  \phi_j &=& \sum_{j'=0,J-1}\hat{\phi}_{j'} \,{\rm e}^{\,{\rm i}\,j\,j'\,2\,\pi/J},\\[0.5ex]  \rho_j &=& \sum_{j'=0,J-1}\hat{\rho}_{j'} \,{\rm e}^{\,{\rm i}\,j\,j'\,2\,\pi/J},  \end{array}

vốn tự thỏa mãn các điều kiện biên tuần hoàn φJ = φ0ρJ = ρ0. Lưu ý rằng ρ̂0 = 0, vì ∫ 0Ln(x) dx = n0. Các giá trị ρ̂j khác thu được từ

\hat{\rho}_j = \frac{1}{J}\sum_{j'=0,J-1}\rho_{j'} \,{\rm e}^{-{\rm i}\,j\,j'\,2\,\pi/J},

với j = 1, J − 1. Dạng biến đổi Fourier của phương trình Poisson cho ta

φ̂0 = 0

\displaystyle{  \hat{\phi}_j = - \frac{\hat{\rho}_j}{j^2\,\kappa^2} }

với j = 1, J / 2, trong đó κ = 2π / L. Sau cùng,

φ̂j = φ̂ * J − j

với j = J / 2 + 1 đến J − 1; điều này đảm bảo rằng φj là số thực. Dạng rời rạc của PT ([eqm4]) là

Ej = (φj − 1 − φj + 1) / (2 δx). 

Dĩ nhiên, j = 0j = J − 1 là các trường hợp đặc biệt mà có thể giải được từ các điều kiện biên tuần hoàn. Cuối cùng, giả sử rằng tọa độ của electron thứ i nằm giữa các điểm nút thứ j(j + 1), tức là xj < ri < xj + 1. Ta có thể nội suy tuyến tính để tìm điện trường chiếu theo electron thứ i này:

\displaystyle{  E(r_i) = \left(\frac{x_{j+1}-r_i}{x_{j+1}-x_j}\right)E_j+\left(\frac{r_i - x_j}{x_{j+1}-x_j}\right)  E_{j+1}. }

Ví dụ chương trình PIC 1 chiều

Chương trình sau thực hiện các ý tưởng đã trình bày ở trên.

Hàm chính của chương trình đọc vào các tham số tính toán, kiểm tra xem chúng có nghĩa không, đặt tọa độ ban đầu của các electron, rồi tính theo các phương trình chuyển động của electron từ t = 0 đến một giá trị tmax cho trước, bằng một chương trình con RK4 với bước tính cố định cho trước, δt. Thông tin về các tọa độ pha-không gian của electron cùng điện trường được ghi ra các file một cách thường xuyên.

// 1-d PIC code to solve plasma two-stream instability problem.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <blitz/array.h>
#include <fftw.h>

using namespace blitz;

void Output (char* fn1, char* fn2, double t, 
         Array<double,1> r, Array<double,1> v);
void Density (Array<double,1> r, Array<double,1>& n);
void Electric (Array<double,1> phi, Array<double,1>& E);
void Poisson1D (Array<double,1>& u, Array<double,1> v, double kappa);
void rk4_fixed (double& x, Array<double,1>& y, 
              void (*rhs_eval)(double, Array<double,1>, Array<double,1>&), 
              double h);
void rhs_eval (double t, Array<double,1> y, Array<double,1>& dydt);
void Load (Array<double,1> r, Array<double,1> v, Array<double,1>& y);
void UnLoad (Array<double,1> y, Array<double,1>& r, Array<double,1>& v);
double distribution (double vb);

double L; int N, J;

int main()
{
  // Parameters
  L;            // Domain of solution 0 <= x <= L (in Debye lengths)
  N;            // Number of electrons
  J;            // Number of grid points
  double vb;    // Beam velocity
  double dt;    // Time-step (in inverse plasma frequencies)
  double tmax;  // Simulation run from t = 0. to t = tmax

  // Get parameters
  printf ("Please input N:  "); scanf ("%d", &N);
  printf ("Please input vb:  "); scanf ("%lf", &vb); 
  printf ("Please input L:  "); scanf ("%lf", &L);
  printf ("Please input J:  "); scanf ("%d", &J);
  printf ("Please input dt:  "); scanf ("%lf", &dt); 
  printf ("Please input tmax:  "); scanf ("%lf", &tmax);
  int skip = int (tmax / dt) / 10;
  if ((N < 1) || (J < 2) || (L <= 0.) || (vb <= 0.) 
      || (dt <= 0.) || (tmax <= 0.) || (skip < 1))
    {
      printf ("Error - invalid input parameters\n");
      exit (1);
    }

  // Set names of output files
  char* phase[11]; char* data[11];
  phase[0] = "phase0.out";phase[1] = "phase1.out";phase[2] = "phase2.out"; 
  phase[3] = "phase3.out";phase[4] = "phase4.out";phase[5] = "phase5.out"; 
  phase[6] = "phase6.out";phase[7] = "phase7.out";phase[8] = "phase8.out";
  phase[9] = "phase9.out";phase[10] = "phase10.out";data[0] = "data0.out";  
  data[1] = "data1.out"; data[2] = "data2.out"; data[3] = "data3.out"; 
  data[4] = "data4.out"; data[5] = "data5.out"; data[6] = "data6.out"; 
  data[7] = "data7.out"; data[8] = "data8.out"; data[9] = "data9.out"; 
  data[10] = "data10.out";

  // Initialize solution
  double t = 0.;
  int seed = time (NULL); srand (seed);
  Array<double,1> r(N), v(N);
  for (int i = 0; i < N; i++)
    {
      r(i) = L * double (rand ()) / double (RAND_MAX);
      v(i) = distribution (vb);
    }
  Output (phase[0], data[0], t, r, v);

  // Evolve solution
  Array<double,1> y(2*N);
  Load (r, v, y);
  for (int k = 1; k <= 10; k++)
    {
      for (int kk = 0; kk < skip; kk++)
        {
           // Take time-step
           rk4_fixed (t, y, rhs_eval, dt);

           // Make sure all coordinates in range 0 to L.
           for (int i = 0; i < N; i++)
             {
               if (y(i) < 0.) y(i) += L;
               if (y(i) > L) y(i) -= L;
             }

           printf ("t = %11.4e\n", t);
        }
      printf ("Plot %3d\n", k);

      // Output data
      UnLoad (y, r, v);
      Output(phase[k], data[k], t, r, v);
    }

  return 0;
}

Đoạn chương trình dưới đây ghi dữ liệu kết quả mô phỏng ra nhiều file dữ liệu.

// Write data to output files

void Output (char* fn1, char* fn2, double t,
         Array<double,1> r, Array<double,1> v)
{  
  // Write phase-space data
  FILE* file = fopen (fn1, "w");
  for (int i = 0; i < N; i++)
    fprintf (file, "%e %e\n", r(i), v(i));
  fclose (file);

  // Write electric field data
  Array<double,1> ne(J), n(J), phi(J), E(J);  
  Density (r, ne);
  for (int j = 0; j < J; j++)
    n(j) = double (J) * ne(j) / double (N) - 1.;
  double kappa = 2. * M_PI / L; 
  Poisson1D (phi, n, kappa);
  Electric (phi, E);

  file = fopen (fn2, "w");
  for (int j = 0; j < J; j++)
    {
      double x = double (j) * L / double (J);
      fprintf (file, "%e %e %e %e\n", x, ne(j), n(j), E(j));
    }
  double x = L;
  fprintf (file, "%e %e %e %e\n", x, ne(0), n(0), E(0));
  fclose (file);
}

Đoạn chương trình dưới đây trả lại một vận tốc ngẫu nhiên tuân theo hàm phân bố Maxwell kép trong trường hợp hai luồng truyền ngược nhau. Thuật toán đã dùng ở đây được gọi là phương pháp bác bỏ, và sẽ được đề cập ở phần sau của khóa học.

// Function to distribute electron velocities randomly so as 
// to generate two counter propagating warm beams of thermal
// velocities unity and mean velocities +/- vb.
// Uses rejection method.

double distribution (double vb)
{ 
  // Initialize random number generator
  static int flag = 0;
  if (flag == 0)
    {
      int seed = time (NULL);
      srand (seed);
      flag = 1;
    }

  // Generate random v value
  double fmax = 0.5 * (1. + exp (-2. * vb * vb));
  double vmin = - 5. * vb;
  double vmax = + 5. * vb;
  double v = vmin + (vmax - vmin) * double (rand ()) / double (RAND_MAX);

  // Accept/reject value
  double f = 0.5 * (exp (-(v - vb) * (v - vb) / 2.) +
            exp (-(v + vb) * (v + vb) / 2.));
  double x = fmax * double (rand ()) / double (RAND_MAX);
  if (x > f) return distribution (vb);
  else return v;
}

Đoạn chương trình dưới ước tính mật độ electron trên một lưới giãn cách đều, cho trước các tọa độ tức thời của electron.

// Evaluates electron number density n(0:J-1) from 
// array r(0:N-1) of electron coordinates.

void Density (Array<double,1> r, Array<double,1>& n)
{
  // Initialize 
  double dx = L / double (J);
  n = 0.;

  // Evaluate number density.
  for (int i = 0; i < N; i++)
    {
      int j = int (r(i) / dx);
      double y = r(i) / dx - double (j);
      n(j) += (1. - y) / dx;
      if (j+1 == J) n(0) += y / dx;
      else n(j+1) += y / dx;
    }
}

Các hàm sau để bọc thư viện fftw với các hàm tuần hoàn.

// Functions to calculate Fourier transforms of real data 
// using fftw Fast-Fourier-Transform routine.
// Input/ouput arrays are assumed to be of extent J.

// Calculates Fourier transform of array f in arrays Fr and Fi
void fft_forward (Array<double,1>f, Array<double,1>&Fr, 
      Array<double,1>& Fi)
{
  fftw_complex ff[J], FF[J];

  // Load data
  for (int j = 0; j < J; j++)
    {
      c_re (ff[j]) = f(j); c_im (ff[j]) = 0.;
    }

  // Call fftw routine
  fftw_plan p = fftw_create_plan (J, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_one (p, ff, FF);
  fftw_destroy_plan (p); 

  // Unload data
  for (int j = 0; j < J; j++)
    {
      Fr(j) = c_re (FF[j]); Fi(j) = c_im (FF[j]);
    }

  // Normalize data
  Fr /= double (J);
  Fi /= double (J);
}

// Calculates inverse Fourier transform of arrays Fr and Fi in array f
void fft_backward (Array<double,1> Fr, Array<double,1> Fi, 
      Array<double,1>& f)
{
  fftw_complex ff[J], FF[J];

  // Load data
  for (int j = 0; j < J; j++)
    {
      c_re (FF[j]) = Fr(j); c_im (FF[j]) = Fi(j);
    }

  // Call fftw routine
  fftw_plan p = fftw_create_plan (J, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_one (p, FF, ff);
  fftw_destroy_plan (p); 

  // Unload data
  for (int j = 0; j < J; j++)
      f(j) = c_re (ff[j]); 
}

Đoạn chương trình sau giải phương trình Poisson 1 chiều để tìm điện thế tức thời trên một lưới đều.

// Solves 1-d Poisson equation:
//    d^u / dx^2 = v   for  0 <= x <= L
// Periodic boundary conditions:
//    u(x + L) = u(x),  v(x + L) = v(x)
// Arrays u and v assumed to be of length J.
// Now, jth grid point corresponds to
//    x_j = j dx  for j = 0,J-1
// where dx = L / J.
// Also,
//    kappa = 2 pi / L

void Poisson1D (Array<double,1>& u, Array<double,1> v, double kappa)
{
  // Declare local arrays.
  Array<double,1> Vr(J), Vi(J), Ur(J), Ui(J);

  // Fourier transform source term
  fft_forward (v, Vr, Vi);

  // Calculate Fourier transform of u
  Ur(0) = Ui(0) = 0.;
  for (int j = 1; j <= J/2; j++)
    {
      Ur(j) = - Vr(j) / double (j * j) / kappa / kappa;
      Ui(j) = - Vi(j) / double (j * j) / kappa / kappa;
    } 
  for (int j = J/2; j < J; j++)
    {
      Ur(j) = Ur(J-j);
      Ui(j) = - Ui(J-j);
    }

  // Inverse Fourier transform to obtain u
  fft_backward (Ur, Ui, u);
}

Hàm sau đây để tính điện trường trên một lưới đều, từ điện thế.

// Calculate electric field from potential

void Electric (Array<double,1> phi, Array<double,1>& E)
{
  double dx = L / double (J);

  for (int j = 1; j < J-1; j++)
    E(j) = (phi(j-1) - phi(j+1)) / 2. / dx;
  E(0) = (phi(J-1) - phi(1)) / 2. / dx;
  E(J-1) = (phi(J-2) - phi(0)) / 2. / dx;
}

Đoạn chương trình sau tính vế phải trong các phương trình chuyển động của electron. Nó được viết để dùng với chương trình giải RK4 bước cố định đã được đề cập đến từ phần trước của khóa học.

// Electron equations of motion:
//    y(0:N-1)  = r_i
//    y(N:2N-1) = dr_i/dt

void rhs_eval (double t, Array<double,1> y, Array<double,1>& dydt)
{
  // Declare local arrays
  Array<double,1> r(N), v(N), rdot(N), vdot(N), r0(N);
  Array<double,1> ne(J), rho(J), phi(J), E(J);

  // Unload data from y
  UnLoad (y, r, v);

  // Make sure all coordinates in range 0 to L
  r0 = r;
  for (int i = 0; i < N; i++)
    {
      if (r0(i) < 0.) r0(i) += L;
      if (r0(i) > L) r0(i) -= L;
    }

  // Calculate electron number density
  Density (r0, ne);

  // Solve Poisson's equation
  double n0 = double (N) / L;
  for (int j = 0; j < J; j++)
    rho(j) = ne(j) / n0 - 1.;
  double kappa = 2. * M_PI / L; 
  Poisson1D (phi, rho, kappa);

  // Calculate electric field
  Electric (phi, E);

  // Equations of motion
  for (int i = 0; i < N; i++)
    {
      double dx = L / double (J);
      int j = int (r0(i) / dx);
      double y = r0(i) / dx - double (j);

      double Efield;
      if (j+1 == J)
         Efield = E(j) * (1. - y) + E(0) * y;
      else
         Efield = E(j) * (1. - y) + E(j+1) * y;

      rdot(i) = v(i);
      vdot(i) = - Efield;
    }

  // Load data into dydt
  Load (rdot, vdot, dydt);
}

Các hàm sau đây nạp vào và gỡ các tọa độ pha-không gian của electron từ véc-tơ nghiệm y được dùng cho đoạn chương trình RK4.

// Load particle coordinates into solution vector

void Load (Array<double,1> r, Array<double,1> v, Array<double,1>& y)
{
  for (int i = 0; i < N; i++)
    {
      y(i) = r(i);
      y(N+i) = v(i);
    }
}

// Unload particle coordinates from solution vector

void UnLoad (Array<double,1> y, Array<double,1>& r, Array<double,1>& v)
{
  for (int i = 0; i < N; i++)
    {
      r(i) = y(i);
      v(i) = y(N+i);
    }
}

Phân bố pha-không gian của electron được ước tính tại các thời điểm khác nhau trong một phép tính 1 chiều bất ổn định của 2 dòng với N = 20000, J = 1000, L = 100, vb = 3, và δt = 0,1.[pic1]

Kết quả

Các Hình [pic1] và [pic2] cho thấy những phân bố pha-không gian của electron được xác định tại những thời điểm cách đều nhau trong phép tính bất ổn định hai dòng với 2 × 104 electron. Có thể thấy rằng sự phân bố ban đầu có dạng hai dải đều, ứng với hai luồng electron đối nhau. Tuy nhiên, theo thời gian, hai dải hình thành một cách tự phát những cấu trúc giãn nở ra và cuối cùng biến đổi cả dạng phân bố pha-không gian thành một tập hợp các xoáy cuộn nối với nhau. Ở trạng thái cuối cùng này, về cơ bản các electron đều nảy qua nảy lại trong một điện trường giả-tuần toàn được phát sinh từ sự bất đồng đều của mật độ electron. Nói cách khác, sự bất ổn định đã phá hủy một cách hiệu quả hai luồng. Vì vậy, sự bất ổn định hai luồng là vấn đề chính được xem xét đến trong các máy gia tốc hạt, vốn có các luồng hạt mang điện được phun đối đầu nhau.

Sự phân bố pha-không gian của electron được ước tính tại các thời điểm khác nhau trong bài toán bất ổn định hai luồng 1 chiều, với N = 20000, J = 1000, L = 100, vb = 3, và δt = 0,1.[pic2]

Thảo luận

Hiển nhiên là những ý kiến được bàn luận ở trên có thể được khái quát hóa một cách khá dễ dàng, để tính với các phân bố hạt mang điện tích trong hai và ba chiều. Các chương trình PIC có ưu điểm là tương đối dễ viết. Không may là những chương trình kiểu này thường chịu nhiễu động thống kê rõ rệt, vì chúng nói chung chỉ giải được với một số tương đối ít các hạt (điển hình là  ≤ 106). Những hệ vật lý thực tế không có biểu hiện giống sự nhiễu động thống kê như vậy, vì chúng thường cấu thành từ số lớn (cỡ số Avogadro  ∼ 1024) các hạt tương tác với nhau. Một vấn đề khác đối với các chương trình PIC là chúng không xử lý tốt sự va chạm giữa các hạt mang điện. Lý do của điều này là ở chỗ thường có rất nhiều hạt trong mỗi ô (để hợp với thực tế), và tầm ngắn của trường Coulomb từ những hạt này có xu hướng khử lẫn nhau (nhớ lại là điện trường chỉ được tính ở đỉnh, hay biên giới của ô).


  1. T.H. Stix, The theory of plasma waves, 1st Ed. (McGraw-Hill, New York NY, 1962).
Advertisements

%(count) bình luận

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

One response to “Chương 8: Chương trình tính dùng phương pháp hạt-trong-ô (particle-in-cell)

  1. Pingback: Nhập môn Vật lý tính toán | Blog của Chiến

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 Đăng xuất / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Đăng xuất / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Đăng xuất / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Đăng xuất / Thay đổi )

Connecting to %s