CC BY
34Математика й Математическое
моделирование
Ссылка на статью: // Математика и Математическое моделирование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 06. С. 61-77.
DOI: 10.7463/шаШш.0615.0823537
Представлена в редакцию: 08.11.2015 Исправлена: 22.11.2015
© МГТУ им. Н.Э. Баумана
УДК 519.634
Численное сравнение решений кинетических модельных уравнений
Фролова А. А.1' '' а а &о!о уа@у аш1ехди
вычислительный центр им. А.А. Дородницына РАН Федерального исследовательского центра «Информатика и управление» РАН,
Москва, Россия
Модель Шахова, дающая правильное число Прандтля при переходе к гидродинамическому режиму, является обобщением модели Бхатнагара - Гросса - Крука (БГК, BGK). Она получена разложением в ряд интеграла обратных столкновений по полиномам Эрмита до третьего порядка. Использование этого же разложения при другом выборе свободного параметра приводит к линеаризованной эллипсоидальной статистической модели (ESL), дающей также правильное число Прандтля. Точность ESL- модели при решении задач течения разреженного газа ранее не исследовалась. В работе проводится численное сравнение решений модельных уравнений с решением полного уравнения Больцмана для трех тестовых задач распада разрыва. При вычислениях используется консервативный метод дискретных ординат. Сравнение показывает, что модельные уравнения могут существенно искажать профили макропараметров перед фронтом ударной волны, при этом ESL-модель дает значительно большее отклонение от решения уравнения Больцмана, чем модель Шахова.
Ключевые слова: уравнение Больцмана, модельные уравнения, эллипсоидальная статистическая модель (ES) , модель Шахова ^-модель)
уравнение БГК,
Введение
В настоящее время значительно возрос интерес к моделированию различных физических явлений с использованием кинетического описания. Это связано как с необходимостью более детального изучения газодинамических течений в областях сильной неравновесности газа (например, в задачах микроэлектроники или аэрокосмических проблемах), так и с развитием вычислительных средств, позволяющих проводить расчеты широкого класса задач на базе кинетических уравнений.
Для решения полного уравнения Больцмана наиболее развитыми являются методы статистического моделирования (DSMC)[1] и прямого решения уравнения Больцмана с использованием дискретных ординат [2]. Каждый из этих методов имеет свои преимущества и свои недостатки. Несмотря на развитие и усовершенствование в последнее время
этих подходов [3 - 9], они остаются весьма трудоемкими. Особенно объем вычислений возрастает при решении методом дискретных ординат трехмерных задач в сверхзвуковых и гиперзвуковых режимах, а также при исследовании течений многокомпонентного газа.
В связи с этим представляется важным поиск и развитие различных альтернативных подходов, основанных на тех или иных упрощениях кинетического уравнения. Одним из таких упрощений является аппроксимация интеграла столкновений модельным релаксационным членом. В настоящее время разработано много различных моделей для интеграла столкновений. Наиболее популярными являются модель Бхатнагара - Гросса - Крука (БГК, BGK) [10], эллипсоидальная статистическая модель (ES) [11] , модель Шахова модель) [12] и модели с частотой, зависящей от молекулярной скорости [13-14].
Отметим, что если S-модель широко применялась с 70 годов и достаточно хорошо зарекомендовала себя (см. например, список литературы в [15]), то ES - модель стала широко использоваться только после доказательства для нее Н-теоремы [16], то есть сравнительно недавно, и все возможности ее еще не до конца исследованы.
В работе [17] для задач течения одноатомного газа и смесей газов предложена линеаризация ES- модели (в дальнейшем ESL). Указывается, что данная модель может быть более точной, чем S- модель, так как аккуратнее описывает "хвосты" функции распределения. Так как численных расчетов, подтверждающих это, не приводится, представляется важным провести тестирование ESL модели. В данной работе мы сравним решения, получаемые по модели Шахова и ESL- модели с решениями полного кинетического уравнения для некоторых задач газовой динамики.
1. Кинетическое уравнение
Для лучшего понимания излагаемого материала кратко опишем кинетическое уравнение и некоторые модельные уравнения.
Кинетическое уравнение Больцмана без учета действия силового члена может быть записано в виде [18]:
/ + 4■(!*/) = I(г,г,£) ,
где /(г, г) - функция распределения по скоростям, зависящая от пространственного вектора г = (х, у,е), вектора скорости £ = ) и времени Интеграл столкновений
I(г, г,£) в случае бинарных столкновений представляет собой интегральный локальный оператор в физическом пространстве и для упругих соударений одноатомного одноком-понентного газа записывается следующим образом:
I(4) =\ ^ | (/£)/(£' ) - /&)/($)) I я I *(I я I , . (1)
К3 52
Здесь (£, ^) - скорости частиц до соударений, (£ , £ 1) - скорости после соударений, связанные законами сохранения импульса и энергии, g = - вектор относительной ско-
рости сталкивающихся частиц и = $лпдс10(1£ - элемент телесного угла. Дифференциальное сечение столкновений а(g,в) зависит от относительной скорости и угла рассеивания в между векторами g и g = ^ . Интеграл (1) можно записать также в виде:
I(& = 0(& (& ,
где О(£) - интеграл обратных столкновений, а у(£) - частота столкновений.
Используя функцию распределения, макроскопические переменные, плотность числа частиц п, средняя скорость и, температура Т, компоненты тензора напряжений Р , вектор теплового потока q , давление р и компоненты тензора неравновесных напряжений рг/ определяются по формулам:
т
п =\ М, и = 1 (■ /С^, Т = -—- | с 2м,
: пJ 3квп2
с т с 1
Ру = т^с^М, ч = - ¡с 2сМ, р = ~(Рп + Р22 + Р33),
т
2 - - - з
Ру = Р — Р ■
Здесь с = \ — и - вектор относительной скорости молекул газа, т - масса молекул, кв -константа Больцмана. Нижние индексы г и ] принимают значения от 1 до 3 и обозначают компоненты величин вдоль осей х,у,7 соответственно.
Интеграла столкновений обладает рядом важных свойств. Первое из них - это ортогональность к сумматорным инвариантам щ(^) = (1, 2 У,
¡ I = 0 ■
я3
Данное свойство является следствием симметрии интеграла и выполнения законов сохранения (массы, импульса и энергии) при упругом столкновении частиц. Второе важное свойство интеграла столкновений состоит в выполнении Н-теоремы Больцмана, основное содержание которой, утверждает, что для замкнутого пространства функция
Больцмана Я = |/ ,п/С, _ряет неравенству ^ , 0, что С(_ву(, возрае-
танию энтропии замкнутой системы. Условие равенства интеграла столкновений нулю достигается в равновесии, когда функция распределения становится Максвеллианом. Кроме этого, интеграл столкновений на этапе релаксации сохраняет положительность функции распределения.
Основные черты кинетических решений могут быть получены решением более простых модельных уравнений, общий вид которых записывается в виде
§ + ^ (2)
где Е+ - аппроксимация интеграла обратных столкновений, деленного на частоту столкновений О(%) / у(£) , а 1/т - аппроксимация частоты столкновений у(£).
Одной из широко используемых моделей является модель BGK, в которой функция
¥+ равна Максвеллиану, М (п, иТ) = п(———)3/2 ехр(-—с2 / 2квТ), а частота столкновений
2жквТ
не зависит от скорости частиц и выражается через коэффициент вязкости л и давление p , 1 / т = р / и . Данная модель удовлетворяет выше приведенным условиям для интеграла столкновений и является аппроксимирующей моделью первого приближения для псевдо-максвелловских молекул. Одним из существенных недостатков данного приближения является неправильная величина числа Прандтля Рг при переходе к гидродинамическому режиму. Тем не менее, эта очень важная модель, дающая правильное качественное поведение кинетического решения и на ее основе был создан эффективный метод для расчета течений многокомпонентного и химически реагирующего газа [19-20]. Для получения правильного числа Прандтля в [11] была предложена эллипсоидальная статистическая модель ES, в которой функция ¥ + и частота столкновений 1/ т равны соответственно
¥+ = . п ехр(-1 сТ-1с) , 1 = ■ Р
^е<2яТ) 2 ' т (1 -ЛЕ8)/'
Здесь тензор Т выражается через температуру Т и тензор напряжения Р согласно формуле,
Т = (1 - ЛЕ5 ) КТ1 + ЛЕ5 Р / р,
где р - плотность, Я - газовая постоянная, а параметр ЛЕ8 может принимать значения -0.5 < ЛЕ8 < 1. ES-модель, также как и BGK-модель, обладает основными свойствами интеграла столкновений, сохраняя положительность функции распределения.
В [21] предложена аппроксимация с помощью разложения по ортого-
нальной системе полиномов Эрмита до третьего порядка, которая имеет вид
¥+ = М (п, и, Т)(1 + (1 - Рт)-Рс с + —(1 - Рг —)а с (—— 5)).
и 2КТр 1] 5КТр ¡л 2КТ 2
Если время релаксации т = ¡и / р, то функция
¥+ = М (п, и, Т) (1 + -^(1 - Рг)чс (£- - 5) 1 . (3)
^ 5КТр 2КТ 2 J
Релаксационный член с такой функцией ¥ + известен в литературе как модель Шахова модель). Если т = ¡и/Рг р, то функция
¥ + = М(п,и, Т)(1 + Ц-^-^сс 1 . (4)
^ Рг 2КТр ' )
Выражение (4) соответствует линеаризации ES - модели (далее ESL) , которое было рассмотрено в [17] и предложено для расчетов течений, как одноатомного газа, так и смесей газов. Обе модели (как S- модель, так и ESL- модель) могут давать отрицательные значения функции распределения, так как ¥ + не всюду положительна. Это является существенным недостатком, поэтому Н-теорема для этих моделей доказана только для случая
близкого к равновесию, когда д или р^ малы и ¥ + > 0 [21,17]. Используя на этапе релаксации в уравнении (2) функцию ¥ + в виде (3) и (4) и интегрируя с весовыми функциями
2
сс и сс , получим правильные скорости релаксации теплового потока и компонент тензора вязких напряжений для газа с псевдомаксвелловскими молекулами
Жр 1} р Жц 2 р —у = - р , =— Ч . Жг л 4 Жг 3л
Несмотря на то, что порядки аппроксимаций S-модели и ESL - модели одинаковы, получаемые решения могут иметь отличия. В [17] указывается, что ESL-модель более правильно описывает "хвосты" функции распределения и потому должна давать более точные значения параметров течения. Однако это утверждение требует численной проверки. Проведем сравнение решений указанных модельных уравнений с решением полного уравнения Больцмана.
2. Описание численного метода
В данной работе для численного решения кинетического уравнения используется метод дискретных ординат. В этом случае вычисления осуществляются в области скоростного пространства ^, ограниченного сферой, центр и радиус которой определяются по характерным параметрам задачи. В данной области вводится Декартова сетка скоростей . Функция распределения задается в центрах ячеек скоростной сетки, а кинетическое
уравнение, записанное для каждого скоростного узла, образует систему гиперболических уравнений в физическом пространстве с нелинейным источником [22]:
/ дг
Узлы сетки выбираются так, чтобы описать функцию распределения. Число скоростных узлов Ыа в трехмерном скоростном пространстве обычно достаточно большое, от нескольких тысяч до сотни тысяч. Система уравнений с источником в виде полного интеграла столкновений численно решается по явной схеме методом конечного объема, при этом значения на гранях ячеек для второго порядка точности вычисляются с использованием монотонной реконструкции. Для модельных уравнений возможно применение неявных схем [23-24], но в данной работе такой подход не используется.
В силу высокой размерности интеграла столкновений при его вычислении используют, как правило, метод Монте-Карло. Так как скорость сходимости обычного метода
Монте-Карло достаточно низкая 1/^ (N -число испытаний), то применяются процедуры понижения дисперсии, позволяющие улучшить точность вычисления, не увеличивая числа испытаний. Известными методами понижения дисперсии являются выбор случайных узлов по значимой выборки [25], выделение главной части, т.е. разбиение функции распределения на равновесную функцию и поправку к ней [4,7], а также задание квазиравномерных сеток [3,5]. Выбор скоростей столкновений в последнем случае осуществля-
- + Ч ■{1а/а) =1 а, а = Ма.
ется по последовательностям Коробова[26]. Этот метод, хотя и требует достаточно большого количества испытаний, оказывается достаточно эффективным для расчетов на равномерной сетке в скоростном пространстве. Так как последовательность случайных узлов не зависит от вида функции распределения, то расчет столкновений (наиболее трудоемкая часть) проводится только один раз для всех ячеек в физическом пространстве. Однако если функция распределения имеет сложную структуру, локализованную в узкой области скоростного пространства, или расчетная область очень протяженная, например, как при гиперзвуковых течениях, метод квазиравномерных сеток теряет свою эффективность.
В наших расчетах при умеренных числах Маха используется выборка случайных узлов по последовательности Коробова, а при больших числах Маха (М>8) скорости сталкивающихся частиц находятся по методу значимой выборки. При использовании данного
подхода предварительно строится кумулятивное распределение Е(Лп), равномерно распределенное на интервале (0,1) и имеющее вид:
п
Е(Л) = ТР/(ЛУа ,
а=0
где Vа - объем скоростной ячейки с центром Ла, а Р[ - распределение, которое в каждой ячейке физического пространства близко к интегрируемой функции /(Л),
^ (Л)=
Т/ЛаУа
а=1
Случайные узлы Л, I = 1,■■■, с заданным распределением р/определяются из решения уравнения, Е(Л1) = Г, где Г последовательность равномерно распределенных случайных чисел. То есть, для каждого г1 находится интервал скоростных узлов [Лп ,Лп+1\ такой, чтобы Е(Лп—1) < Е(Л) < Е(Лп) ■ Искомый узел случайной последовательности полагается равным Лп ■
Аналогично работам [3-7] для построения кубатурной формулы используется слабая форма интеграла Больцмана в восьмимерном пространстве О х О х 2ж х Ът
1 Ът 2Л
1(Л)=1 ¡ ¡ ¡ ¡[¿(Л—л—л)—¿(Л—Л)—¿(Л—ЛШЛ^ъсъсЕСЛСЛ,
2 0 0 О О
где КЛ,Л) = /(Л)/(Л;) вклад от прямых столкновений, ¿(Л —Л) - дельта-функция Дирака, Ь- прицельное расстояние и Ът - максимальное значение Ь, зависящее от закона взаимодействия частиц. Используя метод Монте-Карло с выбором узлов по распределению р/ и с равномерным выбором параметров столкновений е, ъ, для аппроксимации интеграла получается выражение
77-/) Ыс /■ ^
1 (5*)I— (д.+д -д* - д , ^ 1=1 -11-
где Д - аппроксимация дельта-функции на дискретной скоростной сетке
\И V для || ||ш< й./2,
д.-.) = ] 7 Ь/"" .
[0 иначе.
Здесь || .||да=max{| .х |,| .у |,| . |}, И. -шаг скоростной сетки.
Для выполнения законов сохранения на каждом столкновении проводится двухточечное проектирование скоростей обратных столкновений (5,5) на узлы скоростной сетки [3].
При расчете модельных уравнений для выполнения численно законов сохранения и получения правильной скорости релаксации неравновесных моментов (теплового потока и компонент тензора неравновесных напряжений) проводится коррекция параметров функции Г+, так чтобы
^ (Г +- / )Уа= 0,
и для Б- модели
а для ЕБЬ- модели
т- Е (г - /)с2сауа=- ч ,
2 „ 3 Л
тЕ (Га+ - /а) С шС]аУа =- -Р у . а " Л '
Решение этих нелинейных систем находится методом Ньютона.
3. Численные результаты сравнения решений модельных уравнений
Приведем для некоторых задач газовой динамики численные решения, полученные с использованием модельных уравнений и уравнения Больцмана. Рассмотренные нами ниже одномерные тесты, несмотря на свою простоту, являются элементами сложных газодинамических течений. Анализ и сравнение результатов в этих задачах позволяет увидеть реальную точность используемых моделей.
Расчеты для различных форм интеграла столкновений проводятся на одинаковых равномерных скоростных сетках и одинаковых равномерных сетках в физическом пространстве. Во всех ниже приведенных результатах используются безразмерные величины, которые нормализуются на следующие характерные параметры задачи, температуру Тю,
скорость ых = у!квТх / т , плотность числа частиц пда, длину Ь, время ^ = Ь / мю. Число Кнудсена, определяемое как Лт / Ь, где Аю - длина свободного пробега невозмущенного потока, в расчетах полагается равным 1. Принято, что взаимодействие молекул происхо-
а
дит по модели твердых сфер. Тогда безразмерная частота столкновений для Б-модели
В качестве первого примера рассмотрим одномерную нестационарную задачу разлета газа. Начальные условия задаются следующим образом, слева параметры течения газа п = 1, и = —2, Т = 1, а справа п = 1, и = 2, Т = 1.
Сильное разрежение, образующееся в центре области, приводит к увеличению локальной длины свободного пробега X << 1/п, что соответствует кинетическому режиму. Однако из-за малых значений плотности и, соответственно, частоты столкновений решения слабо чувствительны к форме интеграла и практически совпадают для модельных уравнений (Б, БББ, БОК) и для уравнения Больцмана (рис.1). Небольшое отличие в решениях имеют лишь тепловые потоки, при этом максимальное значение отклонения не превышает 4%. Таким образом, на волнах разрежения все модельные уравнения хорошо аппроксимируют полное уравнение Больцмана.
Во втором примере рассмотрим задачу о контактном разрыве. Начальные условия, слева п = 9, и = 0, Т = 1, справа п = 1, и = 0, Т = 9 ■
Расчеты показывают, что плотность числа частиц и температура, полученные из решения модельных и полного кинетического уравнения, практически совпадают, а отличия значений скорости и теплового потока не превышают 3% и 5% соответственно. Наиболее чувствительным макропараметром является ри. Максимальные отличия от значений, получаемых из уравнения Больцмана, видны в области начального разрыва (рис.2) и достигают 30%. В областях распространения волн ошибки значений рхх уменьшаются, при
этом ББЬ-модель дает более точные значения экстремумов.
В третьем расчете представлены результаты задачи о распаде разрыва, образованного встречными потоками, движущимися с одинаковыми абсолютными значениями скоростей. Начальные условия слева п = 1, и = 4, Т = 1, справа п = 1, и = —4, Т = 1 ■
Решением данной задачи являются две ударные волны, движущиеся в противоположные стороны от центра и ограничивающие область покоящегося газа. В начальные моменты времени в центре области формируется неравновесная функция распределения. На рис.3 представлены профили компоненты тензора неравновесных напряжений и теплового потока при 1=0^75. Экстремальные значения теплового потока модельные уравнения воспроизводят правильно, а компоненты тензора неравновесных напряжений с ошибкой порядка 20%. При этом решения модельных уравнений близки между собой. На рис.4 приведены профили плотности числа частиц и температуры в центре области (при х=15) в зависимости от времена Видно, что решение по БББ-модели ближе к полному кинетическому решению, чем решения по другим модельным уравнениям.
При больших временах экстремальные значения теплового потока и компоненты тензора неравновесных напряжений, посчитанные по модельным уравнениям (Б , БББ) и по уравнению Больцмана, практически совпадают (рис. 5). Основные отличия от полного кинетического решения появляются перед движущимися фронтами ударных волн. Ошиб-
ка по тепловому потоку для S- модели (ESL- модели), отнесенная к максимальному значению величины, составляет 20% (25 %), а для компоненты тензора неравновесных напряжений 25% (35%). Из графиков, представленных при t= 1.6, видны искажения профилей теплового потока и компоненты тензора неравновесных напряжений, особенно для ESL-модели.
В следующем (четвертом) примере рассмотрим решение стационарной задачи о структуре ударной волны. На рис.(6-7) показаны результаты расчетов ударной волны с числом Маха =10. Также как и в предыдущем примере, видно, что модельные уравнения дают правильные значения экстремумов q и pxx, а основные отличия наблюдаются перед фронтом ударной волны. При этом отклонение решения, получаемое по ESL- модели, значительно больше.
Изменение поведения неравновесных макропараметров и паразитный нагрев газа перед фронтом волны объясняется наличием потока быстрых молекул из дозвуковой области в область сверхзвука. Если частота не зависит от скорости, то для быстрых частиц длина свободного пробега Л = / v ^ « и релаксация функции происходит медленно, в отличие от решения уравнения Больцмана, в котором при взаимодействии молекул по закону твердых сфер, <^х / v ^ const. Для модели ESL затухание функции еще более медленное, так как частота столкновений содержит дополнительный множитель Pr .
На рис.8(а) показана функция f (£x, Щ /2, Щ / 2) перед фронтом ударной волны при
x=5. На рис.8(б) приведена та же функция для < 0 . Видно, что "всплеск" функции, полученный решением по ESL-модель в области отрицательных скоростей больше, чем по другим модельным уравнениям (BGK, S). Так как значения функции распределения в области отрицательных скоростей малы, то полученные отклонения функции распределения мало влияют на плотность газа, однако ошибка в высоких моментах оказывается значительной.
Из приведенных численных примеров видно, что в областях сильной неравновесности газа отклонения решений модельных уравнений (S , ESL) одного порядка и проявляются в основном в поведении продольной компоненты тензора неравновесных напряжений, а остальные макропараметры, полученные из решения модельных уравнений и уравнения Больцмана, достаточно близки. Однако значительную погрешность модельные уравнения дают перед фронтами ударных волн, где отсутствие зависимости частоты столкновений от скорости, приводит к избыточным потокам горячих молекул в область холодного газа. При этом ESL-модель дает больше отклонение из-за меньшей частоты столкновений. Положение температурного фронта ударной волны также более точное для модели Шахова. Этот факт дает основания полагать, что в пространственных задачах обтекания тел положение ударной волны и значения макропараметров за фронтом будет передаваться S-моделью корректно. Однако в пристеночных областях (аналог задачи о контактном разрыве) использование ESL- модели может дать меньшую ошибку величины теплового потока.
0.0 Рчя
■Я.05
-0.1 а
а
•
V
ю
Рис. 1. Разлет потоков при t=0.45. Продольная компонента тензора неравновесных напряжений (а),
плотность числа частиц (б)
0.8 Рхх 0.4
0.0
-0.4
-0.8
(а) Л
т 3 ъ.
ВоИггпапп - 1 1
3-ллос1е1 • \ N
Е31_-тск1е1 Г
10
15
Рис. 2. Контактный разрыв. Продольная компонента тензора неравновесных напряжений (а), тепловой
поток (б) при 1= 0.35.
Рис. 3. Встречные потоки при t=0.75. Продольная компонента тензора неравновесных напряжений (а),
тепловой поток (б)
Рис. 4. Встречные потоки. Профили плотности числа частиц (а) и температуры (б) в зависимости от времени
в центре области ( х=15).
Рис. 5. Встречные потоки. Сравнение решений модельных уравнений и уравнения Больцмана. Продольная компонента тензора неравновесных напряжений (а), тепловой поток (б) при t=1.6.
Рис. 6. Структура ударной волны М=10. Профили плотности (а), температуры (б).
Рис. 7. Структура ударной волны М=10. Продольная компонента тензора неравновесных напряжений (а),
тепловой поток (б).
Рис. 8. Функция распределения перед фронтом ударной волны (а), функции распределения при отрицательных значениях молекулярной скорости, (крупный план) (б).
Заключение
В данной работе было проведено сравнение решений модельного уравнения Шахова и линеаризованной ББ - модели с решением полного уравнения Больцмана для четырех тестовых задач газовой динамики. Численно показано, что в задаче о контактном разрыве ББЬ-модель более точно передает максимальное значение теплового потока, а также экстремальные значения продольной компоненты тензора неравновесных напряжений для распространяющихся от центра волн. В остальных примерах модельные уравнения дают небольшие отклонения макропараметров (плотности числа частиц, температуры и теплового потока) в областях сильной неравновесности функции распределения, а ошибки значений компоненты тензора неравновесных напряжений для ББЬ- модели и Б- модели одного порядка. Однако перед фронтами ударных волн решения модельных уравнений и уравнения Больцмана отличаются значительно, при этом ББЬ- модель дает большее отклонение макропараметров, чем модель Шахова.
Работа выполнена при финансовой поддержке РФФИ, грант № 15-07-02986А и программы фундаментальных научных исследований, тема №14 по плану НИР ВЦ РАН ФИЦ ИУ РАН , № 01201352394.
Список литературы
1. Bird G. A. Molecular Gas Dynamics and the Direct Simulation of Gas Flow. Oxford: Clarendon Press, 1994. 458 p.
2. Yen S. M. Numerical solution of the nonlinear Boltzmann equation for nonequillibrium gas flow problems // Annual Review of Fluid Mechanics. 1984. V.16. Pp. 67-97.
3. Черемисин Ф.Г. Консервативный метод вычисления интеграла столкновений Больцмана // Доклады Академии Наук. 1997. Т. 357, № 1. С. 53-56.
4. Morris A. B.,Varghese P. L.,Goldstein D. B. Monte Carlo solution of the Boltzmann equation via a discrete velocity model //Journal of Computational Physics. 2011. Vol. 230. Pp 1264-1280. DOI: 10.1016/j.jcp.2010.10.037
5. Черемисин Ф. Г. Решение кинетического уравнения Больцмана для высокоскоростных течений// Ж. вычисл. матем. и матем. физ. 2006. T. 46, № 2. C. 329-343.
6. Arslanbekov R. R., Kolobov V. I., Frolova A. A. Kinetic Solvers with Adaptive Mesh in Phase Space // Physical Review E. 2013. V. 88. 063301.
7. Radtke G. A., Hadjiconstantinou N. G. Variance-reduced particle simulation of Boltzmann transport equation in the relaxation-time approximation // Physical Review E. 2009. V.79. P. 056711. DOI: 10.1103/Phys. Rev. E.79. 056711
8. Radtke G. A., Hadjiconstantinou N. G., Wagner W. Low-noise Monte Carlo simulation of the variable hard sphere gas // Physics of Fluids. 2011. V. 23. P. 030606. DOI: 10.1063/1.3558887
9. Иванов М. С., Коротченко М. А., Михайлов Г. А., Рогазинский С. В. Глобально-весовой метод Монте-Карло для нелинейного уравнения Больцмана // Ж. вычисл. матем. и матем. физ. 2005. T. 45, №10. С. 1860-1870.
10. Bhatnagar P. L. Gross E. P., Krook M. A model for collision process in gases // Physical Review. 1954. V. 94. Pp. 511-525.
11. Holway L. H. New statistical models for kinetic theory: Methods of construction// Physics of Fluids. 1966. V. 9. Pp. 1658-1673.
12. Шахов Е. М. Об обобщении релаксационного кинетического уравнения Крука // Изв. АН. СССР. МЖГ. 1968. №5. С.142-145.
13. Struchtrup H. The BGK model with velocity dependent collision frequency // Continuum Mechanics and Thermodynamics. 1997. V. 9, №1. Pp. 23-31.
14. Zheng Y. , Struchtrup H. Ellipsoidal statistical Bhatnagar-Gross-Krook model with velocity-dependent collision frequency // Physics of Fluids. 2005. V. 17. P. 127103. DOI: 10.1063/1.2140710
15. Титарев В. А. Шахов Е.М. Численный расчет поперечного обтекания холодной пластины гиперзвуковым потоком разреженного газа // Механика жидкости и газа. 2005. №5. С. 152-167.
16. Andries P., Perthame B. The ES-BGK model equation with correct Prandtl numbed/Rarefied Gas Dynamics: 22nd International Symposium: AIP Conf. Proc. 2001, CP 585. Pp. 30-36.
17. Belyi V.V. Derivation of model kinetic equation // EPL. 2015. V. 111. 40011. (DOI: 10.1209/0295-5075/111/4011 )
18. Коган М.Н. Динамика разреженного газа. М.: Наука, 1967. 440 с.
19. Andries P., Aoki K. , Perthame B. A consistent BGK-type model for gas mixture // J. Stat. Phys. 2002. V.106, N. 516. Pp. 993-1113
20. Groppi M., Spiga G. A Bhatnagar-Gross-Krook -type approach for chemically reacting gas mixture // Physics of Fluids. 2004. V. 16, № 12. Pp 4273-4284.
21. Шахов Е.М. Метод исследования движений разреженного газа. М.: Наука, 1974. 203 с.
22. 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 refinrment // Journal of Computational Physics. 2007. V. 223. Pp. 589-608.
23. Титарев В. А. Неявный численный метод расчета пространственных течений разреженного газа на неструктурированных сетках. // Ж. вычисл.матем. и матем. физ. 2010. Т. 50, № 10. С.1811-1826.
24. Титарев В. А. Программный комплекс Несветай-ЗД моделирования пространственных течений одноатомного разреженного газа //Наука и образование. МГТУ им. Н.Э. Баумана. Элект. Журнал. 2014, N. 6. C. 124-154.
25. Tan Z., Varghese P.L. The ^—s method for the Boltzmann equation // Journal of Computational Physics. 1994. V.110. Pp. 327-340.
26. Коробов Н.М. Тригонометрические суммы и их приложение. М.: Наука, 1989. 240 с.
Mathematics I Mathematical Modelling
Mathematics and Mathematical Madelling of the Bauman MSTU, 2015, no. 06, pp. 61-77.
DOI: 10.7463/mathm.0615.0823537
Received: Revised:
08.11.2015 22.11.2015
© Bauman Moscow State Technical Unversity
Numerical Comparison of Solutions of Kinetic Model Equations
A.A. Frolova1*
aafrolovaiffv andexju
institution of Russian Academy of Sciences Dorodnicyn Computing Centre of RAS, Moscow, Russia
Keywords: Boltzmann equation, model equations, the equation of BGK, ellipsoidal statistical model
(ES), Shakhov model (S-model)
The collision integral approximation by different model equations has created a whole new trend in the theory of rarefied gas. One widely used model is the Shakhov model (S-model) obtained by expansion of inverse collisions integral in a series of Hermite polynomials up to the third order. Using the same expansion with another value of free parameters leads to a linearized ellipsoidal statistical model (ESL).
Both model equations (S and ESL) have the same properties, as they give the correct relaxation of non-equilibrium stress tensor components and heat flux vector, the correct Prandtl number at the transition to the hydrodynamic regime and do not guarantee the positivity of the distribution function.
The article presents numerical comparison of solutions of Shakhov equation, ESL- model and full Boltzmann equation in the four Riemann problems for molecules of hard spheres.
We have considered the expansion of two gas flows, contact discontinuity, the problem of the gas counter-flows and the problem of the shock wave structure. For the numerical solution of the kinetic equations the method of discrete ordinates is used.
The comparison shows that solution has a weak sensitivity to the form of collision operator in the problem of expansions of two gas flows and results obtained by the model and the kinetic Boltzmann equations coincide.
In the problem of the contact discontinuity the solution of model equations differs from full kinetic solutions at the point of the initial discontinuity. The non-equilibrium stress tensor has the maximum errors, the error of the heat flux is much smaller, and the ESL - model gives the exact value of the extremum of heat flux.
In the problems of gas counter-flows and shock wave structure the model equations give significant distortion profiles of heat flux and non-equilibrium stress tensor components in front of the shock waves. This behavior is due to fact that in the models under consideration there is no dependency of the collision frequency on the molecular velocity.
As calculations show, the ESL-model describes more accurately the non-equilibrium flow regime, but gives a greater deviation from the Boltzmann equation, than the Shahov model in front of shock waves.
References
1. Bird G.A. Molecular gas dynamics and the direct simulation of gas flow. Oxford: Clarendon Press, 1994. 458 p.
2. Yen S.M. Numerical solution of the nonlinear Boltzmann equation for nonequillibrium gas flow problems. Annual Review of Fluid Mechanics, 1984 vol.16, pp.67-97.
3. Cheremisin F.G. Conservative method of calculating the Boltzmann collision integral. Doklady Akademii Nauk. 1997, vol.357, no.1, pp.53-56.(in Russian)
4. Morris A.B., Varghese P.L., Goldstein D.B. Monte Carlo solution of the Boltzmann equation via a discrete velocity model. Journal of Computational Physics, 2011, vol.230, pp.1264-1280. D0I:10.1016/j.jcp.2010.10.037
5. Tcheremissine F.G. Solution to the Boltzmann kinetic equation for high-speed flows. Computational mathematics and mathematical physics. 2006, vol.46, no.2, pp.315-329. (in Russian). DOI: 10.1134/S0965542506020138
6. Arslanbekov R.R., Kolobov V.I., Frolova A.A. Kinetic solvers with adaptive mesh in phase space. Physical Review E, 2013, vol.88, pp.063301.
7. Radtke G.A., Hadjiconstantinou N.G. Variance-reduced particle simulation of Boltzmann transport equation in the relaxation-time approximation. Physical Review E. 2009, vol.79, pp.056711. DOI: 10.1103/PhysRevE.79.056711
8. Radtke G.A., Hadjiconstantinou N.G., Wagner W. Low-noise Monte Carlo simulation of the variable hard sphere gas. Physics of Fluids, 2011, vol.23, pp.030606. DOI: 10.1063/1.3558887
9. Ivanov M.S., Korotchenko M.A., Mikhailov G.A., Rogasinsky S.V. Global weighted Monte Carlo method for nonlinear Boltzmann equation. Zhurnal vychislitel'noi matematiki i matematicheskoifiziki, 2005, vol.45, no.10, pp.1860-1870. (English version of journal: Computational Mathematics and Mathematical Physics, 2005, vol.45, no.10, pp.1792-1801.)
10. Bhatnagar P.L. Gross E.P., Krook M.A model for collision process in gases. Physical Review, 1954, vol.94, pp.511-525.
11. Holway L.H. New statistical models for kinetic theory: Methods of construction. Physics of Fluids. 1966, vol.9, pp.1658-1673. DOI: 10.1063/1.1761920
12. Shahov Е.М. Generalizatiron of the Krook relaxation kinetic equation. Izvestuya Akademii Nauk SSSR Mekhanika Zhidkosti i Gaza, 1968, no.5, pp.142-145. (in Russian). (English version of journal: Fluid Dynamics, 1968, vol.3, no.5, pp.95-96. DOI: 10.1007/BF01029546 )
13. Struchtrup H. The BGK model with velocity dependent collision frequency. Continuum Mechanics and Thermodynamics, 1997, vol.9, no.1, pp.23-31. DOI: 10.1007/s001610050053
14. Zheng Y., Struchtrup H. Ellipsoidal statistical Bhatnagar-Gross-Krook model with velocity-dependent collision frequency. Physics of Fluids, 2005, vol.17, pp.127103. DOI: 10.1063/1.2140710
15. Titarev V.A. Shakhov E.M. Numerical calculations of hypersonic rarefied gas transverse flow past a cold flat plate. Izvestuya Akademii Nauk. Mekhanika Zhidkosti i Gaza, 2005, no.5, pp.152-167. (in Russian). (English version of journal: Fluid Dynamics, 2005, vol.40, no.5, pp 790-804. DOI: 10.1007/s10697-005-0117-1 )
16. Andries P., Perthame B. The ES-BGK model equation with correct Prandtl number. Rarefied Gas Dynamics: 22nd International Symposium: AIP Conf. Proc., 2001, vol.585, pp.30-36. DOI: 10.1063/1.1407539
17. Belyi V.V. Derivation of model kinetic equation. Europhysics Letters, 2015, vol.111, pp.40011. DOI: 10.1209/0295-5075/111/4011
18. Kogan M.N. Dinamika razrezhennogo gaza [Rarefied gas dynamics]. Moscow, Nauka Publ., 1967. 440 p. (in Russian)
19. Andries P., Aoki K., Perthame B. A consistent BGK-type model for gas mixture. Journal of Statistical Physics, 2002, vol.106, no.5, pp.993-1118. DOI: 10.1023/A:1014033703134
20. Groppi M., Spiga G. A Bhatnagar-Gross-Krook-type approach for chemically reacting gas mixture. Physics of Fluids, 2004, vol.16, no.12, pp.4273-4284. DOI: 10.1063/1.1808651
21. Shakhov E.M. Metod issledovaniya dvizheniy razrezhennogo gaza [The method for analyzing rarefied gas flows]. Moscow, Nauka Publ., 1974, 203 p. (in Russian)
22. 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 Compu-tationalPhysics, 2007, vol.223, no.2, pp.589-608. DOI: 10.1016/j.jcp.2006.09.021
23. Titarev V.A. Implicit numerical method for computing three-dimensional rarefied gas flows using unstructured meshes. Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki, 2010, vol.50, no.10, pp.1811-1826. (in Russian). (English version of journal: Computational Mathematics and Mathematical Physics, 2010, vol.50, no.10, pp.1719-1733. DOI: 10.1134/S0965542510100088 )
24. Titarev V.A. Software package nesvetay-3D for modeling three-dimensional flows of mona-tomic rarefied gas. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2014, no.6, pp.124-154. (in Russian). DOI: 10.7463/0614.0712314
25. Tan Z., Varghese P.L. The A —S method for the Boltzmann equation. Journal of Computational Physics, 1994, vol.110, pp.327-340. DOI: 10.1006/jcph.1994.1030
26. Korobov N.M. Trigonometricheskie summy i ikh prilozhenie [Trogonometric sums and their applications]. Moscow, Nauka Publ., 1989. 240 p. (in Russian).
