Доріжка Кармана: нестійкість потоку, число Струхаля та LBM-симуляція

1. Режими числа Рейнольдса

Число Рейнольдса Re = ρUD/μ = UD/ν (де U — швидкість набігаючого потоку, D — діаметр циліндра, ν = μ/ρ — кінематична в'язкість) керує топологією потоку:

Початок сходу вихорів при Re ≈ 47 — це надкритична біфуркація Гопфа: стаціонарний симетричний розв'язок для сліду втрачає стійкість і народжується граничний цикл (періодичний схід). Параметр біфуркації — Re; біфуркація негістерезисна (надкритична).

2. Число Струхаля та частота сходу вихорів

Безрозмірна частота сходу вихорів — це число Струхаля:

Для кругового циліндра в режимі ламінарного сходу (50 < Re < 1000) Williamson & Brown (1998) дають емпіричну апроксимацію:

За один цикл сходить два вихори протилежного знаку — по одному з кожного боку — тож візерунок сліду має просторовий період 2 × (U/f) × (поперечна відстань). Відстань між однознаковими вихорами вздовж потоку (довжина хвилі вихрової доріжки) дорівнює λ = U/f = D/St.

Щоб вихрова доріжка була стійкою, відношення поперечної відстані до поздовжньої має задовольняти критерій стійкості фон Кармана:

Будь-яке інше відношення спричиняє взаємодію двох рядів і швидко руйнує візерунок. Цю конкретну геометрію вивів фон Карман (1911) з аналізу стійкості подвійного ряду точкових вихорів — елегантне застосування теорії 2D-потенціального потоку.

3. Лінійна стійкість та нестійкість сліду

Перехід від стаціонарного до періодичного потоку при Re ≈ 47 пояснюється аналізом лінійної стійкості рівнянь Нав'є–Стокса, лінеаризованих навколо стаціонарного базового потоку U(x):

При Re < 47 усі власні значення L мають від'ємні дійсні частини — збурення згасають. При Re = 47 комплексно-спряжена пара λ = ±iω₀ перетинає уявну вісь (умова Гопфа), і стаціонарний слід замінюється наростаючим коливанням з частотою ω₀ = 2πSt·U/D. Понад Re = 47 амплітуда насичується через нелінійні ефекти; потік виходить на стійкий граничний цикл (періодичну вихрову доріжку).

Абсолютна проти конвективної нестійкості

Відмінність між абсолютною та конвективною нестійкістю критична для розуміння того, чому вихрова доріжка зберігається:

• Конвективна нестійкість: збурення наростають, але зносяться за потоком (аналогічно радіаційному загасанню). Слід діє як підсилювач, але не як самопідтримний осцилятор.
• Абсолютна нестійкість: збурення наростають у фіксованому місці в просторі. Потік діє як самозбудний осцилятор з чітко визначеною власною частотою.

Кох (1985) показав, що ближній слід циліндра є абсолютно нестійким при Re > 25–35. Ця абсолютна нестійкість «захоплює» частоту сходу, роблячи вихрову доріжку глобальною нестійкістю, а не конвективною хвилею, підсиленою шумом.

4. Опір, підйомна сила та тиск

Почерговий схід вихорів породжує періодично змінну підйомну силу (поперек потоку) з частотою f_s та слабко модульовану силу опору з частотою 2f_s:

Коливання опору з частотою 2f_s виникає через симетричне формування вихорів по черзі з обох боків — кожен новий вихор додає додатній приріст до опору незалежно від того, з якого боку він сходить. Коливання підйомної сили з частотою f_s відповідає почерговому знаку циркуляції в кожному вихорі, що сходить.

Поле тиску розкладається на: (1) середній дефіцит тиску в сліді (відповідальний за опір форми), (2) осцилююче прислідне тиск з періодом 1/f_s та (3) тиск приєднаного вихора (підсмоктування з того боку, де згортається зсувний шар).

5. Вихрозбудні коливання та захоплення частоти

Якщо циліндр закріплено пружно (маса m, пружина k, демпфування c), він може вібрувати в напрямку поперек потоку. Коли частота сходу f_s наближається до власної частоти конструкції f_n = √(k/m) / (2π), виникає явище, що зветься захопленням частоти (lock-in) (або синхронізацією):

• Частота сходу «захоплюється» до f_n у діапазоні зведених швидкостей U_r = U/(f_n D) ≈ 4–8
• У діапазоні захоплення циліндр коливається на f_n незалежно від U
• Амплітуда може сягати 1–2 діаметрів циліндра — набагато більше, ніж коливання ≈0.01D поза захопленням
• Енергія відбирається з потоку й передається в коливання конструкції

Захоплення частоти використовується конструктивно у збирачах енергії потоку (п'єзоелектричні накладки на осцилюючих конструкціях) та в русі риб: форель, що пливе за погано обтічним тілом D-подібної форми, «серфить» доріжку Кармана, відбираючи енергію зі сходячих вихорів через узгоджені в часі хвилеподібні рухи тіла — це експериментально продемонстрували Liao et al. (2003) у Science.

6. LBM D2Q9: теорія

Метод ґраткового Больцмана (LBM) розв'язує дискретне рівняння Больцмана на ґратці. У D2Q9 (2D, 9 швидкостей) кожна комірка зберігає 9 функцій розподілу fᵢ, що представляють популяції частинок, які рухаються в 9 дискретних напрямках:

Час релаксації τ керує в'язкістю: ν = cs²(τ − ½) = (1/3)(τ − ½). Для τ ∈ (0.5, 2) метод стійкий. Число Рейнольдса в одиницях LBM дорівнює Re = U_LBM · D_LBM / ν_LBM.

Граничні умови на поверхні циліндра використовують відбиття назад (bounce-back): будь-яка функція розподілу, що переноситься в тверде тіло, відбивається назад у протилежному напрямку з рівною величиною. Це реалізує умову прилипання з другим порядком точності для нерухомих меж.

7. LBM D2Q9: реалізація

// LBM D2Q9 — вихрова доріжка Кармана
// Сітка: ґратка NX × NY, циліндр у (cx, cy) радіусом R
// Відмальовується завихреність: ω = ∂v/∂x − ∂u/∂y (скінченна різниця)

const W = 320, H = 160;
const CX = 70, CY = H >> 1, R = 12;
const TAU = 0.56;          // час релаксації → ν = (0.56−0.5)/3 = 0.02
const U0  = 0.1;           // швидкість на вході (одиниці LBM, має бути ≪ 1)

// Ваги та вектори швидкостей D2Q9
const W9 = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
const EX = [0, 1,0,-1, 0, 1,-1,-1, 1];
const EY = [0, 0,1, 0,-1, 1, 1,-1,-1];
const OPP= [0, 3,4, 1, 2, 7,8, 5, 6];  // пари для відбиття назад

const N = W * H;
let f    = new Float32Array(9 * N);
let fTmp = new Float32Array(9 * N);
const solid = new Uint8Array(N);

// Позначаємо комірки циліндра
function initSolid() {
  for (let y = 0; y < H; y++) for (let x = 0; x < W; x++) {
    const dx = x - CX, dy = y - CY;
    if (dx*dx + dy*dy <= R*R) solid[y*W+x] = 1;
  }
}

// Рівноважний розподіл
function feq(i, rho, ux, uy) {
  const eu = EX[i]*ux + EY[i]*uy;
  const u2 = ux*ux + uy*uy;
  return W9[i] * rho * (1 + 3*eu + 4.5*eu*eu - 1.5*u2);
}

// Ініціалізуємо однорідний потік
function initFlow() {
  for (let n = 0; n < N; n++) {
    const rho0 = 1.0;
    for (let i = 0; i < 9; i++) f[9*n+i] = feq(i, rho0, U0, 0);
  }
}

function step() {
  const inv_tau = 1 / TAU;

  // Зіткнення
  for (let n = 0; n < N; n++) {
    if (solid[n]) continue;
    let rho = 0, ux = 0, uy = 0;
    for (let i = 0; i < 9; i++) {
      const fi = f[9*n+i];
      rho += fi; ux += EX[i]*fi; uy += EY[i]*fi;
    }
    ux /= rho; uy /= rho;
    for (let i = 0; i < 9; i++) {
      f[9*n+i] += inv_tau * (feq(i, rho, ux, uy) - f[9*n+i]);
    }
  }

  // Перенесення + граничні умови
  for (let y = 0; y < H; y++) for (let x = 0; x < W; x++) {
    const n0 = y*W + x;
    if (solid[n0]) {
      // Відбиття назад для твердих комірок
      for (let i = 0; i < 9; i++) fTmp[9*n0+OPP[i]] = f[9*n0+i];
      continue;
    }
    for (let i = 0; i < 9; i++) {
      let nx = x + EX[i], ny = y + EY[i];
      // Лівий вхід: вхід з постійною швидкістю (Zou-He)
      if (nx < 0) { fTmp[9*n0+i] = feq(i, 1.0, U0, 0); continue; }
      // Правий вихід: нульовий градієнт (копія від сусіда)
      if (nx >= W) { nx = W-1; }
      if (ny < 0)  ny = 0;
      if (ny >= H) ny = H-1;
      const nn = ny*W + nx;
      fTmp[9*nn+i] = f[9*n0+i];
    }
  }
  const tmp = f; f = fTmp; fTmp = tmp;
}

// Зчитуємо завихреність ω = ∂v/∂x − ∂u/∂y (центральні різниці)
function getVorticity(x, y) {
  const getU = (xx, yy) => {
    const n = Math.min(Math.max(yy,0),H-1)*W + Math.min(Math.max(xx,0),W-1);
    let rho=0, ux=0, uy=0;
    for (let i=0;i<9;i++){rho+=f[9*n+i];ux+=EX[i]*f[9*n+i];uy+=EY[i]*f[9*n+i];}
    return rho > 0 ? [ux/rho, uy/rho] : [0,0];
  };
  const [,vR] = getU(x+1, y); const [,vL] = getU(x-1, y);
  const [uT]  = getU(x, y+1); const [uB]  = getU(x, y-1);
  return 0.5 * ((vR - vL) - (uT - uB));
}
    

8. Інтерактивна симуляція вихрової доріжки

Полотно нижче запускає LBM D2Q9-симуляцію на ґратці 320×160 комірок. Завихреність ωz = ∂v/∂x − ∂u/∂y відмальовується за допомогою розбіжної синьо-біло-червоної кольорової мапи. Збільшіть Re, щоб побачити перехід від стаціонарного сліду до повного періодичного сходу вихорів.