Туториал: пишем решатель уравнений Навье-Стокса на Python с нуля для симуляции обтекания крыла

Зачем это нужно? (И почему это не безумие)

Когда слышишь "уравнения Навье-Стокса", первая мысль - бежать. Непрерывность, вязкость, давление, турбулентность. В теории - кошмар. На практике - основа всего, что летает, плавает и едет. От крыла самолета до потока воздуха в серверной.

Готовые CFD-пакеты вроде OpenFOAM или ANSYS Fluent - это здорово. До тех пор, пока не нужно понять, что именно они считают. Черный ящик, который иногда выдает красивую картинку, а иногда - физически невозможную ерунду. Писать свой решатель с нуля - это не академическое упражнение. Это способ заглянуть под капот самой физики течений. Узнать, где и почему все ломается. После этого даже готовые инструменты перестают быть магией.

Что нам понадобится (актуально на март 2026)

Библиотеки. Без них никуда, но и без фанатизма. Мы используем только необходимое, последние стабильные версии.

• Python 3.12+ - на момент 2026 года это стандарт. Особенно важны оптимизации работы с памятью.
• NumPy 2.0+ - основа всех вычислений. В версии 2.0 полностью переработана работа с типами данных и ускорены многие операции.
• SciPy 1.14+ - понадобится для решения систем линейных уравнений. Используем самые новые sparse solvers.
• Matplotlib 3.10+ - визуализация. Последние версии поддерживают аппаратное ускорение рендеринга.

# Установка всего необходимого (для Python 3.12)
pip install numpy==2.0.1 scipy==1.14.0 matplotlib==3.10.0

Уравнения, которые всем управляют

Навье-Стокса в полной красе выглядят пугающе. Но мы возьмем упрощенный случай: несжимаемое течение (плотность постоянна), двумерное. Так легче начать, но физика сохраняется.

В векторном виде уравнение выглядит так:

∂u/∂t + (u·∇)u = -∇p/ρ + ν∇²u

Плюс уравнение непрерывности для несжимаемой жидкости:

∇·u = 0

Где u - вектор скорости, p - давление, ρ - плотность (константа), ν - кинематическая вязкость. Физически: левая часть - это ускорение частицы жидкости, правая - силы давления и вязкости.

Как НЕ надо дискретизировать (чтобы не взорвать все)

Самая частая ошибка новичков - попытка записать разностную схему "в лоб". Берем производные, заменяем разностями, запускаем. Через три итерации значения скорости улетают в бесконечность. Почему? Уравнение Навье-Стокса нелинейное. Термин (u·∇)u - это адвекция, перенос скорости самим же потоком. Неустойчивость растет как снежный ком.

Второй подводный камень - связь скорости и давления. В каждый момент времени поле скорости должно быть бездивергентным (∇·u = 0). Но давление явно не входит в уравнение непрерывности. Получается скрытое ограничение. Решать напрямую - все равно что собрать пазл вслепую.

1 Выбираем метод: проекция для чайников

Мы используем метод проекции (fractional step method). Идея гениальна в простоте: разделить решение на два этапа. Сначала находим промежуточную скорость, игнорируя давление. Затем "проецируем" эту скорость на пространство бездивергентных полей, попутно находя давление.

Алгоритм на одном временном шаге:

1. Решаем уравнение для промежуточной скорости u* (без градиента давления)
2. Решаем уравнение Пуассона для давления p
3. Корректируем скорость: u^{n+1} = u* - Δt/ρ * ∇p

2 Строим сетку: где живут наши переменные

Используем размеченную сетку (staggered grid), еще ее называют сеткой Маркуса-Вольфа. Почему не обычная сетка? Потому что на обычной сетке возникает чётно-нечётная дискретизация, которая не чувствует колебаний давления. Это известная проблема, которая ломала многие первые решатели.

На размеченной сетке:

• Давление p хранится в центрах ячеек
• Скорость u (x-компонента) хранится на серединах вертикальных граней
• Скорость v (y-компонента) хранится на серединах горизонтальных граней

Такой хитрый трюк автоматически обеспечивает выполнение дискретного аналога уравнения непрерывности. Физически это значит, что поток через ячейку считается корректно.

import numpy as np

# Параметры задачи
nx, ny = 101, 101  # Количество узлов по x и y
lx, ly = 2.0, 2.0  # Размер области
dx = lx / (nx - 1)
dy = ly / (ny - 1)

# Инициализация полей на размеченной сетке
p = np.zeros((ny, nx))          # Давление в центрах ячеек
u = np.zeros((ny, nx + 1))      # x-скорость на вертикальных гранях
v = np.zeros((ny + 1, nx))      # y-скорость на горизонтальных гранях

# Временные параметры
dt = 0.001        # Шаг по времени
nt = 500          # Количество шагов
rho = 1.0         # Плотность
nu = 0.01         # Кинематическая вязкость

3 Пишем ядро решателя: от уравнений к коду

Теперь самая мякотка - реализация шага по времени. Распишем по формулам, что именно мы делаем на каждой итерации.

def solve_navier_stokes(u, v, p, dt, dx, dy, rho, nu, nt):
    """Основной цикл решения уравнений Навье-Стокса"""
    
    # Предварительный расчет констант
    dtdx = dt / dx
    dtdy = dt / dy
    dx2 = dx * dx
    dy2 = dy * dy
    
    for n in range(nt):
        # Сохраняем предыдущие значения скорости
        un = u.copy()
        vn = v.copy()
        
        # 1. Промежуточная скорость (без учета давления)
        # x-компонента
        u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
                         un[1:-1, 1:-1] * dtdx * (un[1:-1, 1:-1] - un[1:-1, :-2]) -
                         vn[1:-1, 1:-1] * dtdy * (un[1:-1, 1:-1] - un[:-2, 1:-1]) +
                         nu * dt * (
                             (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, :-2]) / dx2 +
                             (un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[:-2, 1:-1]) / dy2
                         ))
        
        # y-компонента
        v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
                         un[1:-1, 1:-1] * dtdx * (vn[1:-1, 1:-1] - vn[1:-1, :-2]) -
                         vn[1:-1, 1:-1] * dtdy * (vn[1:-1, 1:-1] - vn[:-2, 1:-1]) +
                         nu * dt * (
                             (vn[1:-1, 2:] - 2*vn[1:-1, 1:-1] + vn[1:-1, :-2]) / dx2 +
                             (vn[2:, 1:-1] - 2*vn[1:-1, 1:-1] + vn[:-2, 1:-1]) / dy2
                         ))
        
        # 2. Уравнение Пуассона для давления
        # Дискретный лапласиан
        b = np.zeros((ny, nx))
        b[1:-1, 1:-1] = rho * (
            (u[1:-1, 2:] - u[1:-1, 1:-1]) / dt / dx -
            (u[1:-1, 1:-1] - u[1:-1, :-2]) / dt / dx +
            (v[2:, 1:-1] - v[1:-1, 1:-1]) / dt / dy -
            (v[1:-1, 1:-1] - v[:-2, 1:-1]) / dt / dy
        )
        
        # Решаем уравнение Пуассона методом релаксации
        for _ in range(50):  # Итерации Якоби
            pn = p.copy()
            p[1:-1, 1:-1] = (
                (pn[1:-1, 2:] + pn[1:-1, :-2]) * dy2 +
                (pn[2:, 1:-1] + pn[:-2, 1:-1]) * dx2 -
                b[1:-1, 1:-1] * dx2 * dy2
            ) / (2 * (dx2 + dy2))
            
            # Граничные условия для давления (Неймана)
            p[:, -1] = p[:, -2]  # Правая граница
            p[0, :] = p[1, :]    # Верхняя граница
            p[:, 0] = p[:, 1]    # Левая граница
            p[-1, :] = p[-2, :]  # Нижняя граница
        
        # 3. Коррекция скорости с учетом давления
        u[1:-1, 1:-1] -= dt / rho * (p[1:-1, 1:-1] - p[1:-1, :-2]) / dx
        v[1:-1, 1:-1] -= dt / rho * (p[1:-1, 1:-1] - p[:-2, 1:-1]) / dy
        
        # Граничные условия для скорости (no-slip)
        u[0, :] = 0
        u[-1, :] = 0
        u[:, 0] = 0
        u[:, -1] = 0
        
        v[0, :] = 0
        v[-1, :] = 0
        v[:, 0] = 0
        v[:, -1] = 0
        
        # Входящий поток слева (например, постоянная скорость)
        u[:, 0] = 1.0
        
    return u, v, p

4 Добавляем крыло: почему это сложнее, чем кажется

До сих пор у нас была пустая область. Теперь поместим в нее профиль крыла. Самый простой вариант - прямоугольное препятствие. Но даже это создает проблемы.

Крыло - это внутренняя граница. Жидкость не может проходить сквозь него. Значит, на границе крыла нужно выставить условия no-slip (u=0, v=0). Но наша сетка не адаптирована под сложную геометрию. Прямоугольное крыло, расположенное по линиям сетки, - это единственный случай, который легко реализовать.

def add_wing(u, v, wing_x, wing_y, wing_w, wing_h):
    """Добавляем прямоугольное крыло в поле скорости"""
    # Преобразуем физические координаты в индексы сетки
    i_start = int(wing_x / dx)
    j_start = int(wing_y / dy)
    i_end = i_start + int(wing_w / dx)
    j_end = j_start + int(wing_h / dy)
    
    # Обнуляем скорость внутри крыла и на его границах
    u[j_start:j_end, i_start:i_end+1] = 0
    v[j_start:j_end+1, i_start:i_end] = 0
    
    # Также обнуляем скорость в непосредственной близости от крыла
    # для более стабильного решения
    u[j_start-1:j_end+1, i_start-1:i_end+2] = 0
    v[j_start-1:j_end+2, i_start-1:i_end+1] = 0
    
    return u, v

# Создаем крыло в центре области
wing_center_x = lx / 2
wing_center_y = ly / 2
wing_width = 0.3
wing_height = 0.1

u, v = add_wing(u, v, wing_center_x - wing_width/2, 
                wing_center_y - wing_height/2, wing_width, wing_height)

5 Визуализируем: без картинки это просто числа

Рассчитали поля - теперь нужно их увидеть. Иначе как понять, что все работает? Визуализация в CFD - это не просто красота. По картине течения опытный глаз сразу видит проблемы: области отрыва, нефизичные колебания, неустойчивости.

import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

def visualize(u, v, p, wing_x, wing_y, wing_w, wing_h):
    """Визуализация поля скорости и давления"""
    
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # Создаем сетку для визуализации
    x = np.linspace(0, lx, nx)
    y = np.linspace(0, ly, ny)
    X, Y = np.meshgrid(x, y)
    
    # Интерполируем скорости на одну сетку (для стрелок)
    # u и v на разных сетках, нужно усреднить
    u_center = 0.5 * (u[:, :-1] + u[:, 1:])
    v_center = 0.5 * (v[:-1, :] + v[1:, :])
    
    # 1. Поле давления с контурами
    contour = ax1.contourf(X, Y, p, levels=50, cmap='RdBu_r')
    plt.colorbar(contour, ax=ax1, label='Давление')
    
    # Добавляем крыло
    ax1.add_patch(Rectangle(
        (wing_x, wing_y), wing_w, wing_h,
        facecolor='gray', edgecolor='black', alpha=0.7
    ))
    
    # 2. Векторное поле скорости (каждые 5 точек для наглядности)
    skip = 5
    ax2.quiver(
        X[::skip, ::skip], Y[::skip, ::skip],
        u_center[::skip, ::skip], v_center[::skip, ::skip],
        scale=30, color='black'
    )
    
    # Линии тока
    ax2.streamplot(X, Y, u_center, v_center, color='blue', linewidth=0.5, density=2)
    
    # Крыло на втором графике
    ax2.add_patch(Rectangle(
        (wing_x, wing_y), wing_w, wing_h,
        facecolor='gray', edgecolor='black', alpha=0.7
    ))
    
    ax1.set_title('Поле давления')
    ax2.set_title('Поле скорости (векторы и линии тока)')
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    ax2.set_xlabel('x')
    ax2.set_ylabel('y')
    
    plt.tight_layout()
    plt.show()

# Запускаем симуляцию
u, v, p = solve_navier_stokes(u, v, p, dt, dx, dy, rho, nu, nt)

# Визуализируем результаты
visualize(u, v, p, 
          wing_center_x - wing_width/2, 
          wing_center_y - wing_height/2, 
          wing_width, wing_height)

Где собака зарыта: ошибки, которые сломают вашу симуляцию

Код работает? Отлично. Но работает ли он правильно? Вот список того, что может пойти не так даже при правильной реализации.

А что дальше? Куда развивать этот решатель

Рабочий прототип - это только начало. Если хочется двигаться дальше, вот направления, куда можно копать.

• Трехмерность - добавить ось z. Объем вычислений вырастет в кубе, но увидите настоящие вихревые шнуры.
• Турбулентность - для высоких чисел Рейнольдса нужно моделировать мелкомасштабные пульсации. LES или RANS модели.
• Подвижные границы - крыло, которое меняет угол атаки или форму. Методы погруженных границ или ALE.
• Параллельные вычисления - распараллелить циклы с помощью Numba или переписать на Cython. Или сразу использовать GPU через CuPy.

Главное, что теперь уравнения Навье-Стокса для вас - не страшные символы в учебнике, а конкретные строки кода, которые можно менять, ломать и чинить. Вы увидели, как давление и скорость танцуют друг с другом, как вязкость гасит вихри, как граничные условия управляют всем течением. Это понимание дороже любой готовой библиотеки.

Когда в следующий раз будете смотреть на полет самолета, вы уже будете знать: за подъемную силу крыла отвечает не какая-то абстрактная "аэродинамика", а конкретное распределение давления, которое можно посчитать вот такими вот циклами for и массивами NumPy. Магия превращается в ремесло. А ремесло - в инженерию.

Подписаться на канал