Новый подход к построению и анализу алгоритмов метода частиц для уравнения Больцмана

Более чем столетняя история кинетического уравнения Больцмана [9], [126] подобна величественному айсбергу. Вершину его украшают блестящие математически строгие результаты по существованию, единственности, асимптотическим свойствам решения, начиная с Д. Гильберта [148], С. Чэпмена [35], Т. Карлемана [18] и Г. Греда [140] и вплоть до ярких результатов последних лет, принадлежащих А. А. Арсеньеву [3], А. В. Бобылеву [8], Н. Б. Масловой [40], Р. Кэфлишу [132], О. Э. Ланфорду [87], К. Бардошу и С. Укаи [118], Р. Ди Перне и П. Лионсу [135]. Другая, не столь известная, но наиболее значительная, часть этого айсберга, олицетворяет собой огромный труд сотен исследователей по численному и аналитическому решению уравнения Больцмана и его обобщений во многочисленных приложениях. Именно эти приложения, число которых постоянно растет, и обеспечивают фундаментальность, значимость и непотопляемость в мировом океане науки такого неординарного явления как уравнение Больцмана.

Хорошо известно, что многолетний опыт по аналитическому и асимптотическому исследованию уравнения Больцмана обобщен в классических книгах С. Чэпмена и Г. Каулинга [35], Дж. Ферцигера и Г. Капера [33], М. Н. Когана [20], К. Черчиньяни [36] и современных монографиях [3], [40], [39]. Что же касается опыта и методов численного решения уравнения Больцмана, построения вычислительных алгоритмов на основе строгих с математической точки зрения следствий уравнения и его аналогов, исследования свойств сходимости, то эта проблематика не так полно отражена в литературе и мы не можем назвать какой-нибудь одной монографии, в которой с единой точки зрения достаточно полно и строго рассматривались бы наиболее известные подходы к численному решению кинетического уравнения Больцмана, от регулярных сеточных методов до алгоритмов метода частиц и алгоритмов прямого статистического моделирования. По всей видимости, единственной монографией, посвященной регулярным и полурегулярным методам численного решения кинетического уравнения Больцмана, сочетающим сеточную аппроксимацию дифференциального оператора переноса со специальной техникой вычисления многомерных интегралов столкновений, является монография В. В. Аристова и Ф. Г. Черемисина [2]. В настоящее время этот подход активно развивается также в работах Ж. Шнайдера и Ф. Роже [183] с использованием фундаментального результата [179], A.B. Бобылева и С. Рязанова [123] и других.

Более подробно отражен в литературе подход к численному решению уравнения Больцмана, связанный с использованием так называемого метода частиц. Метод частиц является одним из основных алгоритмов математического моделирования кинетических процессов. Существует две точки зрения на метод макрочастиц в кинетической теории. Согласно первой точке зрения метод макрочастиц — это метод имитационного моделирования процессов на микроуровне [6, 5, 16, 34]- согласно второй точке зрения — это метод численного решения кинетических уравнений [3, 178]. Кинетические уравнения являются математическим аппаратом кинетической теории, описывающей эволюцию макроскопических систем, состоящей из большого «10(10) — -10(26)) числа частиц. При классическом подходе основным объектом в кинетической теории является од-ночастичная функция распределения f (x, v, t) фазовых координат {ж-г>}, где V — скорость и времени t. Утверждается, что функция распределения / при определенных условиях, в частности при Л^физ. —> оо, удовлетворяет кинетическому уравнению.

— Cif).

Оператор С в типичных приложениях представляет собой нелинейный интегро-дифференциальный оператор. Существенная особенность кинетических уравнений состоит в том, что их решения в любой фиксированный момент времени можно интерпретировать как плотность некоторой меры (меры Радона) ?it{ da- (g) dv) относительно меры Лебега в фазовом пространстве.

Лг (с1сс<8> с1г>) = /(ж, г", с! зг с! г".

С практической точки зрения интересны имеющие физический смысл моменты этой меры по переменной V (макропараметры).

Мф (х, ь) = а| ^(у)$?у, ж3 связанных с молекулярными признаками ф{у). Метод частиц как имитационный метод прямого статистического моделирования (метод Монте-Карло) имеет богатую историю. По всей видимости, впервые он был предложен Н. Метрополисом и С. Уламом в 1949 году [171]. В дальнейшем метод прямого статистического моделирования успешно использовался для решения важнейших прикладных задач, в частности, задачи переноса нейтронов в ядерных реакторах [15, 32], переноса носителей заряда в полупроводниках и полупроводниковых приборах [17, 34] и других задач, приводящих к кинетическим уравнениям с линейными интегралами столкновений. Были созданы теоретические основы метода статистического моделирования для решения широкого класса задач, разработаны эффективные способы построения схем метода частиц решения линейных интегральных уравнений, основанный на использовании ряда Неймана и соответствующей итерационной процедуры [14]. Позднее были разработаны методы частиц решения нелинейных уравнений газовой динамики [5, 144] и уравнения Власова [29, 3].

Метод частиц постепенно превратился в «промышленный» алгоритм решения широкого класса задач. К достоинствам этого подхода следует отнести его относительную «нечувствительность» к размерности задачи, ясную интерпретируемость и физическую наглядность, простоту программирования, возможность эффективного распараллеливания. Все это обусловило стремление исследователей использовать преимущества метода частиц при моделировании сред, описываемых на основе кинетического уравнения Больцмана.

Основной трудностью при разработке алгоритмов прямого статистического моделирования течений разреженного газа являлась эффективная организация учета парных столкновений. Здесь под эффективным учетом парных столкновений мы понимаем такую процедуру, которая в единицу времени требовала бы числа арифметических действий, пропорционального числу макрочастиц. С точки зрения имитационного эвристического подхода в общем случае разработать такую процедуру оказалось достаточно сложно, поскольку она не отвечала естественным комбинаторным представлениям о числе пар. Этим объясняется развитие в шестидесятые годы методов типа метода пробных частиц [147, 68], основанных на приближенной итерационной процедуре для нелинейного уравнения Больц-мана и позволяющая в определенном смысле избежать перебора пар сталкивающихся частиц. К сожалению, не приводили к эффективным вычислительным алгоритмам и методы, основанные на точном представлении решения кинетического уравнения Больцмана с помощью ветвящихся случайных процессов специального вида [15], так как требовали значительного числа арифметических действий. Следует отметить, что в последнее время достигнут заметный прогресс в оптимизации соответствующих вычислительных схем [72].

В 1963;1970 гг. был предложен первый алгоритм прямого статистического моделирования течений разреженного газа, известный как схема Дж. Берда [6], основанная на методе расщепления по физическим процессам и позволяющая эффективно учитывать влияние парных столкновений. Метод основывался на эвристических соображениях и его связь с уравнением Больцмана была не ясна. Тем не менее подход Берда явился этапом в развитии алгоритмов метода частиц, с его помощью был решен ряд принципиальных задач динамики разреженного газа, имеющих большое практическое значение.

Метод прямого статистического моделирования постоянно развивался и совершенствовался как самим Дж. Бердом, так и другими авторами. При этом особое внимание уделялось формальному вероятностному обоснованию методов типа Берда и их анализу с позиций общей теории методов Монте-Карло [14, 31], что позволило не только обратиться к известным способам повышения их эффективности, но и создать принципиально новые подходы к прямому статистическому моделирования течений разреженного газа. Здесь прежде всего следует назвать фундаментальную работу О. М. Белоцерковского и В. Е. Яницкого [55, 56], положившую начало моделированию течений разреженного газа на основе схемы испытаний.

Бернулли. В [54, 113] был проведен теоретический анализ схемы Берда и схемы испытаний Бернулли, построен аналог ТУ-частичного основного кинетического уравнения.

В динамике разреженного газа при исследовании методов прямого статистического моделирования с позиций общей теории методов Монте-Карло принято основываться не на уравнении Больцмана, а на линейном уравнении Леонтовича-Каца [19, 88] относительно^{-частичной} функции распределения. Это означает, что в качестве базовой математической модели рассматривается не уравнение Больцмана, а некоторый случайный процесс, уравнением Колмогорова для которого является Л/^-частичное уравнение. Уравнение Леонтовича сыграло большую роль в развитии и анализе методов частиц моделирования течений разреженного газа, на основе уравнения Леонтовича-Каца с использованием общей теории методов Монте-Карло М. С. Ивановым и С. В. Рогазинским в 1988;1990гг. был разработан один из наиболее эффективных алгоритм этого класса — метод мажорантной частоты [16], так же, как и метод Берда, предполагающий построение всей траектории соответствующего случайного процеса.

Вместе с тем, связь уравнения Леонтовича-Каца с уравнением Больцмана установлена в простейших ситуациях на физическом уровне строгости. Решающим предположением при этом является требование «молекулярного хаоса» [20, 26]. Это условие не выполнено при прямом статистическом моделировании посредством конечного числа частиц [16]. Отметим также, что не существует аналога уравнения Леонтовича-Каца как математической модели для исследования классических кинетических процессов в полупроводниковой плазме, в то время как моделирование методом частиц играет играет при изучении полупроводников и полупроводниковых приборов ведущую роль [34].

Кинетическое уравнение Больцмана служит основой математических моделей в различных приложениях. В то же время ведущим алгоритмом математического моделирования в этих приложениях является метод частиц, который в каждом конкретном случае развивался самостоятельно.

55 55 и независимо от других приложении как правило на основе своих эвристических представлений. Это объясняет актуальность и целесообразность разработки алгоритмов построения схем метода частиц решения непосредственно уравнения Больцмана, справедливых для большинства приложений, причем доказательство сходимости таких методов не должно базироваться на гипотезе «молекулярного хаоса» .

Настоящим прорывом в создании эффективных алгоритмов решения уравнения Больцмана и моделирования течений разреженного газа ознаменовалась вторая половина восьмидесятых годов. Этот прорыв был обусловлен появлением работ К. Нанбу [176], предложившим новый подход к построению схем метода частиц непосредственно из уравнения Больцмана, работ А. В. Скорохода [30] и Х. Танаки [192], о приближении к решению уравнения Больцмана для псевдомаксвелловских молекул, фундаментальной работы А. Шнитмана [190] о свойствах пространственно-однородного уравнения Больцмана в случае произвольного закона межмолекулярного взаимодействия. Особое значение в этом ряду имеет работа А. А. Арсеньева [45, 46], где он впервые высказал идею о возможности построения алгоритмов стохастического моделирования решения уравнения Больцмана как алгоритмов, реализующих систему стохастических дифференциальных уравнений по случайной мере Пуассона.

Основной целью диссертации является разработка новых способов построения алгоритмов метода частиц как численного метода математической теории кинетического уравнения Больцмана. В качестве приложений рассматриваются динамика многокомпонентного разреженного газа и перенос заряда в многодолинных полупроводниках.

Предлагаемый в диссертации подход к построению алгоритмов метода частиц восходит к идеям А. В. Скорохода [30] и А. А. Арсеньева [45, 46, 47] о приближении обобщенного решения уравнения Больцмана решением некоторой системы стохастических уравнений по мере Пуассона и к идее А. А. Арсеньева о возможности использования таких систем для построения алгоритмов типа схемы Берда, при этом речь идет уже не о эвристической схеме прямого моделирования, а о численных алгоритмах решения непосредственно нелинейного уравнения Больцмана. Таким образом разрабатываемый в диссертации подход заключается в следующем. Кинетическое уравнение рассматривается как уравнение для плотности некоторой меры. Процесс построения метода частиц разбивается на два этапа. Во-первых, следует указать метод построения разностного аналога кинетического уравнения. В качестве разностного аналога уравнения выбирается система разностных стохастических уравнений по мере Пуассона на случайной временной сетке относительно узлов кубатурной формулы, аппроксимирующих меру — решение исходного уравнения. При этом разностным параметром является величина где N — число узлов кубатурной формулы, интерпретируемых как макрочастицы. Во-вторых следует указать эффективный метод решения (точного или приближенного) аппроксимирующей системы стохастических уравнений. При этом различаются потраекторное сильное решение и слабое решение в смысле распределений [13]. Напомним, что под эффективным методом решения мы понимаем численный алгоритм, требующей для своей реализации О (А/") арифметических операций. Такая точка зрения отличается от подхода «прямого статистического моделирования» и более близка классическому подходу теории разностным схем [27]. Она позволяет также по-другому взглянуть на выбор числа частиц и поставить вопрос об оптимизации вычислительного процесса за счет выбора аппроксимирующей системы и численного метода ее реализации. В качестве способа построения аппроксимирующей системы стохастических уравнений для узлов в диссертации предлагается метод ветвей функции рассеяния.

Проиллюстрируем основную идею нашего подхода на примере классического уравнения Больцмана динамики разреженного газа в пространственно — однородной ситуации. Первым шагом к построению метода частиц решения кинетического уравнения является переход к обобщенной (для мер) форме этого уравнения. Рассмотрим уравнение Больцмана в обобщенной форме (безразмерные переменные) t.

Ф, к) = (ф, мо> + ?({А[н, у. а)) а з, (0.1) о где оператор, А: й (К3) С6(Ж6).

5(2).

Здесь сг (с1£) — стандартная инвариантная относительно вращений мера на единичной сфере.

Функция скачка ф (у, у 1,?), п Е у, у Е И3, характеризует скачкообразное изменение скоростей частиц у и у при парном столкновении, отвечающем геометрическому параметру п, причем предполагаются выполненными «законы сохранения», имеющие в нашем случае вид.

— ф 0,^1,0 = -ф (у, уи?), уъу е И3, п е 5(2) (закон сохранения импульса при элементарном акте взаимодействия) и.

У1-У + ф (^1,0 (^,^1,0) = 0 закон сохранения энергии при элементарном акте взаимодействия). Непрерывная неотрицательная функция ?) зависит от потенциала взаимодействия частиц и удовлетворяет условиям типа.

И (у + ,-п) = г (у, уи?) = Н (уиу,?), уиу в И3, п? 5(2).

Функция /г (г>, г>1,?) характеризует тенденцию взаимодействия частиц со скоростями у и г>1, находящихся в малой окрестности точки х и чье взаимное расположение задается посредством вектора п. Поэтому функцию /г,(г>, г>1, ?) называют также функцией рассеяния.

Наибольшую трудность при численном решении уравнения Больцмана представляет собой интеграл столкновений. При решении прямыми методами с использованием квадратурных формул основная сложность заключается в необходимости использования интерполяции в скоростном пространстве, обусловленной функцией скачка ф. В свою очередь, одной из основных трудностей при разработке эффективных алгоритмов метода частиц решения уравнения Больцмана является зависимость функции рассеяния к от модуля относительной скорости у — что имеет место практически для всех моделей взаимодействия. В случае, когда это не так (псевдомаксвелловские молекулы), эффективные алгоритмы статистического моделирования течения разреженного газа были известны ранее, хотя их строгое обоснование появилось лишь в последнее время. Напомним, что под эффективным алгоритмом метода частиц мы понимаем алгоритм, требующий для реализации арифметических действий, число которых пропорционально количеству частиц.

Для преодоления этой трудности в диссертации предлагается представление интеграла столкновений в специальном каноническом виде. Возможность такого представления связана с очевидными тождествами 1 = }%(а>17) ?77, о справедливого для любых 0 < а < 1, и у (у + % • Ф) — у (у) = [у (х + г) — у (ж)] • % справедливого для любых функций у (у), любых скачков Ф и любых индикаторов т. е. функций, принимающих только два значения, = Ц-Предположим для простоты, что т (О > 0.

Введем в рассмотрение функцию-индикатор

Ь (У, уъ€) < Нтахуо (?) К.

Нтахи) (?) а и соответствующую модифицированную функцию скачка.

Ф (^ь?, а) = | (у, уи^а)ф (У, Уь?) и новое пространство «геометрических» параметров столкновений с мерой в = [0,1], 0 = {?, а}е&, т (<10) = Нтаяг"{?)(т (<1?) ¿-а.

Предлагаемый в диссертации алгоритм построения схем метода частиц состоит из двух основных этапов. На первом этапе за счет специального выбора модифицированной функции скачка Ф и нового пространства геометрических параметров Хг с конечной мерой у (с1£) исходное уравнение для мер преобразуется к виду (0.1) с оператором, А в канонической форме.

А[Н, ф]{у, у1) = Л[ф](у, у1) =.

I [ф (у + Ъ (у, уи?))-ф (у)]1((10.

Отметим, что в каноническом представлении выбор пространства XI с конечной мерой не зависит от значений г-, г>1, и определяется типом аппроксимации функции /г, а модифицированная функция скачка является функционалом от функции рассеяния, т. е. г>1, ?) =.

Рз1[К (у, у 1,?). Пространство Х% с конечной мерой 7(а£) мы в дальнейшем будем называть геометрическим шаблоном (или просто шаблоном) динамической схемы стохастически взаимодействующих частиц для уравнения Больцмана. Наряду с выбором фазового пространства частиц выбор геометрического шаблона и модифицированных функций скачка в каноническом представлении фактически определяет динамическую схему стохастически взаимодействующих частиц, которая выбирается в качестве разностного аналога уравнения Больцмана как система стохастических дифференциальных уравнений по мере Пуассона, задающих изменение фазовых координат частиц во времени (уравнения движения). В случае пространственно-однородного уравнения Больцмана динамическая схема с пространством скоростей в качестве фазового пространства имеет вид: уг {t) = vi (0) + У£l<i<j<N г.

УI (?в х а.

О 0 где символами р^, 1 < I < з < N обозначено семейство независимых в совокупности случайных мер Пуассона на каноническом шаблоне с мерой интенсивности вх (И) = -^т (с!0) аг,.

Рч = № 1 < з < М.

На основании формулы Ито и оценок решения нелинейного уравнения Больцмана доказывается, что при выполнении стандартных условий на функции рассеяния и скачка эмпирические меры.

А = ~ Е «М»)} 1<�г<�Ы в смысле слабой сходимости (для мер) сходятся по вероятности к неслучайной мере, являющейся решением уравнения Больцмана :

Дг /Д N -«¦ оо, 0, если только в начальный момент времени.

Р) /Д ЛГ оо.

В диссертации динамические схемы системы стохастически взаимодействующих частиц построены и обоснованы:

• Для системы уравнения Больцмана (глава 2);

• Для пространственно-неоднородного уравнения Больцмана (глава 3);

• Для уравнения Больцмана динамики многокомпонентного разреженного газа с учетом физико-химических процессов (главы 2 и 4).

• Для уравнения Больцмана теории полупроводников (глава 5).

В главе 4 разрабатываемый в диссертации подход применяется для моделирования течений многокомпонентного разреженного газа с учетом физико-химических процессов. Как известно, математические модели различных процессов динамики многокомпонентного разреженного газа при определенных условиях могут быть описаны посредством систем интегро-дифференциальных уравнений относительно одночастичных функций распределения компонент. С физической точки зрения эти задачи характеризуются широкими пределами типичных значений макропараметров, наличием значительного числа кинетических процессов, протекающих с различными скоростями. С математической точки зрения речь идет о решении систем нелинейных интегро-дифференциальных уравнений с интегралами столкновений сложного вида, разнообразием характерных параметров, что затрудняет использование классических аналитических, численно-асимптотических и разностных подходов. Регулярные и полурегулярные методы, основанные на разностной аппроксимации кинетических уравнений, детерминированных или стохастических квадратурных формулах для вычисления интегралов, а также разложении неизвестных функций в ряд по малому параметру могут быть эффективно использованы в относительно простых ситуациях [24, 39, 35, 25]. В [1] рассмотрены способы решения линеаризованных систем для смеси газов, находящейся вблизи детального химического равновесия.

Известно, что для решения задач с физико-химическими превращениями хорошо зарекомендовали себя алгоритмы метода частиц. Как и в случае простого газа в рамках подхода прямого статистического моделирования этот метод состоит в конструировании, исходя из физической картины явления, такого случайного процесса, что оценки математического ожидания некоторых случайных функций соответствовали в определенном смысле характеристикам реальной системы. Одним из наиболее широко используемых вариантов метода частиц для моделирования смесей газов является процедура Берда [194, 6], позволяющая получить физически приемлемые результаты при достаточно большом числе модельных частиц [187].

Другие схемы статистического моделирования основаны на использовании случайного процесса Пуассона и уравнений типа Леонтовича-Каца[55, 78, 103, 76, 100, 77, 109] (подробный обзор представлен в [81]).

Особенности моделирования методом частиц поведения многокомпонентной смеси с физико-химическими взаимодействиями приводят к идее использования весовых алгоритмов, основные черты которых были изложены в работах [6, 14]. Впервые обоснованный подход к ведению зависящих от времени случайных весовых множителей (дополнительной фазовой характеристики частицы) для этого класса задач был предложен в работе [85]. В приложениях к задачам физико-химической кинетики метод частиц прямого статистического моделирования состоит в следующем: многокомпонентный реагирующий газ, состоящий из s компонент, моделируется системой частиц, координаты и скорости которых в начальный момент времени распределены в соответствии с заданными функциями распределения компонент. Частицы взаимодействуют друг с другом и перемещаются по расчетной области в соответствии с некоторыми принятыми соглашениями. Эти соглашения-модели различаются прежде всего по способам реализации пространственных неоднородностей рассматриваемых течений. Кроме традиционного метода «частиц в ячейках» [5, 6, 55], предполагающего постоянство свойств газа внутри пространственной ячейки и только в ней взаимное рассеяние частиц, в диссертации предлагаются полевые и структурно-полевые алгоритмы, см. также [80], не требующие введения расчетной пространственной сетки и допускающие взаимодействия частиц из поля течения или из структуры с учетом дополнительного параметра столкновения — столкновительной зоны.

Важное различие между алгоритмами статистического моделирования состоит в способе организации физико-химических взаимодействий в смеси газов. Первая группа методов при моделировании каждого компонента смеси использует число частиц, пропорциональное его концентрации [76, 156]. При неупругом акте частицы переходят в новое состояние с выполнением всех законов сохранения. Общее число частиц в сорте меняется, во всей системе остается постоянным. Использование подобных моделей (пропорционального представительства) оправдано при расчете смесей с близкими концентрациями компонент в течении всего времени решения или при достаточных ресурсах ЭВМ. С другой стороны ясно, что зависимость числа моделирующих частиц от концентрации не позволяет эффективно с точки зрения вычислительных затрат моделировать смеси с большими перепадами концентраций компонентов [111]. Эту проблему можно решить с помощью введения дополнительной характеристики частицывесового множителя. В [62, 110] веса частиц одного сорта полагаются равными и в начальный момент времени, и вдоль всей расчетной траектории процесса/ число моделирующих частиц не зависит от концентраций. При неупругом взаимодействии частица с некоторой вероятностью остается в своем сорте или переходит в другой, меняя при этом и свой вес. Законы сохранения выполняются в среднем. Число частиц в системе и в каждом компоненте не сохраняется. Достоинства схемы" независимость числа частиц от концентрации сорта и отсутствие дисперсий весов для частиц одного сорта.

Независимость числа частиц от концентраций была достигнута в работе [85] с помощью введения для каждой частицы ее индивидуального весового множителя. Законы сохранения выполняются в среднем. Число частиц в сорте во время релаксации остается постоянным, при неупругом взаимодействии происходит перераспределение между частицами определенной части их веса. При использовании подобных алгоритмов, однако, появляется особенность, состоящая в том, что дисперсии оценок вычисляемых величин могут расти и достигать неприемлемо больших значений, что существенно влияет на точность получаемых результатов и на время решения задачи на ЭВМ [111]. В [97] разработана схема моделирования неравновесных химических реакций, имеющая меньшие по сравнению с [85] дисперсии оценок макропараметров, и являющаяся модификацией алгоритма из [85]. В [70] подходы [55, 85, 97] обобщены на пространственно-неоднородный случай, построен весовой алгоритм, соответствующий интегро-конечноразностной схеме, аппроксимирующей уравнение Больцмана с точностью до величины первого порядка малости по характерному пространственному размеру и временному интервалу расщепления процессов при предположении статистической независимости частиц.

В отдельную ветвь можно выделить методы структурного стохастического моделирования [99, 102], основанные на детализации кинетических моментов релаксационных процессов и успешно применяемых для расчетов некоторых течений разреженных многокомпонентных газовых сред.

Разработка и обоснование новых эффективных алгоритмов метода частиц решения системы уравнения Больцмана на основе предлагаемого в диссертации подхода является целью второй главы диссертации.

Весовые алгоритмы метода частиц также рассматриваются в диссер-тационой работе как численные алгоритмы реализации динамических схем. При таком подходе проблема выбора весового алгоритма сводится к построению соответствующей аппроксимации уравнения Больцмана его разностным аналогом — системой стохастических уравнений. За счет удачного выбора аппроксимации, которую, вообще говоря, трудно угадать исходя только из физических предпосылок и здравого смысла, оказалось возможным при фиксированном числе частиц и интервале времени достигать любого наперед заданного уровня дисперсии плотности.

В рассматриваемых стохастических алгоритмах получила дальнейшее развитие идея из [97], состоящая в том, что изменение весов частиц при неупругих взаимодействиях оптимальнее (с точки зрения эффективности и точности вычислений) производить почти наверное. В отличие от [97] рассмотренные алгоритмы позволяют также изменять долю «отщепляемого» веса частиц в зависимости от характера решаемой задачи и требуемой точности получаемых результатов. При конечном времени расчета дисперсию плотности компонентов смеси можно сделать как угодно малой. Соответствующие теоремы сходимости приведены в главе 2. В главе 3 алгоритмы метода частиц для системы уравнений Больцмана обобщены на пространственно-неоднородный случай.

Практическому применению стохастических алгоритмов моделирования течений многокомпонентного разреженного газа с учетом физико-химических процессов посвящена глава 4, где также рассматривается вопрос об учете внутренних степеней свободы молекул в рамках модели Ларсена-Боргнакке. Стохастические алгоритмы тестировались на примере задачи о структуре плоской стационарной ударной волны в одноком-понентном газе, смеси нейтральных и реагирующих газов. Метод частиц для решения данной задачи в случае нейтральной смеси использовался в [120, 186, 70, 86]. В [70, 53] разработана методика моделирования ударной волны с учетом химических реакций. Подробный анализ данных подходов изложен в [69].

В главе 4 приведены также решения задач пространственно-однородной и неоднородной релаксации для смеси газов, о сохранении максвеллов-ского распределения, позволившие провести численное сравнение эффективности и точности различных способов аппроксимации сечения рассеянияпроведено сравнение с другими алгоритмами метода частиц: схемами [6, 55, 78, 170]. Результаты численных расчетов экспериментально подтвердили эффективность предложенных алгоритмов.

Известно, что решение задач обтекания тел разреженным газом в двух-и трехмерном приближении с учетом физико-химических процессов является достаточно сложной проблемой. Такие задачи решались автором на основе нового подхода в рамках общеевропейского проекта «Гермес» [168], некоторые результаты этой работы приведены в диссертации, в частности приведены результаты расчета двумерной задачи обтекания эллипса разреженным газом. В случае реагирующей смеси в несколько раз возрастает трудоемкость расчетов на ЭВМ, что значительно «удорожает» систематическое исследование влияния неупругих взаимодействий на структуру возникающих течений /с^еваск.

Для иллюстрации возможностей развиваемого подхода для решения указанного выше класса задач в диссертации рассматривается модель [122]. В [122] показано, что в случае осесимметричного затупленного тела при наличии некоторых дополнительных (возможно, эмпирических) данных.

173, 174] достаточно рассчитать течение вдоль оси симметрии потока. В работе [122] решена задача расчета одномерного течения, моделирующего обтекание тормозного экрана межорбитального транспортного аппарата. Исследовано влияние числа компонент в потоке и химических реакций на характеристики обтекания. Получающиеся при этом результаты хорошо согласуются с экспериментом и не требуют чрезмерных затрат ресурсов ЭВМ. Пример решения такой задачи приведен в главе 4. В качестве исходной рассмотрена задача обтекания эллиптического цилиндра под нулевым углом атаки с учетом вращательных степеней свободы молекул (модель Ларсена-Боргнакке). Учет же бимолекулярных химических реакций проводился в рамках постановки Берда. Это позволило на качественном уровне проиллюстрировать влияние физико-химических превращений на структуру течения.

Важной областью приложения предлагаемого в диссертации подхода является теория полупроводников. Рассмотрению проблем моделирования процессов переноса заряда в многодолинном полупроводнике посвящена глава 5.

Разработка эффективных численных алгоритмов решения систем кинетических уравнений теории полупроводников является сегодня одной из актуальных проблем математического моделирования в физике. Известно, что математические модели процессов переноса заряда в полупроводниках при определенных условиях могут быть описаны системами интегро-дифференциальных уравнений относительно одночастичных функций распределения носителей заряда, принадлежащих к различным долинам или подзонам зоны проводимости (или валентной зоны) полупроводника. Обращение к кинетическому уровню описания позволяет корректно моделировать процессы переноса в субмикронных элементах интегральной электроники, где в силу сильной неравновесности электронной плазмы неприменимы (или ограниченно применимы) модели, основанные на уравнениях для макроскопических характеристик [43, 104, 105].

На практике значительное распространение получили методы частиц (методы Монте-Карло) численного моделирования кинетических явлений в полупр ов одниках.

Как правило, такие методы строятся как методы имитационного моделирования процессов на микроуровне [49, 17, 34, 188]. Другой подход к построению методов частиц основан на приближении обобщенного (слабого) решения кинетического уравнения и использует предлагаемый в диссертации метод динамических схем. стохастически взаимодействующих частиц.

Метод частиц (как метод имитационного моделирования решения кинетического уравнения) заключается в моделировании движения электронов в координатном и импульсном пространстве, в процессе которого частицы движутся под действием сил, учитываемых уравнением, и рассеиваются (перебрасываются в пространстве импульсов и переходят из долины в долину) в соответствии с вероятностями рассеяния, входящими в кинетическое уравнение, с помощью процедур Монте-Карло [22, 137]. Различают одночастичные и многочастичные методы Монте-Карло. Одночас-тичные методы применяются в основном для решения линейных стационарных задач и заключаются в моделировании движения одного электрона достаточно длительное время. Для вычисления средних используется утверждение, согласно которому средние по времени и ансамблю совпадают. Многочастичные методы позволяют рассматривать переходные процессы и учитывать многочастичные эффекты (например, электрон-электронное рассеяние), в них используется усреднение по ансамблю. Известно [160], что впервые метод частиц для моделирования явлений переноса в полупроводниках применил Куросава. Детальное описание методов частиц моделирования решения кинетического уравнения приведено в обзорах [49, 17, 154, 166, 43]. Известно их обоснование на физическом уровне строгости [137, 182].

В настоящее время накоплен богатый опыт расчета свойств полупроводников и полупроводниковых приборов в рамках кинетического подхода [50, 22, 129, 137, 138, 166]. Точность и достоверность результатов численных экспериментов определяется как методом решения кинетического уравнения, так и степенью соответствия математической модели реальности. При этом большое значение имеет полнота и корректность учета в математической модели различных механизмов рассеяния электронов в полупроводниках. К основным механизмам рассеяния носителей заряда обычно относят рассеяние на деформационно-акустическом потенциале, на полярных оптических и междолинных фононах и рассеяние на ионизированных примесях. Межэлектронное взаимодействие (е-е рассеяние), наряду с другими нелинейными механизмами (электрон-дырка, дырка-дырка), обычно не учитывается в математической модели [17, 154]. Это связано в основном с алгоритмическими и вычислительными трудностями, возникающими при попытках учета е — е взаимодействия в рамках стандартного подхода методов частиц, и является следствием принципиальной нелинейности соответствующих интегралов столкновений.

Кроме того, среди специалистов по вычислительной физике известно мнение, что парные взаимодействия электронов не должны оказывать существенного влияния на основные транспортные характеристики электронов в полупроводниках [145, 146, 152]. Это заблуждение основано, по всей видимости, на том факте, что, в силу известного свойства сохранения, интеграл парных взаимодействий не дает непосредственного вклада в макроскопическое уравнение для средней дрейфовой скорости в однодо-линном приближении. В многодолинном полупроводнике с неэквивалентными долинами дрейфовая скорость не является инвариантом. В любом случае в силу нелинейности задачи и сложного характера взаимодействия линейного и парного рассеяний, влияние последнего может быть весьма значительным.

В последнее время ряд авторов указывал на необходимость разработки эффективных методов учета электрон-электронного рассеяния (в различных приближениях) при численном исследовании свойств горячих электронов [130, 142, 149, 153, 167, 169].

Первые алгоритмы учета парных взаимодействий носителей заряда в рамках кинетического подхода опирались на упрощающие предположения относительно предполагаемого вида функции распределения электронов, позволяющие линеаризовать интеграл е —е взаимодействия. Так, в 1971 г. предложен метод итераций [117], который состоял из следующих шагов. Первый шаг совпадал со стандартным методом статистического моделирования [154] (без учета е — е рассеяния). На втором шаге осуществлялся ряд последовательных итераций алгоритма Монте-Карло с учетом е — е рассеяния. При этом вероятность электрон-электронного рассеяния считалась постоянной, а функция распределения электронов — смещенной максвелловской, соответствующей дрейфовой скорости и средней энергии, полученным на предыдущем шаге итераций. Процедура заканчивалась по достижении стационарного решения. Похожие алгоритмы учета парных столкновений, носящие модельный характер и отражающие общую тенденцию к максвеллизации функции распределения в результате пар ных столкновений, использовались в работах [134, 152, 172, 177, 191].

Подход, не опирающийся на какие-либо априорные предположения относительно формы функции распределения электронов, был предложен в работе [169]. Авторы рассматривали движение двух модельных частиц, каждая из которых представляла часть системы электронов. Для получения значения сечения е — е рассеяния для одной из частиц использовалась «память» другой частицы, хранившая информацию о ее состоянии в предыдущие моменты времени (что позволяло оценить функцию распределения). В случае объявления парного взаимодействия, «память» рассматриваемой в данный момент частицы обновлялась, в то время как «память» партнера оставалась неизменной. Затем частицы менялись местами, и теперь уже первая частица моделировала окружающий вторую частицу электронный газ. Процедура повторялась рекуррентно до достижения стационарного решения.

Вышеизложенные методы, построенные исходя из физической интуиции исследователей, пригодны для решения только стационарных задач. Метод, применимый для исследования переходных процессов, описан в [167, 184]. Авторы использовали многочастичный алгоритм метода статистического моделирования. Частота электрон-электронного рассеяния вычислялась по эмпирической формуле.

Попытки построить более экономичные (в смысле затрат машинного вре.

V* о мени) алгоритмы учета парных взаимодеиствии предпринимались и ранее. Так, в [143, 145, 146] предложен метод, авторы которого моделировали одновременно поведение двух модельных электронов («электронных состояний», по терминологии авторов), рассматривая их совместную частоту рассеяния. При этом они предполагали, что состояние половины ансамбля электронов совпадает с состоянием первой модельной частицы, а состояние второй половины — с состоянием второй модельной частицы. Для получения достаточной статистики параллельно моделировалось несколько тысяч пар частиц. Число арифметических операций на шаге по времени в таком методе пропорционально числу пар. К сожалению, управляющие уравнения соответствующего марковского процесса никакого отношения не имеют к рассматриваемой физической ситуации.

Упоминание о другой возможности учета электрон-электронного рассеяния можно найти в [130]. Суть подхода заключается в введении индивидуального саморассеяния на механизме парных взаимодействий путем замены соответствующей функции сечения рассеяния ее максимальным значением. Механизм электрон-электронного рассеяния с введенной таким образом вероятностью рассеяния рассматривается в рамках стандартной процедуры метода Монте-Карло наряду с другими механизмами рассеяния. Вычислительные затраты метода пропорциональны числу частиц. Однако достаточного обоснования метода даже с физической точки зрения предложено не было.

Таким образом, несмотря на то, что моделирование парных столкновений имеет уже немалую историю, признана необходимость создания эффективных, строго обоснованных алгоритмов учета парных столкновений носителей заряда в численных моделях.

Предметом исследования пятой главы настоящей диссертации являются эффективные численные алгоритмы решения уравнения Больцмана, описывающего на кинетическом уровне эволюцию многодолинной полупроводниковой плазмы. Особое внимание уделяется учету нелинейных процессов (электронэлектронных взаимодействий). С математической точки зрения речь идет о решении системы нелинейных интегро-дифференциальных уравнений с больцмановскими интегралами столкновений, отвечающими различным механизмам линейного и нелинейного рассеяний.

При рассмотрении этого класса задач общий подход развивается на случай системы кинетических уравнений Больцмана для многодолинного полупроводника с функциями рассеяния достаточно общего вида. Обоснование алгоритмов метода стохастически взаимодействующих частиц для уравнения Больцмана теории полупроводников посвящены разделы 5.3, 5.4, 5.5.

Особенность моделирования методом частиц задач в много долинном приближении заключается в необходимости введения дополнительной характеристики сорта частиц (в смысле принадлежности к различным долинам) и учете их перехода из сорта в сорт. Аналогичная задача решалась в главе 2 для системы уравнений Больцмана, описывающей поведение многокомпонентного реагирующего газа. Наиболее естественный подход здесь заключается в простом изменении сорта частиц в соответствии с принципом пропорционального представительства. Так как заселенность центральных и боковых долин зоны проводимости может сильно различаться, такой подход требует большого объема статистической информации (числа частиц и числа испытаний), особенно при учете нелинейных взаимодействий. Другая возможность учета многодолинности заключается в построении весовых схем, в основе которых лежит введение для различных частиц зависящих от времени случайных весовых множителей [6, 85, 23]. Разделе 5.6 посвящен численным алгоритмам решения соответствующих рассматриваемому классу задач динамических схем на основе подходов, сформулированных в главе 1 с учетом особенностей моделирования полупроводниковой плазмы.

В разделе 5.7 в рамках многодолинной (TLX) модели проведено исследование изменения под действием электрон-электронных взаимодействий зависимостей стационарных дрейфовой скорости, средней энергии и заселенности долин от напряженности внешнего электрического поля в п-GaAs при различных концентрациях электронов. Проведенные численные расчеты позволили сделать вывод о существенности влияния электрон-электронных взаимодействий на перенос горячих электронов в n-GaAs в рассмотренном диапазоне электрических полей (3−15kBNocm). При этом, установлено немонотонное по концентрации электронов влияние е — е рассеяния на среднюю дрейфовую скорость ансамбля электронов (из линейных механизмов рассеяния учитывалось рассеяние на акустическом деформационном потенциале, полярных оптических и меж долинных фоно-нах).

Основные результаты диссертации.

1. Разработан новый подход к построению и анализу алгоритмов метода частиц для уравнения Больцмана, использующий специальное каноническое представление для слабой формы кинетического уравнения, основанное на введении модифицированных функций скачка, и заключающийся в выборе разностного аналога уравнения Больцмана в виде динамических схем стохастически взаимодействующих частиц, а также численного алгоритма решения соответствующей системы стохастических дифференциальных уравнений по мере Пуассона. На основе нового подхода впервые предложен математически обоснованный метод частиц решения уравнения Больцмана.

2. Предложен метод ветвей функции рассеяния построения динамических схем стохастически взаимодействующих частиц для уравнения Больцмана, в том числе схем с индивидуальными весовыми множителями. Метод ветвей обобщен на пространственно-неоднородный случай. Построены полевые и структурно-полевые динамические схемы для пространственно-неоднородного уравнения Больцмана, не требующие разбиения расчетной области на пространственные ячейки.

3. Разработаны эффективные алгоритмы численной реализации динамических схем, не требующие построения всей траектории системы стохастически взаимодействующих частиц и позволяющие определять приближенное (слабое) решение только в узлах некоторой неслучайной сетки по времени с вычислительными затратами, пропорциональными числу частиц. Доказана сходимость метода расщепления решения системы стохастических дифференциальных уравнений.

4. На основе предложенного подхода разработаны и обоснованы новые эффективные алгоритмы решения уравнения Больцмана для многокомпонентного разреженного газа с учетом внутренних степеней свободы молекул и химических реакций. С помощью построенных вычислительных схем решен ряд задач динамики разреженного газа, в том числе течения, моделирующего при определенных условиях обтекание межорбитального транспортного корабля при входе в атмосферу. Качественно описано влияние химической кинетики на структуру возникающего течения.

5. В рамках предложенного подхода разработаны новые стохастические алгоритмы решения нелинейного уравнения Больцмана для многодолинного полупроводника. Построены и обоснованы схемы с индивидуальными весовыми множителями. Предложены алгоритмы численной реализации динамических схем для уравнения Больцмана теории полупроводников с вычислительными затратами, пропорциональными числу моделирующих частиц. Проведено исследование влияния электрон-электронного рассеяния на транспорт носителей заряда в п — GaAs (ГЪХ-модель)с учетом рассеяния на акустическом деформационном потенциале, полярных оптических и междолинных фоно-нах. Для принятой модели установлено немонотонное по концентрации электронов влияние парных столкновений на среднюю дрейфовую скорость ансамбля электронов.

В заключение приношу огромную благодарность моему Учителю академику A.A.Самарскому за постоянную помощь и поддержку, академику О. М. Белоцерковскому, профессору А. А. Арсеньеву за интерес к работе и сделанные замечания, моим коллегам М. Н. Бурцевой, О. И. Зваровой, И. Л. Цветковой. И. Е. Юферову за большую помощь.

Я глубоко благодарен также моим немецким коллегам, особенно профессорам Х. Нойнцерту и Х. И. Кенигу за все, что они для меня сделали.