Вычислительная гидродинамика на Python: моделируем 2D-течение несжимаемой жидкости

4 часа назад Прокомментировать Просмотры: 10

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

Первое соотношение представляет собой уравнение Навье — Стокса (закон сохранения импульса), второе — уравнение неразрывности [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^2\mathbf{u}^n,\quad \nabla^2 p^{n+1} = \frac{\rho}{\Delta t}\nabla\cdot\mathbf{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])/(2dx) +
(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] для заданной конфигурации достигает 10 000 — это многократно превышает критический порог и переводит течение в турбулентный режим.
При этом явные конечно-разностные схемы обладают строгими ограничениями по устойчивости: малейшая корректировка шага интегрирования по времени существенно трансформирует итоговую динамику.
Рассмотрим, к примеру, поведение потока при увеличении шага по времени до dt = 0.0019:
Шаг по времени 0.0019
При шаге dt = 0.002 критерий Куранта нарушается, и численная схема катастрофически расходится:
шаг по времени 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 не был огромным)
В данном тесте кинематическая вязкость снижена на порядок, а сетка задана резко анизотропной (20 узлов по горизонтали против 150 по вертикали).
При этом граничная скорость движения крышки также уменьшена в 10 раз, благодаря чему инвариантное число Рейнольдса [16] и общий гидродинамический режим остаются неизменными.
Эволюция потока представлена ниже:
Уменьшаем вязкость
Снижение диссипации импульса приводит к росту мелкомасштабных вихрей и выраженным нестационарным пульсациям.
Для демонстрации крайнего случая анизотропии предельно сожмем сетку вдоль оси X (всего 5 ячеек) и детализируем по оси Y (200 ячеек):
# --- Параметры задачи ---
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) моделированию связан с рядом серьезных сложностей:


Размерность расчетной сетки. Сетка 100 × 100 в 2D содержит 10 000 узлов, тогда как эквивалентная трехмерная модель требует уже 1 000 000 ячеек, а качественная дискретизация — десятки миллионов.


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


Условия устойчивости (CFL). Критерии Куранта [17] в 3D требуют предельно малого временного шага dt, что делает прямой расчет ресурсоемким даже на мощных рабочих станциях.


Рендеринг и постпроцессинг. Вместо плоских векторных диаграмм требуются алгоритмы построения изоповерхностей, объемный рендеринг (volumetric rendering) и методы прореживания векторных полей.
По этой причине промышленное моделирование 3D гидродинамики выполняется в специализированных CFD-пакетах (ANSYS Fluent, OpenFOAM, COMSOL и др.) [17] на базе высокопроизводительных кластеров.
В таких системах преимущественно используются метод конечных объемов (FVM) и метод конечных элементов (FEM).
Тем не менее, на обычном ПК удалось получить базовую псевдо-3D визуализацию:


Псевдо-3D
Важно подчеркнуть, что фундаментальный вопрос о глобальном существовании и гладкости решений уравнений Навье — Стокса в трехмерном пространстве остается нерешенным и входит в список семи Задач тысячелетия [11], учрежденных Математическим институтом Клэя [12] с наградой в $1 000 000.
Для двумерного случая корректность постановки строго доказана. В 1969 году выдающийся советский математик Ольга Александровна Ладыженская [10] в фундаментальном труде «Математические вопросы динамики вязкой несжимаемой жидкости» [9] установила глобальную однозначную разрешимость начально-краевой задачи для двумерной системы уравнений Навье — Стокса [1], [4], [15].
Статья носит научно-популярный и демонстрационный характер.
Автор делится инженерным опытом реализации и не претендует на статус академического эксперта в области теоретической физики.

Если подобные задачи появятся в школьных экзаменах — это станет серьезным испытанием для нервной системы любого абитуриента!

                                                                        Литература:
1.https://ru.wikipedia.org/wiki/Уравнения_Навье_—_Стокса
2.https://cyberleninka.ru/article/n/metod-postroeniya-resheniy-uravneniy-navie-stoksa
3.https://cyberleninka.ru/article/n/ob-ustanovivshihsya-resheniyah-uravneniya-navie-stoksa
4.https://djvu.online/file/8VpV5lDf4CKnl (Ландау и Лишпиц 6 том Гидродинамика).
5.https://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=zvmmf&paperid=9339&option_lang=rus&ysclid=mt14ysrr17927377606
6.Роуч П. Вычислительная гидродинамика. М.: Мир, 1980.
7.http://www.unn.ru/pages/issues/vestnik/99999999_West_2013_1(3)/47.pdf?ysclid=mt150iq918618064055
8.https://cfd-education.ru/wp-content/uploads/2026/05/BOOK_rus_final.pdf
9.https://reallib.org/reader?file=505102&pg=22
10.https://ru.wikipedia.org/wiki/Ладыженская,_Ольга_Александровна
11.https://ru.wikipedia.org/wiki/Задачи_тысячелетия
12.https://ru.wikipedia.org/wiki/Математический_институт_Клэя
13.https://en.wikipedia.org/wiki/Projection_method_(fluid_dynamics)
14.https://en.wikipedia.org/wiki/Helmholtz_decomposition
15.https://habr.com/ru/articles/171327/
16.https://ru.wikipedia.org/wiki/Число_Рейнольдса
17.https://ru.wikipedia.org/wiki/Вычислительная_гидродинамика
18.https://ru.wikipedia.org/wiki/Метод_прямоугольников#Составные_формулы_для_равномерных_сеток




 
Источник