В этой статье в качестве эксперимента я постараюсь максимально доступно рассказать, как работает новый метод расчёта гидродинамики, основанный на решении уравнения Шрёдингера.
Всем привет. В этой статье я хотел бы рассказать о новом методе расчёта гидродинамики, основанном на решении уравнения Шрёдингера вместо уравнений, типично используемых для гидродинамики вроде Навье-Стокса. Сам метод очень подробно и полно раскрыт в диссертации Albert Chern'а, названной "Fluid Dynamics with Incompressible Schrödinger Flow". Однако, статья Chern'а кому-то может показаться написанной на не самом доступном языке, поэтому своей статьёй я бы хотел в первую очередь если не объяснить в деталях, как работает этот метод, то хотя бы объяснить, какими интересными свойствами он обладает, и что же именно скрывается за его математикой. Попутно я кратко расскажу о том, как устроены классические методы расчёта гидродинаимики и как новый подход от них отличается. В качестве эксперимента я бы хотел попробовать написать статью так, чтобы каждый, кто отдалённо интересуется программированием физики, нашёл в ней что-то интересное, понятное, и новое для себя — от начинающего программиста до бывалых расчётчиков.
Уравнение Навье-Стокса — это уравнение, описывающее течение жидкости. Уравнение Шрёдингера — это уравнение, описывающее эволюцию квантовой системы. И вдруг выясняется, что эти уравнения на самом деле при некоторых условиях являются тождественными. Не знаю, как у вас, но у меня после прочтения этого факта возник только один вопрос: "o___O?". Я достаточно давно занимаюсь гидродинамикой в качестве хобби и работы, в институте с переменным успехом проходил квантмех и я без понятия, но я понятия не имел, что эти уравнения (а значит и соответствующие им физические явления) настолько близки. Сначала я подумал, что парень, который этот факт открыл, уже готовит шампанское к Нобелевке, но нет — как выяснилось, факт эквивалентности этих уравнений был открыт ещё в начале века неким Madelung'ом, который пытался найти физическую аналогию для объяснения уравнения Шрёдингера.
Почему это важно? В первую очередь потому, что это обозначает глубинное родство квантовомеханических и гидродинамических систем. В диссертации того паренька больше сотни страниц уделено тому, как это вообще так получилось. С участием явления сверхтекучести, которая является загадочным связующим звеном, так как проявляет очевидные свойства идеальной жидкости, являющиеся исключительно следствием квантовой механики. Я же в этой статье далее я рассмотрю только некоторые из параллелей, которые из этого, простите, вытекают.
Следующее очень важное следствие эквивалентности уравнения Шрёдингера и Навье-Стокса — это что решение одного из них эквивалентно решению другого. Так вот уравнение Навье-Стокса — нелиненое, его очень неудобно и неэффективно в общем случае решать, в то время как уравнение Шрёдингера — линеное и его решать гораздо проще. Чтобы составить представление, насколько же неудобным по сей день считается уравнение Навье-Стокса, могу сообщить, что существует целый международный фонд грантов для исследователей, которым хоть какую-то базу под них подстроит, так как(цитата):
Even basic properties of the solutions to Navier–Stokes have never been proven.
Уравнение Шрёдингера же, хоть и описывает мутную квантовую физику, поддаётся решению гораздо легче и эффективнее. Короче, я могу очень долго гудеть про то, как это невероятно и офигенно, но давайте уже перейдём к чему-то более конкретному.
Решение классической гидродинамики на пальцах
Всю эту страницу я буду рассказывать, как работает стандартный солвер Навье-Стокса, который изобрели давным-давно. В основном — для сравнения, для некоторого бэграунда и, возможно, кому-то будет просто интересно почитать. Если вы открыли статью только из интереса именно к уравнению Шрёдингера в гидродинамике, то листайте сразу на вторую страницу.
Что вообще такое — уравнение гидродинамики? Что такое уравнение Навье-Стокса и как его понять? С ответом на этот вопрос гораздо лучше меня справились миллионы авторов статей по этому делу, например, классическая статья от нвидии, по которой многие начинали: https://developer.download.nvidia.com/books/HTML/gpugems/gpugems_ch38.html Однако, я попробую написать очень сжато и на пальцах, что это всё значит и что с этим обычно делают.
Уравнение Навье-Стокса описывает закон, которому обязана подчиняться скорость каждой точки пространства, заполненного равномерной несжимаемой жидкостью. Представьте себе, например, бассейн с водой, в котором выделили некоторый куб, достаточно далеко от стенок, поверхности и дна, в котором нет ничего кроме воды. Вода в нём может как угодно течь, но не может ни образовывать пузырей, ни с чем-то сталкиваться (мы для простоты опустим эти эффекты). Тогда само уравнение Навье-Стокса описывает закон, которому будет подчиняться скорость каждой точки воды в этом кубе:
+ Показать
− Скрыть
\(\frac{\partial \vec u}{\partial t}=-(\vec u \cdot \vec \nabla)\vec u-\frac{1}{\rho}\vec \nabla p + \nu {\vec \nabla}^2 \vec u + \vec F\) \(\vec \nabla \vec u = 0\)
прежде чем вообще смотреть на это уравнение, предлагаю сразу из него выбросить ненужное — то, что нам всё равно не пригодится для понимания и только место занимает. Это член, отвечающий за диффузию \(\nu {\vec \nabla}^2 \vec u\) (у идеальной жидкости один фиг диффузии нет), и за внешнюю силу \(\vec F\) (так как мы обойдёмся без неё). Остаётся система:
\(\frac{\partial \vec u}{\partial t}=-(\vec u \cdot \vec \nabla)\vec u-\frac{1}{\rho}\vec \nabla p\) \(\vec \nabla \vec u = 0\)
Здесь перевёрнутый треугольник называется оператором Набла, который обозначает дифференцирование. Причём смысл этого оператора меняется в зависимости от того, где именно он стоит (например, перед вектором или скаляром). Я постараюсь объяснить смысл каждого его вхождения по порядку. На пальцах смысл всей формулы в следующем. \(\vec u(\vec x)\) — это значение скорости жидкости, которое определяется в каждой точке пространства \(\vec x\). Уравнение описывает закономерности, которым обязана подчиняться эта величина, если она описывает поведение несжимаемой жидкости. Работает хоть для двумерного, хоть для трёхмерного случая. В левой части первого уравнения стоит \(\frac{\partial \vec u(\vec x)}{\partial t}\) — это величина называется производной по времени и показывает, как быстро и куда(это вектор) изменится скорость в точке \(\vec x\) в момент времени \(t\).
+ Как представить производную по времени
− Скрыть
Нулевой вектор производной по времени обозначает, что скорость в этой точке сейчас не меняется, а, например, вектор (10, 0)[м/c2] обозначает, что за следующую секунду скорость вырастет на 10[м/с] по оси x(если сама производная не поменяется).
Слагаемое вида \(-(\vec v \cdot \vec \nabla)\vec u\) называется адвекцией и говорит, что поле скоростей \(\vec u\) в этой точке утекает в направлении \(\vec v\). В нашем же случае \(\vec u = \vec v\), то есть поле скоростей сносит само себя. Это, кстати, и называется нелинейностью и из-за этого возникает миллион проблем при решении этого уравнения.
+ Как представить производную по направлению
− Скрыть
В принципе, смысл этого члена достаточно интуитивно можно представить именно как утекание каждой точки воды по вектору её скорости. Однако, в общем случае производная векторного поля \(\vec u\) по направлению \(\vec v\) обозначается как \((\vec v \cdot \vec \nabla)\vec u\) и обозначает, как меняется функция \(\vec u\) в направлении \(\vec v\) для этой точки.
Слагаемое же \(-\frac{1}{\rho}\vec \nabla p\) является ускорением, которое получает жидкость в точке из-за градиента давления.
+ Как представить градиент
− Скрыть
Оператор \(\vec \nabla\), действующий на скалярное поле(например, давление), называется градиентом. Если слева от некоторой точки давление больше, чем справа, то градиент в ней будет направлен вправо и будет увлекать за собой жидкость в этом направлении. Например, ветер всегда дует в направлении, обратном градиенту давления воздуха (отсюда и минус). Электрический ток течёт в направлении градиента электрического потенциала:
\(E=\vec \nabla \phi\)
Второе уравнение \(\vec \nabla \vec u = 0\) называется уравнение непрерывности, а оператор \(\vec \nabla\) здесь действует на вектор и называется дивергенцией.
+ Как представить дивергенцию
− Скрыть
Оператор дифференцирования, действующий на вектор, называется дивергенцией. Дивергенция, равная нулю, говорит, что для каждого маленького кубика сколько в него жидкости втекает, столько и вытекает. А так как любой объём можно разбить на маленькие кубики, то свойство будет справедливо и для объёма любой формы. Это свойство называют также условием несжимаемости, так как если бы в какой-то объём втекало больше жидкости, чем вытекало, это бы означало, что жидкость в объёме накапливается, сжимаясь. Другой случай применения дивергенции, который может помочь её представить — это теорема Гаусса:
\(\vec \nabla E=\rho\)
Эта теорема говорит, что напряжённость электрического поля, которая "вытекает" из некоторого объёма, всегда вызвана электрическим зарядом плотности \(\rho\) внутри этого объёма. Если в объёме заряда нет, то и дивергенция нулевая.
То есть, одним предложением уравнение Навье-Стокса можно описать так: темп изменения скорости определяется течением и градиентом давления, но жидкость при этом не может сжиматься.
Классическое решения уравнения Навье-Стокса
Как я уже говорил раньше, аналитически это уравнение решить нельзя. По крайней мере, науке способы его решения на бумаге с ручкой неизвестны. Но можно решить численно, то есть с помощью программы с некоторой погрешностью. Для начала определимся, над какими вообще данными будет работать наша программа и в каком представлении. Так как уравнение описывает закон, которому подчиняется скорость жидкости, то проще всего так и поступить — представить скорость каждой точки жидкости в некотором условном кубе как трёхмерный массив векторов. Для двумерного случая — то же самое, только массив двумерный. Важно понимать, что сетка остаётся на месте, а жидкость течёт вдоль скорости, которая хранится в этой сетке. Такой подход называется подходом Эйлера (альтернативный подход, в котором сама сетка съезжает по течению, называется Лагранжевым). Так как жидкость несжимаемая, то плотность её постоянна, для простоты примем, что она равна единице.
Посмотрим теперь, как это уравнение можно программно решить. Для этого можно использовать подход, который называется расщеплением — разбить сложный физический процесс, состоящий из нескольких элементарных, на отдельные чередующиеся стадии и считать, что на каждой стадии работает только один элементарный процесс, а остальные выключены. Как ни странно, можно доказать (см. статью выше), что это — на самом деле математически обоснованная стратегия. Поэтому будем считать, что состояние скоростей для каждой точки в текущий момент времени \(\vec u(\vec x, t)\) нам известно. А для расчёта состояния в следующий момент времени \(t+dt\), разобьём сложный процесс гидродинамической эволюции на простые стадии:
1) снесём поле скоростей по течению. это может немного "сжать" жидкость.
2) найдём такое давление, чтобы жидкость "расжалась".
Первый шаг называется адвекцией, второй — проекцией.
Адвекция
Адвекция, или течение, можно приближённо посчитать достаточно легко — если известно, что в точке \(\vec x\), в момент времени \(t\) скорость равна \(\vec u(\vec x, t)\), то в момент времени \(t+dt\) скорость в неё притечёт жидкость из точки \(\vec x - \vec u(x, t)\cdot dt\).
\(\vec u^*(\vec x, t+dt)=\vec u(\vec x - \vec u(\vec x, t), t)\)
То есть мы получили промежуточное значение скорости, котороже уже утекло по течению, но теперь в нём нарушено условие непрерывности.
Это особенно удобно программируется на GPU, так как это можно посчитать, если хранить скорость в текстуре и её обновлять, просто читая тексели со смещением \(- \vec u(x, t)\cdot dt\) и используя стандартную аппаратную линейную интерполяцию.
Проекция
Проекция берёт скорость, для которой нарушено условие непрерывности \(\vec u^*\) и ищет такое давление, которое её "выправит" до нормальной скорости \(\vec u\). Умные мужики доказали, что такое поле можно найти единственным образом и оно всегда будет градиентом некоторого скалярного поля (давления, в нашем случае):
\(\vec u(\vec x, t+dt)=\vec u^*(\vec x, t+dt) + \vec \nabla p\)
Помножим обе стороны этого равенства на оператор дифференцирования:
\(\vec \nabla \vec u(\vec x, t+dt)=\vec \nabla \vec u^*(\vec x, t+dt) + \vec \nabla^2 p\)
"ПОГОДИ-КА СУСЕЛ, ЭТО ЕЩЁ ЧТО" — можете меня спросить вы. Всё по порядку, но на самом деле отсюда для общего понимания достаточно знать, что если \(\vec u^*(\vec x)\) известно(а оно известно), то отсюда можно найти давление \(p(\vec x)\). Если вспомнить, что в нашем случае дивергенция скорости равна нулю, то остаётся вот такое выражение.
\(\vec \nabla^2 p=-\vec \nabla u^*\)
В правой части этого равенства стоит дивергенция скорости, которую можно легко приблизительно посчитать, если известна скорость \(\vec u^*\)(а она известна). В левой части стоит штука, которая называется лапласианом давления.
+ Как представить лапласиан
− Скрыть
Лапласиан — это оператор дифференцирования (ещё называется оператор набла) в квадрате, то есть применённый дважды к скалярному полю. Первый раз применяем оператор дифференцирования — получаем градиент. Второй раз — получаем дивергенцию. Таким образом оператор лапласа — это дивергенция градиента скалярного поля. Его можно представить как изменение потока скорости через маленький кубик, которое будет вызвано давлением в точке. Ещё одна аналогия — как поменяется дивергенция электрического поля в объёмчике, если в него положить заряд плотностью \(\rho\) (опять же, теорема Гаусса):
\(\vec nabla \vec E = \rho\), \(\vec nabla \phi=\vec E\) => \(\vec nabla^2 \vec phi = \rho\)
Уравнение вида "лапласиан чего-то неизвестного равен чему-то известному" называется уравнением Пуассона. Что бы это ни значило, существует стандартный итеративный алгоритм, который позволяет его решить, то есть найти такое давление, чтобы его лапласиан был равен чему угодно. "Что угодно" мы знаем — это дивергенция промежуточной скорости, поэтому считаем по ней давление. Далее для давления считаем градиент и вычитаем результат из промежуточной скорости, чтобы получить окончательную скорость для следующего шага по времени:
\(\vec u=\vec u^* + \vec \nabla p\)
Шаги адвекции и проекции повторяем до посинения, рассчитывая всё дальше и дальше эволюцию поля течений по времени. Для визуализации можно, например, напускать частиц, которые могу сноситься этим полем скоростей. Результат выглядит так:
Важно понять, что в этом видосе, равно как и во всех остальных гифках этой статьи, жидкость на самом деле находится в большом кубе (границы которого не показаны), а не только там, где видны частицы. Частицы только уносятся полем скоростей, как, например, частицы дыма уносятся полем скоростей воздухе. Сами частицы никакой роли в физике процесса не играют и только позволяют относительно наглядно его продемонстрировать. Частицы обычно добавляются заранее туда, где ожидаются какие-то интересные турбулентности.
Важные особенности классического подхода
"Всё здорово, сусел, но в названии статьи ты написал что-то там про Шрёдингера! Он вообще где? Зачем нам это всё?" — спросите вы. Вопрос резонный. Но всю крутость подхода со Шрёдингером можно осознать, только если иметь представление о слабых сторонах классического солвера, который мы рассмотрели в предыдущей главе. В чём же они заключаются? Давайте об этом поговорим.
Основа любого расчётного метода — это то, как в нём представлены моделируемые данные. В рассмотренном нами подходе мы храним значение скорости для каждой точки. Например, в текселях двумерной или трёхмерной текстуры. Этот способ здорово работает, если требуется описать ровное поле течений, в котором нет особенностей (так называются завихрения и разные другие неоднородности). Неоднородностей обычно нет в вязких жидкостях вроде мёда или майонеза, поэтому метод очень здорово подходит, чтобы моделировать майонез. Но более текучие среды (например, вода, воздух и дым) отличаются тем, что в них существенную роль играют злополучные турбулентные течения — мелкие завихрения, имеющие очень сложную и нерегулярную структуру, даже образующие фракталы, которые очень неудобно описывать просто их значениями в каждой точке текстуры/массива. Если попытаться их моделировать, то все мелкие особенности быстро смазываются и расплываются, что соответствует поведению вязкой жидкости. Такое поведение называется численной вязкостью — это вязкость жидкости, которая появляется не потому что она является частью уравнения, которое мы решаем, а это паразитная вязкость, всплывающая как паразитное следствие нашего метода решения. Более того, напомню, что первое, что мы сделали, не успев взглянуть на уравнение Навье-Стокса — выкинули из него вязкость, так в ней недостатка точно не будет.
А вот избавиться от вязкости гораздо труднее, чем случайно её посчитать. Один из способов — это измельчать расчётную сетку. Чтобы таким методом получить что-то хоть как-то похожее на дым, понадобится сетка минимум 1024x1024x1024, то есть как минимум гигабайт памяти, если хранить по 1 байту на узел. А хранить захочется как минимум трёхкомпонентную скорость, то есть, скорее всего, 32 гигабайта в сумме. Это не только не разумно с точки зрения затрат памяти, это ещё и очень медленно. Другой способ — это представлять скорость не её направлением в каждой точке, а как сумму маленьких элементарных вихрей. Этот метод называется также методом дискретных вихрей. В нём вообще всё не так просто с процессами порождения новых вихрей и удаления старых, с поддержанием нужной плотности (так как вихри друг друга уносят, как частицы) и ещё миллион проблем, можете сами почитать, если интересно. Другой подход основан на том, что в реальных течениях вихри имеют свойство образовывать вращающиеся нити. Представьте медленно движущийся жгут, вокруг которого быстро вращается жидкость. Если такой жгут замыкается в кольцо, получается тороидальный вихрь, образующий знакомое кольцо дыма:
Существуют подходы, которые вместо хранения величины скорости в точках, хранят именно параметры таких жгутов. Но такие методы полагаются на топологию, поэтому в них необходимо считать, как жгуты взаимодействуют, сливаются, распадаются и вообще происходящее быстро теряет простоту и наглядность.
Однако, у классического метода есть одно очень важное положительное свойство — в нём вообще нет параметров. Обратите внимание, что для расчёта используется только скорость и больше вообще ничего — ни вязкости, ни даже плотности. В уравнении Навье-Стокса без вязкости есть плотность, но её можно "спрятать" в нормировку давления, поэтому можно сказать, что в исходном уравнении параметров также нет. Забегая вперёд, замему, что в солвере на уравнении Шрёдинге будет параметр. Загадочный.
На следующей странице мы рассмотрим, как же применить уравнения Шрёдингера, чтобы смоделировать тот же самый процесс, и какой в этом профит. Будет много картинок.
