Гидродинамика на python: Пишем CFD симуляцию плоского течения несжимаемой жидкости

2026-08-22 в 4:06, admin, рубрики: python, вычислительная математика, гидродинамика, конечно-разностные методы, ландавшиц, механика сплошных сред, несжимаемая жидкость, симуляция, течения, уравнение навье-стокса

В этой статье мы решим классическую задачу из курса гидродинамики [4] : сделаем симуляцию плоского течения в канале. Для этого мы напишем на Python код, численно решающий систему дифференциальных уравнений в частных производных, описывающих поведение несжимаемой жидкости:

$frac{partial mathbf{u}}{partial t} + (mathbf{u} cdot nabla) mathbf{u}=-frac{1}{rho} nabla p + nu nabla^2 mathbf{u},quad nabla cdot mathbf{u}=0$

Первое уравнение- уравнение Навье-Стокса, а второе- уравнение неразрывности [1], [4].

Вывод этих уравнений можно прочитать в статье [15].

Применим дискретизацию и напишем конечно-разностную схему для симуляции (равномерная сетка 100 на 100) [18]:

Производная по x (центральная разность):

$left(frac{partial u}{partial x}right)_{i,j} approx frac{u_{i+1,j} - u_{i-1,j}}{2,Delta x}$

Вторая производная по x:

$left(frac{partial^2 u}{partial x^2}right)_{i,j} approx frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{(Delta x)^2}$

Дивергенция скорости:

$(nabla cdot mathbf{u})_{i,j} approx frac{u_{i+1,j} - u_{i-1,j}}{2,Delta x} + frac{v_{i,j+1} - v_{i,j-1}}{2,Delta y}$

Лапласиан скорости:

$(nabla^2 v)_{i,j} approx frac{v_{i+1,j} + v_{i-1,j} - 2v_{i,j}}{(Delta x)^2} + frac{v_{i,j+1} + v_{i,j-1} - 2v_{i,j}}{(Delta y)^2}$

Применим метод проекций [13], основанный на разложении Гельмгольца [14]:

$mathbf{u}^*=mathbf{u}^n - Delta t (mathbf{u}^n cdot nabla)mathbf{u}^n + Delta t,nu nabla^2mathbf{u}^n,quad nabla^2 p^{n+1}=frac{rho}{Delta t}nablacdotmathbf{u}^*,quad mathbf{u}^{n+1}=mathbf{u}^* - frac{Delta t}{rho}nabla p^{n+1}$

Здесь используется несколько чётких шагов:

Сначала система продвигается во времени до положения в середине временного шага, при этом решаются приведенные выше уравнения переноса массы и импульса с использованием подходящего метода адвекции. Этот этап называется прогнозирующим.

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

Затем выполняется корректирующая часть алгоритма. В ней используются центрированные по времени оценки скорости, плотности и т. д. для формирования состояния на конечном временном шаге.

Затем применяется окончательная проекция, обеспечивающая соблюдение ограничения на расходимость поля скоростей. Теперь система полностью обновлена в соответствии с новым временем.

Источник: [13]

Напишем на языке программирования python код, реализующий данную вычислительную схему:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import os

# --- Параметры задачи ---
L = 1.0
Nx, Ny = 100, 100
dx = L / (Nx - 1)
dy = L / (Ny - 1)
dt = 0.001
nu = 0.01
U_lid = 10.0

nt = 2000          # сколько шагов по времени (можно уменьшить для теста)
save_every = 10    # сохранять каждый N-й кадр (чтобы GIF не был огромным)

# --- Инициализация полей ---
u = np.zeros((Nx, Ny))
v = np.zeros((Nx, Ny))
p = np.zeros((Nx, Ny))

# Вспомогательные массивы
un = np.zeros_like(u)
vn = np.zeros_like(v)
pn = np.zeros_like(p)

# --- Функция решения уравнения Пуассона для давления ---
def solve_poisson(p, u, v, dx, dy, dt, nit=50):
    b = np.zeros_like(p)
    b[1:-1, 1:-1] = (1/dt) * ((u[2:, 1:-1] - u[:-2, 1:-1])/(2*dx) +
                              (v[1:-1, 2:] - v[1:-1, :-2])/(2*dy))

    for _ in range(nit):
        pn[:] = p
        p[1:-1, 1:-1] = ((pn[2:, 1:-1] + pn[:-2, 1:-1])*dy**2 +
                         (pn[1:-1, 2:] + pn[1:-1, :-2])*dx**2 -
                         b[1:-1, 1:-1]*dx**2*dy**2) / (2*(dx**2 + dy**2))
        p[:, 0] = 0; p[:, -1] = 0
        p[0, :] = 0; p[-1, :] = 0
    return p

# --- Подготовка сетки для визуализации ---
x = np.linspace(0, L, Nx)
y = np.linspace(0, L, Ny)
X, Y = np.meshgrid(x, y)

# --- Массив для хранения кадров (для GIF) ---
frames = []

print("Запуск симуляции и сбор кадров...")
for n in range(nt):
    un[:] = u
    vn[:] = v

    # Предсказание скорости (явная схема)
    u[1:-1, 1:-1] = (un[1:-1, 1:-1]
                     - un[1:-1, 1:-1] * (un[2:, 1:-1] - un[:-2, 1:-1])/(2*dx)
                     - vn[1:-1, 1:-1] * (un[1:-1, 2:] - un[1:-1, :-2])/(2*dy)
                     + nu * ((un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[:-2, 1:-1])/dx**2
                             + (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, :-2])/dy**2)) * dt

    v[1:-1, 1:-1] = (vn[1:-1, 1:-1]
                     - un[1:-1, 1:-1] * (vn[2:, 1:-1] - vn[:-2, 1:-1])/(2*dx)
                     - vn[1:-1, 1:-1] * (vn[1:-1, 2:] - vn[1:-1, :-2])/(2*dy)
                     + nu * ((vn[2:, 1:-1] - 2*vn[1:-1, 1:-1] + vn[:-2, 1:-1])/dx**2
                             + (vn[1:-1, 2:] - 2*vn[1:-1, 1:-1] + vn[1:-1, :-2])/dy**2)) * dt

    # Граничные условия для скорости
    u[:, 0] = 0; u[:, -1] = 0
    u[0, :] = 0; u[-1, :] = U_lid
    v[:, 0] = 0; v[:, -1] = 0
    v[0, :] = 0; v[-1, :] = 0

    # Давление
    p = solve_poisson(p, u, v, dx, dy, dt)

    # Коррекция скорости градиентом давления
    u[1:-1, 1:-1] -= (dt / 1.0) * (p[2:, 1:-1] - p[:-2, 1:-1]) / (2*dx)
    v[1:-1, 1:-1] -= (dt / 1.0) * (p[1:-1, 2:] - p[1:-1, :-2]) / (2*dy)

    # Повторное применение граничных условий
    u[:, 0] = 0; u[:, -1] = 0
    u[0, :] = 0; u[-1, :] = U_lid
    v[:, 0] = 0; v[:, -1] = 0
    v[0, :] = 0; v[-1, :] = 0

    # Сбор кадра для анимации
    if n % save_every == 0:
        speed = np.sqrt(u**2 + v**2)
        fig, ax = plt.subplots(figsize=(6, 5))
        # Фон: модуль скорости
        cf = ax.contourf(X, Y, speed.T, levels=40, cmap='viridis', alpha=0.7)
        # Стрелки: поле скорости (прорежем, чтобы не было каши)
        ax.quiver(X[::4, ::4], Y[::4, ::4],
                  u[::4, ::4].T, v[::4, ::4].T,
                  scale=25, headwidth=3, headlength=4, color='white', alpha=0.8, linewidth=0.4)
        ax.set_title(f'Шаг по времени: {n * dt:.3f} с')
        ax.set_xlabel('x (м)')
        ax.set_ylabel('y (м)')
        ax.axis('equal')
        ax.axis('off')
        plt.tight_layout(pad=0)
        fig.canvas.draw()
        image = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8')
        image = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))
        frames.append(image)
        plt.close(fig)

print(f"Собрано кадров: {len(frames)}")

# --- Сохранение в GIF ---
from PIL import Image

images = [Image.fromarray(frame) for frame in frames]
gif_path = 'navier_stokes_cavity.gif'
images[0].save(
    gif_path,
    save_all=True,
    append_images=images[1:],
    duration=100,  # мс на кадр
    loop=0
)
print(f"GIF сохранён: {os.path.abspath(gif_path)}")

Видна завихрённость течения и кривое поле скоростей. Действительно, посчитав число Рейнольдса [16] для данного течения, получим 10000, что значительно больше критического значения, значит оно является турбулентным.

Стоит отметить вычислительную неустойчивость данной схемы: если мы меняем шаг по времени, то полученная гифка будет довольно сильно отличаться от этой.

Например, делаем шаг по времени 0.0019:

При шаге 0.002 численная схема разлетается:

Теперь проследим, как изменение параметров жидкости и сетки влияет на полученную визуализацию. Например, зададим такие параметры:

# --- Параметры задачи ---
L = 1.0
Nx, Ny = 20, 150
dx = L / (Nx - 1)
dy = L / (Ny - 1)
dt = 0.01
nu = 0.001
U_lid = 1

nt = 1000          # сколько шагов по времени (можно уменьшить для теста)
save_every = 10    # сохранять каждый N-й кадр (чтобы GIF не был огромным)

То есть мы сделали жидкость слишком жидкой (простите за тавтологию), уменьшив вязкость в 10 раз, и ужали сетку до 20 на 150 (сделали большую детализацию по оси y и маленькую по оси x).

Следует отметить, что начальную скорость мы уменьшили тоже в 10 раз, поэтому число Рейнольдса [16] осталось неизменным, как и характер течения.

Результат видно на экране:

Видно, что вихрей стало значительно больше, наблюдается турбулентность.

Чтобы всё увидеть ещё лучше, пойдём на крайние меры: сузим сетку по X до 5 и растянем до 200 по Y. То есть зададим следующие параметры:

# --- Параметры задачи ---
L = 2.0
Nx, Ny = 5, 200
dx = L / (Nx - 1)
dy = L / (Ny - 1)
dt = 0.0001
nu = 0.05
U_lid = 10

nt = 1000          # сколько шагов по времени (можно уменьшить для теста)
save_every = 10    # сохранять каждый N-й кадр (чтобы GIF не был огромным)

Получим следующую красивую анимацию:

Вывод:таким образом, в данной статье получена симуляция на Python течения несжимаемой жидкости: при помощи конечно-разностной схемы численно решено уравнение Навье-Стокса [1], [4], [15] на плоскости.

Написав код, мы меняли параметры системы: меняли временной шаг, изменяли вязкость среды, корректировали сетку и смотрели, как при этом меняется финальная анимация. Предоставляю читателям, которые интересуются гидродинамикой, возможность подвигать параметры самостоятельно и получить конструкции, которых нет в данной статье.

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

Примечание: в 3d моделировать намного сложнее по следующим причинам:

Размер сетки. В 2D у нас 100 на 100= 10000 ячеек, а в 3D при той же сетке 1000000,

а при нормальной детализации намного больше.
Время счёта. Каждый шаг по времени становится в разы тяжелее: больше операций, больше памяти, медленнее сходимость.
Устойчивость. Условие CFD [17] в 3D жёстче: шаг по времени dt придётся уменьшать, и симуляция будет идти долго даже на хорошей машине.
Визуализация. В 3D «просто стрелки» уже не работают: нужны изоповерхности, объёмные рендеры, векторные поля с прорежением.

Поэтому уравнение Навье-Стокса в 3D моделируют с использованием специальных CFD-симуляторов (ANSYS, OpenFOAM, COMSOL и другие) [17] и на мощных суперкомпьютерах.

При этом часто применяют метод конечных объёмов и метод конечных элементов.

Тем не менее у меня получилось сделать псевдо-симуляцию в 3D на обычном компьютере:

Следует добавить, что вопрос о существовании и единственности решения уравнения Навье-Стокса в трёхмерном пространстве является открытой математической проблемой и входит в список нерешённых Задач тысячелетия [11], за решение которых Математический институт Клэя [12] выплатит премию в миллион долларов.

Для двумерного потока задача о существовании и единственности решения уравнения Навье-Стокса решена положительно. В 1969 году советский и российский математик Ольга Александровна Ладыженская [10] предоставила книгу «Математические вопросы динамики вязкой несжимаемой жидкости» [9] ,в которой доказана глобальная однозначная разрешимость (то есть существование и гладкость решений на любом промежутке времени) начально-краевой задачи для двумерной системы уравнений Навье — Стокса [1], [4], [15].

Статья носит чисто демонстрационный и иллюстративный характер.

Автор не имеет специального физико-математического образования и не является профессиональным математиком!