CC BY
53
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 - 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Программный комплекс <<Несветай-3Д>>
моделирования пространственных течений
одноатомного разреженного газа
# 06, июнь 2014
DOI: 10.7463/0614.0712314
Титарев В. А.
УДК 519.677
Россия, Вычислительный центр им. A.A. Дородницына РАН
titarev@ccas.ru
Введение
Моделирование процессов в микроустройствах (микроканалах, микронасосах и т.д.) и анализ обтекания спускаемых аппаратов произвольной (реалистичной) формы требуют создания методов компьютерного моделирования пространственных течений разреженного газа. Одним из таких методов является численное решение кинетического уравнения Больц-мана для функции распределения молекул по скоростям. В виду высокой размерности уравнения и сложной структуре интеграла столкновений аналитическое решение кинетического уравнения возможно лишь в специальных случаях, в то время как численное решение требует значительных затрат машинного времени. Заметим, что в приложениях основной интерес представляет не функция распределения, а некоторое число ее первых моментов — макроскопических параметров газа, таких как плотность, скорость, давление. Поэтому для практических расчетов можно заменить точное уравнение Больцмана некоторым приближенным уравнением, таким чтобы макропараметры газа, вычисленные с помощью этого приближенного уравнения, совпадали или были бы близки с макропараметрами, вычисленными с использованием точного уравнения.
В настоящее время для моделирования течений одноатомного разреженного газа широкое распространение получило модельное уравнение неполного третьего приближения, получившие название модели Шахова, или S-модели [1,2]. Сравнения с расчетами по точному уравнению Больцмана, методу Монте-Карло и с экспериментальными данными показали хорошую точность S-модели [3, 4, 5, 6, 7, 8, 9]. Успешное применение S-модели стимулировало ее обобщение на течения двухатомного газа с учетом вращательных степеней свободы, и развитие соответствующего численного метода [10, 11, 12]. Тем не менее, несмотря на
относительную простоту, S-модель является сложным интегро-дифференциальным уравнением высокой размерности, численное решение которого требует развития высокоточных параллельных численных методов решения.
В большинстве имеющихся работ для численного интегрирования кинетического уравнения с точным либо модельным интегралами столкновений используются явные методы дискретизации по времени и расщепление уравнения на процессы переноса и релаксации. Используемые разностные схемы реализованы на многоблочных структурированных [13], неструктурированных декартовых [14] и тетраэдральных [15] сетках. Созданные на основе данных методов комплексы программ позволяют проводить численное моделирование течений разреженного газа для объектов достаточно сложной пространственной формы. Однако использование явных методов дискретизации по времени и пространственных сеток без адаптации к поверхности тела ведет к большим затратам машинного времени. Проведение вычислений с гарантированной точностью становится затруднительным даже для относительно простых задач, таких как истечение газа через круглый канал. Таким образом, актуальной является задача разработки более эффективных и универсальных методов решения пространственных задач теории разреженных газов на основе численного решения кинетических уравнений.
В настоящей работе приводится обзор последних результатов по разработке [16, 17, 18] и применению нового численного метода решения кинетического уравнения Больцмана с Б-модельным интегралом столкновений [1, 2] и реализующего его программного комплекса «Несветай-3Д>>. Основными блоками метода являются квазимонотонная схема второго порядка точности на гибридных неструктурированных сетках, консервативная процедура расчета макроскопических параметров газа и полностью неявный метод дискретизации по времени без итераций на верхнем слое. В качестве примеров применения метода рассматриваются течения разреженного газа в каналах и обтекание модели спускаемого аппарата. Приводится сравнение с известными в литературе результатами. Исследуется влияние различных параметров счета (разрешения пространственной сетки, значения числа Куранта, числа процессоров и т.д.) на точность счета и скорость сходимости к стационарному решению.
1. Расчетные уравнения
Размерный вид уравнения. Состояние разреженного газа в точке х = (х1,х2,х3) в момент времени £ определяется функцией распределения молекул по скоростям / = = /(£, х, £), где £ = (£ь£2,£3) — компоненты вектора молекулярной скорости по направлениям (х1,х2,х3) соответственно. Макроскопические переменные, такие как числовая плотность п, температура Т, давление р, средняя скорость газа и = (и1,и2, и3) и поток тепла Я = (я1 ,Я2, q3), выражаются через функцию распределения в виде интегралов по простран-
ству молекулярных скоростей:
п = J /С£, пи = J £/<!>£, зтпЯдТ +1 т'пи2 = ^т^ £2/С£,
Я = 1 т ^ юь2/ ¿£, V = £ — и, р = тп, р = рЕдТ, и2 = иаиа, V2 = ьаьа, £2 = £а£а, С£ = .
Здесь т — масса молекулы; Яд — газовая постоянная; предполагается суммирование по повторяющимся греческим индексам в пределах от 1 до 3. Кинетическое уравнение Больцмана для функции распределения / молекул по скоростям с Б-модельным интегралом столкновений [1, 2] имеет вид
%+ьд/=*]=-(/+—/), /+=/м
т дха —
1 + 4(1 — Рг)Б,Са( С2 — 2)
п 1С V
/м = (2пП9Т)з/2 ехр(—с2), Б = п]Сгс2/С£, С = -ЩТ, С2 = СвСв,
где — = — (Т) — вязкость; число Прандтля для одноатомного газа Рг = 2/3; Рг = 1 соответствует модели БГК [19].
В работе рассматриваются несколько вариантов задания граничных условий на поверхности тела. Первое условие соответствует диффузному отражению с полной тепловой аккомодацией к температуре поверхности Т8иг. При этом функция распределения отраженных молекул записывается в виде:
/- = (ВД/2 ехр(—2дУ ' п- Ч* = —/с£■ (1)
' ?П<0
где £п — проекция вектора скорости молекулы £ на внешнюю нормаль к поверхности тела го8иг. Плотность отраженных молекул находится из условия непротекания.
Второе условие соответствует испарению или конденсации с заданными плотностью п8иг и температурой Т8иг:
/- = ОЛК-Т^ехр (—ЩТъ) ' £п > 0 (2)
Безразмерный вид уравнения. Для проведения расчетов удобнее использовать безразмерные переменные. Пусть п*, Т*, I* — некоторые характерные значения числовой плотности, температуры и пространственной координаты. Введем характерные значения давления и скорости р* = тп*Я9Т*, V* = \]2КдТ*. Перейдем к безразмерным переменным по формулам
/ х ./ / п / Р гр/ Т
х = -¡-, Ь = у t, п = —, Р = —, Т =
1 * 1 * п* р* Т *
и = £/ = 1 Я/ = // = /
V* V* тп*ь3 п*ь3
Степень разреженности газа характеризуется параметром разреженности 8, который обратно пропорционален числу Кнудсена Кп, определенному по длине свободного пробега Л*:
„ 1*р* 8 1 Л*
8 =-= —^тт", Кп = —.
Ъу/п Кп 1*
Далее будем обозначать безразмерные величины теми же буквами, что и размерные, поскольку это не приводит к неоднозначности. Кинетическое уравнение принимает вид:
д/ + = 3 (3)
81 + ^ дхг = 3 (3)
-а
где
3 = V(/+ - /), V = 8р, /+ = /м (1 + |(1 - Рг^Цс2 - |)) ,
5 = ПТТзз/2, /м = (ЛТрехр(-с2).
Для макроскопических величин получаем:
п,пи, 3пТ + пи2, ^ = ^ ^1, Щ. (4)
В безразмерных переменных граничное условие диффузного отражения (1) переписывается в виде
/- =(пТ^(-£)■ п- = « = - ^ (5)
?П<0
Для условия испарения/конденсации (2) имеем
/" = (пТ^- £) • & >0-
В дальнейшем принималось ц ~ л/Т, что соответствует межмолекулярному взаимодействию по закону жестких сфер.
2. Численный метод
Общий вид метода. Стационарное решение кинетического уравнения (3) находится методом установления по времени. В пространстве молекулярных скоростей несобственные интегралы от функции распределения заменим собственными интегралами по некоторой ограниченной области. Вид области и ее размеры определяются решаемой задачей. Введем в расчетной области расчетную (как правило конечно-объемную) сетку с узлами (центрами ячеек) ^ = (^ ). Общее количество узлов сетки обозначим через N. Пусть 3к — вектор, компонентами которого являются к-компонента молекулярной скорости во всех узлах сетки:
3к = (£кЪ ^k2, £>k3, . . . )Т.
Функции f, f будем задавать в центрах ячеек скоростной сетки, интерпретируя их как векторы длины с компонентами
/ = / (*, X, Ь), /^ = /<*>(*, X, Ь).
Кинетическое уравнение (3) переписывается в виде системы из N5 уравнений, записанной в векторной форме
д д
аf + ах" <Е°°f ' = 7 • 7 ="<f- f(6)
Здесь через ° обозначена операция покомпонентного умножения двух векторов: с = а ° Ь — вектор с компонентами сг = агЬг. Отметим, что в существующей литературе как правило используется скалярная запись кинетического уравнения для заданного значения ^. Однако для реализации на современных многоядерных процессорах с векторными регистрами удобнее использовать запись кинетического уравнения как векторного закона сохранения (6).
Помимо удобства программной реализации, векторная запись (6) дает возможность новой по сравнению ранними работами [16] интерпретации неявной разностной схемы для кинетического уравнения.
Расчетная область в физических переменных разбивается на контрольные объемы (ячейки) V тетраэдальной, пирамидальной, гексаэдральной или призматической формы. Каждая из ячеек образована несколькими треугольными или четырехугольными гранями Лц, где I — номер грани. Общее количество ячеек равняется ^зрасе. Для получения разностной схемы уравнение (6) интегрируется по четырехмерному контрольному объему области V х [¿га, ¿""+1], объемный интеграл от конвективных производных преобразуется в сумму поверхностных интегралов по границам граням расчетной ячейки от потоков ^к, спроецированных на нормали иц к граням. Стандартная аппроксимация интегралов от потоков и правой части (модельного интеграла столкновений) на верхнем слое приводит к неявной одношаговой схеме вида
Д/г Рга+1 рга+1 _ 1 \ Л Л>п+1 I тга+1 (п\
— = Яг , Яг Ф + ^ , (7)
где f П = f (¿п, — среднее значение функции распределения по ячейке V в момент времени ¿га; Д/г = f"+1 - /П — приращение функции распределения за шаг по времени. Нормальный численный поток через грань Лц определяются по формуле
Ф^1 = J ° f ^ йиИ = пИ31 + п2132 + п313з. (8)
Аа
Вектор состоит из проекций узлов скоростной сетки на внешнюю единичную нормаль иц грани I ячейки V*.
Проведем приближенную линеаризацию (7) по времени. Для модельного интеграла столкновений будем опускать достаточно сложную зависимость макроскопических величин от Г I так что
« ^ - иПДГг.
Формула для численного потока (8) содержит интеграл от функции распределения по грани ячейки. Как обычно в противопоточных схемах типа Годунова [20, 21, 22], на каждой грани есть два значения решения. Первой значение f- получается экстраполяцией на грань в рамках ячейки V. Второе значение f + соответствует экстраполяции на грани Лц снаружи ячейки. Численный поток через грань I ячейки V вычисляется по формуле
ФП = , Гы )1М (9)
где ц — номер ячейки, соседней к грани I (если такая существует); ^ — номер соответствующей грани ячейки ц. Функция двух векторных аргументов С(ГГ+) по аналогии с вычислительной газодинамикой является решателем задачи о распаде разрыва с начальными данными Г-, Г+. В рамках настоящей работы будут рассматриваться два вида функции С: точный решатель
1 21
G(f-, f+)™1 = 1 ΄„ ◦ [(f- + f +) - sign(î„„) ◦ (f + - f-)] (10)
и решатель типа Русанова
G(ff+)Rus = 2 [înu ◦ (f- + f+) - £Г(/ + - f-)], = max ^ |. (11)
2 3
В общем случае схемы высокого порядка значения на гранях ячейки находятся с помощью процедуры реконструкции решения и определяются по значениям f в нескольких соседних ячейках, количество которых определяется типом и порядком пространственной аппроксимации разностной схемы. Проводя линеаризацию численного потока, при вычислении производных будем учитывать зависимость потока только от значений f в ячееке i и ее непосредственном соседе i по грани l, что соответствует противопоточной схеме первого порядка аппроксимации по пространству. Так же учтем, что j -я компонента ФП+1 зависит только от j -х компонент f ь f, . Получаем
Здесь учтено, что производные от Ф^ по Гь Г^ являются диагональными матрицами и могут быть записаны в виде векторов. Численный поток на нижнем слое ФП вычисляется с полным порядком аппроксимации пространству, используя необходимое количество соседних ячеек. Группируя члены, получаем линеаризованный вид разностной схемы (7):
(1++ Д £ Г) ◦ дг ,+Д £ (Г) ◦ ДГ „ = дшг. (12)
Здесь /5 = diag(1,..., 1)т — единичная матрица размерности N5. Соотношения (12) связывают приращения решения в ячейке V и ее соседях для всех ячеек пространственной сетки г = 1,..., Д,расе. Описание численного метода будет законченным, если будут указаны следующие части численного алгоритма: способ вычисления численных потоков в правой части схемы (12) c высоким порядком аппроксимации, метод вычисления макроскопических величин, входящих в модельный интеграл столкновений, и способ решения линейно системы уравнений для приращений решения Д/¿. Подробное описание этих элементов численного алгоритма дается ниже.
Аппроксимация оператора переноса на произвольной сетке. Численный поток через грань I ячейки V вычисляется по формуле (9) c использованием точного решения задачи о распаде разрыва — функции ^ехас\ Значение f- = f " получается экстраполяцией на грань в рамках ячейки Второй значение f + = fсоответствует экстраполяции на грани Лц снаружи ячейки. Для схемы первого порядка аппроксимации достаточно положить = f Однако для расчета течений газа при умеренных и малых числах Кнудсена необходимо использовать схемы как минимум второго порядка аппроксимации. В настоящем методе используется кусочно-линейное представление функции распределения в пространственной ячейке, обеспечивающее второй порядок аппроксимации по пространству
Здесь егк — базисные функции для ячейки Коэффициенты еЩ, с"2, с" данного представления вычисляются методом наименьших квадратов с использованием значений / в ячейках шаблона реконструкции, который формируется путем рекурсивного добавления к ячейке V необходимого числа ее непосредственных соседей. Также следует отметить, что метод наименьших квадратов применяется в локальной системе координат, связанной с ячейкой что позволяет улучшить число обусловленности соответствующей матрицы. Переход к локальной системе координат показан на рис. 1, на котором представлен типичный шаблон схемы второго порядка на гексаэдральной сетке.
Р" = + е"1ег1(х,У,^) + сП2ег2(х,У,^) + ^¿з^у,^.
(13)
I
Рис. 1. Шаблон метода наименьших квадратов в физической (слева) и локальной (справа) системах координат для гексаэдральной сетки
Подробное описание процедуры построения кусочно-полиномиального представления решения описано в [23] для тетраэдральных ячеек и в [16, 24] для трехмерных ячеек произвольной формы. Для линейной схемы можно положить средние значения функции распределения Г™ по каждой грани I ячейки равными средним значениям рП полученного реконструкционного многочлена р™(х). При этом рП представляются в виде линейной комбинации значений функции распределения в ячейках шаблона, коэффициенты которой определяются только геометрией пространственной сетки. Как известно, линейные схемы высокого порядка точности немонотонны в областях быстрого изменения решения [20]. Для подавления паразитных осцилляций в блок реконструкции вводится так называемый ограничитель наклонов [21]. Для расчета стационарных течений удобен ограничитель наклонов, предложенный в [25].
Известно, что точность аппроксимации методов наименьших квадратов падает для сильно вытянутых ячеек [26]. Также, для некоторых классов задач, например, при расчете течений в длинных каналах, течение может иметь выделенное направление, градиенты вдоль которого значительно меньше поперечных производных, что не используется в полностью трехмерном методе. Поэтому для гексаэдральных ячеек предусмотрена возможность переключения на локально-одномерный (структурированный) способ вычисления значений на гранях ячеек [27]:
где Аг — расстояние от центра ячейки i до центра грани l; SL и SR — левая и правая оценки производной (наклона) решения; ^ld(a, b) — одномерный ограничитель наклонов, применяемый покомпонентно. В настоящем методе могут использовать разные ограничители, включая классический ограничитель Колгана [21] (типа minmod), monotonized central (MC) ограничитель [28], а также ограничитель ван Албада [22].
Далее для всех граней ячеек сетки, соответствующих граничному условию, формула для численного потока (9) видоизменяется так, чтобы аппроксимировать граничное условие. При этом макроскопические параметры, входящие в условие диффузного отражения (5), вычисляются таким образом, что условие непротекания выполняется тождественно.
Нахождение макроскопических величин. Вычисление модельного интеграла столкновений в правой части схемы требует знания макроскопических величин: плотности, температуры, компонент скорости и вектора теплового потока. Данные величины определяются как интегралы от функции распределения по молекулярной скорости (4). Прямая аппроксимация выражений (4) посредством квадратурной формулы с весами Uj (индексы i и n опущены для простоты)
f П = f ? + ^id(SnL, SR)Ai,
(14)
\ / 1 \
(15)
приводит к нарушению в численной схеме законов сохранения массы, импульса и энергии [29]. Появляющийся при численном счете паразитный источник прямо пропорционален параметру разреженности 8. В результате моделирование переходных и сплошносред-ных течений газа может требовать неприемлемо большого числа узлов N5 в скоростной сетке.
Короткий обзор существующих способов построения консервативных схем для модельного кинетического уравнения приведен в [17]. В настоящей работе используется относительно новый метод построения консервативной аппроксимации интеграла столкновений [29, 30, 31]. Макропараметры газа находятся как решение системы уравнений, используемой для построения Б-модельного уравнения:
N
Н (Ш) = £
¿=1
1
£2
2
V2/
(
(/Г - /? V? +
0 0
0
\
0, Ш
у 2 Рг д у
п и
Т д
(16)
Для каждой пространственной ячейки г система из восьми уравнений (16) решается методом Ньютона
М(Ш*-1 )(Ш - Ш5-1) = -Н(Ш*-1), 5 = 1, 2, .
М =
а н
В качестве начального приближения используются значения, полученные из приближения (15).
В частном случае Рг = 1 (модель БГК [19]) функция /^ не содержит значений вектора поток тепла, так что последние три уравнения системы (16) становятся избыточными. В данном случае система (16) совпадает с методом, предложенным в [32, 33] для БГК-модели.
Медленной частью метода Ньютона является вычисление матрицы Якоби М системы уравнений (16), представляющей из себя дискретные суммы от производных /^ по макроскопическим переменным. Точное вычисление М требует существенных затрат машинного времени, особенно для более сложных кинетических моделей, таких как модель Рыкова двухатомного газа [10]. Следуя [31], перейдем в выражении для матрицы М от численного интегрирования к точному (аналитическому). В результате М выражается в явном виде через макропараметры газа. Как правило, использование приближенной формулы для М существенно уменьшает время счета для 8 < 1000, несмотря на потерю квадратичной сходимости итераций.
Решение уравнений на верхнем слое. При решении системы (12) будет полагать, что в левой части численный поток аппроксимируется противопоточной схемой первого порядка. Возможны два варианта метода. В первом варианте метода для получения явных выражений
£
3
для производных от численного потока используется функция Русанова (11) с аргументами f ¿, f ¿г. Неявный конечно-объемный метод (12) переписывается в виде
((1 + Л^Г)/« + щ + )|Л«|) ° Afг +
+ Щ\ - Я?"1«)1А*1 ° ЛД = Л^П (17)
Более аккуратным способом дискретизации левой части является использование точной противопоточной формулы, выражаемой функцией (10). Неявная схема (12) переписывается в виде
((1 + Л^Г)/? + щ Е& ° (I« + ^¿г )\4*|) ° ЛЛ +
+ Щ\ Е ° (I« - 81§п^) IАгг | ° Лfгl = Л1Щ. (18)
Схема (17) так как не требует дорогостоящего многократного вычисления величины шгг для каждой ячейки и поэтому является менее трудозатратной с вычислительной точки зрения. С другой стороны, форма записи метода (18) является более надежной при расчете сильно неравновесных течений газа, например истечения газа в вакуум. Вариант (17) является новым и предлагается в настоящей работе, в то время как (18) использовался во всех предыдущих работах автора.
Очевидно, что обе системы (17), (18) обладают свойством диагонального преобладания. Группируя члены и выделяя в явном виде коэффициент перед Лf ¿, получаем
* ° Лf г + Л1 Е Сгг ° Лf г1 = ЛШ™ (19)
г
Прямое решение системы уравнений (19) является довольно затратной операцией. Поэтому производится приближенная факторизация разреженной матрицы системы на основе подхода, предложенного в [34]. Система решается в два этапа. Сначала производится обратный ход метода, и находятся промежуточные значения Лf *:
^ ° Лf * = - Е Сгг ° Лf * + ЛШ™, г = Ж^се, ..., 1. (20)
I: ¿1<г
Прямой ход метода выдает окончательные значения для приращения решения:
^ ° Лf г = Лf * - Е Сгг ° Лf ¿1, г =1,..., Ж^. (21)
I: ¿1 >г
Использование формул (20), (21) не требует хранения в памяти матрицы системы уравнений; выполняемое число операций линейно пропорционально числу пространственных ячеек Ж^. В результате основные во время расчета одного шага по времени основные затраты машинного времени приходятся на процедуру реконструкции и вычисления макроскопических величин. Тестовые расчеты показывают, что один шаг неявного метода требует лишь на 30% больше машинного времени, чем соответствующего явного.
3. Программный комплекс <<Несветай-3Д>>
Структура пакета. Комплекс <<Несветай-3Д>> состоит из вычислительного ядра (базовой библиотеки), непосредственно кинетического решателя и препроцессора. Исходный код пакета написан на языке Fortran 90. При разработке кода принимались во внимания требования переносимости и совместимости с различными операционными системами (Windows 7, Linux) и возможности работы на вычислительных системах различных типов, от рабочих станций до суперкомпьютеров. Переносимость и совместимость программной реализации означают, что описываемый пакет программ не использует никаких программных решений, зависящих от типа ОС и конкретной архитектуры вычислительных узлов. Параллельная конфигурация использует MPI (версия стандарта 2) для модели с распределенной па- мятью и OpenMP (стандарт 3.0) для модели с общей памятью. Работа с кодом ведется в интегрированной среде разработки (IDE) Microsoft Visual Studio 2010 с использованием компилятора Intel Fortran 14. Для управление версиями используется пакет Git.
Вычислительное ядро. Ядро пакета представляет собой набор модулей, реализующих базовые операции, необходимые для проведения расчетов. Ядро реализует процедуры чтения пространственных сеток в форматах Neutral и StarCD и построение информации о связности сетки (соседи граней, ячейке и вершин). Для аппроксимации оператора переноса используются алгоритмы реконструкции скалярных функций методом наименьших квадратов на произвольной сетке (до 5-го порядка аппроксимации включительно) и алгоритмы нелинейного ограничения решения типа Total Variation dimishing (TVD). Для проведения параллельных вычислений с помощью MPI реализовано создание дополнительных структур данных и алгоритмы параллельного ввода-вывода решения в файл перезапуска на основе процедур MPI I/O. В ядре также входят процедуры вывода пространственных и поверхностных данных в формате Tecplot.
Остановимся на особенностях представления данных в пакете. Расчетная область в физическом пространстве на неструктурированной гибридной сетке определяется списком координат узлов сетки, списком элементов и списком внешних граней сетки. Общее количество узлов, элементов и граней равняется Nvert, Nspace и Nface, соответственно. В пакете допускаются элементы четырех типов: гексаэдры, тетраэдры, призмы и пирамиды. Внешние грани этих элементов являются треугольниками или четырехугольниками. Отметим, что в общем случае четырехугольные грани могут быть неплоскими. Каждая грань задается списком вершин. Для элементов для удобства работы хранятся как список вершин, так и список граней. Топология сетки ставит в соответствие сеточным элементам (элементам и граням) набор составляющих их узлов (вершин) сетки. Обратная топология для каждого узла сетки хранит содержащие его элементы, что позволяет организовать быстрый доступ к соседним узлам выбранного узла и ко всем элементам, содержащим узел.
В скоростном пространстве расчетная область задается списком узлов и весами квадратурной формулы. Как правило, используется разбиение области интегрирования на гекса-
эдральные ячейки, однако, при использовании цилиндрической системы координат в сетку могут входить призмы.
Препроцессор. Препроцессор пакета используется для разбиения пространственных сеток на блоки для проведения параллельных вычислений на многоблочных сетках. Распределение ячеек исходной одноблочной неструктурированной сетки между процессорами производится с помощью программы Metis [35]. Используя это разбиение, препроцессор распределяет ячейки сетки между процессорами, проводит перенумерацию всех элементов сетки (вершин, граней и ядер), обеспечивает необходимое для проведения расчетов перекрытие между блоками и создает списки ячеек для обмена данными между всеми процессорами. Вся необходимая информация записывается в виде отдельных блоков так, что при проведении параллельных расчетов каждая нить считывает свою часть сетки. Благодаря этому параллельная и последовательная версии решателя практически совпадают.
Кинетический решатель. Кинетический решатель является надстройкой над ядром и реализует неявную разностную схему решения кинетического уравнения. Основные расчетные величины (макроскопические параметры и функция распределения) задаются в центрах ячеек пространственной сетки. Данные представляются в виде блочных векторов размерности Nspace по M элементов в блоке, где N — число элементов и M — число величин, заданных в элементе. Например, поля физических переменных представляются в виде блочного вектора с размером блока, равным 8: плотность, три компоненты скорости, температура и три компоненты вектора потока тепла. Для функции распределения длина блока равна числу узлов скоростной сетки Ng.
Помимо значений макроскопических величин и функции распределения в центрах ячеек пространственной сетки, кинетический решатель хранит в оперативной памяти значения функции распределения на гранях сетки (по два значения на грань) и правую часть R^. Таким образом, основные затраты оперативной памяти для хранения решения в байтах оцениваются как
TotalMemory = (2NspaceN + 2NfeCeNf + 8Nspace) RealSize,
где для вычислений с двойной точностью RealSize = 8.
Заметим, что для стационарных течений возможен также экономичный способ организации расчетов, при котором значения функции распределения на гранях не сохраняются в памяти ЭВМ.
В программной реализации разностной схемы алгоритм реконструкции применяется в каждой ячейке сразу для всего вектора f; результат записывается в массив значений f+, f- размера Nface. Далее в одномерном цикле по списку граней находится численный поток через каждую грань iface. Окончательно вклад в правую часть схемы R™ находится с учетом ориентации грани относительно ячейки. Для первого варианта неявной схемы (17)
величины emax находятся на этапе инициализации счета и хранятся в оперативной памяти, в то время как sign в (18) должны вычисляться для каждой ячейки в процессе решения линейной системы уравнений.
Организация параллельных вычислений с помощью MPI. Численное моделирование пространственных течений разреженного газа требует использования шестимерных расчетных сеток с числом узлов порядка нескольких миллиардов. Для ускорения вычислений метод решения должен быть реализован на многопроцессорных ЭВМ с разбиением расчетной сетки на блоки. При этом желательным является сохранение одного из основных достоинств метода — быстрой сходимости к стационарному режиму при использовании неявной дискретизации по времени.
Для кинетического уравнения возможны два подхода к построению параллельной версии численного метода. Первый подход является традиционным и предполагает разбиение пространственной сетки на блоки с некоторым перекрытием с учетом шаблона разностной схемы. Реализация явной разностной схемы на многоблочных неструктурированных сетках хорошо описана в литературе, см., например, [36]. Прямое обобщение неявного алгоритма (12) на многоблочные сетки сталкивается с серьезными трудностями как с алгоритмической точки зрения, так и в обеспечении хорошей масштабируемости метода по числу процессоров. С практической точки зрения наиболее предпочтительной является программная реализация метода без обмена данными между блоками на этапе решения линейной системы (20), (21). Получающийся при этом параллельный многоблочный алгоритм решения задачи не является строго эквивалентным одноблочному методу, но сходится к тому же стационарному решению.
Во втором подходе, впервые примененном для решения точного уравнения Больцмана с помощью явной разностной схемы в работе [37], используется разбиение скоростной сетки. При этом каждому процессору ставится в соответствие набор точек скоростного пространства, для которого строится решение кинетического уравнения с разностной схемы. На каждом шаге по времени вклады в интегральные суммы для вычисления макроскопических величин и плотности отраженных молекул в граничном условии складываются и рассылаются с помощью команды MPLAllReduce на каждый из процессоров. Получающийся параллельный алгоритм решения задачи полностью эквивалентен последовательному методу счета и для модельного уравнения прост в программной реализации.
В текущей версии «Несветай-3Д» реализованы оба подхода к организации расчетов с MPI, с некоторыми модификациями по сравнению с ранними версиями пакета. Для первого подхода используется более широкое перекрытие блоков сетки, с тем, чтобы пространственная реконструкция данных могла проводится и в непосредственных соседях к границе блока. В более ранних версиях пакета [ 17] использовалось минимально возможное перекрытие блоков, но при этом дополнительно требовался обмен данными на гранях граничных ячеек, что существенно усложняло программный код. Для второго подхода реализована гибкая
процедура обмена данными на этапе расчеты макроскопических величин: MPI_AllReduce вызывается только для тех ячеек, в которых итерационный процесс (16) не сошелся.
Сравнение двух подходов к организации параллельных вычислений проведено в [18, 38]. Показана хорошая масштабируемость кода на 4-х ядерных процессорах Intel Xeon при использовании при использовании до 512 вычислительных ядер.
4. Приложения
С помощью вычислительного пакета «Несветай-3Д» могут моделироваться как простые течения разреженного газа, так и течения в сложных геометрических областях. Благодаря реализованным в пакете алгоритмам повышенной точности и неявной дискретизации по времени, «Несветай-ЗД» позволяет строить численные решения задач динамики разреженного газа с хорошей точностью за приемлемое время счета. В настоящем разделе рассматриваются примеры расчетов течений в микроканалах и решения задач внешнего обтекания тел сложной формы.
Течения газа в микроканалах. Исследование течений разреженных газов в каналах и трубах важно для многих приложений. Современное состояние вопроса отражено в обзорах [7, 8]. К настоящему времени достаточно полно разработаны методы расчета течений в трубах и каналах бесконечной длины и произвольного поперечного сечения на основе линеаризованных кинетических моделей [39]. Одной из актуальных современных задач является нелинейная задача об истечении в вакуум. В статье [40] методом Монте-Карло изучены особенности такого течения для относительно коротких круглых труб. С увеличением относительной длины трубы проблема усложняется и требует все больших и больших затрат счетного времени. В недавней работе [41] истечение газа через короткую круглую трубу предложено в виде одной из стандартных верификационных задач, на которой можно сравнивать результаты счета по разным методам и комплексам программ.
Детальная постановка задачи о течении газа через трубу конечной длины изложена в работах [27, 40, 43]. Рассматривается стационарное течение разреженного газа из резервуара 1 в резервуар 2 через трубу длины L и кругового поперечного сечения с радиусом входного течения а. В общем случае радиус трубы может меняться по длине. Введем декартову систему координат с осью z, расположенной по оси трубы. Начало координат z = 0 выберем в центре входного сечения трубы, так что резервуары 1 и 2 занимают полупространства z < 0 и z > L соответственно. Давление и температуру в резервуарах 1 и 2 обозначим p10 и p20, T10 и T20; при этом p20 = 0. Боковые стенки резервуаров и трубы поддерживаются при постоянной температуре T10 = T20 = Tw, на них происходит полная аккомодация импульса и энергии падающих молекул и диффузное отражение их c равновесной функцией распределения при той же температуре Tw и с плотностью , определяемой условием непротекания. В расчетах значения давления, температуры и вязкости в резервуаре 1 используются как характерные значения p*, T*, в то время как радиус входного сечения
трубы а служит масштабом координаты Основной расчетной величиной является расход массы М через сечение трубы, отнесенный к тп10@1а2. При анализе результатов удобнее пользоваться величиной Q, представляющей собой расход М газа через сечение трубы, нормированный на величину расхода М°тв при свободно-молекулярном истечении в вакуум через отверстие [44, 45]:
Q = М/М 0отв, М 0отв = /п/2. (22)
Основным является случай круглой трубы постоянного радиуса, равного радиусу входного сечения а. В цикле работ [17, 27, 43, 46] решение задачи было построено пакетом «Несветай-3Д» для малого отношения давлений |р10/р20| — 1 ^ 1 (линеаризованная задача), а также двух нелинейных режимов р10/р20 = 2 (умеренный перепад давлений) и р10/р20 = то (течение в вакуум). Рассматривались трубы разной относительной длины 1 < Ь/а < 50.
Наиболее интересным является случай истечения в вакуум. В работе [18] была проведена серия верификационных расчетов с целью продемонстрировать точность и сходимость метода решения для короткой трубы Ь/а = 1. В расчетах использовалась последовательность из трех пространственных сеток. Первая сетка состояла из 5600 гексаэдров; при этом вдоль трубы использовалось 10 ячеек, в поперечном сечении — 116 ячеек. Шаг сетки по нормали к поверхности трубы равнялся 0. 06 единицам. Вторая и третья сетки получались последовательным удвоением числа ячеек по каждому направлению и состояли из 40761 и 349 905 ячеек соответственно. На рис. 2 представлены расчетная область в пространственных переменных и разрешение сетки. В скоростном пространстве использовалась расчетная область в цилиндрической системе координат с числом узлов 25 х 32 х 32 для 6 < 1. Таким образом, число узлов в шестимерной расчетной сетке может достигать 4.5 х 109.
I
Рис. 2. Пространственные сетки для течения в короткой круглой трубе
В табл. 1 приведены результаты расчетов пакетом «Несветай-3Д», методом ББМС [40] и экспериментальные данные [42]. Диапазон изменения параметра разреженности 0 < 8 < < 500 соответствует переходу от свободномолекулярного режима течения до практически континуального. Видна быстрая сходимость расчетов расхода массы по пространственной сетке. Приемлемая точность для всех 8 обеспечена уже на второй, относительно грубой, сетке. Отклонение от результатов счета по методу Монте-Карло не превышает 2% и достигается при 8 = 1; для остальных значений параметра разреженности согласие лучше. Расхождение с экспериментальными данными составляет 2-3% при оценке ошибки измерений в 2%.
Т а б л и ц а 1
Сходимость по пространственной сетке для нормированного расхода массы
5 «Несветай-3Д» Б8МС [40] Эксперимент [42]
(1) (11) (111)
0.0 0.666 0.670 0.672 0.672
0.1 0.678 0.683 0.684 0.680 0.675
1.0 0.758 0.766 0.768 0.754 0.743
10.0 1.035 1.061 1.066 1.062 1.06
100.0 1.290 1.351 1.367 1.358 1.33
200.0 1.331 1.406 1.425 1.412
500.0 1.331 1.454 1.474 1.449
Примечание. Нормированный расход массы Q определятся согласно (22); используется схема (14); варианты (1), (и) и (ш) соответствуют сеткам с 5600, 40761 и 349905 гексаэдрами соответственно.
В цитируемых работах также было построено решение задачи для Ь/а = 10, 20, 50. С точки зрения верификации метода наиболее важным является случай Ь/а = 10, для которого есть решение методом ББМС. В цитируемых работа показано отличное совпадение результатов счета по Б-модельному кинетическому уравнению и метода БМБС; расхождение между расчетами составляет не более 1%.
Помимо трубы постоянного радиуса, решения задачи построено для конической трубы для различных значений радиусов входного и выходного сечений [47], а также для составной круглой трубы [48], состоящей из двух частей одинаковой длины постоянного радиуса; при этом отношение радиусов частей равно 2. Наиболее интересным с точки зрения возникающей картины течения является случай расширяющейся составной трубы, в которой вторая секция, примыкающая к вакуумной камере, имеет вдвое больший радиус по сравнению с первой. На рис. 3 показана плоскость симметрии пространственной сетки для данного случая. Интересной особенностью течения в такой трубе является формирование при 8 ^ 1 во второй секции диска Маха и зоны возвратного течения. На рис. 4 показаны линии уровня числа Маха и линии тока для 8 = 1000. В области расширяющегося течения газа на оси трубы макроскопические величины следуют по универсальной кривой, слабо зависящей от 8. При этом в течении присутствует волна разрежения Прандтля — Майера, которая отражается от
Рис. 3. Сечение пространственной сетки для течения в составной трубе длины Ь = 10а
0 5 2 10
в
Рис. 4. Картина течения в расширяющейся составной трубе для 5 = 1000 и Ь/а = 10: а — числовая плотность; б — число Маха; в — линии тока
второй части трубы вниз по течению. С ростом значения параметра разреженности (уменьшения числа Кнудсена) структура ударной волны становится более выраженной и позиция маховского диска устанавливается около г ~ 9.
Задачи внешнего обтекания тел сложной формы. В качестве иллюстрации применения пакета программ к задачам внешней аэродинамики построим поле течения вокруг двух геометрических моделей [ 18, 49].
Первая модель представляет собой модель ракеты с двумя крыльями и четырьмя стабилизаторами. Геометрия модели подробно описана в [50]. Диаметр ракеты и значения давления и температуры в набегающем потоке выбираются в качестве масштабов длины /*, р*, Т*. Температура поверхности принималась равной температуре набегающего потока (холодное тело). Расчеты проводились для числа Маха набегающего потока М= 3, параметра разреженности 8 =1 и угла атаки 10 градусов. В пространственных переменных использовалась гибридная тетра-призматической сетка, что позволило аккуратно разрешить все основные особенности течения при минимальных трудозатратах на ее построение. Общее число пространственных ячеек Ж,расе = 3.3 х 105, из которых 7.4 х 103 ячеек — призмы. Скоростная сетка состояла из 25 х 16 х 32 узлов. Таким образом, общее число ячеек в 6-мерной расчетной сетке равнялось « 4.2 х 109.
Пространственная сетка разбивалась на 128 блоков; в каждом блоке было по 2580 ячеек. На рис. 5 показаны детали сетки вокруг тела, включая призматический слой на поверхности и разбиение на блоки. В расчетах использовалась схема второго порядка с полностью трехмерным шаблоном (13). Число Куранта принималось равным 25. Рис. 6 дает общее представление о картине течения. Представлены изолинии плотности в плоскостях х-г и у-г. Видна типичная картина обтекания холодного тела сверхзвуковым потоком разреженного
Рис. 5. Поверхностная сетка из треугольников с разбиением на блоки (слева) и сечение объемной сетки с призматическим слоем (справа) вокруг модели ракеты
Рис. 6. Линии уровня плотности в плоскостях х — г (слева) и у — г (справа)
газа. Течение характеризуется гладко меняющимися профилями плотности и температуры на линии торможения перед телом, а так же резким падение плотности и давление в следе. Вдоль поверхности все величины меняются достаточно медленно за исключением области крыльев и стабилизаторов. Следует так же отметить, что небольшая потеря симметрии течения в плоскости y-z является следствием использование достаточно грубой и несимметричной тетраэдральной сетки.
Вторая модель представляется собой Воздушно-Космический Аппарат (ВКА) ЦАГИ [51, 52]. Модель ВКА имеет сложную форму и состоит из фюзеляжа с затупленным носком, двух крыльев, вертикального киля и щитка. Общая длина тела со щитком — 10 м. Общая длина ВКА и значения давления и температуры в набегающем потоке выбираются в качестве масштабов длины /*, p*, T*. Температура поверхности принималась равной температуре набегающего потока (холодное тело). Расчеты проводились для числа Маха набегающего потока M^ = 2, параметра разреженности принимался равным ö = 1000, что приближенно соответствует высоте полета 100 км (длина свободного пробега 1 см).
В расчетах использовалась гибридная пространственная сетка со сгущением к носу тела и поверхностям большой кривизны. Вблизи поверхности вводился призматический слой толщиной 4 ячейки, высота каждой ячейки — 20 мм. Остальная часть расчетной области заполнялась тетраэдрами. Сетка разбивалась на 256 блоков. На рис. 7 приведена поверхностная сетка. Общее число пространственных ячеек равнялось Nspace ~ 533 х 103, из которых 386 х 103 — тетраэдры, 147 х 103 — призмы и 551 — пирамиды. Каждый из 256 блоков состоял из 2100 внутренних ячеек и 900. .. 2800 фиктивных ячеек, которые требовались для ТВД-схемы (13). Скоростная сетка состояла из 243 узлов. Общее количество ячеек в 6-мерной сетке равнялось « 6.9 х 109.
На рис. 8 показаны уровни плотности и давление в плоскости симметрии. Распределение давления по поверхности тела представлено на рис. 9. Как и в случае обтекания ракеты, видна типичная картина течения вокруг затупленного холодного тела. Вблизи носка тела начинает
Рис. 7. Расчетная сетка на поверхности ВКА с границами блоков
Рис. 8. Линии уровня плотности (слева) и давления (справа) в плоскости х-у
У
Рис. 9. Распределение давления по поверхности ВКА
формироваться отошедшая ударная волна. Поскольку локальный параметр разреженности, посчитанный по радиусу закругления носка достаточно мал 8П08е ~ 45, фронт волны сильно размыт. Для выбранного режима обтекания за тело образуется зона возвратного течения (не показана).
Заключение
В работе приведен обзор программного комплекса «Несветай-3Д», предназначенного для моделирования пространственных течений одноатомного разреженного газа. Основной комплекса является неявная конечно-объемная ТВД схема типа Годунова, реализованная на произвольных неструктурированных сетках. С помощью «Несветай-3Д»можно решать как исследовательские фундаментальные задачи механики разреженных газов, так и моделировать течения разреженного газа в реальных промышленных приложениях.
В настоящее время ведется работа по реализации в пакете неявной разностной схемы для нестационарных течений разреженного газа, адаптация для эффективного моделирования задач гиперзвукового обтекания тел сложной пространственной формы и обобщение на случай двухатомного газа на основе модели Рыкова [10, 11, 12].
Работа выполнена при поддержке РФФИ (проекты 13-01-00522-а и 14-08-00604-а). Для расчетов использовались суперкомпьютеры Университета Кренфилда (Великобритания), НИВЦ МГУ им. М.В. Ломоносова и МФТИ.
Список литературы
1. Шахов Е.М. Об обобщении релаксационного кинетического уравнения Крука // Известия АН СССР. Механика жидкости и газа. 1968. № 5. С. 142-145.
2. Шахов Е.М. Метод исследования движений разреженного газа. М.: Наука, 1974. 207 с.
3. Жук В.И., Рыков В.А., Шахов Е.М. Кинетические модели и задача о структуре ударной волны//Известия АН СССР. Механика жидкости и газа. 1973. №4. С. 135-141.
4. Рыков В.А., Черемисин Ф.Г., Шахов Е.М. Численные исследования по динамике разреженных газов // Журнал вычислительной математики и математической физики. 1980. Т. 20, №5. С. 1266-1283.
5. Шахов Е.М. Поперечное обтекание пластины разреженным газом // Известия АН СССР. Механика жидкости и газа. 1972. № 6. С. 106-113.
6. Шахов Е.М. Численные методы решения аппроксимирующих кинетических уравнений //Численные методы в динамике разреженных газов: сб. Вып. 2. М.: ВЦ АН СССР, 1975. С. 35-76.
7. Sharipov F., Seleznev V. Data on internal rarefied gas flows // J. Phys. Chem. Ref. Data. 1998. V. 27, no. 3. P. 657-706.
8. Шарипов Ф.М., Селезнев В.Д. Движение разреженных газов в каналах и микроканалах. Екатеринбург: УРО РАН, 2008. 232 с.
9. Graur I., Polikarpov A.Ph. Comparison of different kinetic models for the heat transfer problem // Heat and Mass Transfer. 2009. Vol. 46. P. 237-244. DOI: 10.1007/s00231-009-0558-x
10. Рыков В.А. Модельное кинетическое уравнение для газа с вращательными степенями свободы // Известия АН СССР. Механика жидкости и газа. 1975. № 6. С. 107-115.
11. Рыков В.А., Титарев В.А., Шахов Е.М. Структура ударной волны в двухатомном газена основе кинетической модели // Известия РАН. Механика жидкости и газа. 2008. №2. С.171-182.
12. Ларина И.Н., Рыков В.А. Расчет течений разреженного двухатомного газа через плоский микроканал // Журнал вычислительной математики и математической физики. 2012. Т. 52, №4. С. 720-732.
13. Li Z.-H., Zhang H.-X. Numerical investigation from rarefied flow to continuum by solving the Boltzmann model equation // International Journal for Numerical Methods in Fluids. 2003. Vol.42, no. 4. P. 361-382.
14. Kolobov V.I., Arslanbekov R.R., Aristov V.V., Frolova A.A., Zabelok S.A. Unified solver for rarefied and continuum flows with adaptive mesh and algorithm refinement // Journal of Computational Physics. 2007. Vol. 223, no. 2. P. 589-608. DOI: 10.1016/j.jcp.2006.09.021
15. Аникин Ю.А., Клосс Ю.Ю., Мартынов Д.В., Черемисин Ф.Г. Компьютерное моделирование и анализ эксперимента Кнудсена 1910 года // Нано- и микросистемная техника. 2010, №8. С. 6-14.
16. Титарев В.А. Неявный численный метод расчета пространственных течений разреженного газа на неструктурированных сетках // Журнал вычислительной математики и математической физики. 2010. Т. 50, № 10. С. 1811-1826.
17. Titarev V.A. Efficient deterministic modelling of three-dimensional rarefied gas flows // Commun. Comput. Phys. 2012. Vol. 12, no. 1. P. 161-192.
18. Titarev V.A., Dumbser M., Utyuzhnikov S.V. Construction and comparison of parallel implicit kinetic solvers in three spatial dimensions // Journal of Computational Physics. 2014. Vol. 256, no. 1. P. 17-33. DOI: 10.1016/j.jcp.2013.08.051
19. Bhatnagar P.L., Gross E.P., Krook M. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems //Phys. Rev. 1954. Vol. 94, no. 3. P. 511-525.
20. Годунов С.К. Разностный метод численного расчета разрывных решений уравнений гидродинамики//Математический сборник. 1959. Т. 47, №3. С. 271-306.
21. Колган В.П. Применение принципа минимальных значений производной к построению конечно-разностных схем для расчета разрывных течений газовой динамики // Ученые записки ЦАГИ. 1972. Т. 3, №6. С. 68-77.
22. Toro E.F. Riemann solvers and numerical methods for fluid dynamics. 3rd ed. Springer Berlin Heidelberg, 2009. 724 p. DOI: 10.1007/b79761
23. Dumbser M., Kaser M. Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems // Journal of Computational Physics. 2007. Vol. 221, no. 2. P. 693-723. DOI: 10.1016/j.jcp.2006.06.043
24. Tsoutsanis P., Titarev V.A., Drikakis D. WENO schemes on arbitrary mixed-element unstructured meshes in three space dimensions// Journal of Computational Physics. 2010. Vol. 230, no. 4. P. 1585-1601. DOI: 10.1016/j.jcp.2010.11.023
25. Venkatakrishnan V. On the accuracy of limiters and convergence to steady-state solutions // 31st Aerospace Science Meeting and Exhibit, Reno, NV, January 11-14, 1993. AIAA paper 93-0880. DOI: 10.2514/6.1993-880
26. Петровская Н.Б., Волков А.В. Влияние геометрии сетки на точность реконструкции решения в конечно-объемных и конечно-элементных схемах высокого порядка // Математическое моделирование. 2010. T. 22, №3. С. 145-160.
27. Titarev V.A., Shakhov E.M. Computational study of a rarefied gas flow through a long circular pipe into vacuum // Vacuum. Spec. iss. "Vacuum Gas Dynamics. Theory, experiments and practical applications". 2012. Vol.86, no. 11. P. 1709-1716. DOI: 10.1016/j.vacuum.2012.02.026
28. Leer B. van. Towards the ultimate conservative difference scheme. V. A second order sequel to Godunov's method // Journal of Computational Physics. 1979. Vol. 32, iss. 1. P. 101-136. DOI: 10.1016/0021-9991(79)90145-1
29. Titarev V.A. Conservative numerical methods for model kinetic equations // Computers and Fluids. 2007. Vol. 36, no. 9. P. 1446-1459.
30. Титарев В.А., Шахов Е.М. Численный расчет поперечного обтекания холодной пластины гиперзвуковым потоком разреженного газа // Известия РАН. Механика жидкости и газа. 2005. №5. C. 139-154.
31. Титарев В.А. Численный метод расчета двухмерных нестационарных течений разреженного газа в областях произвольной формы // Журнал вычислительной математики и математической физики. 2009. T. 49, №7. С. 1255-1270.
32. Mieussens L. Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics // Math. Models and Meth. Appl. Sci. 2000. Vol. 8, no. 10. P. 1121-1149.
33. Gusarov A.V., Smurov I. Gas-dynamic boundary conditions of evaporation and condensation. С. Numerical analysis oftheKnudsenlayer//Phys. Fluids. 2002. Vol. 14, no. 12. P. 4242-4255.
34. Men'shov I.S., Nakamura Y. An implicit advection upwind splitting scheme for hypersonic air flows in thermochemical nonequilibrium // A Collection of Technical Papers of 6th Int. Symp. on CFD, Lake Tahoe, Nevada, 1995. Vol. 2. 1995. P. 815.
35. Karypis G., Kumar V. Multilevel k-way partitioning scheme for irregular graphs // J. Parallel Distrib. Comput. 1998. Vol. 48. P. 96-129.
36. Dumbser M., Kaser M., Titarev V. A., Toro E.F. Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems // Journal of Computational Physics. 2007. Vol. 226, iss. 1. P. 204-243. DOI: 10.1016/j.jcp.2007.04.004
37. Аристов В.В., Забелок C.A. Детерминистический метод решения уравнения Больцмана с параллельными вычислениями // Журнал вычислительной математики и математической физики. 2002. Т. 42, №3. С. 425-437.
38. Думбсер М., Титарев В.А., Утюжников С.В. Неявный многоблочный метод решения кинетического уравнения на неструктурированных сетках // Журнал вычислительной математики и математической физики. 2013. Т. 53, №5. С. 762-782. DOI: 10.7868/S0044466913050141
39. Титарев В.А., Шахов Е.М. Неизотермическое течение газа в длинном канале на основе кинетической S-модели // Журнал вычислительной математики и математической физики. 2010. Т. 50, № 12. С. 2246-2260.
40. Varoutis S., Valougeorgis D., Sazhin O., Sharipov F. Rarefied gas flow through short tubes into vacuum // Journal of Vacuum Science and Technology A. 2008. Vol. 26, no. 1. P. 228-238. DOI: 10.1116/1.2830639
41. Sharipov F. Benchmark problems in rarefied gas dynamics // Vacuum. Spec. iss. "Vacuum Gas Dynamics. Theory, experiments and practical applications". 2012. Vol.86, no. 11. P. 16971700. DOI: 10.1016/j.vacuum.2012.02.048
42. Fujimoto T., Usami M. Rarefied gas flow through a circular orifice and short tubes // Journal of Fluids Engineering. 1984. Vol. 106, no. 4. P. 367-373. DOI: 10.1115/1.3243132
43. Титарев В.А., Шахов Е.М. Концевые эффекты при истечении разреженного газа через длинную трубу в вакуум // Известия РАН. Механика жидкости и газа. 2013. № 5. С. 146158.
44. Шахов Е.М. Осесимметричная нелинейная задача о стационарном течении разреженного газа в трубе кругового сечения // Журнал вычислительной математики и математической физики. 1996. Т. 36, №8. С. 169-179.
45. Varoutis S., Valougeorgis D., Sharipov F. Simulation of gas flow through tubes of finite length over the whole range of rarefaction for various pressure drop ratios // Journal of Vacuum Science and Technology A. 2009. Vol. 27, no. 6. P. 1377-1391. DOI: 10.1116/1.3248273
46. Titarev V.A. Rarefied gas flow in a circular pipe of finite length // DOI: 10.1016/j.vacuum.2013.01.012
47. Titarev V.A., Shakhov E.M., Utyuzhnikov S.V. Rarefied gas flow through a diverging conical pipe into vacuum // Vacuum. 2014. Vol. 101. P. 10-17. DOI: 10.1016/j.vacuum.2013.07.030
48. Titarev V.A., Shakhov E.V. Rarefied gas flow into vacuum through a pipe composed of two circular sections of different radii // Vacuum. 2014. DOI: 10.1016/j.vacuum.2014.02.019 (in press).
49. Titarev V.A. Direct numerical solution of model kinetic equations for flows in arbitrary three-dimensional geometries // Proc. 28th Int. Symp. on Rarefied Gas Dynamics 2012. AIP Conference Proceedings 1501, Amer. Inst. Physics, 2012. P. 262-271. DOI: 10.1063/1.4769518
50. Garanzha V.A., Kudryavtseva L.N., Utyuzhnikov S.V. Variational method for untangling and optimization of spatial meshes // J. Comp. and Appl. Math. 2014. Vol. 269. P. 24-41.
51. Ваганов А.В., Дроздов С.М., Косых А.П., Нерсесов Г.Г., Челышева И.Ф., Юмашев В.Л. Численное моделирование аэродинамики крылатого возвращаемого космического аппарата // Ученые записки ЦАГИ. 2009. Т. 40, № 2. С. 3-15.
52. Ваганов А.В., Дроздов С.М., Задонский С.М., ЧелышеваИ.Ф., Косых А.П., Нерсесов Г.Г., Юмашев В.Л. Исследование аэродинамики крылатого воздушно-космического аппарата с отклоненным балансировочным щитком // Ученые записки ЦАГИ. 2009. Т. 40, №5. С. 3-15.
SCIENTIFIC PERIODICAL OF THE BAUMAN MSTU
SCIENCE and EDUCATION
EL № FS77 - 48211. №0421200025. ISSN 1994-0408
electronic scientific and technical journal
Software Package \Nesvetay-3D"
for modeling three-dimensional flows
of monatomic rarefied gas
# 06, June 2014
DOI: 10.7463/0614.0712314
Titarev V. A.
Dorodnycin Computing Centre of RAS 105005, Moscow, Russian Federation
titarev@ccas.ru
Analysis of three-dimensional rarefied gas flowsin microdevices (micropipes, micropumps etc) and over re-entry vehicles requires development of methods of computational modelling. One of such methods is the direct numerical solution of the Boltzmann kinetic equation for the velocity distribution function with either exact or approximate (model) collision integral. At present, for flows of monatomic rarefied gas the Shakhov model kinetic equation, also called S-model, has gained wide-spread use. The equation can be regarded as a model equation of the incomplete third-order approximation. Despite its relative simplicity, the S-model is still a complicated integro-differential equation of high dimension. The numerical solution of such an equation requires high-accuracy parallel methods.
The present work is a review of recent results concerning the development and application of three-dimensional computer package Nesvetay-3D intended for modelling of rarefied gas flows. The package solves Boltzmann kinetic equation with the BGK (Krook) and Shakhov model collision integrals using the discrete velocity approach. Calculations are carried out in non-dimensional variables. A finite integration domain and a mesh are introduced in the molecular velocity space. Next, the kinetic equation is re-written as a system of kinetic equations for each of the discrete velocities. The system is solved using an implicit finite-volume method of Godunov type. The steady-state solution is computed by a time marching method. High order of spatial accuracy is achieved by using a piece-wise linear representation of the distribution function in each spatial cell. In general, the coefficients of such an approximation are found using the least-square method. Arbitrary unstructured meshes in the physical space can be used in calculations, which allow considering flows over objects of general geometrical shape. Conservative property of the method with respect to the model collision integral is guaranteed by means of a special procedure to calculate macroscopic variables. Another important part of the numerical method is the fast solution of the
linear system for time increments of the distribution function. The solution is based on the LU-SGS approach so that the number of operations is linearly proportional to the number of cells in the spatial mesh. Large problems can be solved on hundreds of CPU cores using the Message Passing Interface (MPI).
Performance and robustness of the numerical method and computer code are illustrated on a number of problems, including rarefied gas flows through microchannels into vacuum and external flows over re-entry vehicles on the high altitude of flight. Rarefied gas flows through simple and composite channels of circular cross sectional area are considered. Comparisons with the results of other authors and experimental data are shown.Good spatial mesh convergence of the method is demonstrated. For the flow in a composite channel the formation of a Mach disk is shown. One of the examples of external flow calculation is the analysis of rarefied gas flow over model winged re-entry space vehicle (RSV), which is a three-dimensional object of a very complex shape. The RSV is proposed by Central Aerohydrodynamic Institute (TsAGI). Its geometry includes a blunt fuselage, sweptwings, keel and flap. The flow pattern over the RSV and surface pressure distribution are shown.
At present the work is under way to incorporate an implicit numerical scheme for time-dependent rarefied gas flows, an adaptation algorithm for the efficient modelling of hypersonic flows over complex objects and extension to the diatomic gas modelling on the bases of the Rykov kinetic model.
Publications with keywords: parallelism, rarefied gas, S-model, kinetic equation, implicit TVD method
References
1. Shakhov E.M. [Generalization of the Krook kinetic relaxation equation]. Izvestiya AN SSSR Mekhanika Zhidkosti i Gaza, 1968, no. 5, pp. 142-145. (English translation: Fluid Dynamics, 1968, vol. 3, iss. 5, pp. 95-96. DOI: 10.1007/BF01029546).
2. Shakhov E.M. Metod issledovaniya dvizheniy razrezhennogo gaza [Method of investigation of rarefied gas flows]. Moscow, Nauka Publ., 1974. 207 p. (in Russian).
3. Zhuk V.I., Rykov V.A., Shakhov E.M. [Kinetic models and the shock structure problem]. Izvestiya ANSSSR. Mekhanika Zhidkosti i Gaza, 1973, no. 4, pp. 135-141. (English translation: Fluid Dynamics, 1973, vol. 8, iss. 4, pp. 620-625. DOI: 10.1007/BF01013101).
4. Rykov V.A., Cheremisin F.G., Shakhov E.M. [Numerical investigations in rarefied gas dynamics]. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki, 1980, vol.20, no. 5, pp. 1266-1283. (English translation: USSR Computational Mathematics and Mathematical Physics, 1980, vol. 20, iss. 5, pp. 168-185. DOI: 10.1016/0041-5553(80)90094-4).
5. Shakhov E.M. [Transverse flow of a rarefield gas around a plate]. Izvestiya ANSSSR. Mekhani-ka Zhidkosti i Gaza, 1972, no. 6, pp. 106-113. (English translation: Fluid Dynamics, 1972, vol. 7, iss. 6, pp. 961-966. DOI: 10.1007/BF01176114).
6. Shakhov E.M. [Numerical methods of solution of approximating kinetic equations]. Chislen-nye metody v dinamike razrezhennykh gazov: sb. Vyp. 2 [Numerical methods in rarefied gas dynamics: coll. papers. Iss. 2]. Moscow, Computing Center of the USSR Academy of Sciences Publ., 1975, pp. 35-76. (in Russian).
7. Sharipov F., Seleznev V. Data on internal rarefied gas flows. J.Phys. Chem. Ref. Data, 1998, vol. 27, no. 3, pp. 657-706.
8. Sharipov F.M., Seleznev V.D. Dvizhenie razrezhennykh gazov v kanalakh i mikrokanalakh [Flows of rarefied gases in channels and microchannels]. Ekaterinburg, Russian Academy of Science, Ural Branch, Institute of Thermal Physics Publ., 2008. 232 p. (in Russian).
9. Graur I., Polikarpov A.Ph. Comparison of different kinetic models for the heat transfer problem. Heat and Mass Transfer, 2009, vol. 46, iss. 2, pp. 237-244. DOI: 10.1007/s00231-009-0558-x
10. Rykov V.A. [A model kinetic equation for a gas with rotational degrees of freedom]. Izvestiya ANSSSR Mekhanika Zhidkosti i Gaza, 1975, no. 6, pp. 107-115. (English translation: Fluid Dynamics, 1975, vol. 10, iss. 6, pp. 959-966. DOI: 10.1007/BF01023275).
11. Rykov V.A., Titarev V.A., Shakhov E.M. [Shock wave structure in a diatomic gas based on a kinetic model]. Izvestiya RAN. Mekhanika Zhidkosti i Gaza, 2008, no. 2, pp. 171182. (English translation: Fluid Dynamics, 2008, vol.43, iss. 2, pp. 316-326. DOI: 10.1134/S0015462808020178).
12. Larina I.N., Rykov V.A. [Computation of rarefied diatomic gas flows through a plane microchannel]. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki, 2012, vol. 52, no. 4, pp. 720-732. (English translation: Computational Mathematics and Mathematical Physics, 2012, vol. 52, no. 4, pp. 637-648. DOI: 10.1134/S0965542512040112).
13. Li Z.-H., Zhang H.-X. Numerical investigation from rarefied flow to continuum by solving the Boltzmann model equation. International Journal for Numerical Methods in Fluids, 2003, vol.42, no. 4, pp. 361-382.
14. Kolobov V.I., Arslanbekov R.R., Aristov V.V., Frolova A.A., Zabelok S.A. Unified solver for rarefied and continuum flows with adaptive mesh and algorithm refinement. Journal of Computational Physics, 2007, vol. 223, no. 2, pp. 589-608. DOI: 10.1016/j.jcp.2006.09.021
15. Anikin Yu.A., Kloss Yu.Yu., Martynov D.V., Cheremisin F.G. [Computer Simulation and Analysis of the Knudsen experiment of the 1910 year]. Nano- i mikrosistemnaya tekhnika { Journal of NANO and MICROSYSTEM TECHNIQUE, 2010, no. 8, pp. 6-14.
16. Titarev V.A. [Implicit numerical method for computing three-dimensional rarefied gas flows on unstructured meshes]. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki, 2010, vol. 50, no. 10, pp. 1811-1826. (English translation: Computational Mathematics and Mathematical Physics, 2010, vol. 50, no. 10, pp. 1719-1733. DOI: 10.1134/S0965542510100088).
17. Titarev V.A. Efficient deterministic modelling of three-dimensional rarefied gas flows. Commun. Comput. Phys., 2012, vol. 12, no. 1, pp. 161-192.
18. Titarev V.A., Dumbser M., Utyuzhnikov S.V. Construction and comparison of parallel implicit kinetic solvers in three spatial dimensions. Journal of Computational Physics, 2014, vol. 256, no. 1, pp. 17-33. DOI: 10.1016/j.jcp.2013.08.051
19. Bhatnagar P.L., Gross E.P., Krook M. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev., 1954, vol. 94, no. 3, pp. 511-525.
20. Godunov S.K. [A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics]. Matematicheskiy sbornik, 1959, vol. 47, no. 3, pp. 271-306. (English translation: Sbornik: Mathematics, 1959, vol. 47, pp. 357-393).
21. Kolgan V.P. [Application of the minimum-derivative principle in the construction of finite-difference schemes for numerical analysis of discontinuous solutions in gas dynamics]. Uchenye zapiski TsAGI, 1972, vol. 3, no. 6, pp. 68-77. (in Russian).
22. Toro E.F. Riemann solvers and numerical methods for fluid dynamics. 3rd ed. Springer Berlin Heidelberg, 2009. 724 p. DOI: 10.1007/b79761
23. Dumbser M., Kaser M. Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems. Journal of Computational Physics, 2007, vol. 221, no. 2, pp. 693-723. DOI: 10.1016/j.jcp.2006.06.043
24. Tsoutsanis P., Titarev V.A., Drikakis D. WENO schemes on arbitrary mixed-element unstructured meshes in three space dimensions. Journal of Computational Physics, 2010, vol.230, no. 4, pp. 1585-1601. DOI: 10.1016/j.jcp.2010.11.023
25. Venkatakrishnan V. On the accuracy of limiters and convergence to steady-state solutions. 31st Aerospace Science Meeting and Exhibit, Reno, NV, January 11-14, 1993. AIAA paper 93-0880. DOI: 10.2514/6.1993-880
26. Petrovskaya N.B., Volkov A.V. [The impact of grid geometry on the accuracy of higher order finite-volume and finite-element schemes]. Matematicheskoe modelirovanie, 2010, vol.22, no. 3, pp. 145-160. (in Russian).
27. Titarev V. A., Shakhov E.M. Computational study of a rarefied gas flow through a long circular pipe into vacuum. Vacuum. Spec. iss. "Vacuum Gas Dynamics. Theory, experiments and practical applications", 2012, vol. 86, no. 11, pp. 1709-1716. DOI: 10.1016/j.vacuum.2012.02.026
28. Leer B. van. Towards the ultimate conservative difference scheme V. A second order sequel to Godunov's method. Journal of Computational Physics, 1979, vol. 32, iss. 1, pp. 101-136. DOI: 10.1016/0021-9991(79)90145-1
29. Titarev V.A. Conservative numerical methods for model kinetic equations. Computers and Fluids, 2007, vol. 36, no. 9, pp. 1446-1459.
30. Titarev V.A., Shakhov E.M. [Numerical Calculation of the Hypersonic Rarefied Transverse Flow Past a Cold Flat Plate]. Izvestiya RAN. Mekhanika Zhidkosti i Gaza, 2005, no. 5, pp. 139-154. (English translation: Fluid Dynamics, 2005, vol.40, iss. 5, pp. 790-804. DOI: 10.1007/s10697-005-0117-1).
31. Titarev V.A. [Numerical method for computing two-dimensional unsteady rarefied gas flows in arbitrarily shaped domains]. Zhurnal vychislitel'noy matematiki i matematicheskoyfiziki, 2009, vol.49, no. 7, pp. 1255-1270. (English translation: Computational Mathematics and Mathematical Physics, 2009, vol. 49, iss. 7, pp. 1197-1211. DOI: 10.1134/S0965542509070112).
32. Mieussens L. Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics. Math. Models andMeth. Appl. Sci., 2000, vol. 8, no. 10, pp. 1121-1149.
33. Gusarov A.V., Smurov I. Gas-dynamic boundary conditions of evaporation and condensation: Numerical analysis of the Knudsen layer. Phys. Fluids, 2002, vol. 14, no. 12, pp. 4242-4255. DOI: 10.1063/1.1516211
34. Men'shov I.S., Nakamura Y. An implicit advection upwind splitting scheme for hypersonic air flows in thermochemical nonequilibrium. A Collection of Technical Papers of 6th Int. Symp. on CFD, Lake Tahoe, Nevada, 1995. Vol. 2. 1995, p. 815.
35. Karypis G., Kumar V. Multilevel k-way partitioning scheme for irregular graphs. J. Parallel Distrib. Comput., 1998, vol.48, pp. 96-129.
36. Dumbser M., Kaser M., Titarev V. A., Toro E.F. Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems. Journal of Computational Physics, 2007, vol. 226, iss. 1, pp. 204-243. DOI: 10.1016/j.jcp.2007.04.004
37. Aristov V.V., Zabelok C.A. [A deterministic method for solving the Boltzmann equation with parallel computations]. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki, 2002, vol. 42, no. 3, pp. 425-437. (English translation: Computational Mathematics and Mathematical Physics, 2002, vol. 42, iss. 3, pp. 406-418).
38. Dumbser M., Titarev V.A., Utyuzhnikov S.V. [Implicit multiblock method for solving a kinetic equation on unstructured meshes]. Zhurnal vychislitel'noy matematiki i matematicheskoyfiziki, 2013, vol. 53, no. 5, pp. 762-782. DOI: 10.7868/S0044466913050141 (English translation: Computational Mathematics and Mathematical Physics, 2013, vol.53, iss. 5, pp. 601-615. DOI: 10.1134/S0965542513050126).
39. Titarev V.A., Shakhov E.M. [Nonisothermal gas flow in a long channel analyzed on the basis of the kinetic S-Model]. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki, 2010, vol. 50, no. 12, pp. 2246-2260. (English translation: Computational Mathematics and Mathematical Physics, 2010, vol. 50, iss. 12, pp. 2131-2144. DOI: 10.1134/S0965542510120146).
40. Varoutis S., Valougeorgis D., Sazhin O., Sharipov F. Rarefied gas flow through short tubes into vacuum. Journal of Vacuum Science and Technology A, 2008, vol. 26, no. 1, pp. 228-238. DOI 10.1116/1.2830639
41. Sharipov F. Benchmark problems in rarefied gas dynamics. Vacuum. Spec. iss. "Vacuum Gas Dynamics. Theory, experiments and practical applications", 2012, vol. 86, no. 11, pp. 16971700. DOI: 10.1016/j.vacuum.2012.02.048
42. Fujimoto T., Usami M. Rarefied gas flow through a circular orifice and short tubes. Journal of Fluids Engineering, 1984, vol. 106, no. 4, pp. 367-373. DOI: 10.1115/1.3243132
43. Titarev V.A., Shakhov E.M. [End effects in rarefied gas outflow into vacuum through a long tube]. Izvestiya RAN. Mekhanika Zhidkosti i Gaza, 2013, no. 5, pp. 146-158. (English translation: Fluid Dynamics, 2013, vol. 48, iss. 5, pp. 697-707. DOI: 10.1134/S001546281305013X).
44. Shakhov E.M. [The axisymmetric nonlinear steady flow of a rarefied gas in a pipe of circular cross-section]. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki, 1996, vol.36, no. 8, pp. 169-179. (English translation: Computational Mathematics and Mathematical Physics, 1996, vol.36, no. 8, pp. 1123-1131.).
45. Varoutis S., Valougeorgis D., Sharipov F. Simulation of gas flow through tubes of finite length over the whole range of rarefaction for various pressure drop ratios. Journal of Vacuum Science and Technology A, 2009, vol. 27, no. 6, pp. 1377-1391. DOI: 10.1116/1.3248273
46. Titarev V.A. Rarefied gas flow in a circular pipe of finite length. Vacuum, 2013, vol.94, pp. 92-103. DOI: 10.1016/j.vacuum.2013.01.012
47. Titarev V.A., Shakhov E.M., Utyuzhnikov S.V. Rarefied gas flow through a diverging conical pipe into vacuum. Vacuum, 2014, vol. 101, pp. 10-17. DOI: 10.1016/j.vacuum.2013.07.030
48. Titarev V.A., Shakhov E.V. Rarefied gas flow into vacuum through a pipe composed of two circular sections of different radii. Vacuum, 2014. DOI: 10.1016/j.vacuum.2014.02.019 (in press).
49. Titarev V.A. Direct numerical solution of model kinetic equations for flows in arbitrary three-dimensional geometries. Proc. 28th Int. Symp. on Rarefied Gas Dynamics 2012. AIP Conference Proceedings 1501, Amer. Inst. Physics, 2012, pp. 262-271. DOI: 10.1063/1.4769518
50. Garanzha V.A., Kudryavtseva L.N., Utyuzhnikov S.V. Variational method for untangling and optimization of spatial meshes. J. Comp. andAppl. Math., 2014, vol. 269, pp. 24-41.
51. Vaganov A.V., Drozdov S.M., Kosykh A.P., Nersesov G.G., Chelysheva I.F., Yumashev V.L. [Numerical simulation of aerodymics winged re-entry space vehicle]. Uchenye zapiski TsAGI, 2009, vol.40, no.2, pp. 3-15. (English translation: TsAGI Science Journal, 2009, vol.40, iss. 2, pp. 131-149. DOI: 10.1615/TsAGISciJ.v40.i2.10).
52. Vaganov A.V., Drozdov S.M., Zadonskiy S.M., Chelysheva I.F., Kosykh A.P., Nersesov G.G., Yumashev V.L. [Study of aerodynamics of the aerospace vehicle with a deflected balancing flap]. Uchenye zapiski TsAGI, 2009, vol.40, no. 5, pp. 3-15. (English translation: TsAGI Science Journal, 2009, vol. 40, iss. 5, pp. 517-534. DOI: 10.1615/TsAGISciJ.v40.i5.10).
