CC BY
41
Вычислительные технологии
Том 16, № 2, 2011
Модификации метода PIC-MCC для моделирования многокомпонентной плазмы емкостного ВЧ-разряда*
Д. А. Арискин, И. В. Швейгерт Институт теоретической и прикладной механики СО РАН, Новосибирск, Россия
e-mail: dmitry.ariskin@gmail.com
Представлены кинетическая и гибридная модели многокомпонентной плазмы высокочастотного емкостного разряда в химически активных газах. Предложены модификации метода частиц в ячейках, позволяющие сократить время расчета.
Ключевые слова: ВЧЕ-разряд, многокомпонентная плазма, кинетическая модель, гибридная модель, метод частиц в ячейках, распараллеливание.
Введение
Газовый разряд в смесях молекулярных газов, таких как SiH3, СН4, С2Н2, широко используется для создания различного рода пленок и покрытий. Пленки, получаемые в углеводородных смесях, благодаря их исключительным свойствам — высокой прочности, химической инертности, низкому значению абразивного коэффициента, прозрачности в инфракрасном диапазоне — имеют огромное значение для многих приложений, включая создание трибологических материалов, пассивирующих слоев и катодов с низкой работой выхода электронов [1-3]. Для выращивания таких пленок широко используется высокочастотный емкостной (ВЧЕ) разряд.
В результате плазмохимичееких реакций в ВЧЕ-разряде в смеси, содержащей молекулярные соединения, образуются стабильные и активные компоненты, которые осаждаются на подложку и стенки камеры и формируют наноразмерные частицы в объеме разряда. Плазма в этом случае является многокомпонентной, поскольку кроме электронов в разряде присутствуют несколько различных типов как положительных, так и отрицательных ионов. Химические реакции приводят также к изменению состава фонового газа. Таким образом, задача моделирования такой системы усложняется в силу необходимости расчета большого количества компонентов и взаимодействий между ними.
Для моделирования многокомпонентной плазмы ВЧЕ-разряда используются кинетические, гидродинамические и гибридные методы. При кинетическом подходе для описания электронного и ионного компонентов плазмы применяется уравнение Больцмана, которое может решаться различными методами [4], например, с помощью двучленного приближения [5]. Кроме того для решения уравнения Больцмана широко используется статистический метод частиц в ячейках с розыгрышем столкновений методом Монте-Карло (Particle in cell Monte Carlo collisions method (PIC-MCC)) [6, 7]. В настоящее время метод PIC-MCC является наиболее мощным инструментом теоретического исследования газоразрядной плазмы и позволяет получать решение ab initio. Однако
* Работа выполнена при финансовой поддержке РФФИ (грант № 08-02-00833а) и Интеграционного проекта СО РАН № 113-2009.
недостатком кинетических методов является их вычислительная трудоемкость. Применительно к моделированию многокомпонентной плазмы затраты растут пропорционально количеству рассматриваемых ионных компонентов. Для решения этой проблемы в данном случае, как правило, пренебрегают изменениями в составе фонового газа [8], что позволяет рассматривать фиксированный набор столкновительных процессов, В ряде работ для упрощения задачи рассматриваются только основные типы ионов (см., например, [9]) или система ионов заменяется одним эффективным ионом с усредненными характеристиками [10]. Такие приближения неприменимы при детальном моделировании плазмохимичееких процессов.
Гидродинамические методы описывают систему с помощью одного или нескольких моментов уравнения Больцмана для каждого типа частиц, рассматриваемых в задаче [11]. При таком подходе задача сводится к решению системы дифференциальных уравнений в частных производных. Достоинством этого метода является его вычислительная эффективность. Вместе с тем область применения данного метода ограничена системами, в которых длина свободного пробега ионов и электронов много меньше характерных размеров задачи. Гидродинамические модели не позволяют также получать информацию о функции распределения заряженных частиц по энергии.
Кроме описанных подходов широко используются гибридные методы, объединяющие кинетический и гидродинамический подходы. Так, в работах [8, 12] метод Монте-Карло применялся для расчета скоростей ионизации и диссоциации, а также транспортных коэффициентов для электронов. Эти коэффициенты в дальнейшем использовались в гидродинамической модели. Недостаток такого подхода в том, что распределение электронов по энергии считается функцией локального электрического поля, и, таким образом, модель не является полностью самосогласованной. Как следствие, применимость данного метода ограничена областью более высоких давлений газа, поскольку при низких давлениях неравновесность функции распределения электронов по энергии играет существенную роль. Полностью самосогласованная гибридная модель для описания ВЧЕ-разряда была предложена в работе [13], где кинетические уравнения решаются одновременно с гидродинамическими, С использованием этой модели удалось избежать численного нагрева электронов на статистических флуктуациях электрического поля,
В настоящей работе дано описание двух моделей многокомпонентной плазмы, кинетической и гидродинамической, использованных ранее [14, 15] для моделирования ВЧЕ-разряда в смеси ацетилена и аргона, В первой кинетической модели электроны и все типы ионов описываются уравнениями Больцмана, которые решаются методом PIC-MCC, Во второй гибридной модели электроны описываются в кинетическом приближении, а ионы — в гидродинамическом, что позволило существенно сократить время расчетов. Обе модели учитывали химические процессы в разряде и процессы прокачки и перемешивания. Предложены также модификации метода PIC-MCC для многокомпонентных систем и методы оптимизации алгоритма расчета. Проводится сравнение результатов расчетов, полученных с использованием различных моделей.
Рассматривается одномерная модель ВЧЕ-разряда, но все предложенные модификации метода частиц в ячейках применимы и к двумерным моделям.
Приведенные результаты расчетов получены на двухпроцессорной (восьмнядерной) машине 2 х Intel(R) Xeon(R) Е54208 (4 cores) х 2.43 ГГц, 16 GB RAM.
Материал статьи изложен следующим образом. В разделе 1 представлена кинетическая модель многокомпонентной плазмы, в разделе 2 дано описание гибридного ме-
тода, позволяющего существенно уменьшить время расчетов при моделировании многокомпонентной плазмы по сравнению с кинетическим алгоритмом. Раздел 3 посвящен проблеме учета химических процессов в рамках кинетической модели, в разделе 4 описаны предложенные модификации кинетического алгоритма, необходимые для ускорения расчетов и повышения точности результатов, в разделе 5 предложен эффективный метод расчетов задач с учетом вторичной электронной эмиссии, в разделе 6 описаны результаты распараллеливания системы с помощью библиотеки ОрепМР.
1. Кинетическая модель многокомпонентной плазмы
Рассмотрим одномерную модель ВЧЕ-разряда, Для применимости одномерного приближения электроды должны быть симметричными и их характерные размеры существенно большими межэлектродного расстояния. Один электрод считается заземленным, ко второму приложено напряжение с заданными амплитудой и частотой. Внешняя цепь не учитывается. Электроны и все типы ионов, как положительные так и отрицательные, описываются с помощью уравнений Больцмана для функции распределения
дг дг ту дv
■, (1)
где г V ту — радиус-вектор, скорость и масса соответствующей частицы, Е — электрическое поле, — интеграл столкновений, для электронов включающий упругие и неупругие столкновения с атомами и молекулами смеси, в которой горит разряд. Для ионов учитываются только упругие столкновения и резонансная перезарядка. Парные кулоновекие взаимодействия не рассматриваются, так как для емкостного ВЧ-разряда характерна низкая степень ионизации, и, как следствие, частота таких столкновений много меньше частоты столкновений с нейтральными компонентами.
Неупругие взаимодействия различных типов ионов с нейтральными компонентами и другими ионами описываются уравнениями химической кинетики
= ^ ^ кт>1ПтП1 — У ^ кр^ПрЩ, (2)
где кт,1 — скорости реакций, продуктом которых является ион с индексом г, кр^ — скорости реакций, для которых ион с индексом г является реагентом. Эти уравнения включают, в частности, реакции рекомбинации между положительными и отрицательными ионами.
Электрическое поле описывается уравнением Пуассона
Аф = 47г I епе — У^ Чгп ■• (3)
Транспорт нейтральных компонентов описывается с помощью уравнений диффузии дпг,
дг
У(£гаУ(п„)) — ^ ] кт,1птп1 "У ] кр,ппрпп + ^геас^у, (4)
где кт,1 — скорости реакций, продуктом которых является нейтральный компонент с индексом п, кр,п — скорости реакций, для которых нейтральный компонент с индексом п является реагентом, Оп — коэффициент диффузии. Коэффициенты диффузии
для различных компонентов смеси определяются из парных коэффициентов диффузии для каждого компонента фонового газа, Бгеас1;а- описывает изменение концентрации компонента за счет процессов прокачки и перемешивания:
Бгеасt,j — + Брцт\>,] + Бт1х^ , (5)
где БйБритр,г, Бт;x,j — скорости изменения концентрации соответственно за счет входного потока газа, оттока и перемешивания газа между областью реактора, запятой разрядом, и областью, не запятой им. Задача решается до установления равновесия между всеми компонентами системы. Под равновесием здесь понимается состояние, в котором с продолжением счета параметры плазмы не меняются при сравнении распределений но различным периодам ВЧ-разряда,
В методе Р1С-МСС обычно используются расщепление но физическим процессам и шаг но времени дня ионов, гораздо больший шага дня электронов |6, 161, В результате при моделировании двухкомнопептпой плазмы можно практически пренебречь затратами па расчет ионного компонента. Однако при моделировании многокомпонентной плазмы, когда число различного тина ионов увеличивается, затраты на расчет ионного компонента становятся значительными. На рис, 1 приведены зависимость времени расчета 1000 периодов разряда дня различных вариантов задачи от количества рассматриваемых типов ионов, а также доля времени, затраченная па расчет движения электронов и их столкновительных процессов. Данные получены дня ВЧЕ-разряда в смеси аргона с ацетиленом при давлении Р — 70 мТорр, амплитуде приложенного напряжения и0 — 90 В и межэлектродном растоянии й — 7 см. Отношение между шагом по времени для ионов и электронов в этих расчетах Д^/ ДЬе — 30, количество псевдочастиц па каждый компонент плазмы — 20 000, Следует также отметить, что включение в модель дополнительных типов ионов и химических реакций ведет к увеличению времени установления равновесия в системе и, как следствие, к увеличению времени расчетов, по это является свойством физической системы, а не алгоритма. В данном случае пас интересует зависимость трудоемкости временного шага от количества рассматриваемых типов ионов Л^.
г, с 6000
5000
4000
3000
2000
1000
2 4 6 8 10 12 14 16 18 20 22 Ц
Рис. 1. Время расчета 1000 периодов разряда в зависимости от количества рассматриваемых типов ионов N¿5 серые колонки — кинетическая модель, белые — гибридная модель; заштрихованная область время, затраченное на расчет электронного компонента
С увеличением N затраты на расчет динамики ионов линейно растут, и, например, при расчете параметров разряда в смеси аргона с ацетиленом [14], где рассматривалось 22 различных типов ионов, затраты на расчет ионного компонента составили около 60%, На рис, 1 также приведены данные для гибридной модели, которая будет обсуждаться в следующем разделе.
Точность метода Р1С-МСС ограничена прежде всего погрешностями, определяемыми статистическими флуктуациями, связанными с дискретным представлением электронов и ионов. Последние заменяются псевдочастицами, концентрация которых мала по сравнению с реальным количеством электронов и ионов. Уменьшить эту погрешность можно только увеличением количества псевдочастиц. Для проверки точности расчетов методом Р1С-МСС проводились тестовые расчеты со значительно большим количеством псевдочастиц Для рассматриваемых задач оптимальным является
соотношение Мрагк — 50т. е, число псевдочастиц в каждой ячейки должно равняться 50, Отклонение решения, получаемого при увеличении числа Мраг1; в четыре раза, не превышает 3 %. Другой источник погрешности связан с пространственной сеткой. Наличие сетки приводит к погрешностям аппроксимации уравнения Пуассона и интерполяции силы, что в свою очередь может вызывать нефизический разогрев электронов. Снизить погрешность можно уменьшением шага пространственной сетки. Для оценки ее величины были проведены расчеты с вдвое меньшим шагом пространственной сетки (при этом шаг по времени пропорционально уменьшался, а количество псевдочастиц увеличивалось). Отклонение решения от более точного не превышало 7%.
Таким образом, разработана кинетическая модель многокомпонентной плазмы ВЧЕ-разряда в химически активной смеси. Эта модель в дальнейшем используется для проверки точности гибридных моделей.
2. Гибридная модель многокомпонентной плазмы
Ионы имеют существенно меньшую длину свободного пробега, чем электроны, и в ВЧЕ-разряде распределение ионов по энергии существенно отличается от максвелловского только при очень низких давлениях. Таким образом, в широком диапазоне параметров для ионов можно использовать гидродинамическое описание. При низком давлении для первичного установления равновесия возможно применение гидродинамического приближения с последующим уточнением с помощью кинетического моделирования. Кроме того, некоторые типы ионов можно рассматривать в кинетическом приближении (если, например, интерес представляют их функции распределения), а остальные типы ионов — в гидродинамическом,
В гидродинамическом приближении для положительных и отрицательных ионов уравнение Больцмана заменяется законами сохранения — уравнениями неразрывности и импульса
д'Щ, д ПгУг ,
/ у кт,1птп1 — У, кр,гпрпг, (6)
дг дх дьг дьг qiE к в Т дп
дг дх шшп дх
^^то пМ; (7)
где ^тот — эффективная частота передачи импульса. Как и в предыдущей модели, самосогласованно с (6), (7) решаются уравнения (1) для электронов, а также уравнения (3)-(5).
п, см
10'
10
-з
0.5
1.0
1.5
X, см
Рис. 2. Распределения ионов п и электронов пе в ВЧЕ-разряде (13.56 МГц) в чистом аргоне, полученные с помощью кинетической (штрих) и гибридной (штрихпунктир) моделей; параметры разряда: Р = 100 мТорр, д = 2 см, и0 = 250 В
Дня проверки применимости гибридной модели были выполнены расчеты ВЧЕ-нлазмы в аргоне с одним типом ионов. На рис. 2 приведены результаты расчета разряда в аргоне при Р = 100 мТорр с помощью двух описанных выше моделей — кинетической и гибридной. Параметры задачи были взяты из |17|, При указанном достаточно высоком давлении обе модели дают близкие результаты, которые хорошо совпадают с результатами работы |17|, Отклонения таких важных параметров разряда как шири-па приэлектродного слоя, напряженность электрического ноля в нриэлектродном слое, концентрация плазмы в квазииейтралыюй области не превышают 5%.
Как и ожидалось, применение гибридной модели не дает существенного выигрыша в производительности, поскольку ионы в этом случае рассчитываются с существенно большим шагом но времени, чем электроны. Выигрыш во времени составил 40 с при общем времени расчета на установление 50 мин.
Преимущество использования гибридного метода проявляется при увеличении количества рассматриваемых типов ионов. На рис. 1 приведена зависимость времени расчета для кинетического и гибридного методов от Л^. Видно, что для задачи, в которой моделируется 22 тина ионов |14|, затраты на расчет ионного компонента составляют порядка 60% от общего времени расчетов, и применение гибридного метода дает выигрыш в производительности больше чем в 2 раза.
Таким образом, гибридный метод эффективен дня расчета много компонентной плазмы. Достоинством этого метода является также возможность выбора модели (кинетической или гидродинамической) отдельно для каждого тина ионов.
3. Учет химических процессов
При решении уравнений Больцмана (1) дня ионов методом Р1С-МСС и уравнений химической кинетики (2) возникает проблема пересчета распределения ионов дня получения самосогласованного решения. В процессе счета новые веса рассчитывались но формуле
Шк — п/щ Я(хк - ху),
3 = 1
(8)
|хк — ху| > к
(9)
где п3 и п^ — новое и старое значения концентрации в ]-ш узле пространственной сетки, В [9] был использован аналогичный подход для учета реакций рекомбинации в силановой плазме, В этой работе показано, что применение такого метода позволяет уменьшить количество нефизических шумов, вызванных плохой статистикой при моделировании столкновптельных процессов. Отметим также, что нами этот метод был разработан независимо и использовался в работе [15],
Другой важной проблемой, возникающей при совместном решении уравнений Больц-мана и уравнений химической кинетики, является большое различие во временах установления равновесия в плазме и химического равновесия. Если для установления равновесия в плазме с одним типом ионов достаточно примерно 1000 периодов разряда, то установление химического равновесия между различными типами ионов требует 10 ООО и более периодов, особенно в электроотрицательных смесях, где происходит накопление отрицательных ионов. Кроме того, за счет плазмохимических процессов изменяется состав фонового газа, что в свою очередь приводит к изменению функций распределений электронов и ионов, а это также значительно увеличивает время установления равновесия.
Для решения данной проблемы был использован гибридный (с усреднением по электронному компоненту) метод, который является модификацией метода, изложенного в разделе 2, Суть метода состоит в использовании усредненных параметров для электронного компонента плазмы, полученных из предварительных расчетов.
Схема расчета с помощью этого метода следующая:
— сначала осуществляется расчет параметров разряда до установления равновесия с помощью гибридного метода, описанного в разделе 2, Моментом установления равновесия считается период разряда, после которого скорость относительного изменения параметров становится меньше заданной величины;
— затем в течение заданного количества периодов разряда собирается статистика о концентрации электронов, а также вычисляются средние скорости различных столк-новительных процессов (ионизация, диссоциация, прилипание) в те моменты времени периода разряда, когда производится расчет ионного компонента;
— после этого расчет осуществляется с шагом по времени, характерным для ионного компонента, с использованием рассчитанных распределений концентраций электронов и скоростей реакций.
Данный процесс повторяется несколько раз до тех пор, пока изменение в составе фонового газа за один цикл не станет меньше заданного.
При использовании этого метода затраты на расчет электронного компонента уменьшаются примерно в 500 раз и становятся малыми по отношению к затратам на расчет ионного компонента. Расчет ионного компонента можно производить с помощью как кинетической, так и гидродинамической модели.
На рис, 3 приведены зависимость времени расчета 1000 периодов разряда для различных вариантов задачи от количества рассматриваемых типов ионов, а также доля
Рис. 3. Зависимость времени расчета 1000 периодов разряда от количества рассматриваемых типов ионов N1, полученная с использованием гибридной модели с усреднением по электронному компоненту; заштрихованная область время, затраченное на расчет электронного компонента
времени, затраченная на расчет движения электронов и их столкповительных процессов. Расчет производился с помощью гибридного алгоритма с усреднением но электронному компоненту. Данные получены дня задачи, описанной в раздело 1. Видно, что применение кинетической модели дня всех типов ионов является неэффективным, так как при большом количество ионов основное время расчетов тратится на расчет ионного компонента и выигрыш от исключения электронного компонента становится незначительным. При использовании же гидродинамического описания дня ионов, напротив, удается получить значительный выигрыш в производительности (более чем в 10 раз).
Химические процессы имеют, как правило, большие характерные времена установления равновесия но сравнению с временем установления равновесия в двухкомнонептной смеси (электроны, ионы одного тина). Таким образом, включение химических процессов в модель разряда ведет к увеличению времени установления равновесия и времени расчетов в целом. Кроме того, учет химических процессов в рамках кинетической модели требует модификации Р1С-МСС-алгоритма, что увеличивает его трудоемкость. Дня решения этой проблемы можно эффективно использовать гибридный метод с усреднением но электронному компоненту.
4. Переменные веса частиц
Как было показано в предыдущем раздано, при совместном решении уравнений Больц-мапа методом Р1С-МСС и уравнений химической кинетики возникает необходимость пересчета весов нсевдочастиц, что в итоге ведет к формированию набора нсевдочастиц с различными весами. В ряде случаев это может приводить к накоплению большого количества частиц с маленькими весами, которые не вносят заметного вклада в расчет функции распределения. Данная проблема особенно актуальна дня расчета электроотрицательных газов, так как отрицательные ионы заперты в центре разрядного
промежутка и представляющие их псевдочастицы не уходят на электроды. При этом количество псевдочастиц постоянно увеличивается, а их веса вследствие реакций нейтрализации уменьшаются. Чтобы избежать накопления частиц с маленькими весами, производят склеивание этих частиц [17, 18] либо в момент их образования, либо периодически в заданные моменты времени. При нехватке частиц проводится их дробление, В нашей модели был реализован подход, при котором склеивание частиц осуществляется периодически с достаточно большим интервалом времени, В силу того что веса псевдочастиц одного типа могут существенно различаться, склеивать эффективнее всего самые легкие частицы, а дробить — самые тяжелые. При этом в приведенной последовательности выполняются следующие операции:
— частицы сортируются по весу;
— выбирается определенное количество частиц с наименьшим весом для удаления (они расположены в конце массива);
— выбранные частицы сортируются по номеру ячейки сетки;
— для каждой удаляемой частицы ищется "пара" — частица, принадлежащая той же ячейке сетки. Поиск осуществляется с конца массива (по мере возрастания веса). При этом не допускается, чтобы "пары" различных удаляемых частиц совпадали;
— наконец, параметры каждой "пары" заменяются на суммарные параметры обеих частиц и удаляемые частицы исключаются их массива частиц.
Трудоемкость данной операции составляет 0(п 1пп)+к0(п), где к, п — соответственно количество удаляемых и оставшихся частиц. Эта операция выполняется достаточно редко, и затраты на нее пренебрежимо малы.
Представленный алгоритм применялся для всех типов частиц, включая электроны, В результате описанных выше операций разброс весов частиц относительно среднего уменьшается. Само среднее значение также при изменении параметров плазмы меняется, что позволяет поддерживать количество частиц, а следовательно, трудоемкость и точность расчета на заданных уровнях.
Для сильно электроотрицательных газов значительное укрупнение псевдочастиц может приводить к ухудшению точности расчета в приэлектродном слое, где концентрация ионов существенно меньше их концентрации в квазинейтральной области. Для повышения точности расчета в приэлектродном слое применяется также метод разделения частиц (см., например, [9]),
Использование переменных весов частиц, необходимое для включения в кинетическую модель уравнений химической кинетики, может приводить к увеличению времени расчетов без повышения точности. Склеивание и дробление частиц позволяет поддерживать их количество на заданном уровне,
5. Переменные шаги по времени для частиц одного типа
В данном разделе рассматривается оптимизация алгоритма для задач, в которых важную роль играет вторичная эмиссия электронов с электродов. Простейший способ увеличения производительности — увеличение шага пространственной сетки Ах и шага по времени Аг. Пределом для такого увеличения являются условия стабильности и точности алгоритма. Подробное обсуждение этих ограничений дано в работах [6, 7, 17], Приведем ограничения, характерные для явной схемы расчета движения частиц в методе Р1С-МСС, которую мы используем:
WpeAte < 2,
ADe/Ax > 0.3, vsAts/Ах < 1, Ats < 0.1
(10)
(И) (12)
(13)
где wpe — плазменная частота, Аве _ электронный дебаевекий радиус, vs — характерная скорость частиц типа s, ^coii,max — максимальная частота столкновений. Наилучшее соотношение между точностью и скоростью расчетов достигается вблизи границ неравенств (10)—(13).
Особенностью ВЧЕ-разряда в молекулярных газах является большая частота возбуждения колебательных уровней фонового газа. Эти процессы, как правило, имеют значительные сечения и невысокие пороговые энергии, В результате возбуждения колебательных состояний в молекулярных газах температура электронов гораздо ниже, чем в атомарных. Электроны в области малых полей (в квазинейтральной плазме) термализуютея и не участвуют в процессах ионизации и генерации радикалов [15]. Для таких электронов ограничение на шаг по времени, связанное с условием Куранта vsAts/Ax < 1
Поэтому, если позволяют прочие ограничения (10), (11), (13), можно увеличить шаг по времени для группы термализованных электронов. Данный подход использовался в работе [19] для моделирования тлеющего разряда,
В [20, 21] представлен метод, в котором электроны делятся на несколько групп в соответствии с их энергией таким образом, что каждая следующая группа рассчитывается с вдвое большим шагом по времени, чем предыдущая. Этот подход был разработан для неявных схем расчета, позволяющих преодолевать ограничение, связанное с плазменной частотой wpeAte < 2. Неявные схемы используются, как правило, если есть необходимость расчета больших пространственных областей, при этом они являются более трудоемкими по сравнению с явными схемами при равных шагах по времени. Для ВЧЕ-разряда не характерна большая пространственная протяженность, поэтому
vsAts/Ax < 1
является более жестким по сравнению с условием wpeAte < 2, в силу чего мы попытались применить метод группировки к явной схеме, К сожалению, оказалось, что для рассматриваемых задач вычислительные затраты на логику перегруппировки электронов по энергиям достаточно велики и существенного выигрыша в производительности не получено. Кроме того, применение этого метода уменьшает эффективность распараллеливания, о котором пойдет речь в следующем разделе.
Более эффективным оказался метод, суть которого заключается в пересчете уравнений движения частицы с меньшим шагом по времени, когда нарушается условие Куранта, Если частица проходит за один шаг по времени больше одной ячейки сетки, то производится пересчет ее движения с меньшим шагом по времени Ate,new = Ate/(N +1), где N — количество пройденных частицей ячеек пространственной сетки. На рис, 4 приведены результаты тестовых расчетов разряда в смеси аргона с ацетиленом, где существенную роль играет вторичная эмиссия электронов с электродов. Давление газа P = 70 мТорр, U0 = 200 В, d =2 см, коэффициент вторичной электронной эмиссии Y = 0.2, Состав фонового газа считался постоянным. Расчеты проводились с различными временными шагами для электронов, причем в базовом варианте Ato = 2 х 10-11 с. Видно, что увеличение шага по времени для электронов вдвое Ate = 2At0 мало влияет
а б
Рис. 4. Распределения концентрации электронов (а) и средней энергии электронов (б)
на параметры разряда (отклонение ключевых параметров разряда не превышает 1 %), но при Аte = 4Аt0 ошибка становится существенной (отклонение параметров разряда достигает 20%). При использовании описанной модификации метода удается получить более высокую точность (отклонение ключевых параметров разряда не превышает 5%). При этом время расчета 1000 периодов разряда немодифицированным алгоритмом с Аte = 2At0 составило 10 110 с, а время расчета модифицированным алгоритмом с Аte = 4Ato — 5570 с. Таким образом, получен практически двукратный рост производительности при малых потерях точности. Значительный прирост производительности объясняется тем, что затраты времени на расчет быстрых электронов пренебрежимо малы по сравнению с затратами на расчет термализованных электронов.
Использование различных шагов по времени является привлекательным методом ускорения расчетов, особенно для неявных схем. Для явной схемы расчета более эффективным оказался метод пересчета уравнений движения для быстрых частиц, что позволяет получать более высокую точность с меньшим шагом по времени для задач, в которых вторичная эмиссия электронов с электродов играет существенную роль.
6. Распараллеливание с помощью ОрепМР
Для уменьшения времени расчетов было также проведено распараллеливание алгоритма. Существуют два ключевых метода декомпозиции алгоритма для распараллеливания — на эйлеровой сетке и на лагранжевой сетке. В первом случае расчетная область разбивается на подобласти и каждый процессор осуществляет расчет для одной подобласти. При этом распараллеливаются все части алгоритма, что делает его хорошо масштабируемым. Этот метод особенно эффективен, когда нужно провести расчет для большой пространственной области, особенно в многомерном случае.
В наших задачах вычисления на пространственной сетке при решении уравнения Пуассона и расчете уравнений гидродинамики занимают лишь небольшую часть времени. Это позволяет достигать хороших результатов, используя метод декомпозиции на лагранжевой сетке. При этом параллельно проводится расчет движения псевдочастиц, столкиовительиых процессов, а также пересчет плотности частиц. Затем плотности, рассчитанные разными процессорами, складываются. Масштабируемость этого метода
ограничена долой не распараллеленной части алгоритма. К достоинствам метода можно отнести хорошую балансировку загрузки без применения специальных алгоритмических приемов. Кроме того, следует отметить, что оба метода могут быть использованы совместно.
В случае одномерной модели ВЧЕ-разряда межэлектродпые расстояния невелики и количество ячеек пространственной сетки составляет примерно 100-1000, поэтому декомпозиция на лагранжевой сетке выглядит предпочтительней.
Дня распараллеливания был использован интерфейс ОрепМР, так как on позволяет реализовать распараллеливание с гораздо меньшими затратами по сравнению, например, с MPI. Кроме того, применение ОрепМР не создает проблем переносимости, поскольку программа с использованием ОрепМР может быть скомпилирована на системе, не имеющей ОрепМР.
Недостатком такого подхода является то, что программа применяется только на системах с общей памятью. Это ограничение может оказаться существенным, если необходимо рассчитывать большое количество нсевдочастиц дня повышения точности или в случае большого межэлектродного расстояния. При этом дня большинства рассматриваемых нами задач хорошая точность достигается при количество нсевдочастиц 10 000 100 000, когда производительность перестает существенно расти при использовании 10 процессоров и более. Последнее есть следствие того, что при данном число процессоров затраты времени на не распараллеленную часть алгоритма и на синхронизацию становятся больше затрат на распараллеленную часть алгоритма. Такие же результаты получены в |17| при реализации этого алгоритма с помощью MPI.
Дня сравнения полученных данных с результатами работы |17| были проведены тестовые расчеты дня чистого аргона.
На рис. 5 показан прирост производительности в зависимости от количества процессоров для расчета ВЧЕ-разряда (13.56 МГц) в аргоне. P = 10 мТорр, d = 5 см, U0 = 800 В. В расчетах использовалась равномерная пространственная сетка с количеством ячеек ^gnd = 400. Шаг по времени для электронов Ate = 1/(2000frf), где frf — частота приложенного напряжения. Шаг по времени для ионов At = 20Ate, Время расчета 5000 периодов (характерное время установления) па одном ядре — 170 мин.
Рис. 5. Выигрыш в производительности (С) для разного количества частиц в зависимости от количества процессоров (п); сплошная линия — 40 000 псевдочастиц, штрих — 20 000, пунктир 40 000 [171
Из графиков на рис, 5 видно, что прирост производительности в значительной степени зависит от количества рассчитываемых псевдочастиц. При этом даже для небольшого количества частиц — 20 ООО можно получить пятикратный прирост производительности на восьмиядерпой машине. На рисунке приведена также зависимость прироста производительности, полученная в работе [17] с использованием MPI, К сожалению, в [17] не приведено достаточно данных, чтобы воспроизвести расчеты. Однако мы использовали те же ограничения на шаг сетки и шаг по времени, что и в этой работе, и позволяет сравнивать результаты. Как видно, в нашем случае удалось добиться лучшего прироста производительности, что указывает на большую эффективность ОрепМР по сравнению с MPI,
Заключение
В работе представлены кинетическая и гибридная модели многокомпонентной плазмы высокочастотного емкостного разряда в химически активных газах. Кинетическая модель описывает электроны и все типы ионов с помощью уравнений Больцмана, которые решаются с помощью метода PIC-MCC, Редкие неупругие столкновптельные процессы для ионов, в том числе реакции нейтрализации между положительными и отрицательными ионами, описываются уравнениями химической кинетики. Помимо этого модель включает уравнение Пуассона для расчета электрического потенциала и уравнение диффузии для нейтральных компонентов. Для кинетической модели представлены модификации Р1С-МСС-алгоритма, необходимые для описания многокомпонентной плазмы, такие как пересчет весов частиц при решении уравнений химической кинетики и удаление излишних частиц для поддержания точности и скорости расчетов на заданном уровне.
Кинетическая модель позволяет получить решения ab initio, однако требует больших вычислительных затрат. Для ускорения расчетов были использованы гибридные методы. При увеличении количества рассчитываемых типов ионов затраты на расчет их траекторий движения кинетическим методом занимают большую часть времени расчета, Для уменьшения времени, затрачиваемого на моделирование ионного компонента, использована гибридная модель, в которой электроны описываются кинетически, а ионы — в гидродинамическом приближении. Применение этого метода позволило рассмотреть большое количество различных типов ионов без существенного увеличения расчетного времени по сравнению с разрядом с одним типом ионов.
Включение химических процессов в модель разряда ведет к увеличению времени установления равновесия и, как следствие, к увеличению времени расчетов. При длительном процессе установления химического равновесия функция распределения электронов меняется медленно, что позволяет в определенном интервале времени заменить параметры электронного компонента (плотность, скорости столкповительпых процессов) их средними значениями. Гибридный метод с усреднением по электронному компоненту позволяет сократить время расчета для систем с большим количеством компонентов и химических процессов.
Для ускорения расчета с помощью метода PIC-MCC были применены два подхода — переменных шагов по времени для электронов и распараллеливание, В первом случае при нарушении условия Куранта осуществляется пересчет уравнения движения частицы с меньшим шагом. Этот метод продемонстрировал высокую эффективность в задачах, где вторичная эмиссия электронов с электродов играет существенную роль.
Распараллеливание с декомпозицией на лагранжевой сетке является эффективным методом ускорения расчетов методом PIC-MCC, Использование для такого распараллеливания интерфейса ОрепМР позволяет добиться большего прироста производительности на системах с общей памятью. Тем не менее при расчете больших пространственных областей или при решении многомерных задач более эффективным может оказаться распараллеливание с помощью MPI или применение декомпозиции на эйлеровой сетке за счет того, что количество процессоров на системах с общей памятью, как правило, не очень велико, В дальнейшем планируется реализовать данное распараллеливание с использованием ресурсов современных графических процессоров и библиотеки CUDA, Описанные модели и методы составляют основу комплекса программ для моделирования ВЧЕ-разряда в химически активных смесях, К достоинствам этого комплекса относится возможность применения гибридных методов (в частности, возможность выбора модели для каждого типа ионов в отдельности), а также использование методов ускорения расчета (методы переменных весов частиц и переменных шагов по времени и распараллеливание).
Список литературы
fl] Grill A. Review of the tribologv of diamond-like carbon // Wear. 1993. Vol. 168, No. 1-2. P. 145-153.
[2] Robertson J. Diamond-like amorphous carbon // Materials Sci. and Engineering. R: Reports. 2002. Vol. 37, No. 4-6. P. 129.
[31 Obraztsov A.N., Volkov A.P., Nagovitsyn K.S. et al. CVD growth and field emission properties of nanostructured carbon films //J. Phvs. D: Appl. Phvs. 2002. Vol. 35, No. 4. P. 357.
[4] Aristov V. Direct Methods of Solving the Boltzmann Equation and Study of Nonequilibrium Flows. Dordrecht: Kluwer Acad. Publ., 2001.
[5] Pitchford L., Phelps A. Comparative calculations of electron-swarm properties in N2 at moderate E/N values // Phvs. Rev. A. 1982. Vol. 25, No. 1. P. 540.
[61 Бэдсе л Ч., Ленгдон А. Физика плазмы и численное моделирование. \!.: Энергоатом-издат, 1989.
[7] Хокни Р., Иствуд Д. Численное моделирование методом частиц. М.: Мир, 1987.
[8] kljshner m.J. A model for the discharge kinetics and plasma chemistry during plasma enhanced chemical vapor deposition of amorphous silicon //J. Appl. Phvs. 1988. Vol. 63, No. 8. P. 2532.
[9] Вшивков В.А., Лазарева Г.Г., Снытников А.В. Адаптивное изменение массы модельных частиц при моделировании тлеющего ВЧ-разряда в силановой плазме // Вычисл. технологии. 2008. Т. 13, № 1. С. 22-30.
[10] Schweigert I., Alexandrov A., Ariskin D. et al. Effect of transport of growing nanopar-ticles on capacitivelv coupled rf discharge dynamics // Phvs. Rev. E. 2008. Vol. 78, No. 2. P. 026410.
[11] Belenguer P., Boeuf J. Transition between different regimes of rf glow discharges // Phvs. Rev. A. 1990. Vol. 41, No. 8. P. 4447.
[12] Nienhuis G.J., Goedheer W.J., Hamers E.A.G. et al. A self-consistent fluid model for radio-frequency discharges in SiH4-H2 compared to experiments //J. Appl. Phvs. 1997. Vol. 82, No. 5. P. 2060.
[13] Schweigert I.v., Schweigert V.A. Combined PIOMCC approach for fast simulation of a radio frequency discharge at a low gas pressure // Plasma Sources Sci. and Technology. 2004. Vol. 13, No. 2." P. 315.
[14] Ariskin D.A., Schweigert I.V., Alexandrov A.L. et al. Modeling of chemical processes in the low pressure capacitive radio frequency discharges in a mixture of Ar/C2H2 //J- Appl. Phvs. 2009. Vol. 105, No. 6. P. 063305.
[15] Ariskin D.A., Schweigert I.V. Influence of an acetylene impurity on the properties of a radio-frequency gas discharge in argon //J. Experimental and Theoretical Phvs. 2009. Vol. 109, No. 4. P. 707.
[16] Brackbill J.U., Cohen B.I. Multiple Time Scales. Acad. Press, 1985.
[17] Kawamura E., Birdsall C.K., Vahedi V. Physical and numerical methods of speeding up particle codes and paralleling as applied to RF discharges // Plasma Sources Sci. and Technology. 2000. Vol. 9, No. 3. P. 413.
[18] Sommerer T.J., kljshner M.J. Numerical investigation of the kinetics and chemistry of rf glow discharge plasmas sustained in He, N2, 02, He/N2/02, He/CF4/02, and SiH4/NH3 using a Monte Carlo-fluid hybrid model // J. Appl. Phvs. 1992. Vol. 71, No. 4. P. 1654.
[19] Швейгерт В.А., Швейгерт И.В. Математическое моделирование прикатодной области стационарного тлеющего самостоятельного разряда // ПМТФ. 1988. № 4. С. 16-23.
[20] Friedman A., Parker S., Ray S., Birdsall С. Multi-scale particle-in-cell plasma simulation // J. Comput. Phvs. 1991. Vol. 96, No. 1. P. 54.
[21] Parker S.E., Friedman A., Ray S.L., Birdsall C.K. Bounded multi-scale plasma simulation: application to sheath problems //J. Comput. Phvs. 1993. Vol. 107, No. 2. P. 388-402.
Поступила в редакцию 2 июня 2010 г., с доработки — 8 ноября 2010 г.
