ТЕПЛОФИЗИКА ВЫСОКИХ ТЕМПЕРАТУР, 2008, том 46, № 5, с. 664-673
ТЕПЛОФИЗИЧЕСКИЕ СВОЙСТВА ВЕЩЕСТВА
УДК 533.1
МОДЕЛИРОВАНИЕ СКОРОСТИ ДИССОЦИАЦИИ И ВРЕМЕНИ КОЛЕБАТЕЛЬНОЙ РЕЛАКСАЦИИ ДВУХАТОМНЫХ МОЛЕКУЛ МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ
© 2008 г. Д. В. Котов, С. Т. Суржиков
Институт проблем механики РАН, Москва Поступила в редакцию 23.07.2007 г.
Излагается метод расчета скорости диссоциации и времени колебательной релаксации двухатомных молекул N2 и O2 методом молекулярной динамики, в котором реализуется одновременное моделирование всего ансамбля частиц на основе классических представлений о динамике молекул с внутренними степенями свободы. Обсуждается соответствие результатов расчета экспериментальным данным.
PACS: 52.20.Hv,52.65.Rr
ВВЕДЕНИЕ
Во многих явлениях физической газовой динамики при поступательных температурах частиц, достигающих десятков тысяч градусов, существенную роль играют неравновесные процессы физико-химической кинетики, для которых характерно сильное различие между равновесными распределениями атомов и молекул по поступательным и внутренним степеням свободы, а также отсутствие равновесных распределений в пределах отдельно взятых степеней свободы. Такая ситуация часто встречается при изучении структуры сильных ударных волн, при гиперзвуковом обтекании летательных аппаратов, в газовых разрядах и при взаимодействии мощного лазерного излучения с веществом.
Расчетно-теоретическое описание неравновесных физико-химических процессов во многом базируется на моделях физической и химической кинетики газофазных химических реакций и моделях колебательной и электронной релаксации молекулярных компонент газовых смесей.
Разработано значительное число теоретических моделей, позволяющих достичь удовлетворительного описания имеющихся экспериментальных данных. Такие модели различаются прежде всего степенью подробности учета кинетических процессов (однотемпературная и многотемпературная кинетика, поуровневая кинетика и т.п.), а также различными способами описания элементарных энергообменных процессов в молекулах. К общему недостатку большинства теоретических моделей относится их полуэмпирический характер, что затрудняет процедуру экстраполяции кинетических данных в область высоких температур. Напомним, что подавляющее число
экспериментальных данных ограничено сверху температурами T = 8000-10000 К, а практика аэрофизических расчетов сталкивается с температурами в несколько десятков тысяч градусов.
Принимая во внимание исключительную сложность получения экспериментальных данных и большое разнообразие процессов, подлежащих учету при высоких температурах, следует подчеркнуть актуальность разработки расчетно-теорети-ческих моделей, основанных на ab-initio принципах моделирования, т.е. использующих лишь исходную спектроскопическую информацию о структуре частиц.
К таким методам относится развитый в данной работе метод молекулярной динамики. За более чем сорокалетнюю историю своего развития он стал достаточно надежным инструментом исследования микроскопической структуры материалов и физико-химических процессов в молекулярной физике. В последнее время этот метод является одним из основных подходов к изучению веществ на наноуровне [1].
Применительно к задачам физико-химической механики и химической кинетики можно выделить два основных направления развития молеку-лярно-динамического подхода. В работах [2-5] развит метод классических траекторий (МКТ), с использованием которого вычисляются скорости элементарных химических превращений и константы поуровневой кинетики. В основе МКТ лежит вычислительный эксперимент по изучению статистических закономерностей бомбардировки частицы-мишени потоком пробных частиц. И те и другие частицы подразумеваются структурными (т.е. обладающими внутренними степенями свободы), так что в процессе взаимодействия ча-
стиц происходит перераспределение их кинетической энергии по внутренним степеням свободы.
Результатом таких расчетов [2, 4, 5] является, как правило, определение сечения реакции а(Ег, V, I), которое представляет собой многомерный интеграл по варьируемым параметрам. Например, для трехатомной системы А + ВС это пятимерный интеграл
СТ(Е„ V, I) = ПЬ^х < Р(Е„ V, I)>,
< P(Et, v, J )> =
(2 п) 3bL
J bdb J si
sin 6 dQ x
2 п 2п
2п
x
J dфJ dnJPr(Et, v, J, b, 6, d^,
где Ь - прицельный параметр; углы ф и б задают ориентацию молекулы ВС; - начальная фаза колебаний молекулы ВС; п - угол, определяющий ориентацию вектора углового момента молекулы ВС; V и I - соответственно колебательные и вращательные квантовые числа молекулы ВС; ЕI - относительная энергия столкновения; РГ(Е, V, I, Ь, б, ф, п, - вероятность реакции при заданных параметрах, которая принимает значение 1, если реакция произошла, и 0 - в противном случае. Ввиду большой размерности интеграла для его численного расчета используется метод Монте-Карло. Уровневые константы скорости к(Т, V) получаются усреднением сечения реакции по максвелловскому распределению для поступательной энергии и по больцмановскому распределению для вращательных состояний.
Второй из упомянутых подходов, т.е. прямое моделирование динамики молекул в их ансамбле, применяется в данной работе. Имеющийся опыт использования такого подхода [6-8] свидетельствует о его перспективности, в особенности с учетом непрерывно возрастающей производительности компьютеров.
Типичная постановка задачи, решаемой с использованием второго подхода, состоит в следующем. Рассматривается термостат, в котором заключено N взаимодействующих частиц. Химические реакции и обмен энергией между поступательными и внутренними степенями свободы моделируются одновременным интегрированием траекторий этих частиц. Результатом расчетов является зависимость параметров смеси (степень протекания реакции, распределение молекул по квантовым состояниям и др.) от времени. При необходимости от этих данных можно прийти к константам скоростей реакций, времени релаксации и другим характеристикам процесса.
Если конечной целью является определение констант скорости реакций, то эффективность
данного метода по затратам компьютерного времени значительно ниже, чем при вычислении сечений процессов посредством MKT. Однако рассматриваемая модель более приближена к реальной ситуации в газовой смеси, поскольку не использует максвелл-больцмановское усреднение и не предполагает вычисления числа столкновений, а непосредственно моделирует каждое из них. Поэтому модель лучше описывает химическую кинетику в системах, где нет максвеллов-ского распределения по скоростям, например, в разрядах, газовых лазерах и верхних слоях атмосферы, а также способна учитывать различные коллективные эффекты.
1. ОПИСАНИЕ ИСПОЛЬЗУЕМОГО МЕТОДА МОЛЕКУЛЯРНОЙ ДИНАМИКИ
Дана газовая смесь, состоящая из N частиц (атомов и двухатомных молекул), которые могут вступать в химические реакции. Граничные условия - периодические, т.е. вылет частицы через одну из границ компенсируется появлением аналогичной частицы на симметричной границе. Рассматриваются только парные взаимодействия: взаимодействие атомов внутри одной молекулы и взаимодействие атомов, принадлежащих разным молекулам. Задача решается в рамках квазиклассического описания динамики частиц. Считается, что внутримолекулярное взаимодействие описывается потенциалом Морзе
TU ) т-ч -2 в( r - re) -e(r - г.)
U (r) = Dee -2 Dee , (1)
где De - энергия диссоциации, r - расстояние между атомами в молекуле, re - равновесное расстояние между ними.
Межмолекулярное взаимодействие описывается потенциалом Леннарда-Джонса
U( Гц) = 4£;
12
6-,
(2)
где £у, - параметры этого потенциала, гц - расстояние между атомами г и ц, не связанными в молекулу.
Требуется определить распределение молекул по колебательным уровням энергии, степень и скорость диссоциации а(г) и другие параметры газовой смеси, которые будут обсуждаться ниже.
В данной модели сразу все частицы одновременно взаимодействуют друг с другом, поэтому поверхность потенциальной энергии (ППЭ) зависит от N координат. В связи со сложностью более или менее точного построения такого потенциала в модели выбран самый простой вид ППЭ - суперпозиция потенциалов Морзе и Леннарда-Джонса. Поэтому следует ожидать, что один из существенных вкладов в погрешность расчетов связан
b
п
0
0
0
0
0
с неточностью задания ППЭ. Впрочем, отметим, что проблема задания потенциальной поверхности для взаимодействующих частиц в методах молекулярной динамики является одной из ключевых и находящихся в процессе разработки. Кроме этого, на точность расчетов будет сильно влиять статистическая погрешность, так как число моделируемых частиц N остается по-прежнему не столь большим даже при многопроцессорных расчетах. Еще один вклад в погрешность будет давать численное интегрирование траекторий, но эту ошибку всегда можно контролировать, меняя шаг интегрирования.
Перед началом компьютерного моделирования динамики ансамбля частиц производится инициализация системы. Начальные координаты каждой частицы определяются случайным образом, так чтобы расстояние межу ними превосходило некоторый наперед заданный минимальный радиус межмолекулярного взаимодействия
Скорость центра масс частицы уст разыгрывается в соответствии с максвелловским распределением так, чтобы среднеквадратичная скорость
была равна vct =
3 kT M
, где M - масса частицы,
т.е. свободного атома или молекулы. В последнем случае М = m1 + m2, где m1 и m2 - массы атомов, из которых состоит данная молекула; Т - начальная поступательная температура системы.
Направление движения частицы выбирается в соответствии с равновероятным распределением:
V с
f \
cos ф cos б sin ф cos б sin б
; ф = 2
П = аге^ (2%2-1 )е [-п/2; п/2],
где и %2 - случайные числа, равномерно распределенные на отрезке [0, 1]; ф, б - углы, характеризующие направление скорости центра масс частицы. Для двухатомных молекул выбирается также случайное (равновероятное) направление, вдоль которого ориентируется ось молекулы.
Для каждой молекулы в соответствии с больц-мановским распределением определяются электронное состояние, колебательная еуШ и вращате
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.