Актуальность исследования и развития математических моделей гистерезиса и методов моделирования физических полей связана с широким использованием в технических устройствах элементов, выполненных из магнитных материалов и диэлектриков, имеющих неоднозначные зависимости характеристик. Современный уровень развития техники предполагает все более активное использование компонентов и материалов на пределе их возможностей, что требует широкого применения компьютерного моделирования при проектировании и совершенствовании математических моделей с целью повышения их соответствия исследуемым объектам [1, 2].
Целью исследования является усовершенствование инструментов моделирования магнитного поля с учетом гистерезиса. В рамках поставленной цели выполнено объединение в единый метод ранее разработанных модификаций модели Джилса – Атертона [3] и метода моделирования магнитного поля, построенного на основе граничных интегральных уравнений [4].
Материалы и методы исследования
Согласно модели Джилса – Атертона, в объеме ферромагнитного материала векторы магнитной индукции 
, напряженности магнитного поля 
и намагниченности 
коллинеарны, а их проекции на общую ось удовлетворяют соотношению
(1)
Связь M(H) задается дифференциальными соотношениями, основанными на безгистерезисной кривой Man(H), для описания которой была выбрана функция Ланжевена
где Ms – намагниченность насыщения, a – параметр формы. Также в модели используются коэффициент сцепления k, параметр c, характеризующий выпуклость границ домена, и коэффициент междоменной связи α. Перечисленные параметры можно найти по измеренной петле гистерезиса, используя метод Нелдера – Мида [5] безусловной оптимизации.
В случае, когда междоменным взаимодействием можно пренебречь и положить α = 0, для описания кривой M(H), выходящей из известной начальной точки (H0, M0), можно использовать явные выражения, устойчивые к ошибкам численной реализации [3].
При 
рабочая точка будет двигаться вправо, т.е. в сторону возрастания H. Возможны два случая. Если 
, то рабочая точка всегда будет находиться под Man(H), а ее движение будет описываться кривой
(2)
где 
. Если 
, то траектория M(H) и кривая Man(H) пересекутся в точке (Hcr, Mcr), определяемой из системы уравнений
, (3)
. (4)
При 
траектория будет описываться кривой
, (5)
а при H ≥ Hcr – выражением (4), где вместо (H0, M0) необходимо подставить (Hcr, Mcr).
При 
рабочая точка будет двигаться влево, т.е. в сторону убывания H. Возможны два случая. Если 
, то рабочая точка всегда будет находиться над Man(H), а ее движение будет описываться кривой
, (6)
где 
. Если 
, то при 
траектория определяется согласно (5), а при H ≤ Hcr – выражением (6), где вместо (H0, M0) необходимо подставить (Hcr, Mcr), определяемую из уравнений (3),(4).
Для общего случая при α ≠ 0, для определения M(H) можно использовать итерационный процесс:
, 
, (7)
, (8)
где функции M(Hj) определяются согласно (2)–(6). Результатом является функция 
, которая может быть подставлена в (1) для определения индукции.
Вычислительные правила (1)–(8) позволяют для любой начальной точки (H0, B0) выделить из петли гистерезиса однозначную зависимость B(H), представляющую собой непрерывную, возрастающую, гладкую всюду, кроме точки (H0, B0), функцию (рис. 1, а). Такая зависимость соответствует характеристике ферромагнитного материала при любом монотонном изменении значения напряженности H магнитного поля относительно начального заданного значения H0 и может быть использована для моделирования магнитного поля постоянных токов и постоянных магнитов. Каждая выделенная линия B(H) имеет одну точку пересечения с осью H – остаточную индукцию Br, и с осью B – коэрцитивную силу Hc.
а) б)
Рис. 1. Модель Джилса – Атертона магнитного гистерезиса: 1 – начальная кривая намагничивания, 2 – петля гистерезиса, 3 – выделенная однозначная зависимость, 4 – модель Джилса – Атертона: а) однозначная характеристика B(H) для начальной точки (H0, B0), б) петля гистерезиса и начальная кривая для модели Джилса – Атертона
В доступных программах анализа физических полей Femm и Elcut для моделирования кривой размагничивания B(H), не проходящей через начало координат, используется подход, предполагающий смещение этой кривой вправо на величину коэрцитивной силы [6, 7]. Магнитная проницаемость μ в этом случае определяется по формуле
. (9)
Смещение кривой размагничивания в [6, 7] достигается размещением по границам магнита эффективных поверхностных токов с поверхностной плотностью ±Hc. Установлено, что аналогичный эффект достигается заданием остаточной намагниченности, значение которой 
. Такой подход использован в работе для расчета поля для ферромагнетика с выделенной однозначной характеристикой B(H).
Для анализа распределения плоскопараллельного магнитного поля была использована модификация метода граничных интегральных уравнений, предполагающая одновременное использование скалярных потенциала простого слоя (ППС) и потенциала двойного слоя (ПДС) [4]. Согласно выбранному методу магнитная индукция ищется в виде суперпозиции:
где 
– безвихревое поле реакции ферромагнетиков, а 
– поле постоянных магнитов, создаваемое в однородной среде с магнитной проницаемостью μ0.
Постоянный магнит представляется в виде совокупности элементарных однородно намагниченных постоянных магнитов, границы каждого из которых покрываются слоем фиктивных магнитных зарядов простого слоя с известной плотностью σ = μ0(Mr)n, где (Mr)n – нормальная составляющая остаточной намагниченности. Индукция, наводимая магнитом с индексом PM в точке Q, определяется по формуле
,
где 
вычисляется интегрированием по границам 
магнита:
,
,
– ППС, 
– расстояние от точки интегрирования P до точки наблюдения Q, 
вне магнита PM, а внутри 
,
где 
, 
– проекции вектора 
на оси координат x и y, соответственно.
Каждый ферромагнетик разбивается на элементарные области однородности магнитной проницаемости μ, а поле реакции 
моделируется введением скалярного ПДС, двойной слой зарядов которого размещается на границах 
областей однородности:
,
,
,
где 
– нормаль в точке P, τ – неизвестная плотность ПДС, которая определяется из интегрального уравнения второго рода
, (10)
где точка Q пробегает границы всех элементарных областей разбиения ферромагнетиков 
, 
, 
и 
– магнитные проницаемости со стороны положительного и отрицательного направлений нормали 
в точке Q соответственно.
На основе приведенных формул в работе получен алгоритм моделирования магнитного поля, предполагающий выполнение следующих этапов:
1) замена рассматриваемой области расчета кусочно-однородной средой, в том числе разбиение ферромагнетиков на области однородности SE и задание в каждой из них начальных значений магнитной проницаемости;
2) размещение первичных источников магнитного поля (плотностей ППС на границах элементарных постоянных магнитов);
3) определение плотностей скалярного ПДС на границах элементарных областей SE путем решения интегрального уравнения (10);
4) расчет усредненного значения магнитной индукции в каждой области SE:
,
где 
– внешняя нормаль, а 
и 
– магнитные проницаемости с внутренней и внешней сторон границы LE в точке 
;
5) корректировка магнитной проницаемости в каждой элементарной области SE с применением параметра релаксации ν:
,
где μk и μk+1 – старое и новое значения магнитной проницаемости соответственно, μ(Hcp) рассчитывается согласно (9), а значение Hcp определяется из уравнения 
с помощью метода половинного деления отрезка (дихотомии). Выполнение алгоритма завершается, когда для каждой элементарной области однородности SE выполняется 
, где ε – заранее заданное максимально допустимое значение погрешности определения магнитной проницаемости. Иначе переход на шаг в).
Результаты исследования и их обсуждение
Рассмотрена задача моделирования плоскопараллельного магнитного поля уединенного ферромагнитного образца C45 прямоугольной формы с остаточной намагниченностью, направленной вдоль оси y (рис. 2, а). Результаты сравнивались с решением, полученным в Femm 4.2. Для материала C45 на тороидальном образце с использованием прибора MagHyst была получена петля гистерезиса с начальной кривой намагничивания (рис. 1, б), определены значения остаточной индукции и коэрцитивной силы: Br = 0,949 Тл, Hc = –682,9 А/м. В Femm была экспортирована кривая размагничивания со смещением вправо на величину 
.
По имеющейся петле гистерезиса были получены параметры модели Джилса – Атертона: Ms = 1460000 А/м, a = 1230 А/м, k = 870 А/м, c = 0,17, α = 0,0023, количество повторений n в (8) принято равным 6 (рис. 1, б). Построенная модель характеризовалась следующими значениями остаточной индукции и коэрцитивной силы: Br = 0,939 Тл, Hc = –687,2 А/м. В качестве однозначной характеристики B(H) была выбрана кривая размагничивания модели Джилса – Атертона. В модели была принята кусочно-постоянная аппроксимация плотности ПДС, количество граничных элементов оказалось равным 336.
а) б)
Рис. 2. Расчет магнитного поля уединенного образца C45: 1 – результат, полученный предложенным методом, 2 – график, импортированный из FEMM 4.2: а) область расчета: ферромагнетик, разбитый на 21 подобласть, б) графики y – составляющей магнитной индукции вдоль оси y
Результаты расчетов представлены на рис. 2, б. Хорошее совпадение графиков демонстрирует корректность предложенного метода и объясняется хорошим приближением значения коэрцитивной силы при заданных параметрах модели Джилса – Атертона.
Заключение
В работе предложена модификация метода граничных интегральных уравнений, предназначенная для моделирования стационарного плоскопараллельного магнитного поля в нелинейных средах с гистерезисом и основанная на выделении из области, ограниченной петлей гистерезиса, однозначной характеристики, задаваемой с использованием явных выражений модели Джилса – Атертона. Комбинированное использование скалярных ППС и ПДС совмещает преимущества обоих способов размещения зарядов, а именно:
- устойчивость решения в области, занимаемой ферромагнетиком, обеспечиваемую использованием ПДС;
- непрерывность и конечность значений плотности фиктивных магнитных зарядов ПДС, что упрощает ее аппроксимацию;
- простоту моделирования поля постоянных магнитов, обеспечиваемую применением ППС.
Несмотря на то, что в работе рассмотрено только плоскопараллельное поле, предложенный метод может быть легко перенесен на случай трехмерной задачи. При этом использование скалярных потенциалов не вызовет избыточного увеличения размерности, в отличие от векторных. Метод может быть применен к моделированию медленно меняющегося магнитного поля. Для этого достаточно выполнить фиксацию состояния ферромагнитного материала в каждой элементарной области и, используя эти точки в качестве начальных, построить новые зависимости B(H).
Результаты работы получены при поддержке ФГБОУ «Фонд содействия развитию малых форм предприятий в научно-технической сфере» в рамках НИР «Разработка программных модулей моделирования размагничивания постоянных магнитов для проектирования высокоэффективных электроприводов мотор-колес» (11690ГУ/2017).