CC BY
21Математика к Математическое
моделирование
И Ир://та! hfnelpLtb.ru 1531М 2412-5911
Ссылка на статью:
// Математика и математическое моделирование. 2018. № 05. С. 1-16
Б01: 10.24108/шаШш.0518.0000145
Представлена в редакцию: 26.09.2018
© НП «НЕИКОН»
УДК 539.3; 539.5
Конечно-элементное моделирование необратимого процесса поляризации сегнетоэлектрических керамик
Скалиух А.С.
а -Е-ька1длк11@ дтаД.сот
1 Институт математики, механики и компьютерных наук им. И.И.Воровича Южного федерального университета,
Ростов-на-Дону, Россия
В работе представлена конечно-элементная модель, описывающая необратимые процессы деформирования и поляризации. С использованием инкрементальной теории, принципа возможной работы и техники метода конечных элементов, строится система линейных алгебраических уравнений относительно приращений узловых значений основных переменных задачи: вектора перемещений и электрического потенциала. Матрица и правые части системы зависят от остаточных параметров и определяются в каждом равновесном состоянии. Как следствие, нелинейная задачи описывается последовательностью линейных задач. Наряду с геометрической информацией о конечных элементах, внешних нагрузках и свойствах тела, приходится дополнительно хранить информацию об элементных значениях остаточных параметров. Разработанная модель имплантирована в конечно-элементный комплекс, позволяющий находить поля остаточных деформаций и поляризации, физические характеристики частично поляризованного тела и его локальную анизотропию.
Ключевые слова: конечные элементы, итерационный процесс, остаточная деформация, остаточная поляризация
1. Введение
Необратимые процессы в поликристаллических сегнетоэлектрических средах, или керамиках, появляются, когда внешние механические или электрические поля достигают пороговых значений, например, при предварительной поляризации. Математические модели таких процессов позволяют оценить внутренние поля остаточной поляризации и деформации, определить локальную анизотропию, найти физические модули материала и физические характеристики конструкции в целом. Сложность задач связана в основном с поликристаллическим строением материала и доменной структурой сегнетоэлектрика, а также нелинейными зависимостями между механическими и электрическими параметрами, которые представляют собой диэлектрические и деформационные петли гистерезиса [1, 2], [3]. Как правило, модели строятся в квазистатической постановке с использованием
следующих гипотез и допущений: рассматривается безмоментная, геометрически линейная теория; процессы необратимы и описываются гистерезисными зависимостями; наблюдается фазовый переход "тело - тело", при котором меняется класс изотропии; значения механических напряжений и электрических полей не приводят к разрушению или электрическому пробою. Основными уравнениями являются полевые уравнения механики сплошных сред и электростатики диэлектриков, геометрические соотношения для деформаций и электрического поля, и определяющие соотношения гистерезисного типа для сопряженных электромеханических полей. Наиболее сложной частью являются определяющие соотношения, поскольку они требуют «настройки» модели на конкретный тип материала, что является сложной задачей. В отличие от линейных задач, где достаточно знать только упругие, пьезоэлектрические и диэлектрические постоянные материала, здесь необходимо определять функциональные зависимости, которые при возрастающих нагрузках одни, а при убывающих - другие. Подобные задачи возникают в явлениях пластичности, магнетизма, концентрации газа в пористых телах и т.д. Поэтому моделирование гис-терезисных зависимостей в таких случаях опирается на некоторые общие принципы и закономерности. Тонкости моделирования одномерных задач можно найти в [4]. Все модели, описывающие гистерезисные зависимости в сегнетокерамических материалах, можно условно разделить на феноменологические [5, 6, 7, 6, 9] и микромеханические [10, 11, 12, 13, 14, 15, 16]. Отдельно отметим подход [17], в котором используются статистические методы для описания необратимых частей. Однако для всех моделей общим является положение о представлении тензора деформаций и вектора поляризации в виде обратимых (упругих) и необратимых (пластических) составляющих. Отличительной особенностью является конкретное представление необратимых частей, которые в феноменологическом подходе формируются для представительного объема в целом, а в микромеханических моделях вначале, на микроуровне, определяется распределение всех векторов спонтанной поляризации, а затем путем усреднения находятся интегральные составляющие, которые присущи представительному объему. Тем не менее, в каждом из подходов для остаточных частей формулируются либо дифференциальные уравнения с производными по времени, либо уравнения в дифференциалах. Физические характеристики материала зависят от остаточных параметров, поэтому их определение также необходимо и для конкретизации определяющих соотношения для обратимых составляющих.
Задачей данного исследования является построение конечно-элементной модели для определения механических и электрических искомых параметров в необратимом процессе деформирования и поляризации поликристаллического сегнетоэлектрического тела произвольного объема. Чтобы значительно не увеличивать объем работы здесь не рассматриваются вспомогательные модели по выводу уравнений в дифференциалах для определения остаточной деформации и остаточной поляризации, а приводятся лишь ссылки на соответствующие работы. Рассматриваемая модель имплантирована в конечно-элементный комплекс ЛСБЬЛК, а результаты ее работы представлены в виде иллюстративного материала.
2. Постановка задачи
Рассматривается тело из термически деполяризованной керамики, к части поверхности которого приложены нагрузки силового характера. На другой части поверхности заданы условия кинематического характера. Кроме механических имеются и электрические граничные условия, выражающиеся в том, что на некоторой части поверхности нанесены электроды, соединенные с генератором напряжений, и имеется часть поверхности без электродов, как показано на рис. 1.
Рис. 1. Общий случай нагружения и Рис. 2. Структура материала: а) тело;
закрепления тела Ь) частица; с) домены; ф ячейка
Внешние механические и электрические нагрузки медленно меняются во времени, достигая таких значений, при которых в теле начинаются необратимые процессы. Требуется построить математическую модель, описывающую необратимый процесс, связанный с изменением внутренней структуры.
Под представительным объемом понимаем малый по сравнению с общим телом объем, содержащий в себе огромное количество кристаллитов, в которых имеется множество доменов [3], как показано на рис. 2. Каждый домен характеризуется вектором спонтанной поляризации р и тензором спонтанной деформации г8. Внутреннее строение сегнето-электриков таково, что имеются механизмы закрепления доменных стенок, препятствующие свободному переключению доменов из одного положения в другое. Поэтому механические напряжения и электрическое поле небольшой интенсивности лишь деформируют стенки доменов, вызывая обратимые параметры деформации и поляризации Ре, ге, являющиеся таким образом, параметрами состояния. При увеличении интенсивности нагрузок, по достижении ими пороговых значений, механизмы закрепления доменных стенок разрушаются, домены меняют свою ориентацию. Это проявляется в том, что меняются направления главных осей тензоров спонтанной деформации и направления векторов спонтанной поляризации, приводя к появлению остаточной (пластической) деформации и остаточной поляризации, которые определяются посредством усреднения по представительному объему спонтанных величин:
^ N 1 N
£о =—Е (£. )*' ро= ^ Е (р.)к •
N к=1 N к=1
Эти параметры являются параметрами процесса, они будут рассматриваться в дальнейшем, как необратимые параметры. Очевидно, что для случая термически деполяризованной керамики векторы спонтанной поляризации ориентированы хаотически в представительном объеме, их огромное количество, поэтому остаточные параметры такого состояния равны нулю. В необратимых процессах полная поляризация и деформация складываются из обратимых и остаточных частей: Р = Ре + Р0, £ = ее + £0, но определяющие соотношения формулируются для них раздельно. Кроме того, в уравнения Максвелла входит вектор электрической индукции, который определяется следующим образом Б = £оЕ + Ре + Р0.
С макроскопической точки зрения математическая постановка задачи включает в себя следующие уравнения [3]:
- полевые уравнения
V-о + р Г = 0; V-Б = де, (1)
где о, Б - тензор напряжений и вектор электрической индукции, р, Г, де - массовая плотность, плотность массовых сил и плотность заряда соответственно, V - набла оператор Гамильтона;
- геометрические соотношения
£ = (Уи + Vu т)/2, Е = -Уф , (2)
где £, Е - тензор деформаций и вектор напряженности электрического поля, и, ф - вектор перемещений и электрический потенциал соответственно;
- определяющие соотношения для обратимых составляющих [18]
£ - £ 0 = 8(£ о): о + а т(Ро) - Е,
Б - Ро = а(Ро): о + э(£о) - Е,
где 8(£0), а(Р0), э(£0) - тензор упругих податливостей, тензор пьезоэлектрических констант, тензор диэлектрических проницаемостей; эти физические характеристики материала зависят линейным образом от остаточных параметров
в(£ о) = §о + К : £ о, а(Ро) = N - Ро, э(£ о) = э о + М : £ о ; - определяющие соотношения для необратимых составляющих [3]
dP0 = т(Рш- Р0 )| <Ее |,
<£ о = п(£ £ о) 1 <о £- <о е(I где Р = Р (Р0, £0), £ш = £ш (Р0, £0) векторная и тензорная функции предельных значений поляризации и деформации, выраженные в виде поверхностных интегралов по единичной сфере, которые в силу громоздкости здесь не выписываются; Ее = Е + а Р0, ое = о + Р£0 - внутренние эффективные поля Вейсса; <о е{и, йо е -главные значения приращения тензора эффективных напряжений; и, наконец, т, п, а, р -известные постоянные.
Замыкают постановку задачи граничные условия на S, где имеет место два независимых разбиения: £ = ^ = ^ :
п • а = р* (^ на ; и = и* (^ на £и; п • Б = д* на ; ф = ф*(/) на £, (5)
и начальные условия для необратимых составляющих неполяризованной керамики:
е„и = 0, Ро !г=о = 0 • (6)
Требуется определить решения уравнений (1)-(4), удовлетворяющие условиям (6),
(7).
3. Метод решения
В квазистатическом процессе граничные условия представлены последовательностью значений: (р*, и*, ф} е (р*г., и*, фг}: р* = р* ), и*г. = и* (¿г.), ф*г. = ф* (¿г.), I = 1,..., п, вследствие чего имеем последовательность бесконечно близких равновесных состояний С(0), С(1), ...,С(г-1), С(г),..., С(п), где С(0), С(г) и С(п) - состояния начальное, текущее и конечное соответственно. Начальные условия (6) задаются в состоянии С(0), и все параметры, однозначно определяющие систему, считаются известными. Задача заключается в том, чтобы по известным параметрам предыдущего состояния С(г) найти все параметры состояния С(г+1).
Для решения воспользуемся принципом возможных перемещений
| (а : ¿г -р? • ¿и - Б • ¿Е + де5 ф)^а - | р • - | д^фйБ = 0, (7)
I* ииии I ^
а ^
где ¿и = 0 на £м, ¿ф = 0 на £ .
При переходе от С(г) к С(г+1) состоянию все полевые характеристики равно как и граничные функции главных и естественных условий получают приращения, поэтому в состоянии С(г+1) могут быть представлены в виде:
о(г) + Ло(г), г (г) + Лг (г), г (г) + Лг (г), Е(г) + ЛЕ(г), Ре(г) + ЛРе(г), р(г) + Лр(г), ?(г) + л? (г), и(г) + Ли(г), ф() + Лф(г), и*г) + Ли*г), р*г) +Лр*г), ф*г) + Лф*г). Величины со знаком « Л» играют роль малых приращений. Принцип возможных перемещений (7) применим для каждого равновесного состояния, в частности для С(г+1) состояния:
I[(о(г) + Ло(г)) : ¿(г(г) + Лг(г)) - (Б(г) + ЛБ(г)) • ¿(Е(г) + ЛЕ(г)) -а
р(?(г) + л?(г)) • ¿(и(г) + Ли(г))]^а - I (р(г) + Ар• ¿(и(г) + Ли(г))^£ = 0.
Варьированию подлежат приращения функции перемещения и электрического потенциала, т.е. ¿(и(г) + Ли(г)) = ¿Ли(г), ¿(ф(г) + Лф(г)) = ¿Лф(г), откуда получаем
¿(е(0 + Ае(г )) = ¿Ае(г), ¿(Е(г) + АЕ(г)) = ¿АЕ(г). Учитывая, что механические напряжения, электрическая индукция и функции внешних естественных граничных условий считаются неизменными в каждом состоянии, и отбрасывая выражение
| (о(0 : ¿(Ае() - Ае(г-1)) - Б(г) • ¿(АЕ(г) - АЕ(г-1)) -
Q
— pf(0 • S(Au(0 - Au('—1:i))dQ — Jpi° • S(Au(i) - Au('—1))dS,
содержащие члены высшего порядка малости, получим:
J(Ao(i) : SAs(i) — AD(i) • SAE(i) — pAf(i) • SAu(i) +Ag(°SАФ)/0 —
Q
— jAp• SAu(i)dS — JAq• 5^(0)dS = 0 (8)
Оценим, как изменяются определяющие соотношения для обратимых составляющих (3) при переходе от одного состояния к другому. С этой целью отметим, что неполяризо-ванная керамика является изотропным телом и не обладает пьезоэлектрическими свойствами, в то время как поляризованная до состояния насыщения имеет поляризацию I Pol= Psat, и главные значения деформации So: =-S sat /2, 8 0П =—S sat /2,, 8ош =8 sat, причем третья главная ось совпадает с направлением вектора остаточной поляризации. В неполяризованном состоянии керамика представляет собой изотропное тело, тензоры упругих податливостей, пьезоэлектрических модулей и диэлектрических проницаемостей которой имеют значения:
S = S0, d = 0, э = э0,
а в состоянии поляризации до насыщения керамика относится к классу трансверсально-изотропного тела, у которой физические модули принимают значения:
S = S sat, d = d sat , э = э sat •
Введем поправочные тензоры
S1 = Ssat — S0, d1 = dsat , э1 = эsat — э0 •
Тогда можно построить линейные функции физических модулей, зависящих от остаточных параметров следующего вида [18]:
S(s 0) = S0 + ^ Si, d = 1^1 di, э1 = э 0 + ^ э1 (9)
8 sat P sat 8 sat
Соотношения (3) справедливы для любого равновесного состояния. Запишем их для состояний C('} и C('+1) и вычтем одно из другого, получим
As= S(s('+1)): о('+1) — S(s: о(i) + dT(P0(i+1)) • E('+1) — dT(P0(O) • E(i), AD(° = d(P0('+1)): о('+1) — d(P0(O): о(i) + a(s('+1)) • E('+1) — a(s• E(i). После несложных преобразований с учетом (9) получаем
S
S
S
D
Аг« =
Др(г)
: До(г) +
АБ(г) =
)) +
А I Ро(г) I
а т(Р0г)) +
ДI Ро(г) I йт
• ДЕ(г) +
Ав
(О
А I Р(г) I 0111 ^ : о(° ат • Е(г),
: Ао (° +
э(г 0)) +
Ав
а)
АI Ро(г) I
а : о(г) +
Ав
а)
• АЕ(г) +
-Э • Е(°.
В отличие от линейных соотношений (3) определяющие соотношения для приращений (10) имеют добавочные слагаемые, которые можно рассматривать как дополнительные приращения обратимой части деформации и поляризации вследствие перестройки структуры материала.
Примем следующие допущения.
1. Т.к. для бесконечно близких равновесных состояний выполняются оценки АI Р0(г) I << I Р0(г) I , I Ав(щ I << I в од I, то можно пренебречь слагаемыми в прямоугольных скобках.
2. Последние слагаемые являются малыми величинами, вследствие чего ими также можно пренебречь.
3. Чтобы приращения вектора остаточной поляризации и остаточной деформации при переходе от С() к С (г+:) состоянию не содержали неизвестных значений, заменим их приращением предыдущего состояния, т.е.
АI Ро(г) I ^ АI Ро('-1) I , АI в0щ I ^ АI В0Щ4 I.
В рамках сделанных допущений, определяющие соотношения для приращений обратимых составляющих принимают вид
Аг= 8(г: Ао(0 + ат(Р0(°) • АЕ(г),
(11)
АБ(г) = а(Р0(г)) : Ао(г) + э(г(,)) • АЕ(г).
Прежде, чем использовать эти соотношения, необходимо определить тензор остаточных деформаций и вектор остаточной поляризации для С0) состояния. Для этого воспользуемся уравнениями (4). Заметим, что эти уравнения содержат приращение остаточных параметров как в левой, так и в правой частях за счет эффективных полей, что осложняет задачу. Поэтому для определения приращений остаточных параметров предлагается использовать метод последовательных приближений. Дифференциалы заменяются на приращения, после чего итерационный процесс можно представить следующими соотно-
шениями:
(Аг0-1)), = п(гГ - г0-1)),^Аа/0-1)),- - (Аа/ (г 0°)* =г0г-1) + (Аг 0-1)),, (АРо(г-1)), = - Р0г-1)),-1 I (АЕе/('-1)),-1 I ,
(Ро(г)), =Р<'-1) + (АР0г-1)),,
(12)
в
в
э
в
в
Очевидно, что итерационный процесс возможен только для монотонных функций нагрузок. В силу этого для первого шага итерации принимается
(АР0( '-1) )0 = 0, (АеГ )0 = 0, и, кроме этого (Р0( г-1))^ = Р0('-1), (е^-4)^ = е('-1). На практике обычно число итераций не превышает десяти для достижения приемлемой точности.
В практических задачах удобно использовать матричное представление Фойхта, когда тензорам напряжений и деформаций ставятся в соответствие векторы:
0 ^ 0 = {оп, а22, азз, а2з, а^ а^^
е ^ е = {811, 822, 833, 2823, 2813, 2812} ;
тензорам 8, d, э ставятся в соответствие матрицы 18, {!, э, элементы которых, определяются по элементам тензоров по следующему правилу
о , о _ Г и ' = ] а ^ о I9 - (' + ^ * ]
иЫ] ^ ика I- ■ ■
э ^э Н'- р+]),'
] ]
а векторы не меняются.
Тогда определяющие соотношения (11) можно записать в матричной форме:
Ае = 8(е : А0(0 + d т(Р0(О) • АЕ(0, АБ= d(Р0(°): А0(0 + э (е • АЕ(0, и преобразуем их к другой, более удобной для дальнейшего форме, решив относительно А0^ . Получим соотношения
А0« = С(е(0) • Ае« -ет(Р0('),е0')) • ДЕ('),
АБ е) = е (Ро(0, е 0')) • Ае «) + Э (Р00, е (')) • АЕ «
где введены обозначения
С = 8-1, е = с1 • 8-1, Э = э-с • 8-1: dт. В дальнейшем будет рассматриваться только матрично-векторные обозначения, поэтому обозначения «шапочка» сверху опускаем.
Выберем согласованную конечно-элементную сетку, задаваемую в = , аппроксимирующую область О. На каждом конечном элементе От неизвестные полевые функции и, ф аппроксимируем в форме:
и - N (х) • и, ф» N (х) • Ф , где N - матрица функций формы для поля перемещений и ; NТ - вектор-строка функций формы для поля электрического потенциала ф; и, Ф - локальные векторы соответствующие узловым степеням свободы, размерность которых зависит от выбранного типа конечного элемента. Приращение искомых функций находим по приращениям искомых узловых значений, которые в С ^ состоянии определяются следующим образом:
Ли(0 = МТ • ли(0, Лф = МТ • ЛФ(0,
и ' ф , (14)
Л£(0 = В • ли(0, ЛЕ(0 = -В0 • ЛФ(0
где матрицы В = Б • N и В0 = Бф • N, а Бх, Бф " матрицы дифференциальных операций, связывающие деформации с перемещениями и электрическое поле с потенциалом. Конкретный вид матриц форм и дифференциальных операций также связан с выбором вида конечного элемента.
Интегральное выражение (10) можно переписать в векторном виде (символ Л означает приращение, 5 вариацию):
£ |[(Ло(0)т • 5(Лг(0) - (ЛГ(0)т • 5(Ли(°) - (ЛБ(0)т -5(ЛЕ(0) +
т О
Лд()5(Лф(0)]^О - £ |(Лр(0)т • 5(Ли- £ |Лч.(05(Лф= 0. (15)
п Б, Ч
Суммирование в объемных интегралах распространяется на все конечные элементы, а в поверхностных только на грани тех конечных элементов, которые выходят на соответствующие границы. Основными неизвестными являются приращения вектора перемещения Ли(г) и электрического потенциала Лф(г).
Дальнейшие действия стандартные для метода конечных элементов. Подставляя (13), (14) в (15), и используя независимость соответствующих вариаций функций в объеме и на поверхности, получим локальную систему линейных алгебраических уравнений:
Кии -Ли« + Киф • ЛФ(0 = г« + г«,
КТф^Ли« - Кфф • ЛФ(г) = г2(г) + гфг),
(16)
где
Кии =| ВТ • С • В; ао, Киф = | ВТ • ет • В0 dО, Кфф=| Вт • Э • В0ао,
От От От
Т00 = |N • аР* )¿Б + |N • Лf(0ао, Г^ = |ВТ • c • Л£(;-1)ао,
г2(г) = |NфЛд*г)аз + |Nф Аао, Гффг) = |Вт • е• Л^¿О- |Вт • АР0(г-1)ао.
От От От
Проводя операцию сборки матрицы и правых частей по значениям локальных матриц и правых частей, получаем общую систему уравнений, после решения которой находятся неизвестные узловые значения приращений перемещений и электрического потенциала Ли(г), ЛФ(г).
4. Результаты численных экспериментов
Схематично алгоритм решения задачи можно описать следующим образом. Пусть в состоянии С(г) все параметры известны, т.е. ЛР0( г-1), Лг(г-1), Р0( г-1), г(г-1), а также
Ли(г -1), ЛФ(г-1).
Переход к следующему состоянию осуществляется заданием временного параметра t'+1, после чего определяются граничные функции р('+1), ф('+1) и находятся их приращения
Ар(0, Аф!0 (рассматривается наиболее часто встречающийся случай, когда = 0, д = 0).
Дальнейшие действия связаны с каждым конечным элементом, для которого находятся значения физических модулей С, е, Э, элементов локальных матриц и правых частей Кии, Киф, Кфф, Г^, Г^, Гсф'). Далее собирается общая матрица и правые части. Из решения системы находятся приращения Аим, АФМ, по которым определяются узловые значения и('+1) = и(0 +Аи(0, Ф('+1) = Ф(0 + АФ(0 . Далее определяются элементные значения электрического поля, полной деформации и упругой части деформации Е( '+1), е('+1), е('+1). Из определяющих соотношений находятся элементные значения механических напряжений о(г). И, наконец, определяются элементные значения приращения остаточных параметров АР^, Ае^, и сами остаточные параметры Р^, е(0 . Таким образом,
все параметры С('+1) состояния определены. Далее осуществляется переход С('+2) состоянию, или, если внешние нагрузки достигли конечных значений, осуществляется выход из итерационного цикла.
Данный алгоритм был имплантирован в конечно-элементный комплекс ЛСБЬЛК. В качестве примера приводится поле остаточной поляризации цилиндрического преобразователя с центральным и кольцевыми электродами, как показано на рис. 3. Включая в электрическую схему электроды цилиндра различным образом, можем сформировать совершенно разные поля остаточной поляризации. Например, на рис. 4 показана половина осевого сечения для случая, когда электрический потенциал на верхних электродах одного знака, а на нижнем электроде - другого.
Рис. 3. Общий вид керамического Рис. 4. Поле остаточной
диска с электродами поляризации в осевом сечении
По полю остаточной поляризации и деформации определяется локальная анизотропия материала и вычисляются физические характеристики материала. Эта информация необходима для дальнейшего анализа, например при проведении модального или гармонического анализа неоднородно поляризованных пьезокерамических тел.
5. Обсуждение результатов
Сложность моделирования необратимых процессов поляризации и деформирования поликристаллических сегнетоэлектрических материалов заключается в том, что наравне с упругими составляющими появляются остаточные характеристики, влияющие на структу-
ру материала и, как следствие, на сами упругие деформации и поляризацию. Для квазистатических процессов необратимый процесс можно рассматривать пошагово и разделить на части, в которых последовательно вначале находятся искомые характеристики в виде перемещений и электрического потенциала, а затем и приращение остаточных параметров. Использование техники конечных элементов позволяет нелинейную задачу свести к решению последовательности систем линейных алгебраических уравнений, у которых как матрица, так и правые части меняются, при переходе от одного состояния к другому. Для каждого конечного элемента, который рассматривается как представительный объем, необходимо определять необратимые составляющие. С этой целью приходится использовать элементные значения электрического поля и механических напряжений. При определении же основных неизвестных, таких как перемещения и электрический потенциал, используются узловые значения. Выбор той или иной вспомогательной модели, позволяющей определять остаточные параметры, является чисто субъективной проблемой, однако при конечно-элементном подходе приходится хранить информацию об остаточных параметрах на каждом конечном элементе. Поэтому желательно, чтобы такие модели были связаны с минимальной хранимой информацией. Например, для модели Джила-Атертона [17], достаточно хранить информацию только об остаточных параметрах, к которым относятся элементные значения компонент тензора деформации и вектора поляризации. Однозначная разрешимость возникающих линейных систем уравнений основана на том, что элементы матрицы, связанные с частично поляризованной керамикой, принимают промежуточные значения, между значениями деполяризованного состояния и состояния полного насыщения. И в том и в другом случае задача решается однозначно. В процессе поляризации материал керамики приобретает локальную анизотропию. Вычисление физических характеристик в глобальных осях в каждом состоянии осуществляется с помощью матриц перехода, которые однозначно определяются по хранимой информации об остаточных параметрах, но эти элементарные, громоздкие формулы здесь не приводятся. Описанный алгоритм был имплантирован в конечно-элементный комплекс ACELAN, позволяющий не только определять остаточные параметры и находить физические характеристики локально анизотропного материала, но и определять многие интересующие пользователей физические поля.
6. Заключение
Предложена конечно-элементная модель, описывающая необратимые процессы в поликристаллических сегнетоэлектрических средах. Необратимый процесс рассмотрен в квазистатической постановке для множества значений внешних нагрузок, в виде последовательности равновесных состояний, для каждого из которых определяются приращения искомых величин по приращениям граничных условий. Предложенная конечно-элементная техника позволяет рассматривать задачи как в трехмерной, так и в двумерной постановке путем выбора соответствующих конечных элементов. Нелинейность задачи заменяется решением последовательности линейных задач, для систем линейных алгеб-
раических соотношений, в которых учитываются изменяющиеся остаточные параметры. Предложенная техника имплантирована в конечно-элементный комплекс, позволяющий находить поля распределения остаточных величин, физические характеристики частично поляризованного тела и локальную анизотропию, возникающую вследствие необратимости процесса деформирования и поляризации.
Работа поддержана РФФИ, грант 17-08-00860-a.
Список литературы
1. Lynch C.S. The effect of uniaxial stress on the electro-mechanical response of 8/65/35 PLZT // Acta Materialia. 1996. Vol. 44. No. 10. Pp. 4137-4148. DOI: 10.1016/S1359-6454(96)00062-6
2. Dayu Zhou, Kamlah M. Dielectric and piezoelectric performance of soft PZT piezoceramics under simultaneous alternating electromechanical loading // J. of the European Ceramic Society. 2005. Vol. 25. No. 12. Pp. 2415-2420. DOI: 10.1016/j .jeurceramsoc.2005.03.073
3. Белоконь А.В., Скалиух А.С. Математическое моделирование необратимых процессов поляризации. М.: Физматлит, 2010. 328 с.
4. Skaliukh A.S. About mathematical models of irreversible polarization processes of a ferroelectric and ferroelastic polycrystals // Ferroelectrics and their applications / Ed. by Husein Irzaman. IntechOpen, 2018. Ch. 4. Pp. 39-69. DOI: 10.5772/intechopen.78262
5. Huber J.E., Fleck N.A. Multi-axial electrical switching of a ferroelectric: theory versus experiment // J. of the Mechanics and Physics of Solids. 2001. Vol. 49. No. 4. Pp. 785-811. DOI: 10.1016/S0022-5096(00)00052-1
6. Kamlah M., Bohle U. Finite element analysis of piezoceramic components taking into account ferroelectric hysteresis behavior // Intern. J. of Solids and Structures. 2001. Vol. 38. No. 4. Pp. 605-633. DOI: 10.1016/S0020-7683(00)00055-X
7. McMeeking R.M., Landis C.M. A phenomenological multi-axial constitutive law for switching in polycrystalline ferroelectric ceramics // Intern. J. of Engineering Science. 2002. Vol. 40. No. 14. Pp. 1553-1577. DOI: 10.1016/S0020-7225(02)00033 -2
8. Landis C.M. Fully coupled, multi-axial, symmetric constitutive laws for polycrystalline ferroelectric ceramics // J. of the Mechanics and Physics of Solids. 2002. Vol. 50. No. 1. Pp. 127-152. DOI: 10.1016/S0022-5096(01)00021 -7
9. Haug A., Knoblauch V., McMeeking R.M. Combined isotropic and kinematic hardening in phenomenological switching models for ferroelectric ceramics // Intern. J. of Engineering Science. 2003. Vol. 41. No.8. Pp. 887-901. DOI: 10.1016/S0020-7225(02)00320-8
10. Huber J.E., Fleck N.A. Ferroelectric switching: a micromechanics model versus measured behaviour // Eur. J. of Mechanics A/Solids. 2004. Vol. 23. No. 2. Pp. 203-217. DOI: 10.1016/j.euromechsol .2003.11.006
11. Kamlah M., Liskowsky A.C., McMeeking R.M., Balke H. Finite element simulation of a polycrystalline ferroelectric based on a multidomain single crystal switching model // Intern. J. of Solids and Structures. 2005. Vol. 42. No. 9-10. Pp. 2949-2964.
DOI: 10.1016/j .ijsolstr.2004.09.045
12. Haug A., Onck P.R., Van der Giessen E. Development of inter- and intragranular stresses during switching of ferroelectric polycrystals // Intern. J. of Solids and Structures. 2007. Vol. 44. No. 6. Pp. 2066-2078. DOI: 10.1016/j.ijsolstr.2006.07.024
13. Pane I., Fleck N.A., Chu D.P., Huber J.E. The influence of mechanical constraint upon the switching of a ferroelectric memory capacitor // Eur. J. of Mechanics A/Solids. 2009. Vol. 28. No. 2. Pp. 195-201. DOI: 10.1016/j.euromechsol.2008.09.002
14. Jayabal K., Menzel A., Arockiarajan A., Srinivasan S.M. Micromechanical modelling of switching phenomena in polycrystalline piezoceramics: application of a polygonal finite element approach // Computational Mechanics. 2011. Vol. 48. No. 4. Pp. 421-435.
DOI: 10.1007/s00466-011-0595-4
15. Daniel L., Hall D.A., Withers P.J. A multiscale model for reversible ferroelectric behaviour of polycrystalline ceramics // Mechanics of Materials. 2014. Vol. 71. Pp. 85-100.
DOI: 10.1016/j.mechmat.2014.01.006
16. Осипова Н.Г., Семенов А.С. Моделирование нелинейного поведения пьезокерамики и тетрагональной структуры методами конечно-элементной гомогенизации // Науч-техн. ведомости С.-Петербург. гос. политехн. ун-та. Физ.-матем. науки. 2011. № 4(134). С. 56-64.
17. Smith R.C., Zoubeida Ounaies. A domain wall model for hysteresis in piezoelectric materials // J. of Intelligent Material Systems and Structures. 2000. Vol. 11. No. 1. Pp. 62-79. DOI: 10.1106/HPHJ-UJ4D-E9D0-2MDY
18. Скалиух А.С. Моделирование необратимых процессов поляризации и деформации в керамике // 14-й Российско-СНГ-балтийско-японский симп. по сегнетоэлектричеству (С.-Петербург, Россия, 14-18 мая 2018): Тез. докл. СПб., 2018. С. 71.
Mathematics I Mathematical Modelling
f.Vjvf
rtiFJli' }OHT*TU1l
h tip:/ ArraihfnelpLtb.ru
ISSN 2412-591 i
Mathematics and Mathematical Modeling, 2018, no. 05, pp. 1-16.
DOI: 10.24108/mathm.0518.0000145
Received: 26.09.2018
© NP "NEICON"
Finite-element Modeling Irreversible Polarization Process of Ferroelectric Ceramics
A.S. Skaliukh1' a -5-5kaliukh:g gmail.com
1Institute of Mathematics, Mechanics and Computer Sciences n.a. I.I.Vorovich of the Southern Federal University,
Rostov on Don, Russia
Keywords: finite elements, iterative process, irreversible strain, irreversible polarization
A finite-element model developed for quasi-static processes describes irreversible processes of deformation and polarization occurring in polycrystalline ferroelectric media due to the effect of intense electric fields and mechanical loads. The paper presents external parameters such as strain and polarization as a sum of residual and reversible parts. Using the incremental theory, the virtual work law, and the constitutive relations for reversible and irreversible components, a system of linear algebraic equations was built for the increments of nodal values of the main variables, namely the displacement vector and the electric potential, during the transition from one equilibrium state to another.
The constructed constitutive relations connect the reversible parts of the strain and polarization with the stresses and the electric field in the form of linear tensor equations. It is shown that the physical characteristics depend on the residual parameters so that the coefficients of elastic compliance and dielectric constant linearly depend on the principal values of residual strain, and the piezoelectric modules depend linearly on the module of residual polarization. The constitutive relations for the increments of the residual parameters are determined as element values for each finite element from the equations in differentials. Ultimately, the task is reduced to a system of linear algebraic equations, the matrix and right sides of which depend on the residual parameters and are determined at each equilibrium state. As a result, the non-linearity of the problem is replaced by solving a sequence of linear problems until the external loads reach their final values.
The model is implanted into a finite-element complex, which allows us to determine the fields of residual strain and polarization, the physical characteristics of a partially polarized body, and local anisotropy for the case of complete and partial polarization.
References
1. Lynch C.S. The effect of uniaxial stress on the electro-mechanical response of 8/65/35 PLZT. Acta Materialia, 1996, vol. 44, no. 10, pp. 4137-4148. DOI: 10.1016/S1359-6454(96)00062-6
2. Dayu Zhou, Kamlah M. Dielectric and piezoelectric performance of soft PZT piezoceramics under simultaneous alternating electromechanical loading. J. of the European Ceramic Society,, 2005, vol. 25, no. 12, pp. 2415-2420. DOI: 10.1016/j.jeurceramsoc.2005.03.073
3. Belokon' A.V., Skaliukh A.S. Matematicheskoe modelirovanie neobratimykh protsessov poliarizatsii [Mathematical modeling of irreversible processes of polarization]. Moscow: Fizmatlit Publ., 2010. 328 p. (in Russian).
4. Skaliukh A.S. About mathematical models of irreversible polarization processes of a ferroelectric and ferroelastic polycrystals. Ferroelectrics and their applications / Ed. by Husein Irzaman. IntechOpen, 2018. Ch. 4. Pp. 39-69. DOI: 10.5772/intechopen.78262
5. Huber J.E., Fleck N.A. Multi-axial electrical switching of a ferroelectric: theory versus experiment. J. of the Mechanics and Physics of Solids, 2001, vol. 49, no. 4, pp. 785-811. DOI: 10.1016/S0022-5096(00)00052-1
6. Kamlah M., Bohle U. Finite element analysis of piezoceramic components taking into account ferroelectric hysteresis behavior. Intern. J. of Solids and Structures, 2001, vol. 38, no. 4, pp. 605-633. DOI: 10.1016/S0020-7683(00)00055-X
7. McMeeking R.M., Landis C.M. A phenomenological multi-axial constitutive law for switching in polycrystalline ferroelectric ceramics. Intern. J. of Engineering Science, 2002, vol. 40, no. 14, pp. 1553-1577. DOI: 10.1016/S0020-7225(02)00033 -2
8. Landis C.M. Fully coupled, multi-axial, symmetric constitutive laws for polycrystalline ferroelectric ceramics. J. of the Mechanics and Physics of Solids, 2002, vol. 50, no. 1,
pp. 127-152. DOI: 10.1016/S0022-5096(01)00021 -7
9. Haug A., Knoblauch V., McMeeking R.M. Combined isotropic and kinematic hardening in phenomenological switching models for ferroelectric ceramics. Intern. J. of Engineering Science, 2003, vol. 41, no. 8, pp. 887-901. DOI: 10.1016/S0020-7225(02)00320-8
10. Huber J.E., Fleck N.A. Ferroelectric switching: a micromechanics model versus measured behaviour. Eur. J. of Mechanics A/Solids, 2004, vol. 23, no. 2, pp. 203-217.
DOI: 10.1016/j.euromechsol .2003.11.006
11. Kamlah M., Liskowsky A.C., McMeeking R.M., Balke H. Finite element simulation of a polycrystalline ferroelectric based on a multidomain single crystal switching model. Intern. J. of Solids and Structures, 2005, vol. 42, no. 9-10, pp. 2949-2964.
DOI: 10.1016/j .ijsolstr.2004.09.045
12. Haug A., Onck P.R., Van der Giessen E. Development of inter- and intragranular stresses during switching of ferroelectric polycrystals. Intern. J. of Solids and Structures, 2007, vol. 44, no. 6, pp. 2066-2078. DOI: 10.1016/j.ijsolstr.2006.07.024
13. Pane I., Fleck N.A., Chu D.P., Huber J.E. The influence of mechanical constraint upon the switching of a ferroelectric memory capacitor. Eur. J. of Mechanics A/Solids, 2009, vol. 28, no. 2, pp. 195-201. DOI: 10.1016/j.euromechsol.2008.09.002
14. Jayabal K., Menzel A., Arockiarajan A., Srinivasan S.M. Micromechanical modelling of switching phenomena in polycrystalline piezoceramics: application of a polygonal finite element approach. Computational Mechanics, 2011, vol. 48, no. 4, pp. 421-435.
DOI: 10.1007/s00466-011-0595-4
15. Daniel L., Hall D.A., Withers P.J. A multiscale model for reversible ferroelectric behaviour of polycrystalline ceramics. Mechanics of Materials, 2014, vol. 71, pp. 85-100.
DOI: 10.1016/j.mechmat.2014.01.006
16. Osipova N.G., Semyonov A.S. The simulation of nonlinear behavior for ferroelectric ceramics with tetragonal structure by finite-element homogenizing. Nauchno-tekhnicheskie vedomosti S.-Peterburgskogo gosudarstvennogo politekhnicheskogo universiteta. Fiziko-matematicheskie nauki [St.-Peterburg Polytechnic Univ. J.: Physics and Mathematics], 2011, no. 4(134), pp. 56-64 (in Russian).
17. Smith R.C., Zoubeida Ounaies. A domain wall model for hysteresis in piezoelectric materials. J. of Intelligent Material Systems and Structures, 2000, vol. 11, no. 1, pp. 62-79.
DOI: 10.1106/HPHJ-UJ4D-E9D0-2MDY
18. Skaliukh A.S. Modelirovanie neobratimykh protsessov poliarizatsii i deformatsii v keramike[Modeling of the irreversible phenomena of polarization and deformation in ceramics]. 14 Rossijsko-SNG-baltijsko-iaponskij simpozium po segnetoelektrichestvu [14th Russia/CIS/Baltic/Japan Symp. on Ferroelectricity: RCBJSF-2018 (St. Petersburg, Russia, May 14-18, 2018)]: Abstracts. St. Petersburg, 2018. P. 71 (in Russian).
