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

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

Механика

УДК 51-7:519.6

Л.И. Мороз, А.Г. Масловская

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

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

Ключевые слова: модель переключения поляризации сегнетоэлектрика, дифференциальное уравнение дробного порядка, формула Грюнвальда — Летникова, предикт-корректорная схема, вычислительный эксперимент.

L.I. Moroz, A.G. Maslovskaya

PREDICTOR-CORRECTOR METHOD FOR THE NUMERICAL MODELING DOMAIN BOUNDARIES FRACTAL DYNAMICS IN FERROELECTRICS

The present study suggests computational scheme of predictor-corrector method for solving fractional-order differential equation applied to modelling polarization switching current in ferroelectrics. The mathematical model considers the fractal nature of domain structure dynamics and includes an initial problem for the fractional order differential equation. Predictor-corrector scheme was constructed with use of Gryunvald - Letnikov formula to solve the problem numerically. The method verification as well as program implementation of the model were performed. The computing experiments results were described to estimate electron-induced polarization switching current in ferroelectrics.

Key words: model of polarization switching in ferroelectrics, fractional order differential equation, Gryunvald - Letnikov formula, predict-corrector scheme, computing experiment.

Введение

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

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

Математическая постановка задачи и вычислительная схема

Пусть математическая постановка задачи моделирования фрактальной динамики некоторой физической системы в общем случае определяется задачей Коши для обыкновенного дифференциального уравнения дробного порядка:

(1)

с начальным условием: х(^0) = х0,

где х(0 - решение уравнения; / (х,/) - некоторая заданная функция правой части уравнения; а - порядок дробного дифференцирования: 0 < а < 1.

Вычислительную схему для реализации модели в постановке (1) будем строить с использованием формулы Грюнвальда - Летникова [6, 8, 9, 10] для численной аппроксимации производной дробного порядка на временной сетке соА = = /7г, /' = :

^=_!т1Ш,((4 (2)

где к - шаг по времени; Г (а) - гамма-функция Эйлера.

Подставляя аппроксимацию (2) в (1), получим явную конечно-разностную схему в виде:

( «-1 Л

(3)

=r(i)

V /=

Г(,7-g)

х

J I-_/+1 /=! )

где / = 0,Г -1; Л™ = ■

Г(-а)-Г(7 + 1)

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

метода одинакового или смежного порядка. Таким образом, на первом шаге рассчитаем «пробное» приближение Л"г+1, а при помощи неявной схемы уточним это значение, подставив его в правую часть:

(

\

j+1

i -0.N,

>-j+l

i - 0, N.

(4)

(5)

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

Данный метод легко распространяется и на системы дифференциальных уравнений дробного порядка [11].

Тест-пример

Для реализации математических моделей фрактальной динамики в постановке (1) с использованием вычислительной схемы (4)-(5) предложена прикладная программа в 111111 Matlab. С целью верификации работы программы построим численное решение следующего дифференциального уравнения половинного порядка с начальным условием х(0) = 0:

df

0<f<l.

Аналитическое решение задачи (6) имеет вид [10]: 28

315л/тг

Г

(6)

(7)

Результат сравнения точного и приближенного решений тест-задачи (6) приведен на рис. 1 а. На графике представлены: аналитическое решение, численное решение, полученное с помощью явной формулы (3), и численное решение, полученное на основе метода прогноза и коррекции (4)-(5).

Рис. 1. Сравнение численных решений, полученных с помощью явной схемы и схемы «предиктор-корректор», с аналитическим решением задачи (6) - я, абсолютная величина невязок численных решений - б.

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

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

Практическая оценка абсолютной погрешности численного решения задачи (6)

Метод Оценка погрешности \ численного решения

Норма Евклида 0.055 Равномерная норма £ = пих|Д|

Явная схема (3) 0.015

Схема «предиктор-корректор» (4)-(5) 0.019 0.005

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

Решение прикладной задачи

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

Математическая модель предложена для описания динамики 180° доменной структуры и формирования тока переключения поляризации сегнетоэлектрика. Ток переключения поляризации выражался с использованием соотношения для скорости движения доменной стенки. Принимая во внимание особенности переключения поляризации в режиме инжекции электронов, в модель вводили механизм клинообразного роста доменов: расширение области переключения путем роста клиньев с фиксированной шириной основания 2у и меняющимся углом наклона доменной границы ср. Процесс переключения поляризации заканчивается при переходе серединой клина z(t} = I + x(f)/2 границы кристалла. Выражение для вычисления поляризационного тока имеет вид:

I(w) =

4'Ps-d-y das

4-Ps-d-y

s2 -т.

dw°

L2

при 0 < s< 1 -

/

1 —

M

das dwa

(8)

где Ь - толщина кристалла, м; - площадь электрода, м2; .V = х/ /, - безразмерное расстояние;

м/—--безразмерное время; а - динамическая фрактальная размерность, варьируемая в диапазоне

т,

0.5 < а < 1; соэ° ср = у" / (х" + у ); тх=Ыул - характерное время пробега толщины кристалла доменной стенкой со скоростью V«,; / - глубина инжектирования пучка электронов в образец; М - число дискретных интервалов, на которые разбивается основание клина у при организации движения стен-

ки к - 0,М.

Динамика движения границы описывается уравнением дробного порядка:

(

das

: ехр

Л

x,w + —(l-1 2/V

2 \Ls -cos ф— ; 2

, *(wo) = 0, 0<s(w)<2(\-l/L).

(9)

где т = 21',. /_/ - характерный параметр - «время», в течение которого ток с плотностью у инициирует появление зарядов с поверхностной плотностью а = 21\..

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

Представляя уравнение (9) в виде

das idw~

приходим к общей постановке задачи (1), для решения которой воспользуемся сконструированной схемой (4)-(5).

Для проведения вычислительного эксперимента использовался следующий набор параметров, соответствующих данным физического эксперимента [14]. Моделирование проводилось для сегнето-электрического кристалла триглицинсульфата (ТГС). Геометрические параметры образца: толщина кристалла ¿=10° м; облучаемая площадь пятна на верхней грани кристалла равна .V 0.25-КГ' м2; линейные размеры облучаемой области кристалла d=2y=0.5-10° м (число областей, на которые разделено основание области переключения и=20). Физические параметры: диэлектрическая проницаемость образца 8=50, спонтанная поляризация Р5=2-10~~ Кл/м~. Электрическая постоянная в0 =8.85-101" Ф/м.

Характерные параметры времени: Т1=13.5 с, т2=37 с, Тз=16.5 с, которые соответствуют значению поля активации 5=6-105 В/м и плотности тока инжекции у—2.4-10° А/м2 (при токе зонда 1= 6-Ю10 А).

На рис. 2 представлен результат моделирования импульса тока переключения поляризации для кристалла ТГС в сравнении с экспериментальной кривой [14] (расчет проведен при а = 0.73). Зависимость поляризационного тока от времени характеризуется наличием максимума и повторяет особенности экспериментальных данных. До момента, соответствующего началу переключения поляризации, проходит время, необходимое для накопления стартового заряда, создающего в необлученной части кристалла поле, равное коэрцитивному.

.........I -2 -

- [ ь i / ' J* А ¡г VA ; Т

у

\\

\ « ч

J vx ..........i..........1.........

20 25 /, С

Рис. 2. Модельное представление тока переключения поляризации кристалла ТГС (2) и данные эксперимента (1) [14].

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

Заключение

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

1. Scott, J.F. Fractal dimensions in switching kinetics of ferroelectrics - Cambridge: University of Cambridge Press, 1998. - 9 p.

2. Galiyarova, N.M., Bey, A.B., Kuznetzov, E.A., Korchmariyuk, Y.I. Fractal dimensionalities and microstructural parameters of piezoceramics PZTNB-1 // Ferroelectrics. - 2004. - V. 307. - P. 205-211.

3. Мейланов, P.П., Садыков, С.А. Фрактальная модель кинетики переключения поляризации в сегнетоэлек-триках // Журнал технической физики. - 1999. - Т. 69, № 5. - С. 128-129.

4. Tadic, В. Switching current noise and relaxation of ferroelectric domains // Eur. Phys. J. B. - 2002. - V. 28. -P. 81-89.

5. Maslovskaya, A.G., Barabash, Т.К. Multifractal analysis of electron beam stimulated process of polarization reversal in ferroelectrics // Physics Procedia. - 2012. - № 23. - P. 81-85.

6. Самко, С.Г., Килбас, А.А., Маричев, О.И. Интегралы и производные дробного порядка и некоторые их приложения - Минск: Наука и техника, 1987. - 688 с.

7. Maslovskaya, A.G., Barabash, Т.К. Dynamic simulation of polarization reversal processes in ferroelectric crystals under electron beam irradiation // Ferroelectrics. - 2013. - V. 442. - P. 18-26.

8. Scherera, R., Kallab, S.L., Tangc, Y., Huang, J. The Griinwald - Letnikov method for fractional differential equations // Computers & Mathematics with Applications. - 2011. - V. 62. - P. 902-917.

9. Петухов, A.A., Ревизников, Д.Л. Алгоритмы численного решения дробно-дифференциальных уравнений // Вестник МАИ. - 2009. - Т. 16. - С. 228-234.

10. Xiao, К., Zhou, S.-B., Zhang, W.-W. Numerical Solution for Fractional-Order Differential Systems with Time-Domain and Frequency-Domain Methods // Journal of Electronic Science and Technology of China. - 2008. - V. 6 (3). - P. 294-298.

11. Биссенова, K.B., Масловская, А.Г. Предиктор-корректорная схема и численная реализация фрактальной модели динамики финансовых систем // Вестник АмГУ. Серия «Естественные и экономические науки». - 2018. -Вып. 81.-С. 11-18.

12. Bhalekar, S., Gejji, V.D. A predictor-corrector scheme for solving nonlinear delay differential equations of fractional order // Journal of fractional calculus and applications. - 2011. - V. 1, № 5. - P. 1-9.

13. Lekdee, N., Sirisubtawee, S., Koonprasert, S. Exact solutions and numerical comparison of methods for solving fractional-order differential systems // Proc. of the International MultiConference of Engineering and Computer Scientists (Hong Kong). - 2018. - Vol. 1. - P. 1-10.

14. Масловская, А.Г., Копылова, И.Б. Взаимодействие электронных пучков средних энергий с сегнетоэлек-трическими материалами - Владивосток: Дальнаука, 2010. - 204 с.