МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ПРОЦЕССА ЗАРЯДКИ ДИЭЛЕКТРИКОВ, ХАРАКТЕРИЗУЮЩЕГОСЯ ЭФФЕКТОМ ЗАПАЗДЫВАНИЯ

Математика . Прикладная математика.

Механика

УДК 519.6:51-73:004.942

А.В. Павельчук, Н.Л. Габрелян, А.Г. Масловская

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ПРОЦЕССА ЗАРЯДКИ ДИЭЛЕКТРИКОВ, ХАРАКТЕРИЗУЮЩЕГОСЯ ЭФФЕКТОМ ЗАПАЗДЫВАНИЯ

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

MATHEMATICAL MODEL OF DIELECTRIC CHARGING PROCESS CHARACTERIZED BY HYSTERESIS PHENOMENON

A mathematical statement of the simulation problem was proposed to describe evolutionary process of electron-induced charging in polar dielectrics with hysteresis phenomenon. The applied problem consists of modified representation of nonlinear "reaction-diffusion" type equation, model distribution of electron energy loses in irradiated target and set of the equations specifying coordinate dependences of an electric field-vector as well as polarization vector component stimulated by electron irradiation. There was formulated a hypothesis of computational scheme construction to realize the model using the grid splitting method.

Введение

В настоящее время реакционно-диффузионные системы имеют широкий спектр приложений для модельного представления процессов и явлений в гидродинамике, биологии, химии, физике, теории массо- и теплопереноса и других областях. Среди многочисленных моделей можно указать отдельный класс систем, формализуемых с помощью уравнений с частными производными в присутствии эффекта запаздывания или наследственности [1-3]. Уравнения такого типа также носят название функционально-дифференциальных уравнений [1]. При решении прикладных задач физический смысл запаздывания часто связывают с конечной скоростью распространения возмущений или инерционной природой самой системы, которая формирует отклик на внешнее воздействие не мгновенно, а с некоторым временным лагом. Исследованиям подобных систем посвящен значительный круг работ как зарубежных, так и отечественных авторов [1-8]. Следует заметить, что только для достаточно ограниченного класса задач в постановке уравнений с частными производными параболического типа с эффектом запаздывания можно построить аналитические решения [5]. Многочисленные исследования посвящены анализу свойств решений - устойчивости, асимптотичности, периодичности, осцилляции и пр. Вместе с тем для изучения объекта на основе постановки и проведения вычислительного эксперимента требуется построить гибкую и эффективную вычислительную схему. В связи с этим широкое распространение на практике получили численные методы, в том числе конечно-разностные [6-8].

Серия предшествующих работ [9-13] была посвящена развитию подходов к построению и реализации математических моделей стационарного и эволюционного процессов зарядки полярных диэлектрических материалов в неравновесных условиях электронного облучения. Специфика формулируемых задач соответствует теоретическому описанию процессов взаимодействия электронного зонда с сегнетоэлектрическими материалами, которые подлежат анализу и модификации с помощью методов растровой электронной микроскопии. В этом случае в качестве одного из аспектов полевого воздействия на полярный материал рассматривается эффект зарядки инжектированными в образец электронами. Подобное воздействие, в частности, используется для создания контролируемых доменных структур субмикронного размера [14]. Авторская модификация модели динамической зарядки была введена в рассмотрение с учетом собственной радиационно-стимулированной проводимости образца и на основе предварительной симуляции транспорта электронов методом Монте-Карло [9, 12]. В математической постановке модель описывается как начально-граничная задача для системы уравнений с частными производными, причем базовое соотношение представляет собой нелинейное реакционно-диффузионное уравнение. В указанных работах сконструированы различные варианты вычислительных схем, предназначенные для программной реализации модели и базирующиеся на использовании конечно-разностных схем расщепления (явной схемы расщепления, метода дробных шагов Яненко [15], метода переменных направлений [16]). Принимая во внимание физические законы и механизмы, определяющие диффузионный характер модели процесса зарядки, можно предположить, что эффект запаздывания реакции системы на внешнее воздействие может быть учтен в соответствующем соотношении математической постановки задачи. Поэтому данная работа направлена на разработку физико-математической модели эволюционных процессов зарядки полярных диэлектриков, описываемых уравнениями типа «реакция - диффузия» с учетом эффекта запаздывания по времени.

Концептуальная постановка задачи моделирования зарядки диэлектриков

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

Рассмотрим систему, включающую уравнение непрерывности, или уравнение сохранение заряда, и уравнение Пуассона. Дополним систему выражением, определяющим связь между потенциалом и напряженностью поля:

дР = О - Лу 1, дг

р=-р, (1) ££0

Е = -grad р,

где р — объемная плотность заряда, Кл/м3; \\ — плотность тока, А/м2; е — диэлектрическая проницаемость материала; е0 — электрическая постоянная, Ф/м; О - генерационное слагаемое, отвечающее за действие объемного источника в объекте, Кл/(м3с); р - потенциал, В; Е - напряженность поля, В/м.

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

= - В ■ grad р, (2)

где В = V ■ I - коэффициент диффузии электронов, м2/с; I - средняя длина свободного пробега, м; V - средняя тепловая скорость, м/с.

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

т — = -еЕ, (3)

е Ж У &

те - масса электрона, кг; е - заряд электрона, К.

Скорость движения электрона в конце свободного пробега (при начальной скорости -0 = 0),

до остановки при столкновении электрона с каким-либо атомом, равна: у1 =—г. Усредняя последнее

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

_ ёЕ- ёЕ .„.

V = — г =—г, (4)

где г - среднее время свободного пробега электрона между столкновениями с атомами.

Для г из последнего следует: г = . Средняя длина свободного пробега электрона может

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

У _ те-2 кТ

I = -г = —— = — . (5)

еЕ еЕ

В последнем выражении использована формула кинетической энергии теплового движения для одной степени свободы: т—2 = кТ 2 = 2 ,

где к - константа Больцмана, Дж/К; Т - температура, К. Таким образом, для коэффициента диффузии В с помощью (5) получим: г\\ — т —кТ кТ кТ

В=v ■1 = =тпЕ--=тп—. (6)

еЕ еЕ е

/ип = -— = —- дрейфовая подвижность электронов, м2/(Вс).

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

Г = аЕ, (7)

где а = тп Р - электрическая проводимость, См/м.

Таким образом, полный ток представляется суммой диффузионного и дрейфового токов:

j = jdf + Г . (8)

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

у jdr = а(\\у Е + (Е^а(а)

при а(ИуЕ = аР- = тР и (Е^га(а) = (Е^гаё(тпр)) = тп (E,gradр).

Диффузионная компонента для равновесного случая представлена выражением: Шу jdlf = - В Шу (grad р) = - О Ар .

Математическая постановка задачи моделирования электронно-индуцированной зарядки

диэлектриков с запаздыванием

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

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

(E, grad p) =

Er dp + Ez Ъ-Р

r dr z dz

d_2p+1 dp+d^p

dr 2 r dr dz2

Рис. 1. Схематическое представление геометрии образца и внутреннего источника.

Модифицируем представление (2) с учетом временного запаздывания:

f (r,z,t + tj ) = -D • gradp(r,z,t + tp), (10)

где tj и tp - лаги диффузионного тока и градиента плотности зарядов.

Учитывая, что t = tj -tp> 0, выражение (10) будет

означать, что система реагирует на изменение градиента плотности инжектированных зарядов изменением диффузионного тока не мгновенно, а с некоторым временным лагом t. Далее потребуем, чтобы время запаздывания t было учтено в (9) и при преобразовании уравнения непрерывности системы (1) с учетом суперпозиции вкладов диффузионного и дрейфового токов (8).

В данном случае можно провести некоторую аналогию с математическим моделированием процесса теплопроводности, характеризуемого запаздыванием по времени и представляющего собой диффузионный по природе процесс. В литературе описана так называемая модель двухфазного запаздывания (DPL - «dual-phase-lagging» model) для формализации закона Фурье и модификации

классического уравнения теплопроводности [17-20]. Это также приводит авторов к математической модели в виде начально-граничной задачи для функционально-дифференциального уравнения параболического типа.

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

др( г, г, г)

= О(г,7)-т р(г,7,()2 -т

( др( г, 2, г) др( г, 7, г)

Е & "

+ Е,

Г д2р( г, г, г-т) 1 др( г, г, г-т) д2р( г, г, г-т)

д 2р( г, г, г) 1 др( г, г, г) д 2р( г, г, г) р( г, г, г)

дг2 г дг дг2 о

Е (г, г, г ) = -§га<1 р( г, г, г),

где 0 < г < Я, 0 < г < Z - геометрические размеры объекта, т < г < Т - период действия источника, с начальным условием:

р( г, г, г )=р0 (г, г, г) при 0 < г < т, 0 < г < Я, 0 < г < 2 (12)

и соответствующими физическому смыслу задачи граничными условиями:

др дг др э2

= 0, др дг

= 0, р 2 = 0, р\\ Я = 0, (

= 0, Р, 2 = 0,р( Я = 0

при г > 0.

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

Спецификация внутреннего источника. Для определения геометрии внутреннего источника и введения аналитической аппроксимации функциональной зависимости распределения потерь энергии в облученном материале требуется предварительная симуляция транспорта электронов. Для моделирования распределения электронных траекторий может быть эффективно применен метод статистических испытаний или Монте-Карло [9, 11]. В вычислительной схеме положим, что электрон с энергией старта падает перпендикулярно плоскости поверхности образца в некоторую точку под определенным углом. Позиция электрона в точке задается с использованием значений углов рассеяния: азимутальным углом и углом отклонения. Значения углов и вид взаимодействия (упругое и неупругое) определяются с помощью генератора случайных чисел. Длина свободного пробега на каждом шаге рассчитывается с использованием модельного сечения рассеяния Мотта. Расчет изменения энергии при неупругом рассеянии электронов проводился на основе модифицированного закона Бете для многокомпонентных материалов. Изменение траектории и потерь энергии для каждого электрона рассчитывается до тех пор, пока величина энергии не уменьшится до некоторого порогового значения. Согласно концепции метода Монте-Карло для достижения статистической достоверности требуется моделирование для достаточно большого числа историй электронов (1000-10000). Моделирование транспорта электронов позволяет задать начальное распределение плотности зарядов в образце р0 (г, г, г) при решении задачи о моделировании релаксационных процессов или функцию, опреде+

ляющую генерационное слагаемое О (г, z), для задачи моделирования динамики зарядки диэлектрика.

Расчет координатных зависимостей вектора напряженности и индуцируемой электронным зондом компоненты вектора поляризации. В отличие от потенциала ф напряженность электрического поля Е является векторной функцией, которая в каждой точке пространства характеризуется величиной поля и направлением. Связь между напряженностью и потенциалом Е = — grad^ и применение

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

значение модуля вектора напряженности |Е| =^(Ег )2 + (Е2 )2 . Степень поляризации диэлектрика, индуцируемая инжекцией электронного пучка, характеризуется вектором поляризации Р, для вычисления которого можно воспользоваться связью между вектором напряженности и вектором поляризации: Р = (е — 1)е0 Е [20].

Комментарий к вопросу о конструировании вычислительной схемы реализации модели. Реализация детерминированной модели в постановке начально-граничной задачи для системы уравнений с частными производными, включая реакционно-диффузионное функционально-дифференциальное уравнение, потенциально может быть проведена на основе конечно-разностного метода переменных направлений [16]. Ранее алгоритм на основе указанного метода был построен для предшествующей локально-равновесной модели. Подобная вычислительная схема является абсолютно устойчивой для двумерной по пространственным координатам задачи и имеет второй порядок аппроксимации по времени и координатам г и z. Следует заметить, что в работах [6] и [7] авторами приводятся результаты конструирования вычислительных схем на основе конечно-разностного метода переменных направлений для отдельных классов функционально-дифференциальных уравнений параболического типа. Тем не менее вопрос построения, анализа и программной реализации вычислительной схемы для сформулированной постановки задачи (11)-(13) требует дальнейшего рассмотрения.

Заключение

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