CC BY
7_ДОКЛАДЫ АН ВШ РФ_
2016_апрель-июнь_№ 2 (31)
- ТЕХНИЧЕСКИЕ НАУКИ -
УДК 531.19:533.72
ПРИМЕНЕНИЕ МЕТОДА ПРОБНЫХ ЧАСТИЦ ДЛЯ ВЫЧИСЛЕНИЯ КОЭФФИЦИЕНТОВ ПЕРЕНОСА В РАЗРЕЖЕННЫХ ГАЗАХ
К.А. Финников
Сибирский федеральный университет
Традиционные методы расчета коэффициентов переноса в разреженных газах основаны на разложении первого приближения функции распределения частиц в ряд по какому-либо базису функций. В ряде случаев могут оказаться полезными альтернативные методы, основанные на статистическом моделировании движения частиц. В качестве подобного метода в данной статье предлагается вариант метода пробных частиц, модифицированный таким образом, чтобы плотность пробных частиц аппроксимировала первое приближение к функции распределения в методе Чепмена-Энскога. Основными особенностями метода являются следующие. Во-первых, это введение дополнительной характеристики пробной частицы - статистического веса, учитываемого в виде множителя при определении вклада пробных частиц в функцию распределения. Во-вторых, это специальный алгоритм моделирования столкновений пробных частиц с фоновыми, отражающий структуру линеаризованного столкновительного члена в уравнении для первого приближения функции распределения. В данном алгоритме столкновение приводит к исчезновению пробной частицы, вступившей в столкновение, и рождению трех новых, в результате чего вместо траектории единичной пробной частицы рассчитывается тернарное дерево траекторий. Траектории частиц зависимы друг от друга в пределах одного дерева, но различные деревья траекторий могут рассчитываться независимо. Это обеспечивает неограниченную возможность масштабирования алгоритма в многопроцессорных системах, характерную для других методов пробных частиц. Рассмотрен случай однокомпонентного газа с центральным потенциалом межчастичного взаимодействия (потенциалы твердых сфер и Леннард-Джонса). В случае газа твердых сфер результаты, полученные с использованием разработанного метода, идентичны известным эталонным аналитическим результатам. Расчеты коэффициентов переноса в инертном газе с использованием потенциала Леннард-Джонса сопоставлялись с данными эксперимента. Расхождение составляет около 1 %, что соответствует точности приближения потенциала межатомного взаимодействия формулой Леннард-Джонса. Разработанный метод может использоваться как параллельно с традиционным, в целях контроля вычислительных погрешностей и отсутствия ошибок в аналитических вычислениях, так и самостоятельно.
Ключевые слова: коэффициенты переноса, разреженный газ, приближение Чепмена-Энскога, метод пробных частиц, вязкость, теплопроводность, статистическое моделирование.
Б01: 10.17212/1727-2769-2016-2-108-122
Введение
Определение коэффициентов переноса разреженных газов на основании данных о потенциалах межчастичных взаимодействий является одним из основных приложений кинетической теории. Соотношения, описывающие явления переноса, появляются в результате применения приближения малых чисел Кнудсена к уравнению Больцмана, в первом порядке такого приближения [1]. Основным результатом здесь является уравнение для первого приближения функции рас-
Работа выполнена при частичном финансировании РНФ (соглашение № 14-19-00312).
© 2016 К.А. Финников
пределения /1 или, в случае многокомпонентного газа, система таких уравнений, в которой фигурируют /1 всех компонент. В случае однокомпонентного
газа уравнение для /1 содержит известные слагаемые, пропорциональные пространственным и временным производным газодинамических параметров, которые в совокупности мы обозначим как Я(у), и линейный оператор, действующий на /1 -столкновительный член:
Я (у) + Я(/1) = 0; Я(/1) = = п\[ /0 (у') /1 (V') + /0 (V') /1 (у') - /0 (у) /1 (V) - /0 (V) /1 (у) ] х
х|у - VI(! у - V |)(1)
Нахождение /1 и ее моментов - потоков частиц, импульса и энергии и в дальнейшем коэффициентов переноса является задачей скорее технической, но в то же время достаточно сложной и вычислительно затратной. Наиболее употребителен способ решения этой задачи путем разложения /1 в ряд с использованием полиномов Сонина. В результате такого разложения уравнение (1) сводится к системе линейных уравнений и, в итоге, выводятся выражения для коэффициентов переноса, включающие интегралы столкновений а1 (Т). Последние определяются законами парного взаимодействия частиц. Расчет интегралов столкновений может быть достаточно трудоемким, особенно в случае нецентрального характера взаимодействия частиц, и зачастую требует привлечения методов Монте-Карло. Получаемые таким способом значения коэффициентов переноса имеют как статистическую погрешность, так и погрешность, связанную с обрывом разложения /1 в ряд. Если первая может быть снижена путем простого увеличения объема статистики, то уменьшение второй требует существенного усложнения расчетных выражений и увеличения числа необходимых интегралов столкновений. В некоторых случаях разложение по полиномам Со-нина сходится достаточно медленно. Так, расчет коэффициентов переноса электронной компоненты в плазме может потребовать не менее 5-10 слагаемых ряда [2]. Аналогичной особенностью обладают все смеси сильно различающихся по молекулярной массе компонент, в частности гелий-ксеноновая смесь, имеющая важное практическое значение [3]. Медленная сходимость разложения связана со сложной зависимостью времени релаксации скорости более легкой молекулы от ее энергии. Вычисления коэффициентов переноса существенно усложняются в случае многокомпонентного газа, при рассмотрении газа частиц, имеющих внутренние степени свободы.
Таким образом, представляет интерес разработка альтернативных методов расчета коэффициентов переноса в разреженных газах, основанных на прямом моделировании движения частиц. Примером такого подхода является работа [4], в которой развивается расчетный метод, основанный на моделировании равновесного состояния газа и применении формул Грина-Кубо. Особенностью данного метода является то, что в нем необходимо моделировать коллектив достаточно большого (порядка нескольких тысяч) числа молекул. В настоящее время метод опробован только для случая газа твердых сфер.
Методы пробных частиц являются мощным инструментом, позволяющим напрямую решать уравнение Больцмана путем расчета множества траекторий частиц. Этим методам присуща простота алгоритма, сохраняющаяся при необходимости учитывать большое количество столкновительных процессов, возможность применения параллельных расчетов на неограниченном числе процессоров, возможность получения решения уравнения Больцмана, погрешность которого обусловлена практически только конечностью статистической выборки [5]. Отмеченные особенности относятся к тем случаям использования методов пробных частиц, в которых искомой является функция распределения частиц, сталкивающихся только с частицами других сортов, неподвижными или имеющими известную функцию распределения; последние называются фоновыми или полевыми. В таких случаях уравнение Больцмана линейно, из чего вытекает возможность независимого расчета траекторий частиц.
Уравнение (1), имеющее структуру, отличную от структуры уравнения Больцмана для пробных частиц, в то же время схоже с последним в том, что оно линейно и по своей природе является уравнением переноса частиц в пространстве скоростей. Это дает основания полагать, что некоторая модификация метода пробных частиц может быть использована для получения решения (1). Подобный метод может использоваться для расчета коэффициентов переноса как совместно с традиционными методами, с целью взаимного контроля расчетных результатов, так и независимо от них, в тех случаях, когда он будет более удобным и вычислительно экономичным.
Целью настоящей работы является адаптация метода пробных частиц к задаче расчета коэффициентов переноса. Будет рассмотрен случай однокомпонент-ного газа и обсуждена возможность обобщения метода на случай многокомпонентных смесей.
1. Адаптация метода пробных частиц
Общей чертой методов пробных частиц является то, что решаемые с помощью них уравнения (уравнение Больцмана, родственные ему уравнения переноса нейтронов и радиационного переноса тепла) выражают сохранение числа частиц и их перенос в фазовом пространстве. Уравнение (1) можно интерпретировать аналогичным образом, но применительно не к полной концентрации частиц в фазовом пространстве, а к ее возмущению. Тогда Я( у) в (1) получает смысл источника, обусловленного наличием возмущающих факторов - неоднородности газодинамических параметров. Интеграл столкновений выражает перераспределение возмущения концентрации по пространству скоростей. В результате перераспределения возмущения последнее затухает, поскольку источник всегда дает нуль при интегрировании по направлениям. Процесс генерирования и затухания возмущения может быть смоделирован процессом рождения и гибели пробных частиц. При этом алгоритм моделирования рождения и гибели пробных частиц нуждается в следующих модификациях.
I. Поскольку как источник возмущения, так и /1 знакопеременны, необходимо учесть возможность отрицательного вклада пробной частицы в /1. Для этого пробным частицам присваивается специальная характеристика которая в дальнейшем будет именоваться весом. Вес частицы используется в качестве дополни-
тельного множителя при определении вклада пробной частицы в возмущение функции распределения и может иметь как положительную, так и отрицательную величину.
II. Поскольку частицы, которые рассматриваются как пробные, принадлежат к тому же сорту, что и фоновые частицы, процесс столкновения между ними требует специального описания. Изменение скорости фоновой частицы в результате столкновения вносит дополнительное возмущение, которое можно описать, используя введенную выше дополнительную характеристику пробной частицы -вес. Пусть пробная частица с весом V и скоростью у сталкивается с фоновой частицей, имеющей скорость V , и в результате первая приобретает скорость у', а вторая - скорость V'. Тогда результат столкновения можно представить как совокупность следующих событий:
1) исчезновение пробной частицы;
2) возникновение трех новых пробных частиц: со скоростью у' и весом V , со скоростью V и весом -V, со скоростью V' и весом V . Такое описание процесса столкновений соответствует структуре подынтегрального выражения в столкно-вительном члене (1), содержащего 4 слагаемых. Аналогично, в обычном методе пробных частиц изменение скорости пробной частицы, идентичное исчезновению частицы с прежним значением скорости и появлению - с новым, отражает структуру столкновительного члена в уравнении Больцмана, содержащего два слагаемых.
Таким образом, особенностью метода пробных частиц в применении к рассматриваемой задаче является то, что столкновение пробной частицы с фоновой должно моделироваться как исчезновение пробной частицы, вступившей в процесс столкновения, и появление трех новых пробных частиц. Все множество частиц, порождаемых одной, может быть представлено в виде тернарного дерева, ветви которого соответствуют частицам, а точки ветвления - столкновениям. Время существования каждой пробной частицы ограничено ее свободным пробегом и ее скорость постоянна. В дальнейшем мы будем пользоваться понятиями первичная частица для первоначально инициируемых пробных частиц и дерево частиц для множеств частиц, порождаемых первичными. Будет удобным считать, что первичные частицы сами принадлежат к деревьям, порождаемым ими.
С формальной точки зрения, дерево частиц должно развиваться неограниченно. На практике приходится использовать те или иные критерии прерывания расчета дерева. В качестве условия прекращения расчета выберем истечение определенного периода времени /1гее с момента старта первичной частицы дерева. Уточним, что каждая частица дерева либо исчезает в результате столкновения, либо достигает времени /1гее. Конечность /1гее порождает погрешность результата, в добавление к обычной погрешности статистических методов, связанной с конечностью количества рассчитанных деревьев частиц.
Алгоритм расчета дерева частиц наиболее удобно построить на основе принципа рекурсии, поскольку часть дерева, порождаемая любой из его частиц, сама является тернарным деревом. Общая структура рекуррентно вызываемой процедуры выглядит следующим образом.
Процедура Расчет Дерева (параметры: скорость у, вес V , время от момента старта первичной частицы дерева /).
Начало.
1. Выбрать следующие случайные величины, взаимосвязанные и зависящие от у, плотности и температуры газа и потенциала взаимодействия частиц: т - время движения частицы до столкновения, V - скорость частицы фона, с которой произойдет столкновение, у' и V' - скорости частиц после столкновения.
2. Если т > /1гее -t, то т := /1гее - /.
3. Учесть вклад частицы в возмущение, описываемый величинами у, т, м .
4. Если т < /1гее -/, то:
4.1. Выполнить РасчетДерева (у', м, / + т);
4.2. Выполнить РасчетДерева (V', м, / + т);
4.3. Выполнить РасчетДерева (V, - м, / + т).
Конец процедуры.
Дерево частиц, порождаемое первичной пробной частицей со скоростью у0 и весом м, рассчитывается в процессе выполнения описанной процедуры, вызываемой с параметрами (у0 , м, 0). Расчеты деревьев, порожденных разными первичными частицами, могут проводиться параллельно.
Выбор времени движения частицы до столкновения и других параметров, которые должны быть определены в п. 1 вышеприведенного алгоритма, может производиться различными способами. Здесь возможно применение прямого метода определения искомого вектора случайных величин с использованием обратной функции распределения [6]. В настоящем исследовании использовался метод, более отвечающий его целям и являющийся, возможно, более вычислительно затратным, но в то же время алгоритмически простым. Свободный пробег частицы рассматривается как совокупность шагов достаточно малой временной продолжительности. Если длительность шага достаточно мала для того, чтобы вероятность столкновения даже в течение достаточно большого количества шагов А81ер была
малой, возможные столкновения на А§1ер последовательных шагов могут считаться статистически независимыми. Поэтому на каждом шаге скорость фоновых частиц может определяться в соответствии с их собственной функцией распределения (максвелловской), а вероятность столкновения - рассчитываться исходя из величины относительной скорости пробной и фоновой частицы и полного сечения столкновения, зависящего от этой скорости. Использование этого метода вносит погрешность, связанную с конечностью вероятности столкновения в течение единичного временного шага. Влияние этой погрешности будет видно из расчетных результатов. Отметим, что и использование метода, использующего обратную функцию распределения, также вносит погрешность, связанную с конечностью массива значений обратной функции.
Возмущение, создаваемое произвольным источником К( у), может быть определено в результате усреднения по деревьям, первичные частицы которых имеют распределение по скоростям и веса, соответствующие форме Я(у). Сформулируем следующее утверждение, доказательство которого выходит за рамки настоящего исследования: если первичные частицы имеют случайные скорости, подчиняющиеся функции распределения / (у) и им присваиваются веса, описываемые такой зависимостью м( у), что
/ (у)м( у) = Я( у),
то для решения уравнения для /1 (у) действует следующее слабое предельное соотношение
—УъЩ з( у - уР)-> /1(у), (2)
где N - количество первичных частиц, суммирование производится по всем частицам всех деревьев, хг-, м>{, ур - времена существования частиц, их веса и скорости. Поясним, что понятие слабого предела подразумевает равенство между любыми моментами правой и левой частей (2), достигаемое в указанном пределе; это понятие используется здесь по той причине, что левая часть (2) в отличие от
/1 (у) - неограниченная функция.
Соотношение (2) используется для вывода выражений для моментов /1 (у), в том числе - потоков импульса и энергии. Так, компонента потока импульса определяется следующим образом:
тху = тп\УхУу/1(у) ду и тпN^т^ урхуру, (3)
точное равенство достигается в пределе N ^да, /1гее ^ да .
Рассмотрим случай возмущения, создаваемого сдвигом средней скорости газа. Пусть средняя скорость газа и = (их (у) ,0,0). Тогда в локально-сопровождающей системе отсчета
ди тугу,, 0
Я( у) = -ди^-гТУ/ 0( у). (4)
ду кТ
Для такой формы источника возможным (не единственным) вариантом является следующий выбор функции распределения и весовой функции первичных частиц:
/ (у) = / (у); *(у) = - дх у. (5)
ду кТ
После расчета необходимого числа деревьев частиц, скорости и веса которых определяются в соответствии с (5), данные расчета используются для определения тху согласно (3); далее очевидным образом находится коэффициент вязкости.
Аналогичным образом рассматривается случай неоднородности температуры и рассчитывается коэффициент теплопроводности. Описанный порядок действий является простейшим методом рассмотрения явлений переноса с помощью адаптированного метода пробных частиц.
Результаты расчетов коэффициентов переноса, полученные с использованием описанного метода для газа твердых сфер, показаны на рисунке 1. Расчетные коэффициенты вязкости и теплопроводности приведены в отношении к аналитическим, полученным методом Чепмена-Энскога в приближении одного полинома Сонина [1]. В расчетах варьировалась величина /1гее с целью наблюдения динамики сходимости результатов по этой величине. Для удобства анализа результаты
приведены в зависимости от параметра, производного от /1гее - среднего числа поколений частиц в дереве, приблизительно пропорционального . Количество рассчитываемых деревьев выбиралось таким, чтобы погрешность статистического усреднения не превышала 1 %, и составляло от 10 до 20 тысяч. Величина временного шага при расчете свободного пробега частиц выбиралась такой, чтобы вероятность столкновения в течение шага не превышала в одном случае, 0,1, в другом случае - 0,2. Полученные в этих двух случаях результаты отличаются друг от друга не более чем на 0,5 %.
Как можно заметить, сходимость результатов по является достаточно медленной, особенно для коэффициента теплопроводности. При этом увеличение числа поколений частиц на единицу приводит к увеличению расчетных затрат втрое. Для 17 поколений частиц получение приведенных результатов требует около 300 часов расчетного времени на процессорном ядре с частотой 2 ГГц. В следующем разделе будет показано, какие усовершенствования возможны в описанном методе применительно к частным задачам - задачам расчета коэффициентов вязкости и теплопроводности.
1,1
0,4 -I-i-i-i-1
0 5 10 15 20
Среднее количество поколений в деревьях частиц
Рис. 1 - Отношение расчетных коэффициентов переноса в газе
твердых сфер к аналитическим Fig. 1 - Calculated transport coefficients in a hard sphere gas in relation to analytical ones
2. Итерационный алгоритм
Медленная сходимость результата при увеличении ^ee показывает, что нецелесообразно улучшать точность расчетов путем простого увеличения этого параметра. Погрешность конечного ttree связана с пренебрежением вкладом, который
вносят в f1 частицы дерева после истечения этого интервала времени от начала старта первичной частицы. Учесть их вклад можно итерационным способом. Будем рассматривать каждую частицу дерева, существующую на момент обрыва расчета t = ttree, как первичную, порождающую собственное дерево частиц, и
определять вклад этого дерева в /1 на основе уже имеющейся информации. Для максимально эффективной организации итерационного процесса необходимо учесть аксиальную симметрию рассеяния при столкновении частиц. Рассмотрим фундаментальное решение уравнения (1), получаемое для точечного источника возмущения:
5(у - Уо) + Щ/1) = 0: /(у) = К (у, у<,). Для К (у, Уо), как частного случая решения (1), действует соотношение (2), при
/ (у) = 5(у - Уо), ^(у) = 1.
По причине аксиальной симметрии процесса рассеяния К (у, уо) зависит от угла между векторами у и уо , но не от направлений каждого из них по отдельности. Для произвольной формы источника искомое возмущение определяется как
/(у) = | К (у, уо) Я( уо)^о.
Подставим сюда источник (4), а полученное выражение используем для определения компоненты потока импульса. Выполняя интегрирование по направлениям, приходим к следующему выражению для коэффициента вязкости:
ц = 110n Jmv2 J rn^fo (v0)J K(v, v0 )(cos2 6 -1/ 3)dQe dv0dv ;
о
где 6 - угол между v и vo . Вводя функцию
Кц (vo) = J mv2 —1o J K (v, vo )(cos2 6 - 1/3)d (5)
10 o
приведем выражение для коэффициента вязкости к следующему виду:
2
ц = nJ^ fo (vo )Кц (vo )dvo. (6)
Функция Кц(vo) подлежит непосредственному определению из расчетов деревьев пробных частиц. Подставляя соотношение (2) в (5), получаем
Кц(vo) * 11oN^T,W,m(vP)2 (cos2 6, - V3), (7)
где данные в правой части получены для деревьев частиц, первичные частицы которых имеют скорость vo (направление этого вектора не имеет значения, но должно быть фиксировано в процессе расчета конкретного дерева) и вес 1; 6, -
углы между vp и vo . Функцию Кц (vo) следует рассчитывать для дискретного набора скоростей vo = hv,2hv... Nvhv, в котором максимальное значение должно в несколько раз превышать тепловую скорость. Значения Кц (vo) обновляются в
ходе итераций. На каждой итерации ее значения, полученные на предыдущей итерации, используются для учета вклада в сумму в выражении (7) от частиц, существующих на момент обрыва расчета дерева, t = ttree. Пусть на момент обрыва
расчета дерева существует пробная частица со скоростью vp и весом w . Вклад ее и тех частиц, которые она в дальнейшем может породить, в сумму (7) равен
2 w ^ cos2 6-3 j Kц (vp),
где 6 - угол между vp и v0 , а значение K^ (vp) определяется путем интерполяции по массиву значений K^, взятых с предыдущей итерации.
Аналогичным образом рассматривается случай неоднородности температуры и определяется коэффициент теплопроводности; для него
( 2
v0 f0(v0) KX (v0)dv0 ,
x = ^ J
5
2kT 2
v у
1 ^ m (vp )
Kx(vo) «n^cos6t , (8)
а вклад частицы со скоростью vp и весом w вклад в сумму (8) равен
cos 6 wKx (vp).
Отметим, что коэффициенты вязкости и теплопроводности могут рассчитываться совместно, с использованием одних и тех же данных расчетов деревьев частиц.
В итерационном варианте метода исключается погрешность, связанная с конечностью ttree. Величина ttree влияет только на скорость сходимости итерационного процесса. Последняя растет с увеличением ttree, однако вместе с этим растут и затраты на выполнение итерации. Оптимальной с точки зрения вычислительной скорости является такая величина ttree, при которой в среднем в дереве
происходит одно разветвление.
Итерационный вариант метода отличается наличием погрешности, связанной с конечностью массивов значений K^, Kx . В ходе тестирования алгоритма было
выяснено, что при использовании кубического сплайна для интерполяции значений K^, Kx результат практически перестает зависеть от количества элементов
массива уже при количестве элементов выше 3o. При этом увеличение количества элементов массивов не вызывает роста длительности расчета, если критерием окончания расчета выбрать достижение требуемой погрешности статистического усреднения. Дело в том, что с ростом количества элементов растет и объем статистики, набираемый в процессе выполнения единичной итерации. Поэтому были зафиксированы следующие параметры алгоритма: количество элементов массивов K^, Kx - 100, диапазон скоростей для определения этих функций - от 0
до 5>/2kT / m .
Результаты расчетов коэффициентов переноса в газе твердых сфер с помощью итерационного алгоритма приведены в табл. 1. Расчетные результаты приведены в отношении к аналитическим, полученным в приближении Чепме-на-Энскога при пяти слагаемых в разложении по полиномам Сонина [7], рассматриваемым как эталонные. Расчеты проводились до достижения погрешности статистического усреднения 0,02 % при доверительной вероятности 99 %. Приведены результаты, полученные для различных величин временного шага, используемых при расчете свободного пробега. Длительность временного шага характеризуется средней вероятностью столкновения на протяжении шага. Результаты демонстрируют очень быструю сходимость по величине временного шага. Фактически, уже при максимальной из рассмотренных величин шага погрешность составляет доли процента. Возможным объяснением этого явления является то, что корректное распределение по длительности свободного пробега и скоростям фоновой частицы реализуется не для каждой конкретной частицы в отдельности, а для большого количества моделируемых частиц. Расчетные результаты свидетельствуют о возможности использования метода для высокоточных расчетов коэффициентов переноса.
Таблица 1 / Table 1
Коэффициенты переноса в газе твердых сфер, в отношении к эталонным значениям. Рсоц - приведенный временной шаг (средняя вероятность столкновения в течение временного шага) Transport coefficients in a hard sphere gas in relation to the reference values. Pcon is the reduced time step (mean probability of collision
during a time step)
Pcoll ^ расч/^эксп ^расч/ ^эксп
Приближение Чепмена-Энскога, один полином Сонина
- 0,984 0,977
Метод пробных частиц
0,6 1,0018 1,0025
0,4 1,0009 1,0021
0,2 0,9997 1,0014
0,1 0,9989 1,0003
0,05 0,9989 1,0004
Получение результатов для минимального временного шага, с погрешностью статистического усреднения 0,02 %, потребовало 14 часов расчетного времени на процессорном ядре с частотой 2 ГГц. Получение результата с погрешностью усреднения 0,2 %, что примерно на порядок меньше погрешности метода Чепме-на-Энскога при использовании одного полинома Сонина, возможно при максимальной величине временного шага и в 100 раз меньшем объеме статистики, что соответствует расчетному времени 42 секунды.
Описанный метод применим к случаю газа с более сложным потенциалом межчастичного взаимодействия. В сравнении с потенциалом твердых сфер случай потенциалов взаимодействия, не обращающихся в нуль при конечных расстояниях, отличается необходимостью искусственно ограничивать радиус взаимодействия некоторой величиной ртах. При определении вероятности столкно-
вения на временном шаге используется полное сечение столкновения, равное лРшах • Конечность pmax порождает дополнительную погрешность результата.
Полный список источников погрешности результата включает следующее.
1. Конечность статистической выборки.
2. Конечность количества временных шагов за время свободного пробега.
3. Конечность pmax.
4. Погрешность численного интегрирования при расчете угла рассеяния.
5. Конечность массивов значений K^, .
Как уже отмечалось, вклад п. 5 может быть сделан сколь угодно малым без увеличения объемов вычислений. Удобным способом оценки вклада п. 3, 4 является использование расчетных процедур для определения угла рассеяния в традиционном методе расчета коэффициентов переноса - путем расчета интегралов
l s
столкновений Q ' . Как показывают расчетные результаты, вклад этих источников легко может быть снижен до величин менее 0,01 %. В случае потенциала Леннард-Джонса для этого достаточно выбрать pmax > 4ст и количество интервалов численного интегрирования не менее 40.
В табл. 2 приведены результаты расчета коэффициентов переноса в аргоне для температуры 300 К. Межатомное взаимодействие описывалось потенциалом Лен-
нард-Джонсом с параметрами ст = 3,42-10"10 м, е/k = 117,8 К [8]. Погрешность статистического усреднения составляет 0,1 % при доверительной вероятности 99 %. Расчетные результаты приведены в отношении к экспериментальным значениям коэффициентов переноса: ц = 2,272 -10"5 Па ■ с, Х = 0,01781 Вт/м ■ К [9]. Как и в случае с потенциалом твердых сфер, наблюдается достаточно быстрая сходимость по параметру точности метода - средней вероятности столкновения на временном шаге. Отклонение расчетного результата от экспериментальных данных составляет около 1 %. Как показали результаты расчетов для случая потенциала твердых сфер, погрешность метода достаточно мала и, следовательно, отмеченное отклонение можно связать в первую очередь с неточностью выражения потенциала Леннард-Джонса.
Таблица 2/Table 2
Расчетные коэффициенты переноса в отношении к экспериментальным (аргон, в расчетах используется потенциал
Леннард-Джонса) Transport coefficients in a hard sphere gas in relation to the measured values (argon, Lennard-Jones potential was used in calculations)
Pcoll ц расч/цэксп ^■расч/ ^эксп
Приближение Чепмена-Энскога, один полином Сонина
- 1,017 1,008
Метод пробных частиц
0,4 1,017 1,012
0,2 1,014 1,010
0,1 1,013 1,010
0,05 1,013 1,009
Заключение
Метод расчета коэффициентов переноса на основе модифицированного алгоритма пробных частиц продемонстрировал свою работоспособность. Как показал опыт использования метода, с его помощью возможно рассчитывать коэффициенты вязкости и теплопроводности однокомпонентного газа за приемлемое время с точностью, сопоставимой с точностью традиционного метода или превышающей ее. Аналогично известным и широко используемым методам пробных частиц разработанный метод допускает применение параллельных вычислений на неограниченном числе процессоров. Сравнительно небольшой необходимый объем памяти для реализации алгоритма - порядка нескольких сотен байт - делает возможным использование графических ускорителей (GPU).
Обсудим возможные направления обобщения метода и его приложения. Обобщение метода на случай произвольного числа компонент может быть выполнено на основе следующих соображений:
- наличие неоднородности является источником возмущения функции распределения в каждой из компонент;
- источник возмущения в одной из компонент вызывает возмущение в функциях распределения всех компонент, и каждое дерево частиц может включать частицы всех сортов. Соответственно вместо единичной функции K^ необходимо рассчитывать массив таких функций размерностью Nc х Nc, где Nc - количество компонент; то же относится и к K-k ;
- необходимо рассматривать дополнительный тип возмущения - градиент массовой доли компоненты, приводящий к возникновению потоков массы и энергии.
Обобщение на случай частиц, в столкновениях которых происходит изменение внутреннего состояния (вращательного и колебательного), также не создает принципиальных трудностей, особенно в более простом варианте метода, который позволяет рассчитывать коэффициенты переноса с точностью порядка 10 %. Для этого будет достаточным задавать случайное внутреннее состояние первичных частиц деревьев и частиц фона, с которыми происходят столкновения пробных частиц. В более точном варианте метода, итерационном, необходимо определять K^, K^ как функции не только модуля скорости, но и внутреннего состояния первичной частицы.
В качестве возможных приложений метода следует, в первую очередь, назвать задачи определения коэффициентов переноса в многокомпонентных смесях. Аналитические выражения для определения коэффициентов переноса, выведенные в теории Чепмена-Энскога, существенно усложняются уже при переходе от одно-компонентного газа к двухкомпонентному. Вследствие этого вместо непосредственного применения теории Чепмена-Энскога широко используются разнообразные правила смесей, связывающих значения коэффициентов переноса в смеси с их значениями в чистых газах [10]. Представленный метод, после его обобщения на случай многокомпонентных смесей, можно рекомендовать для независимой про-
f 1
верки результатов, получаемых как в методах, основанными на разложении j в
ряд, так и с использованием правил смесей.
Еще одной возможной областью приложения метода является определение параметров межчастичного взаимодействия на основе данных о коэффициентах
переноса. При решении этой задачи необходимо как можно точнее рассчитывать коэффициенты переноса при заданном потенциале взаимодействия. Метод пробных частиц позволяет достичь любой требуемой точности путем увеличения объема вычислений; наибольший вклад в погрешность имеет статистическую природу.
ЛИТЕРАТУРА
1. Рудяк В.Я. Статистическая теория диссипативных процессов в газах и жидкостях. -Новосибирск: Наука, 1987. - 272 с.
2. Митчнер М., Кругер Ч. Частично ионизованные газы. - М.: Мир, 1976. - 497 с.
3. Космические ядерные энергетические установки суб- и мегаваттного класса. Ч. 1 - концепции реакторов (обзор) / Ю.Г. Драгунов, Б.А. Габараев, В.В. Ужанова, М.С. Беляков, М.М. Селиверстов // Проблемы машиностроения и автоматизации. - 2014. - № 2. -С. 95-107.
4. Рудяк В.Я., Лежнев Е.В. Стохастическое моделирование вязкости разреженных газов // Доклады АН ВШ РФ. - 2015. - № 3 (28). - С. 99-108. - doi: 10.17212/1727-2769-2015-399-108.
5. Белоцерковский О.М., Хлопков Ю.И. Методы Монте-Карло в прикладной математике и вычислительной аэродинамике // Журнал вычислительной математики и математической физики. - 2006. - Т. 46, № 8. - С. 1494-1518.
6. Котов В.Э. Методика учета движения фона при моделировании транспорта частиц методом Монте-Карло // Вопросы атомной науки и техники. Серия: Термоядерный синтез. - 2002. - Вып. 3/4. - С. 169-175.
7. Haviland J.K. The solution of two molecular flow problems by the Monte Carlo method // Methods of Computational Physics. - 1965. - Vol. 4. - P. 109-209.
8. White J.A. Lennard-Jones as a model for argon and test of extended renormalization group calculations // Journal of Chemical Physics. - 1999. - Vol. 111, N 20. - P. 9352-9356. -doi: 10.1063/1.479848.
9. Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermo-physical property Library CoolProp / I.H. Bell, J. Wronski, S. Quoilin, V. Lemort // Industrial & Engineering Chemistry Research. - 2014. - Vol. 53, N 6. - P. 2498-2508. -doi: 10.1021/ie4033999.
10. Reid R.C., Prausnitz J.M., Sherwood T.K. The properties of gases and liquids. -New York: McGraw-Hill, 1977. - 688 p.
APPLICATION OF THE MONTE CARLO TEST PARTICLE METHOD TO CALCTLATE TRANSPORT COEFFICIENS
IN DILUTED GASES
Finnikov K.A.
Siberian Federal University, Russia
Conventional methods of calculating transport coefficients in diluted gases are based on theexpansion of the first-order Chapman-Enskog velocity distribution function. Alternative methods based on statistical simulation of particle motion can be useful in some cases. An example of such a method is proposed in the present paper. It is a version of the Monte-Carlo test particle method modified in such a way that the density of test particles approximates to the first-order velocity distribution function. The main features of the method are as follows. First, a statistical weight as an additional characteristic of a test particle is introduced. This statistical weight is taken into account as a multiplier in determining the contribution of a particle to the velocity distribution function. Second, it is a special algorithm of simulating collisions between test particles and background particles that reflects the structure of the linearized collision term in the equation for the first approximation of the velocity distribution function. In this algorithm a collision leads
to colliding test particle vanishing and the appearance of three new test particles. As a result, instead of a single test particle trajectory, a ternary tree of trajectories is calculated. Trajectories of test particles are interdependent within a tree of trajectories, but different trees are calculated separately. This provides an unlimited scalability of the algorithm in multiprocessor systems similarly to other test particle methods. The case of a monocomponent gas with a central potential of interparticle interaction (hard sphere and Lennard-Jones potentials) is examined. The results obtained with use of the proposed method for hard sphere gas are identical to the known analytical results. Calculated transport coefficients in an inert gas described by the Lennard-Jones interatomic potential are compared with experimental data. The disagreement is about 1 %, which corresponds to the accuracy of the approximation of interatomic interaction by the Lennard-Jones expression. The proposed method can be used both independently and in combination with conventional methods to cross-check calculation errors.
Keywords: transport coefficients; diluted gas; test particle method; Chapman-Enskog expansion; viscosity; thermal conductivity; Monte-Carlo simulation.
DOI: 10.17212/1727-2769-2016-2-108-122
REFERENCES
1. Rudyak V.Ya. Statisticheskaya teoriya dissipativnykh protsessov v gazakh i zhidkostyakh [Statistical theory of dissipative processes in gases and liquids]. Novosibirsk, Nauka Publ., 1987. 272 p.
2. Mitchner M., Kruger Ch.H. Partially ionized gases. New York, A Wiley-Interscience Publ., 1973. 518 p. (Russ. ed.: Mitchner M., Kruger Ch. Chastichno ionizovannye gazy. Moscow, Mir Publ., 1976. 497 p.).
3. Dragunov Yu.G., Gabaraev B.A., Uzhanova V.V., Belyakov M.S., Seliverstov M.M. Kos-micheskie yadernye energeticheskie ustanovki sub- i megavattnogo klassa. Ch. 1 -kontseptsii reaktorov (obzor) [Sub-megawatt and megawatt class space nuclear power facilities. Pt. 1 - conceptions of reactors (review)]. Problemy mashinostroeniia i avtomatizatsii -Engineering & Automation Problems, 2014, no. 2, pp. 95-107.
4. Rudyak V.Ya., Lezhnev E.V. Stokhasticheskoe modelirovanie vyazkosti raz-rezhennykh gazov [Stochastic simulation of rarefied gas viscosity]. Doklady Akademii nauk vysshei shko-ly Rossiiskoi Federatsii - Proceedings of the Russian higher school Academy of sciences, 2015, no. 3 (28), pp. 99-108. doi: 10.17212/1727-2769-2015-3-99-108
5. Belotserkovskii O.M., Khlopkov Yu.I. Monte Carlo methods in applied mathematics and computational aerodynamics. Computational Mathematics and Mathematical Physics, 2006, vol. 46, no. 8, pp. 1418-1441. doi: 10.1134/S0965542506080124. Translated from Zhurnal vysshei matematiki i matematicheskoi fiziki, 2006, vol. 46, no. 8, pp. 1494-1518.
6. Kotov V.E. Metodika ucheta dvizheniya fona pri modelirovanii transporta chastits metodom Monte-Karlo [Method for taking into account the motion of background in Monte-Carlo simulation of particle transport]. Voprosy atomnoi nauki i tekhniki. Seriya: Termoyadernyi sintez - Problems of Atomic Science and Engineering. Series "Thermonuclear Fusion", 2002, iss. 3-4, pp. 169-175.
7. Haviland J.K. The solution of two molecular flow problems by the Monte-Carlo method. Methods of Computational Physics, 1965, vol. 4, pp. 109-209.
8. White J.A. Lennard-Jones as a model for argon and test of extended renormalization group calculations. Journal of Chemical Physics, 1999, vol. 111, no. 20, pp. 9352-9356. doi: 10.1063/1.479848
9. Bell I.H., Wronski J., Quoilin S., Lemort V. Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property Library CoolProp. Industrial & Engineering Chemistry Research, 2014, vol. 53, no. 6, pp. 2498-2508. doi: 10.1021/ie4033999
10. Reid R.C., Prausnitz J.M., Sherwood T.K. The properties of gases and liquids. New York, McGraw-Hill, 1977. 688 p.
СВЕДЕНИЯ ОБ АВТОРАХ
Финников Константин Андреевич - родился в 1974 году, кандидат физико-математических наук, доцент кафедры теплофизики Сибирского федерального университета. Область научных интересов: физика низкотемпературной плазмы, вычислительная гидродинамика, методы статистического моделирования. Опубликовано 48 научных работ. (Адрес: 660041, Красноярск, пр. Свободный, 79. Email: [email protected]).
Finnikov Konstantin Andreevich (b. 1974) - PhD, Associated Professor, Thermophysics department, Siberian Federal University. His research interests are currently focused on low-temperature plasma physics, computational fluid dynamics, Monte Carlo simulation methods. He is author of 48 scientific papers. (Address: 79 Svobodny ave., Krasnoyarsk, 660041, Russia. Email: [email protected]).
Статья поступила 21 апреля 2016 г.
Received April 21, 2016
To Reference:
Finnikov K.A. Primenenie metoda probnykh chastits dlya vychisleniya koeffitsientov perenosa v razrezhennykh gazakh [Application of test particle Monte Carlo method to calculation of transport coefficients of diluted gases]. Doklady Akademii nauk vysshei shkoly Rossiiskoi Fe-deratsii - Proceedings of the Russian higher school Academy of sciences, 2016, no. 2 (31),
pp. 108-122. doi: 10.17212/1727-2769-2016-2-108-122
