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

Георгий Чащин

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

Проделан первый шаг в освоении нового в науке численного метода, применительно к задачам акустики

Математика

Дипломы

Вуз: Московский физико-технический институт (государственный университет) (МФТИ)

ID: 5f39b5c5cd3d3e0001bb2197

UUID: 5ab12850-c23f-0138-1b68-0242ac180006

Язык: Русский

Опубликовано: около 4 лет назад

Просмотры: 8

10.66

Георгий Чащин

Московский физико-технический институт (государственный университет) (МФТИ)

Комментировать 11

Рецензировать 0

Скачать - 1,6 МБ

Поделиться работой

Федеральное государственное автономное образовательное учреждение высшего образования «Московский физико-технический институт (национальный исследовательский университет)» Физтех-школа Аэрокосмических Технологий Кафедра газовой динамики, горения и теплообмена Направление подготовки / специальность: 24.04.02 Системы управления движением и навигация (магистратура) Направленность (профиль) подготовки: Баллистика, аэродинамика и процессы управления летательных аппаратов РАЗРАБОТКА БЫСТРОГО ЧИСЛЕННОГО МЕТОДА РАСЧЕТА РАСПРОСТРАНЕНИЯ АКУСТИЧЕСКИХ ВОЗМУЩЕНИЯ НА ОСНОВЕ РЕШЕТОЧНЫХ УРАВНЕНИЙ БОЛЬЦМАНА (магистерская диссертация) Студент: Чащин Георгий Сергеевич _________________________ (подпись студента) Научный руководитель: Макаров Владимир Евгеньевич, канд. физ.-мат. наук _________________________ (подпись научного руководителя) Консультант (при наличии): _________________________ _________________________ (подпись консультанта) Москва 2019

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

Содержание Введение .................................................................................................................................... 3 Теоретический раздел ........................................................................................................ 5 I. 1. Основы кинетической теории ................................................................................................... 5 2. Решёточные уравнения Больцмана ........................................................................................... 9 3. Газовая динамика решёточной теории ................................................................................... 15 4. Акустика в решёточном газе ................................................................................................... 17 II. Практический раздел .................................................................................................... 21 Заключение.............................................................................................................................. 26 Список использованных источников ..................................................................................... 27 2

Введение Метод решеточных уравнений Больцмана (МРУБ) или LBM (LatticeBoltzmannMethod) – технология вычислительной газовой динамики, в основе которого лежат численные решения системы дискретизованных кинетических уравнений. Метод решеточных уравнений Больцмана прост в своей реализации и эффективен при масштабировании на массивно-параллельных системах [1621], что и обеспечило ему популярность. В работе [22] представлены возможности использования МРУБ для расчета акустических характеристик современных дозвуковых гражданских самолетов, типичный пример которого приведен на рисунке В.1. 3 1 2 4 3 4 Рисунок В.1 – Типичный современный дозвуковой гражданский самолет во взлетнопосадочной конфигурации с основными зонами генерации шума: 1 – вентилятор при наличии звукопоглощающих конструкций в его канале; 2 –выхлопная струя; 3 – элементы механизации крыла; 4 – шасси Расчеты акустических характеристик всего самолета проводились с учетом шума, излучаемого вентилятором [23-27], шума выхлопной струи при ее взаимодействии с планером [28-34], шума элементов механизации [35-39] и шума шасси [40-45]. 3

В настоящей работе сделан шаг в освоении теории МРУБ и рассмотрен один из простейших примеров ее применения в задаче распространения 2D аэродинамических возмущений. Автор благодарит сотрудника ВЦ РАН Ильина О.В. за проведенный курс лекций по теории МРУБ и сотрудника отдела 017 ЦИАМ Шорстова В.А. за помощь в построении алгоритма решения модельной задачи. . 4

I. Теоретический раздел 1. Основы кинетической теории В кинетической теории любой газ состоит из огромного числа частиц, движущихся с различными скоростями и в различных направлениях, причём размер каждой из частиц должен быть много меньше как характерного размера занимаемой газам области (или обтекаемого тела), так расстояния между самими частицами. Макроскопические параметры газа (как плотность и скорость) определяются суммой микроскопических параметров всех частиц, находящихся в этой точке. Столкновения частиц попарные, упругие и мгновенные. Для описания такой системы используется понятие функции распределения 𝑓(𝑡, 𝑥̅ , 𝜉 ̅), с помощью которой можно найти мгновенные значения плотности, плотности импульса и энергий, а также температура как мера хаотического движения: 𝜌(𝑡, 𝑥̅ ) = ∫ 𝑓(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ, 𝜌(𝑡, 𝑥̅ )𝑢̅ (𝑡, 𝑥̅ ) = ∫ 𝜉 ̅𝑓(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ, 2 𝜌(𝑡, 𝑥̅ )𝐸 (𝑡, 𝑥̅ ) = ∫|𝜉 ̅| 𝑓(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ (1.1) 3𝜌(𝑡, 𝑥̅ )𝑅𝑇 = ∫|𝑣̅ |2𝑓(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ, где (𝑣̅ (𝑡, 𝑥̅ ) = 𝜉 ̅ − 𝑢̅(𝑡, 𝑥̅ ) - скорость относительно потока. В отсутствии объёмных сил изменение такой системы описывается уравнением Больцмана 𝑑𝑓 𝑑𝑡 = 𝜕𝑓 𝜕𝑡 𝜕𝑓 + 𝜉 ̅ = Ω(𝑓), 𝜕𝑥̅ (1.2) где Ω(𝑓) - интеграл столкновений, обладающий свойством: ∫ Ω(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ = 0, ∫ 𝜉 ̅Ω(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ = 0, 2 ∫|𝜉 ̅| 𝜉Ω(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ = 0 ∫|𝑣̅ |2𝜉Ω(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ = 0. 5 (1.3)

В модели Батнагара-Гросса-Крука интеграл столкновений имеет вид Ω(𝑓 ) = 𝑓𝑒𝑞 − 𝑓 𝜏 (1.4) , где 𝑓 𝑒𝑞 (𝑡, 𝑥̅ , 𝑣̅ ) – равновесная функция распределения, определяемая как функция распределения Максвелла 𝑓 𝑒𝑞 = ρ √2π𝑅𝑇 v2 e −2RT (1.5) , а 𝜏 - время релаксации. То есть газ находиться в состоянии локального равновесия, а изменение картины течения связано с отклонением от равновесного состояния. Отметим важное свойство равновесной функции распределения – сепарабельность: 𝑒𝑞 𝑒𝑞 𝑓 𝑒𝑞 (𝑣̅ ) = 𝑓1 (𝑣𝑥 )𝑓1 (𝑣𝑦 ). 𝑒𝑞 Под 𝑓1 (1.6) следует понимать равновесную функцию в одномерном пространстве. Ещё отметим, что интегралы (1) при замене в них функции распределения равновесной не изменяются, что видно из формул (3). Важную роль в кинетической теории занимают следующие тензоры, называемые моментами: 𝑀 = ∫ 𝑓𝑑 2𝜉 = 𝜌 𝑀𝛼 = ∫ 𝜉𝛼 𝑓𝑑 2𝜉 = 𝜌𝑢𝛼 (1.7) 𝑀𝛼𝛽 = ∫ 𝜉𝛼 𝜉𝛽 𝑓𝑑 2 𝜉 𝑀𝛼𝛽𝛾 = ∫ 𝜉𝛼 𝜉𝛽 𝜉𝛾 𝑓𝑑 2𝜉 Третий тензор связан ещё с одним тензором: 𝑀𝛼𝛽 = 𝜌𝑢𝛼 𝑢𝛽 + ∫ 𝑣𝛼 𝑣𝛽 𝑓𝑑 2𝜉, (1.8) а 𝜎𝛼𝛽 = − ∫ 𝑣𝛼 𝑣𝛽 𝑓𝑑 2𝜉 – тензор напряжений, получаемый с помощью замены из (1.1). 6

Особое место занимает свёртка последнего тензора: 2 𝑀𝛼𝛽𝛽 = ∫ 𝜉𝛼 |𝜉 ̅| 𝑓𝑑 2𝜉 = 𝜌|𝑢̅|2𝑢𝛼 + 𝜌𝑅𝑇 + 2𝜎𝛼𝛽 𝑢𝛽 + ∫ 𝑣𝛼 |𝑣̅ |2𝑓𝑑 2𝜉, (1.9) 𝑞𝛼 = ∫ 𝑣𝛼 |𝑣̅ |2 𝑓𝑑 2𝜉 - вектор теплового потока. Проинтегрируем уравнение Больцмана 𝜕𝑓 𝜕𝑓 ∫ 𝜕𝑡 𝑑 2 𝜉 + ∫ 𝜉𝛾 𝜕𝑥 𝑑 2𝜉 = ∫ Ω(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ . 𝛾 (1.10) Произведя замену порядков интегрирования и дифференцирования и используя свойство интеграла столкновений, получим уравнение неразрывности 𝜕𝑀 𝜕𝑡 + 𝜕𝑀𝛾 𝜕𝑥𝛾 = 0 или 𝜕𝜌 𝜕𝑡 + 𝜕(𝜌𝑢𝛾 ) 𝜕𝑥𝛾 = 0. (1.11) Проинтегрировав уравнение Больцмана, умноженное на 𝜉𝛼 ∫ 𝜉𝛼 𝜕𝑓 2 𝜕𝑓 2 𝑑 𝜉 + ∫ 𝜉𝛼 𝜉𝛾 𝑑 𝜉 = ∫ 𝜉𝛼 Ω(𝑡, 𝑥̅ , 𝜉 ̅)d2 ξ 𝜕𝑡 𝜕𝑥𝛾 (1.12) После аналогичных преобразований получим уравнение Навье-Стокса 𝜕𝑀𝛼 𝜕𝑡 + 𝜕𝑀𝛼𝛾 𝜕𝑥𝛾 = 0 или 𝜌 𝜕𝑢𝛼 𝜕𝑡 + 𝜌𝑢𝛽 𝜕𝑢𝛼 𝜕𝑥𝛽 = 𝜕𝜎𝛼𝛾 𝜕𝑥𝛾 (1.13) По тем же соображениям имеем из уравнения ∫ 𝜉𝛼 𝜉𝛼 𝜕𝑓 2 𝜕𝑓 2 𝑑 𝜉 + ∫ 𝜉𝛼𝜉𝛼 𝜉𝛾 𝑑 𝜉 = ∫ 𝜉𝛼 𝜉𝛼 Ω(𝑡, 𝑥̅ , 𝜉 ̅)d2ξ 𝜕𝑡 𝜕𝑥𝛾 (1.14) закон сохранения энергии в двух видах: 𝜕(𝜌𝐸) 𝜕𝑀𝛼𝛼𝛾 𝜕(𝜌𝐸) 𝜕(𝜌𝐸𝑢𝛽 ) + = 0 или + 𝜕𝑡 𝜕𝑥𝛾 𝜕𝑡 𝜕𝑥𝛽 𝜕(𝑢𝛼 𝜎𝛼𝛽 ) 𝜕𝑞𝛽 = + 𝜕𝑥𝛽 𝜕𝑥𝛽 (1.15) Точный вид уравнения Навье-Стокса определяется разложением Чепмена-Энскога, когда предполагаем справедливость разложения ∞ 𝑓 = ∑ 𝜀 𝑛 𝑓 (𝑛), 𝑛=0 где 𝜀 ≪ 1 − число Кнудсена. 7 (1.16)

Нулевой член ряда есть равновесная функция, ограничивая разложение нулевым порядком, получим уравнения Эйлера. Рассматривая нулевой и первый порядок – Уравнение Навье-Стокса для вязкого газа. 8

2. Решёточные уравнения Больцмана В методе решёточных уравнений Больцмана частицы имеют конечный набор возможных скоростей. Расчётная область покрывается квадратной сеткой с постоянным шагом Δℎ. Производится дискретизация по времени с постоянным шагом Δ𝑡. Получаются следующие скорости: 𝑐̅0 = (0; 0), 𝑐̅3 = (− 𝑐̅6 = (− 𝑐̅1 = ( Δℎ ; 0), Δ𝑡 Δℎ Δℎ Δ𝑡 ; Δ𝑡 Δℎ ; 0), Δ𝑡 𝑐̅4 = (0; − ), 𝑐̅7 = (− Δℎ Δ𝑡 Δℎ Δ𝑡 𝑐̅2 = (0; Δℎ Δ𝑡 Δ𝑡 ), Δℎ Δℎ ; ), Δ𝑡 Δ𝑡 Δℎ Δℎ 𝑐̅8 = ( ; ) Δ𝑡 Δ𝑡 ), ;− Δℎ 𝑐̅5 = ( ), (2.1) Такая модель носит название D2Q9, где сначала указана размерность пространства, а затем число скоростей. Эта модель образуется как прямое (декартово) произведение модели D1Q3 самой на себя. Одномерная модель представлена следующими скоростями: с0 = 0, 𝑐1 = Δℎ Δ𝑡 , 𝑐−1 = − Δℎ Δ𝑡 (2.2) В связи с этим допущением преобразуется и уравнение Больцмана к виду: 𝑑𝑓𝑖 𝑑𝑡 = 𝜕𝑓𝑖 𝜕𝑡 + 𝑐̅𝑖 𝜕𝑓𝑖 𝜕𝑥̅ 𝑒𝑞 = 𝑓𝑖 − 𝑓𝑖 𝜏 для 𝑖 = 0,1, … ,8 , (2.3) Здесь индекс нумерует компоненты функций распределения для соответствующей скорости частиц. Также меняются определения моментов: 𝑀 = ∑𝑖 𝑓𝑖 = 𝜌, 𝑀𝛼 = ∑𝑖 𝑓𝑖 𝑐𝑖𝛼 = 𝜌𝑢𝛼 , 𝑀𝛼𝛽 = ∑𝑖 𝑓𝑖 𝑐𝑖𝛼 𝑐𝑖𝛽 , 𝑀𝛼𝛽𝛾 = ∑𝑖 𝑓𝑖 𝑐𝑖𝛼 𝑐𝑖𝛽 𝑐𝑖𝛾 . 9 (2.4)

Замена интегрирования суммированием необходимо также и в свойствах интеграла столкновений, однако: ∑𝑖 Ω(𝑓𝑖 ) = 0, ∑𝑖 Ω(𝑓𝑖 ) 𝑐̅𝑖 = 0, (2.5) ∑𝑖 Ω(𝑓𝑖 )|𝑐̅𝑖 |2 = 0 , ∑𝑖 Ω(𝑓𝑖 )|𝑐̅𝑖 − 𝑢̅|2 𝑐̅𝑖 = 0 . Равновесная функция такой статистики тоже отличается, однако введя в рассмотрение решёточную скорость звука (с помощью которой можно ввести и температуру решётки), разложим функцию Максвелла в ряд Тейлора по скорости потока (что вводит ограничение на число Маха, квадрат которого должен быть много меньше единицы): 𝑓𝑖 𝑒𝑞 = |𝑢 ̅|3 𝑂( где 𝑤𝑖 ρ 𝑐𝑠3 − e 2 √2π𝑐𝑠 ̅ )2 (c̅i − u 2 2cs )] = 𝑤𝑖 𝜌 [1 + = ̅ 𝑐̅𝑖 𝑢 𝑐𝑠2 − ρ e 2 √2π𝑐𝑠 + (|c̅i |)2 2c2 s ̅)2 (𝑐̅𝑖 𝑢 2𝑐𝑠4 − [1 + |𝑢 ̅|2 2𝑐𝑠2 ̅ 𝑐̅𝑖 𝑢 𝑐𝑠2 + |𝑢 ̅|3 + 𝑂( 𝑐𝑠3 ̅)2 (𝑐̅𝑖 𝑢 2𝑐𝑠4 − |𝑢 ̅|2 2𝑐𝑠2 + (2.6) )], - весовые коэффициенты, которые просто определить для одномерного случая. Равновесной функцией распределения в этом случае будет: 𝑒𝑞 𝑓0 = 𝑤0𝜌 [1 − 𝑒𝑞 𝑓1 𝑢2 2𝑐𝑠2 ], 𝑐𝑢 (𝑐𝑢)2 𝑢2 = 𝑤1𝜌 [1 + 2 + − ] 𝑐𝑠 2𝑐𝑠2 2𝑐𝑠2 𝑒𝑞 𝑓−1 = 𝑤−1𝜌 [1 − 𝑐𝑢 (𝑐𝑢)2 𝑐𝑠 2𝑐𝑠2 2 + − 𝑢2 2𝑐𝑠2 (2.7) ] Однако из неё должны получаться моменты, то есть: 𝑒𝑞 𝑒𝑞 𝑒𝑞 𝑓−1 + 𝑓0 + 𝑓1 𝑒𝑞 𝑒𝑞 = 𝜌𝑢, 𝑒𝑞 = 𝜌𝑢2 + 𝜌𝑐𝑠2 , −𝑐𝑓−1 + 𝑐𝑓1 𝑒𝑞 𝑐 2𝑓−1 + 𝑐 2 𝑓1 где 𝑐 = Δℎ Δ𝑡 = 𝜌 , (2.8) – модуль скорости. 10

Решая эту систему получим, что 𝑒𝑞 𝑓0 𝑒𝑞 𝑓1 = 𝜌(1 − = 𝜌( 𝑢 𝑐𝑠2 𝑐 𝑐2 + 2𝑐 𝑢 𝑒𝑞 𝑓−1 = 𝜌(− 𝑢2 2𝑐 − 2 ), 𝑢2 𝑐𝑠2 2𝑐 2𝑐 2 + − 2 ), 𝑢2 𝑐𝑠2 2𝑐 2𝑐 2 − 2 (2.9) ) , Однако, очевидно, что свёртка третьего момента будет воспроизведена неточно, и подстановкой получим равенство: 𝜌𝑢𝑐 2 = 3𝜌𝑐𝑠2 + 𝜌𝑢3 . (2.10) Пренебрегая кубической ошибкой, получим соотношения для скорости звука решётки, необходимого для наилучшего определения моментов: 𝜌𝑐𝑠 = 1 √3 Δℎ 𝑐= √3Δ𝑡 . (2.11) Подставив (10) в формулу (8) и сравнив с выражением (6), получим 2 1 3 6 весовые коэффициенты в одномерном случае:𝑤0 = , 𝑤−1 = 𝑤1 = . Из свойства сепарабельности равновесной функции её весовые коэффициенты получим перемножением коэффициентов одномерной. Итого: 4 1 9 9 𝑤0 = , 𝑤1 = 𝑤2 = 𝑤3 = 𝑤4 = , 𝑤5 = 𝑤6 = 𝑤7 = 𝑤8 = 1 36 (2.12) Получим конечно-разностную схему для уравнения Больцмана. Для этого в уравнении (3) произведём замену на новую переменную 𝜁 и получим 𝑑𝑓𝑖 𝑑𝜁 С = 𝜕𝑓𝑖 𝑑𝑡 𝜕𝑡 𝑑𝜁 производными + 𝑑𝑡 𝑑𝜁 𝜕𝑓𝑖 𝑑𝑥̅ 𝜕𝑥̅ 𝑑𝜁 =1 и 𝑒𝑞 = 𝑓𝑖 − 𝑓𝑖 𝜏 𝑑𝑥̅ 𝑑𝜁 (2.13) . = 𝑐̅𝑖 , проинтегрировав которые получим 𝑡 = 𝜁 + 𝑡0 и 𝑥̅ = 𝑐̅𝑖 + 𝑥̅0 - уравнение характеристики из точки 𝑥̅0 в момент времени 𝑡0 . Таким образом получаем обыкновенное дифференциальное уравнение: 𝑒𝑞 𝑑𝑓𝑖 1 𝑓 = − 𝑓𝑖 + 𝑖 𝑑𝜁 𝜏 𝜏 11 (2.14)

Решение этого уравнение может быть приведено в квадратурах следующим образом: 𝑓𝑖 (𝜁) = 𝑒 − 𝜁 𝜏 (𝐶 + 𝑒𝑞 𝜁 𝑓𝑖 (𝜁′) 𝜁′ 𝑒 𝜏 𝑑𝜁′) , ∫𝜁 𝜏 0 𝐶 = 𝑐𝑜𝑛𝑠𝑡, 𝜁0 = 𝑐𝑜𝑛𝑠𝑡 (2.15) Произвольную константу определим, приняв 𝜁0 за начальную точку: 𝜁0 𝑓𝑖 (𝜁0) = 𝑒 − 𝜏 𝐶. Выразив С и проведя несложные алгебраические преобразования получим 𝑓𝑖 (𝜁) = 𝑒 − 𝜁− 𝜁0 𝜏 𝑒𝑞 (𝑓𝑖 (𝜁0) + 𝜁 𝑓 (𝜁′) 𝜁 ∫𝜁 𝑖 𝜏 𝑒 0 ′− 𝜁 𝜏 0 (2.16) 𝑑𝜁′). Вернёмся к старым переменным с учётом произвольности выбора начального момента и проинтегрируем уравнение за один шаг по времени, то есть возьмём 𝜁0 = 𝑡, 𝜁 = 𝑡 + Δ𝑡, 𝜁 ′ = 𝑡 ′ : 𝑓𝑖 (𝑡 + Δ𝑡, 𝑥̅ + 𝑐̅𝑖 Δ𝑡) = 𝑒 ′ 𝑒𝑞 Δ𝑡 −𝜏 (𝑓𝑖 (𝑡, 𝑥̅ ) + 𝑡+Δ𝑡 𝑓𝑖 (𝑡′) 𝑡 −𝑡 𝑒 𝜏 𝑑𝑡′). ∫𝑡 𝜏 (2.17) Разложив экспоненциальную функцию по Тейлору и аппроксимировав интеграл по формуле трапеции получим: 𝑓𝑖 (𝑡 + Δ𝑡, 𝑥̅ + 𝑐̅𝑖 Δ𝑡) Δ𝑡 Δ𝑡 2 Δ𝑡 Δ𝑡 Δ𝑡 2 𝑒𝑞 = (1 − + 2) 𝑓𝑖 (𝑡, 𝑥̅ ) + ((1 − + 2 ) 𝑓𝑖 (𝑡, 𝑥̅ ) 𝜏 2𝜏 2𝜏 𝜏 2𝜏 (2.18) 𝑒𝑞 + 𝑓𝑖 (𝑡 + Δ𝑡, 𝑥̅ + 𝑐̅𝑖 Δ𝑡)) Δ𝑡 Произведём замену 𝑓̃𝑖 = 𝑓𝑖 + Ω(𝑓𝑖 ) (или то же самое, чтои 𝑓𝑖̅ = 2 𝑒𝑞Δ𝑡 𝑓𝑖 + 𝑓𝑖 2𝜏 Δ𝑡 1+ 2𝜏 ) после алгебраических упрощений, получим: 𝑒𝑞 𝑓 𝑓𝑖̅ (𝑡 + Δ𝑡, 𝑥̅ + 𝑐̅𝑖 Δ𝑡) = 𝑓𝑖̅ (𝑡, 𝑥̅ ) + Δ𝑡 𝑖 (𝑡,𝑥̅ )− 𝑓̅𝑖 (𝑡,𝑥̅ ) 𝜏+ Δ𝑡 2 + 𝑂(Δt 2 . (2.19) Заметим, что из свойств интеграла столкновения, что такая замена не изменяет значения моментов от функции распределения, поэтому различий между двумя функций можно не делать, а тождество (18) есть конечноразностная схема второго порядка для шага по времени для уравнения Больцмана. 12

Подводя итог, имеем алгоритм решения задачи: 1. На старом временном слое в каждом узле решётки вычисляется 𝑒𝑞 из значений 𝜌 и 𝑢̅ равновесная функция 𝑓𝑖 . 2. Определяем эффект столкновений каждого вида частиц в решётке по формуле :𝑓𝑖∗(𝑡, 𝑥̅ ) = 𝑓𝑖 (𝑡, 𝑥̅ ) + Δ𝑡 𝑒𝑞 𝑓𝑖 (𝑡,𝑥̅ )− 𝑓𝑖 (𝑡,𝑥̅ ) 𝜏+ . Δ𝑡 2 3. Переходим на новый временной слой, находя распространение частиц для каждой скорости по узлам решётки: 𝑓𝑖 (𝑡 + Δ𝑡, 𝑥̅ + 𝑐̅𝑖 Δ𝑡) = 𝑓𝑖∗ (𝑡, 𝑥̅ ). 4. Уже на новом временном слое определяем газодинамические параметры с помощью первых моментов: 𝜌 = ∑𝑖 𝑓𝑖 и 𝑢̅ = 1 𝜌 Важное ∑𝑖 𝑓𝑖 𝑐̅𝑖 . замечание: в начальный момент времени 𝑒𝑞 распределения считается равной равновесной, то есть 𝑓𝑖 (0, 𝑥̅ ) = 𝑓𝑖 . 14 функция

3. Газовая динамика решёточной теории Для получения уравнений газовой динамики решётки воспользуемся уравнением Больцмана с введённым в него числом Кнудсена, являющегося малым параметром: 𝑒𝑞 𝑑𝑓𝑖 𝜕𝑓𝑖 𝜕𝑓𝑖 𝑓 − 𝑓𝑖 = + 𝑐̅𝑖 = 𝑖 , 𝑑𝑡 𝜕𝑡 𝜕𝑥̅ 𝜏0 𝜀 (3.1) где 𝜏0 − характерное время задачи. И разложением функции распределения Чепмена-Энскога до первого порядка по числу Кнудсена: 𝑓𝑖 = 𝑓𝑖 𝑒𝑞 + 𝜀𝑓𝑖 (1) , в котором за нулевой член разложения принята равновесная функция распределения, а первый член определяется подстановкой выражения в (19) и получается: .𝑓𝑖 (1) 𝑒𝑞 = −𝜏0 ( 𝜕𝑓𝑖 𝜕𝑡 𝑒𝑞 + 𝑐̅𝑖 𝜕𝑓𝑖 𝜕𝑥̅ 𝑒𝑞 ) = −𝜏0 𝑑𝑓𝑖 𝑑𝑡 . (3.2) Заметим также из свойств интеграла столкновений (2.5) следуют дополнительные свойства первого приближения: ∑𝑖 𝑓i (1) = 0 и ∑𝑖 𝑓i (1) 𝑐̅𝑖 = 0 . (3.3) Опираясь на (3.3) просуммируем уравнение (3.2) по всем скоростям решётки, умноженное на 1, 𝑐̅𝛼 и 𝑐̅𝛼 𝑐̅𝛽 , тогда получим соотношения для моментов от равновесной функции первого приближения: 𝜕𝑀𝑒𝑞 𝜕𝑡 𝑒𝑞 + 𝜕𝑡 + 𝑒𝑞 𝜕𝑀𝛼𝛽 𝜕𝑡 𝜕𝑥𝛾 =0, 𝑒𝑞 𝑒𝑞 𝜕𝑀𝛼 𝜕𝑀𝛾 𝜕𝑀𝛼𝛽 𝜕𝑥𝛽 = 0, 𝑒𝑞 + 𝜕𝑀𝛼𝛽𝛾 𝜕𝑥𝛾 = − (3.4) 1 𝜏0 (1) 𝑀𝛼𝛽 . Моменты от равновесной функции вычисляются и самой функции (2,6) с определёнными весовыми коэффициентами (2,12): 𝑒𝑞 𝑀𝑒𝑞 = ∑𝑖 𝑓𝑖 = 𝜌, 𝑒𝑞 𝑒𝑞 𝑀𝛼 = ∑𝑖 𝑓𝑖 𝑐𝑖𝛼 = 𝜌𝑢𝛼 , (3.5) 𝑒𝑞 𝑒𝑞 𝑀𝛼𝛽 = ∑𝑖 𝑓𝑖 𝑐𝑖𝛼 𝑐𝑖𝛽 = 𝜌𝑢𝛼 𝑢𝛽 + 𝜌𝑐𝑠2𝛿𝛼𝛽 , 15

𝑒𝑞 𝑒𝑞 𝑀𝛼𝛽𝛾 = ∑𝑖 𝑓𝑖 𝑐𝑖𝛼 𝑐𝑖𝛽 𝑐𝑖𝛾 = 𝜌𝑐𝑠2(𝑢𝛼 𝛿𝛽𝛾 + 𝑢𝛽 𝛿𝛼𝛾 + 𝑢𝛾 𝛿𝛼𝛽 ). Теперь видно, что первое соотношение из (3,4) есть уравнение неразрывности, а второе – уравнение Эйлера: 𝜕𝜌 𝜕(𝜌𝑢𝛾 ) 𝜕(𝜌𝑢𝛼 ) 𝜕(𝜌𝑢𝛼 𝑢𝛽 ) 𝜕𝑝 + =0 и + = − , 𝜕𝑡 𝜕𝑥𝛾 𝜕𝑡 𝜕𝑥𝛽 𝜕𝑥𝛼 (3.6) где 𝑝 = 𝜌𝑐𝑠2 . Выражение же последнего момента из (3.5) связано с трудоёмкими преобразованиями производных и группировкой слагаемых, но в конечном итоге получается следующее соотношение: (1) 𝑀𝛼𝛽 = ∑𝑖 𝑓𝑖 (1) 𝑐𝑖𝛼 𝑐𝑖𝛽 = − 𝜌𝑐𝑠2 𝜏0 ( 𝜕𝑢𝛼 𝜕𝑥𝛽 + 𝜕𝑢𝛽 𝜕𝑥𝛼 ) − 𝜏0 𝜕(𝜌𝑢𝛼 𝑢𝛽 𝑢𝛾 ) 𝜕𝑥𝛾 . (3.7) Умножив уравнение (3.1) на 𝑐̅𝛼 и просуммировав по скоростям, используя разложение Чепмена-Энскога и его свойства увидим уравнение Навье-Стокса в моментной форме: 𝑒𝑞 𝑒𝑞 𝜕𝑀𝛼 𝜕𝑡 + Использовав 𝜕𝑀𝛼𝛽 𝜕𝑥𝛽 ранее (1) + 𝜀 𝜕𝑀𝛼𝛽 𝜕𝑥𝛽 (3.8) =0. полученные выражения напишем его в газодинамических переменных: 𝜕(𝜌𝑢𝛼 ) 𝜕(𝜌𝑢𝛼 𝑢𝛽 ) + 𝜕𝑡 𝜕𝑥𝛽 𝜕(𝜌𝑢𝛼 𝑢𝛽 𝑢𝛾 ) 𝜕𝑝 𝜕 𝜕𝑢𝛼 𝜕𝑢𝛽 = − + 𝜈 (𝜌 ( + )) + 𝜏 , 𝜕𝑥𝛼 𝜕𝑥𝛽 𝜕𝑥𝛽 𝜕𝑥𝛼 𝜕𝑥𝛾 (3.9) где 𝑝 = 𝜌𝑐𝑠2 𝜈 = 𝑐𝑠2𝜏. Как видно из последней формулы, закон сохранения импульса газа в решёточной теории содержит лишнее слагаемое с кубом скорости, но учитывая ограничение по скорости газа и малый коэффициент перед ним, можно пренебречь данной ошибкой. Также в вязкой составляющей отсутствует слагаемое с дивергенцией скорости, тем ни менее при сделанных ограничениях сжимаемость газа довольно слабая, и эта ошибка тоже вносит незначительный вклад. Стоит обратить внимание, что уравнение Эйлера 16

воспроизводятся верно, в решёточной теории вязкость обязана быть ненулевой. 4. Акустика в решёточном газе Рассмотрим распространение малых колебаний в газе. Для этого в уравнении неразрывности из (3.6) и в уравнении Навье-Стокса (3.9) примем следующие представления для плотности и компонент скорости 𝜌 = 𝜌0 + 𝜌 ′ и 𝑢𝛼 = 𝑢′𝛼 , причем должны быть верны оценки: 𝜌 ′ ≪ 𝜌0 = 𝑐𝑜𝑛𝑠𝑡 и √𝑢′𝛼 𝑢′𝛼 ≪ 𝑐𝑠 . Кроме того, проведём линеризацию уравнений, отбросив слагаемые более высокого порядка малости, чем первый, в результате чего они примут следующий вид 𝜕𝜌′ 𝜕𝑡 + 𝜌0 𝜕𝑢′𝛾 𝜕𝑥𝛾 = 0 и 𝜌0 𝜕𝑢′𝛼 𝜕𝑡 = − 𝑐𝑠2 𝜕𝜌′ 𝜕𝑥𝛼 + 𝜈𝜌0 𝜕 𝜕𝑥𝛽 ( 𝜕𝑢′𝛼 𝜕𝑥𝛽 + 𝜕𝑢′𝛽 𝜕𝑥𝛼 ) (4.1) Для того чтобы выразить неизвестные переменные через одну, удобно ввести потенциал скорости согласно 𝑢𝛼 = 𝜕𝜑 𝜕𝑥𝛼 . Теперь можно выразить перепад плотности и второго уравнения: 𝜌 = − 𝜌0 𝜕𝜑 𝑐𝑠2 𝜕𝑡 + 2𝜈𝜌0 𝑐𝑠2 Δx φ . (4.2) Подставляя полученное выражение в уравнение неразрывности из (4.1) получим: Δx φ − Полученное 1 ∂2 φ 𝑐𝑠2 ∂𝑡 2 + 2𝜈 𝜕Δx φ 𝑐𝑠2 уравнение 𝜕𝑡 =0. называется (4.3) волновым и описывает распространение звука в газе. Рассмотрим движение звуковой волны постоянной частоты, то есть подставим 𝜑(𝑡, 𝑥̅ ) = 𝜓(𝑥̅ )𝑒 −𝑖𝜔𝑡 :, где под i необходимо понимать мнимую единицу, а не номер скорости в решётке, но поскольку мнимая единица никогда не стоит в индексе, а номер скорости только в индексе, то и путаницы в обозначениях возникнуть не должно. Δx ψ + ω2 𝑐𝑠2 ψ− i 2𝜈𝜔 𝑐𝑠2 Δx ψ = 0. 17 (4.4)

Данное уравнение является уравнением Гельмгольца, которое можно привести к стандартному виду используя обозначение волнового числа𝑘 = 𝜔 𝑐𝑠 . За одно сделаем и допущение, что Δx ψ + (𝑘 √1 + i 2𝜈𝜔 2 ) ψ 𝑐𝑠2 2𝜈𝜔 𝑐𝑠2 ≪ 1. = 0. (4.5) Теперь видно, что распределение в пространстве амплитуды звуковой волны постоянной частоты в вязкой среде описывается уравнением Гельмгольца с комплексным волновым числом. С его учётом можно преобразовать выражение для пульсации плотности: 𝜌′ = 𝜌0 𝜔 𝑐𝑠2 (− 𝜔 2𝜈𝜔2 𝑐𝑠 𝑐𝑠4 2 + i(1 − ))𝜓𝑒 −𝑖𝜔𝑡 . (4.6) В полярных координатах (под 𝑟 = |𝑥̅ | следует понимать модуль радиус-вектора, 𝜃 - полярный угол) в открытом пространстве решение уравнения Гельмгольца может быть получено методом разделения переменных для любых значений волнового вектора и записано c использованием функций Ганкеля первого вида в конструкции: 𝜓 = ∑∞ 𝑛= −∞ 𝐶𝑛 √1+ 2𝜈𝜔 i 2 𝑐𝑠 𝑛 𝐻𝑛 (𝑘𝑟√1 + i 2𝜈𝜔 𝑐𝑠2 ) 𝑒 𝑖𝑛𝜃 , где 𝐶𝑛 = 𝑐𝑜𝑛𝑠𝑡. (4.7) Подставляя значение вязкости равной нулю, получается решение в невязкой среде, обозначаемое индексом «0»: 𝑖𝑛𝜃 𝜓0 = ∑∞ . 𝑛= −∞ 𝐶𝑛 𝐻𝑛 (𝑘𝑟) 𝑒 Введём в рассмотрение коэффициент поглощения 𝜇 = (4.8) 𝜈𝜔2 𝑐𝑠3 . С его помощью можно показать, что решение (4.7) стремиться к виду 𝜓 = 𝑒 −𝜇𝑟 𝜓0c удалением от источников. Сравнив выражения (4.7) и (4.8) достаточно показать, что 1 || √1 + i 2𝜈𝜔 𝑛 𝐻𝑛 (𝑘𝑟√1 + i 2𝜈𝜔 ) − 𝑒 −𝜇𝑟 𝐻𝑛 (𝑘𝑟) || → 0 2 𝑐𝑠 𝑐𝑠2 при 𝑟 → ∞. 18 (4.9)

Вспомним некоторые разложения Ломмеля для функций Бесселя и Неймана: 1 ∞ 𝑛 (−1)𝑚 ( ℎ𝑧) 2 𝐽𝑛 (𝑧√1 + ℎ) = √1 + ℎ ∑ 𝑚! 𝑚=0 ∞ 𝑛 𝑌(𝑧√1 + ℎ) = √1 + ℎ ∑ 𝑚 𝑚 1 (−1)𝑚 ( ℎ𝑧) 2 𝑚! 𝑚=0 𝐽𝑛+𝑚 (𝑧) (4.10) 𝑌𝑛+𝑚 (𝑧) С их использованием и разложением экспоненты в ряд Тейлора предел преобразуется ∞ (−1)𝑚 (𝜇𝑟)𝑚 𝑚 |∑ (𝑖 𝐻𝑛+𝑚 (𝑘𝑟) − 𝐻𝑛 (𝑘𝑟)) | → 0 𝑚! (4.11) 𝑚=0 при 𝑟 → ∞. Для функций Ганкеля справедливо следующая асимптотическая формула: 𝜋 𝐻𝑛 (𝑧) = √ 𝑒 𝑖(𝑧− 2𝑧 𝑛𝜋 𝜋 − 4) 2 (∑𝑀−1 𝑙=0 (−1)𝑚 (𝑛,𝑙) (2𝑖𝑧)𝑙 + 𝑂 (𝑧 −𝑀 )) , (4.12) 1 где (𝑛, 𝑙) = Γ(𝑛+𝑙+2) 1 . 𝑙!Γ(𝑛−𝑙+2) Окончательно вид оценки представляется в форме 𝜋 √2𝑘𝑟 |∑∞ 𝑚=0 (−1)𝑚 (𝜇𝑟)𝑚 𝑚! 𝜋 (𝑖 𝑚 − 1)| = √ |𝑒 𝑖𝜇𝑟 − 𝑒 −𝜇𝑟 | → 0 , 2𝑘𝑟 (4.13) при 𝑟 → ∞. Подводя итог, приходим к выводу, что решёточной теорией возможно описание распространения звуковых волн в газе, но чтобы картина течения соответствовала физической необходимо подобрать основные параметры схемы так, чтобы скорость и коэффициент поглощения звука в решётке были равны их физическим значениям. 20

II. Практический раздел Исследование возможностей моделирования акустических возмущений методом решёточных уравнений Больцмана проводится для серии частот (f) в квадратной расчётной области с сеткой, состоящей из квадратных ячеек, и точечного источника. Основными параметрами такой решётки являются: число ячеек, приходящееся на сторону расчётной области (N), число ячеек, приходящееся на длину волны (n), время релаксации (𝜏). Также используются и другие параметры решётки. Скорость звука приравнена к физической, все размерные величины выражены в СИ. Число ячеек всегда чётное, источник расположен в центре правой верхней ячейки, являющейся ближайшей к центру квадрата. Мощность источника подобрана так, чтобы амплитуда перепада плотности в узлах ячейки, содержащей источник. Используемые параметры задачи приведены в таблице 1. Таблица 1 – Параметры рассмотренной задачи f, Гц N n 𝜏 𝜐 100 300 15 1.5*10−8 0.1089 500 300 15 1.5*10−8 0.1089 1000 300 15 1.5*10−8 0.1089 5000 300 15 1.5*10−8 0.1089 10000 300 15 1.5*10−8 0.1089 21

На рисунке 1 представлены примеры распространения волн от точечного источники с частотой f=500 Гц на сетке, содержащей 15 ячеек на дину волны, для разных значениях вязкости. Минимальная вязкость Максимальная вязкость Рисунок 1 – Распространения волн от точечного источники с частотой f=500 Гц на сетке, содержащей 15 ячеек на дину волны, для разных значениях вязкости 22

На рисунках 2 – 5 представлено сопоставление полученных результатов расчетов с различными вариантами аналитического решения. Рисунок 2 – Сопоставление полученных результатов расчетов с различными вариантами аналитического решения (частота f=100 Гц, сетка –15 ячеек на длину волны, минимальная вязкость) Рисунок 3 – Сопоставление полученных результатов расчетов с различными вариантами аналитического решения (частотаf=500 Гц,сетка – 15 ячеек на длину волны, минимальная вязкость) 23

Рисунок 4 – Сопоставление полученных результатов расчетов с различными вариантами аналитического решения (частота f=5000 Гц, сетка –15 ячеек на длину волны, минимальная вязкость) Рисунок 5 – Сопоставление полученных результатов расчетов с различными вариантами аналитического решения (частота f=10000 Гц, сетка –15 ячеек на длину волны, минимальная вязкость) 24

На рисунке 6 представлено сопоставление численных решений, полученных на 2-х различных сетках с аналитическим экспоненциальным приближением. Весь график Фрагмент графика Рисунок 6 – Сопоставление полученных результатов расчетов с с аналитическим экспоненциальным приближением (частота f=500 Гц, сетки 1 и 2 –15 и 20 ячеек соответственно на длину волны, минимальная 25

вязкость) Заключение 1. Освоена теория перспективного численного метода, позволяющего строить чрезвычайно экономичные и эффективные алгоритмы решения широкого круга аэродинамических и акустических задач для дозвуковых течений газа на декартовых сетках. 2. В качестве примера использования метода создано программное приложение, позволяющее рассчитывать распространение двумерных акустических возмущения и получены различные аналитические решения по распространению акустических волн от точечного источника, необходимыедля оценки точностиполучаемых численных решений. 3. Выполненные расчеты продемонстрировали эффективность созданного алгоритма в широком диапазоне параметров метода и задачи и его высокую точность: так, в частности, сопоставление численного и аналитического решений задачи о распространении акустических возмущений от точечного источника, свидетельствуют о том, что наблюдаемые отличие для сеток, содержащих 15-20 ячеек на длину волны, не превосходят величины 0.8 дБ. 4. Выполненная работа дает возможность перейти непосредственно к созданию быстрого метода распространения пространственных акустических возмущения и к анализу детальных свойств исходного метода с целью реализации возможности его применения для расчета турбулентных течений с применением вихреразрешающих подходов, что, в свою очередь, открывает перспективны расчета шума сложных практически важных конфигураций. 26

Список использованных источников 1. D. Casalino, A. Hazir. “Building a Numerical Turbofan Demonstrator Exa Corporation. DGLR DEGA X-Noise Workshop, Braunschweig”, November 11th (2015). 2. Mann, A., Pérot, F., Kim, M.-S., Casalino, D. and Fares, E., “Ad- vanced Noise Control Fan Direct Aeroacoustics Predictions using a LatticeBoltzmann Method”, AIAA-2012- 2287, 2012. 3. Pérot, F., Freed, D., Mann, A., “Acoustic absorption of porous materi- als using LBM”, AIAA 2013-2070. 4. Pérot F., Mann, A., Kim, M.-S., Casalino, D. and Fares, E., “Investi- gation of inflow condition effects on the ANCF aeroacoustics radiation using LBM”, InterNoise, New York, 2012. 5. Mann, A., Pérot F., Kim, M.-S and Casalino, “Characterization of Acoustic Liners Absorption using a Lattice-Boltzmann Method”, AIAA2013- 2271, 2013. 6. Casalino, D. Ribeiro A.F.P.et al. “Towards Lattice-Boltzmann Predic- tion of Turbofan Engine Noise”, AIAA-2014-3101. 7. Lew, P., Lyrintzis, A. S., and Mongeau, L., “Noise Prediction of a Subsonic Turbulent Round Jet using the Lattice-Boltzmann Method,” Journal of the Acoustical Society of America, Vol. 128, No. 3, September 2010, pp. 1118–1127. 8. Lew, P., Najafi-Yazdi, A., and Mongeau, L., “Unsteady Numerical Simulation of a Round Jet with Impinging Microjets for Noise Suppression,” Journal of the Acoustical Society of America, Vol. 133, No. 6, June 2013. 9. Habibi, K., Gong, H., Najafiyazdi, A. and Mongeau, L., “Prediction of the Sound Radiated from Low-Mach Internal Mixing Nozzles with Forced Mixers using the Lattice Boltzmann Method”, AIAA paper 2013-2143, 2013. 27

10. Habibi, K., Gong H, and Mongeau L. "A robust numerical approach for the prediction of turbofan engine noise." Proceedings of Meetings on Acoustics. Vol. 19. 2013. 11. K. Habibi, L. Mongeau, D. Casalino, P.-T. Lew, “Aeroacoustic Study of Internal Mixing Nozzles with Forced Lobed Mixers using a High-Mach Subsonic Lattice Boltzmann scheme”, AIAA-2014-3313 12. P.-T. Lew et al, K. Habibi and L. Mongeau, “An extended Lattice Boltzmann Methodology for High Subsonic Jet Noise Prediction”, AIAA2014-2755. 13. Casalino, D. and Lele, S. “Lattice-Boltzmann Simulation of Coaxial Jet Noise Generation”, Center for Turbulence Research, Proceedings of the Summer Program 2014. 14. M. Murayama, Y. Yokokawa, T. Imamura, K. Yamamoto, H. Ura, T. Hirai, Numerical investigation on change of airframe noise by flap side-edge shape, AIAA Paper 2013-2067, 2013. 15. L.G.C. Simoes, D.S. Souza, M.A.F. Medeiros, On the small effect of boundary layer thicknesses on slat noise, AIAA Paper 2011-2906, 2011. 16. P.-T. Lew, R. Shock, A. Najafi-Yazdi, L. Mongeau, S. Colavincenzo, R. Lapointe, G.Waller, Noise prediction from a partially closed slat junction, AIAA Paper 2013-2161, 2013. 17. E. Fares, S. Noelting, Unsteady flow simulation of a high-lift configu- ration using a lattice Boltzmann approach, AIAA Paper 2011-869, 2011. 18. B. Koenig, E. Fares, S. Noelting, Lattice-Boltzmann flow simulations for the HILIFTPW-2, AIAA Paper 2014-0911, 2014. 19. Seror, C., Sagaut, P., and Belanger, A., “A Numerical Aeroacoustics Analysis of a Detailed Landing Gear,” AIAA Paper No. 2004-2884, May 2004. 20. G.A. Brès, D.M. Freed, M.Wessels, S. Noelting, F. Pérot, Flow and noise predictions for the tandem cylinder aeroacoustic benchmark, Physics of Fluids 24 (3) (2012) 036101. 28

21. A. Keating, P. Dethioux, R. Satti, S. Noelting, J. Louis, T. Van de Ven, R. Vieito, Computational aeroacoustics validation and analysis of a nose landing gear, AIAA Paper 2009-3154, 2009. 22. Casalino, D., Ribeiro, A. F. P and Fares, E. “Facing rim cavities fluc- tuation modes”, Journal of Sound and Vibration. Vol. 333, pp. 2812–2830, 2014. 23. Casalino, D., Ribeiro, A. F. P., Fares, E., and Noelting, S., “Lattice- Boltzmann Aeroacoustic Analysis of the LAGOON Landing Gear Configuration,” AIAA Journal , Vol. 52, No. 6, 2014, pp. 1232–1248. 24. Bouvy, Q., Rougier, T., Ghouali, A., Casalino, D., Appelbaum, J., and Kleinclaus, C. “Design of Quieter Landing Gears through Lattice-Boltzmann CFD Simulations”, AIAA 2015-3259. 29 Powered by TCPDF (www.tcpdf.org)

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв

10.66

Георгий Чащин

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

6315.12

Игорь Маслеников

и хорошего настроения

6315.12

Игорь Маслеников

удачи

6315.12

Игорь Маслеников

успехов в конкурсе

6315.12

Игорь Маслеников

Наверное было затрачено много времени и труда на работу