«МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ С ПОМОЩЬЮ МНОГОПРОЦЕССОРНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМ ПРОЦЕССОВ ЭЛЕКТРОННОГО ТРАНСПОРТА В ВАКУУМНЫХ И ТВЕРДОТЕЛЬНЫХ МИКРО- И НАНОСТРУКТУРАХ ( ...»

ИНСТИТУТ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ

РОССИЙСКОЙ АКАДЕМИИ НАУК

На правах рукописи

ПОЛЯКОВ СЕРГЕЙ ВЛАДИМИРОВИЧ

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

С ПОМОЩЬЮ МНОГОПРОЦЕССОРНЫХ ВЫЧИСЛИТЕЛЬНЫХ

СИСТЕМ ПРОЦЕССОВ ЭЛЕКТРОННОГО ТРАНСПОРТА

В ВАКУУМНЫХ И ТВЕРДОТЕЛЬНЫХ МИКРО- И

НАНОСТРУКТУРАХ

(05.13.18 - Математическое моделирование, численные методы и комплексы программ)

ДИССЕРТАЦИЯ

на соискание ученой степени доктора физико-математических наук

Москва –

СОДЕРЖАНИЕ

Наименование Стр.

Введение Глава 1. Вычислительные основы решения задач электронного транспорта в микро- и наноструктурах 1.1. Математические модели электронного транспорта в твёрдотельных микро- и наноструктурах 1.1.1. Диффузионно-дрейфовая и квазигидродинамическая модели 1.1.2. Квантово-механические модели 1.2. Численные методы для анализа ДД и КГД моделей 1.2.1. Проблема численного решения краевых и начально-краевых задач для эволюционного уравнения общего вида 1.2.2. Схемы экспоненциальной подгонки на декартовых сетках 1.2.3. Схемы экспоненциальной подгонки на нерегулярных треугольных сетках 1.2.4. Схемы экспоненциальной подгонки на нерегулярных тетраэдральных сетках 1.3. Численные методы решения одномерных уравнений Фоккера-Планка и Шрёдингера 1.4 Разрешение проблемы некорректности одномерных краевых задач для уравнений Фоккера-Планка и Шрёдингера Глава 2. Параллельные алгоритмы и технологии 2.1. Параллельные алгоритмы на основе преобразования Фурье 2.2. Параллельные алгоритмы на основе метода прогонки 2.2.1. Базовый алгоритм распараллеливания 2.2.2. Обобщения базового параллельного алгоритма 2.3. Параллельные итерационные методы решения уравнения Пуассона и стационарных схем экспоненциальной подгонки 2.4. Параллельная реализация нестационарных схем экспоненциальной подгонки 2.4.1 Параллельные алгоритмы решения нестационарных задач на ортогональных сетках 2.4.2. Технология решения задач на нерегулярных сетках 2.5. Распараллеливание по группам и балансировка загрузки 2.6. Гибридная технология параллельного программирования Глава 3. Моделирование низкотемпературного примесного пробоя в полупроводниковых структурах Глава 4. Моделирование процессов латерального переноса фотоиндуцированных носителей заряда в гетероструктурах с двумерным электронным 4.4. Одномерная задача в условиях однородного освещения 4.5. Латеральный перенос в случае неоднородного освещения 4.5.2.1. Полуаналитический подход и его результаты Глава 5. Моделирование электронного транспорта в квантовых каналах Глава 6. Моделирование неравновесных процессов в ячейках автокатодных дисплеев и других автоэмиссионных микро- и наноструктур 6.5.1. Результаты моделирования в случае заданного распределения 6.5.2. Результаты моделирования реальной двумерной структуры 6.5.3. Результаты моделирования реальной трёхмерной структуры Глава 7. Моделирование процессов образования и миграции пор Введение Диссертационная работа посвящена созданию математических методов исследования, параллельных численных алгоритмов и комплексов программ вычислительных систем (МВС) неравновесных электронных процессов в эксперимента в относительно новой прикладной области, фундаментальный характер предполагаемых исследований, ориентация на использование вычислительной сложности анализируемых математических моделей.

Рассмотрим эти вопросы подробнее.

Исследование физических процессов в электронных приборах является достаточно развитой отраслью науки и опирается на успехи в различных разделах физики в конце XIX и начале XX веков. Расцвет этой отрасли начался в 1950-х годах с массового внедрения радио, телевидения и электронных вычислительных машин в повседневную жизнь. Развитие минимального энергопотребления. Номенклатура электронных устройств расширялась как по их назначению, так и по принципам действия. К концу моделирования как основной методологии научных исследований [1] появились первые технологии компьютерного моделирования, причём в самых высокотехнологичных областях науки, таких как атомная и ядерная физика, аэро- и гидродинамика и других. Поскольку эти технологии были изначально ориентированы на использование электронной вычислительной техники, они естественным образом проникли и в область исследований, связанную с развитием элементной базы компьютеров и систем связи.

К началу 1990-х годов (когда автор начал работу в данной области), в электронике сложилось уже большое самостоятельное направление – математическое моделирование электронных приборов [2-16], охватывающее множество фундаментальных и прикладных проблем из различных разделов физики и системотехники, использующее различные типы математических моделей, как на схемном, так и на физическом уровне. Поскольку основной технологией производства электронных схем и систем была твёрдотельная кремниевая технология (см., например, [16]), то большая часть исследований была связана с изучением различных физических процессов в кремниевых и других полупроводниковых структурах. При этом в фундаментальных работах по физике полупроводников появились разделы, связанные с приложениями в электронике, и наоборот, в приборных исследованиях – обширные разделы, связанные с физикой полупроводниковых систем. Это обстоятельство позволяет не противопоставлять одни работы в данной области другим.

Практически одновременно стали складываться несколько подходов к теоретическому и численному исследованию электронных процессов. Один из них базировался на достижениях классической механики, другой вобрал в себя статистические и квантово-механические методы. Иллюстрацией достижений в рамках первого подхода могут служить такие известные работы как [17-45], в рамках второго – [46-58]. До некоторого момента подобные исследования велись относительно независимо и объединялись лишь в рамках асимптотических методов многомасштабного моделирования и теории возмущений [59-65]. Однако, когда активные элементы электронных приборов преодолели микронный масштаб и стали перемещаться в область наноразмеров (середина 1990-х годов), возникла необходимость объединения всех подходов. Исследования, учитывающие одновременно классические и квантовые свойства субмикронной структуры появились как раз в это время и бурно развиваются до настоящего момента.

Со временем популярными стали работы по физике конденсированного состояния низкоразмерных (мезоскопических) и наноструктур. При этом спектр исследуемых материалов не ограничился полупроводниками.

Приведём список работ [66-113], который далеко не полон, однако показывает, как за последние 15 лет традиционные классические подходы сменялись неклассическими и смешанными, микроструктуры превратились в мезо- и наноструктуры, а теоретический анализ все более ориентировался на высокопроизводительные компьютерные вычисления.

Анализ эволюции методов математического моделирования за последние 15 лет показывает (и это подверждают, например, работы [109, 111-113]), что при моделировании реальных электронных систем низкой размерности необходим учёт различных геометрических и временных масштабов и сильно разнородных физических процессов в рамках одной задачи. Поэтому сегодня востребованы смешанные математические модели, включающие одновременно классические, полуклассические и квантовые описания электронных процессов в микро- и наноструктурах. При этом базовым подходом на макроуровне является чаще всего приближение механики сплошной среды (МСС), используемое как в классическом, так и в квантовом представлении. А на микро- и наноуровнях используются модели, учитывающие атомно-молекулярную структуру вещества.

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

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

статистического или имитационного моделирования, записанные как в классической, так и квантово-механической форме.

Учитывая особое внимание к микро- и наноструктурам с одномерным или двумерным электронным газом в связи с перспективами их применения в качестве базовых элементов новой электроники, выделим из общего ряда задачи моделирования одномерного электронного транспорта в квантовых структурах, которые чаще всего базируются на решении квазистационарной задачи туннелирования в приближении Хартри-Фока с учетом электронэлектронных взаимодействий и корреляционных эффектов. Именно этим моделям в диссертации было уделено особое внимание.

Сочетание в одной задаче классических и неклассических описаний приводит к существенному усложнению математической модели.

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

Поэтому в дисертации были развиты методы и подходы изначально ориентированные на параллельные вычисления.

Примеров успешного решения сложных физических проблем в выбранной прикладной области теперь уже достаточно [109, 112]. Однако в каждом конкретном случае могут понадобиться новые математические методы и их компьютерные и суперкомпьютерные реализации. Не является исключением и настоящая диссертационная работа, в которой большинство из использованных математических моделей были относительно новыми и не имели развитой вычислительной базы для их детального анализа и соответственно результатов моделирования. Экспериментальные данные в выбранных приложениях до сих пор остаются неполными и имеют фрагментарный характер, не позволяющий сделать однозначные выводы о природе исследуемых физических процессов.

Конкретные цели и задачи диссертации. Основной целью данной работы было создание математических основ для моделирования с помощью современных многопроцессорных вычислительных систем неравновесных электронных процессов в микро- и наноструктурах твердотельной и вакуумной электроники и проведение вычислительных экспериментов для ряда актуальных проблем. В качестве конкретных прикладных задач были рассмотрены проблемы моделирования процессов примесного пробоя и одномерного электронного транспорта в квантовых каналах наноструктур на основе AlGaAs, электронной эмиссии из кремниевых автоэмиттеров субмикронных размеров, образования и миграции пор в межсоединениях электронных схем в современных чипах в результате электрических и термомеханических воздействий. Все эти задачи объединяет сильная нелинейность и многомасштабность моделируемых процессов, а также повышенная сложность при численной, а затем и компьютерной реализации.

Следует подчеркнуть, что имеющиеся численные подходы и доступные на сегодняшний день комплексы программ до сих пор не позволяют решать рассмотренные в диссертации задачи в полном объёме и с необходимой точностью. Более того, многие из них появились под влиянием исследований аналогичных представляемому в данной работе и проведенных в последние 15 лет.

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

• на основе предварительного анализа разработать или выбрать наиболее адекватные математические модели изучаемых процессов;

• разработать или адаптировать к конкретным условиям эффективные • реализовать численные методики в виде комплексов последовательных и/или параллельных программ;

• провести детальное компьютерное моделирование исследуемых процессов и выработать рекомендации для специалистов из выбранных предметных областей;

• обобщить полученные математические результаты на случай более общих постановок задач из выбранных классов.

Методы исследования. В диссертационной работе применялся практически весь аппарат методов математического моделирования. На уровне моделей в основу работы легли модели механики сплошной среды и полуклассические квантово-механические модели (мотивация такого подхода подробно рассмотрена в п. 1.1 гл. 1). Основу численных алгоритмов составили сеточные методы решения обыкновенных дифференциальных уравнений и уравнений с частными производными. При построении пространственно-временных аппроксимаций дифференциальных уравнений использовались методы конечных разностей и конечных объёмов, применяемые как на ортогональных, так и на нерегулярных треугольных и тетраэдральных сетках (детали этого выбора обсуждаются в п. 1.2, 1.3 гл. 1).

При построении нерегулярных треугольных и тетраэдральных сеток использовались прямые методы типа Делоне. Для некоторых из построенных численных схем проводился анализ устойчивости и сходимости с помощью принципа максимума и/или метода энергетических неравенств.

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

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

Основные результаты работы можно сформулировать следующим образом:

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

2. Развиты математические модели электронной эмиссии с поверхности кремниевых микро- и наноструктур и электро-, термо- и стрессмиграции пор в медных межсоединениях электронных схем.

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

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

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

6. Созданы комплексы параллельных программ для моделирования процессов электронной эмиссии с поверхности кремниевых автокатодов и электро-, термо- и стресс- миграции пор в медных межсоединениях электронных схем.

7. С помощью разработанных вычислительных методик и комплексов программ получены новые результаты в исследовании процессов низкотемпературного примесного пробоя в полупроводниках типа GaAs, одномерного электронного транспорта в наноструктурах на основе AlGaAs, автоэлектронной эмиссии из кремниевого микрокатода с учетом реальной геометрии структуры.

Научная новизна полученных в диссертации результатов состоит в следующем.

1. В диссертации исследуются новые математические модели, описывающие электронные процессы в микро- и наноструктурах и разработанные автором совместно со специалистами и коллегами из Фрязинского отделения ИРЭ РАН, МГТУ «СТАНКИН», LSI Logic Incorporation. В качестве таковых можно указать модели примесного пробоя и одномерного электронного транспорта в наноструктурах на основе AlGaAs, модель неравновесного электронного транспорта в кремниевом автоэмиттере субмикронных размеров, модель процессов электро-, термо- и стрессмиграции в медных межсоединениях электронных схем.

2. В диссертации разработаны оригинальные численные методы анализа изучаемых процессов на базе как традиционных, так и новых математических моделей. Среди развитых численных подходов отметим конечно-объемные схемы экспоненциальной подгонки для решения начально-краевых задач для эволюционных уравнений на ортогональных и регуляризации и численные методы решения пространственно одномерных нелокально нелинейных квантовых задач в приближении Хартри-Фока.

3. В процессе работы над диссертацией автором созданы новые разработанные численные методы, а именно, параллельные реализации конечно-объемных схем экспоненциальной подгонки на ортогональных и нерегулярных сетках в одномерном, двумерном и трехмерном случаях, а также параллельный алгоритм решения нелокально нелинейной системы моделирования одномерного электронного транспорта в квантовом канале вычислителей.

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

• проведено численное моделирование задачи о примесном пробое в наноструктуре на основе AlGaAs, для которой в литературе имелись только теоретические оценки;

электронного газа наноструктуры на основе AlGaAs с целью определения возможностей оптической диагностики таких структур на этапе роста (ранее для данной задачи известны были только результаты нескольких натурных экспериментов, проведенных во Фрязинском отделении ИРЭ РАН);

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

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

Теоретическая и практическая ценность результатов диссертации заключается в ниже следующем.

• Совместно со специалистами в области твердотельной и вакуумной микро- и наноэлектроники разработаны новые математические модели, более адекватно воспроизводящие моделируемые процессы и явления.

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

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

• На основе предложенных математических моделей и разработанных численных методов созданы комплексы программ для персональных и суперкомпьютеров, позволяющие проводить детальные исследования задач из выбранных прикладных областей. Один из комплексов внедрен в промышленную систему моделирования.

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

Личный вклад соискателя. В диссертационной работе представлены результаты, полученные при решающем вкладе соискателя. Основные результаты диссертации получены лично соискателем. Исключение составляют формулировки математических моделей для рассмотренных в работе приложений и физические интерпретации полученных результатов моделирования. В частности, физико-математические модели низкотемпературного примесного пробоя, фотоиндуцированного классического и неравновесного квантового электронного транспорта в селективнолегированных гетероструктурах с двумерным электронным газом (см. гл. 1, 3, 4) разработаны профессором В.А. Сабликовым из Фрязинского отделения кремниевых автоэмиттеров субмикронных размеров (см. гл. 1, 6) разработана математическая модель образования и миграции пор под действием разработана совместно с профессором В.Я. Сухаревым и доктором Jun Ho Choy из LSI Logic Incorporation.

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

диссертационной работе, были представлены и обсуждались на всероссийских и международных конференциях:

• The Third International Congress on Industrial and Applied Mathematics (ICIAM'95), Hamburg (Germany), July, 3-7, 1995.

Зеленогорск, 26 февраля - 1 марта 1996 г.

• 23-rd International Simposium on Compound Semiconductors, 23- September, 1996, S.-Peterburg, Russia.

Технология - 96", С.-Петербург, 23-27 сентября 1996 г.

Технология - 97", С.-Петербург, 23-27 июня 1997 г.

• Третья международная научная конференция "Математические модели конденсированных системах и других средах", Тверь, 29 июня - 3 июля • Fourth Int. Conf. on Numerical Methods and Applications – NMA'98, Sofia (Bulgaria), 19-23 August 1998.

Технология", С.-Петербург, 22-26 июня 1998 г.

Технология", С.-Петербург, 14-18 июня 1999 г.

• Четвертый Всероссийский семинар "Проблемы теоретической и прикладной электронной оптики", Москва, 21-22 октября 1999 г.

"Полупроводники 99", Новосибирск, 25-29 октября 1999 г.

• 12-th Internatinal Vacuum Microelectronics Conference, IVMC'99, Darmstadt (Germany), July 6-9, 1999.

математическое моделирование и приближенные методы", посвященная памяти академика А.Н. Тихонова, Обнинск, 15-19 мая 2000 г.

• 4-ая Международная научная конференция "Математические модели конденсированных системах и других средах", г. Москва, 27 июня - июля 2000 г.

Технология", С.-Петербург, 19-23 июня 2000 г.

• 3-rd International Conference "Finite Difference Schemes: Theory and Applications (FDS-2000)", Palanga, Lithuania, Sept. 1-4, 2000.

• Всероссийская научная конференция "Высокопроизводительные вычисления и их приложения", Черноголовка, 30 октября - 2 ноября • 13-th Internatinal Vacuum Microelectronics Conference, IVMC'00, Darmstadt, Germany, July 6-9, 2000.

• 6th International Computational Accelerator Physics Conference (ICAP 2000), Darmstadt (Germany), Sept. 11-14, 2000.

COMPUTATIONAL PHYSICS – MTCP2000”, Dubna (Russia), July 24Int. Conf. "Displays and Vacuum Electronics (DVE 2001)", GarmishParteanKirche (Germany), May 2-3, 2001.

• Int. Conf. "Dynamical systems modelling and stability investigation", Kyiv (Ukraine), May 22-25, 2001.

• 4th International Vacuum Electron Sources Conference (IVESC'2002), Saratov (Russia), July 15-19, 2002.

• V International Conference on Numerical Methods and Applications – NM & A 02, Borovets (Bulgaria), August 20-24, 2002.

• Четвёртый Всероссийский семинар "Сеточные методы для краевых задач и приложения", Казань, 13-16 сентября 2002 г.

• 5th International Congress on Mathematical Modeling, Dubna (Russia), September 30 - October 6, 2002.

• Международная конференция "Математические идеи П.Л. Чебышова и их приложение к современным проблемам естествознания", Обнинск, 14-18 мая 2002 г.

• Пятый Всероссийский семинар "Сеточные методы для краевых задач и приложения", посвященый 200-летию Казанского государственного университета, Казань, 17-21 сентября 2004 г.

• IV International Congress on Mathematical Modeling, Nizhny Novgorod (Russia), September 20-26, 2004.

• Всероссийская научная конференция «Научный сервис в сети:

технологии параллельного программирования», Новороссийск, 18- сентября 2006 г.

• Международная научная конференция «Параллельные вычислительные технологии (ПаВТ'2007)», Челябинск, 29 января - 2 февраля 2007 г.

• Восьмой всероссийский семинар "Проблемы теоретической и прикладной электронной и ионной оптики", Москва, 29-31 мая 2007 г.

• Всероссийская научная конференция "Научный сервис в сети Интернет: многоядерный компьютерный мир. 15 лет РФФИ", Новороссийск, 24-29 сентября 2007 г.

• Седьмой Всероссийский семинар "Сеточные методы для краевых задач и приложения", Казань, 21-24 сентября 2007 г.

• XIV научно-техническая конференция с участием зарубежных специалистов «Вакуумная наука и техника», Сочи, 8-15 октября 2007 г.

• Всероссийская научная конференция "Научный сервис в сети Интернет: решение больших задач", Новороссийск, 22-27 сентября • Международная научная конференция "Моделирование нелинейных процессов и систем", Москва, 14-18 октября 2008 г.

• Девятый Всероссийский семинар "Проблемы теоретической и прикладной электронной и ионной оптики", Москва, 27-29 мая 2009 г.

• Internatilonal Conference «Mathematical Modeling and Computational Physics (MMCP'2009)», Dubna, July 7-11, 2009.

• Всероссийская суперкомпьютерная конференция "Научный сервис в сети ИНТЕРНЕТ. Масштабируемость, параллельность, эффективность", Абрау-Дюрсо, 21-26 сентября 2009 г.

Результаты работы обсуждались на рабочих семинарах ИММ РАН, НИВЦ МГУ, МСЦ РАН, РНЦ «Курчатовский институт».

Реализация и внедрение результатов работы.

Работа выполнена в Институте математического моделирования РАН и включает результаты пятнадцатилетних исследований соискателя по созданию математических моделей, численных методов, параллельных алгоритмов и комплексов программ для моделирования актуальных научно-технических задач в области твердотельной и вакуумной микро- и наноэлектроники.

математического моделирования РАН, проектов Программ фундаментальных исследований Президиума и Отделения математических наук РАН, проектов Российского фонда фундаментальных исследований, проектов ИНТАС, проекта Научно-технической программы Союзного государства «Разработка и использование программно-аппаратных средств Грид-технологий перспективных высокопроизводительных (суперкомпьютерных) вычислительных систем семейства «СКИФ», проектов Центра математического моделирования ИММ РАН – МГТУ «СТАНКИН», а также в рамках научного сотрудничества с компанией LSI Logic Incorporation (США) – производителем чипов для персональных и промышленных компьютеров.

Научные положения диссертации и разработанные на их основе методики, алгоритмы и программные комплексы использовались для совместных исследований в следующих организациях: Фрязинское отделение Института радиотехники и электроники им. В.А. Кательникова РАН, ФГУП «НИИФП им. Ф.В. Лукина», Центр математического моделирования ИММ РАН – МГТУ «СТАНКИН», LSI Logic Incorporation.

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

Основные публикации. По теме диссертации опубликовано 54 работы, из них 21 – статьи в ведущих отечественных и зарубежных рецензируемых журналах, в том числе 14 – статьи в российских рецензируемых журналах из списка ВАК. Основные публикации приведены в конце работы в виде отдельного списка [A1-A48].

Благодарности. Автор выражает глубокую благодарность руководителям Института математического моделирования РАН А.А.

Самарскому и Б.Н. Четверушкину, создавшим поистине творческую атмосферу в ИММ РАН и обративших внимание многих молодых исследователей на такое перспективное направление математического моделирования как высокопроизводительные вычисления в естественных и гуманитарных науках. Отдельно хочется поблагодарить Б.Н. Четверушкина за его постоянное внимание и поддержку в работе. Также автор выражает искреннюю благодарность своим учителям Ю.Н. Карамзину и В.А.

Трофимову за формирование научного мировоззрения и постоянное внимание к работе. Автор выражает глубокую благодарность своим старшим коллегам и соавторам В.А. Сабликову и В.А. Федирко за инициацию интереса к фундаментальным и прикладным проблемам микро- и наноэлектроники и неоценимую помощь в развитии физико-математических моделей, использовавшихся в диссертации. Также автор благодарен своим коллегам и соавторам И.В. Абалакину, Е.Н. Аристовой, В.Г. Бобкову, А.С.

Болдареву, П.Н. Вабищевичу, В.В. Вьюркову, В.А. Гасилову, Е.Н.

Головченко, И.А. Граур, Т.Г. Елизаровой, И.Г. Захаровой, Е.Л. Карташевой, Г.М. Кобелькову, С.Г. Кобелькову, Т.К. Козубской, Э.М. Кононову, О.А.

Косолапову, П.С. Кринову, Т.А. Кудряшовой, С.И. Мартыненко, О.Ю.

Милюковой, В.А. Николаевой, О.Г. Ольховской, И.В. Попову, А.А.

Свердлину, С.А. Сукову, Л.Ю. Тремсиной, И.В. Фрязинову, А.Г. Чурбанову, Е.В. Шильникову, М.В. Якобовскому за многочисленные обсуждения широкого спектра проблем вычислительной математики, математического моделирования, параллельных вычислений и программирования, затронутых в диссертации. Отдельную благодарность хочется выразить руководителям и сотрудникам Межведомственного суперкомпьютерного центра РАН Г.И.

Савину, Б.М. Шабанову, Ж.Е. Вершининой, О.С. Аладышеву, П.Н. Телегину и другим за многолетнее сотрудничество и поддержку при проведении вычислительных экспериментов на параллельных вычислительных системах МСЦ РАН. Хочется также поблагодарить свою семью за долготерпение и помощь при подготовке диссертации.

ГЛАВА Вычислительные основы решения задач электронного транспорта в микро- и наноструктурах моделирования электронных процессов в микро- и наноструктурах твердотельной и вакуумной электроники. В начале главы перечисляются используемые математические модели. Затем представляются разработанные математических моделей.

твёрдотельных микро- и наноструктурах Рассмотрим математические модели, которые будут использоваться в диссертации для описания процессов электронного транспорта в микро- и наноструктурах, применяющихся в твердотельной электронике. Как показывает анализ публикаций, посвященных теоретическим и численным исследованиям в данной предметной области, большинство используемых математических моделей основываются на детерминистическом подходе и лежат в рамках приближения механики сплошной среды (МСС) [114-117].

Объяснением этому факту может служить как широта возможностей для описания электронных процессов в рамках моделей МСС, так и большая работоспособность.

вариант метода Монте-Карло [118-122]), а также методы имитационного моделирования (например, методы классической молекулярной динамики [123-135]), имеют меньшее распространение вследствие их высокой вычислительной емкости и сложности анализа получаемых результатов. К тому же, очень часто статистические подходы объединяются с моделями сплошной среды. Например, широко используется метод частиц [136], квантовый метод Монте-Карло [137, 138] (объединяющий вариационный и диффузный варианты), квантовый метод молекулярной динамики [139, 140].

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

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

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

Современная электроника оперирует с приборами, в которых размеры активных элементов лежат в субмикронном и частично нанометровом диапазоне (условно от 500 до 10 нм). Например, силовые элементы современных чипов и питающие их линии по толщине имеют размеры от нескольких микрон до нескольких сотен нанометров. В то же время затворы процессорных транзисторов и переключателей имеют размеры 32 нм в производстве, 15-18 нм – в технологических экспериментах и дизайне, 5- нм – в научных разработках. Для такого рода объектов в литературе появился термин – мезоскопические структуры [141]. Под мезоскопическими структурами понимаются структуры, размеры которых на несколько порядков больше атомных, однако в них сохраняется квантовая когерентность частиц. В результате, в таких структурах реализуются квантовые эффекты, существенно влияющие на поведение структуры в целом. Для описания электронных процессов в мезоскопических структурах используются как классические, так и квантовые модели. Однако наиболее эффективными оказываются смешанные описания.

Учитывая сказанное, в работе используются все три подхода. При этом для описания классической динамики заряженных частиц в мезоскопических структурах выбрана гидродинамическая модель в диффузионно-дрейфововом [3, 26, 37] и квазигидродинамическом [6, 8, 12, 14, 15, 18, 142-155] приближениях. Для учета квантовых эффектов используются стационарные и квазистационарные квантово-механические модели, основанные на кинетическом уравнении Фоккера-Планка [33, 156-158] и уравнениях Шрёдингера в линейном приближении [159] и в приближении Хартри-Фока [160-165].

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

1.1.1 Диффузионно-дрейфовая и квазигидродинамическая модели Диффузионно-дрейфовая модель (ДДМ) выводится из кинетического уравнения Больцмана с помощью метода моментов в предположении о том, что при внешнем воздействии на изучаемую твердотельную структуру локально изменяются только концентрации заряженных и/или электрически нейтральных частиц, и при этом сохраняется полная энергия всей системы частиц. В такой ситуации ДДМ включает объемные уравнения неразрывности для концентраций свободных частиц и массовые балансные уравнения для концентраций связанных частиц. ДДМ дополняет уравнение Пуассона для потенциала самосогласованного электрического поля, если хотя бы одна группа частиц – заряженные. Если по какой-то причине изменяется энергия всей рассматриваемой системы (например, по структуре течет ток или она нагревается вследствие другого внешнего воздействия), то ДДМ дополняется уравнением энергии, которое в большинстве случаев можно записать в виде уравнения теплопроводности.

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

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

Здесь использованы следующие общепринятые обозначения: n и p – объемные концентрации свободных неравновесных электронов и дырок, nD и p A – объемные концентрации связанных электронов на донорах и дырок на акцепторах, jn и j p – плотности объемных токов свободных электронов и дырок, e – элементарный заряд, равный заряду электрона, n, p и Dn, D p – объемные коэффициенты подвижности и диффузии электронов и дырок, связанные соотношениями [31]:

( k = n, p, D, A ) – темпы генерации и рекомбинации соответствующих носителей заряда, удовлетворяющие принципу детального равновесия [31] и – напряженность и потенциал электрического поля, – теплоемкость при постоянном давлении, QT – плотность потока энергии, T – коэффициент теплопроводности, GT и RT – темпы генерации и релаксации энергии, div и – операторы дивергенции и градиента в декартовых координатах ( x, y, z ), t – время. Уравнения (1.1)-(1.4) записываются в ограниченной области с границей для t > 0.

Начальные условия для уравнений (1.1), (1.2) и (1.4) можно, например, принять в виде Условия на границе могут иметь весьма разнообразный вид. Для примера возьмем ситуацию, когда на границе заданы нормальные распределение потенциала 0 и температура окружающей среды T0 :

Здесь n – внешняя нормаль к поверхности области.

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

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

Здесь k – вектор-столбцы, составляющие тензор термоупругих напряжений, u – вектор смещений,, и K – безразмерные коэффициенты Ламэ и модуль всестороннего сжатия, и – функции нагрузки, возникающей при тепловом расширении структуры и изменении её массового состава, зависящие от температуры и массовых долей компонент среды. На границе структуры задаются либо компоненты вектора смещения, либо напряжения.

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

Для этого можно записать либо нестационарную систему уравнений Максвелла для основных характеристик поля, либо уравнений ЛифшицаЛандау для моментов, либо воспользоваться различными приближениями электромагнитостатики [167]. В диссертации этот вариант задач не рассматривается, однако разработанный численный подход применялся в [168] к решению системы одномерных нестационарных уравнений ЛифшицаЛандау при численном анализе процессов перемагничивания многослойной металлической пленки.

При рассмотрении процессов в многослойных структурах, содержащих металлы, полупроводники и диэлектрики, в каждом слое (подобласти) рассматривается конкретный вариант ДДМ и ставятся соответствующие условия сопряжения.

В диссертации ДДМ применяется для анализа процессов переноса фотовозбужденных носителей заряда в полупроводниковой гетероструктуре с двумерным электронным газом (гл. 4), а также для моделирования процессов образования и миграции пор в межсоединениях электрических схем (гл. 7). В последнем случае рассматривается многослойная структура металл-диэлектрик и учитываются электро- и термомеханические изменения ее геометрии.

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

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

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

В гл. 6 диссертации исследуется двух- и трехтемпературный варианты КГД модели для кремниевой структуры. Поэтому в дополнение к уравнениям (1.1), (1.3), (1.4) в этом случае добавляются одно или два уравнения для плотностей энергии электронов и дырок wn = nk BTn и wp = pk BTp ( Tn и Tp – температуры неравновесных электронов и дырок):

Здесь Q n и Q p – плотности потоков энергии электронов и дырок, n, p и Dn, D p – объемные коэффициенты дрейфа и диффузии плотности энергии электронов и дырок, Gn, G p и Rn, R p – темпы генерации и релаксации плотности энергии электронов и дырок. Начальные условия для уравнений (1.10) вытекают из (1.7). Граничные условия аналогичны (1.8):

В заключение пункта отметим, что ДД и КГД модели рассматриваются в пространственной геометрии, имеющей от одного до трех измерений, как в стационарном, так и в нестационарном случаях.

1.1.2 Квантово-механические модели Для описания электронного транспорта в мезоскопических структурах с учетом квантовых эффектов в диссертации используются три модели.

Первая модель базируется на кинетическом уравнении Фоккера-Планка [33, 156-158] для симметричной части функции распределения электронов по энергии f ( t, ) в нуль-мерном приближении:

Здесь A( ), B( ) – нелинейные интегральные коэффициенты, зависящие от энергии явным образом, G { f } и R { f } – генерационное и релаксационное столкновений, t > 0 – время, > 0 – энергия, отсчитываемая от дна зоны проводимости.

Стационарное уравнение (1.12) применяется в гл. 3 для моделирования явления низко-температурного примесного пробоя в полупроводниках типа GaAs. Оно учитывает диффузию и дрейф электронов в самосогласованном электрическом поле, процессы ударной ионизации, возбуждения и рассеяния взаимодействия. Для стационарного уравнения (1.12) формулируется нелокальная краевая задача в ограниченной области энергий p, M с условиями:

Здесь F { f } – нелинейный функционал от f.

Стационарную краевую задачу (1.12), (1.13) назовем моделью ФоккераПланка (МФП).

Вторая модель описывает одномерный стационарный квантовый перенос электронов через заданный потенциальный барьер. Она базируется на линейном стационарном уравнении Шрёдингера и применяется в гл. 6 для решения задачи туннелирования электронов через потенциальный барьер на границе твердое тело – вакуум. Условно назовем эту модель линейной моделью электронного туннельного транспорта (ЛМЭТТ). В размерном виде ее можно записать следующим образом:

Здесь = ( x) – волновая функция электронов с энергией, V ( x) – Здесь амплитуда падающей на барьер волны принята равной единице, что не ограничивает общности модели. В общей постановке задачи, когда барьер нельзя считать сосредоточенным на некотором конечном интервале, асимптотикой (1.15).

Смешанные краевые условия определяются непрерывностью волновой функции и ее первой производной на границах барьера:

отношению плотности потока в прошедшей волне к плотности потока в падающей волне:

Третья модель описывает квазистационарный электронный транспорт в квантовом канале полупроводниковой наноструктуры. Возникающая здесь задача туннелирования является нелинейной и многочастичной. Однако при использовании подхода Слэтера [163] и приближения Хартри-Фока [162] данная модель также формулируется в терминах одночастичных волновых функций и учитывает изменение эффективного потенциального барьера за взаимодействия. Также модель учитывает разделение электронных потоков по направлению и спину. В диссертации данная модель рассмотрена в одномерном случае для цилиндрически симметричной геометрии канала. Для данной геометрии потенциал самосогласованного электрического поля можно учесть в уравнениях Шрёдингера с помощью одномерной функции Грина.

Данную модель назовем условно нелинейной моделью электронного туннельного транспорта Шрёдингера в операторной форме без подробной расшифровки конкретных слагаемых:

Здесь ( ) ( x) – волновая функция электронов с энергией, спином s и направлением движения, U (s ) ( ) – нелинейный интегральный оператор, зависящий от всего набора волновых функций. Краевые задачи для уравнения (1.18) формулируются аналогично (1.15), (1.16).

пространственно нульмерном или одномерном квазистационарном случае и имели дополнительную координату в энергетическом пространстве. В этом электронного туннельного транспорта – квазидвумерные.

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

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

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

1.2 Численные методы для анализа ДД и КГД моделей В данном пункте рассмотрим основные численные алгоритмы, применявшиеся для решения поставленных в диссертации прикладных задач, обсуждающиеся ниже численные методы базируются в основном на методе сеток (по которому имеется обширнейшая библиография [170-236]) и рассматриваются в данной главе на примере решения отдельных уравнений, составляющих типовые компоненты используемых математических моделей.

При этом часть модельных уравнений анализируется в стационарном случае (рассматриваются краевые задачи), часть – в нестационарном (начальнокраевые задачи).

При решении дифференциальных уравнений с помощью метода сеток используются два его варианта: метод конечных разностей (МКР) (см., например, [170-212]) и метод конечных элементов (МКЭ) (см., например, [213-226]). Обобщением метода конечных разностей на нерегулярных сетках является метод конечных объёмов (МКО) (см. [188, 198, 200, 225, 227-236]).

Метод конечных элементов является вариантом вариационных методов Рэлея-Рица и Бубнова-Галёркина, в котором используются конечномерные пространства с базисом, состоящим из финитных функций. Методы конечных разностей и конечных объёмов совпадают с методом БубноваГалёркина в случае линейных базисных функций и медианных контрольных объёмов (см. [237-246]). По существу в МКЭ дискретизируется только расчетная область, в МКР и МКО дополнительно дискретизируются дифференциальные операторы.

В настоящее время метод конечных элементов наиболее широко используется при решении стационарных задач механики и теории упругости, поскольку позволяет использовать аппроксимации очень высокого порядка. В тоже время все большее число нестационарных эволюционных задач в областях со сложной геометрией решается методами конечных разностей и конечных объёмов ввиду устойчивости, простоты и экономичности этих методов. Эти обстоятельства позволяют ограничиться в данной работе рассмотрением схем, полученных в рамках МКР и МКО.

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

1.2.1 Проблема численного решения краевых и начально-краевых задач для эволюционного уравнения общего вида Традиционной проблемой математического моделирования неравновесных процессов переноса заряженных частиц в нелинейных средах является сильный перепад электрического и/или магнитного полей в отдельных подобластях расчётной области. В результате этого перепада концентрации частиц изменяются на порядки и даже на десятки порядков.

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

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

Последние два способа (использование подробной адаптивной сетки и параллельных вычислений) представляются наиболее универсальными, вычислительных ресурсов. Однако они предполагают, что вычисления проводятся по устойчивым сеточным схемам.

В случае расчётных областей прямоугольной формы известно множество конечно-разностных методов решения уравнений типа конвекциядиффузия. Основной задачей при построении таких методов является дифференциальное решение. То есть разностная схема должна быть не только устойчивой и сходящейся к дифференциальному решению, но и проявлять такие фундаментальные свойства как консервативность и монотонность (в сильном или слабом смысле в зависимости от вида параболического типа и их стационарных вариантов примером таких схем могут служить однородные консервативные монотонные схемы А.А.

Самарского [247], а также слабо монотонные консервативные схемы Е.И.

Голанта [248] и Н.В. Кареткиной [249]. В настоящей работе предложены их эволюционных уравнений, а также применения сеток различного типа.

евклидова пространства модельное эволюционное уравнение общего вида матричные, векторные и скалярные коэффициенты которого K = { K }, =1,2,3, координат r = ( x1, x2, x3 ), времени t и неизвестной функции u. В (1.19) div и – операторы дивергенции и градиента в декартовых координатах, (, ) – скалярное произведение в трехмерном пространстве. Будем предполагать, что уравнение (1.19) сохраняет параболический тип в области D (0, tmax ] ).

Для этого достаточно равномерной эллиптичности оператора L (см., например, [250, 251]). В частности, если выполнены условия (где k0 – некоторая константа, которая не зависит от координат, времени и значений искомой функции), то равномерная эллиптичность L имеет место, и уравнение (1.19) сохраняет параболичность во всей области.

Для уравнения (1.19) ставится следующая начально-краевая задача:

g12 + g 2 > 0, r D ; D – кусочно-гладкая граница области D, n – вектор внешней нормали к границе. В случае прямоугольной области D вместо условий (1.221) на границе области могут использоваться периодические условия (обозначим их условно (1.222)) или обобщенные дополнительные условия вида где ( r, t ; u ) – некоторый линейный или нелинейный функционал от u, а g – функция координат, времени и решения, определенная на границе области D. Примером условий (1.223) могут служить так называемые интегральные условия:

Основным требованием к задаче (1.19)-(1.22) является существование классического единственного решения на временном интервале (0, tmax ]. Для этого достаточно кусочной непрерывности по совокупности переменных всех исходных данных задачи [250, 251].

Перед обсуждением дискретизации уравнения (1.19) и для учёта экспоненциального характера его решения, определяющегося полями E и F, преобразуем пространственный оператор L. Для этого введем следующие векторные, матричные и скалярные переменные:

где ( x10, x2, x3 ) – произвольная конечная точка пространства, в том числе одна из точек области D, – -функция Кронекера. Тогда оператор (1.201) можно переписать в следующем виде:

Это представление будет использоваться ниже при дискретизации задачи (1.19)-(1.22). Нетрудно показать, что преобразования (1.23) не накладывают дополнительных ограничений на коэффициенты задачи.

При разработке сеточных численных методов решения задачи (1.19)независимо от пространственной размерности и наличия или отсутствия нелинейности используется интегро-интерполяционный метод [188, 200], который применяется как на ортогональных декартовых сетках (в том числе неравномерных), так и на нерегулярных треугольных и тетраэдральных сетках (как квазиравномерных, так и локальносгущающихся). Рассмотрим последовательно все варианты предлагаемого численного алгоритма.

1.2.2 Схемы экспоненциальной подгонки на декартовых сетках Как уже отмечалось, сеточные аппроксимации уравнений типа «конвекция-диффузия», к которым относится (1.19), известны давно. Так, в работе А.А. Самарского [247] рассмотрена методика построения консервативных монотонных разностных схем на ортогональной декартовой сетке для стационарного и нестационарного вариантов уравнения (1.19) для диагонального тензора K и E = 0. В этой работе предложена так называемая регуляризированная монотонная схема, имеющая второй порядок аппроксимации по пространственным переменным. С помощью принципа регуляризированной схемы к дифференциальному решению со вторым порядком точности по пространству. В [188, 200] регуляризированная схема обобщена на случай цилиндрически и сферически симметричной геометрии задачи.

регуляризированной схемы А.А. Самарского на ортогональной сетке для диагонального тензора K и двух вариантов одномерного уравнения (1.19):

когда E 0, F 0 и когда E 0, F 0. Для этих вариантов были предложены варианты схемы А.А. Самарского в исходном и сопряженном пространствах.

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

Отметим далее, что подход, разработанный А.А. Самарским и Е.И.

Голантом не использовал преобразования (1.23). В отличие от этих работ в статье Н.В. Кареткиной [249] такое преобразование проводилось. В ней был рассмотрен случай нестационарного уравнения (1.19) с диагональным тензором K и F 0. Предложенные в работе Н.В. Кареткиной разностные схемы фактически являются обобщением на нестационарный и многомерный случай схем экспоненциальной подгонки [252], широко применяющихся при решении краевых задач для обыкновенных дифференциальных уравнений.

Построенные в [249] аппроксимации отличаются от [247, 248] и достаточно, чтобы доказать существование и единственность разностного решения, а также устойчивость и сходимость схем к дифференциальному Дополнительно в [249] был рассмотрен вопрос реализации построенных схем с помощью метода немонотонной прогонки (он будет обсуждаться ниже в п.1.3 в связи с параллельной реализацией).

В настоящей диссертации было предложено естественное и очевидное объединение схем А.А. Самарского и Н.В. Кареткиной на случай E 0, F (см. [A27]). При этом аппроксимация экспоненциальных членов несколько отличалась от предложенной в [249]. Кроме того, было предложено в случае диагонального тензора K использовать локально-одномерные схемы обобщенной экспоненциальной подгонки. В работе [A42] было также показано, что схемы обобщенной экспоненциальной подгонки можно реализовать и на нерегулярных сетках, например, треугольных и тетраэдральных (см. п. 1.2.1.2), и этому не препятствует наличие смешанных производных.

Для иллюстрации разработанного численного подхода в случае ортогональных декартовых сеток рассмотрим стационарное уравнение (1.19) в области D = {( x1, x2, x3 ) : 0 < x < 1, = 1,2,3} :

В качестве граничных условий рассмотрим условия Дирихле:

представление оператора L в форме (1.202). Его можно записать подробно:

производных, методику [247, 248] для аппроксимации недивергентных представлений, можно получить сеточную задачу вида Здесь yh, f h и Lh – сеточные аналоги функций u, f и оператора L на, rh – узлы сетки, – граничные узлы сетки. Для Lh на полном 27-ми точечном нерегулярном шаблоне получаем следующее выражение, записанное в безиндексной форме:

где K h,, Vh, = Gh, yh и qh – разностные аналоги соответствующих непрерывных функций, K h, = K h,, = (1 + R ) к диагональным компонентам тензора, полученные по аналогии с [247, 248] в которых R = 0.5 Fh+ h+1) Fh h, Fh± = 0.5 Fh, ± Fh,. Функции Gh, в отличие от [249] аппроксимируются по формуле трапеций:

Заметим, что в качестве точки ( x10, x2, x3 ) выбрано начало координат.

Сделаем далее замечание относительно аппроксимации на границе в случае других типов граничных условий. Условия Неймана и Ньютона (входят в (1.221)) обычно аппроксимируются со вторым порядком по пространственным координатам. Для этого на границе записывается аналог (1.29), учитывающий граничные условия. Если задача периодическая, то никаких отличий от задачи Дирихле не возникает. Если используются обобщенные граничные условия (1.223), содержащие функционалы, то здесь возникает вопрос об аппроксимации функционалов со вторым порядком точности по пространственным переменным. Например, если эти условия можно записать в виде интегралов (1.224), то для их аппроксимации используются квадратурные формулы необходимого порядка точности.

Реализация схемы (1.27), (1.28) в линейном случае осуществляется методами решения линейных алгебраических систем, обсуждаемыми в гл. 2.

В нелинейном случае для реализации схемы (1.27), (1.28) предлагается использовать метод простой итерации по нелинейности:

Здесь s – параметры метода, Lh y h – оператор Lh на s-ой итерации. В качестве начального приближения можно выбрать любую ограниченную функцию, удовлетворяющую граничным условиям задачи. На каждой итерации (1.30) возникает линейная алгебраическая задача, решение которой также рассматривается в гл. 2.

Из общей теории решения систем нелинейных алгебраических уравнений следует (см., например, [253, 254]) утверждение:

Утверждение 1.1. При наличии у исходной задачи (1.27), (1.28) единственного ограниченного решения и выполнения условий итерационный процесс (1.30) сходится в норме C ( ) к решению задачи (1.27), (1.28) со скоростью геометрической прогрессии со знаменателем q.

В (1.31) y h – точное решение сеточной задачи, производных от Lh [ yh ], Ah – положительно определенный оператор вблизи корня y h, h – вектор, обращающий приближенное разложение в ряд Тейлора в равенство.

Доказательство утверждения 1.1 следует из принципа сжимающих отображений [253, 254].

Заметим, что использование метода Ньютона при решении задач рассматриваемого класса не всегда оправдано, поскольку он очень трудоемок в реализации и может иметь более жесткие по сравнению с (1.31) условия сходимости.

Рассмотрим вопрос о порядке аппроксимации схемы. Действуя аналогично [198] можно показать, что разностная схема (1.27), (1.28) с сеточным оператором (1.29) аппроксимирует задачу (1.24), (1.25) с сказанное про методику аппроксимации при других типах граничных условий, то легко понять, что и в этих случаях порядок аппроксимации не изменится.

Рассмотрим далее вопрос об устойчивости и сходимости построенной разностной схемы. Прежде всего, отметим, что устойчивость и сходимость схемы (1.27), (1.28) зависит от вида нелинейности. В частности, если исходная задача (1.24), (1.25) линейная, и все коэффициенты являются ограниченными кусочно-непрерывными функциями своих аргументов, то с помощью объединения методик [198] и [248] можно показать, что разностная схема (1.27), (1.28) устойчива и сходится к решению дифференциальной доказать такую же скорость сходимости в норме C ( ), если использовать теоремы вложения [198] или принцип максимума для сопряженного уравнения [248]. Результат можно распространить на случай граничных условий других типов (Неймана, Ньютона, смешанных, периодических, интегральных).

дифференциальной задачи (1.24), (1.25) результат предыдущих рассуждений не изменится, однако доказать его в каждом конкретном случае будет достаточно сложно. Здесь не будем рассматривать данный вопрос, поскольку ниже он исследуется для более общего случая нерегулярных сеток.

Суммируя сказанное можно сформулировать следующую теорему.

Теорема 1.1. Разностная схема (1.27) с граничными условиями общего вида и сеточным оператором (1.29) устойчива и сходится к точному решению дифференциальной задачи (1.24), (1.21) с пространственным выполнены следующие условия:

1) коэффициенты исходной дифференциальной задачи (1.24), (1.21) рассматриваемой области D, а их производные в области непрерывности ограничены;

2) решение дифференциальной задачи (1.24), (1.21) существует, единственно и является достаточно гладким в области D ;

3) шаги сетки удовлетворяют неравенствам:

где h – положительные константы, независящие от шагов сетки.

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

По указанным выше причинам доказательство теоремы проводить не будем. Отметим только, что оно проводится с помощью комбинации априорных оценок погрешности разностного решения и погрешности аппроксимации в сеточных нормах W21 ( ) и C ( ) [198]. Ограничения на K = { K, / G }, =1,2,3 должна быть равномерно положительно определена в узлах сетки.

Перейдем к конкретным применениям построенной схемы.

подгонки. Для этого возьмем в качестве D отрезок [0,1] и запишем краевую задачу для одномерного стационарного уравнения конвекции-диффузии:

с тремя типами граничных условий:

Здесь k, E, F, q, f – функции, которые могут зависеть от x и u, m, m, m – параметры, которые могут зависеть от решения, m (u ) – линейные или нелинейные функционалы по u.

подгонки во внутренних узлах неравномерной сетки x имеет следующий вид:

Здесь используются следующие стандартные обозначения: i – индекс узлов уравнений (1.35) имеют следующий вид:

где значения соответствующих функций на сетке с учетом нелинейности.

Отметим, что при реализации вычисление функций Gi не требуется. Вместо этого в расчетных формулах фигурируют отношения Gi +1 / Gi, которые легко вычисляются:

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

Далее рассмотрим нестационарный вариант задачи (1.19)-(1.22) с нестационарную схему экспоненциальной подгонки с весами yh = yh (t + ), yh yh (t ) (то же относится и к другим функциям), [0,1] – вес схемы, оператор Lh определен в (1.29).

Для реализации нелинейного варианта схемы с весами (также как и для стационарной задачи) предлагается использовать метод простой итерации:

где в качестве параметра метода выступает шаг по времени.

Исследование сходимости итераций (1.38) приводит к следующему утверждению.

Утверждение 1.2. При наличии у задачи (1.37) единственного ограниченного решения на каждом временном слое и выполнения условий итерационный процесс (1.38) сходится в норме C ( ) к решению задачи (1.37) на каждом временном слое со скоростью геометрической прогрессии со знаменателем q.

В (1.39) y h – точное решение задачи (1.37) на верхнем слое, h – вектор, обращающий приближенное разложение в ряд Тейлора в равенство. Доказательство сходимости, как и выше, опирается на принцип сжимающих отображений [253, 254].

Для схемы (1.37) можно доказать следующую теорему.

Теорема 1.2. Разностная схема (1.37) с граничными условиями общего вида и сеточным оператором (1.29) устойчива и сходится в норме L2 ( ) C (t ) к точному решению дифференциальной задачи (1.19)-(1.22) с (где = 1, если 0.5, и = 2, если = 0.5 ), если выполнены следующие условия:

1) коэффициенты исходной дифференциальной задачи (1.18), (1.20), (1.21) рассматриваемой области D (0, tmax ], а их производные в области непрерывности ограничены;

2) решение дифференциальной задачи (1.19)-(1.22) существует, единственно и является достаточно гладким в области D (0, tmax ] ;

3) шаги сетки удовлетворяют неравенствам:

где h, C – положительные константы, независящие от шагов сетки.

Доказательство теоремы в линейном случае можно провести методом энергетических неравенств [188, 198, 200]. Для этого необходимо записать задачу для погрешности разностного решения и провести необходимые оценки, используя ограниченность разностного и дифференциального решения, а также соответствующих производных последнего.

В нелинейном случае для явной схемы применяется методика, использующаяся в линейном случае. Для значений > 0 необходимо сначала доказать ограниченность разностного решения на новом слое, а затем рассмотреть задачу для погрешности разностного решения с учетом нелинейности. При этом также можно использовать энергетический метод.

Такая комбинированная методика применялась автором в работах [255-259].

В случае диагонального тензора K можно предложить использовать локально-одномерные схемы экспоненциальной подгонки.

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

где функции находятся с помощью решения серии одномерных задач, операторы Lh, [] имеют вид (1.34) по соответствующей координате и в сумме составляют оператор Lh [], вес = 0 или 1. Начальные условия для схемы (1.41) совпадают с условиями для схемы (1.37).

Схема с весом = 0 – линейная и реализуется с помощью комбинации прогонок. Схема с весом = 1 – нелинейная. Для каждого из этапов утверждением, аналогичным Утв. 1.2. Заметим также, что для обеспечения лучшей устойчивости и сходимости порядок применения одномерных операторов следует менять от шага к шагу (например, 1-2-3 на нечетном шаге и 3-2-1 – на четном) [210].

Теорема 1.3. Разностные схемы (1.41) с граничными условиями общего вида и сеточным оператором (1.29) устойчивы и сходятся в норме L2 ( ) C (t ) к точному решению дифференциальной задачи (1.18), (1.20), 1) коэффициенты исходной дифференциальной задачи (1.18), (1.20), (1.21) рассматриваемой области D (0, tmax ], а их производные в области непрерывности ограничены;

2) решение дифференциальной задачи (1.18), (1.20), (1.21) существует, единственно и является достаточно гладким в области D (0, tmax ] ;

3) шаги сетки удовлетворяют неравенствам:

где h, C – положительные константы, независящие от шагов сетки.

Если необходимо рассчитать существенно нестационарный процесс, вместо локально одномерной схемы (1.41), имеющей первый порядок двуциклического расщепления [210]:

yhj +5/8 yhj +3/ имеющую второй порядок аппроксимации по времени и безусловно устойчивую в норме L2 ( t ). Реализация и анализ сходимости схемы (1.43) проводятся также, как и для схемы (1.41).

Рассмотрим кратко случай криволинейной области D и применения аппроксимаций на нерегулярной сетке, не являющейся произведением сеток по отдельным пространственным направлениям. В этой ситуации нерегулярном шаблоне, предложенную в [188, 198, 200]. В итоге мы получим схемы экспоненциальной подгонки на нерегулярной сетке. При этом необязательно выбирать точку ( x1, x2, x3 ) внутри области D, поскольку соответствующая часть интеграла в реальных вычислениях не понадобится.

Альтернативным подходом является метод фиктивных областей (см., например, [209]). Для реализации схем экспоненциальной подгонки в рамках данного подхода достаточно доопределить в фиктивных областях уравнение диагональными компонентами тензора K.

В заключение пункта отметим, что предложенные обобщенные схемы экспоненциальной подгонки имеют большие перспективы при решении параболическом приближении. Преимуществом схем является то, что ограничение сверху на величину = max ( h E ) (максимум произведения электрического или магнитного поля), которое имеют любые разностные схемы, для схем экспоненциальной подгонки существенно ниже. Например, при увеличении до единицы и более большинство схем начинает быстро Экспоненциальные же схемы устойчиво работают и при больших (до 150единиц при вычислениях с двойной точностью) ввиду их продемонстрировано в работе [A7] при решении задачи о динамике фотоиндуцированных зарядов в полупроводниковой гетероструктуре (см. гл.

3). Указанный предел величины можно повысить до 30000 и более, если реализована в рамках новой 64-битной платформы вычислительных систем.

Недостатком представленных схем экспоненциальной подгонки было отсутствие обобщения на случай нерегулярных областей и нерегулярной сетки. В частности, автору не были известны их аналоги для треугольных и тетраэдральных сеток. Восполнить этот пробел удалось в работе [A42], в произвольной нерегулярной треугольной сетки и отмечена возможность распространения результата на случай тетраэдральной сетки, а также проведён теоретический анализ схем и определены условия их устойчивости и сходимости. Эти результаты рассматриваются в следующем пункте.

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

Принципиальным моментом при этом является свойство ортогональности используемых сеток.

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

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

Рассмотрим двумерный вариант уравнения (1.24) в произвольной замкнутой двумерной области D с кусочно-гладкой границей D. Для упрощения описания предлагаемого численного подхода удобно перейти от координат ( x1, x2 ) к координатам ( x, y ). Оператор L в форме (1.26) в этом случае будет иметь вид Компоненты вектора V в (1.44) имеют вид:

Граничное условие зададим с помощью вектора потока W с компонентами в следующем виде Здесь Wn – некоторая функция, зависящая от координат и решения u.

Как и выше предполагается, что двумерная краевая задача (1.24), (1.44), (1.47) имеет единственное решение, обладающее в D достаточной гладкостью. Искомая функция u ( x, y ) удовлетворяет принципу максимума и следующему интегральному тождеству которое легко получить, проинтегрировав уравнение (1.24) по области D и применив формулу Остроградского (см., например, [260]).

Перейдем к построению конечно-объемной схемы экспоненциальной подгонки на треугольной сетке для задачи (1.24), (1.44), (1.47), предполагая, что область D может иметь произвольно сложную форму, в том числе невыпуклую и/или многосвязную, с конечным числом точек разрыва гладкости границы.

содержащая как внутренние точки области D, так и точки ее границы D.

Множеством P будем обозначать все внутренние точки P. На P построена триангуляция про которую известно следующее:

1) T (P ) содержит все узлы P ;

2) образующие T ( P ) треугольники Tm имеют ненулевую площадь и пересекаются не более чем по образующим их вершинам или ребрам;

3) объединение треугольников Dh = Tm имеет ту же связность, что и 5) в случае > 0 величина стремится к 1 при бесконечном измельчении T (P ) с учетом криволинейности границы;

6) центр масс каждого треугольника проектируется на каждую его сторону (см. рис.1.1).

Рисунок 1.1. Пример тупоугольного треугольника, центр масс которого проектируется на все три его стороны.

Вопрос построения сетки, удовлетворяющей указанным требованиям, в данной работе подробно не рассматривается, хотя один из возможных подходов предложен автором в работе [261].

В общем случае триангуляция T (P ) может не удовлетворять критерию Делоне, и в ней могут быть треугольники с тупыми углами.

Однако такие треугольники являются граничными.

P = P \ P, а также множество граничных ребер треугольников E (P ).

Количество граничных ребер обозначим L.

Для построения конечно-объёмной схемы введем в области D дуальную к P сетку V = {Vi = V ( Pi ), i = 1,..., N }, состоящую из контрольных объемов. Для этого в каждой точке сетки Pi P определим множество всех треугольников, вершинами которых она является, H i = { Tm T (P ) : Pi Tm }, и назовем его шаблоном в точке Pi. Пусть количество точек в шаблоне равно N i + 1 (включая точку Pi ), а количество треугольников равно Mi.

Пронумеруем точки и треугольники в смежном порядке в направлении против часовой стрелки: Pi, Pi1,..., PiN, Tm1,..., TmM. При этом каждый треугольник Tm j образован точками Pi, Pi j, Pi j +1 (см. рис. 1.2), которые удобно также переименовать в Pi, Pm ), Pm+ ) (см. рис. 1.3). Заметим также, если точка Pi – внутренняя, то шаблон H i полный, и M i = N i, при этом N i 3. Если точка Pi – граничная, то шаблон H i неполный, и M i < N i, при этом Ni 2.

В каждом из треугольников шаблона определим точку пересечения медиан (центр масс треугольника) M m j ( j = 1,..., M i ). Соединим эти точки последовательно и получим замкнутую ломаную Li, окружающую точку Pi.

Фигура, ограниченная ломаной Li (зеленые линии на рис. 1.2), составляет контрольный объем Vi (заштрихованная область на рис. 1.2) в точке Pi с площадью Vi. Если точка Pi является граничной, то ломаную Li замыкается через проекции точек M j на соответствующие граничные ребра и саму точку Pi (см. рис. 1.2б). Последнее можно сделать ввиду выполнения условия 6) для триангуляции.

Введенный контрольный объём называется медианным контрольным объёмом. Его достоинство состоит в том, что он всегда существует, является выпуклым и пересекается с соседними объёмами только по границе. Такие контрольные объёмы используются при решении задач МСС методом конечных элементов. Недостаток медианных объёмов в рамках МКО – несимметричная аппроксимация самосопряженных дифференциальных операторов. Поэтому в рамках МКО используется барицентрический контрольный объём. Он отличается от медианного тем, что все звенья ломаных разбиваются на две части, каждая из которых лежит на медиане соответствующего треугольника. Новая ломаная Li состоит из 2 M i звеньев (отмечены красным на рис. 1.2), каждое из которых соединяет центр масс M m j (синие точки на рис. 1.2) одного из треугольников с серединой Qik (красные точки на рис. 1.2) одной из сторон этого треугольника.

Рисунок 1.2. Схема определения медианного и барицентрического контрольных объёмов для внутренних (а) и граничных точек (б). Синие точки – центры масс треугольников, красные – середины ребер треугольников, зеленые – пересечения ребер и отрезков, соединяющих смежные центры масс. Зелеными линиями и штриховкой выделен медианный объём, красными линиями – барицентрический объём.

Рисунок 1.3. Составные части барицентрического контрольного объема в треугольнике Tm j и направления нормалей и конормалей на их границах.

Заметим, что в качестве контрольных объёмов можно также рассмотреть ячейки Дирихле, в которых вместо центров масс берутся точки пересечения серединных перпендикуляров. Однако использование таких контрольных объёмов существенно ограничивает класс триангуляций, на которых проводится решение задачи. А именно, все треугольники триангуляции должны быть строго остроугольные. Иначе контрольные объёмы могут пересекаться или вовсе не существуют. Построенная на них схема в случае слабо тупоугольных треугольников может быть не консервативной.

В дальнейшем предполагается, что используются барицентрические контрольные объёмы Vi. Их преимущество состоит в том, что они могут быть построены на широком классе триангуляций, удовлетворяющих условиям 1)и попарно не пересекаются.

контрольного объема, представленные на рис. 1.2, 1.3. Как видно из рисунков, пересечение контрольного объема Vi и любого треугольника Tm j шаблона H i дает четырехугольник Tm j Vi Tm j = Tm ) Tm+ ), состоящий из двух смежных треугольных частей Tm ) и Tm+ ).

В дальнейшем используются обозначения: Pi = ( xi, yi ), Pm ) = xmj ), ymj ), Pm+ ) = xm+j ), ym+j ) – вершины треугольника Tm j, упорядоченные в направлении против часовой стрелки; Qmj ), Qm+j ) – середины сторон, исходящих из точки Pi ;

S m j, S m±j ) и S m j – площади Tm j, Tm± ) и Tm j (см. рис. 1.3). На границе Tm j, состоящей из четырех направленных отрезков определим векторы конормалей и нормалей к этим отрезкам по следующим формулам:

(m±j) = l (m±j) / | l (m±j) |, (m±j) = lm±j ) / | lm±j ) |, n (m±j) = (m±j), n (m±j) = (m±j), где = – оператор поворота на 90 градусов по часовой стрелке.

Также удобно использовать векторы и следующие соотношения:

независимую локальную нумерацию вершин, рёбер и их длин, а также обозначения сеточной функции в вершинах:

lmin = min lm +1), lmax = max lm +1) ;

Нумерация вершин производится в направлении против часовой стрелки.

Увеличение параметра на единицу производится циклически в рамках множества {1,2,3}, то есть 3 + 1 = 1.

Определения и соотношения (1.49)-(1.53) будут использованы для компактной записи разрабатываемой схемы.

Методика построения конечно-объёмной схемы состоит в том, что исходное уравнение разделяется на два с помощью введения вектора потока W, и каждое из них аппроксимируется отдельно путем интегрирования по например, многочленами Лагранжа. Такой подход применяется обычно для построения разностных схем на ортогональных четырехугольных сетках четырёхугольной ячейки, а значения вектора потока W определяются на сторонах ячейки. В данной работе предлагается использовать эту же тетраэдральных сетках. Причем, как и в [229], предлагается значения неизвестной функции искать в узлах сетки P, а потоки аппроксимировать специальным образом по этим значениям. Такой подход имеет определенное преимущество, в том числе позволяет решать задачу с граничными периодическими и функциональными).

Перейдем к построению схемы. Для этого проинтегрируем уравнение (1.24) по контрольному объему Vi и разделим на его площадь Si все части полученного равенства:

контрольном объёме Vi, умноженными на Vi. В результате получим следующее приближенное соотношение:

где использованы обозначения ui = u ( xi, yi ), f i = f ( xi, yi, ui ), qi = q ( xi, yi, ui ).

Порядок аппроксимации произведенной замены зависит от способа осреднения. Рассмотрим его подробнее на примере произвольной гладкой функции ( x, y ).

В простейшем случае можно положить i = ( Pi ). Тогда в случае произвольной треугольной сетки получаем ошибку порядка O ( l ), где l – величина наибольшего ребра треугольников ( l lmax ), входящих в шаблон H i. Если точка Pi совпадает с центром масс контрольного объёма Vi, а многоугольник Vi – правильный, то порядок аппроксимации повышается до O(l ). Также, порядок аппроксимации можно считать приближенно вторым, если положение точки Pi мало отличается от центра масс, а Vi – почти правильный многоугольник.

Последние два утверждения следуют из того факта, что следующие приближенные формулы

ABC ABC

являются обобщениями формул средних прямоугольников и трапеций для треугольника, S ABC – его площадь.

осреднение. В частности, в каждом треугольнике шаблона Tm j можно приближать подинтегральную функцию плоскостью, построенной по трем её значениям в вершинах треугольника:

Если использовать эти выражения при интегрировании, получим Тогда результирующая аппроксимация будет иметь вид:

Выражение для i в (1.551) можно также записать в виде В дальнейшем удобно ввести параметр, равный либо 0, либо 1, и умножить коэффициенты mj ) и m+j ) на него. Этот приём позволяет выбрать аппроксимацию того или иного типа:

Для аппроксимации правой части будем использовать формулы (1.55):

произведение различных средних:

Дальнейшая разработка схемы связана с аппроксимацией первого и второго слагаемых в (1.54). Рассмотрим сначала первое слагаемое (1.54).

Применим к нему по аналогии с (1.48) формулу Остроградского:

Здесь Vi – граница контрольного объема (то есть ломаная Li ). Контурный интеграл в правой части этого равенства разобьём на составляющие части вдоль звеньев ломаной и учтем, что векторы нормали на них не зависят от переменной интегрирования:

Здесь m±j ) – интегралы, возникающие только в граничных точках сетки. Во внутренних узлах сетки они отсутствуют, то есть m±j ) = 0.

На каждом отрезке ломаной введем средние значения вектора потока Wm±j ). Контурные интегралы в (1.56) будем аппроксимировать по формуле:

Отсюда с учётом граничных условий (1.47) получим:

Средние значения нормальной производной потока Wn(,± )j на границе области возьмем в виде аппроксимировать полусуммой значений Wn на концах ребра PQm±j ).

Средние значения потока Wm±j ) внутри области можно получить, если задать некоторую модель (аппроксимацию) вектора потока W в каждом треугольнике Tm j сетки. Если строить схему второго порядка точности, то достаточно предположить, что вектор потока W является постоянным в каждом треугольнике. Тогда Wmj ) = Wm+j ) Wm j. Конкретное выражение для Wm j в треугольнике Tm j и порядок такой аппроксимации уточним ниже.

Ввиду сказанного, аппроксимация первого слагаемого в (1.54) будет иметь вид построить по аналогии со случаем декартовой сетки. Поскольку W = KDV, то K 1W = DV. Это равенство можно проинтегрировать покомпонентно по треугольнику Tm j и воспользоваться формулой Грина [260]. Тогда получим:

В левой части полученных равенств применим формулы средних для компонент потока:

Контурные интегралы в правых частях (1.59) заменим по формуле трапеций с треугольника:

Значения V,i, V(,±m) j ( = x, y ) аппроксимируются следующим образом:

Аппроксимация значений G,i не требуется, поскольку в последующих расчетных формулах схемы они сокращаются.

Интегральные коэффициенты в левой части (1.59) обозначим как элементы матрицы размерности 2х2 и также заменим с помощью двумерной формулы трапеций в треугольнике:

В правой части (1.62) использована символическая запись.

Приближенные равенства (1.59) с учётом (1.60), (1.61), (1.62) приводят к следующей системе приближенных уравнений для аппроксимации компонент вектора потока W в треугольнике Tm j :

Система уравнений (1.63) совместна в силу выполнения условий равномерной эллиптичности оператора L во всех точках области D (что предполагается при выводе схемы) и при условии, что возмущенный тензор K мало отличается от K. Последнее действительно так, если треугольники сетки достаточно малы (что легко обеспечить при ее построении) и в каждом сеточном уравнении можно выбрать стартовую точку интегрирования ( x0, y0 ) при определении функций G независимо от других (это требование выполняется, поскольку, как и в случае ортогональной сетки, в реальных вычислениях будут участвовать лишь отношения значений функций G в соседних точках шаблона).

В итоге, значения компонент потока W приближаются формулами:

Рассмотрим кратко вопрос о порядке аппроксимации приближенных значений вектора потока. Для этого заметим, что если функцию u заменить в каждом треугольнике Tm j частью плоскости с использованием трех ее значений в вершинах треугольника (как это предложено выше), то погрешность такой аппроксимации будет иметь второй порядок. Градиент от Аппроксимация градиента при интегрировании по треугольнику или по любому контуру, лежащему внутри треугольника, также будет иметь второй порядок. Если теперь рассмотреть вектор-функцию V, то для ее компонент можно получить такой же результат, если интегралы в экспоненте аппроксимировать по формуле трапеций. Наконец, если рассмотреть аппроксимацию вектора W, то и в этом случае получим второй порядок ввиду того, что интегралы от компонент тензора K 1 аппроксимируются по обобщенной формуле трапеций в треугольнике. Таким образом, формулы Рассмотрим далее второе слагаемое (1.54). Заменим в нем интеграл по контрольному объёму на сумму интегралов по его составным частям Tm j, принадлежащим отдельным треугольникам шаблона (см. рис. 1.3):

Мы приняли, что в каждом треугольнике Tm j вектор потока постоянен, и его можно найти из (1.64). Поэтому можем записать цепочку приближенных равенств:

Из последнего приближенного равенства вытекает аппроксимация второго слагаемого в (1.54):

Объединим теперь (1.58) и (1.65) и подставим в (1.54):

Для монотонизации уравнений (1.66) в случае F 0 используем тот же приём, который применялся в случае ортогональных сеток [188, 200]. Для этого заметим, что вектор Fi S m j можно разложить по афинному базису между нормалями больше 0 и меньше ):

Поэтому Очевидно, что коэффициенты m±j ) стремятся к нулю при измельчении сетки со скоростью O ( l ). Схема (1.66) будет монотонной (устойчивой), если все величины 1 + m±j ) > 1, где < 1 – некоторое число, зависящее от параметров сетки. Если задать конкретное значение и подчинить размеры треугольников указанному выше условию, то схема (1.66) будет условно монотонной (условно устойчивой в норме C (P ) ).

Для построения монотонной (устойчивой в норме C (P ) ) схемы введем приближенное преобразование:

Если учесть, что m±j ) Cl ( C = const > 0 не зависит от параметров сетки), то величина в правой части последних равенств оценивается снизу величиной преобразование (1.67), (1.68) позволяет монотонизировать схему (1.66).

Завершая обсуждение аппроксимации уравнений (1.54), заметим, что и второе слагаемое в (1.54) было заменено на сеточный аналог с точностью O(l ). То есть уравнения (1.66) имеют второй порядок аппроксимации.

Далее, если вместо функции u в (1.66) подставить сеточную функцию U h = {U i U ( xi, yi ) U ( Pi ), Pi P } и перейти к равенствам, получим треугольной сетке для двумерной краевой задачи (1.24), (1.44), (1.47):

Сделаем замечание относительно аппроксимации различных типов граничных условий.

Во-первых, граничное условие (1.47) включает в себя условие Неймана Wn = 0. В этом случае величины m±j ) = 0 во всех точках сетки. Эта же ситуация возникает в случае условий Дирихле, поскольку уравнения (1.701) рассматриваются только во внутренних точках сетки.

Во-вторых, (1.47) может быть условием Ньютона, если Wn = u.

Тогда величины m±j ) = (U )i + (U )i + (U )m lm±j ). В случае условий смешанного типа, на каждой части границы получаем соответствующее типу выражение для m±j ).

В-третьих, при аппроксимации интегральных граничных условий величины m±j ) = Cm±j ),kU k, то есть они могут выражаться через значения искомой функции в любых точках сетки (нелокальная аппроксимация).

Похожая ситуация реализуется в случае периодических граничных условий.

Следует также отметить, что в нелинейном случае для реализации схемы (1.70) следует использовать итерационный процесс, аналогичный (1.30). Каждая его итерация реализуется с помощью решения системы линейных алгебраических уравнений с разреженной матрицей нерегулярной структуры. Методы решения таких систем рассматриваются в п. 1.3.

Рассмотрим далее некоторые варианты построенной схемы. Для простоты анализа остановимся на линейных граничных условиях Ньютона Аналог однородной схемы А.А. Самарского. Исследуем случай диагонального тензора K и E 0, F 0. Тогда схема (1.70) переходит в конечно-объёмный аналог на треугольной сетке двумерной однородной консервативной схемы А.А. Самарского:

Схема (1.71) имеет второй (или почти второй, в зависимости от поскольку является частным случаем (1.70).

Для сеточного оператора Lh [U h ], заданного во внутренних точках сетки, выполняется принцип максимума.

Утверждение 1.3. Если сеточная функция U h отлична от константы и удовлетворяет условиям то при l l0 она не может принимать положительного максимального значения во внутренних точках сетки.

Для доказательства этого факта рассмотрим тождества где Во внутренних точках сетки значения m±j ) = 0 и выполнены условия Следует также заметить, что во внутренних точках сетки шаблон схемы будет полным. Поэтому U m+j ) U mj +1, U m+ ) U m ) и суммы составляют внедиагональные коэффициенты схемы. Тогда В силу (1.74) имеем:

положительный максимум функции U h достигается в некоторой внутренней точке сетки Pi0 :

Поскольку функция U h не константа, а сетка P связная, то можно найти такой полный контрольный объём Vi1, в котором Последнее вытекает из того, что искомая функция u ( x, y ), по крайней мере, непрерывна в D. В случае сходимости схемы (которая доказывается ниже без использования оценок принципа максимума) её сеточный аналог U h наследует это свойство на сетке. Поэтому можно построить сетку такого малого размера ( l l0 ), что в окрестности положительного максимума будут находиться только неотрицательные значения U h.