Численные методы решения уравнения больцмана

Моделирование гидродинамики: Lattice Boltzmann Method

Моделирование извержения вулкана
с помощью Lattice Boltzmann Method. (с) Источник

В этой статье я расскажу о численном методе моделирования гидродинамики Lattice Boltzmann Method, LBM. На русском—метод решёточных уравнений Больцмана. Он превосходит другие известные методы (например, finite element method) в легкости распараллеливания, возможности моделирования многофазных потоков, моделировании потоков в пористых средах. Кроме того, вычислительный алгоритм содержит только простейшие арифметические операции. Метод весьма новый, первые коммерческие продукты на его основе стали появляться около 2010 года.

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

Зачем все это нужно? Точнее, в каких отраслях нужно моделировать гидродинамику:

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

Статья будет включать следующие разделы:

• Обзор физики — высокоуровневый обзор основных необходимых уравнений
• Основная идея—описание основного принципа алгоритма
• Технические детали—более подробное объяснение исходных уравнений, вывод вычислительной схемы
  • Уравнение Навье-Стокса
  • Уравнение Больцмана
  • Дискретное уравнение Больцмана
  • Иллюстрация вычислительной схемы
  • Пространственные решетки
  • Равновесные функции распределения
  • Несжимаемость
  • Вязкость и число Рейнольдса
  • Еще раз, все вместе

• Разное
  • Дополнения к алгоритму
  • Подводные камни
  • Существующие решения
  • Что почитать

Обзор физики

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

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

Основная идея

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

Ситуация иллюстрируется картинкой ниже.

Вопросительный знак в ней символизирует тот факт, что я не знаю, какое уравненине описывает поведение плотных жидкостей на микроскопическом уровне. Иногда вместо «микроскопический» говорят «мезоскопический»—в том смысле, что микроскопическое описание—это описание поведения отдельных атомов и молекул, а уравнение Больцмана описывает потоки молекул.

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

Технические детали

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

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

Вариант уравнения для макроскопической динамики несжимаемых жидкостей и газов выглядит так:
(1)

здесь v – скорость потока, ρ – плотность жидкости, p – давление в жидкости, f – внешние силы (например, гравитация).

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

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

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

Если теперь подставить эти величины в уравнение Ньютона F = ma, после несложных преобразований мы получим уравнение выше. Слева находится ma, справа—F.

Уравнение Больцмана

Это уравнение оперирует функцией распределения плотности вероятности частиц по координатам и по скоростям, f(r, v,t). Величина f(x, y, z, v_x, v_y, v_z, t) dx dy dz dv_x dv_y dv_z показывает, какая доля частиц в момент времени t находится в кубе от x до x+dx, от y до y + dy, от z до z + dz со скоростями в диапазоне от v_x до v_x + dv_x, от v_y до v_y + dv_y, от v_z до v_z + dv_z. Ее также можно записать в виде .

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

Основная идея для вывода уравнения похожа на вывод уравнения Навье-Стокса. Давайте мысленно выделим в данный момент времени в данном небольшом объеме пучок молекул, летящих в данном направлении (точнее, узком конусе направлений). Через небольшой промежуток времени dt они окажутся в соседней точке (благодаря наличию скорости), их скорость сама по себе изменится из-за ускорения молекул внешними силами. Кроме того, на этом отрезке пути некоторые молекулы столкнутся с другими, изменят свою скорость, и мы больше не сможем включать их в исходный пучок. С другой стороны, в результате столкновений молекул в том же объеме, летящих в другую сторону, некоторые из них приобретут нужное направление скорости, и мы добавим их в пучок.

Это можно записать в следующем виде:
, (4)
где F—внешняя сила, m—масса молекулы, dN_coll—изменение числа частиц в пучке за счет столкновений.

Как определяется в общем виде величина dN_coll, от читателя останется скрытым. Нам понадобится только ее стандартное приближение—Батнагара-Гросса-Крука (BGK). В этом приближении dN_coll равно
, (5)
где f eq —равновесная функция распределения, распределение Максвелла-Больцмана, а τ – так называемое время релаксации.

В итоге получаем
. (6)
Заметим, что f eq зависит от макроскопических плотности и скорости в данной точке (то есть неявно от координаты и времени). Именно это уравнение нам понадобится в дальнейшем, но обычно его делят на dt и приводят к виду
, (7)
где набла с индексом v – это набла по переменным скорости.

Дискретное уравнение Больцмана

Чтоб получить возможность моделировать динамику непрерывного уравнения Больцмана на компьютере, нам необходимо его дискретизировать. Для этого сначала введем равномерную сетку пространственных координатам—шаг сетки пусть будет одинаковым по всем осям. Поведение жидкости мы будем определять именно в узлах сетки. Фактически, мы разрешаем молекулам находиться только в определенных пространственных узлах. Кроме того, мы должны дискретизировать время—мы будем определять состояние жидкости в равноотстоящие друг от друга моменты времени. Кроме того, мы позволим молекулам иметь только определенные значения скорости—такие, чтоб за шаг по времени они успевали перейти в какой-нибудь соседний узел. Эти разрешенные направления будут одинаковыми для всех пространственных узлов. Очевидно, что по диагональным направлениям скорости частиц будут больше, чем по недиагональным.

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

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

Для простоты ниже будем предполагать, что внешние силы отсутствуют. Мы пронумеруем разрешенные направления скорости от 1 до Q с индексом i. Если теперь массу частиц, пролетающих из данного узла в направлении i за шаг времени обозначить как f_i, то уравнение (6) перейдет в
. (8)
Здесь мы учли, что шаг по времени равен единице, и заменили все dt из (6) на единицу. f_i eq обозначает дискретную равновесную плотность распределения, зависящую от макроскопических массы и скорости в данном узле. Мы не указали, из какого именно узла будем использовать f_i eq —из r + v_i t в момент времени t + 1 или из r в момент времени t. Для вычислительной схемы удобнее оказывается использовать узел r + v_i t в момент времени t + 1. Тогда уравнение выше можно разложить на две составляющие, распространение (streaming step) и столкновение (collision step).

Streaming step:
. (9)
Collision step:
. (10)

Здесь f_i с тильдой обозначает массу частиц, пришедших в узел по направлению i, но еще не столкнувшихся с остальными пришедшими частицами. Streaming step иногда называют advection step.

Для дискретизированных направлений скорости масса и макроскопическая скорость в каждом узле будут рассчитываться как
. (11)

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

Мы покажем, что равновесные функции распределения зависят от массы в узлах и макроскопической скорости. Поэтому после streaming step необходимо пересчитать массы и скорости в каждом узле, пересчитать равновесные функции распределения, и потом совершать collision.

Таким образом, на каждом шаге вычислительной схемы необходимо «распространить» частицы, то есть переместить частицы, летевшие из узла r в направлении i в узел r + v_i (проделать это для всех частиц и направлений). После этого необходимо пересчитать массы, скорости, равновесные функции распределения. Наконец, необходимо «столкнуть» частицы, прилетевшие в данный узел—то есть перераспределить частицы по направлениям.

Иллюстрация вычислительной схемы

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

На следующем рисунке изображена одна итерация пары «streaming/collision». Цветные стрелки изображают потоки летящих молекул. Интенсивность цвета кодирует массу молекул, летящих в данном потоке, длина стрелок примерно соответствует пути, проходимом потоком за шаг по времени (лишь примерно, поскольку стрелки должны идти от центра узла до центра узла).

Пространственные решетки

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

В LBM любая решетка должна содержать нулевой вектор из узла в себя самого—он описывает частицы, которые никуда не летят из данного узла. В LBM решетки обычно обозначаются аббревиатурой DnQm, где n—размерность пространства, m—число векторов в решетке. Например, D2Q9, D3Q19 и т.п.

В двумерном пространстве LBM решетка может состоять, например, из 5 векторов (2 вертикальных, 2 горизонтальных и нулевой вектор из узла в себя самого), а может из 9 векторов, как на иллюстрации выше (2 вертикальных, 2 горизонтальных, 4 диагональных, один нулевой). Это решетки D2Q5 и D2Q9, соответственно.

Очевидными факторами для выбора решеток являются: 1. точность моделирования (интуитивно понятно, что чем больше векторов в решетке, тем точнее моделирование) 2. вычислительные затраты (расчет на решетке D2Q5 будет быстрее расчета на D2Q9). Как ни странно, это не самые важные факторы. Самые важные факторы—это воспроизводимость уравнений Навье-Стокса и симметрия некоторых тензоров, построенных на векторах решетки.

Обычно используются решетки D2Q9, D3Q15, D3Q19. Решетки D2Q9 и D3Q19 изображены ниже. Базисные вектора решетки обычно обозначаются как e_i или c_i (они совпадут со скоростями v_i, введенными ранее при единичном шаге по времени). Ниже мы будем использовать обозначение e_i.

Выпишем базисные вектора для D2Q9:
(11)

и для D3Q19:
(12)

Еще раз напомню, что мы везде предполагаем, что шаг по времени равен единице, поэтому v_i = e_i.

Равновесные функции распределения

В непрерывном случае равновесная функция распределения (распределение Максвелла-Больцмана) равна
. (13)

Здесь есть ранее не встречавшиеся величины: R—универсальная газовая постоянная, T—температура, D—размерность пространства, v—вектор скорости, плотность вероятности для которого мы и находим. Здесь молярная масса газа принимается равной единице (она для нас не важна—важна только макроскопическая плотность). Можно считать это изменением единицы массы в нашем моделировании. Кроме того, функция нормируется на локальную макроскопическую плотность, а не на единицу. Заметим также, что обычно от v не отнимают макроскопическую скорость газа u. Это происходит потому, что обычно распределение исследуют в случае неподвижного газа, нам же необходимо перейти в локальную инерционную систему отсчета, движущуюся с текущей скоростью газа в данной точке в данный момент времени, чтоб использовать распределение Максвелла-Больцмана. Вычитание u и является таким переходом.

Предположим, что в данной точке пространства распределение молекул по скорости оказалось равновесным. Это распределение зависит от макроскопической массы ρ и скорости u. С другой стороны, мы можем вычислить ρ и u по функции распределения (см. формулы (2) и (3)). Очевидно, это вычисление должно дать правильные ρ и u (то есть это в каком-то смысле дополнительное ограничение на функцию распределения):
, . (14)

Такое же требование мы налагаем на вычисление плотности и скорости в дискретном случае:
. (15)

Основным требованием для дискретных равновесных функций распределения является воспроизводство уравнения Навье-Стокса в пределе бесконечно малых шага по времени и шага решетки. Это эквивалентно тому, что если при заданных ρ и u мы попытаемся их снова рассчитать по равновесным функциям распределения в непрерывном и дискретном случае, результаты совпадут (см. замечание выше о том, что в дискретном случае вместо плотности имеется в виду масса), т.е.
. (16)

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

Если просуммировать по направлениям скорости уравнение (10) для collision step, то с учетом формулы (16) можно показать, что collision step не изменяет макроскопические массы и скорости молекул, находящихся в узле.

Оказывается, если просто подставить дискретные векторы скорости e_i=v_i из (11) и (12) в непрерывную равновесную функцию распределения (13), равенства (16) не будут выполняться.

Оказывается также, что мы можем заставить равенства (16) выполняться, если разложим распределение Максвелла-Больцмана (13) в ряд Тейлора до второго порядка макроскопической скорости. Это соответствует тому, что величина u/sqrt(RT) очень мала. Это ограничение всегда должно выполняться в процессе моделирования во всех узлах.

Но даже разложения в ряд Тейлора недостаточно. Мы также должны ввести в дискретизированные функции f_i eq специально подобранные множители w_i (детали вычислений можно найти в этой прекрасной статье—есть и бесплатная версия; конечно же, все окажется немного сложнее, чем в нашем поверхностном описании—на самом деле, расчет происходит через тензоры, вплоть до четвертого ранга, построенные на векторах решетки). Окончательно получим
. (17)

Здесь кроется подвох: мы везде предполагаем, что шаг решетки и времени являются единицами длины и времени, соответственно. Потому мы не можем взять значение R из СИ, а температура здесь не равна температуре моделируемой жидкости в Кельвинах.

Чтоб определить их значения, заметим следующее. Предположим, что в жидкости имеется возмущение—то есть в некоторых узлах есть избыточная масса. Эта масса начнет «растекаться» по пространству, притом в крайнем правом узле в области возмущения она будет двигаться в том числе и по направлению (1, 0, 0) в 3D или (1, 0) в 2D. За единицу времени молекулы вдоль этих направлений пройдут единицу длины, то есть их скорость будет равна единице. Это означает, что скорость звука как скорость распространения возмущений в нашей системе тоже равна единице. С другой стороны, скорость звука равна sqrt(γ R T / μ), где γ—это постоянная адиабаты, μ—молярная масса, которую мы приняли равной единице ранее. Постоянная адиабаты γ равна 1 + 2 / d, где d—число степеней свободы молекул. В идеальном газе оно равно размерности пространства. В нашем газе молекулы могут двигаться только вдоль прямых, соединяющих узлы, потому размерность равна не трем (или двум), а единице. То есть γ=3, а sqrt(3 R T) = 1.

Обычно в литературе по LBM под «скоростью звука» понимают величину
. (18)

Теперь, окончательно,
. (19)

Выпишем значения коэффициентов w_i для самых распространенных решеток.
Для D2Q9:
(20)

Для D3Q19:
(21)

Несжимаемость

Ограничение на малую макроскопическую скорость жидкости теперь можно записать как
. (22)

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

Хорошее максимальное значение для макроскопической скорости в узлах будет, например, 0.01.

Вязкость и число Рейнольдса

Кинематическая вязкость ν в LBM (как обычно, в единицах решетки) рассчитывается как
. (23)
Здесь τ—время релаксации, введенное ранее в формуле (5), c_s=1/sqrt(3)— «скорость звука», введенная в (18).

При моделировании гидродинамики без учета изменений температуры любая система с наперед заданной геометрией (например, труба с квадратным сечением) полностью описывается только одним безразмерным параметром—числом Рейнольдса, равным
, (24)
где v—характеристическая скорость в системе (например, скорость в центре трубы), L—характеристическая длина в системе (например, длина стороны квадрата сечения).

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

1. выбрать характеристическую скорость v. Как мы выяснили, она должна быть намного меньше скорости звука. Например, 0.01.
2. рассчитать необходимую внешнюю силу для такой скорости
3. рассчитать вязкость по (23) так, чтоб получить нужное число Рейнольдса
4. рассчитать время релаксации по (24) так, чтоб получить нужную вязкость

Если задача моделирования составлена в единицах СИ (мол, сторона квадрата сечения трубы 1м, давление на входе трубы X Паскалей, на выходе—Y Паскалей), то сначала надо найти безразмерное число Рейнольдса, и воспользоваться алгоритмом выше.

Еще раз, все вместе

Перед моделированием надо задать начальные макроскопические массы и скорости в каждом узле. Потом задать потоки масс в каждом разрешенном направлении e_i в каждом узле (см. Подводные камни). Самый простой способ—указать потоки из равновесного распределения.

Для моделирования в цикле совершать

1. streaming step по формуле (9)
2. пересчет макроскопических масс и скоростей в каждом узле по формулам (11), пересчет равновесных потоков по всем направлениям по формуле (19)
3. collision step по формуле (10)

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

Разное

Дополнения к алгоритму

Мы не коснулись включения внешних сил в моделирование (например, гравитации). Они приведут к небольшой добавке в уравнение (9) для streaming step.

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

Граничные условия на границе твердых тел в гидродинамике обычно выбираются как no-slip boundary conditions (то есть когда скорость жидкости на поверхности тел равна нулю). Обычно они реализуются через bounce-back conditions (когда потоки молекул отражаются в обратном направлении при пересечении твердой стенки).Об этом можно почитать здесь, раздел 4.6.

Модель столкновений, которую мы представили, называется single relaxation time. Существуют более совершенные модели, например, multiple relaxation time (в частности, double relaxation time).

LBM также поддерживает многофазность (моделирование потока смеси жидкостей или газов с разными параметрами).

Наличие теплопроводности (то есть переноса тепла в системе, изменения температуры в разных точках системы и, как следствие, изменения параметров системы—например, плотности) также поддерживается LBM. Температура моделируется как отдельный «газ», тоже по алгоритму LBM, но скорость этого газа определяется скоростью основной жидкости. Говорят, что температура в этом смысле является passive scalar. Есть немало статей по моделированию через LBM явления Rayleigh–Benard convection.

Мы совершенно не коснулись вопросов эффективной реализации и параллелизации.

Подводные камни

Если вы моделируете систему с теплопроводностью, то она будет описываться двумя безразмерными величинами. В случае Rayleigh-Benard convection обычно выбирают число Прандтля и число Рэлея. Для того, чтобы воспроизвести эту систему в единицах решетки, недостаточно просто воспроизвести обе этих безразмерных величины, правильно задав внутренние параметры системы (характеристическая скорость, внешняя сила, теплопроводность). Дело в том, что между внутренними параметрами будут существовать скрытые зависимости. Подробнее можно почитать тут.

Как мы уже сказали, при выборе характеристической скорости в системе не забывайте о том, что число Маха должно быть много меньше единицы (формула (22)).

LBM может стать нестабильным в некоторых системах при высоких числах Рейнольдса (когда поток, тем не менее, все еще ламинарный).

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

Существующие решения

Есть большой и очень зрелый open-source проект, целиком посвященный LBM—Palabos (PArallel LAttice BOLtzmann). У проекта есть и вики. Разработчики предоставляют платную консультацию по моделированию гидродинамики.

Здесь есть прекрасные обучающие примеры моделирования типичных систем на MATLAB. Например, Rayleigh–Benard convection (в наличии внешняя сила, теплопроводность, граничные условия, правильный перевод величин). Всего 160 строчек в MATLAB.

Есть множество коммерческих решений, например, это или это.

Подробный список коммерческих и некоммерческих пакетов по LBM можно найти на википедии.

Что почитать

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

Статистические методы расчета кинетических коэффициентов переноса газа (стр. 5 )

Из за большого объема этот материал размещен на нескольких страницах: 1 2 3 4 5 6 7 8 9

,

.

Результирующий эффект столкновений с частицами всех сортов получается суммированием (3.147, 3.148) по всем сортам:

.

Приравнивая изменение числа выделенных частиц за время dt, подсчитанные через производные от заданной функции fi(r, ui, t) (3.137) и найденное путем подсчета числа столкновений (3.147, 3.148) и сокращая на drduidt, получим искомые уравнения для одночастичной функции распределения, основанные на выражениях (3.148) и (3.147), соответственно:

. (3.149)

. (3.150)

Уравнение (3.149) называется уравнением Больцмана. Оно было составлено в 1872 году и в течение более 40 лет подвергалось критике с механистических позиций. В течение всей последующей жизни Больцману не удалось, а он опубликовал это уравнение 28‑и лет от роду, найти пути его решения. В 60 лет он покончил жизнь самоубийством. Через 12 лет после его смерти были опубликованы работы, в которых удалось на основе этого уравнения вычислить кинетические коэффициенты и предсказать явление термодиффузии.

Для замкнутости постановки задачи отыскания функции распределения к интегро-дифференциальному уравнению (3.149) требуются граничные условия. Они задаются с помощью ядер рассеяния газовых молекул на поверхности R(u¢®u):

(3.151)

Здесь rw — радиус-вектор точек поверхности твердого тела, ограничивающего газ.

3.3.3. Условия применимости уравнения Больцмана

При выводе уравнения Больцмана несколько раз используется условие разреженности, означающее, что длина свободного пробега l должна быть много больше диаметра молекул d (или, точнее, характерного линейного размера сечения взаимодействия молекул). Во-первых, выбор элементов dr, dt подразумевает разреженность, иначе нельзя было бы удовлетворить неравенствам d3 « |dr| « l3, tст « dt « , которые лежат в основе подсчета числа столкновений. Во-вторых, при расчете числа столкновений учитывались только парные столкновения. Вероятность обнаружить частицу в объеме порядка d3 (d — размер молекулы) из предоставленного ей среднего объема 1/n составляет d3n » d/l. Вероятность застать 2 частицы в объеме

d3 (только тогда произойдет их взаимодействие), очевидно, есть произведение вероятностей для отдельных частиц, т. е. , тогда вероятность тройного столкновения

(d3n)3. Таким образом, пренебрегать тройными столкновениями можно только тогда, когда « 1.

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

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

d2/l2, т. е. опять играет роль параметр разреженности.

Таким образом, если разреженность обеспечена, то уравнение Больцмана должно хорошо «работать». Оно может описывать самые разнообразные ситуации с разреженными газами: моноскоростные молекулярные пучки (когда все молекулы летят с одной и той же скоростью u без столкновений друг с другом), любые по величине градиенты температур, воздействие звука любых частот, ударные волны любых амплитуд, разлет молекул в вакууме и т. д.

3.3.4. Число микросостояний N-частичной газовой системы и его стремление к максимальному значению для функций распределения, удовлетворяющих уравнению Больцмана

Пусть заданы функции распределения газовой смеси fj(r, uj, t). Разобьем объем газа V и пространство скоростей на достаточно большое число элементов (Dr, Duj) таких, что в пределах элемента фазового пространства находилось бы достаточно большое количество как частиц, так и состояний одночастичной задачи. С другой стороны, будем считать, что в пределах элемента фазового объема вероятность того, что состояние занято, является одинаковой (такая вероятность будет различна только при заметной неоднородности плотности частиц в пределах фазового элемента). В этом случае удается рассчитать число состояний N-частичной подсистемы сначала в каждом таком элементе, а затем и во всей системе (см. п.2.3.7). Логарифм этого числа состояний есть энтропия, для плотности которой было получено следующее выражение:

. (3.152)

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

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

3.3.4.1. Абсолютный максвеллиан как решение уравнения Больцмана

Для проведения указанных выше доказательств наиболее удобно использовать уравнение Больцмана в обобщенной форме (3.150):

. (3.153)

Подставляя в (3.153) абсолютный максвеллиан для Fi = 0 видим, что левая часть обращается в нуль, так как f00(ui) не зависит от t и r, а в правой части подынтегральное выражение в скобках даст:

,

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

3.3.4.2. Замечательное свойство моментов интеграла столкновений

Для того, чтобы доказать, что энтропия в изолированной системе, описываемой уравнением Больцмана, не убывает во времени, рассмотрим вначале одно интересное свойство моментов интеграла столкновений. Введем обозначение для интеграла ij-столкновений:

. (3.154)

Заметим, что Jij является функцией (r, ui, t). Поэтому, так же как для функций распределения fi(r, ui, t), удобной величиной, сокращающей информацию, является момент интеграла столкновений:

, (3.155)

где ji(ui) — произвольные функции скорости молекул. Это, в частности, могут быть масса молекулы mi, ее импульс miui, энергия и др. Эти функции называют молекулярными признаками.

Выберем два произвольных значения скорости молекул i-го компонента ui и . Без ограничения общности можно сказать, что скорость может быть получена в результате столкновения молекулы со скоростью ui с молекулой j-го компонента.

В дополнение к интегралу столкновений Jij(ui), даваемого выражением (3.154), рассмотрим интеграл столкновений для второго фиксированного значения скорости :

. (3.156)

Запишем теперь выражения для моментов от интегралов (3.154) и (3.156) с молекулярными признаками ji(ui) и :

, (3.157)

. (3.158)

Так как все переменные в интеграле (3.158) «немые» (по всем ведется интегрирование), то можно провести взаимную замену переменных , без изменения величины . Но такая замена делает (3.157) и (3.158) тождественными, и, следовательно,

. (3.159)

Складывая теперь (3.157) и (3.158) и деля на два, получим:

. (3.160)

Здесь учли соотношение взаимности (3.102) и то, что . Повторяя такие же операции для моментов от интегралов столкновений, являющимися функциями скорости молекул j-того компонента Jij(uj) и , нетрудно получить:

. (3.161)

Если просуммировать (3.160) и (3.161) по всем i и j, то в каждую сумму войдут моменты от всех интегралов, входящих в систему уравнений Больцмана с одним и тем же перечнем молекулярных признаков . Поэтому очевидно следующее равенство:

. (3.162)

Производя суммирование по i и j выражений (3.160) и (3.161), а затем складывая их друг с другом и деля на два, получим:

(3.163)

Если мы имеем дело с одним компонентом, то из всей суммы остается только один член ii, но в обозначениях для (3.163) надо будет вводить следующие индексы: i®1, j®2, указывающие на отнесение характеристики к первой или второй молекуле, участвующей в столкновении, т. е.

(3.164)

Если в качестве молекулярного признака j выбираются так называемые сумматорные инварианты столкновений m, mu, mu 2/2, для которых справедливы законы сохранения, то (3.163, 3.164) обращаются в 0. Это легко понять, исходя из физического смысла момента интеграла столкновений.

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

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

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

3.3.4.3. Теорема о неотрицательности производства энтропии за счет столкновений газ-газ

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

. (3.165)

Найдем производную по времени от энтропии смеси (3.152):

. (3.166)

В изолированной системе ¶s/¶t есть не что иное, как плотность производства энтропии. Подставляя в (3.166) вместо ¶fi/¶t его значение из (3.165), будем иметь:

. (3.167)

Здесь в качестве молекулярного признака выступает величина , но в соответствии с (3.163) сумму (3.167) можно выразить через интеграл от разности сумм молекулярных признаков после и до столкновения:

(3.168)

Таким образом, подставляя (3.168) в (3.167), имеем:

(3.169)

Так как величина типа — неотрицательна при любом неотрицательном значении x и y, то, очевидно, значение правой части (3.169) неотрицательно, т. е.

. (3.170)

Это доказательство впервые было проведено Больцманом, и, так как для энтропии, взятой с обратным знаком, он использовал обозначение H, то с тех пор утверждение (3.170) носит название H-теоремы Больцмана.

Легко понять физический смысл выражения для производства энтропии (3.167). Например, для однокомпонентного газа:

. (3.171)

И, поэтому, s есть не что иное, как суммарное по всем столкновениям (в единице объема и за единицу времени) изменение молекулярного признака в результате этих столкновений. Но — это размер квантовомеханической ячейки фазового пространства, — число частиц, приходящееся на эту ячейку. Для разреженного газа эта величина много меньше 1, поэтому можно считать, что — это вероятность того, что элементарная ячейка, относящаяся к заданному объему фазового пространства, является занятой. Эта вероятность определяет число состояний N-частичной системы, а в соответствии с (2.62), (3.152) — это вклад в энтропию, приходящуюся на одну частицу с радиус-вектором r и скоростью u. Следовательно, производство энтропии есть суммарное по всем столкновениям в единице объема за единицу времени изменение вклада в энтропию каждой частицы за счет перехода u ® u¢.

3.3.4.4. Теорема о неотрицательности производства энтропии за счет столкновений газ-поверхность

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

Баланс энтропии вблизи фазовой границы

Пусть s и S есть энтропия и производство энтропии на единицу площади данного слоя. Уравнение баланса энтропии имеет вид:

(3.172)

Здесь jsг, jsт — плотности потоков энтропии в газе и в твердом теле вблизи границы фаз, соответственно.

Для непроницаемого по отношению к газовым частицам твердого тела величина jsт будет обусловлена плотностью потока тепла jqт:

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

Для стационарного случая имеем:

(3.173)

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

(3.174)

(3.175)

и условие непроницаемости

(3.176)

на основе (3.174) можно получить:

(3.177)

Чтобы показать, что разность (3.177) неотрицательна, рассмотрим подынтегральное выражение второго члена (3.177), деленного на f0, с учетом граничного условия (3.124):

(3.178)

Интеграл из (3.178)

можно рассматривать как среднее значение величины f(u¢)/f0(u¢), рассчитанное с весом

который нормируется на 1 в связи со свойством ядра оставлять максвеллиан без изменения (3.128). Поэтому (3.178) можно переписать в следующем виде:

(3.179)

где С — обозначение функции xln(x).

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

(3.180)

если среднее определяется с весом, нормированным на 1.

Функция

Например, если в качестве C(u) взять u2, а нормированную весовую функцию выбрать в виде f0(u)/n, то неравенство (3.180) будет соответствовать тому, что

В этом нетрудно убедиться, так как

Для нашего случая неравенство (3.180) с учетом (3.178, 3.179) дает:

(3.181)

Умножая (3.181) на |un|f0(u)du и интегрируя в пределах un > 0 с учетом нормировки ядра рассеяния R(u¢®u) на единицу, получим:

(3.182)

Сравнивая это выражение с (3.177), нетрудно видеть, что производство энтропии за счет столкновений газовых молекул с поверхностью, как и следовало ожидать, больше или равно нулю (S ³ 0).

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

3.3.5. Анализ уравнения Больцмана методами теории размерностей и подобия

3.3.5.1. Вывод обезразмеренного уравнения Больцмана

Рассмотрим разреженный неравновесный однокомпонентный газ в отсутствие внешних сил. Функция распределения молекул по скоростям для такого газа удовлетворяет уравнению Больцмана:

(3.183)

Здесь u, u¢ и u1, — скорости пары молекул до и после столкновения, соответственно. В качестве характерной длины возьмем длину свободного пробега l, характерной скорости — среднюю тепловую скорость ut, за характерное сечение столкновений примем s0. Тогда можно ввести следующие безразмерные величины, которые в отличие от размерных будем отмечать символом (*):

Используя безразмерные переменные: уравнение (3.183) можно переписать в виде:

Здесь учитываем, что вероятность перехода в интервалы скоростей есть безразмерная величина, нормированная на 1. Так как , то обезразмеренное уравнение имеет такой же вид, как и уравнение (3.183) (отсутствуют параметры подобия):

(3.184)

Уравнение (3.184) показывает, что производная от функции распределения по траектории для фиксированной скорости u* равна разности прихода частиц в число избранных (обладающих скоростью в интервале (u* ¸ u* + du*) и ухода из этого числа за счет столкновений. Приход обусловлен столкновениями частиц со скоростями , u1, такие частицы будем называть полевыми. Проведенное обезразмеривание переменных приводит к тому, что , , , как правило являются величинами порядка 1. Следует отметить, что в сильно неравновесных ситуациях (r*, u*, t*) может быть заключена между 0 и величиной

1 в зависимости от значения аргументов (некоторые области значений аргументов могут быть «пустыми»). Примером сильной неравновесности может служить стационарный молекулярный пучок с фиксированной скоростью u0, поступающий в вакуумный объем.

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

3.3.5.2. Большие отклонения от равновесия

Пусть моноскоростной стационарный (¶f */¶t* = 0) молекулярный пучок плотностью n поступает в неограниченный объем с равновесным газом того же сорта плотностью n0. Если n0 « n, то величина (), связанная с числом столкновений полевых молекул, будет пропорциональна , и этой величиной в (3.184) можно пренебрегать. Функцию f *(u) можно вынести из-под знака интеграла, а интеграл

будет иметь порядок n0/n. Тогда (3.184) дает:

Найдем характерное расстояние Dr, на котором f * изменяется на свою величину (т. е. изменение Df *

f *) при :

Если плотность полевых молекул n0 приближается к плотности частиц в пучке n, то характерное расстояние, на котором значительно изменяется функция распределения, становится величиной порядка длины свободного пробега l. В противном случае |Dr| может неограниченно увеличиваться (|Dr| ® ¥ при n ® 0).

3.3.5.3. Малые отклонения от равновесия

Рассмотрим газ, равновесное состояние которого слабо возмущено. Например, газ заключен в некоторый объем, температура стенок которого является функцией координат, причем . Для простоты опять будем полагать, что нет явной зависимости функции распределения от времени (). При малых отклонениях от равновесия функция распределения представима в виде:

где — локальномаксвелловская функция распределения.

Разность произведений функций распределения в (3.184) в этом случае можно записать в виде:

(3.185)

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

(3.186)

Возвращаясь к размерным переменным, вместо (3.186) будем иметь:

(3.187)

где введено время релаксации функции распределения .

Уравнение (3.187) показывает, что для фиксированной скорости u отклонение функции распределения f(r, u) от максвеллиана f0 пропорционально производной от f(r, u) вдоль траектории движения молекул со скоростью u.

Уравнение типа (3.187) носит название релаксационного уравнения или релаксационной модели уравнения Больцмана. Часто эту модель называют уравнением Бхатнагара-Гросса-Крука (или сокращенно уравнением БГК) в честь авторов, которые впервые в 1954 году применили эту модель.

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

3.3.5.4. Условия, при которых в уравнении Больцмана возникает малый параметр

Рассмотрим неравновесное состояние газа в задаче с характерным размером L. В этом случае величина L позволяет определить производную в (3.186) следующим способом:

где f *(L) и f *(0) значения функции на траектории движения молекул со скоростью u в точках, отстоящих друг от друга на расстояние

Тогда вместо (3.186) имеем:

(3.188)

Таким образом, отклонение функции распределения от максвеллиана пропорционально числу Кнудсена l/L и характерной разности функций распределения (f *(L) – f *(0)). Последняя задается внешними условиями: скоростью обтекания, различием температуры тела и газа и т. д. Даже если максимальная разность (f *(L) – f *(0)) имеет порядок 1 (т. е. наибольшее свое значение), то отклонение функции распределения от максвелловской в каждой точке области течения может быть много меньше единицы благодаря малости числа Кнудсена.

Оценка (3.188) лежит в основе исторически первого метода решения уравнения Больцмана, предложенного в 1916 году Чепменом и, независимо, Энскогом. Метод развит для течений в вязком режиме (Kn « 1) и соответствует постановке задачи в локальной неравновесной термодинамике. Малость числа Кнудсена позволяет область течения размером L разбить на малые части размером L « Dr « l, в каждой из которых газ находится почти в равновесном состоянии с заданной локальномаксвелловской функцией распределения. Основная задача метода — получение из уравнения Больцмана системы уравнений механики сплошной среды с явными выражениями для коэффициентов переноса локальной неравновесной термодинамики D, h, k, характеризующих соответственно диффузию, вязкость и теплопроводность. В соответствии с идеей квазиравновесия в элементе объема граничные условия для функции распределения не принимаются во внимание.

Однако по мере разрежения газа и роста l метод Чепмена-Энскога, основанный на малости функции возмущения Df *, перестает работать, так как в соответствии с (3.188) при l

L величина Df * может стать

1, а столкновения с поверхностью могут приобрести решающую роль.

В разреженной среде, когда L

l, малость отклонения от равновесия (f * – f0*) можно обеспечить только за счет уменьшения максимальной разности функций распределения (f *(L) – f *(0)), (например, на краях области определения). При такой постановке задачи, в частности, можно рассчитать кинетические коэффициенты нелокальной неравновесной термодинамики, например, так были найдены кинетические коэффициенты для отверстия в свободномолекулярном режиме.

Еще один случай возникновения малого параметра в уравнении Больцмана связан с задачами, в которых область определения имеет два различных по порядку величины характерных размера, например, труба длиной L и радиуса R, причем R « L. Длина свободного пробега молекул газа в таких системах l по мере разрежения стремится к величине 2R. Тогда в соответствии с (3.188):

Таким образом, в этом случае снова имеем возможность существования малых возмущений Df *, несмотря на большую (

1) разность функций распределения на концах области определения при произвольных числах Кнудсена.

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

Рассмотрим теперь методы решения уравнения Больцмана, основанные на строгой кинетической теории.