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



УДК 004.942: 51-7

А.Г. Масловская Благовещенск, Россия

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

A.G. Maslovskaya Blagoveshensk, Russia

THE USE OF VARIATIONAL PRINCIPLES FOR COMPUTER SIMULATION OF FERROELECTRIC DOMAIN BOUNDARY DYNAMICS IN THE TEMPERATURE GRADIENT FIELD

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

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

The paper focuses on computer simulation of domain boundary pattern in nonuniform thermal field, which is based on variational approach. The numerical method of local variations was used to solve the model problem. Simulations allow supervising the ferroelectric domain boundaries dynamics as well as estimating the component of polarization current due to domain boundary movement in the temperature gradient field.

VARIATIONAL PROBLEM. FERROELECTRIC. DOMAIN BOUNDARY THERMAL FIELD. NUMERICAL ALGORITHM. SIMULATION.

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

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

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

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

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

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

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

Интерес к изучению доменных и межфазных границ и их динамики обусловлен определяющей ролью этих явлений в процессах поляризации и переполяризации сегнетоэлектрических материалов. Феноменологическое описание процессов движения доменных границ под действием термодинамических сил использовалось при описании и моделировании процессов переполяризации в ряде работ [3-7]. Так, в работе Алешина и Лучанинова [3] представлен расчет смещения стенок с помощью компьютерного моделирования для модельного полидоменного кристалла, обладающего слоистой структурой типа «паркет». В работе Нестерова [4] представлено детальное исследование механизмов движения доменных и межфазных границ в твердых растворах. Парлинский [5] построил и исследовал компьютерную модель поведения доменной стенки вблизи фазового перехода, в основе которой лежит модель одной ячейки, обладающей потенциальной энергией. Большинство исследований в этой области основывались на рассмотрении кинетики доменных и межфазных границ с макроскопической точки зрения. В макроскопических моделях кроме термодинамических сил вводят феноменологические силы вязкости различной природы (сухого, вязкого трения), иногда - инерции.

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

Следуя термодинамической теории сегнето-электриков Ландау-Гинзбурга-Девоншира допустим, что тепловое, упругое и диэлектрическое поведение однородного диэлектрика полностью описывается шестью величинами: полями температуры Т, энтропии S, напряжения X, деформации х, электрического поля E и электрической индукции Б. Термодинамической функцией, достаточно полно описывающей свойства сегнетоэлектри-ков, может служит свободная энергия - упругая энергия Гиббса Ох(Т, X ,, Б) [8]. В частном случае, когда приложенные механические напряжения отсутствуют, то есть X = 0, неполярная фаза цен-тросимметрична и Б. = Б - индукция направлена вдоль одной из кристаллографических осей, функция О1 имеет вид простого полинома:

О Л в2 +1Б4 + 8 Б6 +... . (1) 1 2 4 6

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

а = Р(Т - Тс), (2)

где в - положительная константа, связанная с константой Кюри-Вейсса С соотношением

CP =-

(е0 - электрическая постоянная).

•Р

В силу того, что О ~ Р при Е = 0, из условия

д°1 п

—1 = 0 нетрудно вычислить спонтанную поля-дО

ризацию Р3 ниже температуры Кюри для сегне-

тоэлектрических кристаллов с фазовыми переходами второго рода [8]:

PS =| (Tc - T).

Y

(3)

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

( Ф ^

ст =

•8-а • P

(4)

где

8 = .-

f

а -s„

- толщина доменной стенки;

/ ~ й2 - коэффициент в выражении вклада корреляционной энергии в свободную энергию кристалла; й - постоянная решетки.

Таким образом, для кристаллов с фазовыми переходами второго рода зависимость удельной поверхностной энергии ст от температуры Т с учетом (2) и (3) примет вид:

ст =

W2 d ß* - - т)3,

3 & Y '

(5)

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

Т (х) = Та + а • г, (6)

где Т0 - начальная температура для исходного

йТ

положения границы; а = — - градиент температуры.

Для нормировки параметров можно задать некоторый уровень отклонения границы от равновесного положения хс, соответствующий температуре Кюри ТС. Тогда

Т (х) = а( х - хс) + Тс. (7)

где ст0 = -

Подстановка последнего в (5) дает ст(z) = ст.(а(zc - z))3/2,

W2

3 л/ё. Y Математическая постановка задачи моделирования конфигурации доменной границы на основе вариационного подхода

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

Задачу о моделировании равновесной конфигурации в поле линейного градиента температуры можно сформулировать как вариационную задачу о минимизации функционала:

Ф( z (х, y)) = Цст( z (х, y ))dS =

'__(9)

= Цст(х, y, z) • у 1 + z'2 + z'^ dxdy ^ min

D

или в двумерном случае, с учетом (8):

X _

ф(z(x)) = |а3/2ст0(zH -z^^/T+zfdx^min, (10)

x0

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

z( x0) = z (X) = 0. (11)

На рис. 1 представлена схема изменения конфигурации доменной границы в тепловом поле.

Для общей постановки прикладной задачи (10), (11) можно ввести обозначения независимых

Рис. 1. Модельное представление конфигурации доменной границы

1

параметров в нормированных единицах: ст0 = 1, Т0 = -1, Тс = 0 , гс = 50. В этом случае, градиент температуры примет значение а = 0,02.

Для того чтобы убедиться, что экстремаль г (х) доставляет минимум функционалу

2+(Л" I

Ф(г(х)) = | а32СТ0 (гл - г)3/2^1 + г'2 ёх,

рассмо-

трим вторую производную подынтегральной функции Г(х, г, г') по г':

32 (гс - г)3/2

Г,, = ст„ - а

(12)

(1 + (гХ)2)32 •

Знаменатель выражения (12) положителен при любых значениях г', числитель также положителен при любых значениях г, поскольку формализация математической модели предполагает выполнение условия г < гс. Из этого следует, что выполняется усиленное условие Лежандра > 0 и функционал (10), согласно теореме о достаточном условии слабого минимума (максимума), на некоторой экстремали г (х) достигает своего минимума.

В общем случае вариационная задача в постановке (10), (11) не разрешается известными аналитическими методами и требует привлечения численных процедур.

Применение метода локальных вариаций для численного решения прикладной задачи

Построим численное решение вариационной задачи (10) с граничными условиями (11), используя метод локальных вариаций [9]. Разобьем отрезок [х0, X] на N равных частей с шагом

X - х _

Ах =-0: х . = х0 + ¡Ах, ; = 0, N . Согласно дан-

N ; 0

ному численному алгоритму введем обозначение для подынтегральной функции Г (х, г, г'): (13)

1,[ г,, г,+1] = Г (х + Ах /2,( г. + г,+1)/2, (г+1 - г)/ Ах) -Ах.

Функции г(х) будем аппроксимировать ломаными, соединяющими точки (х1, г1) в плоскости Охг, а функционал (10), с учетом (13), прибли-

N-1

женно заменим суммой J «I = ^ I [г.., г..+1].

¡=0

Решение задачи (10), (11) будем строить последовательными приближениями г;(п), п = 0, 1, 2, ... , выбирая в качестве начального приближения значения г 0 , удовлетворяющие ограничениям (11). В качестве возможных значений для г.+1

к к 7 к 7

рассмотрим три величины: г ; , г ; + п, г ; + п при

/-К

'■■О----О---0'

Рис. 2. Схема метода локальных вариаций [9]

априорно заданном шаге варьирования п. Далее проводим подсчет сумм:

Ф,= I.-,[ г к-11, г. ] +1,[ г., гкм],

Ф+= I-1[ г*?, г. + П] +1[ г. + П, г.+,],

Ф-= I,-,[г^1,г. -П] +1,[г. -П,г;.+,].

(14)

Искомую величину г 1 + полагаем равной

. . .

тому из трех чисел г.. или г. + п , или ^ + п, которое удовлетворяет ограничениям (11) и которому отвечает наименьшая из сумм (14):

г* при Ф;<Ф+, Ф, <Фг

=

гк+к при Ф+<Ф,, ФГ<Ф". (15) 2] - к при Ф7 < Фг, Ф~ < Ф+

Геометрическая схема метода локальных вариаций представлена на рис. 2.

Для оптимизации вычислений выбор шага дискретизации Ах и шага варьирования п проводился в соответствии со схемой двойного пересчета (т. е. для каждой итерации осуществлялась итерационная обработка с подбором оптимальных значений Ах и п). Данный метод позволяет находить локальный минимум функционала при выполнении условий: Ах ^ 0 П/(Ах)р при р > 2.

Следует отметить, что численное решение задачи (10), (11) может быть найдено и с использованием других известных методов решения задач в вариационной постановке, в частности, с помощью прямого численного решения краевой задачи для обыкновенного дифференциального уравнения, полученного с помощью уравнения Эйлера: г„ = - 3(1 + гх2) 2(гс - г).

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

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

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

Программная реализация модели проведена в ППП МаЙаЬ. На рис. 3 представлен результат численного моделирования динамики доменной границы при достижении равновесного положения.

Оценивая результат в нормированных величинах, можно заключить, что в поле постоянного градиента температуры максимальная величина отклонения доменной стенки от начального положения, соответствующего некоторому значению температуры Т0, составляет Дг = 0,1 усл. ед. при а = 0,02 усл. ед. и при линейном размере домена ДХ = 5 усл. ед. (Дг ~2 %). При увеличении значения градиента будет увеличиваться и максимальное отклонение от равновесного положения. Численное значение величины Дг будет определено параметрами моделирования Т0, Тс, г, а, X и не будет в итоге зависеть от величин, образующих комплекс ст0.

Алгоритм программной реализации

модельной задачи расчета компоненты тока поляризации

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

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

СО

I ($) = ,

(16)

где 5" - площадь грани кристалла, м2; Р5 - значение вектора спонтанной поляризации, Кл/м2; СО

- - изменение отношения переключенного

сИ

объема ДУш к полному объему кристаллического образца У, 1/с.

Графическая интерпретация изменения объема переполяризованной области за счет движения границы для кристалла со 180°-доменными стенками показана на рис. 1. В представленной модели, с учетом ее геометрии, соотношение (16) преобразуется к следующему виду:

I(Г) = 2• Р5= 2Р ^

(17)

У Д 5 Дt ' где - изменение площади грани кристалла, соответствующее изменению объема домена за счет смещения стенки.

Значение вектора спонтанной поляризации Р(Т) кристалла в точке с температурой Т может быть вычислено согласно соотношению (3). Тогда в нормированных единицах получим:

р»

Р(Т )= Р

4(Тс -Т0)

4Т-Т),

(18)

где Р3 - значение спонтанной поляризации сег-нетоэлектрика при температуре, равной значению Т0.

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

2 Р

I(t.) = — • , Ps

J At фс-т0

dx •

1 (zk - 4т - т0-a (zk+o+

, 1 / k k+1\ Lt T7 a / k , k+1\

+ 2 (ZN - ZN MTc - T0 - - (ZN + ZN )

l(

a

(zk - zk+1). т - To--(zk+zk+1)

(19)

где At - шаг по времени, определяемый номером итерации J.

Геометрическая визуализация формы компоненты поляризационного тока, обусловленного движением доменной границы под действием теплового поля, представлена на рис. 4.

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

СПИСОК ЛИТЕРАТУРЫ

1. Sogr, A.A. Advanced modes of imaging of ferroelectric domains in the SEM [Text] / A.A. Sogr, A.G. Maslovskaya, I.B. Kopylova // Ferroelectrics. -2006. -Vol. 341. -P. 29-37.

2. Maslovskaya, A.G. Simulation of ferroelectric domain structure imaging in pyroelectric mode by scanning electron microscopy [Text] / A.G. Maslovskaya // Ferroelectrics. -2010. -Vol. 398. -P. 55-63.

Рис. 4. Компонента поляризационного тока значения параметров моделирования: Р0 = 1, Т0 = -1, Тс = 0

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

3. Алешин, В.И. Моделирование переполяризации кристалла и керамики типа ВаТЮ3 [Текст] / В.И. Алешин, А.Г. Лучанинов // Изв. АН. Сер.физ. -2001. -Т. 65. -№ 8. -С. 1114-1118.

4. Нестеров, В.Н. Моделирование на ЭВМ ла-пласового давления и его роль в некоторых процессах перестройки доменной структуры сегнетоэлектри-ков [Текст] / В.Н. Нестеров, А В. Шильников. -XIV

Всерос. конф. по физике сегнетоэлектриков, 1995. -С. 340.

5. Parlunski, K. Domain pattern formation near phase transition - challenge for computer simulations [Text] / K. Parlunski // Ferroelectrics. -1997. -P. 245-253.

6. Донцова, Л.И. Плотность поверхностной энергии и спонтанное движение доменных стенок в кристаллах ТГС [Текст] / Л.И. Донцова, Э.С. Попов // Изв. АН СССР. Сер. физ. -1975. -Т. 39. -С. 854-856.

7. Нечаев, В.Н. Макроскопическая динамика доменных и межфазных границ в сегнетоэлектриках и сегнетоэластиках-сегнетоэлектриках [Текст] / В.Н. Нечаев, А.М. Рощупкин // Изв. РАН. Сер. Физ.

1995. -Т.59. -№ 9. -С. 11-25.

8. Лайнс, М. Сегнетоэлектрики и родственные им материалы [Текст] / М. Лайнс, А. Гласс. -М.: Мир, 1981. -736 с.

9. Черноусько, Ф.Л. Вариационные задачи механики и управления [Текст] / Ф.Л. Черноусько, Н.В. Ба-ничук. -М.: Наука, 1973. -240 с.

10. Масловская, А.Г. Упорядоченное движение доменных границ сегнетоэлектрических кристаллов в неоднородном тепловом поле [Текст] / А.Г. Масловская, Е.А. Ванина // Научно-технические ведомости СПбГПУ Сер. Физико-математические науки. -2010. -№ 2 (98). -С. 7-11.