В ней векторы и можно выразить через орты осей системы:
, (18)
. (19)
Тогда в стационарном состоянии вектор суммарного поля (6) можно записать в виде:
. (20)
Его абсолютное значение
. (21)
Тогда выражение (17) можно записать в виде:
. (22)
Подставляя (22) в скалярное произведение (16) и используя (18), получим:
. (23)
Если для стационарного состояния поле анизотропии из (8) записать в виде
, (24)
уравнение (23) может быть преобразовано в такое:
, (25)
где .
Для удобства введем обозначения , , и запишем функцию из правой части (25) в виде:
. (26)
Чтобы получить представление о характере функции (26), построим ее графики при различных значениях , , (рис. 3, 4).
Рис. 3. Графики функции при , , (слева) и , , (справа) (штрихом показан график функции )
Рис. 4. Графики функции при , , (слева) и , , (справа) (штрихом показан график функции )
Как видно из рис. 3, 4, в интервале допустимых значений графики функций и пересекаются от одного до трех раз. Это означает, что при различных параметрах , , существует от одного до трех решений уравнения (25), соответствующих возможным положениям вектора статической намагниченности. Однако решение, для которого производная , соответствует неустойчивому положению равновесия, поэтому при расчете петель гистерезиса оно должно исключаться из рассмотрения.
Таким образом, в самом общем случае уравнение (25) имеет два решения и , соответствующие возможным устойчивым положениям вектора статической намагниченности (далее для краткости будем называть их устойчивыми решениями). При любых допустимых значениях параметров , , эти решения имеют разные знаки, что позволяет довольно простым способом учесть гистерезис системы при расчете полевых характеристик. Для этого достаточно при каждом новом значении напряженности внешнего поля из двух решений оставлять то, которое совпадает по знаку с решением, полученным для предыдущего значения . Такой способ позволяет автоматически учесть гистерезис при построении кривых перемагничивания в «прямой» форме, когда независимой переменной является напряженность внешнего поля .
Однако в силу громоздкости выражения (26) аналитическое решение уравнения (25) не представляется возможным. Поэтому при расчете кривых перемагничивания следует воспользоваться численными методами его решения.
Численное решение уравнения в задаче об определении направления вектора статической намагниченности
Прежде чем выбрать какой-либо метод численного решения, отметим одну закономерность, которая остается справедливой при любых допустимых значениях параметров , , и позволяет сформулировать необходимые условия сходимости последовательности приближений () к устойчивому решению . На рис. 3, 4 можно видеть, что если для некоторого значение функции превышает значение аргумента
, (27)
то устойчивое решение находится справа от , поэтому для приближения к нему следует увеличить аргумент на следующей итерации:
. (28)
Если же значение функции меньше значения аргумента
, (29)
то устойчивое решение находится слева от , поэтому для приближения к нему следует уменьшать аргумент на следующей итерации:
. (30)
Указанные закономерности являются необходимым условием сходимости последовательности к устойчивому решению , поэтому численные методы, используемые для решения уравнения (25), должны формировать последовательности, удовлетворяющие условиям (27) - (30).
К настоящему времени известно большое количество методов решения нелинейных уравнений вида (25) [6]. Среди них можно выделить довольно универсальный метод Ньютона, часто называемый методом касательных. Достоинствами метода Ньютона являются сравнительная простота его реализации (метод является одноточечным, т.е. для вычисления приближения к решению достаточно знать приближение и нет необходимости в информации о предыдущих , и др. приближениях) и более высокой скоростью сходимости по сравнению с другими одноточечными методами (например, методом деления пополам).
Согласно методу Ньютона, за приближение к решению принимается абсцисса точки пересечения функции и касательной к функции , проведенной в точке (рис. 5).
Рис. 5. Геометрическая интерпретация метода Ньютона
Для построения алгоритма вычисления приближения по методу Ньютона запишем уравнение касательной к функции , проведенной в точке с координатами :
. (31)
Приравнивая функции и , получаем выражение для определения приближения:
. (32)
гексаферритовый гиромагнитный вектор
Итак, для вычисления приближения по методу Ньютона нужно знать значение функции и ее производной в точке . Это требование принято считать недостатком метода, поскольку вычисление производной в каждой точке последовательности может привести к существенному увеличению объема вычислительных операций.
Для проверки этого опасения найдем производную функции (26). Для этого введем обозначения
, (33)
. (34)
При таких обозначениях функцию (26) можно представить в виде частного , поэтому ее производная вычисляется по правилу дифференцирования дробей:
. (35)
Производная числителя (33):
. (36)
Производная знаменателя (34):
. (37)
Тогда выражение (35) можно преобразовать следующим образом:
. (38)
Таким образом, функция и ее производная связаны, поэтому, если при расчете производной в точке использовать результаты расчета функции в той же точке, можно несколько сэкономить на вычислениях. Таким образом, есть основания полагать, что при правильном построении алгоритма численного решения использование метода Ньютона не должно приводить к значительному увеличению объема вычислительных операций по сравнению с методами, для которых вычисление производной не требуется.
Более существенным недостатком метода Ньютона является возможность выхода за пределы интервала , что может привести к непредсказуемым результатам развития итерационного процесса. Рассмотрим подробнее ситуации, когда это может произойти.
1. Если для текущего приближения значение функции меньше значения аргумента, то по условию (30) -й элемент формируемой последовательности должен быть меньше -го. В этом случае следует опасаться выхода за нижнюю границу допустимого интервала изменения величины :
. (39)
Используя выражение (32), получаем:
. (40)
Считая, что (в противном случае последовательность, формируемая методом Ньютона, будет либо расходящейся, либо сходящейся к неустойчивому решению), перепишем (40):
. (41)
Тогда
, (42)
что равносильно условию для производной функции (26)
. (43)
2. Если для текущего приближения значение функции больше значения аргумента, то по условию (28), -й элемент последовательности должен быть больше -го. В этом случае следует опасаться выхода за верхнюю границу допустимого интервала изменения величины :
. (44)
Используя выражение (32), запишем
. (45)
Считая, что , перепишем это выражение:
. (46)
Тогда
, (47)
что равносильно условию для производной функции (26)
. (48)
Подводя итоги, отметим, что при выполнении условий (29), (43) или (27), (48) получаемое по методу Ньютона приближение выходит за пределы допустимого интервала . В этих случаях метод Ньютона становится неприменим, поэтому для вычисления -го приближения следует воспользоваться каким-то другим методом.
В качестве такого метода предлагается использовать метод простой итерации. Одним из достоинств метода является его простота реализации: в качестве -го приближения к решению используется значение функции (26) от -го приближения (рис. 6):
. (49)
Рис. 6. Геометрическая интерпретация метода простой итерации
Следует отметить, что необходимым для сходимости метода является условие , которое не выполняется в окрестности неустойчивого решения уравнения (25). Таким образом, метод простой итерации автоматически расходится в этой окрестности, поэтому его можно использовать для подстраховки и попадания в область сходимости метода Ньютона. Для иллюстрации эффективности такого способа решения на рис. 7 показан пример использования комбинации методов Ньютона и метода простой итерации.
Рис. 7. Графическая иллюстрация решения уравнения (25) при параметрах , , , комбинацией методов Ньютона и простой итерации
Условием прекращения вычисления элементов последовательности будем считать примерное равенство значений аргумента и функции на -й итерации:
, . (50)
По завершении итерационного процесса решением уравнения (25) будем считать величину :
. (51)
Зная , можно рассчитать по формуле (9).
Таким образом, можно сказать, что для общего случая произвольной напряженности внешнего магнитного поля и произвольной ориентации частицы в нем решена задача о направлении вектора статической намагниченности частицы с известными магнитными свойствами. Изложенное решение использовано для разработки компьютерной процедуры «Статика», алгоритм работы которой представлен на рис. 8; на рис. 9 представлен алгоритм работы используемой в ней процедуры «Численное решение».
Рис. 8. Блок-схема алгоритма работы процедуры «Статика»
Рис. 9. Блок-схема алгоритма процедуры «Численное решение»
Входными параметрами процедуры «Статика» являются константы поля анизотропии , , угол ориентации частицы во внешнем поле , текущее значение напряженности внешнего поля (подразумевая использование процедуры в теле цикла, организованного для расчета кривых перемагничивания, напряженность внешнего поля обозначена квадратными скобками как элемент массива) и начальное приближение к решению . Установочными параметрами процедуры являются точность решения статической задачи и малые константы , (), используемые для обозначения границ применимости частных случаев, для которых решение задачи возможно без привлечения процедуры «Численное решение». Выходными параметрами процедуры являются значения и .
Расчет кривых перемагничивания и петель магнитного гистерезиса
Итерационный процесс, используемый в процедуре «Численное решение», предполагает наличие некоторого начального приближения , с которого начинается формирование последовательности.
Как указывалось выше, в общем случае существует два устойчивых решения уравнения (25), одно из которых характеризуется острым углом (), второе - тупым (). При использовании описанных выше численных методов решения уравнения (25) для поиска положительного корня удобно использовать в качестве начального приближения величину , а для поиска отрицательного корня - величину . Даже если уравнение (25) имеет только одно решение, которое не совпадает по знаку с выбранным начальным приближением , комбинация методов Ньютона и последовательных приближений позволяет сформировать последовательность, сходящуюся к этому решению (для примера можно снова обратиться к рис. 7).
Для учета магнитной предыстории системы при расчете кривых перемагничивания и петель гистерезиса будем выбор начального приближения для текущего значения напряженности внешнего поля проводить на основе знака решения, полученного для предыдущего значения напряженности поля :
. (52)
Если при намагниченность частицы соответствует верхнему участку петель гистерезиса, то в качестве приближения будем использовать
, (53)
если нижнему - приближение
. (54)
На рис. 10 приведены схемы алгоритмов расчета верхнего и нижнего участков петель гистерезиса.
Рис. 10. Блок-схемы алгоритмов расчета верхнего (слева) и нижнего (справа) участков петель гистерезиса
Входными параметрами алгоритмов, помимо упоминавшихся , и , являются границы диапазона изменения напряженности поля , . Установочным параметром является количество точек на графике . Выходными параметрами являются массивы значений напряженности внешнего поля и косинусов соответствующих им углов .
В заключение приведем пример расчета кривых перемагничивания и петель гистерезиса при различной ориентации частицы во внешнем поле. Для наглядности и возможности сопоставления результатов с расчетами, проведенными параметрическим способом в работе [3], положим вторую константу поля анизотропии равной нулю (), а по осям координат будем откладывать нормированные величины напряженности поля и проекции статической намагниченности .
На рис. 11 показан независимый расчет верхнего и нижнего участков петель гистерезиса частицы, гексагональная ось которой ориентирована под углом к внешнему полю.
Рис. 11. Графики кривых перемагничивания частицы при ,
В первом случае нормированная напряженность поля уменьшается от значения до , во втором случае - увеличивается от до . В обоих случаях в графике функции наблюдается разрыв, вызванный скачкообразным переходом вектора статической намагниченности из одного устойчивого положения в другое при превышении напряженностью поля критического значения . Значения различаются для разных направлений изменения внешнего поля (от максимального значения до минимального или наоборот), т.е. зависят от предыстории системы. Если же изменить внешнее поле сначала от до , а потом в обратную сторону, можно получить график петель гистерезиса частицы (рис. 12):
Рис. 12. График петель гистерезиса частицы при ,
Как уже отмечалось, разработанная методика позволяет провести расчет кривых перемагничивания при произвольной ориентации частицы во внешнем поле. Для примера на рис. 13 показано семейство петель гистерезиса, построенных для разных значений угла .