Реализация метода FDTD на PascalABC.NET: от уравнений Максвелла к визуализации волн

7 Ноя в 13:38 Прокомментировать Просмотры: 1

В этой статье мы смоделируем распространение электромагнитных волн, численно решив систему дифференциальных уравнений Максвелла в частных производных.

К задаче мы подойдём основательнее, чем в предыдущих публикациях. Я учёл обоснованные замечания читателей о выборе численного метода, проверил обоснованность схемы и правдоподобие результатов, а также описал условия, необходимые для воспроизведения расчёта.

В системе СИ уравнения Максвелла имеют следующий вид:

Реализация метода FDTD на PascalABC.NET: от уравнений Максвелла к визуализации волн

Первое уравнение — закон электромагнитной индукции Фарадея. Его смысл можно сформулировать так: изменение магнитной индукции вызывает вихревое электрическое поле.

Второе уравнение — теорема о циркуляции магнитного поля. Она утверждает, что электрический ток и изменение электрической индукции порождают вихревое магнитное поле.

Третье уравнение — закон Гаусса: электрический заряд служит источником электрической индукции.

Четвёртое уравнение — закон Гаусса для магнитного поля. Он отражает тот факт, что магнитные заряды не обнаружены.

В уравнениях используются следующие обозначения:

ρ — объёмная плотность стороннего электрического заряда (в системе СИ измеряется в Кл/м³);

j — плотность электрического тока (плотность тока проводимости; в СИ — А/м²);

c — скорость света в вакууме, равная 299 792 458 м/с;

E — напряжённость электрического поля (в СИ — В/м);

H — напряжённость магнитного поля (в СИ — А/м);

D — электрическая индукция (в СИ — Кл/м²);

B — магнитная индукция (в СИ — Тл = Вб/м² = кг•с⁻²•А⁻¹);

∇ — дифференциальный оператор набла.

Источник: [1]

Систему также нужно дополнить материальными уравнениями. Для простейшего случая — однородной изотропной среды — они записываются так:

$\mathbf{D} = \varepsilon \mathbf{E}, \quad \mathbf{B} = \mu \mathbf{H}.$

Здесь коэффициенты

$\varepsilon = \varepsilon_0 \varepsilon_r, \quad \mu = \mu_0 \mu_r.$

(диэлектрическая и магнитная проницаемости) считаются постоянными.

Частный случай этих уравнений — электростатика, при которой электрическое поле не зависит от времени.

Другой частный случай — магнитостатика, когда магнитное поле остаётся неизменным во времени, например, в поле постоянного магнита:

магнитное поле постоянного цилиндрического магнита. Неточная модель. — Магнитное поле постоянного цилиндрического магнита. Приближённая модель.

В этой работе мы рассмотрим динамические электромагнитные процессы, возбуждаемые гауссовым сигналом.

Для расчётов воспользуемся методом FDTD — методом конечных разностей во временной области.

Для моделирования нужны конечно-разностные формулы. Они относятся к классическим и приводятся в научной литературе (https://elibrary.ru/download/elibrary_20799764_92247445.pdf и https://elibrary.ru/vvfyhz):

Для компоненты электрического поля Ez:

$E_z^{n+1}(i,j,k) = E_z^n(i,j,k) + \frac{\Delta t}{\varepsilon} \left[ \frac{H_y^{n+1/2}(i+1,j,k) - H_y^{n+1/2}(i,j,k)}{\Delta x} - \frac{H_x^{n+1/2}(i,j+1,k) - H_x^{n+1/2}(i,j,k)}{\Delta y} \right].$

Для компоненты магнитного поля Hx:

$H_x^{n+1/2}(i,j,k) = H_x^{n-1/2}(i,j,k) - \frac{\Delta t}{\mu} \left[ \frac{E_z^n(i,j,k+1) - E_z^n(i,j,k)}{\Delta z} - \frac{E_y^n(i,j,k+1) - E_y^n(i,j,k)}{\Delta y} \right].$

В качестве источника электромагнитного возмущения зададим классический гауссов импульс:

$S(t) = \exp\!\left(-\frac{(t - t_0)^2}{\sigma_t^2}\right).$

В дискретном времени выражение принимает вид:

$S_n = \exp\!\left(-\frac{(n\,\Delta t - t_0)^2}{\sigma_t^2}\right).$

Уточним, что в программе источник задаётся как плотность тока, изменяющаяся по гауссовому закону.

Параметры импульса: J0 = 10^6 А/м², t0 = 30 шагов и σt = 10 шагов. Источник действует в одной ячейке в центре расчётной области. Возбуждается ячейка Ex с индексами (50, 49, 49), а не строго центральная: так получается из-за смещения компонент на сетке Йи.

Устойчивость схемы обеспечивается выполнением условия Куранта:

$\Delta t \le \frac{1}{c \sqrt{\left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2}\right)}},$

Временной шаг выбран на 5% меньше предельного значения по условию Куранта; число Куранта составляет примерно 0,55.

Чтобы дополнительно проверить корректность работы схемы, через каждые N шагов вычислялась суммарная энергия электромагнитного поля:

$\begin{equation} V(t) = \sum_{i,j,k} \left( |\mathbf{E}_{i,j,k}(t)|^2 + |\mathbf{H}_{i,j,k}(t)|^2 \right), \end{equation}$

Если V(t) не растёт экспоненциально, это свидетельствует о том, что численная ошибка не накапливается и схема остаётся устойчивой в течение всего расчёта.

Кроме того, функция MaxDivergenceE определяет максимальное значение |div E|. Вне области источника дискретная дивергенция E должна быть близка к машинному нулю, тогда как вблизи источника она отлична от нуля.

Граничные условия задаются как PEC на трёх гранях и PMC на остальных трёх: на границах выполняется условие E_tang = 0.

Речь идёт о трёх «нижних» гранях: x = 0, y = 0 и z = 0.

На верхних гранях (x = Nx−1 и т. д.) в массивах нет тангенциальных компонент E, поэтому там задаётся «магнитная» стенка, то есть H = 0.

В идеальном случае такие граничные условия не приводят к потерям энергии.

Расчёт проводился в кубической области размером 100×100×100 мм при шаге сетки Δx = Δy = Δz = 1 мм. Выполнено Nt = 300 временных шагов, что соответствует примерно 550 пс. За это время волновой фронт достигает границ области и отражается от них. Каждые пять шагов визуализировался срез Ez в плоскости z = L/2.

В этом случае обезразмеривание не требуется. Оно становится необходимым, когда величины выходят за удобный для машинной арифметики диапазон: вызывают переполнение или становятся меньше единицы настолько, что теряется точность.

Вычислительную схему реализуем на PascalABC.NET — представителе старого, но по-прежнему удобного семейства Pascal. В Python и C++ у меня возникли проблемы с индексацией, которые не удалось устранить. Ниже приведена версия программы после отладки и исправлений.

program FDTD3D_Animated_PEC_Energy;

uses
  GraphABC, System;

type
  Real = double;

const
  Nx = 100;
  Ny = 100;
  Nz = 100;

  dx = 1.0e-3;
  dy = 1.0e-3;
  dz = 1.0e-3;
  c  = 2.998e8;
  eps0 = 8.854e-12;
  mu0  = 4.0 * Pi * 1.0e-7;

  steps = 300;
  snap_every = 5;

  t0 = 30.0;
  sigma_t = 10.0;

  J0 = 1.0e6;

  // ---------- Раскладка окна ----------
  WinWidth  = 1240;
  WinHeight = 660;
  AnimX     = 10;
  AnimY     = 40;
  AnimSize  = 600;
  PlotX     = 630;
  PlotY     = 40;
  PlotW     = 590;
  PlotH     = 600;

var
  dt: Real;
  Ex, Ey, Ez, Hx, Hy, Hz: array of Real;
  eps_r, mu_r: array of Real;
  src_i, src_j, src_k, slice_z: integer;
  W, H: integer;
  E_hist: array of Real;   // история полной энергии

// ---------- Индексаторы Yee-сетки ----------
function idxEx(i, j, k: integer): integer;
begin Result := (i * (Ny - 1) + j) * (Nz - 1) + k; end;
function idxEy(i, j, k: integer): integer;
begin Result := (i * Ny + j) * (Nz - 1) + k; end;
function idxEz(i, j, k: integer): integer;
begin Result := (i * (Ny - 1) + j) * Nz + k; end;
function idxHx(i, j, k: integer): integer;
begin Result := (i * Ny + j) * Nz + k; end;
function idxHy(i, j, k: integer): integer;
begin Result := (i * (Ny - 1) + j) * Nz + k; end;
function idxHz(i, j, k: integer): integer;
begin Result := (i * Ny + j) * (Nz - 1) + k; end;
function idxEps(i, j, k: integer): integer;
begin Result := (i * Ny + j) * Nz + k; end;

// ---------- Источник ----------
function source(n: integer): Real;
begin
  Result := Exp(-Sqr((n - t0) / sigma_t));
end;

// ---------- PEC: E_tang = 0 на нижних гранях ----------
procedure ApplyPEC;
var i, j, k: integer;
begin
  for j := 0 to Ny - 1 do
    for k := 0 to Nz - 2 do
      Ey[idxEy(0, j, k)] := 0;
  for j := 0 to Ny - 2 do
    for k := 0 to Nz - 1 do
      Ez[idxEz(0, j, k)] := 0;

  for i := 0 to Nx - 1 do
    for k := 0 to Nz - 2 do
      Ex[idxEx(i, 0, k)] := 0;
  for i := 0 to Nx - 2 do
    for k := 0 to Nz - 1 do
      Ez[idxEz(i, 0, k)] := 0;

  for i := 0 to Nx - 1 do
    for j := 0 to Ny - 2 do
      Ex[idxEx(i, j, 0)] := 0;
  for i := 0 to Nx - 2 do
    for j := 0 to Ny - 1 do
      Ey[idxEy(i, j, 0)] := 0;
end;

// ---------- Полная энергия ----------
function TotalEnergy: Real;
var i, j, k: integer; sum: Real;
begin
  sum := 0;
  for i := 0 to Nx - 1 do
    for j := 0 to Ny - 2 do
      for k := 0 to Nz - 2 do
        sum += eps0 * Sqr(Ex[idxEx(i,j,k)]);
  for i := 0 to Nx - 2 do
    for j := 0 to Ny - 1 do
      for k := 0 to Nz - 2 do
        sum += eps0 * Sqr(Ey[idxEy(i,j,k)]);
  for i := 0 to Nx - 2 do
    for j := 0 to Ny - 2 do
      for k := 0 to Nz - 1 do
        sum += eps0 * Sqr(Ez[idxEz(i,j,k)]);
  for i := 0 to Nx - 2 do
    for j := 0 to Ny - 1 do
      for k := 0 to Nz - 1 do
        sum += mu0 * Sqr(Hx[idxHx(i,j,k)]);
  for i := 0 to Nx - 1 do
    for j := 0 to Ny - 2 do
      for k := 0 to Nz - 1 do
        sum += mu0 * Sqr(Hy[idxHy(i,j,k)]);
  for i := 0 to Nx - 1 do
    for j := 0 to Ny - 1 do
      for k := 0 to Nz - 2 do
        sum += mu0 * Sqr(Hz[idxHz(i,j,k)]);
  Result := 0.5 * sum * dx * dy * dz;
end;

// ---------- max |div E| ----------
function MaxDivergenceE: Real;
var i, j, k: integer; d, maxd: Real;
begin
  maxd := 0;
  for i := 1 to Nx - 2 do
    for j := 1 to Ny - 2 do
      for k := 1 to Nz - 2 do
      begin
        d := (Ex[idxEx(i+1,j,k)] - Ex[idxEx(i,j,k)]) / dx
           + (Ey[idxEy(i,j+1,k)] - Ey[idxEy(i,j,k)]) / dy
           + (Ez[idxEz(i,j,k+1)] - Ez[idxEz(i,j,k)]) / dz;
        if Abs(d) > maxd then maxd := Abs(d);
      end;
  Result := maxd;
end;

// ---------- Срез Ez (рисует только свою панель) ----------
procedure DrawSlice(frameArr: array of Real; x0, y0, sz: integer);
var
  i, j, cell: integer;
  val, norm: Real;
  col: Color;
begin
  norm := 0.0;
  for i := 0 to W * H - 1 do
    if Abs(frameArr[i]) > norm then norm := Abs(frameArr[i]);
  if norm = 0 then norm := 1.0;

  cell := sz div W;
  if cell < 1 then cell := 1;

  for i := 0 to W - 1 do
    for j := 0 to H - 1 do
    begin
      val := frameArr[i * H + j] / norm;
      if val < 0 then
        col := RGB(0, 0, 255 + Round(255 * val))
      else
        col := RGB(255, 255 - Round(255 * val), 0);
      Brush.Color := col;
      FillRect(x0 + i * cell, y0 + j * cell,
               x0 + (i + 1) * cell, y0 + (j + 1) * cell);
    end;

  Font.Color := clBlack;
  Font.Size := 11;
  TextOut(x0 + 5, y0 - 22, 'Срез Ez, z = ' + IntToStr(slice_z));
end;

// ---------- График энергии W(t) ----------
procedure DrawEnergyPlot(E_hist: array of Real; nCur, x0, y0, w, h: integer);
var
  i, px, py: integer;
  Emax, frac: Real;
  marginL, marginB, marginT, marginR, plotW, plotH: integer;
begin
  marginL := 60;
  marginB := 40;
  marginT := 25;
  marginR := 20;
  plotW := w - marginL - marginR;
  plotH := h - marginT - marginB;

  Emax := 1e-30;
  for i := 0 to nCur do
    if E_hist[i] > Emax then Emax := E_hist[i];
  if Emax <= 0 then Emax := 1;

  // Оси
  Pen.Width := 1;
  Pen.Color := clBlack;
  Line(x0 + marginL, y0 + marginT, x0 + marginL, y0 + marginT + plotH);
  Line(x0 + marginL, y0 + marginT + plotH,
       x0 + marginL + plotW, y0 + marginT + plotH);

  Font.Color := clBlack;
  Font.Size := 9;

  // Y-метки: 0%, 25%, 50%, 75%, 100%
  for i := 0 to 4 do
  begin
    frac := i / 4;
    py := y0 + marginT + plotH - Round(plotH * frac);
    Line(x0 + marginL - 5, py, x0 + marginL, py);
    TextOut(x0 + 10, py - 7, IntToStr(Round(frac * 100)) + '%');
  end;

  // X-метки
  for i := 0 to 4 do
  begin
    px := x0 + marginL + Round(plotW * i / 4);
    Line(px, y0 + marginT + plotH, px, y0 + marginT + plotH + 5);
    TextOut(px - 12, y0 + marginT + plotH + 8,
            IntToStr(Round((steps - 1) * i / 4)));
  end;

  // Кривая
  Pen.Width := 2;
  Pen.Color := clRed;
  for i := 0 to nCur do
  begin
    px := x0 + marginL + Round(plotW * i / (steps - 1));
    py := y0 + marginT + plotH - Round(plotH * E_hist[i] / Emax);
    if i = 0 then MoveTo(px, py) else LineTo(px, py);
  end;

  // «Бегунок» текущего шага
  if nCur > 0 then
  begin
    Pen.Width := 1;
    Pen.Color := clBlue;
    px := x0 + marginL + Round(plotW * nCur / (steps - 1));
    Line(px, y0 + marginT, px, y0 + marginT + plotH);
  end;

  // Заголовок и абсолютный максимум
  Font.Color := clBlack;
  Font.Size := 11;
  TextOut(x0 + marginL, y0 + 2, 'Полная энергия W(t)');
  Font.Size := 9;
  TextOut(x0 + marginL + 5, y0 + marginT + 5, Format('W_max = {0:E2}', Emax));
  TextOut(x0 + marginL + plotW div 2 - 25, y0 + h - 15, 'номер шага');
end;

// ======================== ОСНОВНАЯ ПРОГРАММА ========================
begin
  SetWindowSize(WinWidth, WinHeight);
  Window.Title := 'FDTD 3D PEC: срез Ez и энергия';

  dt := 0.95 * dx / (Sqrt(3.0) * c);
  WriteLn('CFL: dt=", dt, "  c*dt/dx = ', c*dt/dx,
          '  (предел 1/sqrt(3) = 0.5774)');

  src_i := Nx div 2;
  src_j := Ny div 2;
  src_k := Nz div 2;
  slice_z := Nz div 2;

  SetLength(Ex, Nx * (Ny - 1) * (Nz - 1));
  SetLength(Ey, (Nx - 1) * Ny * (Nz - 1));
  SetLength(Ez, (Nx - 1) * (Ny - 1) * Nz);
  SetLength(Hx, (Nx - 1) * Ny * Nz);
  SetLength(Hy, Nx * (Ny - 1) * Nz);
  SetLength(Hz, Nx * Ny * (Nz - 1));
  SetLength(eps_r, Nx * Ny * Nz);
  SetLength(mu_r,  Nx * Ny * Nz);
  SetLength(E_hist, steps);

  for var i := 0 to High(eps_r) do
  begin
    eps_r[i] := 1.0;
    mu_r[i]  := 1.0;
  end;

  W := Nx - 1;
  H := Ny - 1;
  var frameArr := new Real[W * H];

  // ================= Временной цикл =================
  for var n := 0 to steps - 1 do
  begin
    // ---------- Обновление H ----------
    for var i := 1 to Nx - 2 do
      for var j := 1 to Ny - 2 do
        for var k := 1 to Nz - 2 do
        begin
          var dEz_dy := (Ez[idxEz(i, j, k)] - Ez[idxEz(i, j - 1, k)]) / dy;
          var dEy_dz := (Ey[idxEy(i, j, k)] - Ey[idxEy(i, j, k - 1)]) / dz;
          var mu := mu0 * mu_r[idxEps(i, j, k)];
          Hx[idxHx(i - 1, j, k)] += -(dEz_dy - dEy_dz) / mu * dt;
        end;

    for var i := 1 to Nx - 2 do
      for var j := 1 to Ny - 2 do
        for var k := 1 to Nz - 2 do
        begin
          var dEx_dz := (Ex[idxEx(i, j, k)] - Ex[idxEx(i, j, k - 1)]) / dz;
          var dEz_dx := (Ez[idxEz(i, j, k)] - Ez[idxEz(i - 1, j, k)]) / dx;
          var mu := mu0 * mu_r[idxEps(i, j, k)];
          Hy[idxHy(i, j - 1, k)] += -(dEx_dz - dEz_dx) / mu * dt;
        end;

    for var i := 1 to Nx - 2 do
      for var j := 1 to Ny - 2 do
        for var k := 1 to Nz - 2 do
        begin
          var dEy_dx := (Ey[idxEy(i, j, k)] - Ey[idxEy(i - 1, j, k)]) / dx;
          var dEx_dy := (Ex[idxEx(i, j, k)] - Ex[idxEx(i, j - 1, k)]) / dy;
          var mu := mu0 * mu_r[idxEps(i, j, k)];
          Hz[idxHz(i, j, k - 1)] += -(dEy_dx - dEx_dy) / mu * dt;
        end;

    // ---------- Обновление E ----------
    for var i := 0 to Nx - 1 do
      for var j := 0 to Ny - 2 do
        for var k := 0 to Nz - 2 do
        begin
          var dHz_dy := (Hz[idxHz(i, j + 1, k)] - Hz[idxHz(i, j, k)]) / dy;
          var dHy_dz := (Hy[idxHy(i, j, k + 1)] - Hy[idxHy(i, j, k)]) / dz;
          var eps := eps0 * eps_r[idxEps(i, j, k)];
          Ex[idxEx(i, j, k)] += (dHz_dy - dHy_dz) / eps * dt;
          if (i = src_i) and (j = src_j - 1) and (k = src_k - 1) then
            Ex[idxEx(i, j, k)] += dt * J0 * source(n) / eps;
        end;

    for var i := 0 to Nx - 2 do
      for var j := 0 to Ny - 1 do
        for var k := 0 to Nz - 2 do
        begin
          var dHx_dz := (Hx[idxHx(i, j, k + 1)] - Hx[idxHx(i, j, k)]) / dz;
          var dHz_dx := (Hz[idxHz(i + 1, j, k)] - Hz[idxHz(i, j, k)]) / dx;
          var eps := eps0 * eps_r[idxEps(i, j, k)];
          Ey[idxEy(i, j, k)] += (dHx_dz - dHz_dx) / eps * dt;
        end;

    for var i := 0 to Nx - 2 do
      for var j := 0 to Ny - 2 do
        for var k := 0 to Nz - 1 do
        begin
          var dHy_dx := (Hy[idxHy(i + 1, j, k)] - Hy[idxHy(i, j, k)]) / dx;
          var dHx_dy := (Hx[idxHx(i, j + 1, k)] - Hx[idxHx(i, j, k)]) / dy;
          var eps := eps0 * eps_r[idxEps(i, j, k)];
          Ez[idxEz(i, j, k)] += (dHy_dx - dHx_dy) / eps * dt;
        end;

    ApplyPEC;

    // ---------- Энергия на каждом шаге ----------
    E_hist[n] := TotalEnergy;

    // ---------- Диагностика и отрисовка ----------
    if (n mod snap_every = 0) or (n = steps - 1) then
    begin
      WriteLn('n=', n:4,
              '  W=', E_hist[n]:0:6,
              '  max|divE|=", MaxDivergenceE:0:6);

      for var i := 0 to W - 1 do
        for var j := 0 to H - 1 do
          frameArr[i * H + j] := Ez[idxEz(i, j, slice_z)];

      LockDrawing;
      ClearWindow(clWhite);
      DrawSlice(frameArr, AnimX, AnimY, AnimSize);
      DrawEnergyPlot(E_hist, n, PlotX, PlotY, PlotW, PlotH);
      Redraw;
      Sleep(10);
    end;
  end;

  WriteLn("Расчёт завершён.');
end.

Расчёт выполняется на сетке Йи — специальной смещённой сетке, предназначенной для численного решения уравнений Максвелла методом FDTD. Её предложил Кейн Йи (Kane Yee) в 1966 году.

Идея в том, что в обычной сетке все компоненты поля (Ex, Ey, Ez, Hx, Hy, Hz) располагались бы в одних и тех же узлах. В сетке Йи компоненты разнесены в пространстве на половину шага и, кроме того, сдвинуты относительно друг друга по времени.

Электрическое и магнитное поля обновляются попеременно:

сначала вычисляется En на шаге n;
затем на промежуточном временном шаге находится Hn+1/2;
после этого вычисляется En+1, и цикл повторяется.

Такой алгоритм называется leapfrog — «прыжок через шаг».

Отметим, что схема Йи немонотонна: на границах раздела диэлектриков её фактический порядок точности снижается до первого. Это показано в статье Ж. О. Домбровской и А. Н. Боголюбова, опубликованной в 2017 году: https://elibrary.ru/item.asp?edn=vvfyhz.

В дискретном смысле схема Йи обладает свойством сохранения энергии. Это классический результат, известный с первых работ Йи и подробно рассматриваемый в книгах по FDTD (см., например, Taflove A., Hagness S. C. Computational Electrodynamics: the Finite–Difference Time–Domain Method. Artech House, 2005.). Косвенно проверить это свойство можно, отслеживая энергию на каждом шаге; именно такой контроль я и добавил в программу.

В результате программа строит анимацию распространения электромагнитной волны — срез компоненты Ez:

Срез Ez. Гауссов сигнал. Код до отладки. Не уверен в верности-картинка должна быть симметриной. — Срез Ez для гауссова сигнала. Это версия до отладки; в корректности результата есть сомнения — изображение должно быть симметричным.

Что происходит во времени? На изначально жёлтом фоне из центра постепенно проявляется верхняя треугольная область, заполненная мелкими сине-жёлтыми квадратами, образующими шахматный узор. Внутри неё заметна цельная красная область овальной формы. По мере приближения к левому верхнему углу овал размывается и вытягивается, затем исчезает; к концу анимации остаётся лишь несколько ячеек. Одновременно видна волна в виде извилистой линии, цвет которой чередуется между жёлтым и синим. При движении к левому верхнему углу линия растягивается и выходит за пределы кадра, оставляя участок, направленный к правому верхнему углу. Нижняя треугольная область на протяжении всей анимации остаётся жёлтой.

Сопоставим увиденное с описаниями в научной литературе и разберём возможные артефакты:

Наблюдение	Ожидаемая картина по литературным данным	Физическая интерпретация	Оценка и источник
Волновой фронт расходится от центра	Гауссов импульс, возбуждённый в центре, распространяется наружу. Это типичный результат (Салливан, главы 1 и 3).	Основной лепесток импульса расходится от точечного источника. Наблюдение соответствует ожиданиям.	✓ https://lshoshia.science.tsu.ge/fdtd/Dennissullivan-fixed.pdf
Чередование жёлтых и синих участков	Гауссов импульс включает полуволны противоположных знаков: переднюю и следующую за ней.	В публикациях по FDTD поле часто показывают двухцветной картой: красным обозначают положительные значения Ez, синим — отрицательные.	✓ https://www.researchgate.net/publication/5313575_Amplified_total_internal_reflection_theory_analysis_and_demonstration_of_existence_via_FDTD
Овальная область вытягивается и размывается у угла	Это может быть следствием численной дисперсии: высокочастотные компоненты распространяются с другой скоростью, а сферический фронт становится анизотропным.	Фазовая скорость на сетке Йи зависит от направления. При S ≈ 0,55 фронт может становиться ромбовидным; наибольшая скорость наблюдается вдоль пространственной диагонали.	~ https://www.academia.edu/112873112/Reduction_of_numerical_dispersion_in_FDTD_method_through_artificial_anisotropy
Волна вытягивается и уходит «за угол»	Вблизи частоты Найквиста передний край пакета может двигаться быстрее номинальной скорости; вдоль диагонали сигнал приходит раньше.	При S ≈ 0,55 скорость вдоль диагонали максимальна, поэтому пакет раньше достигает угла. Затем может проявляться мода, распространяющаяся вдоль стенки. Последнее — предположение автора, требующее отдельной проверки.	~ https://www.researchgate.net/publication/3307223_FDTD_dispersion_revisited_faster-than-light_propagation
Красный овал вблизи центра — ближнее поле	Дискретный диполь на сетке Йи может быть анизотропным; ближнее поле сосредоточено в нескольких ячейках.	Источник ориентирован по x, а срез Ez проходит на расстоянии 1,5 ячейки от него. Дальнее поле Ez на таком срезе ослаблено (∝ n·r̂), поэтому преобладает ближнее. Сходная анизотропия точечного диполя на сетке Йи описана в статье 2025 года о резонансных диполях и эффекте Пёрселла.	✓ https://arxiv.org/html/2504.07444
Нижняя половина изображения всё время остаётся жёлтой	До прихода отражённых волн поле должно быть зеркально симметрично относительно источника.	Результат противоречит физическим ожиданиям и данным литературы. Повторная реализация схемы подтвердила, что нулевая область ограничена антидиагональю i + j = 99, проходящей через ячейку источника. Это дефект индексации, а не свойство метода.	✗ —
К концу анимации остаются несколько красных ячеек	Остаточное поле не должно занимать заметную часть цветовой шкалы.	Расчёт показывает, что обрезанный гауссов импульс оставляет паразитный заряд: суммарное значение источника по времени составляет примерно 17,7, а не ноль. Из-за этого сохраняется статическое поле. Нормировка каждого кадра отдельно растягивает его на всю палитру. Возможны также артефакты, связанные с выбором среза и компоненты поля.	~ — оценка по коду

После исправления программы я добавил график изменения энергии во времени:

Отладка кода+ график энергии. Сомневаюсь- должно быть плато энергии, а не экспонента. — Исправленная версия с графиком энергии. Остаются сомнения: ожидалось постоянное значение энергии, а не экспоненциальный рост.

Выводы

Реализация модели FDTD на PascalABC.NET показывает, что даже задачи с достаточно сложной физической моделью можно решать и называть наглядно средствами языка, который многие считают устаревшим. Библиотека GraphABC избавила от необходимости подключать отдельные инструменты визуализации: весь процесс — от постановки задачи до отображения результата — можно выполнить в одной программе.

При разработке особое внимание уделялось правильной индексации компонент на сетке Йи и соблюдению условия устойчивости CFL.

Код получился компактным и понятным, поэтому его можно использовать как основу для дальнейших экспериментов: например, добавить неоднородности среды (для этого потребуется усреднять ε и μ на сетке Йи, иначе схема потеряет консервативность), PML-слои или источники другого типа.

Автор не имеет специального физико-математического образования и не является профессиональным физиком!