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

УДК 53.072; 537.226.4; 548.537.611

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

© 2011 А.А. Лукичев, В.В. Ильина

Институт геологии и природопользования ДВО РАН, г. Благовещенск

Поступила в редакцию 10.03.2011

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

Явление гистерезиса хорошо известно в физике, технике, экономике и других отраслях науки. В физике гистерезисные зависимости встречаются у различных материалов с нелинейным откликом на внешнее воздействие, в частности у ферромагнетиков и нелинейных диэлектриков.

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

Поскольку практический метод получения петли Сойера-Тауэра [4] прост, то можно пред-Лукичёв Александр Александрович, кандидат физико-математических наук, старший научный сотрудник лаборатории конструкционной керамики. E-mail: lukichevaa@mail.ru.

Ильина Виктория Владимировна, кандидат физико-математических наук, младший научный сотрудник лаборатории конструкционной керамики. E-mail: nastja-5051974@mail.ru

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

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

Б(^а) = /и(ю,Н)/и0 Н(0 , (1)

для диэлектриков используется аналогичное выражение:

Р(г,ю) = а(с, Е) ■ Е(г), (2)

здесь ) - магнитная индукция, ¡л0 - магнитная постоянная, /л(о,И)- магнитная проницаемость, Н(Ь)- напряженность магнитного поля, Р(Ь)- вектор поляризации, а(о,Е) - поляризуемость, Е(Ь)-напряженность электрического поля.

Далее будем рассматривать некий абстрактный материал с функцией отклика вида (1), (2): У^,с) = Л(с, X) X^), (3)

где

X (Г) = X 0е (4)

- внешнее гармоническое воздействие, Л(а,Х) -коэффициент, зависящий от частоты изменения внешнего поля сои амплитуды внешнего воздействия Х(Ь), или спектральная функция. С математической точки зрения равенство (3) представляет собой параметрическую зависимость, где t является параметром. Как следует из (3), форма петли гистерезиса определяется видом функции Л(а),Х).

Сначала рассмотрим линейный гистерезис. В этом случае, коэффициент Л(а,Х) не зависит от

амплитуды внешнего воздействия и зависит только от частоты.

Известно, что форма спектров, как магнитной индукции, так и поляризации, во многих случаях близки к спектральной функции линейного осциллятора. Для резонансного режима колебаний [5] это функция Лоренца:

A(<) = 2 2

ю0 -ю + 2¡Ра> '

(5)

для заторможенного режима- укороченная спектральная функция:

A(W):

A)

1

A)

ю)+2ip<j) ю02 1+i<2T2

)

(6)

Y(ю,t) = |A(<)| elcp X(t) = |A(<)|X0

-i(<t + ф)

Для того, чтобы проанализировать зависимость формы петли от частоты и фазы построим график гистерезиса в простейшем виде. Перейдём к мгновенным значениям функции (4): X(t) = X0 sin <о t, аналогично, функция (8) будет: Y(t) = |A(<)|X0 sin(< t + ф) = Y0 sin(< t + ф). (9)

Коэффициент Yg=\A(w)\Xg для простоты пока будем считать постоянным. Петли гистерезиса для функции (9) в координатах (X(t),Y(t)) построены на рис. 2. Как видно из рисунка, петля линейного гистерезиса представляет собой эллипс, ширина и наклон главной оси которого существенно зависит от фазы. Угловой коэффициент главной оси определятся приближенным выражением:

где т = 2в / ®02 - постоянная времени, A0 - значение функции при ®=0. Аналог выражения (6) в физике диэлектриков известен как функция Дебая [6]. Графики модуля функций (5) и (6) в зависимости от частоты показаны на рис. 1а, зависимости фазы колебаний от частоты для этих же функций показаны на рис. 1б. Для описания заторможенных спектров недебаевского вида также применяются эмпирические спектральные функций Коула-Коула и Дэвидсона-Коула [6].

Спектральные функции Лоренца и Дебая могут быть представлены в виде

Л(®) = \Л(ш)\(7) где ] - зависящий от частоты фазовый сдвиг отклика материала на воздействие внешнего поля. Будем считать, что реальные спектры с достаточной точностью могут быть описаны функцией (7).

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

Y

п/2-ф

к =

ю

A(<)

X

п / 2-ф Л sin(n/2 -ф) (10)

ю

здесь в скобках указано значение аргумента t функций Xи Y. При нулевой фазе эллипс вырождается в прямую, с угловым коэффициентом k=Y(/Xg. При р=п/2 главная ось вертикальна, ширина эллипса максимальна, длины осей равны Y( и Х(. При }>п/2 угловой коэффициент меняет знак, и петли с фазой п-убудут зеркальным отражением петель с фазой р.

Учтем частотную зависимость амплитуды колебаний. Построим графики для линейного гистерезиса, используя спектральную функцию Лоренца (5), в резонансном режиме колебаний (Ь=(,1ф(), рис. 3а, и для заторможенного режима (Ь=5ю(), рис. 3б. В резонансном режиме зависимость формы петли от частоты и фазы рез-

1 А(ю),

о.е.

— 0.1

1-0.2

" " - ■ .

3.0 -1—

ол

1.0

*s fil.

—ТЕ-

Рис. 1. а - зависимость от частоты модуля спектральной функции Лоренца (5); б - Зависимость от частоты фазы функции Лоренца. Цифрами на графиках обозначена величина относительного затухания в/а. Пунктирная линия - граница между резонансным и заторможенным режимами колебаний. При затухании выше 3 графики функция Лоренца полностью совпадает с функцией Дебая (6)

0

б

а

Рис. 2. Графики для линейного гистерезиса, функция (9). Цифрами на графиках обозначена величина фазового сдвига

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

1 рез 1 1

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

но частоты резонанса, это следствие асимметрии спектральной функции. При малых затуханиях (Ь<0,1а>д) эти различия будут невелики, при больших затуханиях (Ь>0,1а>0), возникает значительная асимметрия. Также следует учитывать, что а>реззаметно отличается от сод, и срез<п/2 [5], соответственно, амплитуда вертикального эллипса (д>=п/2) будет меньше амплитуды наклонной петли соответствующей с=юрез (рис. 3а).

Если У0>>Х0, что обычно выполняется для нелинейных материалов, то можно считать, что при малых значениях отклика У@,а) боковая образующая петли - прямая, и ее наклон совпадает с наклоном главной оси и определяется выражением (10).

В отличие от резонансного режима, для заторможенных колебаний (Ь>юд) и амплитуда, и фаза меняются плавно в широком диапазоне частот. Как видно из рис. 3б, с ростом частоты петля плавно уширяется, её амплитуда снижается. При нулевой фазе (с0) эллипс вырождается в прямую угловым коэффициентом к=У/Хд, это значение углового коэффициента максимальное. Симметричная петля (д>=п/2) имеет очень маленькую амплитуду.

Как следует из рис. 1б, для резонансного режима колебаний фаза меняется достаточно резко от 0 до -п, в узкой области частот вблизи резонанса. При больших затуханиях (заторможенный режим) фаза будет ненулевой практически во всем интервале рабочих частот, диапазон изменения фазы от 0 до -п/2. Отсюда следует, что

а б

Рис. 3. а - линейный гистерезис для резонансного режима колебаний, спектральная функция (5), в/а =0,1 б - линейный гистерезис для заторможенного режима колебаний, спектральная функция (5), в/а0=5 Цифрами на графиках обозначена величина фазового сдвига Петля с ^=84,2о на рисунке а соответствует максимуму модуля спектральной функции (5)

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

Очевидно, что линейный гистерезис во многом подобен фигурам Лиссажу.

Рассмотрим характерные точки петли и найдем их связь с параметрами входного воздействия и отклика материала. Из рисунков видно, что ширина и наклон петли существенно зависят от фазового сдвига и отношения У0/Х. Найдём ширину петли по точкам пересечения оси ОХ, из равенства: = У0(ю)sin(юt + р) = 0.

Отсюда <=-р+пп, подставляем этот аргумент в функцию внешнего воздействия и получаем значение коэрцитивной силы:

АУ = X>п( р). (11)

Соответственно, ширина петли будет равна 2АХ. Таким образом, ширина петли определяется только величиной входного воздействия и фазой отклика.

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

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

L(a) = ctha - — a

(12)

Для диэлектриков величина a = p0E / Ш, где р0 - элементарный дипольный момент частиц, участвующих в процессе поляризации, к -постоянная Больцмана, Т - температура. В теории магнетизма a = m0Н ^ /Ш , где т0 - элементарный магнитный момент, Ие- эффективное поле [8]. График функции Ланжевена показан на рис. 4.

В теории Ланжевена рассматривается магнитное поле, действующее на отдельно взятый домен внутри материала. Но внутри материала действует эффективное поле, которое является

суперпозицией внешнего и внутреннего (реактивного) полей [8].

^ Fef = Fr + Ft, г (13)

где Fr - реактивное поле материала, Fi - внешнее поле. Реактивное поле появляется как отклик материала на внешнее воздействие и, следовательно, зависит от частоты. В магнитных материалах это молекулярное поле Вейсса, в диэлектриках - локальное поле Лоренца. Для нелинейных материалов, как правило, Fr >> Fi, следовательно Fef « Fr. Очевидно, что реактивное поле пропорционально отклику материала на внешнее воздействие вида (9), т.е.

Fr = F0|Л(<)|sin(<t + р), (14) где ¥д - постоянный коэффициент. Отсюда

то

a =

kT

F0|A(®)|sin(®t + ф), (15)

здесь, m0 - элементарный момент, F - амплитудный коэффициент, определяемый свойством материала.

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

S(X) = sign(Fr(t))(l - e-Fr('}|/, (16)

где £ - некоторый коэффициент, определяемый свойствами материала. График этой функции приведён на рис. 4. Функция (16) также зависит от реактивного поля вида (14). Из рисунка видно, что различие функций L(a) и S(X) невелики.

Теперь рассмотрим вопрос о зависимости амплитуды функции насыщения от частоты. Максимальная амплитуда приведённых выше

Рис. 4. Зависимость функции Ланжевена Ь(а) (12) и экспоненциальной функции насыщения Б(¥) (16) от величины реактивного поля

Рис. 5. Гистерезисная зависимость функции Ланжевена от величины внешнего воздействия.

Цифрами на графиках обозначена величина фазового сдвига

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

Ла(Г,®,X) = У,|Л(®)|, (17)

где Бо -обобщённая функция насыщения, У -амплитудный множитель. Поскольку амплитуда нелинейного гистерезиса будет заведомо меньше амплитуды линейного гистерезиса, здесь введён амплитудный множитель У, равный отношению амплитуд нелинейного и линейного процессов. Очевидно, что произведение У|А(й>)| равно насыщенному значению функции (17). Частотная зависимость реактивного поля определяет форму петли, спектральный множитель влияет на её амплитуду. Построим петлю гистерезиса с помощью равенства (17). В качестве функции насыщения используем функцию Ланжевена. Графики функции (17) показаны на рис. 5.

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

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

Найдём характерные точки нелинейной петли. При малых амплитудах функцию насыщения можно считать линейной, следовательно, ширина петли определяется так же, как и для линейного гистерезиса (11), соответственно коэрцитивная сила равна полуширине петли: Xc = X0 sin р . Найдём значение функции (17) при X(t)=0, для магнитных материалов эта величина называется остаточной намагниченностью. Из условия: rat=0 находим значение реактивного поля для этой точки: Fr0 = F0|^(®)|sin( р), и далее находим соответ-вующее ему значение функции насыщения

S=f(FJ-

Определим наклон боковых образующих нелинейной петли. Разложим функцию Ланжевена в ряд на линейном участке в окрестности нуля L(a))=a/3 [8], откуда, используя (14) и (15) получаем:

Y (t' ®) = YT Fef = YT FoHHsin(® t + р). (19)

Далее, по аналогии с (10), находим угловой коэффициент:

k = YsF0 m0

1

X0 3kT sin(^ /2 - р)

(20)

Если также разложить в ряд экспоненциальную функцию насыщения (16), то полученный угловой коэффициент будет совпадать с (20) с точностью до постоянных множителей. Из равенства (20) следует, что по наклону линейной петли гистерезиса можно определить амплитуду напряженности внутреннего поля ¥.

Для того, чтобы показать применимость разработанной выше модели для описания гистере-зисных зависимостей реальных материалов авторами были сняты гистерезисные зависимости для феррита НМ2000 и трансформаторного железа. Феррит имеет выраженную резонансную зависимость магнитной проницаемости от частоты с резонансной частотой /И/400 кГц. В линейном режиме гистерезисные зависимости для феррита точно повторяли приведённые на рис. 3а. В нелинейном режиме петли имели форму, близкую к приведённым на рис. 5. Для трансформаторного железа, имеющего АЧХ дебаевского вида, как линейные, так и нелинейные петли в основном подчинялись общим закономерностям, показанным на рис. 3б. Отсюда следует, что полученная модель позволяет описать поведение реальных материалов и найти связь формы петли и свойств материала.

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

1. М.А. Красносельский, А.В. Покровский. Системы с гистерезисом. М.: Наука, 1983. 275 с.

2. Hauser H., Fulmek P.L., Grossinger R. Hysteresis modeling and measurement for two-dimensional particle assembles//J. of magnetism and magnetic materials, Vol. 242-245, 2002. Р. 1067-1069.

3. Bottauscio O., Chiampi M., Chiarabaglio D., Repetto M. Preisach-type hysteresis models in magnetic field computation.//Physica B, Vol. 275, 2000.- p.p.24-39.

4. Барфут Дж., Тейлор Дж. Полярные диэлектрики и

их применение. М.: Мир, 1981. 526 с.

5. Ильина В.В., Лукичев А.А. Различные режимы вынужденных колебаний линейного осциллятора с затуханием и исследование соответствующих спектральных функций //Известия Самарского научного центра РАН. 2008. Т. 10. №3. C. 782-790.

6. Поплавко Ю.М. Физика диэлектриков. Киев, Вища школа, 1980. 400 с.

7. Кузнецов А.П., Кузнецов С.П., Рыскин Н.М. Нелинейные колебания. М.: Физматлит, 2002. 292 с.

8. Тамм И.Е. Основы теории электричества. Москва, Наука, 1976. 616 с.

SIMPLE MATHEMATICAL MODEL OF HISTERESIS LOOP FOR NONLINEAR MATERIALS

© 2011 A.A. Lukichev, V.V. Iljina

Institute of Geology and Nature Management of RAS, Blagoveschensk

The paper is devoted to development of simple and easily interpreted mathematical model of a hysteresis loop for nonlinear materials. The analytical expressions are received, allowing to model a hysteresis at various regimes of vibrations and with various functions of saturation. Critical points of a loop are defined and their dependence on properties of material is found. Applicability of the developed model to the description hysteresis dependences of real materials is shown. Keywords: histeresis loop, math simulation, nonlinear dielectrics, ferromagnetic,.

Alexander Lukichev, Candidate of Physics and Mathematics, Senior Staff Scientist of construction ceramics laboratory. E-mail: likichevaa@mail.ru.

Victoria Iljina, Candidate of Physics and Mathematics, Junior Staff Scientist of Construction ceramics laboratory. E-mail: nastja-5051974@mail.ru