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

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

Механика

УДК 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.

Key words: electron irradiation, dielectriccharging, deterministic model, reaction-diffusion equation, parabolic partial derivative equation with hysteresis, computational scheme, finite-difference 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р).

88

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

(9)

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

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

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

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

(E, grad p) =

г

Dp =

Er dp + Ez Ъ-Р

r dr z dz

2

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

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

+Б

дг2

г

дг

д2 2

(11)

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

ееп

дг2 г дг дг2 о

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

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

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

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

др дг др э2

= 0, др дг

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

= 0,

(13)

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

при г > 0.

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

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

+

г=0

г=0

г=0

г=0

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

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

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

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

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

Заключение

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

1. Wu, J. Theory and applications of partial functional differential equations - N.Y.: Springer, 1996. - 368 p.

2. Kolmanovskii, V.B., Myshkis, A.D. Introduction to the theory and Applications of Functional Differential Equations - Kluver: Dordrecht, 1999. - 648 p.

3. Kyrychko, Y.N., Hogan, S.J. On the use of delay equations in engineering applications // Journal of vibration and control. - 2010. - V. 16 (7-8). - P. 943-960.

4. Wu J., Zou, X. Travelling wave fronts of reaction-diffusion systems with delay // Journal of dynamics and differential equations. - 2001. - V. 13. - P. 651-687.

5. Polyanin, A.D., Zhurov, A.I. Exact solutions of linear and nonlinear differential-difference heat and diffusion equations with finite relaxation time // International Journal of Non-Linear Mechanics. - 2013. - V. 54. - P. 115-126.

6. Jin, C.R. The numerical methods for solving some delay differential equations // Ph.D. Paper of Harbin institute of technology. - 2006. - P. 37-59.

7. Лекомцев, А.В. Численные алгоритмы решения некоторых классов эволюционных уравнений с запаздыванием: Автореф. дис. ...канд. физ.-мат. наук. - Екатеринбург, 2010. - 24 с.

8. Лекомцев, А.В., Пименов, В.Г. Сходимость метода переменных направлений численного решения уравнения теплопроводности с запаздыванием // Труды ИММ УрО РАН. - 2010. - Т. 16, № 1. - С. 102-118.

9. Maslovskaya, A.G., Sivunov, A.V. Simulation of electron injection and charging processes in ferroelectrics modified with the SEM-techniques // Solid State Phenomena. - 2014. - V. 213. - P. 119-124.

10. Сивунов, А.В, Масловская, А.Г. Численное моделирование процессов зарядки при диагностике сегнето-электриков методами растровой электронной микроскопии // Компьютерные исследования и моделирование. -2014. - Т. 6, № 1. - С. 107-119.

11. Pavelchuk, A., Barabash, T., Maslovskaya, A.G. Electron injection and polarization reversal processes in ferroelectrics analyzed with SEM: modelling of electron beam-stimulated effects // In: IOP Conf. Series: Materials Science and Engineering, 2016. - V. 110. - 012080 (6).

12. Maslovskaya, A.G., Barabash, T., Veselova, E. Polarization switching response and domain structure dynamics induced in ferroelectrics by incident electron beams // Solid State Phenomena. - 2016. - V. 247. - P. 131-137.

13. Павельчук, А. В., Масловская, А. Г. Конструирование разностной схемы для задачи расчета характеристик полевых эффектов в облученных электронными пучками диэлектриках // Материалы 7-й научно-практ. Internet-конференции «Междисциплинарные исследования в области математического моделирования и информатики», 2016 / отв. ред. Ю.С. Нагорнов. - Ульяновск: ЗЕБРА, 2016. - С. 147-150.

14. He, J., Tang, S.H., Qin, Y.Q., Dong, P., Zhang, H.Z., Kang, C.H., Sun, W.X., Shen, Z.X. Two-dimensional structures of ferroelectric domain inversion in LiNbO3 by direct electron beam lithography // J. Appl. Phys.. - 2003. -V. 93. - P. 9943-9946.

15. Яненко, Н.Н. Метод дробных шагов решения многомерных задач математической физики - Новосибирск: Наука, 1967. - 197 с.

16. Формалев, В.Ф., Ревизников, Д. Л. Численные методы - М.: ФИЗМАТЛИТ, 2004. - 400 с.

17. Reyes, E., Rodrigues, F., Martin, J.A. Analytic-numerical solutions of diffusion mathematical models with delays // Computers and Mathematics with Applications. - 2008. - V. 56. - P. 743-753.

18. Garcia, P., Castro, M.A., Martin, J.A., Sirvent, A. Numerical solutions of diffusion mathematical models with delay // Mathematical and Computer Modelling. - 2009. - V. 50. - P. 860-868.

19. Hong, B-S., Chou, C.-Y. Realization of thermal inertia in frequency domain // Entropy. 2014. - V. 16. -P. 1101-1121.

20. Физика сегнетоэлектриков: современный взгляд / под ред. К.М. Рабе, Ч.Г. Анна, Ж.-М. Трискона; пер. с англ. - М.: БИНОМ. Лаборатория знаний, 2011. - 440 с.