CC BY
28
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
Том 149, кп. 4
Физико-математические пауки
2007
УДК 519.63
ЭКСПОНЕНЦИАЛЬНЫЕ СХЕМЫ ДЛЯ РЕШЕНИЯ ЭВОЛЮЦИОННЫХ УРАВНЕНИЙ НА НЕРЕГУЛЯРНЫХ СЕТКАХ
C.B. Поляков
Аннотация
Рассматривается проблема численного решения эволюционных уравнений па нерегулярных сетках. Для уравнений типа конвекции-диффузии предложены консервативные слабомопотоппые копечпо-объёмпые схемы па треугольных и тетраэдральных сетках. Рассмотрены вопросы устойчивости и сходимости предложенных схем. Подход апробирован па задаче вакуумной микроэлектроники, связанной с расчетами полевой эмиссии электронов с поверхности острийпых кремниевых микрокатодов.
Введение
Традиционной проблемой математического моделирования процессов переноса заряженных частиц в полупроводниках и других материалах с нелинейно изменяющимися свойствами является наличие сильного перепада электрического поля в различных подобластях расчётной области. Оно приводит к тому, что концентрации частиц изменяются на порядки и даже на десятки порядков. В результате многие численные подходы сталкиваются с фундаментальным ограничением на величину приложенного к среде внешнего электрического или магнитного поля, избежать которого практически невозможно. Получение приемлемых по точности результатов достигается несколькими способами. Первый из них состоит в получении аналитического решения в областях простой формы в условиях слабой нелинейности и использовании набора таких решений в качестве базиса для представления искомой функции в более сложных ситуациях. Другой способ связан с построением адаптивной локально сгущающейся пространственной сетки, позволяющей выделить области сильного изменения электрического и магнитного полей. При этом может применяться перестроение сетки при переходе па очередной слой по времени. В дополнение к этим двум способам может применяться увеличение числа узлов сетки в расчёте на возможное распараллеливание вычислений.
Последние два способа (использование подробной адаптивной сетки и параллельных вычислений) представляются наиболее универсальными, поскольку имеют ограничения лишь по величине доступных вычислительных ресурсов. Однако они предполагают, что вычисления проводятся по устойчивым сеточным алгоритмам (схемам). Наиболее распространенными способами построения таких схем являются метод конечных элементов (МКЭ) и метод конечных объёмов (МКО). В настоящее время все большее число эволюционных задач в областях со сложной геометрией решается методом конечных объёмов ввиду его простоты и экономичности [1. 2]. К тому же метод конечных элементов наиболее эффективен при решении стационарных задач. Эти обстоятельства позволяют ограничиться в данной работе рассмотрением схем. полученных в рамках МКО.
В случае расчётных областей прямоугольной формы известно множество разностных методов решения уравнений типа конвекции-диффузии, о которых идет
речь ниже. Основной задачей при построении таких методов было получение на сетке всего спектра свойств дифференциального решения, то есть разностная схема должна быть не только устойчивой и сходящейся к дифференциальному решению, но и проявлять такие фундаментальные свойства, как консервативность и монотонность (в сильном или слабом смысле в зависимости от вида дифференциального уравнения). Для уравнений рассматриваемого типа примерами таких схем могут служить однородные схемы Самарского [3]. слабо монотонные схемы Голанда [4]. схемы Кареткиной [5] и т. д.
Последний из указанных подходов фактически является развитием схем экспоненциальной подгонки [6], применяющихся при решении обыкновенных дифференциальных уравнений. На наш взгляд, он имеет большие перспективы при решении нелинейных нестационарных задач электродинамики. Основным его недостатком на сегодняшний день является ограничение сверху на величину 6 = тах (к |Е|/У0) (максимум произведения шага пространственной сетки на величину электрического поля, отнесенного к величине нормировочного потенциала). Однако это же
6
единицы и выше большинство схем начинает быстро ухудшаться по точности, а некоторые из них теряют устойчивость вычислений.
6
при вычислениях с двойной точностью) и хорошо передают качественные характеристики дифференциального решения, в том числе консервативность и слабую монотонность. Последнее было продемонстрировано в работе [7] при решении задачи о динамике фотоиндуцированных зарядов в полупроводниковой гетероструктуре.
6
вычисления с четвертной точностью, которая аппаратно реализована в рамках новой 64-битной платформы вычислительных систем.
Другим недостатком схем экспоненциальной подгонки является отсутствие обобщения на случай нерегулярной сетки. В частности, автору неизвестны их аналоги для треугольных и тетраэдральных сеток, которые все шире используются при решении практических задач индустриальной направленности. Восполнить этот пробел призвана настоящая работа. В ней построены такие схемы для случая произвольной нерегулярной треугольной сетки и проведено их обобщение на случай тетраэдральной сетки. Проведён теоретический анализ построенных схем. определены условия их устойчивости и сходимости. На примере решения задачи о моделировании электронной эмиссии из кремниевого катодного микроузла продемонстрирована эффективность применения экспоненциальных схем.
1. Постановка задачи в двумерном случае
Для иллюстрации предлагаемого подхода рассмотрим в произвольной двумерной ограниченной области П с кусочно-гладкой границей дП на конечном интервале времени [0, Т] квазилинейное параболическое уравнение вида
ди
— = Ьи + /, Ьи = , IV = К gl•ad и + Ей, (1)
записанное в безразмерных переменных (х, у, Ь) относительно неизвестной функции и(х, у, Ь). Здесь предполагается, что ёгу и grad - обычные операторы дивергенции и градиента, действующие в пространстве переменных (х, у), W - обобщённый вектор потока, Е = Е(х,у,€) = (ЕХ,ЕУ) - кусочно-гладкое заданное векторное поле, К - тензор второго ранга с кусочно-непрерывными компонентами кщ (х, у,Ь, и) {%,] = 1, 2), удовлетворяющими условию равномерной эллиптичности (кц > к0 > 0, 4к11к22 — (к12 + к21)2 > к1 > 0) оператора Ь в П щи Е = 0,
функция / = /(х, у, и) кусочно-непрерывна по совокупности своих аргументов в области Б х (0, Т) х (-те, те).
Уравнение (1) дополняется начальными
и(х,у, 0) = ио(х,у) (2)
и граничными условиями
и(х,у,^= д1(х,у,^, (х,у) € дБ1,
4 > 0, (3)
(^,п) = д2(х,у,г,и), (х,у) € ЗВ2,
Здесь дБ^ (1=1,2) - части границы дБ = дБ1 идБ2 , на которых заданы граничные условия различного типа (Дирихле или Неймана - Ньютона), (•, •) - скалярное произведение векторов в двумерном евклидовом пространстве, п - внешняя нормаль к границе области.
В дальнейшем предполагается существование классического решения начально-краевой задачи (1) (3). Ввиду этого будет выполняться и следующее интегральное тождество
— JJudxdy=j)(W,n)dl+ II $ АхАу, (4)
В дВ В
которое легко получить при интегрировании уравнения (1) по области Б, выноса производной по времени за знак интеграла и применении формулы Гаусса Остроградского в двумерном случае.
2. Построение конечно-объёмной схемы в двумерном случае
Следуя подходу, развитому в [8], построим для примера явную конечно-объёмную схему для задачи (1)-(3). Для этого введем в области D треугольную сетку Qxy, удовлетворяющую критерию Делоне и такую, что каждый её треугольник имеет центр масс, проектирующийся строго на внутренние части всех его сторон. Пусть для данной сетки число вершин треугольников равно N, а число треугольников - М. На отрезке [0, T] введем равномерную сетку Qt с шагом т и числом интервалов Nt. По аналогии со случаем регулярной прямоугольной сетки можем записать следующую конечно-объёмную схему
U j+1 _ Uj
-г = 1 ,...,N, j = 0,...,Nt-l, (5)
т
с начальным условием
U0 = u0(xi,yi), i = 1,..., N, (6)
н граничным условием Дирихле
Uj = gi(xi,yi, tj), (xi,yi) G dDi, j = 0 ,...,Nt. (7)
Здесь Uj - сеточный аналог функции u в узлах треугольной сетки Pi = (xi,yi), Л - сеточный аналог oneратора L, ^j - сеточный аналог правой части уравнения (1) с учетом граничного условия (3) па части границы 3Di. Ниже будут уточнены выражения для Л и у.
( xi , yi , tj )
уравнения (1) по времени на интервале [tj, tj+1] и по пространству по медианаль-ному контрольному объёму Vi с центром в точке Pi и внешней границей, образованной отрезками, соединяющими центры рёбер и центры масс примыкающих
к Рг треугольников. После деления полученного интегрального равенства на величину т|У| (где |У| — площадь контрольного объёма), применения теоремы о первообразной для производной и по времени и формулы Гаусса-Остроградского для слагаемого Ьи, а также формулы о среднем значении функции в окрестности точки (хг ) получаем следующее приближенное равенство (подробнее см. [8])
3+1 3
щ — щ
дVi 1
Здесь и3 = и(Хг,Уг,гз ), W3 = W (х,у,гз), /3 = / (Хг,Уг,гз ,и3), дУц = дУг П дБ<2 , дУг2 = дУ\дУц.
В приближенном равенстве (8) можем учесть, что внешняя граница контрольного объёма У г есть замкнутая ломаная, на каждом звене которой внешняя нормаль есть постоянный вектор. Предположим, что к точке Рг примыкает тг треугольников. Очевидно, если точка Рг лежит внутри области, то число звеньев ломаной равно 2тг, и все они также лежат внутри области. Если точка Рг - граничная, то к 2тг внутренним звеньям ломаной добавляется четное число пг граничных ребер, исходящих из Рг. Па каждом внутреннем звене ломаной зафиксируем некоторое среднее значение вектора потока ^ ^ его через W 3к, и вычислим внешнюю нормаль к звену п¿к и длину звена 1гк ■ Па каждом граничном звене ломаной можно определить среднее значение функции д2 , обозначив его (д2 )3р, а также длину звена 1гр. Если далее применить теорему о среднем значении функции в соответствующих интегралах, то вместо (8) получаем
3+1 з 1 2т? 1
щ ^т^ ^ + т кр + Я- (9) 1 г| к=1 1 г| р=1
Рассмотрим вопрос о вычислении средних значений потока и граничной функции. Приближенные равенства (8) и (9) получены со вторым порядком точности по пространству. Это следует из способа аппроксимации пространственных интегралов от слагаемых щ и / уравнения (1). Данный способ был выбран из соображений минимальности пространственного шаблона схемы. Поэтому будем стремиться сохранить второй порядок пространственной аппроксимации схемы и при вычислении средних значений потока. Для этого достаточно потребовать,
и
и
и
кусочно-постоянным во всей области). Последнее позволяет считать, что и средний вектор потока постоянен во всем треугольнике. Аналогично, на ребрах треуголь-и
и
средние значения граничной функции.
Прежде, чем перейти к определению средних, заметим, что равенство (9) можно упростить, если учесть, что в каждом треугольнике имеется единственный средний вектор потока W3к , а вектор нормали п¿к можно заменить выражением © I
в котором Iгк - ненормированный вектор конормали (то есть вектор, заданный
©
поворота на 90° против часовой стрелки. Можно также учесть, что два последовательных звена ломаной, входящих в одни треугольник, можно «спрямить», взяв вместо суммы двух конормалей результирующий вектор (он будет соединять
середины сторон треугольника, прилежащих к точке Рг, и иметь правую ориентацию), который также можно обозначить I . В результате таких тождественных преобразований (9) переходит в следующее равенство
3 +1 3
иг - иг
1 X © I:,) + ¿>)?р 1гр + Я. (10)
N к=Л|У<| 1
к= 1 р= 1
Найдем теперь выражения для средних значений потока W 3к . Для этого преобразуем сначала выражение для W :
V = grad и + Си (V = К-1 W, С = К-1 Е)
н запишем его покомпонентно:
ди ди
х = :хЩ у = Ъу+ уП"
Далее, выберем произвольную точку Р0 = (х0,у0) облаети Б и произведем в последних равенствах следующие экспоненциальные преобразования:
д
^х 5 У У
У
где
вх = ехр
С х ¿х
еу = ехр
Су ¿у'
1-У0
Затем проинтегрируем получившиеся равенства по произвольному треугольнику Тгк (к = 1,..., тг), примыкающему к Рг, с вершинами в точках Ргк1 (I = 1, 2, 3, нумерация которых осуществляется против часовой стрелки, начиная с точки Ргк1 = Рг) и применим формулы Грина:
j ехУх ¿х ¿у = + ехийу, j еуУу ¿х ¿у = — еуийх.
Тлк дТцс Тлк дТц.
Учтем далее, что Ух = к11Шх + fc12W^ Vy = к21Шх + к22Wy (где кз - элементы К-1
J ехк1^х ¿х ¿у + J ехк12 Wy ¿х ¿у = + ^ ехи ¿у, Тк Тлк ЭТлк
Jeyк2lWx ¿,хЛу + Jеук22 Wy ¿х ¿у = — ^ еуи ¿х.
Тгк Тгк дТлк
Вынесем теперь средние значения компонент потока за знаки интегралов и приблизим двойные интегралы значениями соответствующих функций в центре треугольника, умноженными на его площадь Бгк, а контурные интегралы в левых частях равенств с помощью формулы трапеций. В результате получим следующие приближенные равенства
^х){к + ^ (^у& « С «(21) ^х)1к + «(22) (Wy)3к « ЬЦ', (11)
г (2)
У
где
(1т) $гк
а-, = -
гк з
У^ (вхк
1т
.1=1
(2т) $гк
а-, = -
гк з
(еУ к2г'
.1=1
1к1
т =1, 2,
(1)
(2)
Чк
(угкз - Угк1) (ехи У гк2 - (ехи )У )гк1
(Хгк3 Хгк1) (еУи )У )гк2 - (еУи )У )гк1
гкЗ - (ехи) Ук1
)У >гк3 - (еУи )У >гк1
Система уравнений (11) совместна и имеет единственное решение. Это следует из свойств оператора К и невырожденности треугольников Тгк. Поэтому в конечном итоге можно записать следующие приближенные равенства:
(ж
хУк
(11) (1)
Ь(1) + к(12)Ь(2) (Ж у « к(21)Ь(1) + к(22)Ь(2)
°гк + кгк °гк , (ЖУ)гк ~ кгк °гк + кгк °гк ,
(12)
где к(тп) - элементы соответствующей обратной матрицы (т, п = 1, 2).
Нетрудно заметить, что при определении потоков из равенств (12) понадобится
ех еУ
на Рг. Это обстоятельство позволяет вычислять указанные величины по простым приближенным формулам:
У Чкга
(ех У ) гк 1
(еу> У Чкга
(еУ У )гк1
exp
exp
0.5 (Хгкт - хгк1) ( (^х)\кт + (^х)\к1
т = 2, 3.
°.5 (Угкт - Угк1 ) ( (Оу )Укт + (СУ )ук1
Средние значения граничной функции д2 определим, применив формулу трапеций для приближения граничных интегралов. В результате получим:
ы
3 д 1
гР | (92)^1 + | 1У2;.;,р2
(д2)У
(13)
где Ргр1 = Рг, Ргр2 - концы граничного ребра треугольника Тгр, часть которого входит в контрольный объем Уг, а само ребро принадлежит ЗВ2.
Если теперь в приближенных равенствах (10)—(13) заменить функцию и на ее разностный аналог и, то нетрудно видеть, что (10) переходит в уравнение конечно-объёмной схемы (5), в котором
1
т
Аи? = 4т\ £ ('^ 6 Ч' ^ = # + т крг
к=1
(14)
р=1
и соответствующие средние значения определяются по точным равенствам (11)
(13).
Для иллюстрации разработанного подхода рассмотрим частный случай краевой
К
кц = к22 = ко, к 12 = &21 = 0, вектор Ж = к0gradи + Еи, ана всей границе выполняются условия 3-го рода
(Ж,п)= д2, (х,у) е <9£.
У
Дополнительно предположим, что задача является линейной, то есть ко = = k0(x,y,t), g2 = — ц (и — и0), и пусть ц = const. В этих условиях можно записать построенную схему в более привычном виде, а именно выписать итоговый вид оператора Л и правой части у:
MJi = {Wx)3ik {1у)>к + (Uifc) -
ni
I I p=1 \ /
^ = fi + щ £ (i Ml +\ MQ lip, (is)
p=1
где
. 3 Aylk3 ((e^j — U?) — Aylk2 ((e^j — Uj)
(Wv)3, = H---^-:--:-^-:-L
2 Stk (l/ko)l + (?:c/ko)ik2 + (?x/ko)ik3
(Wy)
3 Axik3 ((ey){k2Ujk2 — j — Axik2 ((ey){k3Ujk3 — j
y>ik
2Sn
(l/ko)j + (ey/ко) jk2 + (ey/ko)jk3
0.5Axikm [(Ex/ko)ikm + (Ex/ko)j 0.5Ayikm {(Ey / ko ) jkm + (Ey/ko)j
(ex)3ikm = eXP
(ey )jkm = eXP
Axikm xikm xi? Ayikm yikm yi? m 2? 3.
Сделаем далее несколько замечаний. Во-первых, аналогичным способом можно построить целое семейство схем по времени
Ui+1 — Uj
= a (wj+1 + yf1) + (1 — а) ^Uj + j
(16)
i = 1,...,N, j = 0,...,Nt — 1,
j
где а - вес схемы, принимающий значения 0, 0.5 и 1.
Во-вторых, используя приведенную процедуру, нетрудно построить схемы для трехмерного уравнения типа (1) на тетраэдральной сетке с барицентрическими контрольными объёмами.
В-третьих, если реализация явных схем вполне очевидна, то в случае неявных нелинейных схем придется использовать итерации по нелинейности по пространственным переменным. В настоящей работе использовались как явные схемы, так и полу неявные (чтобы не использовать итерации по нелинейности). Для реализации последних использовался метод бисопряжеииых градиентов с диагональным предобусловливателем.
Сформулируем теперь следующую теорему:
Теорема 1. Конечно-объш,ные схемы (6), (7), (16) при а = 0 и 1 имеют, порядок аппроксимации Ф = О (Ь2 + т) «Ф = О (Ь2 + т2) пр и а = 0.5. Все три схемы
Рис. 1. Расчетная область
консервативны, и при выполнении условий
Л-шт <Сх, т < С2НЩапп1 а =2, 2,1 для а = 0,0.5,1, (17)
они обладают свойством слабой .монотонности, являются устойчивыми в обычном и асимптотическом смысле и сходятся к дифференциальному решению задачи (1)-(3) с точностью Ф в норме Ь2(Пху) х С (П^ ■ Здесь Ни кШп - средний и минимальный размеры рёбер треугольников пространственной сетки, С1, С2 - положительные константы, не зависящие от параметров пространственной сетки.
Доказательство Теоремы 1 состоит из нескольких этапов. Порядок аппроксимации и консервативность схем вытекают из способа их построения и могут быть проверены непосредственно. Доказательство устойчивости и сходимости вытекает из условий выполнения слабого принципа максимума. При этом в случае неявной и симметричной схем можно построить оценки численного решения на слое по времени в энергетической норме (см., например, [3, 9]).
3. Тестовая задача и результаты расчетов
Конечно-объёмная схема (16), (6), (7) была опробована на примере расчетов динамики электронов в кремниевом автокатоде. Постановка задачи включала следующую систему безразмерных уравнений (подробнее см. [10])
дп
— =¿41 (НУ]п, ]п = пцп Е + D1gY&,d{Dnn), (18)
д'1д
= ¿42 (Ну (¡п + ип, Е)-(гип- 1), С}п = и>п Е+ 02ё1-а.<\(0па1ип), (19) ё1у(еЕ) = -7 (1 - п), Е = - gradу, (20)
где п и тп = п Тп - концентрация и энергия электронов (Тп - электронная температура), ^ п и Ц п - ток и поток энергии эле ктронов, Е и у - вектор напряженности и потенциал электрического поля, а, Бп - нелинейные коэффициенты, зависящие от электронной температуры, ^ , Бк, е, ^ - безразмерные константы.
Уравнения (18)—(19) решаются в подобласти кремниевого катода имеющего форму прямоугольника с небольшим острием в левом верхнем углу (см. рис. 1). Уравнение (20) для электрического поля решается во всей области П = = [0, Ьх] х [0, Ьу], содержащей катод и вакуумный слой П2 , примыкающий сверху
к катоду. Приведенные; уравнения дополняются следующими начальными и граничными условиями
п(х,у, 0) = 1, 'Шп(х,у, 0) = 1, Е(х,у, 0) = 0, (х,у) £ (21)
(3п, п) = зI, (Яп, п) = ЯП (х, у) £ (22)
^(х, 0) = 0, ^(х,Ьу ) = Уа, (Е,п )=0, х = 0,ЬХ. (23)
п
раздела вакуум-кремиий ставятся известные условия сопряжения на нормальную
Е
е1ЕП1) = £2ЕП2), Е« = ЕТ2), (х,у) £ П1 П (24)
Задача (18) (24) возникла в связи с разработкой устройств микровакуумной электроники [10]. использующих в качестве силового элемента кремниевые автокатоды малой размерности (порядка нескольких микрон и меньше). Основной проблемой при их создании оказывается сильный разогрев катода вблизи эмит-тирующей поверхности (в нашем случае раздел областей в силу большой величины электрического поля, который приводит к быстрому выходу прибора из строя.
Для моделирования указанной задачи на основе разработанной выше методики был создан комплекс программ и проведено предварительное исследование электронных процессов в кремниевых катодах микронных размеров. В качестве иллюстрации приведем расчеты для области, изображенной на рис. 1. с размерами 1.5 х 3 мкм. Полуширина анода была равна 1.5 мкм, высота катода -1.5 мкм, высота острия 0.3 мкм. радиус скругления острия 20 им, полуширина основания острия 0.06 мкм. Потенциал на аноде был равен 1000 В. Расчет проводился на треугольной сетке с параметрами: число узлов сетки во всей области составляло 35480. в области катода 17660: число треугольников во всей области 70049. в области катода 34666. Целыо расчетов был поиск области катода, подвергающейся наибольшему воздействию электрического поля.
На рис. 2 и 3 показаны изолинии модуля электрического поля и энергии электронов в области острия во время переходного процесса. Видно, что электрическое поле глубоко проникает в объём острия, разогрев электронов достигает значительной величины в верхней его части. Именно это приводит к наблюдаемым на опыте весьма высоким плотностям тока автоэмиссии у вершины острия, а также, как можно полагать, служит причиной деградации вершины острия в процессе эмиссии [11]. Таким образом, расчеты позволяют найти участки катода, наиболее подверженные деградации. Дальнейшая цель моделирования состоит в оптимизации эмиссионных и прочностных характеристик катодов при их применении в таких устройствах, как автокатодные и нашлемные дисплеи, эмиссионные зонды и др.
Заключение
В работе представлен новый численный подход к решению эволюционных краевых задач со структурой пространственного оператора типа нелинейной конвекции-диффузии. Особенностью численного алгоритма является применение схем экспоненциальной подгонки в случае решения задачи на нерегулярных треугольных и тетраэдральных сетках. В работе проведено подробное построение экспоненциальной конечно-объёмной схемы для модельного уравнения параболического типа общего вида, сформулирована теорема устойчивости и сходимости
Рис. 2. Изолинии модуля электрического поля в области острия катода
Рис. 3. Изолинии энергии электронов в области острия катода
схемы в энергетической норме. Эффективность предложенного подхода продемонстрирована на примере задачи вакуумной микроэлектроники по расчету эмиссионных характеристик кремниевого автокатодного микроузла в случае двумерной реальной геометрии.
Работа выполнена при финансовой поддержке РФФИ (проекты Л*1' 05-07-90054-в. 06-01-00097-а. 07-01-12079-офи).
Summary
S. V. Pulyakuv. Exponential finite volume schemes for solving of evolutional equations 011 irregular grids.
A numerical solution of evolutionary equations 011 irregular grids is considered. The conservative weakly monotone finite volume schemes 011 both triangular and t.et.raliedral grids for convection-diffusion type equations are proposed. The conditions of stability and convergence of these schemes are formulated. The proposed numerical approach are illustrated for one vacuum microelectronics problem connected with computations of electron emission from edge surface of silicon microcat.hode.
Литература
1. Chung T.J. Computational fluid dynamics. Cambridge: Cambridge Univ. Press, 2002. 1022 p.
2. Li R., Chen Zh... Wu W. Generalized difference methods for differential equations. Numerical analysis of finite volume methods. M. Dekker Inc.. 2000. 459 p.
3. Самарский А.А. Теория разностных схем. M.: Наука, 1983. 616 с.
4. Гола,ид Е.И. О сопряженных семействах разностных схем для уравнений параболического типа с младшими членами // Жури, вычисл. матем. и матем. физ. 1978. Т. 18, 5. С. 1162 1169.
5. Каре.ткииа П. В. Безусловно устойчивая разностная схема для параболических уравнений, содержащих первые производные // Жури, вычисл. матем. и матем. физ. 1980. Т. 20, 7. С. 236 240.
6. Дула,и Е., Миллер Док:., Шилдерс. У. Равномерные численные методы решения задач с пограничным слоем. М.: Мир, 1983. 198 с.
7. Поляков С.В., Сабликов В.А. Латеральный перепое фотоипдуцироваппых носителей заряда в гетероструктурах с двумерным электронным газом // Матем. моделирование. 1997. Т. 9, 12. С. 76 86.
8. Карамзин Ю.Н., Попов И.В., Поляков С.В. Разностные методы решения задач механики сплошной среды па неструктурированных треугольных и тетраэдральных сетках // Матем. моделирование. 2003. Т. 15, Л' 11. С. 3 12.
9. Калииичеико М.П., Поляков С.В. Численные методы для двумерной модели распространения лазерного излучения в химически активном газе в случае развитой термодиффузии // Жури, вычисл. матем. и матем. физ. 1997. Т. 37, Л' 3. С. 334 347.
10. Федирко В.А., Поляков С.В. Моделирование па МВС устройств вакуумной микро-электропики // Фундаментальные физико-математические проблемы и моделирование технико-технологических систем / Под ред. Л.А. Уваровой. М.: Япус-К, ИЦ МГТУ «Стапкип», 2004. Вып. 7. С. 138 147.
11. Дюжев П.А., Махиборода М.А., Федирко В.А. Исследование различных режимов автоэлектрошгой эмиссии кремниевого каптилевера // XIV пауч.-техп. копф. с участием зарубежных специалистов «Вакуумная паука и техника»: материалы копф. / Под ред. проф. Д.В. Быкова. М.: МИЭМ, 2007. С. 248 251.
Поступила в редакцию 02.11.07
Поляков Сергей Владимирович кандидат физико-математических паук, заместитель директора Института математического моделирования РАН.
E-mail: polyakoveimamod.ru
