АНАЛИЗ ПРИМЕНИМОСТИ МОДЕЛИ РЕШЁТОЧНОГО ГАЗА D2Q4 ДЛЯ РЕШЕНИЯ ЗАДАЧИ ИНЖЕКЦИИ Текст научной статьи по специальности «Физика»

This work is devoted to the study of the stability of an approximate solution to the problem of restoring the parameters of a mathematical model (continuous dependence of the solution on the initial data), which is found by the Tikhonov regularization method. The purpose of the work is to draw conclusions about the stability of the approximate solution (regularized solution) of the inverse problem with a small change in the measured initial data. To do this, we first consider a mathematical model described by a system of ordinary differential equations. Then, on the basis of the studied mathematical model, the inverse problem of restoring its parameters from the measured initial data was posed. From the inverse problem, the "exact" and "approximate" systems of algebraic equations with respect to the desired parameters will be compiled. Further, the "approximate" system of algebraic equations is solved by Tikhonov's regularization method. As a result, a regularized solution of an "approximate " system of algebraic equations is found, which is an approximation to the desired exact solution of the inverse problem. Finally, the stability of the found regularized solution is investigated. Our work helps to strengthen the theory in the process of constructing a method for solving inverse problems using Tikhonov's regularization method. The results of this work show that the regularized solution is completely stable with a small change in the initial data, i.e., it continuously depends on the initial data. Moreover, this solution is unique and is an approximation to the desired exact solution. Therefore, in practical problems, it can be used in the process of studying behavior, the properties of a mathematical model (solution of a direct problem).

Key words: stability, regularized solution, Tikhonov's regularization method, inverse problem, mathematical model, forecasting.

Le Van Huyen, postgraduate, huyenlevan120193@gmail.com, Russia, St. Petersburg, Peter the Great St. Petersburg Polytechnic University,

Chernenkaya Liudmila Vasilievna, doctor of technical science, professor, Russia, St. Petersburg, Peter the Great St. Petersburg Polytechnic University

УДК 623.463.3

DOI: 10.24412/2071-6168-2023-1-14-22

АНАЛИЗ ПРИМЕНИМОСТИ МОДЕЛИ РЕШЁТОЧНОГО ГАЗА D2Q4 ДЛЯ РЕШЕНИЯ ЗАДАЧИ ИНЖЕКЦИИ

Д.Б. Дмитриенко, А.М. Башанин

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

Ключевые слова: клеточный автомат, решеточный газ, математическое моделирование, ин-жекция, Ansys Fluent, перпендикулярные потоки.

Перпендикулярное столкновение двух потоков жидкости или газа - достаточно распространенная газодинамическая задача, решаемая в различных областях науки техники. Так, например инжек-ция струи газа или топлива в сверхзвуковую часть сопла ракетного двигателя является одним из способов управления его тягой и другими динамическими характеристиками [1]. Также, инжекция, смешение и горение топлива имеют место в рабочем тракте гиперзвуковых прямоточных воздушно-реактивных двигателей, а рассматриваемые процессы изучаются и моделируются на соответствующих установках -поршневых или мембранных гиперзвуковых ударных трубах [2,3].

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

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

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

14

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

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

В данной работе рассматривается модель HPP, представленная впервые в статье [5] и названная по инициалам её авторов - Hardy J., de Pazzis O., Pomeau Y. Она имеет 4 разрешенных направления движения частиц и достаточно успешно применяется, как простая модель для описания движения газов и жидкостей [4, 6].

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

Рис. 1. Модель решеточного газа с ортогональной структурой

Следует пояснить, что для обозначения решеток используется формула типа DnQш, где п -пространственная размерность решетки (п = 1 для одномерной модели, п = 2 для двумерной и т.д.), т -количество допустимых направлений движения частиц, то есть векторов их скоростей. Так например, как уже было приведено выше, решетка изображенная на рис. 1 обозначается как D2Q4. Модель с треугольной структурой имеет 6 возможных направлений и обозначается D2Q6, а для LBM кроме 8 направлений движения частицы, принято еще и состояние покоя, когда V; = (0,0), и соответственно формула решетки - D2Q9.

Частицы газа перемещаются между узлами ортогональной решетки со скоростями V;. Значение скоростей кратно размеру ячейки, что позволяет частицам за один шаг по времени йЬ попасть из одного узла решетки в другой. В простейшем случае вектор V; представляет собой набор единичных скоростей (рис.2):

ух = (1,0); у2 = (0,1); У3 = (-1,0); У4 = (0,-1);

Рис. 2. Разрешенные скорости в HPP

Таким образом новое положение частицы в решетке определяется как:

г.п = г.п-1+у.п-1^м, (1)

где г™ — вектор положения частицы I на текущем шаге по времени (п), ггп_1,угп_1 - вектор положения и скорость, соответственно, на предыдущем шаге по времени (п — 1).

Скорость частицы на текущем и предыдущем шаге не меняется, если не произошло столкновения с другой частицей или со стенкой (рис.3).

Правила описаны всего для трех типов столкновений: с другой частицей при разнонаправленных скоростях (лоб в лоб), с другой частицей при перпендикулярно направленных скоростях и столкновение со стенкой или твердой границей(рис.З).

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

п п+1

Рис. 3. Локальные правила для ортогональной сетки (НРР)

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

Однако следует сказать и о недостатках такой симуляции решеточного газа. Модель часто дает ошибку как раз из-за полного сохранения импульса по вертикальной и горизонтальной оси. Никакая энергия не удаляется из модели ни в результате столкновений, ни в результате движения, поэтому она будет продолжаться бесконечно. Также в модели отсутствует вращательная инвариантность, то есть возможность поворота оси под любым углом, без изменения её свойств, что приводит к значительной анизотропии. Как результат - все неортогональные и несимметричные структуры модели (вихри, размывания и т.п.) имеют квадратную форму [7].

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

П П+1

Рис. 4. Недостатки модели НРР

При определении столкновений, главное условие проверки - это отсутствие частиц с одинаковыми координатами в сетке. То есть в одном узле не должно быть более одной частицы.

На начальном этапе происходит заполнение области частицами с необходимой частотой и структурой (можно задать случайное заполнение, или имитировать слои).

16

Состояние ячейки при этом характеризуется числом заполнения а^, которое является функцией от радиус-вектора положения частицы Г; и текущего шага по времени п, т.е. а^ (г¿,п).

С учетом того, что область моделирования двумерная и дискретная, то для задания наличия/отсутствия частицы достаточно матрицы X на Y, строки и столбцы которой будут содержать логические значения: либо 0 либо 1. X и Y принадлежат множеству целых неотрицательных чисел. Также при заполнении матрицы создается массив скоростей частиц V;, в котором направление скоростей выбирается либо случайно - для хаотического движения газа, либо по выбору - для однонаправленного течения.

Таким образом, наличие частицы на п-ом шаге по времени в узле по радиус-вектору Г; со скоростью V; обозначается как:

а^п) = 1.

Если частица отсутствует в узле с координатами Г; на данном шаге по времени п, то:

а1(г1,п) = 0.

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

Это позволяет хранить результаты вычислений в более сжатом формате и уменьшить количество обрабатываемой информации при расчетах.

При описании эволюции клеточного автомата модели НРР взаимодействие частиц разбивают на две части:

- фазу столкновения, которая описывает число заполнения частиц, входящих в узел - а.

- фазу распространения, которая описывает число заполнения частиц, выходящих из узла, с учетом прописанных локальных правил - а°и£

Именно между этими двумя временными шагами происходит выполнение логики модели за счет работы функции столкновения - П(а1П) = П(г, V).

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

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

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

аош = а1п + п(а*п), (2)

или

ят —ят —атат+2(1 — ат+1)(1 — ат+3) + ят+1ят+з(1 — ят)(1 — ят+2), (3)

где т...т + 3 - индексы возможных направлений скоростей частиц, которых согласно НРР всего 4, при этом вычислительная операция выполняется в цикле, определяя число заполнения для всех направлений распространения.

Ниже, на рис. 5 представлены фазы эволюции клеточного автомата при столкновении частиц.

1 [

)-

а б

в г

Рис. 5. Реализация столкновения: а-б фаза столкновения, в-г фаза распространения

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

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

I

р(г,п) = ^ ат(гт,п), (4)

т=1

17

где I— число возможных направлений скоростей.

В случае потока, где все частицы движутся в одном направлении:

р(г1,п) = а;(г;,п).

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

рСр(г,п) =

ZN1 1

lm=1aí,m(rí,m,n),

где N — число узлов в заданном фрагменте сетки, N1 сетки.

Поток вещества определяется как:

N <5>

число частиц, оказавшихся в заданном фрагменте

у(г,п) =Рср(г,п)и(г,п) = ^ утат(гт,п),

(6)

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

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

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

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

Формула для расчета кинематической вязкости с учетом двумерного пространства и единичной скорости частиц представлена ниже:

_ 1

" 21р(1-2р)' (7) где р = рср(г, п)/1 - средняя плотность частиц в расчете на одно направление скорости.

Для назначения параметров вещества и геометрии, воспользуемся известными данными согласно [8](рис.6).

507,2

щ А

. Я . 228,6 С

Рис. 6. Схема расчетной области

При создании расчетной сетки НРР модели, с учетом дискретности и равномерности узлов, будем ориентироваться на количество ячеек в расчетах с использованием CFD пакетов: от 125 до 300 тыс. ячеек. В таком случае зададим сетку 1000 на 200 ячеек, которая будет имитировать протяженный канал (рис. 7). В нем движется основной поток частиц, частота которых Р1 неоднородна. Снизу на расстоянии Х2 = 460 находится отверстие инжекции шириной 10 ячеек.

При этом количество узлов в области составит 200 тыс., что сопоставимо с размерностью сетки

в Ашу8.

Из отверстия с единичной скоростью движется поток частиц с плотностью Р2. Отношение плотностей частиц Р±/Р2 пропорционально отношению их концентраций.

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

Рис. 7. Расчетная схема задачи

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

4 » _ _ Л_** * Д

n=500 n=700

Рис. 8. Распределение плотности частиц при неравномерном потоке

Также был проведен расчет на 14795 частиц для равномерно заданного потока, и инжектируемой струи более высокой плотности (рис. 9).

n=500 n=700

Рис. 9. Распределение плотности частиц при равномерном потоке

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

Для расчета в программном пакете Ansys была задана расчетная область с линейными размерами 507,2 на 76 мм (рис. 6). Левая граница инициализирована как входную (inlet), а правую как выходную (outlet). Также задается, как вход, граница отверстия, через которое производится вдув потока перпендикулярный оси канала. Используемый в расчете газ - воздух.

Выбрана модель турбулентности k — о> SST, как наиболее удобная и не требующая высокой точности расчетной сетки. Начальные условия реализуемой расчетной схемы (рис. 6) взяты из [8] и приведены в таблице.

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

Газодинамические параметры газов

Моо Рю, Па Тю, K Pi, Па Tf , K

2,61 133758 318,3 1044556 288.3

Ниже представлены поля безразмерных скоростей частиц в области вдуваемой струи, отнесенные к максимальному глобальному значению скорости на соответствующих шагах по времени - 100 и 500 (рис.10). Значения скорости усреднялись по заданному шагу по вертикальной и горизонтальной осям. Также для сравнения приведены поля скоростей, полученные в Ansys Fluent.

Рис. 10. Векторные поля скоростей для 100 и 500 шагов по времени

O.OOJ7iff \ 0.00525

Рис. 11. Поля скоростей полученные в Ansys Fluent

4.0 3.5 3.0 2.5

1.5 1.0

0.5 0.0

-200 -100 0 100 200

X, ММ

Рис. 12. Точки профиля давления, полученные с использованием HPPмодели

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

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

20

"

" , —+-+—

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

3,0

2,5

2,0

р

i

1,0

0,5

0,0

-100 -50 0 50 100

X, мм

Рис. 13. Профили давления, полученные в Ansys и экспериментальные точки

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

HPP модель может хорошо работать при описании течения жидкости или газа в канале с достаточно сложной геометрией, диффузии и распространения вещества (волны) в заданной области или при обтекании различных тел [9].

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

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

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

1. Быков Н.В., Бырдин К.А., Макаренко В.С. Влияние степени нерасчетности инжектируемой струи на силовые характеристики ракетного двигателя. Инженерный журнал: наука и инновации, 2017, вып. 2.

2. Котов М.А. Расчетно-экспериментальные исследования ударно-волновых процессов в гиперзвуковой ударной аэродинамической трубе // дис. .. .канд. физ.мат. наук. - ИПМех РАН, Москва, 2014 - 194с.

3. T. J. McIntyre, A. F. P. Houwing, P. C. Palma, P. A. B. Rabbath, J. S. Fox. Optical and Pressure Measurements in Shock Tunnel Testing of a Model Scramjet Combustor. JOURNAL OF PROPULSION AND POWER Vol. 13, No. 3, May - June 1997. pp. 388-394.

4. Бобков С.П., Соколов В.Л. Дискретное моделирование течения газа при пониженном давлении. Изв. вузов. Химия и хим. технология. 2017. Т. 60. Вып. 2. С. 79-84.

5. Hardy J., de Pazzis O., Pomeau Y. Molecular dynamics of a classical lattice gas: transport properties and time correlation functions // Physical Review. - 1976. - Vol. 13, N 5. - P. 1949-1961.

6. Бобков С.П., Соколов В.Л. Анализ возможностей применения решеточных моделей для исследования процессов в газах при пониженном давлении. Вестник ИГЭУ. 2015. Вып. 4. С. 58-63.

7. R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Physics Reports, vol.222,n.3, p.145, (1992).

SST о эксперимент сетка 1

- - ■сетка 2 сетка 3 сетка 4

i

. J Í

/ jK> -—

8. Яковчук М.С., Численное моделирование динамических процессов вдува струй в сверхзвуковую часть сопла // Вестник СГАУ. 2012. №3-1 (34).

9. Бандман О.Л., Клеточно-автоматные модели естественных процессов и их реализация на современных компьютерах // ПДМ. 2017. №1 (35).

Дмитриенко Даниил Богданович, специалист, оператор, era_1@mil.ru, Россия, Анапа, ФГАУ «ВИТ «ЭРА»,

Башанин Алексей Михайлович, специалист, оператор, Россия, Анапа, ФГАУ «ВИТ «ЭРА»

ANALYSIS OF THE APPLICABILITY OF THE LATTICE GAS MODEL D2Q4 TO SOLVE THE INJECTION

PROBLEM

D.B. Dmitrienko, A.M. Bashanin

The article considers the possibility of using a lattice gas model - a kind of cellular automaton with four degrees of freedom for modeling the injection of perpendicular supersonic flows. The construction and the principle of evolution of a cellular automaton is described. The model was verified on a test problem solved in Ansys Fluent.

Key words: cellular automaton, lattice gas, mathematical modeling, injection, Ansys Fluent, perpendicular flows.

Dmitrienko Daniil Bogdanovich, specialist, operator, era_1@mil.ru, Russia, Anapa, FGAU «MIT

«ERA»,

Bashanin Alexey Mihaylovich, specialist, operator, Russia, Anapa, FGAU «MIT «ERA»

УДК 621.3.046

DOI: 10.24412/2071-6168-2023-1-22-27

ДРОССЕЛЬ БУСТЕРНОЙ СХЕМЫ DC-DC-ПРЕОБРАЗОВАТЕЛЯ

В.В. Воронин

Представлен расчет дросселя бустерной схемы DC-DC преобразователя, на основе его технического задания - габаритные размеры, не более мм - 30*20*20; тип производства -серийный; группа эксплуатации - 2; ток, протекающий через нагрузку - 0,3 А; напряжение на нагрузке - 20 В; напряжение на входе схемы - 7 В; рабочая частота - 20 кГц; уровень помех - низкий. В работе представлены результаты электрического расчета, эскизного проекта конструктивного расчета дросселя на кольцевом сердечнике, включающий в себя: постановку задачи, выбор размеров сердечника, расчет числа витков и выбор материала сердечника, расчет параметров обмотки, расчет габаритных и установочных размеров дросселя.

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

Силовая электроника - стремительно развивающееся направление техники, целью которого являются снижение массы и габаритов устройств питания электронных средств. В устройствах питания, работающих на частоте 50 Гц, размеры трансформатора питания занимали от 20% до 50 % объема и массы устройств питания вследствие очень низкой частоты питающей сети. Увеличение частоты приводит к уменьшению габаритных размеров трансформаторов. Однако передавать электрическую энергию с частотой выше, чем 50 Гц на большие расстояния невозможно вследствие ее излучения в окружающее пространство. Поэтому преобразование промышленной частоты в более высокую производится непосредственно в устройствах питания. Современные преобразовательные устройства работают уже на частотах до 2 МГц.

Основными элементами всех преобразователей являются электронный ключ и дроссель или трансформатор. В качестве электронных ключей обычно используют транзисторы, работающие в ключевом режиме [1]. Последние могут работать только при постоянной полярности напряжения, то есть от источников постоянного напряжения. Этими источниками может быть выпрямленное напряжение питающей сети, химические и другие источники.

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

22