Рівняння Нав'є-Стокса у WebGL — симуляція рідини в реальному часі

1. Рівняння нестисливого Нав'є-Стокса

Рівняння нестисливого Нав'є-Стокса описують рух ньютонівської рідини зі сталою густиною. У векторній формі це два зв'язані диференціальні рівняння з частинними похідними:

Рівняння нерозривності ∇ · u = 0 — це обмеження нестисливості: рідина не може стискатися, тож дивергенція поля швидкостей має скрізь дорівнювати нулю. Це обмеження керує розв'язувачем тиску.

Рівняння імпульсу має чотири члени праворуч:

• Адвекція — −(u · ∇)u: рідина переносить власну швидкість (нелінійна, найскладніша частина)
• Дифузія — ν ∇²u: в'язке згладжування різниць швидкостей
• Градієнт тиску — −(1/ρ)∇p: різниці тиску спричиняють течію
• Об'ємні сили — f: зовнішні впливи (сплеск мишею, гравітація)

2. Розщеплення за операторами — підхід Stable Fluids

Метод Стама просуває розв'язок на один крок за часом Δt за допомогою чотирьох послідовних підкроків, кожен з яких застосовується до сітки поля, збереженої як текстура з плаваючою комою:

3. Напівлагранжева адвекція

Ідея: замість запитувати «куди прямує рідина в комірці сітки x?», запитуй «звідки рідина в x прийшла?» Простеж поле швидкостей назад на Δt і візьми там значення попереднього поля.

У GLSL (один канал текстури з плаваючою комою = компонента швидкості):

// Шейдер адвекції — gl_FragCoord.xy відображається в UV у [0,1]²
uniform sampler2D uVelocity;   // поточне поле швидкостей (vec2, закодований як RG)
uniform sampler2D uSource;     // поле для адвекції (швидкість або барвник)
uniform vec2      uTexelSize;  // 1/роздільна здатність
uniform float      uDt;

void main() {
  vec2 uv   = gl_FragCoord.xy * uTexelSize;
  vec2 vel  = texture(uVelocity, uv).rg;
  vec2 prev = uv - uDt * vel * uTexelSize; // зворотне трасування в просторі текстури
  gl_FragColor = texture(uSource, prev);  // білінійна вибірка (БЕЗКОШТОВНА на GPU)
}

4. Дифузія — ітерація Якобі

Неявне розв'язання (I − ν Δt ∇²) u_new = u_old відносно u_new дає розріджену лінійну систему. На регулярній сітці лапласіан дискретизується так:

// Шейдер Якобі (одна ітерація; ping-pong між двома FBO)
uniform sampler2D uX;        // поточне наближення x_k
uniform sampler2D uB;        // права частина b  (u_old)
uniform float      uAlpha;    // Δx² / (ν·Δt)
uniform float      uBeta;     // 1 / (4 + α)
uniform vec2      uTexelSize;

void main() {
  vec2 uv = gl_FragCoord.xy * uTexelSize;
  vec4 xL = texture(uX, uv - vec2(uTexelSize.x, 0));
  vec4 xR = texture(uX, uv + vec2(uTexelSize.x, 0));
  vec4 xB = texture(uX, uv - vec2(0, uTexelSize.y));
  vec4 xT = texture(uX, uv + vec2(0, uTexelSize.y));
  vec4 b  = texture(uB, uv);
  gl_FragColor = (xL + xR + xB + xT + uAlpha * b) * uBeta;
}

Для рідин із низькою в'язкістю можна взагалі пропустити явну дифузію (ν = 0) — самої числової дифузії від адвекції достатньо для згладжування, і це заощаджує 20–40 проходів Якобі на кадр.

5. Проєкція тиску — розклад Гельмгольца

Після адвекції поле швидкостей u̅ зазвичай дивергентне (∇ · u̅ ≠ 0) — несумісне з реальною поведінкою рідини. Розклад Гельмгольца гарантує, що будь-яке векторне поле можна розкласти на бездивергентну частину та градієнт:

// Шейдер дивергенції
void main() {
  vec2 uv = gl_FragCoord.xy * uTexelSize;
  float vL = texture(uVelocity, uv - vec2(uTexelSize.x,0)).x;
  float vR = texture(uVelocity, uv + vec2(uTexelSize.x,0)).x;
  float vB = texture(uVelocity, uv - vec2(0,uTexelSize.y)).y;
  float vT = texture(uVelocity, uv + vec2(0,uTexelSize.y)).y;
  gl_FragColor.r = 0.5 * (vR - vL + vT - vB);  // центральні різниці
}

// Шейдер віднімання градієнта
void main() {
  vec2 uv = gl_FragCoord.xy * uTexelSize;
  float pL = texture(uPressure, uv - vec2(uTexelSize.x,0)).r;
  float pR = texture(uPressure, uv + vec2(uTexelSize.x,0)).r;
  float pB = texture(uPressure, uv - vec2(0,uTexelSize.y)).r;
  float pT = texture(uPressure, uv + vec2(0,uTexelSize.y)).r;
  vec2  vel = texture(uVelocity, uv).xy;
  vel -= 0.5 * vec2(pR - pL, pT - pB);
  gl_FragColor = vec4(vel, 0, 1);
}

Розв'язання Пуассона для тиску використовує той самий шейдер Якобі, що й дифузія — параметри α = −Δx², β = 1/4. 20–50 ітерацій зазвичай достатньо для плавної течії; розв'язувач Гаусса–Зейделя чи багатосітковий збігався б швидше.

6. Реалізація у WebGL — текстури ping-pong

Розв'язувач рідини у WebGL підтримує кілька цілей рендерингу з плаваючою комою (пари FBO для ping-pong). Мінімальне налаштування потребує:

• швидкість — 2-компонентна (RG) float-текстура
• тиск — 1-компонентна (R) float-текстура
• дивергенція — 1-компонентна (R) float-текстура
• барвник — 4-компонентна (RGBA) float-текстура

// Мінімальний помічник FBO для WebGL2
function createDoubleFBO(gl, w, h, format) {
  const make = () => {
    const tex = gl.createTexture();
    gl.bindTexture(gl.TEXTURE_2D, tex);
    gl.texImage2D(gl.TEXTURE_2D, 0, format.internal,
                   w, h, 0, format.format, gl.FLOAT, null);
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.LINEAR);
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE);
    gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE);
    const fbo = gl.createFramebuffer();
    gl.bindFramebuffer(gl.FRAMEBUFFER, fbo);
    gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, tex, 0);
    return { tex, fbo };
  };
  let read = make(), write = make();
  return {
    get read()  { return read;  },
    get write() { return write; },
    swap() { [read, write] = [write, read]; }
  };
}

// Покадровий цикл оновлення
function update(dt) {
  addForces(velocity, mouseForce);          // крок 1
  for (let i = 0; i < 20; i++)
    jacobi(velocity, velocity, alpha, beta); // крок 2: дифузія
  advect(velocity, velocity, dt);           // крок 3: адвекція швидкості
  computeDivergence(velocity, divergence);  // крок 4a
  pressure.read.clearToZero();               // крок 4b: скинути p
  for (let i = 0; i < 40; i++)
    jacobi(pressure, divergence, -dx*dx, 0.25); // 4c: розв'язання тиску
  subtractGradient(velocity, pressure);     // крок 4d: проєкція
  advect(dye, velocity, dt);                // адвекція барвника
  renderDye(dye);
}

7. Адвекція барвника та візуалізація

Саме поле швидкостей невидиме. Дві популярні техніки візуалізації:

• Адвекція барвника: вприскуй колір у позицію миші й адвектуй його щокадру за допомогою поля швидкостей. Барвник накопичується й виявляє вихори, турбулентність і візерунки перемішування.
• Розфарбовування за ротором/завихреністю: обчисли ω = ∂v/∂x − ∂u/∂y (2D скаляр завихреності) та зістав додатні/від'ємні значення з двома кольорами. Ядра вихорів проявляються як яскраві плями.

// Шейдер ротора — кольоровий прохід утримання завихреності
void main() {
  vec2 uv = gl_FragCoord.xy * uTexelSize;
  float vL = texture(uVelocity, uv - vec2(uTexelSize.x, 0)).y;
  float vR = texture(uVelocity, uv + vec2(uTexelSize.x, 0)).y;
  float uB = texture(uVelocity, uv - vec2(0, uTexelSize.y)).x;
  float uT = texture(uVelocity, uv + vec2(0, uTexelSize.y)).x;
  float curl = 0.5 * ((vR - vL) - (uT - uB));  // ∂v/∂x − ∂u/∂y
  gl_FragColor = vec4(curl, 0.0, 0.0, 1.0);
}

// Сила утримання завихреності — повертає вихорам енергію, протидіючи дисипації
void main() {
  vec2  uv  = gl_FragCoord.xy * uTexelSize;
  float wL  = abs(texture(uCurl, uv - vec2(uTexelSize.x, 0)).r);
  float wR  = abs(texture(uCurl, uv + vec2(uTexelSize.x, 0)).r);
  float wB  = abs(texture(uCurl, uv - vec2(0, uTexelSize.y)).r);
  float wT  = abs(texture(uCurl, uv + vec2(0, uTexelSize.y)).r);
  float w   = texture(uCurl, uv).r;
  vec2  eta = normalize(vec2(wT - wB, wR - wL) + 1e-5);
  vec2  force = eta * w * uConfinement * uTexelSize;
  vec2  vel = texture(uVelocity, uv).xy + force * uDt;
  gl_FragColor = vec4(vel, 0, 1);
}

8. Розширення та подальше читання

• 3D Нав'є-Стокс: розширте кожну текстуру до 3D текстури (TEXTURE_3D у WebGL 2). Споживання пам'яті зростає як N³; сітка 64³ практична в реальному часі.
• Кращі розв'язувачі Пуассона: Якобі потребує багато ітерацій. Гаусс–Зейдель (червоно-чорне впорядкування) збігається у 2 рази швидше; багатосіткові методи збігаються за O(N) роботи незалежно від розміру сітки.
• Методи вихрових частинок: подають завихреність як дискретні точкові вихори, а не сітку — природно адаптивні й вільні від числової дифузії.
• SPH-рідина (на основі частинок): замість ейлерової сітки згладжена гідродинаміка частинок використовує лагранжеві частинки — подробиці див. у статті про SPH.
• Порівняння з LBM: метод ґраткового Больцмана досягає схожих результатів на цілком іншій математичній основі — мезоскопічній кінетичній теорії, а не макроскопічній механіці суцільного середовища.