Рівняння мілкої води: створіть симуляцію цунамі за годину

1. Припущення мілкої води

Цунамі у відкритому Тихому океані має довжину хвилі 100-200 км, а океан глибиною лише ~4 км — відношення глибини до довжини хвилі приблизно 1:30. Коли довжина хвилі набагато більша за глибину рідини, вертикальна швидкість і прискорення всередині стовпа рідини нехтовно малі порівняно з горизонтальними. Це дозволяє звести повні 3D нестисливі рівняння Нав'є-Стокса (див. нашу статтю про Нав'є-Стокса в атмосфері для загального вигляду) до набагато простішої 2D системи, що відстежує лише висоту поверхні води h(x, y, t) та усереднене по глибині горизонтальне поле швидкості u(x, y, t) — без жодного вертикального виміру.

2. Рівняння мілкої води

Відкидання вертикальної структури дає компактну систему з трьох пов'язаних диференціальних рівнянь:

h — повна висота водного стовпа, u і v — усереднені по глибині горизонтальні складові швидкості, а g — прискорення вільного падіння (9.81 м/с²). Праві частини — це сила градієнта тиску: гідростатичний тиск у будь-якій точці пропорційний стовпу води над нею, тож похила поверхня штовхає рідину вниз по схилу точно так само, як сила градієнта тиску в повному рівнянні імпульсу Нав'є-Стокса.

Для хвиль малої амплітуди можна повністю відкинути нелінійні адвективні члени (u∂u/∂x тощо), отримавши лінеаризовані рівняння мілкої води — достатньо добре для поширення цунамі у відкритому океані, хоча повні нелінійні члени стають важливими, коли хвиля посилюється й стає крутішою біля берега.

3. Швидкість хвилі та shoaling

Лінеаризація й поєднання трьох рівнянь дає класичне хвильове рівняння зі швидкістю:

H — незбурена глибина води. В океані глибиною 4 км, c = √(9.81 × 4000) ≈ 198 м/с — понад 700 км/год, швидко як реактивний літак, тому цунамі, викликане землетрусом, може досягти узбережжя за тисячі кілометрів протягом кількох годин. Оскільки c залежить від глибини, коли хвиля наближається до берега й H зменшується, хвиля сповільнюється.

Збереження потоку енергії (приблизно, потік енергії ∝ амплітуда² × групова швидкість має залишатись сталим) означає, що коли швидкість падає, амплітуда має зростати для компенсації — явище, зване shoaling. Закон Гріна дає поширене наближення:

Цунамі, що ледь помітне як хвиля 0.5 м в океані глибиною 4 км, може вирости до кількох метрів, коли H падає до кількох десятків метрів біля берега — саме те, що передбачають рівняння мілкої води, і саме тому землетрус, відчутий у морі, породжує руйнівну хвилю на суходолі.

4. Скінченнорізницева дискретизація

Щоб змоделювати це на комп'ютері, дискретизуємо h, u і v на регулярній 2D сітці з кроком dx, і замінюємо просторові похідні центральними різницями:

та просуваємо сітку вперед у часі простим явним (Ейлер вперед) кроком: обчислюємо всі просторові похідні у поточний час, використовуємо їх для обчислення швидкості зміни h, u, v, а потім просуваємо все на dt.

5. Умова стійкості CFL

Явні схеми стійкі лише якщо інформація не може перетнути більше однієї комірки сітки за часовий крок. Оскільки найшвидший сигнал у цій системі поширюється зі швидкістю хвилі c = √(gH), умова Куранта-Фрідріхса-Леві (CFL) вимагає:

Виберіть dt більший — і симуляція розійдеться в NaN протягом кількох кадрів — це найпоширеніший баг при першій реалізації розв'язувача мілкої води.

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

class ShallowWaterSim {
  constructor(nx, ny, dx) {
    this.nx = nx; this.ny = ny; this.dx = dx;
    this.g = 9.81;
    this.H  = new Float32Array(nx * ny).fill(50);  // глибина у спокої (батиметрія)
    this.h  = new Float32Array(nx * ny);              // збурення висоти поверхні
    this.u  = new Float32Array(nx * ny);
    this.v  = new Float32Array(nx * ny);
  }

  idx(i, j) { return j * this.nx + i; }

  // Додає гаусів горб — симулює початкове зміщення, викликане землетрусом
  addDisturbance(cx, cy, amplitude, radius) {
    for (let j = 0; j < this.ny; j++)
      for (let i = 0; i < this.nx; i++) {
        const d = Math.hypot(i - cx, j - cy);
        this.h[this.idx(i, j)] += amplitude * Math.exp(-(d * d) / (2 * radius * radius));
      }
  }

  step(dt) {
    const { nx, ny, dx, g, H, h, u, v } = this;
    const hNew = new Float32Array(h.length);
    const uNew = new Float32Array(u.length);
    const vNew = new Float32Array(v.length);

    for (let j = 1; j < ny - 1; j++)
      for (let i = 1; i < nx - 1; i++) {
        const k = this.idx(i, j);
        const total = H[k] + h[k]; // повна висота стовпа

        // Центральні різниці для градієнтів
        const dhdx = (h[this.idx(i+1,j)] - h[this.idx(i-1,j)]) / (2 * dx);
        const dhdy = (h[this.idx(i,j+1)] - h[this.idx(i,j-1)]) / (2 * dx);

        // Імпульс: du/dt = -g * dh/dx  (лінеаризовано, без адвекції)
        uNew[k] = u[k] - g * dhdx * dt;
        vNew[k] = v[k] - g * dhdy * dt;

        // Нерозривність: dh/dt = -div(total * u)
        const fluxRight = total * u[this.idx(i+1,j)];
        const fluxLeft  = total * u[this.idx(i-1,j)];
        const fluxUp    = total * v[this.idx(i,j+1)];
        const fluxDown  = total * v[this.idx(i,j-1)];
        const div = (fluxRight - fluxLeft) / (2 * dx) + (fluxUp - fluxDown) / (2 * dx);
        hNew[k] = h[k] - div * dt;
      }

    this.h = hNew; this.u = uNew; this.v = vNew;
    this.applyBoundaries();
  }
}

7. Відбивні межі та узбережжя

Край області потребує граничної умови, інакше скінченнорізницевий шаблон читатиме за межами масиву. Найпростіший фізично обґрунтований вибір — відбивна межа (тверда стіна): обнулити нормальну складову швидкості й віддзеркалити градієнт висоти, щоб хвилі відбивались назад замість того, щоб зникати чи перегортатись на протилежний край.

applyBoundaries() {
  const { nx, ny, u, v, h } = this;
  for (let j = 0; j < ny; j++) {
    u[this.idx(0, j)]      = 0; u[this.idx(nx - 1, j)] = 0;
    h[this.idx(0, j)]      = h[this.idx(1, j)];
    h[this.idx(nx - 1, j)] = h[this.idx(nx - 2, j)];
  }
  // (повторити симетрично для верхнього/нижнього країв, використовуючи v)
}

Щоб симулювати пляж і побачити shoaling у дії, зробіть глибину води у спокої H мілкішою до одного краю сітки замість сталої — лінійна рампа від 50 м до 0.5 м на останніх 20% області достатня, щоб побачити, як хвиля сповільнюється й помітно зростає в амплітуді точно так, як передбачає закон Гріна, прямо перед тим, як вона «ламається», коли лінійне наближення руйнується на дуже малій глибині.