Решеточные методы Больцмана

Решеточные методы Больцмана (LBM) , возникшие на основе метода решеточных газовых автоматов (LGA) (модели Харди- Помо -Пацциса и Фриша - Хасслахера - Помо ), представляют собой класс методов вычислительной гидродинамики (CFD) для моделирования жидкости . Вместо прямого решения уравнений Навье – Стокса плотность жидкости на решетке моделируется с помощью процессов течения и столкновения (релаксации). Метод универсален поскольку модельную жидкость можно легко имитировать обычное поведение жидкости, такое как сосуществование пара и жидкости, и таким образом можно моделировать жидкостные системы, такие как капли жидкости. Кроме того, жидкости в сложных средах, таких как пористые среды, можно напрямую моделировать, тогда как со сложными границами работать с другими методами CFD может быть сложно.

Продолжительность: 20 секунд. 0:20

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

Алгоритм

Схема векторов решетки D2Q9 для двумерной решетки Больцмана

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

В алгоритме есть этапы столкновения и потоковой передачи. Они изменяют плотность жидкости. $\rho ({\vec {x}},t)$, для ${\vec {x}}$ положение и $t$ время. Поскольку жидкость находится на решетке, плотность имеет ряд составляющих: $f_{i},i=0,\ldots ,a$ равно числу векторов решетки, связанных с каждой точкой решетки. В качестве примера здесь показаны векторы решетки для простой решетки, используемой в двумерном моделировании. Эту решетку обычно обозначают D2Q9, для двух измерений и девяти векторов: четыре вектора вдоль севера, востока, юга и запада, плюс четыре вектора по углам единичного квадрата, плюс вектор с обоими нулевыми компонентами. Тогда, например, вектор ${\vec {e}}_{4}=(0,-1)$, т. е. оно указывает строго на юг и поэтому не имеет $x$ компонент, но $y$ компонент $-1$. Итак, одна из девяти составляющих полной плотности в центральной точке решетки, $f_{4}({\vec {x}},t)$, является той частью жидкости в точке ${\vec {x}}$ двигаясь прямо на юг, со скоростью в единицах решетки.

Тогда этапы, которые развивают жидкость во времени, таковы:

Шаг столкновения

$f_{i}({\vec {x}},t+\delta _{t})=f_{i}({\vec {x}},t)+{\frac {f_{i}^{eq}({\vec {x}},t)-f_{i}({\vec {x}},t)}{\tau _{f}}}\,\!$

что такое Бхатнагар Гросс и Крук (BGK) модель релаксации к равновесию посредством столкновений между молекулами жидкости. $f_{i}^{eq}({\vec {x}},t)$ — равновесная плотность в направлении i при данной плотности тока, это можно выразить в приближении Тейлора (см. ниже, в разделе «Математические уравнения для моделирования »):

$f_{i}^{eq}=\omega _{i}\rho \left(1+{\frac {3{\vec {e}}_{i}{\vec {u}}}{c^{2}}}+{\frac {9({\vec {e}}_{i}{\vec {u}})^{2}}{2c^{4}}}-{\frac {3({\vec {u}})^{2}}{2c^{2}}}\right)$

Модель предполагает, что жидкость локально релаксирует до равновесия в течение характерного времени. $\tau _{f}$. Эта временная шкала определяет кинематическую вязкость : чем она больше, тем больше кинематическая вязкость.

Шаг потоковой передачи

$f_{i}({\vec {x}}+{\vec {e}}_{i},t+\delta _{t})=f_{i}({\vec {x}},t)\,\!$

Как $f_{i}({\vec {x}},t)$ это, по определению, плотность жидкости в точке ${\vec {x}}$ во время $t$, который движется со скоростью ${\vec {e}}_{i}$ за шаг по времени, затем на следующий шаг по времени $t+\delta _{t}$ это перетекло в точку ${\vec {x}}+{\vec {e}}_{i}$.

Преимущества

• LBM был разработан с нуля для эффективной работы на архитектурах с массовым параллелизмом , начиная от недорогих встроенных FPGA и DSP и заканчивая графическими процессорами , гетерогенными кластерами и суперкомпьютерами (даже с медленной сетью межсоединений). Это обеспечивает сложную физику и сложные алгоритмы. Эффективность выводит на качественно новый уровень понимания, поскольку позволяет решать задачи, к которым ранее нельзя было подойти (или лишь с недостаточной точностью).
• Этот метод основан на молекулярном описании жидкости и может напрямую включать физические термины, возникающие на основе знаний о взаимодействии между молекулами. Следовательно, это незаменимый инструмент фундаментальных исследований, поскольку он сокращает цикл между разработкой теории и формулировкой соответствующей численной модели.
• Автоматизированная предварительная обработка данных и создание решетки за время, которое составляет небольшую часть всего моделирования.
• Параллельный анализ данных, постобработка и оценка.
• Полностью разрешенный многофазный поток с мелкими каплями и пузырьками.
• Полностью разрешенный поток через сложную геометрию и пористую среду.
• Сложный, связанный поток с теплообменом и химическими реакциями.

и

Как и в случае с CFD на основе Навье-Стокса, методы LBM успешно сочетаются с термоспецифичными решениями, обеспечивая возможность моделирования теплопередачи (проводимость на основе твердых тел, конвекция и излучение). Для многофазных/многокомпонентных моделей толщина границы раздела обычно велика, а соотношение плотностей на границе раздела невелико по сравнению с реальными жидкостями. Недавно эта проблема была решена Юанем и Шефером, которые усовершенствовали модели Шаня и Чена, Свифта и Хэ, Чена и Чжана. Они смогли достичь соотношения плотностей 1000:1, просто изменив уравнение состояния . Было предложено применить преобразование Галилея, чтобы преодолеть ограничение моделирования высокоскоростных потоков жидкости. Быстрое развитие этого метода также позволило успешно смоделировать микрофлюидику . Однако на данный момент LBM все еще ограничен в моделировании с высокими числами Кнудсена течений методы Монте-Карло , где вместо этого используются , а потоки с высокими числами Маха в аэродинамике все еще сложны для LBM, а последовательная термогидродинамическая схема отсутствует.

Развитие по методу LGA

LBM возник из метода автоматов решеточного газа (LGA), который можно рассматривать как упрощенную фиктивную модель молекулярной динамики, в которой пространство, время и скорости частиц дискретны. Например, в двумерной модели FHP каждый узел решетки связан со своими соседями 6 скоростями решетки в треугольной решетке; в узле решетки может находиться либо 0, либо 1 частица, движущаяся с заданной скоростью решетки. Через промежуток времени каждая частица переместится в соседний узел в своем направлении; этот процесс называется этапом распространения или потоковой передачи. Когда более одной частицы прибывает в один и тот же узел с разных направлений, они сталкиваются и меняют свои скорости в соответствии с набором правил столкновений. Шаги потоковой передачи и шаги столкновения чередуются. Подходящие правила столкновений должны сохранять количество частиц (массу), импульс и энергию до и после столкновения. LGA страдает несколькими врожденными недостатками для использования в гидродинамическом моделировании: отсутствие галилеевой инвариантности. для быстрых потоков, статистического шума и плохого масштабирования числа Рейнольдса в зависимости от размера решетки. Однако LGA хорошо подходит для упрощения и расширения возможностей моделей реакционной диффузии и молекулярной динамики .

Основной мотивацией перехода от LGA к LBM было желание убрать статистический шум путем замены булевого числа частиц в направлении решетки на его среднее по ансамблю, так называемую функцию распределения плотности. Вместе с этой заменой правило дискретных коллизий также заменяется непрерывной функцией, известной как оператор коллизий. В разработке LBM важным упрощением является аппроксимация оператора столкновения с помощью релаксационного члена Бхатнагара-Гросса-Крука (BGK). Эта решетчатая модель BGK (LBGK) делает моделирование более эффективным и обеспечивает гибкость транспортных коэффициентов. С другой стороны, было показано, что схему LBM можно рассматривать и как специальную дискретизированную форму непрерывного уравнения Больцмана. Из теории Чепмена-Энскога можно восстановить основную непрерывность и уравнения Навье-Стокса из алгоритма LBM.

и n Q m классификация

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

Популярным способом классификации различных методов по решетке DnQm является схема . Здесь «D n » означает « n размеров», а «Q m » означает « m скоростей». Например, D3Q15 представляет собой трехмерную решетчатую модель Больцмана на кубической сетке с присутствием остальных частиц. Каждый узел имеет кристаллическую форму и может доставлять частицы в 15 узлов: в каждый из 6 соседних узлов, имеющих общую поверхность, в 8 соседних узлов, разделяющих угол, и в себя. (Модель D3Q15 не содержит частиц, движущихся к 12 соседним узлам, имеющим общее ребро; добавление их создаст модель «D3Q27».)

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

Преобразование единиц решетки

В большинстве решетчатых симуляций Больцмана $\delta _{x}\,\!$ является основной единицей шага решетки, поэтому, если область длины $L\,\!$ имеет $N\,\!$ единиц решетки по всей ее длине, пространственная единица определяется просто как $\delta _{x}=L/N\,\!$. Скорости в решеточном моделировании Больцмана обычно выражаются в терминах скорости звука. Таким образом, дискретную единицу времени можно определить как $\delta _{t}={\frac {\delta _{x}}{C_{s}}}\,\!$, где знаменатель $C_{s}$ это физическая скорость звука.

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

Моделирование смесей

Моделирование многофазных/многокомпонентных потоков всегда было проблемой для традиционной CFD из-за движущихся и деформируемых интерфейсов . Более фундаментально, границы раздела между различными фазами (жидкостью и паром) или компонентами (например, нефтью и водой) возникают в результате специфических взаимодействий между молекулами жидкости. Поэтому такие микроскопические взаимодействия сложно реализовать в макроскопическом уравнении Навье–Стокса. Однако в LBM кинетика частиц обеспечивает относительно простой и последовательный способ учета основных микроскопических взаимодействий путем изменения оператора столкновения. Было разработано несколько многофазных/многокомпонентных моделей LBM. Здесь фазовые разделения генерируются автоматически на основе динамики частиц, и для манипулирования границами раздела не требуется специальной обработки, как в традиционных методах CFD. Успешные применения многофазных/многокомпонентных моделей LBM можно найти в различных сложных жидкостных системах, включая нестабильность границы раздела, пузырьковые / капель динамика , смачивание на твердых поверхностях, межфазное скольжение и электрогидродинамические деформации капель.

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

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

Термическая решетка

В настоящее время (2009 г.) метод термической решетки-Больцмана (TLBM) относится к одной из трех категорий: многоскоростной подход, пассивный скалярный подход, и распределение тепловой энергии.

уравнения Навье – Стокса из

См. Также: Вывод уравнений Навье – Стокса.

Начнем с дискретного решеточного уравнения Больцмана (также называемого уравнением ЛБГК из-за используемого оператора столкновений). 2-го порядка Сначала мы разложим ряд Тейлора в левой части LBE. Это выбрано вместо более простого разложения Тейлора 1-го порядка, поскольку дискретный LBE не может быть восстановлен. При выполнении разложения в ряд Тейлора 2-го порядка член нулевой производной и первый член справа сокращаются, оставляя только первую и вторую производные члены разложения Тейлора и оператор столкновения:

$f_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})=f_{i}({\vec {x}},t)+{\frac {\delta _{t}}{\tau _{f}}}(f_{i}^{eq}-f_{i}).$

Для простоты напишите $f_{i}({\vec {x}},t)$ как $f_{i}$. Немного упрощенное разложение ряда Тейлора тогда выглядит следующим образом, где «:» — произведение двоеточия между диадами:

${\frac {\partial f_{i}}{\partial t}}+{\vec {e}}_{i}\cdot \nabla f_{i}+\left({\frac {1}{2}}{\vec {e}}_{i}{\vec {e}}_{i}:\nabla \nabla f_{i}+{\vec {e}}_{i}\cdot \nabla {\frac {\partial f_{i}}{\partial t}}+{\frac {1}{2}}{\frac {\partial ^{2}f_{i}}{\partial t^{2}}}\right)={\frac {1}{\tau }}(f_{i}^{eq}-f_{i}).$

Разложив функцию распределения частиц на равновесную и неравновесную составляющие и воспользовавшись разложением Чепмена-Энскога, где $K$ - число Кнудсена, расширенный Тейлором LBE можно разложить на различные по порядку величины числа Кнудсена, чтобы получить правильные уравнения непрерывной среды:

$f_{i}=f_{i}^{\text{eq}}+Kf_{i}^{\text{neq}},$

$f_{i}^{\text{neq}}=f_{i}^{(1)}+Kf_{i}^{(2)}+O(K^{2}).$

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

$\rho =\sum _{i}f_{i}^{\text{eq}},$

$\rho {\vec {u}}=\sum _{i}f_{i}^{\text{eq}}{\vec {e}}_{i},$

$0=\sum _{i}f_{i}^{(k)}\qquad {\text{for }}k=1,2,$

$0=\sum _{i}f_{i}^{(k)}{\vec {e}}_{i}.$

Тогда расширение Чепмена-Энскога будет:

${\frac {\partial }{\partial t}}=K{\frac {\partial }{\partial t_{1}}}+K^{2}{\frac {\partial }{\partial t_{2}}}\qquad {\text{for }}t_{2}({\text{diffusive time-scale}})\ll t_{1}({\text{convective time-scale}}),$

${\frac {\partial }{\partial x}}=K{\frac {\partial }{\partial x_{1}}}.$

Путем подстановки расширенного равновесия и неравновесия в расширение Тейлора и разделения на разные порядки $K$, уравнения непрерывной среды почти выведены.

Для заказа $K^{0}$:

${\frac {\partial f_{i}^{\text{eq}}}{\partial t_{1}}}+{\vec {e}}_{i}\nabla _{1}f_{i}^{\text{eq}}=-{\frac {f_{i}^{(1)}}{\tau }}.$

Для заказа $K^{1}$:

${\frac {\partial f_{i}^{(1)}}{\partial t_{1}}}+{\frac {\partial f_{i}^{\text{eq}}}{\partial t_{2}}}+{\vec {e}}_{i}\nabla f_{i}^{(1)}+{\frac {1}{2}}{\vec {e}}_{i}{\vec {e}}_{i}:\nabla \nabla f_{i}^{\text{eq}}+{\vec {e}}_{i}\cdot \nabla {\frac {\partial f_{i}^{\text{eq}}}{\partial t_{1}}}+{\frac {1}{2}}{\frac {\partial ^{2}f_{i}^{\text{eq}}}{\partial t_{1}^{2}}}=-{\frac {f_{i}^{(2)}}{\tau }}.$

Затем второе уравнение можно упростить с помощью некоторой алгебры, а первое уравнение — до следующего:

${\frac {\partial f_{i}^{\text{eq}}}{\partial t_{2}}}+\left(1-{\frac {1}{2\tau }}\right)\left[{\frac {\partial f_{i}^{(1)}}{\partial t_{1}}}+{\vec {e}}_{i}\nabla _{1}f_{i}^{(1)}\right]=-{\frac {f_{i}^{(2)}}{\tau }}.$

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

${\frac {\partial \rho }{\partial t}}+\nabla \cdot \rho {\vec {u}}=0,$

${\frac {\partial \rho {\vec {u}}}{\partial t}}+\nabla \cdot \Pi =0.$

Тензор потока импульса $\Pi$ тогда имеет следующий вид:

$\Pi _{xy}=\sum _{i}{\vec {e}}_{ix}{\vec {e}}_{iy}\left[f_{i}^{eq}+\left(1-{\frac {1}{2\tau }}\right)f_{i}^{(1)}\right],$

где ${\vec {e}}_{ix}{\vec {e}}_{iy}$ является сокращением квадрата суммы всех компонентов ${\vec {e}}_{i}$ (т.е. $\textstyle \left(\sum _{x}{\vec {e}}_{ix}\right)^{2}=\sum _{x}\sum _{y}{\vec {e}}_{ix}{\vec {e}}_{iy}$), а равновесное распределение частиц второго порядка, сравнимое с уравнением Навье – Стокса, имеет вид:

$f_{i}^{\text{eq}}=\omega _{i}\rho \left(1+{\frac {{\vec {e}}_{i}{\vec {u}}}{c_{s}^{2}}}+{\frac {({\vec {e}}_{i}{\vec {u}})^{2}}{2c_{s}^{4}}}-{\frac {{\vec {u}}^{2}}{2c_{s}^{2}}}\right).$

Равновесное распределение справедливо только для малых скоростей или малых чисел Маха . Вставка равновесного распределения обратно в тензор потока приводит к:

$\Pi _{xy}^{(0)}=\sum _{i}{\vec {e}}_{ix}{\vec {e}}_{iy}f_{i}^{eq}=p\delta _{xy}+\rho u_{x}u_{y},$

$\Pi _{xy}^{(1)}=\left(1-{\frac {1}{2\tau }}\right)\sum _{i}{\vec {e}}_{ix}{\vec {e}}_{iy}f_{i}^{(1)}=\nu \left(\nabla _{x}\left(\rho {\vec {u}}_{y}\right)+\nabla _{y}\left(\rho {\vec {u}}_{x}\right)\right).$

Наконец, уравнение Навье – Стокса восстанавливается в предположении, что изменение плотности мало:

$\rho \left({\frac {\partial {\vec {u}}_{x}}{\partial t}}+\nabla _{y}\cdot {\vec {u}}_{x}{\vec {u}}_{y}\right)=-\nabla _{x}p+\nu \nabla _{y}\cdot \left(\nabla _{x}\left(\rho {\vec {u}}_{y}\right)+\nabla _{y}\left(\rho {\vec {u}}_{x}\right)\right).$

Этот вывод следует за работой Чена и Дулена.

Математические уравнения для моделирования

Непрерывное уравнение Больцмана представляет собой уравнение эволюции для функции распределения вероятностей одной частицы. $f({\vec {x}},{\vec {e}}_{i},t)$ и функция распределения плотности внутренней энергии $g({\vec {x}},{\vec {e}}_{i},t)$ (Он и др.) соответственно:

$\partial _{t}f+({\vec {e}}\cdot \nabla )f+F\partial _{v}f=\Omega (f),$

$\partial _{t}g+({\vec {e}}\cdot \nabla )g+G\partial _{v}f=\Omega (g),$

где $g({\vec {x}},{\vec {e}}_{i},t)$ связано с $f({\vec {x}},{\vec {e}}_{i},t)$ к

$g({\vec {x}},{\vec {e}}_{i},t)={\frac {({\vec {e}}-{\vec {u}})^{2}}{2}}f({\vec {x}},{\vec {e}}_{i},t),$

$F$ это внешняя сила, $\Omega$ является интегралом столкновений, а ${\vec {e}}$ (также обозначено ${\vec {\xi }}$ в литературе) — микроскопическая скорость. Внешняя сила $F$ связано с температурой внешней силы $G$ по соотношению ниже. Типичным тестом модели является конвекция Рэлея – Бенара для $G$.

$F={\frac {{\vec {G}}\cdot ({\vec {e}}-{\vec {u}})}{RT}}f^{\text{eq}},$

${\vec {G}}=\beta g_{0}(T-T_{avg}){\vec {k}}.$

Макроскопические переменные, такие как плотность $\rho$, скорость ${\vec {u}}$и температура $T$ могут быть рассчитаны как моменты функции распределения плотности:

$\rho =\int f\,d{\vec {e}},$

$\rho {\vec {u}}=\int {\vec {e}}f\,d{\vec {e}},$

${\frac {\rho DRT}{2}}=\rho \epsilon =\int g\,d{\vec {e}}.$

Решеточный метод Больцмана дискретизирует это уравнение, ограничивая пространство решеткой, а пространство скоростей дискретным набором микроскопических скоростей (т. е. ${\vec {e}}_{i}=({\vec {e}}_{ix},{\vec {e}}_{iy})$). Например, микроскопические скорости в D2Q9, D3Q15 и D3Q19 задаются как:

${\vec {e}}_{i}=c\times {\begin{cases}(0,0)&i=0\\(1,0),(0,1),(-1,0),(0,-1)&i=1,2,3,4\\(1,1),(-1,1),(-1,-1),(1,-1)&i=5,6,7,8\\\end{cases}}$

${\vec {e}}_{i}=c\times {\begin{cases}(0,0,0)&i=0\\(\pm 1,0,0),(0,\pm 1,0),(0,0,\pm 1)&i=1,2,...,5,6\\(\pm 1,\pm 1,\pm 1)&i=7,8,...,13,14\\\end{cases}}$

${\vec {e}}_{i}=c\times {\begin{cases}(0,0,0)&i=0\\(\pm 1,0,0),(0,\pm 1,0),(0,0,\pm 1)&i=1,2,...,5,6\\(\pm 1,\pm 1,0),(\pm 1,0,\pm 1),(0,\pm 1,\pm 1)&i=7,8,...,17,18\\\end{cases}}$

Однофазное дискретизированное уравнение Больцмана для плотности массы и плотности внутренней энергии:

$f_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})-f_{i}({\vec {x}},t)+F_{i}=\Omega (f),$

$g_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})-g_{i}({\vec {x}},t)+G_{i}=\Omega (g).$

Оператор столкновения часто аппроксимируется оператором столкновения BGK при условии, что он также удовлетворяет законам сохранения:

$\Omega (f)={\frac {1}{\tau _{f}}}(f_{i}^{\text{eq}}-f_{i}),$

$\Omega (g)={\frac {1}{\tau _{g}}}(g_{i}^{\text{eq}}-g_{i}).$

В операторе столкновения $f_{i}^{\text{eq}}$ — дискретная равновесная функция распределения вероятностей частиц . В D2Q9 и D3Q19 это показано ниже для несжимаемого потока в непрерывной и дискретной форме, где D , R и T — размерность, универсальная газовая постоянная и абсолютная температура соответственно. Частичный вывод от непрерывной формы к дискретной обеспечивается посредством простого вывода со вторым порядком точности.

$f^{\text{eq}}={\frac {\rho }{(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}}-{\vec {u}})^{2}}{2RT}}}$

$={\frac {\rho }{(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}})^{2}}{2RT}}}e^{{\frac {{\vec {e}}{\vec {u}}}{RT}}-{\frac {{\vec {u}}^{2}}{2RT}}}$

$={\frac {\rho }{(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}})^{2}}{2RT}}}\left(1+{\frac {{\vec {e}}{\vec {u}}}{RT}}+{\frac {({\vec {e}}{\vec {u}})^{2}}{2(RT)^{2}}}-{\frac {{\vec {u}}^{2}}{2RT}}+...\right)$

Сдача в аренду $c={\sqrt {3RT}}$ дает окончательный результат:

$f_{i}^{eq}=\omega _{i}\rho \left(1+{\frac {3{\vec {e}}_{i}{\vec {u}}}{c^{2}}}+{\frac {9({\vec {e}}_{i}{\vec {u}})^{2}}{2c^{4}}}-{\frac {3({\vec {u}})^{2}}{2c^{2}}}\right)$

$g^{eq}={\frac {\rho ({\vec {e}}-{\vec {u}})^{2}}{2(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}}-{\vec {u}})^{2}}{2RT}}}$

$\omega _{i}={\begin{cases}4/9&i=0\\1/9&i=1,2,3,4\\1/36&i=5,6,7,8\\\end{cases}}$

$\omega _{i}={\begin{cases}1/3&i=0\\1/18&i=1,2,...,5,6\\1/36&i=7,8,...,17,18\\\end{cases}}$

Поскольку над однокомпонентным потоком уже проделана большая работа, будет обсуждаться следующий TLBM. Многокомпонентный/многофазный TLBM также более интересен и полезен, чем просто один компонент. Чтобы соответствовать текущим исследованиям, определите набор всех компонентов системы (т. е. стенок пористой среды, множества жидкостей/газов и т. д.) $\Psi$ с элементами $\sigma _{j}$.

$f_{i}^{\sigma }({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})-f_{i}^{\sigma }({\vec {x}},t)+F_{i}={\frac {1}{\tau _{f}^{\sigma }}}(f_{i}^{\sigma ,eq}(\rho ^{\sigma },v^{\sigma })-f_{i}^{\sigma })$

Параметр релаксации, $\tau _{f}^{\sigma _{j}}\,\!$, связано с кинематической вязкостью , $\nu _{f}^{\sigma _{j}}\,\!$, по следующему соотношению:

$\nu _{f}^{\sigma _{j}}=(\tau _{f}^{\sigma _{j}}-0.5)c_{s}^{2}\delta _{t}.$

Моменты ​ ​ $f_{i}\,\!$ укажите локальные сохраняющиеся величины. Плотность определяется выражением

$\rho =\sum _{\sigma }\sum _{i}f_{i}\,\!$

$\rho \epsilon =\sum _{i}g_{i}\,\!$

$\rho ^{\sigma }=\sum _{i}f_{i}^{\sigma }\,\!$

и средневзвешенная скорость, ${\vec {u'}}\,\!$, а локальный импульс определяется выражением

${\vec {u'}}=\left(\sum _{\sigma }{\frac {\rho ^{\sigma }{\vec {u^{\sigma }}}}{\tau _{f}^{\sigma }}}\right)/\left(\sum _{\sigma }{\frac {\rho ^{\sigma }}{\tau _{f}^{\sigma }}}\right)$

$\rho ^{\sigma }{\vec {u^{\sigma }}}=\sum _{i}f_{i}^{\sigma }{\vec {e}}_{i}.$

$v^{\sigma }={\vec {u'}}+{\frac {\tau _{f}^{\sigma }}{\rho ^{\sigma }}}{\vec {F}}^{\sigma }$

В приведенном выше уравнении для равновесной скорости $v^{\sigma }\,\!$, ${\vec {F}}^{\sigma }\,\!$ Термин представляет собой силу взаимодействия между компонентом и другими компонентами. Это до сих пор является предметом многочисленных дискуссий, поскольку обычно это параметр настройки, который определяет, как взаимодействуют жидкость-жидкость, жидкость-газ и т. д. Франк и др. перечислите текущие модели для этого термина силы. Обычно используемые выводы - это хромодинамическая модель Ганстенсена, подход Свифта, основанный на свободной энергии как для систем жидкость/пар, так и для бинарных жидкостей, модель Хе, основанная на межмолекулярном взаимодействии, подход Инамуро и подход Ли и Линя.

Ниже приводится общее описание ${\vec {F}}^{\sigma }\,\!$ как указано рядом авторов.

${\vec {F}}^{\sigma }=-\psi ^{\sigma }({\vec {x}})\sum _{\sigma _{j}}H^{\sigma \sigma _{j}}({\vec {x}},{\vec {x}}')\sum _{i}\psi ^{\sigma _{j}}({\vec {x}}+{\vec {e}}_{i}){\vec {e}}_{i}\,\!$

$\psi ({\vec {x}})\,\!$ - эффективная масса и $H({\vec {x}},{\vec {x}}')\,\!$ – функция Грина, представляющая межчастичное взаимодействие с ${\vec {x}}'\,\!$ как соседний сайт. Удовлетворение $H({\vec {x}},{\vec {x}}')=H({\vec {x}}',{\vec {x}})\,\!$ и где $H({\vec {x}},{\vec {x}}')>0\,\!$ представляет собой отталкивающие силы. Для D2Q9 и D3Q19 это приводит к

$H^{\sigma \sigma _{j}}({\vec {x}},{\vec {x}}')={\begin{cases}h^{\sigma \sigma _{j}}&\left|{\vec {x}}-{\vec {x}}'\right|\leq c\\0&\left|{\vec {x}}-{\vec {x}}'\right|>c\\\end{cases}}$

$H^{\sigma \sigma _{j}}({\vec {x}},{\vec {x}}')={\begin{cases}h^{\sigma \sigma _{j}}&\left|{\vec {x}}-{\vec {x}}'\right|=c\\h^{\sigma \sigma _{j}}/2&\left|{\vec {x}}-{\vec {x}}'\right|={\sqrt {2c}}\\0&{\text{otherwise}}\\\end{cases}}$

Эффективная масса, предложенная Шанем и Ченом, использует следующую эффективную массу для однокомпонентной многофазной системы . Уравнение состояния также приведено в условии однокомпонентности и многофазности.

$\psi ({\vec {x}})=\psi (\rho ^{\sigma })=\rho _{0}^{\sigma }\left[1-e^{(-\rho ^{\sigma }/\rho _{0}^{\sigma })}\right]\,\!$

$p=c_{s}^{2}\rho +c_{0}h[\psi ({\vec {x}})]^{2}\,\!$

До сих пор кажется, что $\rho _{0}^{\sigma }\,\!$ и $h^{\sigma \sigma _{j}}\,\!$ системы являются свободными константами для настройки, но после включения в уравнение состояния (EOS) они должны удовлетворять термодинамическим соотношениям в критической точке, так что $(\partial P/\partial {\rho })_{T}=(\partial ^{2}P/\partial {\rho ^{2}})_{T}=0\,\!$ и $p=p_{c}\,\!$. Для ЭОС, $c_{0}\,\!$ составляет 3,0 для D2Q9 и D3Q19, а для D3Q15 он равен 10,0.

Позже это было показано Юанем и Шефером. что эффективную массовую плотность необходимо изменить для более точного моделирования многофазного потока. Они сравнили Шан и Чен (SC), Карнахан-Старлинг (C-S), Ван-дер-Ваальса (vdW), Редлиха-Квонга (R-K), Редлиха-Квонга Соаве (RKS) и Пенга-Робинсона (P- Р) ЭОС. Их результаты показали, что SC EOS было недостаточно и что C–S, P–R, R–K и RKS EOS более точны при моделировании многофазного потока одного компонента.

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

$\rho \theta +\rho uu=\sum _{i}f_{i}{\vec {e}}_{i}{\vec {e}}_{i}.$

Неструктурированные сетки

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

Предполагая $\Omega ^{j}$ - это объем, составленный всеми барицентрами тетраэдров, гранями и ребрами, соединенными с вершиной. ${\boldsymbol {v}}^{j}$, дискретная функция плотности скорости:

$f_{i}({\boldsymbol {v}}^{j},t+\delta t)=f_{i}({\boldsymbol {v}}^{j},t)-\delta t\sum _{k}S_{i}^{jk}f_{i}({\boldsymbol {v}}^{k},t)-{\delta t \over \tau }\sum _{k}C^{jk}(f_{i}({\boldsymbol {v}}^{k},t)-f_{i}^{eq}({\boldsymbol {v}}^{k}))$

где ${\boldsymbol {v}}^{k}$ - положение вершины и ее соседей, а также:

$C^{jk}={1 \over V^{j}}\int _{\Omega ^{j}}w_{k}({\boldsymbol {x}})d\Omega$

$S_{i}^{jk}={1 \over V^{j}}\oint _{\partial \Omega ^{j}}({\vec {e_{i}}}{\vec {n}})w_{k}({\boldsymbol {x}})d\Omega$

где $w_{k}({\boldsymbol {x}})$ представляет собой результат линейной интерполяции ${\boldsymbol {x}}$ вершинами треугольника или тетраэдра, которые ${\boldsymbol {x}}$ лежит внутри.

Приложения

За последние годы LBM зарекомендовал себя как мощный инструмент для решения проблем разной длины и временного масштаба. Некоторые из применений LBM включают в себя:

• Пористые среды
• Биомедицинские потоки
• Науки о Земле (Фильтрация почвы).
• Энергетические науки (Топливные элементы ).

Пример реализации

Это базовая реализация LBM на сетке 100x100 с использованием Python:

#This is a fluid simulator using the lattice Boltzmann method.
#Using D2Q9 and peiodic boundary, and used no external library.
#It generates two ripples at 50,50 and 50,40.
#Reference: Erlend Magnus Viggen's Master thesis, "The Lattice Boltzmann Method with Applications in Acoustics".
#For Wikipedia under CC-BY-SA license.
import math
#Define some utilities
def sum(a):
    s=0
    for e in a:
        s=s+e
    return s
#Weights in D2Q9
Weights=[1/36,1/9,1/36,
   1/9, 4/9,1/9,
   1/36,1/9,1/36]
#Discrete velocity vectors
DiscreteVelocityVectors=[[-1,1],[0,1],[1,1],
     [-1,0],[0,0],[1,0],
     [-1,-1],[0,-1],[1,-1]
]
#A Field2D class
class Field2D():
    def __init__(self,res : int):
        self.field=[]
        for b in range(res):
            fm=[]
            for a in range(res):
                fm.append([0,0,0,
                           0,1,0,
                           0,0,0])
            self.field.append(fm[:])
        self.res = res
    #This visualize the simulation, can only be used in a terminal
    @staticmethod
    def VisualizeField(a,sc,res):
        stringr=""
        for u in range(res):
            row=""
            for v in range(res):
                n=int(u*a.res/res)
                x=int(v*a.res/res)
                flowmomentem=a.Momentum(n,x)
                col="\033[38;2;{0};{1};{2}m██".format(int(127+sc*flowmomentem[0]),int(127+sc*flowmomentem[1]),0)
                row=row+col
            print(row)
            stringr=stringr+row+"\n"
        return stringr
    #Momentum of the field
    def Momentum(self,x,y):
        return velocityField[y][x][0]*sum(self.field[y][x]),velocityField[y][x][1]*sum(self.field[y][x])
#Resolution of the simulation
res=100
a=Field2D(res)
#The velocity field
velocityField=[]
for DummyVariable in range(res):
    DummyList=[]
    for DummyVariable2 in range(res):
        DummyList.append([0,0])
    velocityField.append(DummyList[:])
#The density field
DensityField=[]
for DummyVariable in range(res):
    DummyList=[]
    for DummyVariable2 in range(res):
        DummyList.append(1)
    DensityField.append(DummyList[:])
#Set initial condition
DensityField[50][50]=2
DensityField[40][50]=2
#Maximum solving steps
MaxSteps = 120
#The speed of sound, specifically 1/sqrt(3) ~ 0.57
SpeedOfSound=1/math.sqrt(3)
#time relaxation constant
TimeRelaxationConstant=0.5
#Solve
for s in range(MaxSteps):
    #Collision Step
    df=Field2D(res)
    for y in range(res):
        for x in range(res):
            for v in range(9):
                Velocity=a.field[y][x][v]
                FirstTerm=Velocity
                #The Flow Velocity
                FlowVelocity=velocityField[y][x]
                Dotted=FlowVelocity[0]*DiscreteVelocityVectors[v][0]+FlowVelocity[1]*DiscreteVelocityVectors[v][1]
                # #The taylor expainsion of equilibrium term
                taylor=1+((Dotted)/(SpeedOfSound**2))+((Dotted**2)/(2*SpeedOfSound**4))-((FlowVelocity[0]**2+FlowVelocity[1]**2)/(2*SpeedOfSound**2))
                #The current density
                density=DensityField[y][x]
                #The equilibrium
                equilibrium=density*taylor*Weights[v]
                SecondTerm=(equilibrium-Velocity)/TimeRelaxationConstant
                df.field[y][x][v]=FirstTerm+SecondTerm
    #Streaming Step
    for y in range(0,res):
        for x in range(0,res):
            for v in range(9):
                #Target, the lattice point this iteration is solving
                TargetY=y+DiscreteVelocityVectors[v][1]
                TargetX=x+DiscreteVelocityVectors[v][0]
                # Peiodic Boundary
                if TargetY == res and TargetX == res:
                    a.field[TargetY-res][TargetX-res][v]=df.field[y][x][v]
                elif TargetX == res:
                    a.field[TargetY][TargetX-res][v]=df.field[y][x][v]
                elif TargetY == res:
                    a.field[TargetY-res][TargetX][v]=df.field[y][x][v]
                elif TargetY == -1 and TargetX == -1:
                    a.field[TargetY+res][TargetX+res][v]=df.field[y][x][v]   
                elif TargetX == -1:
                    a.field[TargetY][TargetX+res][v]=df.field[y][x][v]
                elif TargetY == -1:
                    a.field[TargetY+res][TargetX][v]=df.field[y][x][v]
                else:
                    a.field[TargetY][TargetX][v]=df.field[y][x][v]
    #Calulate macroscopic variables
    for y in range(res):
        for x in range(res):
            #Recompute Density Field
            DensityField[y][x]=sum(a.field[y][x])
            #Recompute Flow Velocity
            FlowVelocity=[0,0]
            for DummyVariable in range(9):
                FlowVelocity[0]=FlowVelocity[0]+DiscreteVelocityVectors[DummyVariable][0]*a.field[y][x][DummyVariable]
            for DummyVariable in range(9):
                FlowVelocity[1]=FlowVelocity[1]+DiscreteVelocityVectors[DummyVariable][1]*a.field[y][x][DummyVariable]
            FlowVelocity[0]=FlowVelocity[0]/DensityField[y][x]
            FlowVelocity[1]=FlowVelocity[1]/DensityField[y][x]
            #Insert to Velocity Field
            velocityField[y][x]=FlowVelocity
    #Visualize
    Field2D.VisualizeField(a,128,100)

Внешние ссылки

• Метод LBM
• Энтропийный решеточный метод Больцмана (ELBM)
• dsfd.org: веб-сайт серии ежегодных конференций DSFD (с 1986 г. по настоящее время), на которых обсуждаются достижения в теории и применении решеточного метода Больцмана.
• Сайт ежегодной конференции ICMMES, посвящённой решеточным методам Больцмана и их приложениям.

Дальнейшее чтение

• Дойч, Андреас; Сабина Дорманн (2004). Клеточно-автоматное моделирование формирования биологических паттернов . Биркхойзер Верлаг . ISBN  978-0-8176-4281-5 .
• Суччи, Сауро (2001). Решеточное уравнение Больцмана для гидродинамики и не только . Издательство Оксфордского университета . ISBN  978-0-19-850398-9 .
• Вольф-Гладроу, Дитер (2000). Решеточные газовые клеточные автоматы и решеточные модели Больцмана . Спрингер Верлаг . ISBN  978-3-540-66973-9 .
• Сукоп, Майкл С.; Дэниел Т. Торн-младший (2007). Решётчатое Больцмановское моделирование: введение для геологов и инженеров . Спрингер . ISBN  978-3-540-27981-5 .
• Цзянь Го Чжоу (2004). Решеточные методы Больцмана для мелководных течений . Спрингер . ISBN  978-3-540-40746-1 .
• Хе, X., Чен, С., Дулен, Г. (1998). Новая тепловая модель для решеточного метода Больцмана в пределе несжимаемости . Академическая пресса . {{cite book}}: CS1 maint: несколько имен: список авторов ( ссылка )
• Го, ЗЛ; Шу, К. (2013). Решетчатый метод Больцмана и его применение в технике . Мировое научное издательство .
• Хуанг, Х.; МЦ Сукоп; ХY. Лу (2015). Многофазные решеточные методы Больцмана: теория и применение . Уайли-Блэквелл . ISBN  978-1-118-97133-8 .
• Крюгер, Т.; Кусуматмая, Х.; Кузьмин А.; Шардт, О.; Сильва, Г.; Вигген, Э.М. (2017). Решетчатый метод Больцмана: принципы и практика . Спрингер Издательство ISBN  978-3-319-44647-9 .