Аэродинамическая труба

Как это устроено

Струи дыма наверху никто не рисовал. Пространство вокруг болида разбито на 2,4 миллиона ячеек, и видеокарта несколько сотен раз в секунду пересчитывает, как в каждой из них движется воздух. Дым просто летит по этому посчитанному ветру.

Зачем Формуле-1 аэродинамическая труба

Современный болид едет быстро не столько из-за мотора, сколько из-за воздуха. Крылья и днище прижимают машину к трассе силой, которая на высокой скорости превышает её собственный вес. Чем сильнее прижим, тем больше сцепление шин и тем быстрее можно пройти поворот. Расплата — сопротивление: воздух тормозит машину, и часть мощности мотора уходит только на то, чтобы его растолкать.

Команды проверяют свои идеи двумя способами: расчётом на суперкомпьютерах и продувкой в аэродинамической трубе. По правилам FIA модель в трубе — не больше 60 % настоящего размера, а скорость потока — не больше 50 м/с. Пол под моделью — бегущая лента, колёса крутятся. Число часов в трубе и объём расчётов ограничены, и отстающим командам разрешено больше. Наша труба — та же идея, только машина полноразмерная, а вместо вентилятора уравнения.

Уравнения Навье–Стокса

Движение воздуха при скоростях заметно ниже звуковой описывают уравнения Навье–Стокса. Это второй закон Ньютона для каждой капли воздуха плюс условие, что воздух не сжимается:

$$\rho\Big(\frac{\partial \mathbf u}{\partial t} + (\mathbf u\cdot\nabla)\,\mathbf u\Big) = -\nabla p + \mu\,\nabla^2 \mathbf u,\qquad \nabla\cdot\mathbf u = 0.$$

Слева — масса, умноженная на ускорение: капля ускоряется, если меняется поток в данном месте, и если её саму переносит в место, где поток другой. Справа — силы: перепад давления $\nabla p$ толкает воздух туда, где давление ниже, а вязкость $\mu$ сглаживает разницу скоростей между соседними слоями.

Насколько важна вязкость, говорит число Рейнольдса — отношение инерции к трению:

$$\mathrm{Re} = \frac{U L}{\nu} = \frac{69{,}4\ \text{м/с}\times 5{,}3\ \text{м}}{1{,}5\cdot 10^{-5}\ \text{м}^2/\text{с}} \approx 2{,}5\cdot 10^{7}.$$

Двадцать пять миллионов: вязкость важна только в тонком пограничном слое у кузова толщиной в миллиметры, а за машиной поток рвётся на вихри всех размеров — это турбулентность. Готовой формулы для такого течения нет, его можно только считать.

Решётка Больцмана

Решать эти уравнения можно по-разному. Мы взяли решёточный метод Больцмана — на нём работают и учёные (FluidX3D), и промышленные программы для автомобилей (PowerFLOW). Вместо скорости и давления в каждой ячейке хранятся 19 чисел $f_i$. Каждое число — сколько воздуха летит в одном из 19 направлений $\mathbf c_i$: стоит на месте, к одной из 6 соседних граней или к одному из 12 рёбер (схема D3Q19). Плотность и скорость — просто суммы:

$$\rho = \sum_i f_i,\qquad \rho\,\mathbf u = \sum_i f_i\,\mathbf c_i.$$

Шаг состоит из двух действий. Перенос: каждое $f_i$ перелетает в соседнюю ячейку по своему направлению. Столкновение: в каждой ячейке числа сдвигаются к равновесному распределению — так молекулы, сталкиваясь, приходят к равновесию:

$$f_i(\mathbf x + \mathbf c_i,\ t + 1) = f_i^{\mathrm{eq}} + \Big(1 - \frac{1}{\tau}\Big)\big(f_i - f_i^{\mathrm{eq}}\big),\qquad f_i^{\mathrm{eq}} = w_i\,\rho\Big(1 + 3\,\mathbf c_i\!\cdot\mathbf u + \tfrac92 (\mathbf c_i\!\cdot\mathbf u)^2 - \tfrac32\,\mathbf u^2\Big) + \dots$$

Из этих простых правил получаются в точности уравнения Навье–Стокса — это показывает разложение Чепмена–Энскога. Вязкость задаётся временем релаксации: $\nu = c_s^2\,(\tau - \tfrac12)$, где $c_s^2 = 1/3$ — квадрат скорости звука решётки. Главное удобство метода — всё локально: каждая ячейка смотрит только на соседей. Миллионы ячеек считаются одновременно, видеокарта делает два миллиарда таких обновлений в секунду.

Масштабы. Ячейка — около 3,8 см, труба — 14 × 4,4 × 2,2 м. Поток в решётке идёт со скоростью 0,1 ячейки за шаг, то есть число Маха около 0,17, как у настоящей машины на 250 км/ч. Отсюда шаг по времени: $\Delta t = h\,u_\text{реш}/U \approx 5{,}5\cdot 10^{-5}$ с. За кадр видеокарта делает около десятка шагов, поэтому показ замедлен примерно в 40 раз; точное число видно в табличке справа.

Турбулентность, которую не разрешить

Чтобы честно разрешить все вихри при $\mathrm{Re} = 2{,}5\cdot 10^7$, нужна сетка мельче десятой доли миллиметра — порядка $10^{15}$ ячеек. Поэтому вязкость в расчёте настоящая, воздушная, а влияние вихрей мельче ячейки добавляет подсеточная модель Смагоринского (метод крупных вихрей, LES). Там, где поток сильно сдвигается, вязкость растёт:

$$\nu_t = (C_s\,\Delta)^2\,|S|,\qquad C_s = 0{,}12.$$

В решётке это сводится к своему $\tau$ в каждой ячейке. Его находят по неравновесной части потока импульса $\Pi^{\mathrm{neq}}$, без производных: $\tau = \tfrac12\big(\tau_0 + \sqrt{\tau_0^2 + 18\sqrt2\,C_s^2\,|\Pi^{\mathrm{neq}}|/\rho}\big)$.

Как заставить решётку не взорваться

При настоящей вязкости воздуха $\tau - \tfrac12 \approx 10^{-6}$: схема почти не гасит ничего. Первая версия нашей трубы честно это показала. За машиной рождалась рябь с периодом в две ячейки, и через пару тысяч шагов она заливала всю трубу — даже воздух перед машиной дрожал на 30 % скорости. В пустом ящике с потоком 0,1 и крошечным шумом то же самое: ошибка росла экспоненциально. Так проявляется известная беда схемы: при движущемся потоке у её второго порядка ошибки растут как $u^3$.

Помогли три вещи:

• Рекурсивная регуляризация третьего порядка (O. Malaspinas, 2015; C. Coreixas и др., 2017). Неравновесную часть $f_i$ перестраиваем из $\Pi^{\mathrm{neq}}$, а её третьи моменты вычисляем рекурсивно: $a^{(3)}_{\alpha\beta\gamma} = u_\alpha a^{(2)}_{\beta\gamma} + u_\beta a^{(2)}_{\alpha\gamma} + u_\gamma a^{(2)}_{\alpha\beta}$. Равновесие тоже берём с членами третьего порядка — теми шестью, что решётка D3Q19 умеет хранить. После этого посеянный шум гаснет даже при $\tau - \tfrac12 = 10^{-6}$.
• Объёмная вязкость (Деллар, 2001). След тензора $\Pi^{\mathrm{neq}}$ отвечает за сжатие, то есть за звук; его мы релаксируем отдельно и сильно. Звуковые волны, которыми турбулентный след «звенит» на всю трубу, гаснут, а сдвиговое течение не трогается — для несжимаемого воздуха это законно.
• Поглощающие слои на входе и выходе трубы: волны и след уходят, не отражаясь обратно.

С этим шум в набегающем потоке — 0,03 % скорости, и это при настоящей вязкости воздуха.

Болид в решётке

Модель болида сначала разбивается на кубики по 1 см: поверхность и всё, что внутри неё. Кокпит, воздухозаборник и входы в понтоны закрыты: у настоящей машины там пилот, радиаторы и мотор, и воздух сквозь кузов не проходит. Потом каждая ячейка решётки получает долю твёрдого $\varepsilon$ — среднее по 64 точкам внутри неё. С учётом угла к потоку, высоты кузова и угла закрылка DRS это пересчитывается за миллисекунду, поэтому ручки работают плавно.

В ячейках с твёрдым воздух подтягивается к скорости поверхности. Так устроен метод частично заполненных ячеек (Ноубл и Торчинский, 1998):

$$f_i \leftarrow f_i^{\mathrm{eq}} + \Big(1 - \frac{1 - B}{\tau}\Big) f_i^{\mathrm{neq}} + B\,\Omega^s_i,\qquad \Omega^s_i = f_{\bar\imath} - f_i + f_i^{\mathrm{eq}}(\rho, \mathbf u_s) - f_{\bar\imath}^{\mathrm{eq}}(\rho, \mathbf u).$$

Вес $B$ растёт от 0 до 1 вместе с $\varepsilon$. Элементы крыльев тоньше ячейки, поэтому они считаются сплошными: иначе воздух проходил бы сквозь них. У колёс скорость поверхности — вращение $\mathbf u_s = \boldsymbol\omega \times \mathbf r$, причём низ колеса едет вместе с дорогой. Сама дорога — движущаяся стенка: отражённое от неё $f$ получает добавку $6\,w_i\,\rho\,\mathbf c_i\cdot\mathbf u_w$.

Как считается сила

Всё, что воздух потерял на кузове, кузов получил. Импульс, который столкновение отбирает у воздуха в каждой твёрдой ячейке, суммируется по частям машины: переднее крыло, днище, заднее крыло, колёса, остальной кузов.

$$\mathbf F = -\sum_{\text{ячейки}} B \sum_i \Omega^s_i\,\mathbf c_i,\qquad F_\text{Н} = F_\text{реш}\cdot\frac{\rho_\text{возд}\,h^4}{\Delta t^2}.$$

Силы скачут вместе с турбулентностью, поэтому показываем среднее за последние ~0,08 секунды течения. Килограмм-сила — это вес одного килограмма, $9{,}81$ Н.

Что упрощено и где числа расходятся с правдой

Течение здесь честное, но сетка грубая. Поэтому абсолютные силы отличаются от настоящих, а вот изменения от DRS, клиренса и грязного воздуха показывают правильное направление. Вот что даёт расчёт на 250 км/ч при ячейке 3,8 см, а рядом — оценки для настоящего болида по открытым данным:

	Наш расчёт	Настоящий болид (оценка)
Сопротивление	≈ 800–1000 кгс	≈ 350 кгс ($C_D A \approx 1{,}2\ \text{м}^2$)
Вертикальная сила	≈ 500 кгс вверх	≈ 1200 кгс вниз ($C_L A \approx 4\ \text{м}^2$)
Переднее + заднее крыло	≈ 150 кгс вниз	≈ 600 кгс вниз

Откуда расхождение:

• Пограничный слой не разрешён. У настоящего кузова он тоньше миллиметра, у нас ячейка 3,8 см. Без пристеночной модели скруглённые тела — шины, понтоны — ведут себя как при маленьком числе Рейнольдса: поток отрывается раньше, след шире, сопротивление выше. Проверка на шаре это подтверждает: $C_D \approx 0{,}33$ вместо 0,1–0,2 у настоящего шара на таком Re.
• Щель под днищем — две ячейки. Граунд-эффект, главный источник прижима, держится на течении в зазоре в несколько сантиметров. Его нужно разрешать в десятки ячеек, иначе днище недорабатывает. Поэтому по умолчанию машина приподнята на 3,5 см.
• Упрощённое днище. У бесплатной модели оно плоское, без тоннелей Вентури, как у машин 2009–2021 годов.
• Тесная труба. Машина закрывает около 15 % сечения. Это обычно и для настоящих труб, но там вносят поправку, а мы показываем «сырые» цифры.
• Нет радиаторов, тормозов и выхлопа; шины жёсткие, подвеска не работает, поэтому нет и «дельфинирования».
• Модель — приблизительная. Это любительская 3D-модель, похожая на W14, а не настоящая геометрия машины.

С измельчением сетки цифры движутся в нужную сторону. На ячейке 3,0 см (качество «Максимум») крылья прижимают примерно в полтора раза сильнее, а суммарная подъёмная сила падает впятеро. Командам для правдивых чисел нужны сотни миллионов ячеек, пристеночные модели и проверка в трубе.

Как проверено

• Тот же шаг расчёта написан второй раз — на JavaScript, в двойной точности. На маленькой трубе с вращающимся шаром видеокарта совпадает с ним до $7\cdot10^{-7}$ в f32 и до $5\cdot10^{-5}$ в f16; силы совпадают до $10^{-4}$.
• Течение Куэтта (движущийся пол под неподвижной стенкой) даёт точный линейный профиль; это проверка движущейся стенки.
• Затухание вихря Тейлора–Грина совпадает с $e^{-2\nu k^2 t}$ с точностью 1 % — значит, вязкость задана верно.
• Обтекание цилиндра при Re = 100 даёт число Струхаля 0,16–0,17 (справочное 0,165).
• Пустая труба остаётся идеально однородной, масса сохраняется до $10^{-13}$.

Источники

Krüger T. et al. The Lattice Boltzmann Method: Principles and Practice. Springer, 2017. — Latt J., Chopard B. Lattice Boltzmann method with regularized pre-collision distribution functions // Math. Comput. Simul., 2006. — Malaspinas O. Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization // arXiv:1505.06900, 2015. — Coreixas C. et al. Recursive regularization step for high-order lattice Boltzmann methods // Phys. Rev. E 96, 033306, 2017. — Dellar P. Bulk and shear viscosities in lattice Boltzmann equations // Phys. Rev. E 64, 031203, 2001. — Noble D., Torczynski J. A lattice-Boltzmann method for partially saturated computational cells // Int. J. Mod. Phys. C 9, 1998. — Hou S. et al. A lattice Boltzmann subgrid model for high Reynolds number flows // Fields Inst. Commun. 6, 1996. — Lehmann M. FluidX3D — решётка Больцмана на GPU, приём FP16S. — Технический регламент FIA Formula 1 и Aerodynamic Testing Restrictions.