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

На правах рукописи

Шершнёв Антон Алексеевич

МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ ГАЗА В ПЕРЕХОДНОМ РЕЖИМЕ НА ОСНОВЕ РЕШЕНИЯ МОДЕЛЬНЫХ КИНЕТИЧЕСКИХ УРАВНЕНИЙ

05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

г 5 ш

Новосибирск - 2013

005531671

Работа выполнена в Федеральном государственном бюджетном учреждении науки Институте теоретической и прикладной механики им. С. А. Христиановича Сибирского отделения Российской академии наук

Научный руководитель: кандидат физико-математических наук

Кудрявцев Алексей Николаевич

Официальные оппоненты: Григорьев Юрий Николаевич,

доктор физико-математических наук, профессор,

ИВТ СО РАН, главный научный сотрудник

Сковородко Петр Алексеевич, кандидат физико-математических наук, ИТ СО РАН, ведущий научный сотрудник

Ведущая организация: Федеральное государственное бюджетное

учреждение науки Институт прикладной математики им. М. В. Келдыша Российской академии наук, г. Москва

Защита состоится 23 октября 2013 г. в 14:00 на заседании диссертационного совета ДМ003.046.01 на базе Федерального государственного бюджетного учреждения науки Института вычислительных технологий Сибирского отделения Российской академии наук по адресу: 630090, Новосибирск, пр. Ак. Лаврентьева 6, конференц-зал ИВТ СО РАН.

С диссертацией можно ознакомиться в библиотеке ИВТ СО РАН. Автореферат разослан « » июля 2013 г.

Ученый секретарь диссертационного совета, к.ф.-м.н., доцент

А. С. Лебедев

Общая характеристика работы

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

Обычно в микроустройствах встречаются течения в околоконтинуальном или переходном режимах, что подразумевает учет эффектов разреженности, поскольку средняя длина свободного пробега молекул А в этом случае не может считаться пренебрежимо малой по сравнению с характерным размером задачи L. Хорошо известно, что континуальный подход, основанный на решении уравнений Навье-Стокса с граничными условиями прилипания на твердой стенке, применим для случаев, когда параметр Kn = A/L, называемый числом Кнудсена, удовлетворяет условию Кп < 0,01. Использование условий скольжения и скачка температуры позволяет расширить область применимости до Кп = 0,1. Однако в переходном режиме (0,1 < Кп < 10) уравнения Навье-Стокса не способны описать течения с достаточной точностью даже при использовании граничных условий скольжения.

Для описания течений в переходном режиме активно развиваются как континуальные, так и кинетические подходы. Континуальные подходы основаны на системах уравнений, являющихся приближениями более высоких порядков по числу Кнудсена, чем уравнения Навье-Стокса, таких как «дополненные» (augmented) уравнения Барнетта и уравнения БГК-Барнетта (Zhong, 1993, Agarwal, 2001), стабилизированные моментные уравнения (система R13, Struchtrup, Torrilhon, 2003, Иванов, Крюков, Тимохин и др., 2012), на некоторых других подходах (квазигазодинамические уравнения, Елизарова, Четверушкин, 1984, Шеретов, 1996).

Несомненно, наиболее универсальным решением является использование кинетического подхода, который основан на представлении газа в виде ансамбля частиц, чья функция распределения по скоростям подчиняется уравнению Больц-мана. В настоящее время наиболее эффективным численным методом решения уравнения Больцмана является метод прямого статистического моделирования (ПСМ), интенсивно развивающийся в течение почти полувека (Bird, 1965, Бело-церковский, Яницкий, 1975, Иванов, Рогазинский, 1988). Трудности с применением ПСМ для моделирования медленных и нестационарных течений связаны с его стохастическим характером, приводящим к статистическому разбросу данных.

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

скоростей (Григорьев, Михалицын, 1983, Filbet, Russo, 2003), консервативные схемы для аппроксимации интеграла столкновений (Черемисин, 1998), использование дивергентной формы уравнения Больцмана (Савельев, Nanbu, 2002, Малков, Иванов, 2011), прямое решение уравнения Больцмана все еще требует очень больших вычислительных ресурсов.

Более практичным в настоящее время представляется использование подхода, основанного на детерминистическом решении модельных кинетических уравнений, в которых интеграл столкновений заменен членом релаксационного типа. Впервые уравнение подобного типа было рассмотрено Бхатнагаром, Гроссом и Круком (уравнение БГК, 1954). Уравнения Навье-Стокса могут быть выведены из уравнения БГК так же, как и из уравнения Больцмана. Однако это приводит значению числа Прандтля Рг равному 1, в то время как правильное значение для одноатомного газа Рг = 2/3. Для того чтобы устранить этот недостаток, было предложено несколько модификаций уравнения БГК, таких как эллипсоидальная статистическая модель (ЭС модель, Holway, 1966) и модель Шахова (1968).

Для численного решения модельных кинетических уравнений можно ввести конечно-разностную сетку в пространстве скоростей (метод дискретных ординат) и в физическом пространстве. В ранних работах для аппроксимации производных в физическом пространстве использовались простые конечно-разностные схемы, обычно первого порядка точности, в частности метод характеристик (Шахов, 1974). В последние годы для этой цели было предложено использовать современные схемы сквозного счета (Yang, Huang, 1995, Mieussens, 2000, Титарев, 2007). Были разработаны алгоритмы, использующие неструктурированные сетки, развиты эффективные неявные схемы для моделирования стационарных течений

(Титарев, 2009, 2010).

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

Цели диссертационной работы. Целями данной работы являются:

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

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

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

и микросопле.

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

1. Разработан алгоритм прямого конечно-разностного решения модельных кинетических уравнений на основе метода дискретных ординат и WENO схем высокого порядка точности.

2. На основе разработанного алгоритма создан программный код для расчетов на многопроцессорных ЭВМ, параллелизованный с помощью геометрической декомпозиции расчетной области и протокола MPI; с помощью численных экспериментов получены оценки точности квадратурных формул в методе дискретных ординат, позволяющие выбрать параметры дискретизации, необходимые для достижения требуемой точности вычислений.

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

4. Путем сравнения с другими походами (решением уравнений Навье-Стокса, методом ПСМ) получены данные о точности и границах применимости модельных кинетических уравнений при моделировании течений в переходном режиме.

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

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

Представление работы. Результаты докладывались и обсуждались на следующих научных конференциях и семинарах:

— Международная конференция по методам аэрофизических исследований (ICMAR), Новосибирск, 2008, 2011;

— Международный симпозиум по динамике разреженного газа (RGD), Киото, Япония, 2008, Пасифик Гроув, США, 2010;

— Всероссийская конференция молодых ученых «Проблемы механики: теория, эксперимент и новые технологии», Новосибирск, 2009;

_ Европейская конференция по микротечениям (GASMEMS), Эйндховен, Нидерланды, 2009.

— Международная конференция ASME по тепло- и массообмену на микро- и наномасштабах, Шанхай, Китай, 2009;

— Всероссийский семинар «Фундаментальные основы МЭМС- и нанотехноло-гий», Новосибирск, 2010;

_ Всероссийская школа-конференция молодых ученых «Актуальные вопросы

теплофизики и физической гидрогазодинамики», Новосибирск, 2010;

— Международный симпозиум по ударным волнам (ISSW), Манчестер, Великобритания, 2011,

_ Международная конференция по неравновесным процессам в соплах и струях (NPNJ), Алушта, 2012,

а также на научном семинаре ИТПМ СО РАН под руководством академика В.М. Фомина, объединенном научном семинаре ИВТ СО РАН, кафедры математического моделирования НГУ и кафедры вычислительных технологий НГ-ТУ «Информационно-вычислительные технологии» под руководством академика Ю.И. Шокина и проф. В.М. Ковени и на семинаре отдела разреженных газов ИТ СО РАН под руководством академика А.К. Реброва.

Публикации. По материалам диссертации опубликовано 14 работ, в том числе 3 статьи в рецензируемых журналах, рекомендованных ВАК для для представления основных результатов диссертации. Список работ приведен в конце автореферата.

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

Структура и объем диссертации. Диссертационная работа состоит из введения, четырех глав, заключения и списка литературы из 129 наименований. Общий объем диссертации составляет 119 страниц, включая 57 рисунков и 3 таблицы.

Основное содержание работы

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

Глава 1 носит методический характер. В ней представлен общий вид модельного кинетического уравнения релаксационного типа описана связь между функцией распределения и макроскопическими величинами. Модельное кинетическое уравнение для функции распределения / может быть записано в виде

где х = (дг.^г) — вектор пространственных координат, V = (их, иу, vz) = , 1>2, г>з) — вектор молекулярных скоростей, — равновесная функция распределения, а V — частота столкновений, величина обратная к среднему времени между молекулярными столкновениями. В § 1.1 рассмотрен переход к безразмерным переменным, выписаны формулы для вычисления безразмерных макровеличин. В § 1.2 введены укороченные функции распределения, позволяющие уменьшить число независимых переменных и сократить вычислительные затраты при решении одно- и двумерных задач, приведены системы уравнений для укороченных функций распределения и формулы для вычисления газодинамических параметров по укороченным функциям распределения. В § 1.3 описаны три рассматриваемые нами модели столкновительного члена: модель БГК, модель Шахова и ЭС модель. В оригинальной модели БГК равновесная функция распределения является локально-максвелловской функцией

_ _ Р Г с2 }

1 =(2^Г)3/2ехрГ2Лг/- с = (сх,су,се) = у-и. (2)

Здесь р это плотность, и — средняя скорость потока, Т — температура, а И это удельная газовая постоянная, равная отношению универсальной газовой постоянной <% к молярной массе газа. В ЭС модели равновесная функция ищется как наиболее вероятное распределение, у которого заданы моменты до второго порядка включительно. Результат представляется анизотропным гауссианом вида

fN _ гЕ _ Р 1 ^/2^АехР

V \

А = (аи), В = (Ьи) = А~\ (3)

ЯТ _ 1 - Рг рп

где элементы матрицы аи = —8и - , Рг — это число Прандтля, а ри —

тензор давления. Равновесная функция распределения в модели Шахова определяется следующим образом:

14

5 РТ

где р — это давление, ад — вектор теплового потока.

г-50

Рис. 1: Моделирование затухания акустической волны с помощью различных численных схем. На рисунке изображено распределение давления.

Рис. 2: Зависимость обратной толщины ударной волны от числа Маха.

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

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

Раздел 2.1 посвящен дискретизации в пространстве скоростей и методу дискретных ординат. В § 2.1.1 и § 2.1.2 описано вычисление интегралов с помощью метода Симпсона и квадратур Гаусса-Эрмита соответственно. В обоих случаях с помощью численных экспериментов для функций максвелловского типа получены оценки точности интегрирования, предложены процедуры выбора параметров сетки, обеспечивающие вычисление газодинамических величин с требуемой точностью. В § 2.1.3 проведено сравнение эффективности квадратур на примере задачи об обтекании тела произвольной формы с теплоизолированной поверхностью. Показано, что квадратуры Гаусса-Эрмита существенно более эффективны в случае дозвуковых и умеренно сверхзвуковых течений, в то время как для гиперзвуковых течений предпочтительно использовать составное правило Симпсона. Раздел § 2.2 посвящен дискретизации уравнений в координатном пространстве. В первой части раздела описаны используемые конечно-разностные схемы, включая противопоточную схему первого порядка, стандартную ТУИ схему второго порядка точности с ограничителем наклона пипшос!, \VENO схемы 3-го порядка и 5-го порядка и сохраняющую монотонность решения МР5-схему. Переход к общим криволинейным координатам обсуждается в § 2.2.1.

В § 2.3 описано интегрирование по времени системы ОДУ, полученной в результате дискретизации скоростного и координатного пространств. Переход на следующий временной слой проводится с помощью т. н. ТУБ схем Рунге-Кутты 2-го либо 3-го порядка, гарантирующих сохранение монотонности решения. Так-

же обсуждается выбор временного шага. Раздел 2.4 посвящен постановке граничных условий. Обсуждаются условия на входных и выходных границах, твердых стенках (диффузное либо зеркальное отражение), описаны характеристические граничные условия, используемые на дозвуковых границах. В § 2.5 обсуждаются особенности программной реализации численного метода. Расчетный код распараллелен, параллелизация основана на геометрической декомпозиции расчетной области. В серии расчетов измерено параллельное ускорение, показано, что при достаточно большом числе процессоров эффективность параллелизации достигает 86%. Это позволило проводить достаточно требовательные к вычислительным ресурсам расчеты. В частности, при моделировании течения около конечной плоской пластины (§ 4.2) при числе Маха М = 10 общее число узлов в фазовом пространстве достигало 1,86 х 109. В § 2.6 кратко описываются методы решения уравнений Навье-Стокса и уравнения Больцмана; полученные с их помощью данные сравниваются в диссертации с результатами расчетов, выполненных на основе модельных кинетических уравнений.

Глава 3 посвящена численным исследованиям одномерных задач динамики разреженного газа.

На примере простых задач, таких, как задача об одномерной ударной трубе и задача об отражении ударной волны от стенки, проведена верификация численного метода. В § 3.3 рассматривается задача о распространении акустической волны, показано хорошее согласие численного решения модельных кинетических уравнений с аналитическим решением уравнений Навье-Стокса. Кроме того, на примере данной задачи проведено сравнение точности \VENO схемы 5-го порядка с простейшей схемой 1-го порядка и стандартной ТУБ схемой 2-го порядка. Как видно из рис. 1, схемы низких порядков даже на существенно более мелких сетках дают по сравнению с \VENO схемой 5-го порядка худшие результаты в силу большой численной вязкости.

В § 3.4 численно промоделирована теплопередача между параллельными пластинами, получено удовлетворительное согласие с экспериментальными данными и расчетами, выполненными на основе уравнения Больцмана. При этом наблюдается очень хорошее согласие с результатами итерационного решения стационарного модельного кинетического уравнения методом характеристик (Сгаиг, РоНкагроу, 2009). Проверена точность выполнения закона сохранения массы, показано, что изменение полной массы в процессе расчета составило 0,003%.

В § 3.5 рассмотрена задача о структуре ударной волны в аргоне. Расчеты проведены для широкого диапазона чисел Маха, результаты сравниваются с экспериментальными данными Альсмейера и решениями уравнений Навье-Стокса. На рис. 2 представлена т. н. обратная толщина ударной волны, которая определяется как Х/д, где <5 это толщина ударной волны, вычисленная по максимальному наклону профиля плотности. Из трех модельных кинетических уравнений наилучшее совпадение с измерениями показывает ЭС модель. Модель БГК пред-

сказывает несколько заниженную толщину волны, модель Шахова — предсказывает завышенное значение толщины, особенно при больших числах Маха. Уравнения Навье-Стокса, как и следовало ожидать, при высоких числах Маха дают значительно заниженное значение толщины ударной волны.

Кроме того, в расчетах получены данные о распределении макроскопических величин. Продемонстрировано полное совпадение решений модельных кинетических уравнений и аналитического решения Йена (Phys. Fluids 1966). Сравнение функций распределения внутри ударной волны показывает, что они имеют двух-пиковый вид, качественно схожий с бимодальным распределением Мотт-Смита.

Глава 4 посвящена двумерным задачам динамики разреженного газа.

В § 4.1 моделируется конвекция изэнтропического вихря в области с периодическими по обеим координатам граничными условиями. На примере данной задачи изучается, с какой точностью в численном методе выполняются законы сохранения. В процессе расчета вычислялось изменение массы, импульса и энергии в расчетной области. Процентное отношение минимального значения к максимальному для массы составило 99,98%, для импульса 99,93% и для энергии 99,97%.

В § 4.2 было выполнено численное моделирование разреженного течения одноатомного газа около конечной плоской пластины на основе трех различных подходов: уравнений Навье-Стокса с граничными условиями скольжения и скачка температуры, метода ПСМ и конечно-разностного решения модельных кинетических уравнений. Вычисления проведены для числа Кнудсена Кп = 0,01 и чисел Маха М = 2 и М = 10. Показано, что при М = 2 все три подхода дают близкие результаты. Однако в случае М = 10 в обширной области над пластиной отсутствует равновесие между различными поступательными степенями свободы. Это ведет к тому, что уравнения Навье-Стокса не позволяют корректно предсказать

Рис. 5: Отношение продольной температуры Тх к полной температуре Т для чисел Маха Моо = 10 (сверху) и Моо = 2 (снизу).

характеристики течения даже при использовании на поверхности пластины условий скольжения. В частности, на рис. 3 представлены профили давления в сечении поперек пластины, а на рис. 4 — коэффициент трения на поверхности пластины, демонстрирующие существенные различия континуального и кинетического решений. Эти различия объясняются нарушением равновесия между различными поступательными степенями свободы молекул в разреженном гиперзвуковом течении (см. рис. 5).

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

Раздел § 4.3 посвящен численному моделированию гиперзвукового течения разреженного газа у цилиндрически затупленной толстой пластины под нулевым углом атаки с использованием трех различных подходов: метода ПСМ, ЭС модели и уравнений Навье-Стокса. Показано, что хотя при числе Кнудсена по радиусу затупления Кпд = 0,1 уравнения Навье-Стокса с граничными условиями скольжения и температурного скачка на поверхности тела позволяют достаточно успешно описывать начальные эффекты разреженности для вязкого течения за ударной волной, при Кпд = 0,5, когда имеет место существенная анизотропия функции распределения, необходимо применение кинетического подхода, позволяющего получить существенно более точные результаты (рис. 6, 7).

В § 4.4 промоделирован процесс распространения в плоском микроканале ударной волны, возникающей после разрыва диафрагмы, которая в начальный момент разделяет две области с газами с различными давлениями. Проведено сравнение трех подходов: модельных кинетических уравнений, метода ПСМ и уравнений Навье-Стокса с граничными условиями скольжения. Вычисления проведены для чисел Кнудсена Кп = 0,05 и Кп = 0,5, вычисленных по параметрам в области с низким давлением. Показано, что при числе Кнудсена Кп = 0,05 все три подхода дают близкие результаты: хорошее согласие наблюдается как в скорости распространения ударной волны, так и в распределениях макроскопических величин вдоль и

Рис. 6: Коэффициент трения на цилиндрически затупленной пластине, Кпд = 0,5, М = 5.

Рис. 7: Сверхзвуковое обтекание цилиндрически затупленной пластины. Профиль плотности в сечении х = 10 для Кпд = 0,5.

поперек канала. Однако при числе Кундсена Кп = 0,5 существенную роль начинают играть эффекты разреженности. В данном случае хорошо видны преимущества подхода, основанного на использовании модельных кинетических уравнений. По сравнению с уравнениями Навье-Стокса, при Кп = 0,5 ЭС модель более точно описывает распределения продольной скорости поперек канала (рис. 9) и профиль температуры вдоль центральной линии сопла (рис. 10). В то же время, в отличие от данных ПСМ, расчеты на основе ЭС модели не страдают от статистического шума, который особенно заметен в этой нестационарной задаче.

Рис. 8: Распространение ударной волны в микроканале. Изолинии температуры в момент времени { = 80 не, Кп = 0.05 в расчетах на основе метода ПСМ (сверху) и ЭС модели (снизу).

В § 4.5.1 проведены расчеты стационарного течения в плоском клиновидном микросопле со скругленной критической частью и расчетным числом Маха М — 4 на выходе. Число Рейнольдса Ке*, вычисленное по параметрам в критическом сечении, менялось в диапазоне 50 ~ 350. Проведено сравнение полученных результатов с решениями уравнений Навье-Стокса. При Яе» = 350 в целом наблюдается хорошее согласие кинетического и континуального подходов. Наибольшие различия наблюдаются вблизи стенки сопла в профилях давления. При уменьшении числа Рейнольдса расхождения в результатах усиливаются. Заметно опреде-

Рис. 9: Распространение ударной волны в микроканале. Профиль скорости в сечении х/Н = 22,5 в момент времени / = 80 не, Кп = 0,5.

Рис. 10: Распространение ударной волны в микроканале. Распределения температуры вдоль линии симметрии, Кп = 0,5.

ленное отличие температур как на оси сопла, так и, особенно ярко выраженное, вблизи стенки. Это связано с тем, что в пристеночной области функция распределения далека от равновесной из-за влияния отраженных от стенки молекул. Однако сравнение распределений величин вдоль линии симметрии показывает, что кинетические и континуальные расчеты дают очень близкие результаты даже при числе Рейнольдса Ие* = 50. В то же время распределения при различных числах Рейнольдса заметно отличаются в сверхзвуковой части сопла, что свидетельствует о влиянии толстого пограничного слоя на ядро потока.

В § 4.5.2 промоделирован процесс запуска плоского клиновидного микросопла ударной волной. В расчете ударная волна возникала при распаде разрыва в форкамере сопла. Начальный перепад давления был таким, что число Маха возникавшей ударной волны равнялось 1,5. Числа Кнудсена для состояний газа справа и слева от диафрагмы равнялись Кпд = 0,0214 и Кп/. = 0,0026 соответственно. На рис. 11 показаны мгновенные поля течения, полученные в расчетах на основе ЭС модели и уравнений Навье-Стокса, на момент времени, когда первичная ударная волна достигла выхода из сопла. Развивающийся за ударной волной пограничный слой ближе к горлу сопла отрывается от стенки. На рис. 12 показаны распределения вдоль центральной линии сопла, полученные с помощью трех различных подходов. Сравнение с невязким расчетом показывает, что вязкость существенно меняет характер течения. Если интенсивность первичной ударной волны просто несколько уменьшается, то вторичная ударная волна полностью исчезает. Вместо нее наблюдается достаточно плавный подъем давления, начинающийся примерно щтлх/к = 3. В целом, в данной задаче результаты кинетического и континуального расчетов достаточно близки друг к другу.

Рис. 11: Стартовый процесс в клиновидном сопле. Поле числа Маха, полученное помощью ЭС модели и уравнений Навье-Стокса.

Рис. 12: Распределения давления вдоль центральной линии сопла в момент времени г = 11,93, полученные с помощью ЭС модели, уравнений Навье-Стокса для Ие» = 50 и уравнений Эйлера.

В заключении изложены основные результаты работы:

1. На основе метода дискретных ординат и современных методов сквозного счета (WENO схем 3-го и 5-го порядков, МР схемы 5-го порядка) разработан алгоритм прямого численного решения модельных кинетических уравнений (уравнения БГК, модели Шахова и ЭС модели). Получены оценки точности интегрирования в методе дискретных ординат для квадратур Симп-сона и Гаусса-Эрмита. Даны рекомендации по выбору параметров дискретизации, обеспечивающие нужную точность вычисления газодинамических величин.

2. Создан программный код, решающий кинетические уравнения релаксационного типа с тремя моделями столкновительного члена. Проведена верификация программы, на тестовых примерах показана высокая точность численного метода. Программа распараллелена с помощью геометрической декомпозиции расчетной области и протокола MPI. Продемонстрирована достаточно высокая эффективность параллелизации 86% от идеального ускорения при расчете на 64 процессорах).

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

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

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

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

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

Автор выражает глубокую признательность своему научному руководителю А. Н. Кудрявцеву за всестороннюю поддержку. Также автор хотел бы поблагодарить Е. А. Бондаря за предоставленные результаты ПСМ расчетов, Д. В. Хотя-новского за предоставленные результаты моделирования затупленной пластины с помощью уравнений Навье-Стокса и профессора М. С. Иванова за неоднократные полезные обсуждения полученных результатов.

Список работ по теме диссертации

Публикации в периодических изданиях, рекомендованных ВАК

1. Иванов, М. С. Эффекты разреженности при обтекании затупленной передней кромки гиперзвуковым потоком / М. С. Иванов,

Д. В. Хотяновский, А. А. Шершнёв, А. Н. Кудрявцев, А. А. Шевырин, Ш. Ёнемура, Е. А. Бондарь // Теплофизика и аэромеханика. — 2011. — Т. 18, №4, —С. 543-554.

2. Шершнёв, А. А. Численное моделирование сверхзвукового течения газа около плоской пластины на основе кинетических и континуальных моделей / А. А. Шершнёв, А. Н. Кудрявцев, Е. А. Бондарь // Вычислительные Технологии. — 2011. — Т.16, № 6. — С. 93-104.

3. Kudryavtsev, А. N. А numerical method for Simulation of microflows by solving directly kinetic equations with WENO schemes / A. N. Kudryavtsev,

А. A. Shershnev // Journal of Scientific Computing. — 2013. — DOI: 10.1007/s 10915-013-9694-z, 32 p.

Публикации в трудах международных конференций

4. Kudryavtsev, А. N. Numerical simulation of the shock-wave structure with different kinetic and continuum models / A. N. Kudryavtsev, A. A. Shershnev, M. S. Ivanov // Proc. of 14th Int. Conf. on Methods of Aerophysical Research. — 2008. — 1 CDROM. — paper № 14, 7 p.

5. Kudryavtsev, A. N. Comparison of different kinetic and continuum models applied to the shock-wave structure problem / A. N. Kudryavtsev,

A. A. Shershnev, M. S. Ivanov // Proc. of 26th Int. Symp. on Rarefied Gas Dynamics. — Mellville, New York, 2009. — P. 507-512.

6. Shershnev, A. A. Study of rarefaction and non-equilibrium effects using model kinetic equations / A. A. Shershnev, A. N. Kudryavtsev, I. A. Graour // Proc. of 1st GASMEMS Workshop. — 2009. — 1 CDROM. — paper № 10, 8 p.

7. Kudryavtsev, A. N. Numerical simulation of gas microflows by solving relaxation-type kinetic equations / A. N. Kudryavtsev, A. A. Shershnev,

M. S. Ivanov // Proc. of the ASME 2nd Micro/Nanoscale Heat & Mass Transfer Int. Conf. — 2009. — paper № 18520, 10 p.

8. Kudryavtsev, A. N. A study of the finite flat plate problem using various kinetic and continuum models. / A. N. Kudryavtsev, A. A. Shershnev, Ye. A. Bondar // Proc. of 27th Int. Symp. on Rarefied Gas Dynamics. — Mellville, New York, 2011. — V. 2. — P. 934-940.

9. Ivanov, Mikhail. Rarefaction and non-equilibrium effects in hypersonic flows about leading edges of small bluntness / Mikhail Ivanov, Dmitry Khotyanovsky, Alexey Kudryavtsev, Anton Shershnev, Yevgeniy Bondar, Shigeru Yonemura // Proc. of the 27th Int. Symp. on Rarefied Gas Dynamics. — Mellville, New York, 2011. — V. 2. — P. 1295-1301.

10. Shershnev, A. A. Kinetic approach to numerical simulation of gas flows in the transitional regime / A. A. Shershnev, A. N. Kudryavtsev // Proc. of 15th Int. Conf. on Methods of Aerophysical Research. — 2010. — 1 CDROM. — 9 p.

11. Ivanov, M. S. Numerical study of rarefied hypersonic flows about blunted leading edges / M. S. Ivanov, S. Yonemura, D. Khotyanovsky, A. Kudryavtsev, A. Shershnev, Ye. Bondar // Proc. of 15th Int. Conf. on Methods of Aerophysical Research. — 2010. — 1 CDROM.

12. Bondar, Ye. A. Numerical study of hypersonic rarefied flows about leading edges of small bluntness / Ye. A. Bondar, A. A. Shershnev, A. N. Kudryavtsev, D. V. Khotyanovsky, S. Yonemura, M. S. Ivanov // Proc. of 28th Int. Symp. on Shock Waves. — 2011. — 1 USB flash drive — paper № 2561, 6 p.

13. Кудрявцев, A.H. Численное моделирование нестационарных течений разреженного газа путем прямого конечно-разностного решения кинетических уравнений / А. Н. Кудрявцев, А. А. Шершнёв // Материалы IX Международной конференции по неравновесным процессам в соплах и струях. — Москва: изд-во МАИ, 2012. — С. 226-228.

Публикации в тезисах всероссийских конференций

14. Шершнёв, А. А. Метод прямого численного решения кинетических уравнений для моделирования течений газа в переходном режиме /

А. А. Шершнёв // Труды 7-ой всероссийской конференции молодых ученых «Проблемы механики: теория, эксперимент и новые технологии». — 2009. — С. 241-243.

15. Шершнёв, A.A. Метод численного моделирования микротечений на основе модельных кинетических уравнений релаксационного типа /

А. А. Шершнёв, А. Н. Кудрявцев // XI Всероссийская школа-конференция молодых ученых «Актуальные вопросы теплофизики и физической гидрогазодинамики». — 2010. — С. 102.

Подписано в печать 05.07.2013 Тираж 140 экз., Заказ № 037 Отпечатано «Документ-Сервис» 630090, Новосибирск, ул. Институтская, 4/1, тел. (383) 335-66-00

Федеральное государственное бюджетное учреждение науки Институт теоретической и прикладной механики

им. С. А. Христиановича Сибирского отделения Российской академии наук

На правах рукописи

Шершнёв Антон Алексеевич

МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ ГАЗА В ПЕРЕХОДНОМ РЕЖИМЕ НА ОСНОВЕ РЕШЕНИЯ МОДЕЛЬНЫХ КИНЕТИЧЕСКИХ УРАВНЕНИЙ

^ 05.13.18 - Математическое моделирование, численные методы

^ и комплексы программ

см „

т— Диссертация на соискание ученой степени

СО £

^ кандидата физико-математических наук

О ^ СМ 2

^^ Научный руководитель

® кандидат физико-математических наук

Кудрявцев Алексей Николаевич

Новосибирск - 2013

Содержание

Введение............................................................................4

Глава 1. Основные уравнения..................................................16

1.1. Кинетическое уравнение для функции распределения................16

1.2. Приведение к безразмерному виду......................................17

1.3. Укороченные функции распределения ................................18

1.4. Модели столкновительного члена...................20

1.4.1. Уравнение БГК.........................20

1.4.2. Эллипсоидальная статистическая модель..........21

1.4.3. Модель Шахова........................22

Глава 2. Численный метод и программная реализация..........24

2.1. Дискретизация в пространстве скоростей...............24

2.1.1. Метод Симпсона........................25

2.1.2. Квадратурная формула Гаусса-Эрмита...........27

2.1.3. Сравнение двух методов интегрирования..........29

2.1.4. Дискретизация источникового члена.............30

2.2. Дискретизация в координатном пространстве........................31

2.2.1. Обобщенные криволинейные координаты..........39

2.3. Интегрирование по времени......................40

2.4. Граничные условия...........................42

2.4.1. Свободная граница.......................42

2.4.2. Диффузное отражение.....................43

2.4.3. Зеркальное отражение.....................44

2.4.4. Условия на дозвуковой входной границе ..........44

2.5. Параллельная реализация алгоритма.................46

2.6. Уравнения Навье-Стокса и ПСМ...................49

Глава 3. Моделирование одномерных задач динамики разреженного газа 51

3.1. Задача об одномерной ударной трубе..................................51

3.2. Отражение ударной волны от стенки.................53

3.3. Распространение акустической волны.................54

3.4. Теплопередача между пластинами...................57

3.5. Структура ударной волны.......................60

Глава 4. Моделирование двумерных задач динамики разреженного газа 66

4.1. Конвекция изэнтропического вихря..................66

4.2. Конечная плоская пластина......................68

4.3. Цилиндрически затупленная пластина.................77

4.4. Распространение ударной волны в микроканале...........84

4.5. Микросопло...............................93

4.5.1. Стационарное течение в сопле................94

4.5.2. Стартовый процесс в сопле..................99

Заключение....................................103

Литература ....................................105

Введение

Быстрое развитие микротехнологий в течение двух последних десятилетий создало предпосылки к разработке и широкому распространению микро- и на-ноэлектромеханических устройств (NEMS и MEMS). Широкий спектр различных микромашин внедрен и успешно применяется в электронике, космонавтике, в различных отраслях промышленности, в том числе и медицинской. В силу малых размеров устройств, характерные свойства газовых микротечений в них могут существенно отличаться от течений в макроустройствах. Детальное понимание природы течений в микросистемах является одной из наиболее актуальных задач современной механики жидкостей и газов [1,2].

Модель дискретных частиц

Динамика сплошной среды

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

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

Уравнения Эйлера

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

Приближения высших порядков

|-^-1-1-1-1-1-

О --0.01 0.1 1 10 100 -»-00

Невязкий Локальное число Кнудсена кп Свободно-

предел молекулярный предел

Рис. 1: Области применимости математических моделей описания течения газа.

Обычно в микроустройствах встречаются течения в околоконтинуальном или переходном режимах, что подразумевает учет эффектов разреженности, поскольку средняя длина свободного пробега молекул Я в этом случае не может считаться пренебрежимо малой по сравнению с характерным размером задачи Ь. Хорошо известно [3], что континуальный подход, основанный на решении уравнений Навье-Стокса с граничными условиями прилипания на твердой стенке, применим для случаев, когда параметр Кп = Х/Ь, называемый числом Кнудсена, удовлетворяет условию Кп < 0.01. Использование условий скольжения и скачка температуры позволяет расширить область применимости до Кп = 0.1. Численное моделирование течений в микроустройствах в данном режиме со скольжением подробно обсуждается в [4]. Однако в переходном режиме (0.1 < Кп <10) уравне-

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

Одним из возможных способов решения данной проблемы может быть использование континуальных (гидродинамических) уравнений, являющихся приближениями более высоких порядков по числу Кнудсена, чем уравнения На-вье-Стокса. В англоязычной литературе подобные уравнения известны под общим названием «Extended Hydrodynamic Equations». Известно [5], что уравнения Навье-Стокса могут быть выведены из уравнения Больцмана, описывающего эволюцию функции распределения молекул газа по скорости, путем разложения ее по степеням Кп {разложение Чепмена-Энскога) до величин первого порядка малости. Если же рассмотреть члены до второго порядка малости 0(Кп2), то такое разложение приводит к так называемым уравнениям Барнетта, включающим пространственные производные более высоких (выше второго) порядков. Было, однако, показано, что уравнения Барнетта неустойчивы к коротковолновым возмущениям и, при больших числах Кнудсена, могут приводить к нарушению второго начала термодинамики.

Альтернативным разложению Чепмена-Энскога подходом является момент-ный метод, при котором выводятся уравнения для моментов функции распределения. Используя определенные предположения, такую систему уравнений можно замкнуть. Наиболее известной системой моментных уравнений является 13-моментная система уравнений Грэда [5], включающая, кроме уравнений для обычных гидродинамических величин, еще уравнения для компонент тензора вязких напряжений и вектора теплового потока. К сожалению, попытка использования уравнений Грэда также наталкивается на значительные трудности. Эта система уравнений является гиперболической, в результате в ней появляются физически бессмысленные разрывные решения, в частности при расчете внутренней структуры ударных волн.

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

годы были предложены методы регуляризации данных уравнений, позволяющие избавиться от тех нефизических следствий, к которым приводят оригинальные уравнения. Так, уравнения Барнетта можно стабилизировать и согласовать со вторым началом термодинамики, если ввести в них некоторые дополнительные члены из следующего, так называемого супербарнеттовского приближения. Полученные подобным образом «дополненные» (augmented) уравнения Барнетта и т. н. уравнения БГК-Барнетта были довольно успешно использованы для описания некоторых разреженных течений [6, 7].

Подобные методы были также использованы для регуляризации уравнений Трэда, что привело к системе уравнений, известных как уравнения R13 [8]. Было продемонстрировано [9], что полученные на их основе решения согласуются с результатами решения уравнения Больцмана лучше решений уравнений На-вье-Стокса.

Еще одним оригинальным подходом к континуальному описанию разреженных течений является использование так называемых квазигазодинамических уравнений (КГД) [10, 11], которые, в отличие от уравнений Навье-Стокса, выводятся не для мгновенных пространственных средних по малому объему газа, а для величин, полученных пространственно-временным осреднением. То есть, к обычному определению газодинамических величин добавляется сглаживание по некоторому малому временному интервалу. Система КГД уравнений отличается от уравнений Навье-Стокса наличием дополнительных диссипативных членов, имеющих порядок малости (9(Кп2). Вклад этих поправок мал для стационарных и квазистационарных течений при малых числах Кнудсена, но становится значительным для сильно нестационарных течений и течений с числом Кнудсена Kn ~ 1. Как показывают результаты расчетов [12], использование КГД уравнений также позволяет расширить область применения континуального подхода.

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

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

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

В настоящее время наиболее эффективным численным методом решения уравнения Больцмана является метод прямого статистического моделирования (ПСМ) [13]. Являясь мощным инструментом для моделирования стационарных сверх- и гиперзвуковых течений [14], метод ПСМ становится существенно более дорогим в вычислительном плане при расчете медленных дозвуковых течений, когда скорость движения газа много меньше средней скорости молекул. Еще менее эффективен данный метод при моделировании нестационарных течений, для которых статистические ошибки не могут быть уменьшены осреднением по достаточно большому временному интервалу.

От этих трудностей свободен подход, основанный на прямом детерминистическом решении кинетических уравнений для функции распределения. Подобные методы существуют достаточно давно и продолжают активно развиваться [15]. В качестве важных достижений в этой области можно указать применение спектральных методов в пространстве скоростей [16-18], консервативные схемы для аппроксимации интеграла столкновений [19-21], а также использование дивергентной формы уравнения Больцмана [22, 23].

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

Более практичным представляется использование подхода, основанного на детерминистическом решении модельных кинетических уравнений, в которых интеграл столкновений заменен членом релаксационного типа. Впервые уравнение подобного типа было рассмотрено Бхатнагаром, Гроссом и Круком [24] и, независимо, Веландером [25], и общеизвестно как уравнение БГК. Уравнение БГК описывает релаксацию к локально-максвелловскому равновесному распределению. Заметим, что уравнения Навье-Стокса могут выведены из уравнения БГК так же, как и из уравнения Больцмана. Однако это приводит к значению числа Прандт-ля Рг равному 1 [5], в то время как правильное значение для одноатомного газа Рг = 2/3. Для того чтобы устранить этот недостаток, было предложено несколько обобщений оригинального уравнения БГК, таких как эллипсоидальная статистическая модель (ЭС модель) [26] и модель Шахова [27]. Оба уравнения имеют ту же форму, что и исходное уравнение БГК, за исключением того, что равновесная функция распределения отлична от максвелловской. Стоит отметить, что частота столкновений в этих моделях, в отличие от уравнения Больцмана, не зависит от относительной скорости сталкивающихся молекул. Поэтому также были предложены модельные кинетические уравнения с частотой, зависящей от относительной скорости, которые обеспечивают правильное значение числа Прандтля [28, 29].

Использование модельных кинетических уравнений для численного моделирования течений в переходном режиме было подробно рассмотрено Е.М. Шаховым в его хорошо известной книге [30]. В последнее десятилетие вновь возник интерес к этому подходу в связи с появлением такой области применения как моделирование микротечений. Как следствие, в последние годы появилось большое число работ, посвященных модельным кинетическим уравнениям. В [31] для решения модельных кинетических уравнений впервые были использованы современные квазимонотонные схемы сквозного счета. В [32, 33] этот подход был улучшен путем введения процедуры коррекции равновесной функции распределения, которая гарантирует выполнение законов сохранения. Численные эксперименты, проведенные в [33], свидетельствуют о том, что данная процедура позволяет полу-

чить достаточно точные решения на относительно грубых сетках в пространстве скоростей. В [34] с помощью модельных кинетических уравнений были решены некоторые трехмерные задачи. В [35] и [36] были предложены гибридные схемы, основанные на сопряжении модельных кинетических и континуальных уравнений. Метод для решения модельных кинетических уравнений на неструктурированных сетках был разработан в [37]. В [38] подход, основанный на детерминистическом решении модельных кинетических уравнений с помощью схем высокого порядка точности, был применен для исследования скорости генерации энтропии и анализа неравновесности течений разреженного газа. В последние годы В.А. Титаревым с соавторами было опубликовано большое количество работ, посвященных моделированию различных течений разреженного газа на основе модельных кинетических уравнений (в основном, модели Шахова) и TVD схем (см., например, [39-41]).

В настоящей работе для численного решения модельных кинетических уравнений используются так называемые WENO (Weighted Essentially Non-Oscillatory) схемы [42]. WENO схемы, хорошо себя зарекомендовавшие при решении газодинамических уравнений, обеспечивают очень высокий (на гладких решениях — до пятого) порядок аппроксимации, что позволяет уменьшить число точек разностной сетки и повысить экономичность решения сложных многомерных задач. Это свойство становится особенно желательным при решении кинетических уравнений, которые приходится решать в многомерном фазовом пространстве.

Целями диссертационной работы являются:

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

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

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

Научная новизна работы заключается в следующем:

1. Разработан алгоритм прямого конечно-разностного решения модельных кинетических уравнений на основе метода дискретных ординат и WENO схем высокого порядка точности.

2. На основе разработанного алгоритма создан программный код для расчетов на многопроцессорных ЭВМ, параллелизованный с помощью геометрической декомпозиции расчетной области и протокола MPI; с помощью численных экспериментов получены оценки точности квадратурных формул в методе дискретных ординат, позволяющие выбрать параметры дискретизации, необходимые для достижения требуемой точности вычислений.

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

4. Путем сравнения с другими походами (решением уравнений Навье-Стокса, методом ПСМ) получены данные о точности и границах применимости модельных кинетических уравнений при моделировании течений в переходном режиме.

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

и

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

— Международная конференция по методам аэрофизических исслед�