Теплообмін: числові методи скінченних різниць і скінченних елементів

1. Рівняння теплопровідності і навіщо його дискретизувати

Закон Фур'є стверджує, що тепловий потік пропорційний від'ємному градієнту температури, q = −k∇T. Разом із законом збереження енергії це дає рівняння дифузії тепла — параболічне диференціальне рівняння у частинних похідних:

Аналітичні розв'язки існують для нескінченного стрижня, напівнескінченної пластини або простих розкладів у ряди Фур'є на прямокутниках із тривіальними граничними умовами. Реальний об'єкт — радіатор, кристал процесора, шибка вікна — має неправильну геометрію і змішані граничні умови, тож неперервне поле T(x,y,t) замінюють значеннями на дискретній сітці, а похідні — скінченними різницями (або, для неправильних форм, базисними функціями методу скінченних елементів).

2. Метод скінченних різниць: явна схема

Дискретизуємо простір сіткою з кроком Δx, а час — кроками Δt. Друга похідна за простором стає тришаблонним трафаретом, а похідна за часом — прямою різницею: явна схема (FTCS):

Кожне нове значення залежить лише від старих — жодної системи рівнянь розв'язувати не потрібно, легко векторизувати, і саме так власна симуляція сайту «Рівняння теплопровідності 2D» рахує сітку 80×80 щокадру. Проблема — стійкість.

3. Стійкість: умова CFL / фон Неймана

Явна FTCS умовно стійка. Аналіз стійкості фон Неймана (розклад похибки у ряд Фур'є і перевірка, що коефіцієнт підсилення не перевищує 1) дає жорстке обмеження на крок часу:

4. Неявна схема та схема Кранка-Ніколсона

Неявна схема (BTCS) обчислює просторовий трафарет на новому шарі часу, що вимагає розв'язку системи лінійних рівнянь на кожному кроці, але повністю знімає обмеження стійкості:

BTCS безумовно стійка при будь-якому Δt, але лише першого порядку точності за часом; Кранк-Ніколсон — другого порядку точності і теж безумовно стійка ціною розв'язку тридіагональної (1D) або стрічкової (2D, через ADI — метод змінних напрямків) системи на кожному кроці. Для 2D-сіток ADI розбиває один повний крок на два піввкроки, кожен з яких розв'язує дешеву тридіагональну систему вздовж однієї осі, що тримає вартість близькою до явної схеми, знімаючи при цьому обмеження CFL.

5. Метод скінченних елементів: основна ідея

FDM потребує регулярної сітки; метод скінченних елементів працює з довільною геометрією — L-подібним кронштейном, зубцем шестерні, кузовом авто — покриваючи область трикутниками (2D) або тетраедрами (3D) і представляючи температуру як кусково-лінійну функцію на сітці, зважену вузловими значеннями.

Після складання M і K є розрідженими матрицями, і похідна за часом усе одно обробляється явною, неявною чи схемою Кранка-Ніколсона, як у §2–4 — FEM змінює спосіб дискретизації простору, а не кроку в часі. Виграш у тому, що густину сітки можна варіювати локально (дрібнішу біля гострих кутів, грубішу у відкритих ділянках), чого рівномірна сітка FDM зробити без значної переробки не може.

6. Граничні умови на практиці

7. Реалізація на JavaScript

// Явний розв'язувач FTCS для 2D рівняння теплопровідності з граничними умовами Діріхле
function stepHeat2D(T, next, nx, ny, alpha, dt, dx) {
  const r = alpha * dt / (dx * dx);
  if (r > 0.24) throw new Error(`Нестійко: r=${r.toFixed(3)} перевищує 0.25 (2D межа CFL)`);

  for (let j = 1; j < ny - 1; j++) {
    for (let i = 1; i < nx - 1; i++) {
      const idx = j * nx + i;
      const lap = T[idx + 1] + T[idx - 1] + T[idx + nx] + T[idx - nx] - 4 * T[idx];
      next[idx] = T[idx] + r * lap;
    }
  }
  // Граничні умови Діріхле: тримати краєві рядки/стовпці на попередньому значенні
  for (let i = 0; i < nx; i++) { next[i] = T[i]; next[(ny-1)*nx+i] = T[(ny-1)*nx+i]; }
  for (let j = 0; j < ny; j++) { next[j*nx] = T[j*nx]; next[j*nx+nx-1] = T[j*nx+nx-1]; }
  return next;
}

// Алгоритм Томаса — розв'язує тридіагональну систему за O(N), для 1D неявної/Кранка-Ніколсона
function thomasSolve(a, b, c, d) {
  const n = d.length;
  const cp = new Float64Array(n), dp = new Float64Array(n);
  cp[0] = c[0] / b[0]; dp[0] = d[0] / b[0];
  for (let i = 1; i < n; i++) {
    const m = b[i] - a[i] * cp[i-1];
    cp[i] = c[i] / m;
    dp[i] = (d[i] - a[i] * dp[i-1]) / m;
  }
  const x = new Float64Array(n);
  x[n-1] = dp[n-1];
  for (let i = n-2; i >= 0; i--) x[i] = dp[i] - cp[i] * x[i+1];
  return x;
}

8. Застосування та межі

• Охолодження електроніки: дизайн ребер радіатора зазвичай розв'язується 3D-FEM із граничними умовами Робіна, що моделюють конвекцію до рухомого повітря.
• Фізика будівель: теплові моделі цілих будівель використовують грубі неявні сітки FDM з кроком у годину — на цьому масштабі обмежуючим фактором зазвичай є стійкість, а не точність.
• Лиття та зварювання: фазовий перехід (прихована теплота) додає нелінійний член джерела, зазвичай оброблюваний «методом ентальпії», щоб не відстежувати фронт кристалізації явно.
• Прискорення на GPU: явна FDM ідеально паралелізується і напряму лягає на compute-шейдери — саме тому браузерна демонстрація вище працює інтерактивно на ноутбучній GPU.
• Коли неявна схема не варта того: для симуляцій реального часу, де малий крок часу вже потрібен для візуальної плавності, явна FTCS з невеликим запасом зазвичай простіша і достатньо швидка.