Пишем CFD solver для симуляции потока воздуха (часть 1)

Начну издалека. Пару недель назад мы с моим другом решили сделать один очень крутой проект, на мою долю выпала часть аэродинамических подсчетов. Я стал углубляться в эту тему и понял, что на самом деле все не так уж и просто. В зависимости от угла атаки (угла между вектором скорости и нормали), а также от формы тела меняются многие параметры, которые влияют не только непосредственно на тело, но и на окружающий воздух, его давление, температуру, скорость и ‘очень много всего еще’, что в дальнейшем будет снова влиять на наше тело. Получается замкнутый круг из гигантского числа диффуров, который еще можно до бесконечности усложнять, используя все более и более сложные модели. В конце концов я сдался и у меня был выбор: пойти просто посчитать, используя готовые программы, или потратить кучу времени, разобраться во всем этом и под конец сделать все самому. Угадайте что я выбрал? Конечно же второе.

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

Для начала стоит понимать, что конкретно сложного в такой задаче, как симуляция потоков (будь то воды или воздуха). В основном — это количество частиц, которые влияют не только на объект, но и сами на себя. Из — за их количества даже появилось отдельное направление в физике — МКТ (Молекулярно‑кинетическая теория). В моем понимании, она нужна только для того, что не писать 2 закон Ньютона для каждой из частиц, а просто показать, что есть такие параметры как давление, температура и объем, благодаря которым можно усредненно посчитать, то, как газ будет себя вести.

Но тогда встает логичный вопрос, каким образом нам посчитать динамику абсолютно хаотичных потоков. В симуляциях любят использовать сетки, то есть набор упорядоченных точек, в каждую из которых мы будем аппроксимировать какой‑то объем (площадь) и уже писать характеристики для него. Чем больше таких точек у нас будет, тем точнее у нас получится симуляция.

Отлично, у нас уже есть модель, которая позволяет нам ввести количественное описание. Так давайте воспользуемся этим.

Уравнение Навье‑Стокса (без вязкости):

Где:

• V — общая скорость

• (V·∇)V — конвекция

• ρ — плотность воздуха (в нашем случае константа)

• ∇p — градиент давления

Не пугайтесь, здесь все проще, чем могло показаться. Для начала скажем, что изменение общей скорости по времени зависит и от самой скорости, и от давления, тогда напишем, как будет изменяться скорость со временем, только при условии, что на нее влияет сама скорость:

Давайте теперь распишем правую часть, чтобы она не была такой страшной:

при:

где:

• u — проекция скорости на OX

• v — проекция скорости на OY

получим, что cкалярное произведение V и ∇ — это сумма произведений компонент:

тогда:

(V·∇)V = (u·∂u/∂x + v·∂u/∂y, u·∂v/∂x + v·∂v/∂y)

получим, что для горизонтальной оси (см. уравнение Навье‑Стокса):

а для вертикальной:

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

• полностью готовый код можно найти по ссылке: https://github.com/PerKyyling/AirSimulation

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

Теперь от физики давайте перейдем к проге. Для начала простой шаблон python‑скрипта на фреймворке arcade (аналог pygame, но полегче):

import arcadeWIDTH = 1000HEIGHT = 1000TITLE = "Simulation"CELL_SIZE = 10WHITE = (255, 255, 255)LITE_BLUE = (200, 230, 255)# создаем дочерний класс от arcade.Windowclass Window(arcade.Window):    # клонируем init класса    def __init__(self, width, height, title):        super().__init__(width, height, title)    # on_update вызывается каждый раз между отрисовкой    def on_update(self, delta_time):      pass    # метод для рисования    def on_draw(self):      self.clear()def main():    sm = Window(WIDTH, HEIGHT, TITLE)    arcade.run()if __name__ == "__main__":    main()

Теперь давайте предположим, что наше тело — квадрат, оно движется со скоростью V = 10 м/с. Воздух покоится.
Тогда, если связать систему отсчета, связанную с квадратом, то получим, что воздух имеет горизонтальную скорость 10 м/с. Тогда давайте создадим два numpy‑массива, в которые положим горизонтальную скорость (u) и вертикальную скорость (v):

import arcadeimport numpyWIDTH = 1000HEIGHT = 1000TITLE = "Simulation"CELL_SIZE = 10WHITE = (255, 255, 255)LITE_BLUE = (200, 230, 255)class Window(arcade.Window):    def __init__(self, width, height, title):        super().__init__(width, height, title)        self.background_color = arcade.color.WHITE        self.N = 100        self.V = 10  # м * с^-1        self.u_old = np.full((self.N, self.N), self.V, dtype=np.float32)        self.v_old = np.zeros((self.N, self.N), dtype=np.float32)            def on_update(self, delta_time):      pass    def on_draw(self):      self.clear()def main():    sm = Window(WIDTH, HEIGHT, TITLE)    arcade.run()if __name__ == "__main__":    main()

u_old — массив с горизонтальными скоростями (в момент времени t = 0 все u[i][j] = V)

v_old — массив с вертикальными скоростями, поток горизонтален => все 0

Теперь у вас возникнет вопрос, что такое N и откуда индексы old, начну отвечать по порядку.

Для начала давайте разберемся с константами в init, чтобы к ним не возвращаться:

    def __init__(self, width, height, title):        super().__init__(width, height, title)        self.background_color = arcade.color.WHITE        self.N = 100        self.V = 10  # м * с^-1        self.dx = 0.001        self.dy = 0.001        self.dt = 0.0000625        self.u_old = np.full((self.N, self.N), self.V, dtype=np.float32)        self.u_new = np.zeros((self.N, self.N), dtype=np.float32)        self.v_old = np.zeros((self.N, self.N), dtype=np.float32)        self.v_new = np.zeros((self.N, self.N), dtype=np.float32)        # массив для объекта        self.collision = np.zeros((self.N, self.N), dtype=np.float32)        self.p = np.zeros((self.N, self.N), dtype=np.float32)        self.p_old = np.zeros((self.N, self.N), dtype=np.float32)        #сколько времени прошло        self.time_total = 0        # границы объекта        self.left_border = 40   # X_min        self.right_border = 60  # X_max        self.low_border = 40    # Y_min        self.up_border = 60     # Y_max        self.rho = 1.225 / 10 ** 6        # пока нам не нужно        self.epsilon = 1e-3        self.max_disp = 0        # сразу добавляем 1 в места, где есть объект и сбрасываем там все скорости        for j in range(self.low_border, self.up_border):          for i in range(self.left_border, self.right_border):              self.v_old[j][i] = 0              self.u_old[j][i] = 0              self.collision[j][i] = 1        # в дальнейшем проверка, есть ли обьект в клетке (i, j), будет выглядить:        # if self.collision[i][j] == 1

Мы симулируем 1м^3, расстояние между точками — dx, dy (1 см), промежуток времени dt (625 * 10^-7), p — массив давлений, о нем позже, rho — плотность воздуха.

Теперь стоит обсудить, что такое индексы old и new, наша симуляция будет работать очень просто: условно, в момент времени t = t0, мы имели все массивы с индексов old, потом прошло dt и наши массивы по какой то логике поменялись на new в момент времени t = t1. Тогда новый массив после каждой итерации мы записываем как старый и после второго dt делаем то же самое, что и после первого, так мы записываем новые значения в каждый момент времени, зная старый.

Теперь перейдем к реализации конвекции. Для начала разложим наше уравнение для конвекции в численные методы. (i — строка, j — столбец)

вспомним уравнение по оси OX:

его можно разложить в:

то же самое мы можем сделать для вертикальной скорости:

откуда легко можно выразить и

давайте добавим 2 ключевые функции в наш код:

    # после каждой итерации скорость будет проникать сквозь квадрат, ее нужно будет сбрасывать    def zeroingCollisionSpeed(self):        for j in range(self.low_border, self.up_border):            for i in range(self.left_border, self.right_border):                self.u_old[j][i] = 0                self.v_old[j][i] = 0        def speedTransmission(self):        # для стабильности, при  отладке вроде без этого улетало        for j in range(self.low_border, self.up_border):            for i in range(self.left_border, self.right_border):                self.u_new[j][i] = 0                self.v_new[j][i] = 0        # основной цикл для нахождения скорости в след. момент времени        for j in range(1, self.N - 1):            for i in range(1, self.N - 1):                if self.collision[j][i] == 1:                    continue                if self.u_old[j][i] > 0:                    dudx = (self.u_old[j][i] - self.u_old[j][i-1]) / self.dx                else:                    dudx = (self.u_old[j][i+1] - self.u_old[j][i]) / self.dx                if self.v_old[j][i] > 0:                    dudy = (self.u_old[j][i] - self.u_old[j-1][i]) / self.dy                else:                    dudy = (self.u_old[j+1][i] - self.u_old[j][i]) / self.dy                self.u_new[j][i] = self.u_old[j][i] - self.dt * (                        self.u_old[j][i] * dudx + self.v_old[j][i] * dudy)                if self.u_old[j][i] > 0:                    dvdx = (self.v_old[j][i] - self.v_old[j][i-1]) / self.dx                else:                    dvdx = (self.v_old[j][i+1] - self.v_old[j][i]) / self.dx                if self.v_old[j][i] > 0:                    dvdy = (self.v_old[j][i] - self.v_old[j-1][i]) / self.dy                else:                    dvdy = (self.v_old[j+1][i] - self.v_old[j][i]) / self.dy                self.v_new[j][i] = self.v_old[j][i] - self.dt * (                        self.u_old[j][i] * dvdx + self.v_old[j][i] * dvdy)        # граничные условия - воздух влетает слева        self.u_new[:, 0] = self.V        self.v_new[:, 0] = 0        # и уравниваем скорости по остальным границам        self.u_new[:, self.N-1] = self.u_new[:, self.N-2]        self.v_new[:, self.N-1] = self.v_new[:, self.N-2]        self.u_new[0, :] = self.u_new[1, :]        self.v_new[0, :] = 0        self.u_new[self.N-1, :] = self.u_new[self.N-2, :]        self.v_new[self.N-1, :] = 0                # снова сбрасываем скорость внутри квадрата        for j in range(self.low_border, self.up_border):            for i in range(self.left_border, self.right_border):                self.u_new[j][i] = 0                self.v_new[j][i] = 0

Тогда функция on_update будет выглядеть так:

    def on_update(self, delta_time):        self.zeroingCollisionSpeed()        self.speedTransmission()        self.v_old = self.v_new.copy()        self.u_old = self.u_new.copy()

добавим функцию для рисования:

    def on_draw(self):        self.clear()        # находим вектор общей скорости в данных точках        max_speed = np.max(np.sqrt(self.u_new ** 2 + self.v_new ** 2))        if max_speed == 0:            max_speed = 1.0        # т.к. мы симулируем 1 метр, делим его на 100 частей, а в длину программа 1000px        for j in range(self.N):            for i in range(self.N):                x = i * CELL_SIZE                y = j * CELL_SIZE                # если в этой точке объект - рисуй белым                if self.collision[j][i] == 1:                    arcade.draw_lbwh_rectangle_filled(x, y, CELL_SIZE - 1, CELL_SIZE - 1, WHITE)                else:                    # иначе рисуй скорость                    speed = np.sqrt(self.u_new[j][i] ** 2 + self.v_new[j][i] ** 2)                    norm = min(speed / max_speed, 1.0)                    # Синий — медленно, красный — быстро                    color = (int(norm * 255), 0, int((1 - norm) * 255))                    arcade.draw_lbwh_rectangle_filled(x, y, CELL_SIZE - 1, CELL_SIZE - 1, color)

если все успешно совместить, получится:

Теперь добавим давление. Для начала давайте вспомним основное уравнение Навье‑Стокса:

Мы уже знаем относительно момента времени t, какой будет скорость в момент времени t + dt без учета давления. Мы видим, что давление как‑то зависит от скорости, но нам надо выразить скорость через давление, возникает парадокс: одно уравнение — две неизвестных. Но у нас есть один замечательный физических факт, воздух движется из области высокого давление в область низкого, а также есть наша договоренность, что воздух не сжимаемый, то есть его дивергенция равна 0. Теперь давайте вспомним, из‑за чего воздух должен начать изменять свое направление движения. Из‑за препятствия. Внутри препятствия скорость 0, а снаружи кубика V, то есть его дивергенция не 0, ведь создается гигантский градиент. Теперь давайте запишем еще одно очень важное уравнение. Уравнение Пуассона (выводится математически из Навье‑Стокса):

где u* — промежуточная скорость (которая у нас есть), тогда, если:

что в численных:

мы уже в численных методах считали. Но после всего этого остается маленькая проблема, мы знаем u* в каждой точке, но для выражения p[i][j] нам нужны значения p, в четырех соседних точках. Поздравляю, уравнение с 2 неизвестными мы свели к уравнению с 5, но в этом и прелесть. У нас все неизвестные лежат в одной матрице. Допустим, мы ее заполним нулями, тогда при давлении 0 дивергенция скорости в точках не 0, а должна быть 0, тогда мы должны дать нашим скоростям изменить давление так, что:

где epsilon () — еще одна очень маленькая константа в нашем коде. Стоит сказать, почему нам ее нужно вводить, ведь можно просто приравнять к 0, оказывается нельзя. Как я сказал раньше, у нас неизвестная, зависящая от 4 других, что дает 10 000 уравнений. Такую штуку можно решить методом Гаусса‑Зейделя. Кучу раз пробежаться по массиву со случайными величинами и все ближе и ближе приближать их к истинному значению, чем меньше приближение, тем больше времени будет исполняться алгоритм.

Зная все это, давайте исправим on_update и добавим две очень важные функции:

on_update теперь будет выглядеть так:

    def on_update(self, delta_time):        self.zeroingCollisionSpeed()        self.speedTransmission()        iteration = 0 # кол-во итераций        # показывает максимальную разницу от 0        self.max_disp = 1.0 # с дисперсией не имеет ничего общего, но пусть так        # перебираем пустой массив пока не сделаем много итераций или не получим достаточную точность        while abs(self.max_disp) > self.epsilon and iteration < 100:            # функция для уточнения давлений            self.updatePressure()            iteration += 1        # функция для финального исправления скорости в след. момент времени        self.finalUpdateSpeed()        self.v_old = self.v_new.copy()        self.u_old = self.u_new.copy()

напишем updatePressure():

    # буквально создаем промежуточное давление в каждой точке, выражая его из Пуассона    def updatePressure(self):        self.max_disp = 0        for j in range(1, self.N - 1):            for i in range(1, self.N - 1):                if self.collision[j][i] == 0:                    dudx = (self.u_new[j][i+1] - self.u_new[j][i-1]) / (2 * self.dx)                    dvdy = (self.v_new[j+1][i] - self.v_new[j-1][i]) / (2 * self.dy)                    div = dudx + dvdy                    self.p[j][i] = (self.dx**2 * self.dy**2 * (self.rho / self.dt) * div -                                    (self.p[j][i+1] + self.p[j][i-1]) * self.dy**2 -                                    (self.p[j+1][i] + self.p[j-1][i]) * self.dx**2) / \                                   (-2 * (self.dx**2 + self.dy**2))                    # находим отклонение давления в каждой точке относительно предыдущего                    self.max_disp += abs(self.p[j][i] - self.p_old[j][i])        # граничные условия: уравниваем давление по краям        self.p[:, 0] = self.p[:, 1]        self.p[:, self.N-1] = self.p[:, self.N-2]        self.p[0, :] = self.p[1, :]        self.p[self.N-1, :] = self.p[self.N-2, :]        # обнуляем давление внутри квадрата        for j in range(self.low_border, self.up_border):            for i in range(self.left_border, self.right_border):                self.p[j][i] = 0        self.p_old = self.p.copy()

Ну и наконец заканчиваем нашу симуляцию финальным изменением скорости из‑за давления:

  # раньше в коде мы выражали конвекцию из Навье-Стокса, сейчас же просто дописываем его до конца  def finalUpdateSpeed(self):        for j in range(1, self.N - 1):            for i in range(1, self.N - 1):                if self.collision[j][i] == 1:                    continue                self.u_new[j][i] -= (self.dt / self.rho) * (self.p[j][i+1] - self.p[j][i-1]) / (2 * self.dx)                self.v_new[j][i] -= (self.dt / self.rho) * (self.p[j+1][i] - self.p[j-1][i]) / (2 * self.dy)        for j in range(self.low_border, self.up_border):            for i in range(self.left_border, self.right_border):                self.u_new[j][i] = 0                self.v_new[j][i] = 0        # ограничение скорости стоит сделать по 2 вещам:        # 1) нету вязкости, что может добавить резкости при больших перепадах        # 2) чтобы симуляция не взлетела (число Куранта, почитайте, если будет время)        max_allowed = self.V * 5        self.u_new = np.clip(self.u_new, -max_allowed, max_allowed)        self.v_new = np.clip(self.v_new, -max_allowed, max_allowed)

Поздравляю, у нас получилась ОЧЕНЬ медленная симуляция на питоне, а, раз так, то вычисление можно перенести на C++ (это в репозитории, я про логику рассказываю), после всех танцев с бубном получим вот такую картинку:

В заключение, скажу, почему 2д. Как оказалось, для 10^4 точек скорость вычислений уже очень низкая, поэтому в этой части я показал прототип, в который всегда можно добавить любую другую форму объекта и ось, но все равно придется менять язык. В следующей части я найду на чем можно написать 3д солвер (если есть идеи, буду рад их прочитать) и полностью перенесу логику на плюсы.