Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Дифференциальные уравнения движения вязкой жидкости (уравнения Навье – Стокса)

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

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

В вязкой жидкости различают два рода напряжений.

  • 1. Нормальное напряжение Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса– проекция Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксана нормаль п в данной точке поверхности.
  • 2. Касательное напряжение Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса– проекция Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксана касательную плоскость к поверхности в данной точке. Касательные напряжения имеют место лишь при движении вязкой жидкости.

Схема поверхностных сил, действующих в вязкой жидкости, представлена на рис. 4.21. Первый индекс при р указывает нормаль к площадке, на которую действует напряжение, второй – ось, на которую оно спроектировано.

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Рис. 4.21. Схема напряжений, возникающих в вязкой жидкости

Выделим в движущейся жидкости элементарный параллелепипед с ребрами, параллельными осям х, у, z, и рассмотрим поверхностные силы, действующие на сто гранях (рис. 4.22).

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Рис. 4.22. Схема сил, действующих на элементарный параллелепипед

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

В вязкой жидкости, в противоположность невязкой, напряжение зависит от ориентации площадки в данной точке. Однако, как строго доказывается в технической гидромеханике, сумма всех нормальных напряжений Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксав данной точке не зависит от ориентации площадки и, следовательно, эта сумма является скалярной функцией только координат точки и времени, в связи с чем вводится новое понятие о гидромеханическом давлении.

Гидромеханическим давлением в вязкой жидкости называют давление, величина которого равна среднему арифметическому из величин любых трех нормальных напряжений в данной точке:

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Знак «минус» берется потому, что рxx, руу, рzz, направленные по внешней нормали, всегда отрицательны, а р должно быть положительным, как это обычно принимают в гидравлике. Таким образом, понятия гидромеханического давления в вязкой жидкости и гидродинамического давления в невязкой идеальной жидкости существенно различны.

Дадим упрощенный вывод уравнений движения вязкой жидкости применительно лишь к частному случаю несжимаемой жидкости. Рассмотрим вначале одномерное движение жидкости в направлении, параллельном оси Ох.

Выделим в потоке движущейся жидкости элементарный параллелепипед с ребрами dxy dy, dz (рис. 4.23). Соотношение (4.16) для сил, действующих на элементарный параллелепипед, в данном случае будет

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса(4.17)

где Fтр – сила трения, определяемая по формуле

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Рис. 4.23. Схема сил к выводу уравнений Навье – Стокса

По закону Ньютона для касательного напряжения т имеем

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

где Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса– динамическая вязкость.

Отсюда сила трения будет равна

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Формулы для силУравнения движения вязкой несжимаемой жидкости уравнения навье стоксабыли указаны в параграфе 4.8. Подставляя их и выражение для силы трения в формулу (4.17), получим

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

где Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

В общем случае движения в трехмерном пространстве, когда vx изменяется по всем направлениям, а не только в направлении оси z, проекция силы трения на ось х определится более сложным выражением

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Тогда уравнение движения в проекции на ось х будет

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Для всех трех осей х, у, z получим в развернутом виде

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Последние три уравнения представляют систему дифференциальных уравнений Навье – Стокса для вязкой несжимаемой жидкости.

В векторной форме можно записать

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

i, j, k – орты координатных осей (см. параграф 3.3).

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

Добавим к полученным уравнениям движения уравнение неразрывности для несжимаемой жидкости

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Полагая, что массовые силы X, Y, Z заданы, получим замкнутую систему четырех уравнений с четырьмя неизвестными функциямиУравнения движения вязкой несжимаемой жидкости уравнения навье стокса

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

Видео:Science show. Выпуск 51. Уравнение Навье - СтоксаСкачать

Science show. Выпуск 51. Уравнение Навье - Стокса

Уравнение Навье-Стокса и симуляция жидкостей на CUDA

Привет, Хабр. В этой статье мы разберемся с уравнением Навье-Стокса для несжимаемой жидкости, численно его решим и сделаем красивую симуляцию, работающую за счет параллельного вычисления на CUDA. Основная цель — показать, как можно применить математику, лежащую в основе уравнения, на практике при решении задачи моделирования жидкостей и газов.

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Видео:Уравнение Навье-Стокса на пальцах. МЛФ#2Скачать

Уравнение Навье-Стокса на пальцах. МЛФ#2

Уравнение Навье-Стокса для несжимаемой жидкости

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

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

Условно уравнение Навье-Стокса можно разделить на пять частей:

  • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— обозначает скорость изменения скорости жидкости в точке (его мы и будем считать для каждой частицы в нашей симуляции).
  • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— перемещение жидкости в пространстве.
  • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— давление, оказываемое на частицу (здесь Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— коэффициент плотности жидкости).
  • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— вязкость среды (чем она больше, тем сильнее жидкость сопротивляется силе, применяемой к ее части), Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— коэффициент вязкости).
  • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— внешние силы, которые мы применяем к жидкости (в нашем случае сила будет играть вполне конкретную роль — она будет отражать действия, совершаемые пользователем.

Также, так как мы будем рассматривать случай несжимаемой и однородной жидкости, мы имеем еще одно уравнение: Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса. Энергия в среде постоянна, никуда не уходит, ниоткуда не приходит.

Будет неправильно обделить всех читателей, которые не знакомы с векторным анализом, поэтому заодно и бегло пройдемся по всем операторам, присутствующим в уравнении (однако, настоятельно рекомендую вспомнить, что такое производная, дифференциал и вектор, так как они лежат в основе всего того, о чем пойдет речь ниже).

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

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Оператор набла представляет из себя векторный дифференциальный оператор и может быть применен как к скалярной функции, так и к векторной. В случае скаляра мы получаем градиент функции (вектор ее частных производных), а в случае вектора — сумму частых производных по осям. Главная особенность данного оператора в том, что через него можно выразить основные операции векторного анализа — grad (градиент), div (дивергенция), rot (ротор) и Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса(оператор Лапласа). Стоит сразу же отметить, что выражение Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксане равносильно Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— оператор набла не обладает коммутативностью.

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

Видео:Вывод уравнений Навье-Стокса - Лекция 3Скачать

Вывод уравнений Навье-Стокса - Лекция 3

Численное решение уравнения Навье-Стокса

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

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

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

Теперь выразим градиент, дивергенцию и оператор Лапласа с учетом нашей координатной сетки ( Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— индексы в массиве, Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— взятие соответствующих компонентов у вектора):

ОператорОпределениеДискретный аналог
gradУравнения движения вязкой несжимаемой жидкости уравнения навье стоксаУравнения движения вязкой несжимаемой жидкости уравнения навье стокса
divУравнения движения вязкой несжимаемой жидкости уравнения навье стокса Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса
Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксаУравнения движения вязкой несжимаемой жидкости уравнения навье стоксаУравнения движения вязкой несжимаемой жидкости уравнения навье стокса
rotУравнения движения вязкой несжимаемой жидкости уравнения навье стоксаУравнения движения вязкой несжимаемой жидкости уравнения навье стокса

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

Перемещение частиц

Данные утверждения работают только в том случае, если мы можем найти ближайшие частицы относительно рассматриваемой на данный момент. Чтобы свести на нет все возможные издержки, связанные с поиском таковых, мы будет отслеживать не их перемещение, а то, откуда приходят частицы в начале итерации путем проекции траектории движения назад во времени (проще говоря, вычитать вектор скорости, помноженный на изменение времени, из текущей позиции). Используя этот прием для каждого элемента массива, мы будем точно уверены, что у любой частицы будут «соседи»:

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Положив, что Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— элемент массива, хранящий состояния частицы, получаем следующую формулу для вычисления ее состояния через время Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса(мы полагаем, что все необходимые параметры в виде ускорения и давления уже рассчитаны):

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Заметим сразу же, что при достаточно малом Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксаи скорости мы можем так и не выйти за пределы ячейки, поэтому очень важно правильно подобрать ту силу импульса, которую пользователь будет придавать частицам.

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

Вязкость

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

. В таком случае итеративное уравнение для скорости примет следующий вид:

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Мы несколько преобразуем данное равенство, приведя его к виду Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса(стандартный вид системы линейных уравнений):

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

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

Внешние силы

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Импульс-вектор легко посчитать как разность между предыдущей позицией мыши и текущей (если такая имелась), и здесь как раз-таки можно проявить креативность. Именно в этой части алгоритма мы можем внедрить добавление цветов в жидкость, ее подсветку и т.п. К внешним силам также можно отнести гравитацию и температуру, и хоть реализовать такие параметры несложно, в данной статье рассматривать их мы не будем.

Давление

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

Назовем Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксавекторное поле, полученное после расчета перемещения, внешних сил и вязкости. Оно будет иметь ненулевую дивергенцию, что противоречит условию несжимаемости жидкости (Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса), и чтобы это исправить, необходимо рассчитать давление. Согласно теореме разложения Гельмгольца, Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксаможно представить как сумму двух полей:

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

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

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

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

Граничные и начальные условия

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

Начальные условия будут весьма простыми — изначально жидкость неподвижна (скорость частиц равна нулю), и давление также равно нулю. Граничные условия будут задаваться для скорости и давления приведенными формулами:

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Тем самым, скорость частиц на краях будет противоположна скорости у краев (тем самым они будут отталкиваться от края), а давление равно значению непосредственно рядом с границей. Данные операции следует применить ко всем ограничивающим элементам массива (к примеру, есть размер сетки Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса, то алгоритм мы применим для клеток, отмеченных на рисунке синим):

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Краситель

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

В формуле Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксаотвечает за пополнение красителем области (возможно, в зависимости от того, куда нажмет пользователь), Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксанепосредственно является количество красителя в точке, а Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— коэффициент диффузии. Решить его не составляет большого труда, так как вся основная работа по выводу формул уже проведена, и достаточно лишь сделает несколько подстановок. Краску можно реализовать в коде как цвет в формате RGB, и в таком случае задача сводится к операциям с несколькими вещественными величинами.

Завихренность

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

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

Метод Якоби для решения систем линейных уравнений

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

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

Для расчет вязкости подставляем: Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса, Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса, Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса, здесь параметр Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса— сумма весов. Таким образом, нам необходимо хранить как минимум два векторных поля скоростей, чтобы независимо считать значения одного поля и записывать их в другое. В среднем, для расчета поля скорости методом Якоби необходимо провести 20-50 итераций, что весьма много, если бы мы выполняли вычисления на CPU.

Для уравнения давления мы сделаем следующую подстановку: Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса, Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса, Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса, Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса. В результате мы получим значение Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксав точке. Но так как оно используется только для расчета градиента, вычитаемого из поля скорости, дополнительные преобразования можно не выполнять. Для поля давления лучше всего выполнять 40-80 итераций, потому что при меньших числах расхождение становится заметным.

Видео:Уравнение Навье — Стокса для чайниковСкачать

Уравнение Навье — Стокса для чайников

Реализация алгоритма

Реализовывать алгоритм мы будем на C++, также нам потребуется Cuda Toolkit (как его установить вы можете прочитать на сайте Nvidia), а также SFML. CUDA нам потребуется для распараллеливания алгоритма, а SFML будет использоваться только для создания окна и отображения картинки на экране (В принципе, это вполне можно написать на OpenGL, но разница в производительности будет несущественна, а вот код увеличится еще строк на 200).

Cuda Toolkit

Сначала мы немного поговорим о том, как использовать Cuda Toolkit для распараллеливания задач. Более подробный гайд предоставляется самой Nvidia, поэтому здесь мы ограничимся только самым необходимым. Также предполагается, что вы смогли установить компилятор, и у вас получилось собрать тестовый проект без ошибок.

Чтобы создать функцию, исполняющуюся на GPU, для начала необходимо объявить, сколько ядер мы хотим использовать, и сколько блоков ядер нужно выделить. Для этого Cuda Toolkit предоставляет нам специальную структуру — dim3, по умолчанию устанавливающую все свои значения x, y, z равными 1. Указывая ее как аргумент при вызове функции, мы можем управлять количеством выделяемых ядер. Так как работаем мы с двумерным массивом, то в конструкторе необходимо установить только два поля: x и y:

где size_x и size_y — размер обрабатываемого массива. Сигнатура и вызов функции выглядят следующим образом (тройные угловые скобки обрабатываются компилятором Cuda):

В самой функции можно восстановить индексы двумерного массива через номер блока и номер ядра в этом блоке по следующей формуле:

Следует отметить, что функция, исполняемая на видеокарте, должна быть обязательно помечена тегом __global__ , а также возвращать void, поэтому чаще всего результаты вычислений записываются в передаваемый как аргумент и заранее выделенный в памяти видеокарты массив.

За освобождение и выделение памяти на видеокарте отвечают функции CudaMalloc и CudaFree. Мы можем оперировать указателями на область памяти, которые они возвращают, но получить доступ к данным из основного кода не можем. Самый простой способ вернуть результаты вычислений — воспользоваться cudaMemcpy, схожей со стандартным memcpy, но умеющей копировать данные с видеокарты в основную память и наоборот.

SFML и рендер окна

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

строка в начале функции main

создает изображение формата RGBA в виде одномерного массива с константной длиной. Его мы будем передавать вместе с другими параметрами (позиция мыши, разница между кадрами) в функцию computeField. Последняя, как и несколько других функций, объявлены в kernel.cu и вызывают код, исполняемый на GPU. Документацию по любой из функций вы можете найти на сайте SFML, в коде файла не происходит ничего сверхинтересного, поэтому мы не будем надолго на нем останавливаться.

Вычисления на GPU

Чтобы начать писать код под gpu, для начала создадим файл kernel.cu и определим в нем несколько вспомогательных классов: Color3f, Vec2, Config, SystemConfig:

Атрибут __host__ перед именем метода означает, что код может исполнятся на CPU, __device__ , наоборот, обязует компилятор собирать код под GPU. В коде объявляются примитивы для работы с двухкомпонентными векторами, цветом, конфиги с параметрами, которые можно менять в рантайме, а также несколько статических указателей на массивы, которые мы будем использовать как буферы для вычислений.

cudaInit и cudaExit также определяеются достаточно тривиально:

В функции инициализации мы выделяем память под двумерные массивы, задаем массив цветов, которые мы будем использовать для раскраски жидкости, а также устанавливаем в конфиг значения по умолчанию. В cudaExit мы просто освобождаем все буферы. Как бы это парадоксально не звучало, для хранения двумерных массивов выгоднее всего использовать одномерные, обращение к которым будет осуществляться таким выражением:

Начнем реализацию непосредственного алгоритма с функции перемещения частиц. В advect передаются поля oldField и newField (то поле, откуда берутся данные и то, куда они записываются), размер массива, а также дельта времени и коэффициент плотности (используется для того, чтобы ускорить растворение красителя в жидкости и сделать среду не сильно чувствительной к действиям пользователя). Функция билинейной интерполяции реализована классическим образом через вычисление промежуточных значений:

Функцию диффузии вязкости было решено разделить на несколько частей: из главного кода вызывается computeDiffusion, которая вызывает diffuse и computeColor заранее указанное число раз, а затем меняет местами массив, откуда мы берем данные, и тот, куда мы их записываем. Это самый простой способ реализовать параллельную обработку данных, но мы расходует в два раза больше памяти.

Обе функции вызывают вариации метода Якоби. В теле jacobiColor и jacobiVelocity сразу же идет проверка, что текущие элементы не находятся на границе — в этом случае мы должны установить их в соответствии с формулами, изложенными в разделе Граничные и начальные условия.

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

Расчет завихренности представляет из себя уже более сложный процесс, поэтому его мы реализуем в computeVorticity и applyVorticity, заметим также, что для них необходимо определить два таких векторных оператора, как curl (ротор) и absGradient (градиент абсолютных значений поля). Чтобы задать дополнительные эффекты вихря, мы умножаем Уравнения движения вязкой несжимаемой жидкости уравнения навье стоксакомпоненту вектора градиента на Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса, а затем нормализируем его, разделив на длину (не забыв при этом проверить, что вектор ненулевой):

Следующим этапом алгоритма будет вычисление скалярного поля давления и его проекция на поле скорости. Для этого нам потребуется реализовать 4 функции: divergency, которая будет считать дивергенцию скорости, jacobiPressure, реализующую метод Якоби для давления, и computePressure c computePressureImpl, проводящие итеративные вычисления поля:

Проекция умещается в две небольшие функции — project и вызываемой ей gradient для давления. Это, можно сказать, последний этап нашего алгоритма симуляции:

После проекции мы смело можем перейти к отрисовке изображения в буфер и различным пост-эффектам. В функции paint выполняется копирование цветов из поля частиц в массив RGBA. Также была реализована функция applyBloom, которая подсвечивает жидкость, когда на нее наведен курсор и нажата клавиша мыши. Из опыта, такой прием делает картину более приятной и интересной для глаз пользователя, но он вовсе не обязателен.

В постобработке также можно подсвечивать места, в которых жидкость имеет наибольшую скорость, менять цвет в зависимости от вектора движения, добавлять различные эффекты и прочее, но в нашем случае мы ограничимся своеобразным минимумом, ведь даже с ним изображения получаются весьма завораживающими (особенно в динамике):

И под конец у нас осталась одна главная функция, которую мы вызываем из main.cppcomputeField. Она сцепляет воедино все кусочки алгоритма, вызывая код на видеокарте, а также копирует данные с gpu на cpu. В ней же находится и расчет вектора импульса и выбор цвета красителя, которые мы передаем в applyForce:

Видео:Программа решения уравнений Навье-Стокса. Лекция 1. Введение.Скачать

Программа решения уравнений Навье-Стокса. Лекция 1. Введение.

Заключение

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

Как автор данной статьи, я буду искренне признателен любым комментариям и дополнениям, и постараюсь ответить на все возникшие у вас вопросы под этим постом.

Дополнительный материал

Весь исходный код, приведенный в данной статье, вы можете найти в моем Github-репозитории. Любые предложения по улучшению приветствуются.

Оригинальный материал, послуживший основой для данной статьи, вы можете прочесть на официальном сайте Nvidia (англ). В нем также представлены примеры реализации частей алгоритма на языке шейдеров:
developer.download.nvidia.com/books/HTML/gpugems/gpugems_ch38.html

Доказательство теоремы разложения Гельмгольца и огромное количество дополнительного материала про механику жидкостей можно найти в данной книге (англ, см. раздел 1.2):
Chorin, A.J., and J.E. Marsden. 1993. A Mathematical Introduction to Fluid Mechanics. 3rd ed. Springer.

Канал одного англоязычного ютубера, делающего качественный контент, связанной с математикой, и решением дифференциальных уравнений в частности (англ). Очень наглядные ролики, помогающие понять суть многих вещей в математике и физике:
3Blue1Brown — YouTube
Differential Equations (3Blue1Brown)

Также выражаю благодарность WhiteBlackGoose за помощь в подготовке материала для статьи.

И под конец небольшой бонус — несколько красивых скриншотов, снятых в программе:

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса
Прямой поток (дефолтные настройки)

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса
Водоворот (большой радиус в applyForce)

Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса
Волна (высокая завихренность + диффузия)

Также по многочисленным просьбам добавил видео с работой симуляции:

Видео:Математика в неожиданных местах: Уравнения Навье-Стокса - 3Скачать

Математика в неожиданных местах: Уравнения Навье-Стокса - 3
  • В книжной версии

    Том 21. Москва, 2012, стр. 650

    Скопировать библиографическую ссылку:

    • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса
    • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса
    • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса
    • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса
    • Уравнения движения вязкой несжимаемой жидкости уравнения навье стокса

    НАВЬЕ́ – СТ О́КСА УРАВНЕ́НИЯ, диф­фе­рен­ци­аль­ные урав­не­ния дви­же­ния сплош­ной сре­ды (жид­ко­сти или га­за), учи­ты­ваю­щие её вяз­кость. Вы­ве­де­ны Л. На­вье в 1822 (опубл. в 1827) на ос­но­ве уп­ро­щён­ной мо­де­ли мо­ле­ку­ляр­ных взаи­мо­дей­ствий. В 1845 Дж. Стокс в ре­зуль­та­те изу­че­ния ста­цио­нар­но­го дви­же­ния не­сжи­мае­мой жид­ко­сти по­лу­чил эти урав­не­ния в совр. фор­ме с ис­поль­зо­ва­ни­ем за­ко­нов со­хра­не­ния мас­сы и им­пуль­са для сплош­ной сре­ды.

    🎥 Видео

    Вязкость. Ламинарное и турбулентное течения жидкостей. 10 класс.Скачать

    Вязкость. Ламинарное и турбулентное течения жидкостей. 10 класс.

    Вывод уравнений движения идеальной жидкости - Лекция 2Скачать

    Вывод уравнений движения идеальной жидкости - Лекция 2

    Уравнение навье стоксаСкачать

    Уравнение навье стокса

    Гладкое решение уравнения Навье — СтоксаСкачать

    Гладкое решение уравнения Навье — Стокса

    Р.В. Шамин. Лекция № 2 Уравнения гидродинамики идеальной жидкости со свободной поверхностьюСкачать

    Р.В. Шамин. Лекция № 2 Уравнения гидродинамики идеальной жидкости со свободной поверхностью

    Visual C++. Программа решения уравнений Навье-Стокса. Урок 6. Решение модельного уравнения энергииСкачать

    Visual C++. Программа решения уравнений Навье-Стокса. Урок 6. Решение модельного уравнения энергии

    лекция 5 гидромеханика веснаСкачать

    лекция 5 гидромеханика весна

    Лекция 2. Уравнение неразрывностиСкачать

    Лекция 2.  Уравнение неразрывности

    Механика - Движение вязкой жидкости. Эффект МагнусаСкачать

    Механика - Движение вязкой жидкости. Эффект Магнуса

    Вывод уравнения неразрывности - Лекция 1Скачать

    Вывод уравнения неразрывности - Лекция 1

    Горбачев В. И. - Основы механики сплошных сред. Часть 2 - Теория течения вязкой жидкости. Часть 1Скачать

    Горбачев В. И. - Основы механики сплошных сред. Часть 2 - Теория течения вязкой жидкости. Часть 1

    Математика в неожиданных местах: Уравнения Навье-Стокса - 1Скачать

    Математика в неожиданных местах: Уравнения Навье-Стокса - 1

    Урок 137. Движение тела в жидкости и газе.Скачать

    Урок 137. Движение тела в жидкости и газе.

    Халилов В. Р. - Теоретическая механика - Гидродинамика вязкой среды (Лекция 15)Скачать

    Халилов В. Р.  -  Теоретическая механика  -  Гидродинамика вязкой среды (Лекция 15)
    Поделиться или сохранить к себе: