1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 63, № 12, 2023

Журнал вычислительной математики и математической физики, 2023, T. 63, № 12, стр. 2051-2065

Трехмерное моделирование высокоскоростного движения тела в трубе с разреженным газом

А. Н. Якунчиков 1, *, А. Р. Юлдашева 1, **

1 Механико-математический факультет
119991 Москва, Ленинские горы, МГУ, 1, Россия

* E-mail: art-ya@mail.ru
** E-mail: a.r.iuldasheva@gmail.com

Поступила в редакцию 17.07.2023
После доработки 25.07.2023
Принята к публикации 22.08.2023

Полный текст (PDF)

Аннотация

В работе исследуется обтекание тела, которое движется с высокой дозвуковой скоростью в трубе, заполненной разреженным газом. Данная аэродинамическая задача рассматривается применительно к проблеме создания высокоскоростного вакуумного транспорта в условиях конечных чисел Кнудсена. Выбраны параметры, максимально приближенные к целевым характеристикам проектов таких систем, а именно: скорость порядка 1000 км/ч, значительный относительный поперечный размер тела, смесь азота и кислорода (воздух) в качестве газа. Задача решалась в трехмерной постановке. Библ. 31. Фиг. 7. Табл. 3.

Ключевые слова: разреженный газ, вакуумный транспорт, аэродинамическое сопротивление, событийное молекулярно-динамическое моделирование, Hyperloop, EDMD.

ВВЕДЕНИЕ

Задача настоящего исследования тесно связана с идеей высокоскоростного вакуумного транспорта, которая впервые была сформулирована томским ученым Вейнбергом [1] еще в 1914 г. Данная идея была развита Маском [2] в его концепции вакуумной транспортной системы, которая получила название Hyperloop. В ее основе лежит создание сети наземных трубопроводов, внутри которых при пониженном давлении (0.001–0.01 атм) с высокой скоростью (порядка 1000 км/ч) двигаются капсулы с пассажирами или грузами. На данный момент несколько компаний, включая Hyperloop TT и Virgin Hyperloop, продолжают изучение конструкций и технологий, которые бы обеспечили безопасную работу такой транспортной системы. Например, компания Hyperloop One (в настоящее время – Virgin Hyperloop One), основанная в 2014 г., спроектировала полномасштабную пассажирскую капсулу длиной 8.7 м под названием XP-1 и провела эксперимент в трубе длиной 532 м в 2017 г., на тот момент без пассажиров [3]. Максимальная достигнутая скорость в ходе этого эксперимента составила 387 км/ч. Уже в 2020 г. компания провела тестовый запуск с пассажирами внутри обновленной капсулы XP-2. При давлении в трубе порядка 100 Па удалось достичь скорости 107 км/ч [3]. Компания Hyperloop TT также анонсировала полномасштабную пассажирскую капсулу длиной 32 м под названием Quintero One и провела аэродинамические эксперименты в трубе длиной 320 м [4]. Две другие компании, AIRT-IFICIAL в Испании и TransPod в Канаде, также уже несколько лет работают над подобными системами [5].

Были проведены многочисленные теоретические исследования, направленные на изучение аэродинамических характеристик таких систем, а также поиск оптимальной формы капсулы. Исследования [6–12] посвящены моделированию обтекания капсулы в континуальном режиме в широком диапазоне скоростей движения. Расчеты показали, что аэродинамическое сопротивление капсулы значительно возрастает с увеличением скорости и отношения поперечного размера капсулы к диаметру трубы, а также давления в трубе, тогда как влияние длины капсулы и температуры ее поверхности менее существенно. Условия расчетов перечисленных работ собраны в табл. 1.

Таблица 1.  

Условия теоретических и экспериментальных исследований, имеющихся в литературе

Источник Давление, Па Скорость, м/с Поперечное сечение тела Подход
[6, 10] 103–105 417 0.1–0.5 Континуальный
[7] 103 300 0.1–0.4 Континуальный
[8] 102 170–350 0.01–0.3 Континуальный
[9] 102–103 25–350 0.25, 0.35 Континуальный
[11] 103–105 139, 194 0.2–0.75 Континуальный
[12] 102 100–350 0.25, 0.36 Континуальный
[13] 104 14 – Эксперимент
[14] 105 170–520 0.08 Эксперимент
[15] 102 160, 320 0.34 Эксперимент

Также в литературе имеются результаты экспериментальных исследований. В работе [13] удалось продемонстрировать работоспособность транспортной системы на магнитном подвесе при пониженном давлении (0.1 атм), но пока для совсем небольших скоростей движения (около 50 км/ч). В работе [14] проведено экспериментальное исследование движения снаряда, летящего в трубе, в широком диапазоне скоростей. Однако измерения проводились при атмосферном давлении и очень маленьком относительном размере тела (приблизительно 0.08). Было получено, что при этих условиях коэффициент сопротивления имеет максимум (${{C}_{x}} \approx 2$) при числе Маха ${\text{M}} \approx 0.8$. Совсем недавно был проведен эксперимент [15] о движении капсулы в вакуумной трубе (0.001 атм) с двумя скоростями: 160 и 320 м/с, при этом относительный размер капсулы составил 0.34 от поперечного сечения трубы. К сожалению, коэффициенты сопротивления не высчитаны из эксперимента, а внимание было направлено на измерение скорости ударных волн, возникающих при влете капсулы в трубу.

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

1. ПОСТАНОВКА ЗАДАЧИ

Рассматривается задача о движении тела в цилиндрической трубе диаметра $D$. Труба заполнена покоящимся газом (смесью газов) с давлением ${{p}_{0}}$ и температурой ${{T}_{0}}$. Тело движется с постоянной скоростью ${{{v}}_{w}}$, положение тела в сечении трубы фиксировано и определяется смещением ${{y}_{0}}$ (см. фиг. 1). Рассматриваются режимы течения, соответствующие конечным числам Кнудсена и высоким дозвуковым скоростям движения. Температура поверхности трубы и тела предполагалась постоянной и равной ${{T}_{w}} = {{T}_{0}}$. Взаимодействие молекул газа с твердыми поверхностями моделировалось ядром рассеяния с полной аккомодацией энергии и импульса. В этом случае плотность вероятности для скорости после отражения будет

$\rho \left( {\mathbf{u}} \right) = \frac{{2\beta _{w}^{4}{{u}_{n}}}}{\pi }{\text{exp}}\left( { - \beta _{w}^{2}{{u}^{2}}} \right),$
где ${{\beta }_{w}} = \sqrt {m{\text{/}}2k{{T}_{w}}} $, $k$ – константа Больцмана, $m$ – масса молекулы. Были введены следующие безразмерные параметры задачи:

Фиг. 1.

Схема задачи.

• относительный поперечный размер тела $d{\text{/}}D$,

• безразмерное смещение ${{y}_{0}}{\text{/}}D$,

• параметр разрежения $\delta = D{{p}_{0}}{{\beta }_{m}}{\text{/}}{{\mu }_{m}}$,

• безразмерная скорость тела (число Маха) ${\text{M}} = {{{v}}_{w}}{\text{/}}\sqrt {\gamma {{R}_{m}}{{T}_{0}}} $,

где ${{\beta }_{m}} = 1{\text{/}}\sqrt {2{{R}_{m}}{{T}_{0}}} $, ${{\mu }_{m}}$ – вязкость смеси газов при температуре $T = {{T}_{0}}$, ${{R}_{m}}$ – газовая постоянная для смеси газов, $\gamma $ – постоянная адиабаты (бралось значение $\gamma = 1.4$, так как рассматривалась смесь двухатомных газов). Сила сопротивления ${{R}_{x}}$ и подъемная сила ${{R}_{y}}$ (искомые величины) обезразмеривались стандартным образом:

• коэффициент сопротивления ${{C}_{x}} = {{R}_{x}}{\text{/}}\left( {{{\rho }_{0}}{v}_{w}^{2}S{\text{/}}2} \right)$,

• коэффициент подъемной силы ${{C}_{y}} = {{R}_{y}}{\text{/}}\left( {{{\rho }_{0}}{v}_{w}^{2}S{\text{/}}2} \right)$,

где ${{\rho }_{0}} = {{p}_{0}}{\text{/}}\left( {{{R}_{m}}{{T}_{0}}} \right)$ – плотность покоящейся смеси, $S = \pi {{d}^{2}}{\text{/}}4$ – площадь поперечного сечения тела.

2. ЧИСЛЕННЫЙ МЕТОД

Для решения задачи использовался метод событийного молекулярно-динамического моделирования (EDMD) [16–20], который является компромиссом между прямым статистическим моделированием Монте-Карло (DSMC) и классическим молекулярно-динамическим моделированием (MD) и по уровню детализации модели, и по вычислительной нагрузке. Событийный метод является полностью детерминированной и бессеточной (и по пространству, и по времени) реализацией выбранной модели столкновений молекул без стохастического выбора пары для столкновения (как в DSMC) или упрощения интеграла столкновений (как в методах с модельным интегралом столкновений). Таким образом, нет необходимости проводить большие серии тестовых расчетов с различными сетками, шагами по времени, моделями выбора пары для столкновения, алгоритмами измельчения пространственных сеток и т.п., чтобы выбрать корректный набор технических параметров для конкретной постановки задачи. Платой за это преимущество является заметно большая вычислительная нагрузка метода EDMD по сравнению с методом DSMC. При этом EDMD существенно превосходит классическое MD-моделирование по скорости расчета, что позволяет использовать EDMD для решения задач динамики разреженного газа.

Подробное описание метода EDMD, его реализация и верификация на двух классических задачах динамики разреженного газа даны в предыдущих статьях авторов [20]. Также метод EDMD успешно применялся для моделирования течения многокомпонентной смеси газов в области со сложными неизотермическими [21] и движущимися [22–24] границами. Для случая многоатомных газов была предложена модель столкновений [25], которая была успешно протестирована на задаче об истечении азота в вакуум. Ниже описаны основные положения метода EDMD, а также специфические особенности его применения к задаче настоящего исследования.

2.1. Основные положения метода

Основная идея метода состоит в том, что траектории движения молекул газа между столкновениями можно считать прямолинейными. Таким образом, нет необходимости численно интегрировать уравнения движения с малым шагом по времени, как это делается при решении задачи методом классического молекулярно-динамического моделирования. Это предположение позволяет существенно снизить вычислительную нагрузку: характерный размер шага интегрирования при молекулярно-динамическом моделировании – ${{10}^{{ - 15}}}$ с, характерная скорость молекулы порядка тепловой – ${{10}^{3}}$ м/с, а характерная длина свободного пробега в газе при нормальных условиях порядка ${{10}^{{ - 7}}}$  м. Следовательно, даже при атмосферном давлении в среднем при классическом молекулярно-динамическом моделировании пришлось бы сделать около ${{10}^{5}}$ шагов для расчета прямого пролета молекулы между столкновениями. В подходе событийного моделирования эти пролеты предполагаются прямолинейными с постоянной скоростью, которую молекула приобрела при столкновении, и выполняются “за один шаг”. Таким образом, основная задача заключается в нахождении времен столкновений частиц друг с другом. Для этого строится упорядоченный список событий (очередь), которые должны произойти в системе. Событиями здесь называются акты межмолекулярных столкновений, столкновения с поверхностью твердого тела, столкновения с плоскостью симметрии задачи (если имеется), вылет или влет молекулы через открытые границы задачи. События из списка последовательно выполняются, при этом список предстоящих событий частично перестраивается. Например, при столкновении двух частиц меняются их скорости, поэтому часть предстоящих столкновений с другими частицами отменяется, а по новым скоростям считаются новые потенциальные столкновения. Из вышесказанного следует, что время в моделируемой системе меняется дискретно – по временам событий из упорядоченного списка. При этом нет нужды пересчитывать положения всех молекул в системе при каждом таком изменении времени, так как их движение полностью детерминировано, если были известны координата и скорость молекулы на определенный момент времени и список предстоящих событий.

Каждой молекуле приписаны (1) ее сорт, (2) координаты ${{{\mathbf{x}}}_{i}}$, (3) скорость ${{{\mathbf{u}}}_{i}}$, (4) энергия ${{e}_{{in,i}}}$ ее внутренних степеней свободы и (5) время $t$ последнего произошедшего с ней события. Между столкновениями молекулы движутся по прямым линиям. В настоящей работе используется модель столкновения [25], описанная в следующем разделе. В этом случае время столкновения двух молекул определяется из квадратного уравнения:

(2.1)
$\begin{array}{*{20}{c}} {{{{\left( {{{{\mathbf{x}}}_{1}}\left( t \right) - {{{\mathbf{x}}}_{2}}\left( t \right)} \right)}}^{2}} = b_{M}^{2},} \end{array}$
(2.2)
$\begin{array}{*{20}{c}} {{{{\mathbf{x}}}_{k}}\left( t \right) = {\mathbf{x}}_{k}^{'} + {\mathbf{u}}_{k}^{'}(t - t_{k}^{'}),~\quad k = 1,2~,} \end{array}$
где $t_{k}^{'}$ – время последнего события k-й молекулы, ${\mathbf{x}}_{k}^{'}$, ${\mathbf{u}}_{k}^{{\text{'}}}$ – ее координаты и скорость после последнего столкновения, ${{b}_{M}}$ – параметр модели столкновения (см. следующий раздел). Необходимо выбрать минимальный корень с условиями $t > t_{1}^{'}$, $t > t_{2}^{'}$. Если такого корня нет, то эти молекулы не столкнутся.

Аналогичным образом рассчитываются времена столкновений молекулы с поверхностями твердого тела и виртуальными границами (плоскости симметрии, периодические или открытые границы). Эти поверхности задаются набором алгебраических поверхностей с ограничениями $\left( {{{x}_{{{\text{min}}}}},~\;{{x}_{{{\text{max}}}}},~\;{{y}_{{{\text{min}}}}},~\;{{y}_{{{\text{max}}}}},~\;{{z}_{{{\text{min}}}}},~\;{{z}_{{{\text{ma}}x}}}} \right)$ на координаты. В настоящее время реализованы поверхности первого и второго порядка. В задаче настоящего исследования присутствуют плоские границы (входное и выходное сечение трубы), цилиндрическая поверхность (стенка трубы) и поверхность сферы (тело). Например, время столкновения i-й молекулы с k-м плоским элементом границы определяются из следующего линейного уравнения:

${{{\mathbf{x}}}_{i}}\left( t \right) \cdot {{{\mathbf{n}}}_{k}} = {{R}_{k}},$
где ${{{\mathbf{n}}}_{k}}$, ${{R}_{k}}~$ – нормаль и расстояние от начала координат, которые определяют плоскость k-го элемента границы. Аналогичным образом записываются квадратные уравнения для поверхностей второго порядка (цилиндр, сфера, конус).

Для реализации граничного условия на открытых границах задачи необходимо создавать молекулы, которые влетают в расчетную область. В терминах событийного метода это означает, что необходимо с некоторой частотой $\omega $ добавлять событие появления молекулы в определенном месте границы. Частота возникновения молекул $\omega $ и распределение их скоростей определяются давлением, температурой и макроскопической скоростью в резервуаре, из которого они попадают в расчетную область. Координата возникновения частицы выбирается случайным образом с равномерным распределением по площади данного элемента границы. В настоящей работе были реализованы равновесные граничные условия, т.е. когда расчетная область граничит с объемом газа, находящимся в термодинамическом равновесии с числовой плотностью ${{n}_{0}}$, средней скоростью ${{{\mathbf{v}}}_{w}}$ и температурой ${{T}_{0}}$. В этом случае поток молекул, влетающий в расчетную область, будет:

(2.3)
$\begin{gathered} {{J}_{1}} = {{n}_{0}}{{\left( {\frac{\beta }{{\sqrt \pi }}} \right)}^{3}}\mathop \smallint \limits_0^\infty \mathop \iint \limits_{ - \infty }^\infty \,{{u}_{n}}\exp \left( { - {{\beta }^{2}}{{{\left( {{\mathbf{u}} - {{{\mathbf{v}}}_{{\mathbf{w}}}}} \right)}}^{2}}} \right)d{\mathbf{u}} = \\ = \frac{{{{n}_{0}}}}{{2\sqrt \pi \beta }}\left( {{\text{exp}}\left( { - {{\beta }^{2}}\vartheta _{{wn}}^{2}} \right) + \sqrt \pi \beta {{\vartheta }_{{wn}}}\operatorname{erfc} \left( { - \beta {{\vartheta }_{{wn}}}} \right)} \right), \\ \end{gathered} $
где $\beta = \sqrt {m{\text{/}}2k{{T}_{0}}} $, ${\text{erfc}}\left( x \right)$ – дополнительная функция ошибок. Тогда частота создания молекулы будет $\omega = {{J}_{1}}S$, где $S$ – площадь рассматриваемого элемента открытой границы. В соответствии с равновесным распределением Максвелла плотность вероятности для скорости создаваемой молекулы будет

(2.4)
$\rho \left( {\mathbf{u}} \right) = \frac{{{{n}_{0}}}}{{{{J}_{1}}}}{{\left( {\frac{\beta }{{\sqrt \pi }}} \right)}^{3}}{{u}_{n}}\exp \left( { - {{\beta }^{2}}{{{\left( {{\mathbf{u}} - {{{\mathbf{v}}}_{w}}} \right)}}^{2}}} \right) = \frac{{2{{\beta }^{4}}{{u}_{n}}\exp \left( { - {{\beta }^{2}}{{{\left( {{\text{u}} - {{{\mathbf{v}}}_{w}}} \right)}}^{2}}} \right)}}{{\pi \left( {\exp \left( { - {{\beta }^{2}}\vartheta _{{wn}}^{2}} \right) + \sqrt \pi \beta {{\vartheta }_{{wn}}}\operatorname{erfc} \left( { - \beta {{\vartheta }_{{wn}}}} \right)} \right)}}.$

События, которые должны произойти с молекулами системы, размещаются в сортированном списке, который называется очередью событий. Поскольку с каждой молекулой связано как минимум одно предстоящее событие (даже если столкновений с другими молекулами нет, есть столкновение с границами задачи – поверхностями твердого тела или вылет через открытую границу), то размер этой очереди будет как минимум порядка количества частиц в системе $N$. Самая затратная операция, которая выполняется с очередью – поиск места для вставки очередного события. С ростом числа частиц в системе вычислительная сложность поиска в сортированном списке будет расти как логарифм (например, при реализации бинарного дерева). Тестовые запуски показали, что при такой реализации даже при количестве частиц 105 (что неприемлемо мало для задач динамики разреженного газа) счет становится очень затратным. Поэтому для реализации очереди событий был выбран метод с использованием хеш-функции или хеш-таблицы. При такой реализации вычислительная сложность поиска постоянна и не увеличивается с ростом количества элементов в списке (при адекватном подборе параметров). Идея метода состоит в следующем. В памяти выделяется место под ${{N}_{h}}$ элементов и выбирается параметр шага хеш-таблицы ${{\Delta }_{h}}$. Когда требуется расположить событие со временем $t$ в таблице, вычисляется хеш-функция:

$H\left( t \right) = \frac{t}{{{{\Delta }_{h}}}}\bmod {{N}_{h}}.$
Элемент добавляется по данному адресу. Ситуация, когда это место уже занято, называется коллизией хеш-функции. В этом случае элементы располагаются в обычный сортированный список, голова которого находится по данному адресу. Параметры хеш-таблицы ${{\Delta }_{h}}$ и ${{N}_{h}}$ выбираются так (см. [20]), чтобы в среднем хеш-таблица была равномерно заполнена и почти не было коллизий. Данный подход позволил проводить расчеты с нагрузкой по 107 частиц на ядро процессора.

2.2. Модель межмолекулярных столкновений

Актуальной проблемой в динамике разреженного газа является учет многоатомности молекул, так как в части задач это приводит к существенному отличию результатов от случая течения бесструктурных частиц. Особенно показательными являются задачи об истечении газа в вакуум – в экспериментах получено [26, 27], что зависимость поступательных и вращательных температур существенно отличается от адиабатических зависимостей. Поэтому для корректного описания течения в таких задачах требуется учитывать обмен энергией между поступательными и вращательными степенями свободы молекулы, в результате которого распределение энергии по степеням существенно отличается от равновесного. В работах [28, 29] сделан вывод, что необходимо описывать перераспределение энергии в терминах микропараметров: хотя бы в зависимости от поступательной, вращательной энергии до столкновения и прицельного расстояния. В связи с этим используемый авторами настоящего исследования метод событийного молекулярно-динамического моделирования (EDMD) идеально подходит для таких задач, так как в рамках событийного подхода перечисленные микропараметры полностью детерминированы для каждого столкновения (в отличие от DSMС, где часть из них выбрасывается случайным образом в предположении равновесного распределения).

2.2.1. Поступательно-вращательная релаксация. Столкновение двух молекул обычно характеризуется вектором скорости $g$ относительного движения центров масс молекул и энергиями ${{e}_{{in,1}}}$, ${{e}_{{in,2}}}$ внутренних степеней свободы первой и второй молекулы соответственно. В случае термодинамического равновесия можно получить вид плотности вероятности для модуля вектора относительной скорости $g$:

${{\rho }_{g}}\left( g \right) = 2\beta _{r}^{4}{{g}^{3}}{\text{exp}}\left( { - \beta _{r}^{2}{{g}^{2}}} \right),$
где ${{\beta }_{r}} = \sqrt {{{m}_{r}}{\text{/}}2kT} $, ${{m}_{r}} = {{m}_{1}}{{m}_{2}}{\text{/}}\left( {{{m}_{1}} + {{m}_{2}}} \right)$, а ${{m}_{1}}$, ${{m}_{2}}$ – массы молекул. Введем энергию относительного движения центров масс молекул до столкновения ${{e}_{r}} = {{m}_{r}}{{g}^{2}}{\text{/}}2$. Несложно получить выражение для плотности вероятности ${{e}_{r}}$:
${{\rho }_{{tr}}}\left( {{{e}_{r}}} \right) = {{e}_{r}}{{\left( {kT} \right)}^{{ - 2}}}{\text{exp}}\left( { - \frac{{{{e}_{r}}}}{{kT}}} \right)~.$
В случае если у молекулы две вращательных степени, плотность вероятности вращательной энергии ${{e}_{{in,i}}}$ одной молекулы будет
${{\rho }_{{in,i}}}\left( {{{e}_{{in,i}}}} \right) = {{\left( {kT} \right)}^{{ - 1}}}{\text{exp}}\left( { - \frac{{{{e}_{{in,i}}}}}{{kT}}} \right).$
Тогда плотность вероятности для общей энергии вращательных степеней двух молекул ${{e}_{{in}}} = {{e}_{{in,1}}} + {{e}_{{in,2}}}$ будет

${{\rho }_{{in}}}\left( {{{e}_{{in}}}} \right) = {{\left( {kT} \right)}^{{ - 2}}}{{e}_{{in}}}\exp \left( { - \frac{{{{e}_{{in}}}}}{{kT}}} \right).$

Как видно, сформулировать корректную микроскопическую модель релаксации энергии при столкновении молекул в терминах ${{e}_{{in,1}}}$, ${{e}_{{in,2}}}$ и ${{e}_{r}}$ не представляется возможным, так как равновесные распределения данных параметров зависят от температуры, а пара сталкивающихся молекул “не знает” температуру в данной точке пространства (как и другие моменты функции распределения скоростей). Поэтому авторами предложено (см. [25]) перейти к другой тройке параметров:

(1) $\varepsilon = {{e}_{r}} + {{e}_{{in}}}$ – полная энергия движения двух молекул относительно центра масс системы, которая, как следует из законов сохранения, не меняется при столкновении $\varepsilon = \varepsilon {\kern 1pt} '$ (штрихом будем обозначать значения соответствующих параметров после столкновения);

(2) безразмерная величина ${{\gamma }_{r}} = {{e}_{r}}{\text{/}}\varepsilon $, характеризующая долю энергии относительного движения центров масс молекул;

(3) безразмерная величина ${{\gamma }_{1}} = {{e}_{{in,1}}}{\text{/}}{{e}_{{in}}}$, характеризующая долю энергии вращательных степеней первой молекулы в энергии вращательных степеней пары.

Равновесные плотности вероятности для введенных величин $\varepsilon $, ${{\gamma }_{r}}$ и ${{\gamma }_{1}}$ получаются очень удобными для конструирования модели релаксации:

$\begin{gathered} {{\rho }_{\varepsilon }}\left( \varepsilon \right) = \frac{1}{{6{{{\left( {kT} \right)}}^{4}}}}{{\varepsilon }^{3}}{\text{exp}}\left( { - \frac{\varepsilon }{{kT}}} \right), \\ {{\rho }_{r}}\left( {{{\gamma }_{r}}} \right) = 6{{\gamma }_{r}}\left( {1 - {{\gamma }_{r}}} \right),\quad {{\rho }_{1}}\left( {{{\gamma }_{1}}} \right) = 1. \\ \end{gathered} $
Как видно, равновесные плотности вероятности ${{\rho }_{r}}\left( {{{\gamma }_{r}}} \right)$ и ${{\rho }_{1}}\left( {{{\gamma }_{1}}} \right)$ не зависит от температуры $T$ и других моментов функции распределения скоростей, а ${{\rho }_{\varepsilon }}\left( \varepsilon \right)$ – зависит, но это не важно, так как величина $\varepsilon $ является инвариантом столкновения. Таким образом, в предложенных терминах появляется возможность корректно сформулировать микроскопическую модель релаксации энергии при межмолекулярном столкновении. Что и было сделано (см. [25]) на основе анализа траекторий столкновения молекул азота и кислорода. В результате получена модель межмолекулярных столкновений, которая описана ниже.

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

Предлагается следующая модель релаксации для расчета изменения ${{\gamma }_{r}}$ и ${{\gamma }_{1}}$ при столкновении:

$\begin{gathered} R(\varepsilon ,b,{{\gamma }_{r}},{{\gamma }_{1}},\gamma _{r}^{'},\gamma _{1}^{'}) = \left( {\left( {1 - {{\alpha }_{r}}\left( {\varepsilon ,b} \right)} \right)\delta (\gamma _{r}^{'} - {{\gamma }_{r}}) + {{\alpha }_{r}}\left( {\varepsilon ,b} \right){{\rho }_{r}}(\gamma _{r}^{'})} \right) \times \\ \times \;\left( {\left( {1 - {{\alpha }_{1}}\left( {\varepsilon ,b} \right)} \right)\delta (\gamma _{1}^{'} - {{\gamma }_{1}}) + {{\alpha }_{1}}\left( {\varepsilon ,b} \right){{\rho }_{1}}(\gamma _{1}^{'})} \right), \\ \end{gathered} $
где $b$ – прицельное расстояние, ${{\alpha }_{r}}\left( {\varepsilon ,b} \right)$ и ${{\alpha }_{1}}\left( {\varepsilon ,b} \right)$ – коэффициенты релаксации для ${{\gamma }_{r}}$ и ${{\gamma }_{1}}$ соответственно, зависящие от $\varepsilon $ и $b$. Легко показать аналитически, что $R(\varepsilon ,b,{{\gamma }_{r}},{{\gamma }_{1}},\gamma _{r}^{'},\gamma _{1}^{'})$ удовлетворяет условиям детального баланса и нормировки для любых $\varepsilon $ и $b$. Вид функций ${{\alpha }_{r}}\left( {\varepsilon ,b} \right)$ и ${{\alpha }_{1}}\left( {\varepsilon ,b} \right)$ был выбран, исходя из полученных в МД-расчетах результатов [25]:
$\alpha \left( {\varepsilon ,b} \right) = {{\alpha }_{0}}q\left( \varepsilon \right){{w}_{0}}\left( b \right) + {{\alpha }_{\infty }}\left( {1 - q\left( \varepsilon \right)} \right){{w}_{\infty }}\left( b \right),$
где $q\left( \varepsilon \right) = {\text{exp}}\left( { - \varepsilon {\text{/}}{{\varepsilon }_{0}}} \right)$, ${{w}_{0}}\left( b \right) = \sqrt[n]{{1 - {{{\left( {b{\text{/}}{{b}_{0}}} \right)}}^{n}}}}$, ${{w}_{\infty }}\left( b \right) = {\text{exp}}\left( { - {{{\left( {b{\text{/}}{{b}_{M}}} \right)}}^{m}}} \right)$, а величины ${{\alpha }_{0}}$, ${{\alpha }_{\infty }}$, ${{\varepsilon }_{0}}$, ${{b}_{0}}$, ${{b}_{M}}$, $n$, $m$ – набор констант (отличающийся для ${{\gamma }_{r}}$ и ${{\gamma }_{1}}$), значения которых выбирались из условия минимизации невязки с результатами МД расчетов.

Угол отклонения $\chi $ вектора относительной скорости ${\mathbf{g}}$ был аппроксимирован следующим образом:

$\chi \left( {\varepsilon ,b,{{\gamma }_{r}}} \right) = \left\{ \begin{gathered} \pi \left( {1 + {{k}_{1}}\left( {\varepsilon ,{{\gamma }_{r}}} \right){{b}^{{{{n}_{1}}}}} + {{k}_{2}}\left( {\varepsilon ,{{\gamma }_{r}}} \right){{b}^{{{{n}_{2}}}}}} \right),~\quad b \leqslant {{b}_{1}}\left( {{{e}_{r}}} \right), \hfill \\ {{P}_{3}}\left( b \right),~\quad {{b}_{1}}\left( {{{e}_{r}}} \right) < b \leqslant {{b}_{M}}, \hfill \\ \end{gathered} \right.$
где ${{b}_{1}}\left( {{{e}_{r}}} \right) = \left( {{{b}_{M}} - {{b}_{m}}} \right){\text{exp}}\left( { - {{e}_{r}}{\text{/}}{{e}_{0}}} \right) + {{b}_{m}}$, ${{k}_{i}}\left( {\varepsilon ,{{\gamma }_{r}}} \right) = {{k}_{{i,0}}} + {{k}_{{i,1}}}{{\gamma }_{r}} + {{k}_{{i,2}}}\varepsilon $, а ${{P}_{3}}\left( b \right)$ – полином 3 степени, коэффициенты которого однозначно определяются из следующих условий: (1) непрерывности $\chi \left( {\varepsilon ,b,{{\gamma }_{r}}} \right)$ в точке $b = {{b}_{1}}\left( {{{e}_{r}}} \right)$, (2) гладкости $\chi \left( {\varepsilon ,b,{{\gamma }_{r}}} \right)$ в точке $b = {{b}_{1}}\left( {{{e}_{r}}} \right)$, (3) $\chi \left( {\varepsilon ,{{b}_{M}},{{\gamma }_{r}}} \right) = 0$, (4) $\partial \chi {\text{/}}\partial b = 0$ в точке $b = {{b}_{M}}$. Величины ${{b}_{m}}$, ${{e}_{0}}$, ${{n}_{1}}$, ${{n}_{2}}$, ${{k}_{{i,j}}}$ – набор констант, значения которых выбирались из условия минимизации невязки с результатами МД расчетов.

2.2.3. Программная реализация модели.Предложенная микроскопическая модель столкновений не усложняет алгоритм моделирования течения разреженного газа методами EDMD и D-SMC. Наоборот, в предложенной модели упрощается процедура оценки количества необходимых столкновений по сравнению с популярными у исследователей моделями переменных твердых сфер (VHS) и переменных мягких сфер (VSS), так как эффективный размер частицы не зависит от относительной скорости молекул, а полагается равным ${{b}_{M}}$. Поэтому алгоритм в этой части совпадает с наиболее простой моделью твердых сфер (HS). Когда пара молекул для столкновения установлена, алгоритм обработки столкновения выглядит следующим образом:

(1) рассчитываются безразмерные величины ${{\gamma }_{r}} = {{e}_{r}}{\text{/}}\varepsilon $ и ${{\gamma }_{1}} = {{e}_{{in,1}}}{\text{/}}{{e}_{{in}}}$ (это значения до столкновения);

(2) в соответствии с описанным выше ядром $R(\varepsilon ,b,{{\gamma }_{r}},{{\gamma }_{1}},\gamma _{r}^{'},\gamma _{1}^{'})$ разыгрываются значения $\gamma _{r}^{'}$, $\gamma _{1}^{'}$ (значения после столкновения);

(3) рассчитывается угол отклонения $\chi $ вектора относительной скорости по введенному выше выражению $\chi \left( {\varepsilon ,b,{{\gamma }_{r}}} \right)$.

(4) так как полная энергия $\varepsilon $ молекул относительно центра масс системы в ходе столкновения не меняется, то $\varepsilon = \varepsilon {\kern 1pt} '$. По параметрам $(\varepsilon ,\gamma _{r}^{'},\gamma _{1}^{'},\chi )$ однозначно рассчитываются скорости центров масс и вращательные энергии молекул после столкновения.

Авторами настоящей статьи подготовлена и протестирована программная C++ реализация микроскопической модели столкновений (выложена в открытом доступе [30]), которую можно использовать в DSMC и EDMD расчетах. Несмотря на то, что с алгоритмической точки зрения существенных усложнений нет, вычислительные затраты при использовании новой модели возрастают. Это происходит из-за увеличения эффективного размера частицы (он полагается равным ${{b}_{M}}$), что приводит к заметному увеличению количества обрабатываемых столкновений.

2.3. Распараллеливание расчетов по областям

Предварительные расчеты в плоской постановке показали, что случай сверхзвукового движения значительно проще в вычислительном плане, так как ударная волна отходит от тела совсем на небольшое расстояние (1–2 калибра трубы), что позволяет существенно сократить размер расчетной области перед телом. В дозвуковом случае, напротив, возмущения от тела могут быть ощутимы на расстоянии в 5–7 и более калибров трубы перед телом, поэтому такие расчеты требуют значительного увеличения расчетной области. При этом с точки зрения перспектив высокоскоростного вакуумного транспорта наибольший практический интерес представляет трехмерная постановка именно в дозвуковом режиме. Вопрос существенного увеличения расчетной области стоит особенно остро в трехмерной постановке. Поэтому для проведения такого расчета потребовалось распараллелить метод событийного молекулярно-динамического моделирования (EDMD) по областям, что заслуживает отдельного внимания из-за известных проблем с распараллеливанием методов такого типа [31].

В настоящий момент не разработан алгоритм распараллеливания событийного подхода “без потерь” (т.е. без потери событий), который бы при этом оставался эффективным с вычислительной точки зрения [31]. Объясняется это следующим образом. Программная реализация подхода событийного молекулярно-динамического моделирования является бессеточной (нет сетки ни по пространству, ни шага по времени). Время в моделируемой системе меняется по временам происходящих событий. Поэтому введение искусственного временного интервала, на котором расчет на каждом ядре идет независимо, неотвратимо ведет к накоплению событий с молекулами, которые должны были бы произойти на границах областей разных ядер, но остались не обработаны. Если придерживаться консервативного подхода “без потерь”, то после интервала параллельного счета необходимо откатиться назад по времени до первого такого необработанного “граничного события”, обработать его и запустить следующий интервал независимого счета. В итоге, из-за необходимости отката получить ускорение счета не удается. Поэтому приходится выбирать менее консервативный подход с потерей событий. Легко сделать оценку доли таких событий – она определятся долей молекул газа, находящихся в слое границы области, по отношению ко всем молекулам области. Толщина этого слоя пропорциональна длине интервала независимого счета, т.е. легко регулируется. Следует оговориться, что даже при потере этой небольшой части событий расчетная схема всегда остается строго консервативной по массе, импульсу и энергии, поэтому с точки зрения макроскопических параметров результаты для одноядерного расчета и расчета с распараллеливанием получаются идентичными, а различия возникают только в индивидуальных траекториях чрезвычайно малой доли молекул. Данный подход был реализован (технология MPI, в синхронном режиме) и применялся для трехмерного расчета течения в данной работе. Разбивка на области была очень простой – труба нарезалась на куски ее поперечными сечениями (см. фиг. 2). Отношение времени расчета одной и той же задачи на одном ядре и на 5 ядрах (тесты проводились на процессоре Intel Core i7–9700) составило ${{\Delta }_{1}}{\text{/}}{{\Delta }_{5}} = 3.55$, т.е. эффективность текущей параллельной реализации даже на таком небольшом количестве ядер порядка 0.7. Это объясняется тем, что балансировка нагрузки между ядрами осуществлялась в ручном режиме. Поэтому в период времени, когда течение устанавливается, часто возникала ситуация, что часть ядер простаивала, ожидая завершения интервала расчета на остальных. Чтобы повысить эффективность распараллеливания, авторы планируют реализовать автоматическое изменение размеров подобластей в процессе расчета, так как вычислительная нагрузка ядра напрямую зависит от количества молекул в конкретной подобласти, которое заранее не известно и меняется по мере установления течения.

Фиг. 2.

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

3. РАСЧЕТЫ И РЕЗУЛЬТАТЫ

3.1. Проведение расчетов

В качестве тела рассматривалась сфера диаметра $d$, которая движется в трубе с постоянной скоростью ${{{v}}_{w}}$. Температура всех твердых поверхностей бралась постоянной и равной ${{T}_{w}} = {{T}_{0}} = 300$ K. Расчетная область представляла из себя половину трубы (см. фиг. 2), отрезанную по плоскости симметрии Oxy. Длина расчетной области по оси Ox бралась равной $L = 20D$ (20 калибров трубы). Трехмерная постановка задачи очень ресурсоемкая, поэтому вместо полноценного параметрического исследования были рассмотрены конкретные комбинации параметров, которые максимально приближены к требованиям для различных концепций вакуумного транспорта, известных из литературы.

• Рассмотрены следующие относительные размеры тела: $d{\text{/}}D = 0.2$, 0.5, 0.8. Наибольший практический интерес представляют высокие значения $d{\text{/}}D$.

• В качестве газовой среды взят воздух: смесь азота и кислорода с использованием модели столкновений, предложенной авторами в предыдущих работах [25] для учета вращательных степеней свободы в двухатомных газах.

• Целевая скорость для таких транспортных систем порядка 1000 км/час, в случае воздуха и температур около 300 К это соответствует числу Маха ${\text{M}} = 0.8$.

• C точки зрения практических приложений результаты для чисел Кнудсена порядка 1 и более бесполезны, так как при реализации потребуют слишком глубокой откачки трубы. Поэтому расчеты проводились для ${\text{Kn}}$ порядка 0.1 (параметр разрежения $\delta = 10$).

В начальный момент времени труба заполнялась газом с параметрами невозмущенного потока. Далее проводился EDMD расчет до выхода макроскопических параметров на стационарные значения. Выход решения на стационар происходил за время порядка $50D{\text{/}}a$, где $a$ – скорость звука в смеси при $T = {{T}_{0}}$. При этом массовый расход во всех поперечных сечениях расчетной области выравнивался, а суммарное касательное напряжение, действующее на трубу, становилось равно по модулю силе сопротивления тела ${{R}_{x}}$. После установления течения вычислялись стационарные поля основных макропараметров и потоки к твердым границам задачи, значения которых представлены ниже.

3.2. Результаты расчетов

Сначала был рассмотрен случай движения тела с нулевым смещением ${{y}_{0}} = 0$ тела относительно оси трубы. В этом случае задача осесимметрична, но расчеты все равно проводились трехмерным кодом в расчетной области, описанной выше. Было исследовано обтекание тела трех относительных размеров $d{\text{/}}D = 0.2$, 0.5, 0.8 при скорости движения тела ${\text{M}} = 0.8$. Полученные коэффициенты сопротивления приведены в табл. 2.

Таблица 2.  

Коэффициенты сопротивления ${{C}_{x}}$ при ${\text{M}} = 0.8$ для различных размеров тела

$d{\text{/}}D$ 0.2 0.5 0.8
${{C}_{x}}$ 5.1 5.3 10.8

До проведения расчетов у авторов настоящей статьи были опасения, что при заметных размерах тела по отношению к поперечному размеру трубы (например, $d{\text{/}}D = 0.8$) и высоких дозвуковых скоростях (например, М = 0.8) не будет существовать стационарного решения, а будет реализовываться нестационарное течение, напоминающее эффект запирания аэродинамической трубы. Однако проведенные расчеты показали, что стационарное решение существует, соответствующие ему стационарные поля плотности, скорости и давления для него представлены на фиг. 3. Значимым отличием рассматриваемой задачи от аэродинамической трубы является наличие движущихся стенок (в системе координат, связанной с телом, труба набегает вместе с потоком) и вызванных ими касательных напряжений в газе (см. фиг. 4). Эти касательные напряжения работают наподобие насоса, плавно повышая давление перед телом в 5–6 раз (см. фиг. 3). На расстояниях же в 6–7 калибров трубы от тела вверх по потоку течение оказывается почти не возмущено.

Фиг. 3.

Поле плотности, скорости и давления при $d{\text{/}}D = 0.8$ и ${{y}_{0}}{\text{/}}D = 0$ (осесимметричный случай).

Фиг. 4.

Поле касательных напряжений ${{p}_{{xy}}}$ при $d{\text{/}}D = 0.8$ и ${{y}_{0}}{\text{/}}D = 0$ (осесимметричный случай).

Далее были проведены расчеты для размера тела $d{\text{/}}D = 0.8$ и нескольких различных смещений ${{y}_{0}}{\text{/}}D = 0.01$, 0.02, 0.05, 0.09. В этом случае задача перестает быть осесимметричной и появляется подъемная сила ${{R}_{y}}$. Полученные коэффициенты сопротивления ${{C}_{x}}$ и подъемной силы ${{C}_{y}}$ представлены в табл. 3 и на фиг. 5 в зависимости от относительного смещения ${{y}_{0}}{\text{/}}D$. Подъемная сила направлена к оси трубы, т.е. центр поперечного сечения трубы является устойчивым положением равновесия (с точки зрения поперечной силы, действующей на тело со стороны газа, в реальности же есть еще сила тяжести). Характерные для несимметричного случая поля плотности, скорости и давления представлены на фиг. 6. Рассмотрим подробнее поля скорости, соответствующие симметричному и несимметричному случаю (см. фиг. 7). Видно, что перед телом возникает обратная струя, которая, взаимодействуя с набегающим потоком, образует кольцевую струю, обтекающую тело. Еще одной интересной особенностью данного течения является наличие в потоке перед телом контура (на фиг. 7 нанесен штриховой линией), на котором модуль скорости практически нулевой. То есть, если бы этот контур был продолжением тела, то поле течения вне контура практически не изменилось. Более того, несложно получить из интегрального закона сохранения количества движения в стационарном случае, что и сопротивление такого тела равнялось бы сопротивлению тела исходной задачи. Таким образом, основываясь на результатах расчета данной задачи для сферического тела, можно построить более сложную форму тела, которая будет практически эквивалентна сфере с точки зрения аэродинамического сопротивления.

Таблица 3.  

Коэффициент сопротивления ${{C}_{x}}$ и подъемной силы ${{C}_{y}}$ при $d{\text{/}}D = 0.8$ и различных смещениях ${{y}_{0}}{\text{/}}D$ тела относительно оси трубы

${{y}_{0}}{\text{/}}D$ 0 0.01 0.02 0.05 0.09
${{C}_{x}}$ 10.8 10.8 10.8 11.0 12.1
${{C}_{y}}$ 0 0.15 0.31 0.86 2.64
Фиг. 5.

Коэффициент сопротивления и подъемной силы в зависимости от положения тела в сечении трубы.

Фиг. 6.

Поле плотности, скорости и давления при $d{\text{/}}D = 0.8$ и ${{y}_{0}}{\text{/}}D = 0.05$.

Фиг. 7.

Поле скорости перед телом при $d{\text{/}}D = 0.8$ и двух различных смещений ${{y}_{0}}{\text{/}}D = 0$ и ${{y}_{0}}{\text{/}}D = 0.05$. Штриховой линией обозначен контур, на котором модуль скорости практически нулевой.

ЗАКЛЮЧЕНИЕ

Методом событийного молекулярно-динамического моделирования (EDMD) исследовано течение разреженного газа, которое возникает в трубе при движении по ней тела с высокой дозвуковой скоростью. Рассмотрен случай, когда давление в трубе соответствует числам Кнудсена порядка 0.1 (параметр разрежения $\delta = 10$). Расчеты проводились в условиях, максимально приближенных к целевым параметрам концепций высокоскоростного вакуумного транспорта, которые известны из литературы.

• В качестве газа взят воздух (смесь азота и кислорода) с учетом вращательных степеней свободы в двухатомных молекулах.

• Скорость движения тела около 1000 км/ч, что соответствует числу Маха ${\text{M}} = 0.8$ в воздухе при температуре 300 K.

• Рассмотрены поперечные размеры тела вплоть до $d{\text{/}}D = 0.8$, т.е. перекрывающие 64% площади сечения трубы.

Получено, что даже для больших относительных поперечных размеров тела $d{\text{/}}D$ существует стационарное решение. Характерными особенностями этого решения является (1) плавное нарастание давления перед телом за счет значительных касательных напряжений в газе, вызванных движущимися стенками трубы (в системе координат, связанной с телом, стенки трубы набегают вместе с потоком), (2) образование перед телом обратной струи, взаимодействующей с набегающим потоком, за счет чего перед телом образуется контур, на котором модуль скорости практически нулевой. В результате упомянутых касательных напряжений, которые работают наподобие компрессора, и увеличения давления перед телом (в 5–6 раз для $d{\text{/}}D = 0.8$), расчеты предсказывают высокие коэффициенты сопротивления: ${{C}_{x}} = 5.1$, 5.3, 10.8 для поперечных размеров тела $d{\text{/}}D = 0.2$, 0.5, 0.8 соответственно. В случае, если тело движется по прямой, смещенной относительно оси трубы, возникает подъемная сила, направленная к центру поперечного сечения трубы.

Список литературы

  1. Вейнберг Б.П. Движение без трения. СПб.: Естествоиспытатель, 1914.

  2. Musk E. Hyperloop Alpha Documents; SapceX: Hawthorne, CA, USA, 2013.

  3. Hyperloop One Website, https://hyperloop-one.com/.

  4. Hyperloop Transportation Technologies Website, https://www.hyperlooptt.com/.

  5. Kim H., Oh S. Shape optimization of a hyperloop pod’s head and tail using a multi-resolution morphing method// Int. J. Mech. Sci. 2022. V. 223. P. 107227.

  6. Sui Y. et al. An aerothermal study of influence of blockage ratio on a supersonic tube train system// J. Therm. Sci. 2020. V. 1–12.

  7. Chen X. et al. Aerodynamic simulation of evacuated tube maglev trains with different streamlined designs, J. Modern Transp. 2012. V. 20. № 2. P. 115–120.

  8. Zhou P., Zhang J., Li T. Effects of blocking ratio and Mach number on aerodynamic characteristics of the evacuated tube train// Int. J. Rail Transport. 2020. V. 8. № 1. P. 27–44.

  9. Oh J.-S. et al. Numerical Analysis of Aerodynamic Characteristics of Hyperloop System// Energies .2019.V. 12. № 3. P. 518.

  10. Sui Y., Niu J., Ricco P., Yuan Y., Yu Q., Cao X., Yang X. Impact of vacuum degree on the aerodynamics of a high-speed train capsule running in a tube // Int. J. Heat Fluid Flow. 2022. V. 88. P. 108752.

  11. Lluesma-Rodríguez F., González T., Hoyas S. CFD simulation of a hyperloop capsule inside a closed environment// Results Engng. 2021. V. 9. P. 100196.

  12. Le T.T.G. et al. Numerical investigation of aerodynamic drag and pressure waves in hyperloop systems// Mathematics. 2020. V. 8. № 11.

  13. Deng Z., Zhang W., Zheng J., Wang B., Ren Y., Zheng X., Zhang J. A high-temperature superconducting maglev-evacuated tube transport (HTS Maglev-ETT) test system // IEEE Trans. Appl. Supercond. 2017. V. 27. № 6. P. 1–8.

  14. Hruschka R., Klatt D. In-pipe aerodynamic characteristics of a projectile in comparison with free flight for transonic Mach numbers, Shock Waves. 2019. V. 29. № 2. P. 297–306.

  15. Seo Y., Cho M., Kim D.H., Lee T., Ryu J., Lee C. Experimental analysis of aerodynamic characteristics in the Hyperloop system, Aerosp. Sci. Technol., 2023. V. 137. P. 108265.

  16. Donev A., Garcia A.L., Alder B.J. Stochastic Event-Driven Molecular Dynamics // J. Comput. Phys. 2008. V. 227. P. 2644–2665.

  17. Valentini P., Schwartzentruber T.E. A combined Event-Driven/Time-Driven molecular dynamics algorithm for the simulation of shock waves in rarefied gases // J. Comput. Phys. 2009. V. 228. P. 8766–8778.

  18. Bannerman M.N., Sargant R., Lue L., Dynam O. A free O(N) general event-driven molecular-dynamics simulator // J. Comput. Chem. 2011. V. 32. P. 3329–3338.

  19. Akkaya V.R., Kandemir I. Event-driven molecular dynamics simulation of hard-sphere gas flows in microchannels // Math. Probl. Eng. 2015. № 2015.

  20. Yakunchikov A., Kosyanchuk V. Application of event-driven molecular dynamics approach to rarefied gas dynamics problems // Comput. Fluids. 2018. V. 170. P. 121–127.

  21. Yakunchikov A., Kosyanchuk V. Numerical investigation of gas separation in the system of filaments with different temperatures // Int. J. Heat Mass Transf. V. 138. P. 144–151.

  22. Yakunchikov A., Kosyanchuk V. A new principle of separation of gas mixtures in non-stationary transitional flows // Acta Astronaut. 2019.

  23. Artem Yakunchikov. The outflow of gas mixture into vacuum, periodically interrupted by bodies moving towards the jet // Vacuum Volume 209, March 2023, 111778.

  24. Yakunchikov Artem. Heat transfer in a rarefied gas between profiled surfaces moving relative to each other // International Journal of Heat and Mass Transfer Volume 184, March 2022. V. 122339.

  25. Yakunchikov A., Kosyanchuk V., Iuldasheva A. Rotational relaxation model for nitrogen and its application in free jet expansion problem, Phys. Fluids. 2020. V. 32. P. 102006.

  26. Marrone P.V. Temperature and Density Measurements in Free Jets and Shock Waves // Phys. Fluids. 1967. V. 10. P. 521.

  27. Mori H., Niimi T., Akiyama I., Tsuzuki T. Experimental detection of rotational non-Boltzmann distribution in supersonic free molecular nitrogen flows // Phys. Fluids. 2005. V. 17. P. 117103.

  28. Valentini P., Zhang C. and Schwartzentruber T.E. Molecular dynamics simulation of rotational relaxation in nitrogen: Implications for rotational collision number models // Phys. Fluids. 2012. V. 24. P. 106101.

  29. Tokumasu T. and Matsumoto Y. Dynamic molecular collision (DMC) model for rarefied gas flow simulations by the DSMC method // Phys. Fluids. 1999. V. 11. P. 1907–20.

  30. Реализация микроскопической модели столкновений для воздуха, https://multiscale.ru/science/collisionmodel.

  31. Miller S., Luding S. Event-driven molecular dynamics in parallel // J Comput Phys 2004. V. 193. P. 306–16.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал