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

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

Шувалов Павел Вадимович

МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ РАЗРЕЖЕННОГО

ГАЗА НА ОСНОВЕ КИНЕТИЧЕСКОГО УРАВНЕНИЯ БОЛЬЦМАНА НА КЛАСТЕРНЫХ И ГРАФИЧЕСКИХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ

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

программ

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

АВТОРЕФЕРАТ

6 ИЮН 2013

Москва-2013

005060801

005060801

Работа выполнена на кафедре моделирования ядерных процессов и технологий Московского физико-технического института (государственного университета)

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

доцент

Клосс Юрий Юрьевич

Официальные оппоненты: доктор физико-математических наук,

старший научный сотрудник Бишаев Александр Михайлович, кафедра высшей математики Московского физико-технического института (государственного университета), доцент

доктор физико-математических наук, профессор Самохин Александр Борисович,

кафедра прикладной математики Московского государственного технического университета радиотехники, электроники и автоматики, заведующий кафедрой Ведущая организация: Институт физики высоких температур

РАН

Защита состоится « г / » ¿^Ж-у_ 2013 г. в 5 * часов на

заседании диссертационного совета Д 212.156.05 при Московском физико-техническом институте (государственном университете) по адресу 141700, Московская обл., г. Долгопрудный, Институтский пер., д. 9, ауд. 903 КПМ.

Автореферат разослан « / ?» ^ ^ д._2013 г.

Ученый секретарь диссертационного совета

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

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

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

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

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

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

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

Апробация работы. Научные результаты были доложены, обсуждены и получили одобрение специалистов на научных конференциях: International Conference on Computing Science, Амстердам, 2010; Международная конференция по неравновесным процессам в соплах и струях, Алушта, 2010 и 2012; International Shock Interaction Symposium, Москва, 2010; Международная конференция по вычислительной механике и современным прикладным программным системам, Алушта, 2011; 54-я и 55-я научные конференции МФТИ, Москва-Долгопрудный, 2011 и 2012.

В рамках данной работы были получены свидетельства на программы для ЭВМ [9,10].

Публикации. Материалы диссертации опубликованы в 10 работах, из них 6 - статьи в изданиях из списка, рекомендованного ВАК РФ [1-6].

Личный вклад автора в работы с соавторами.

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

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

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

Во введении приводятся постановка задач и цели диссертации. Описана научная новизна, научная и практическая ценность.

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

Для описания процессов в разреженном газе

используется кинетическое уравнение Больцмана.

Применение такого подхода дают наиболее

последовательный метод описания физических процессов

течений разреженного газа в сложных системах. Подобный

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

интегро-дифференциальных уравнений, вычисления

интеграла столкновений, применения суперкомпьютерных

6

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

Кинетическое уравнение Больцмана может быть представлено в следующем виде:

J —оо 0 о

где f = fit, X, р,); = f(t, X, р,); /, • = fit, X, Р,');

f' = f it, х, p;') - функции распределения плотности газов; t - время; р- импульс молекул; т - масса частиц; g -относительная скорость при столкновении; b, е -параметры столкновения; р(, ру, р^', р.' - импульсы частиц до и после столкновения.

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

Левая часть уравнения описывает движение частиц в пространстве без столкновения между собой. Интеграл в правой части уравнения Больцмана называется интегралом

7

столкновений Больцмана и обозначается следующим образом

7/=X111 (Л 1 л .

у' о о

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

В большинстве расчетов для описания взаимодействия молекул использовался потенциал Леннард-Джонса.

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

плотность и0, масса т, диаметр молекул <у.

Для моделирования уравнения Больцмана используется консервативный проекционный метод, разработанный профессором Ф.Г. Черемисиным. Этот метод является ключевым для корректного описания газокинетических процессов в переходной, но очень важной для широкого класса прикладных задач, области по числу Кнудсена 0.01<Кп<10.

Кроме того, в первой главе дано описание основ физики ударных волн, приведены основные положения и термины.

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

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

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

на ранее разработанном проекционном методе вычисления интеграла столкновений.

Высокая степень модульности проблемно-

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

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

Эффективное распараллеливание позволило, с одной стороны, значительно сократить время расчётов (с недель до часов), с другой — обеспечить саму возможность прецизионного расчета на мелких сетках, которые требуют больших объёмов памяти. На сегодняшний день солверы используют две технологии распараллеливания - MPI (Message Passing Interface) применяется для вычислений на

кластерах из обычных процессоров. Nvidia CUD А — для расчетов на видеокартах.

Входные данные

Параметры задачи

■ ......: : ............ ■ ........................

s I Конфигурационный XML файл

Интерактивная оболочка

Расчетная сетка

Генератор сеток прямоугольных

Рис. 1. Схема проблемно-моделирующей среды. Для анализа результатов компьютерного моделирования применяются различные программные продукты. Для работы с прямоугольными сетками применяется система Bkviewer, ориентированная на анализ выходных данных в случае моделирования систем с прямоугольной геометрией.

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

Солверы Визуализация

/1 ••• 1\

ггт СШ ПГШ

Рис. 2. Деление сетки между процессами.

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

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

Для реализации алгоритмы использовалась технология NVIDIA CUDA, которая предоставляет разработчику набор расширения языка С++ для запуска вычислительного кода на

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

Подготовка интегрирующей сетки проводилась на CPU. Для этого использовлся код из солвера RectSolver.

Так как все данные находятся на GPU и скорость обмена данными между CPU и GPU может быть узким местом, то расчет макропараметров происходит так же на GPU. В итоге, копирование производится для нескольких макропараметров, что уменьшает объем данных, копируемых из памяти GPU на несколько порядков.

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

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

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

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

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

тепломассопереноса разреженного газа в сложных системах — ударных волнах в разреженном газе и микронасосе Кнудсена. Все исследования проведены на основе разработанного и изложенного в главе 2 модуля проблемно-моделирующей среды.

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

В процессе работы был проделан ряд расчетов. Варьировались объемы сосудов, размеры трубки и числа Кнудсена. Так же рассматривались различные геометрии нагревания стенок. При большом объеме сосудов процесс развивается слишком медленно, ибо газ в горячем сосуде должен нагреться и излишки должны перетечь в левый сосуд. В проведенных расчетах объем сосудов принимался равным 10 объемам щели. Это позволяло понять качественную картину и получить адекватное время развития процесса. Отношение длины щели к ширине так же имеет большое влияние на скорость развития процесса. Было выбрано отношение длины к ширине 5:2. Были проведены расчеты для широкого диапазона значений чисел Кнудсена — от 0.03 до 5.

На рис. 3 показаны линии тока для установившегося режима для Кп=0.05.

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

что в разреженном газе многие явления протекают по-другому и, кроме того, сильно зависят от степени разреженности газа (которая описывается числом Кнудсена).

0.00« 0.002 0.003 0.004 0.005 0.000 0.007 0008 0.009 0.01 0.011 0.018 0.013

Рис. 3. Линии тока в установившемся режиме.

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

Моделируемая ударная труба показана на рис. 4 и состоит из двух секций: более широкой толкающей секции (секция А) и узкой измерительной секции(секция В). В начальный момент времени давление в толкающей секции в 10 раз больше, чем в измерительной, а температуры в обоих секциях равны.

Распад начального разрыва давления создает ударную волну (УВ), которая движется внутри измерительной секции. Теоретическое значение числа Маха УВ, вычисленное по соотношениям Рэнкино-Гюгонио на основе одномерной теории распада произвольного разрыва, дает М = 1.55. В работе рассматриваются две геометрические размерности задачи: двумерное (21)) плоское течение, что соответствует течению в щели, и трехмерное (ЗО) течение в канале квадратного сечения. Основные результаты получены для случая 30, который требует значительно больших вычислительных ресурсов, но является более реалистичным.

А

Рі

\л/

2Н

Ро

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

Распространение УВ в канале квадратного сечения было детально изучено для Кп = 0.05. Длина канала в расчетах

17

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

Temperature Velocity -Local Mach number

Density Temperature Velocity -Local Mach number

4 \

/7 .4.

20 40 60 80 100 -40 -20 0 20 40 60 00 1

Рис. 5. Параметры газа на центральной линии для Кп=0.05: а) 1=30; б.) 1=60

2D. Кп=0.05 "-

ч \ 3D, Кп=0.05 .........

■ Ч

i V

а.) : 1

20 30 40 50 60 70

90 100 20 40

Рис. 6. Число Маха вдоль ударной трубы: а) 2Б и 30 геометрии, Кп=0.05; б) ЗБ геометрия, Кп=0.025

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

18

максимального абсолютного значения. Чтобы найти эту точку, мы вычисляем вторую производную графика плотности и находим ее нулевое значение. Вычисления производных производятся численно с 3 порядком точности. Определение нулевого значения второй производной происходило в два этапа - сначала определялись узлы, между которыми находится нулевое значение, далее вторая производная(заданная в виде сеточной функции) вокруг этих узлов аппроксимировалась линейной функцией и определялось точное значение. На рис. 6 замедление УВ показано в зависимости от координаты вдоль канала.

В заключении приведены основные результаты диссертации.

Основные результаты диссертации

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

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

3. Проведено комплексное исследование на основе математического и имитационного моделирования стационарного течения разреженного газа в каверне и двумерного плоского течения разреженного газа в насосе Кнудсена.

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

Список публикаций соискателя по теме диссертации

1. Клосс Ю.Ю., Черемисин Ф.Г., Шувалов П.В. Решение уравнения Больцмана для нестационарных течений с ударными волнами в узких каналах // Журнал вычислительной математики и математической физики - М., 2010. - Т.50, №6 - С. 1148-1158.

2. Клосс Ю.Ю., Черемисин Ф.Г., Шувалов П.В. Решение уравнения Больцмана на графических процессорах // Вычислительные методы и программирование. Раздел 1 - М., 2010. - Т.11 - С. 144-152.

3. Додулад О.И., Клосс Ю.Ю., Черемисин Ф.Г., Шувалов П.В. Моделирование распространения ударной волны в микроканале на основе решения уравнения Больцмана // Математическое моделирование - М., 2010. - Т.22, №6 - С. 99-110.

4. Додулад О.К, Клосс Ю.Ю., Мартынов Д.В., Рогозин O.A., Рябченков В.В., Шувалов П.В., Черемисин Ф.Г. Проблемно-моделирующая среда для расчетов и

анализа газокинетических процессов // Нано- и микросистемная техника - М., 2011. - № 2 - С. 1217.

5. Клосс Ю.Ю., Черемисин Ф.Г., Шувалов П.В. Разработка программного солвера для решения задач динамики разреженного газа в кластерной архитектуре // Вестник компьютерных и информационных технологий - М., 2011 - №5 - С. 21-26.

6. Клосс Ю.Ю., Черемисин Ф.Г., Шувалов П.В. Взаимодействие ударной волны с пограничным слоем в узком канале // Математическое моделирование - М., 2011 - Т. 23, № 4 - С. 131-140.

7. К loss Yu.Yu., Shuvalov P.V., Tcheremissine F.G. Solving Boltzmann equation on GPU // Procedía Computer Science - Amsterdam 2010 - V.l, 1.1 - P. 1077-1085.

8. Князев A.H., Сакмаров A.B., Цуриков Д.Ф., Шувалов П.В. Решение уравнения Больцмана на GPU // Материалы XVII международной конференции по вычислительной механике и современным прикладным программным системам (ВМСППС'2011) - Алушта, 2011 - С. 246-248.

9. Свидетельство о государственной регистрации программы для ЭВМ № 2010613639. Программный комплекс для моделирования газокинетических процессов на основе численного решения уравнения Больцмана Rogsolv / Клосс Ю.Ю., Черемисин Ф.Г., Рогозин O.A., Дербакова Е.П., Шувалов П.В. 2010г.

10. Свидетельство о государственной регистрации программы для ЭВМ № 2010613640. Программный комплекс для визуализации результатов моделирования явлений в разреженном газе BKViewer / Клосс Ю.Ю., Черемисин Ф.Г., Хохлов Н.И., Дербакова Е.П., Аникин Ю.А., Рогозин O.A., Шувалов П.В., 2010г.

Шувалов Павел Вадимович

МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ РАЗРЕЖЕННОГО

ГАЗА НА ОСНОВЕ КИНЕТИЧЕСКОГО УРАВНЕНИЯ БОЛЬЦМАНА НА КЛАСТЕРНЫХ И ГРАФИЧЕСКИХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ

Автореферат

Подписано в печать 16.05.2013. Формат 60 х 84 1/16. Усл. печ. л. 1,25. Тираж 120 экз. Заказ номер 895.

Федеральное государственное автономное образовательное учреждение высшего профессионального образования "Московский физико-технический институт (государственный университет)"

Отдел оперативной полиграфии "Физтех-полиграф" 141700, Московская обл., г. Долгопрудный, Институтский

пер., 9

Московский физико-технический институт (государственный университет) Кафедра моделирования ядерных процессов и технлогий

На правах рукописи 04201357930 УДК 001.891.57:53

ШУВАЛОВ Павел Вадимович

Моделирование течений разреженного газа на основе кинетического уравнения Больцмана на кластерных и графических вычислительных системах

Специальность 05.13.18 «Математическое моделирование, численные методы и комплексы программ»

Диссертация на соискание ученой степени кандидата физико-математических наук

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

кандидат физико-математических наук доцент Ю.Ю. Клосс

МОСКВА-2013

СОДЕРЖАНИЕ

ВВЕДЕНИЕ................................................................................4

ГЛАВА 1. МАТЕМАТИЧЕСКИЕ ОСНОВЫ КИНЕТИЧЕСКОЙ ТЕОРИИ ГАЗОВ.........................................................................6

1.1. Основы описания газокинетических процессов..................................................6

1.2.Граничные условия и безразмерные переменные.................................................8

1.3. Основные принципы консервативного проекционного метода вычисления интеграла столкновений.................................................................................................11

1.4. Введение в разностные схемы первого и второго порядков..................................17

1.5. Распространение ударных волн в микроканалах (основные положения)..................23

ГЛАВА 2. АЛГОРИТМЫ, ПРОГРАММНЫЕ СИСТЕМЫ И ИХ

РЕАЛИЗАЦИИ...................................................................36

2.1. Проблемно-моделирующая среда на прямоугольных сетках.................................36

2.2. Солвер для численного решения кинетического уравнения.................................49

2.3. Алгоритмы и програмная реализация вычислений на кластерных архитектурах........64

2.4. Реализация программного солвера на графических процессорах...........................69

2.5. ОсновныИСЬ..........................................................................................87

ГЛАВА 3. МОДЕЛИРОВАНИЕ И АНАЛИЗ ПРИКЛАДНЫХ ФИЗИЧЕСКИХ СИСТЕМ И ЯВЛЕНИЙ.....................................91

3.1.Анализ и компьютерные модели двумерного насоса Кнудсена..............................91

3.2.Изучение и анализ ударных волн в микроканалах на основе программно-моделирующей среды...............................................................................101

3.3. Моделирование индуцированного течения в каверне на графических процессорах..! 17

ЗАКЛЮЧЕНИЕ........................................................................................119

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ.....................................120

ВВЕДЕНИЕ.

Для численного моделирования течений разреженного газа во второй половине XX века были разработаны два основных подхода: метод прямого статистического моделирования и конечно-разностное решение кинетического уравнения Больцмана. В первом подходе моделируется процесс случайных столкновений и перемещения большого числа шестимерных векторов в фазовом пространстве, обозначающих молекулы газа [1]. На их основе вычисляются среднестатистические значения физических величин, отождествляемые с макроскопическими параметрами газа. Данный метод успешно применялся при расчете сверхзвуковых течений разреженного газа, но для медленных течений этот подход может давать недостоверные результаты из-за присущего методам Монте-Карло статистического шума. В настоящих исследованиях развивается второй из названных подходов, который не содержит статистических флуктуаций в решении и позволяет разрешать ничтожно малые изменения параметров течения газа. Последнее обстоятельство важно для разработки компьютерных моделей с высокой точностью, описывающих медленные течения разреженного газа, характерные для условий в микро- и наноустройствах. Следует отметить, что большая размерность задачи, для которой необходимы значительные объемы оперативной памяти и вычислительной мощности, необходимость варьирования геометрических параметров микроустройств, физических характеристик газа требует использования современных суперкомпьютерных систем с различной архитектурой.

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

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

Значительное улучшение метода прямого решения уравнения Больцмана было достигнуто Ф.Черемисиным в 1996г. [2], когда был разработан проекционный метод вычисления интеграла столкновений Больцмана, строго сохраняющий массу, импульс и энергию в процессе молекулярных столкновений. Тем самым был устранен основной недостаток метода прямого решения уравнения Больцмана по сравнению с методом Берда. В дальнейшем консервативный проекционный метод был обобщен на смеси газов и газы с внутренними степенями свободы молекул. Этот метод позволяет экономично рассчитывать сверх- и гиперзвуковые течения, но особенно эффективен при расчете медленных и слабо возмущенных течений, к которым относятся типичные течения в микро-каналах.

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

Глава 1. МАТЕМАТИЧЕСКИЕ ОСНОВЫ КИНЕТИЧЕСКОЙ ТЕОРИИ ГАЗОВ.

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

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

1.1. Основы описания газокинетических процессов.

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

Кинетическое уравнение Больцмана может представлено в следующем виде:

(' о/* 00 2я Ьт

§"+ <= Е / 1 ¡(///Г-Л/^^ейр, (1.Ы)

и1 т ил У -оо О О

где ¿=/(/,х,р() ; /. = /(¿,х,ру) ; // = /(^х,р/) ; ¿^/(^р,.1) - функции распределения газов; / — время; р — импульс молекул; т - масса частиц; g — относительная скорость при столкновении; Ъ, е - параметры столкновения; р,, ру, р;', р,-' - импульсы частиц до и после столкновения.

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

Левая часть уравнения описывает движение частиц в пространстве без столкновения между собой. Интеграл в правой части уравнения (1.1.1)называется интегралом столкновений Больцмана и обозначается следующим образом

00 1я Ът У -00 о О

С учетом введенного обозначения уравнение (1.1.1) переписывается в виде

т' дх ''

Макроскопические параметры газа вычисляются на основе значений функции распределения по формуле

|/(1 р.

Вектор средней скорости

Компоненты тензора давления газа в данной точке Температура как мера средней кинетической энергии

т=

1

¡р1М РЬ •

Ъпкт

при этом учитывается только хаотическая составляющая скоростей молекул Ра =(А.^»Рг) = Р-Рс» ГДе Рс - средний импульс молекул. Вектор потока энергии определяется по формуле

В большинстве расчетов для описания взаимодействия молекул использовался потенциал Леннард-Джонса:

Также возможно описание взаимодействия частиц как твердых сфер диаметром о, степенным потенциалом вида и(г) = / г* , где 4<у<\2 , или любым другим сферически-симметричным потенциалом взаимодействия. 1.2. Граничные условия и безразмерные переменные.

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

При ударе об стенку считается, что часть молекул 1 — а, отражается от стенки зеркально, а другая часть а, - диффузно с максвелловским распределением:

Здесь f¡ - функция распределения падающих на поверхность молекул, аг,пг,Тг -свободные параметры. Параметр а, характеризует долю тангенциального импульса, передаваемого падающими молекулами стенке. При полностью зеркальном отражении а, = 0, а при полностью диффузном отражении а, = 1. Параметр Тг определяет температуру

и (г)=4в[(ст / г)12 — (а / г)6 ],

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

□ ^(хд щк)бог = □ да/с)

Ск<о ск<о

Подставив сюда значение /г из формулы(1.1.1) получим значение пг:

\ М^Ж^Ж

п=2.

т £к< О

' 2 кТл л

В некоторых задачах рассматриваются граничные условия с движущейся стенкой. Примером такой задачи является задача об индуцированном течении в каверне [3]. В таких случаях используется следующая форма функции распределения отраженных молекул (приводится запись для полностью диффузного отражения):

2 кТж

Новый свободный параметр иг задает скорость движения стенки.

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

р / , ♦ 6 » х .. _ те3 т ,1ЛП

Р =-=4 t =—,Ь =-,х =/--,т =—. (1.2.1)

/я0у0 т0 а Л Л0У0 т0

Символом «*» обозначены безразмерные переменные. Здесь введены следующие обозначения:

кТп

у0 = ./ у - характерная скорость,

т,

о

Я = —г=—-—7--длина свободного пробега,

Ы2тт0аэфф

а.. ~ о"л/П(2'2) - эффективное сечение,

эфф

т0 = у - характерное время.

Под эффективным сечением понимается молекулярный диаметр молекул такого газа из твердых сфер, вязкость которого равна вязкости газа с данным молекулярным потенциалом. Коэффициент <Г2(2'2) задает отличие эффективного сечения от молекулярного диаметра. Очевидно, что для потенциала твердых сфер 0(2'2) = 1. Вязкость газа из твердых сфер задается формулой [4; 5]:

5 \ктТ ,„ _

1бег V л

Потенциал Леннарда-Джонса [6]рассчитывается по формуле

Где £ - глубина потенциальной ямы, Т0 - нормировочная температура. Вязкость газа с потенциалом Леннарда-Джонса задается формулой:

где Т' =Т0 / £. Значения П(2,2) для потенциала Леннарда-Джонса приведены в [4; 5]. Из (1.2.2) и (1.2.3) следует, что:

Так, например, для аргона (глубина потенциальной ямы 120К) при300Кзначение поправочного коэффициента равно 1.093.

После подстановки формул (1.2.1) в уравнение Больцмана получим

- V*« Я/с/ ■/; --ть'л\1Р;.

тотоуо ЛтЛ Щ , '"Л '"о1о

5 / V п /

Подставим тп = л/ и поделим на 0

/ /Ял

=Ш«''/; ■-//>Ул-ар1*,

01 т ОХ 1

Подставляем Я = / <— 2 в коэффициент при интеграле, получим / Ч2пщиэфф

"0Ур3у04о~2Я =по.гл= п0а2 _ п0а2 __1

V«, ° Ллп^фф 42яп0(72П{2'2) ЛлП(2'2)

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

дх л/^/и.^ J

Символ «*» в дальнейшем будем опускать. В безразмерном виде уравнение Больцмана принимает вид

д( т' дх

где 7<= №. '/, '-ДЖ,

1.3. Основные принципы консервативного проекционного метода вычисления интеграла столкновений.

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

Интеграл столкновений обладает свойством консервативности по веществу, импульсу и энергии. Проинтегрируем интеграл столкновений с какой либо весовой функцией от скоростей:

I, = ¡рЩ/К-Ж^ЬМей^

он равен интегралу (всего лишь переобозначены переменные)

К = ^Ш-т^аъЧе'м;^

Воспользовавшись соотношениями

£ = = Ъ,ё = в,

получаем Сложив имеем

К = \ ¡(<р - Ч*Ш- Ж са> ^ . Интеграл симметричен относительно % и , поэтому

И наконец, получаем

Я'Х ДЖ)8Ь¿Ъ¿г= сру^)«ц.

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

Теперь, приняв за <р(%) = , получаем искомые соотношения

Однако, при непосредственном численном интегрировании показанная консервативность нарушается. Это требовало разработки другого метода лишенного данного недостатка, что было сделано [2; 7]. Далее опишем его метод.

В скоростном пространстве выбирается область £1 объема V, в которой строится сетка из И* равноотстоящих скоростных узлов Е = . В качестве области в

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

На введенной скоростной сетке функция распределения и интеграл столкновений определяются следующим образом:

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

Интеграл (1.3.1) является 8-ми мерным, что усложняет его численное вычисление. Одним из подходов интегрирования является метод Монте-Карло. Но данный метод

обладает плохой сходимостью при увеличении числа интегрирующих узлов -О^^/^-^.В

связи с этим, применяется метод интегрирования основанный регулярных сетках Коробова [8], ошибка интегрирования при их применении в среднем падает как ^(/^дг) •

В данном случае интегрирующая сетка Коробова строится для , , Ьу, еу из Ыу узлов в области П х £~2 х [0,2лг] х [0, Ьт ] так, что % , совпадают с узлами скоростной сетки Е.

В системе центра инерции относительная скорость сталкивающихся молекул поворачивается в плоскости столкновения на угол рассеяния в. Таким образом, зная угол

рассеяния, нетрудно найти скорости молекул после столкновения. Угол рассеяния 9, в свою очередь, зависит от потенциала парного взаимодействия молекул 17у(гц).

Скорости после столкновения , не лежат в узлах скоростной сетки, поэтому

две последние £-функции в выражении (1.3.1) требуют преобразования:

где узлы , , , _5 и коэффициент ^ выбираются так, чтобы выполнялись закон сохранения импульса и энергии (закон сохранения вещества при данной замене уже выполняется). Коэффициент гу находится из сохранения энергии - £„ = (1 - гу)Е] + гуЕ2. Откуда

где Е0=$1 + Е1 = £ и Е2 = . Выбор

узлов Лу, Лу+ должен обеспечивать выполнение О < гу < 1. Сохранение импульса обеспечивается симметричным выбором аппроксимационных узлов.

Для вычисления f а /д предпочтительнее

использовать следующую интерполяцию

favfpv = ^К-^РК ) " ^У/^+зУц,-1) "' которая обеспечивает равенство нулю интеграла столкновений для максвелловского распределения. Максвелловская функция распределения является равновесной, что означает отсутсвие эволюции функции распределения за счет релаксации, то есть отличие от нуля интеграла столкновений противоречило бы физическому смыслу.

Окончательно оператор интеграла столкновений приобретает вид

1Г+ ¿г'* )-(£;* (1-3.2)

^=1

аппроксимации скоростей на сетке.

2 А 7-2

где В =

УлЪ1ти 4М

А, = <Д/А "(Д

Далее применяется схема «непрерывного счета». Сведем данное уравнение к интегральному уравнению

< -л

Подставляя его в интеграл (1.3.2), получаем

где Ду —V -й член суммы.

Способ вычисления оптимальных коэффициентов сеток Коробова основан на следующей теореме [8].

Теорема 1 Пусть для целых г функция Н(г) определена равенством

<5» р

Рм

1-2 *

\2

1-2<

ч I Р

(1.3.4)

где р - простое число, большее ^ .

Если при г = а достигается минимум функции Н(г) на интервале \<г<р-\, то целые ах =1 ,а2 =а,...,а3 = а"'1 будут оптимальными коэффициентами по модулю р [8]. Минимум функции Н(г) можно установить путем сравнения ее значений �