Содержание к диссертации
Введение
Глава 1. Кинетические модели магнитных материалов 16
1.1. Физико–математическая модель 17
1.2. Вывод уравнения Фоккера–Планка 18
1.3. Многочастичное уравнение Фоккера–Планка для ансамбля из магнитных частиц 22
1.4. Уравнение Брауна для суперпарамагнетиков 28
1.5. Уравнение Фоккера–Планка для ферромагнетика 29
1.6. Безразмерные системы единиц 33
Глава 2. Локальный метод Галеркина на сфере 35
2.1. Сетка и базис 37
2.2. Локальный метод Галеркина 43
2.3. Численные потоки 45
2.4. Ограничитель WENO 47
2.5. Тестирование кода 48
Глава 3. Аналитические оценки и результаты моделирования 51
3.1. Стохастический резонанс в суперпарамагнетиках 51
3.2 . Фазовый переход в ферромагнетике 76
3.3. Перемагничивание ферромагнетика внешним полем 82
Заключение 84
Сокращения и обозначения 86
Список литературы
- Многочастичное уравнение Фоккера–Планка для ансамбля из магнитных частиц
- Безразмерные системы единиц
- Численные потоки
- . Фазовый переход в ферромагнетике
Многочастичное уравнение Фоккера–Планка для ансамбля из магнитных частиц
Как и говорилось выше, в качестве исходной точки построения математической модели мы возьмем систему уравнений Ландау-Лифшица (УЛЛ), записанную для отдельных ферромагнитных частиц где 7 — гиромагнитное отношение, a — параметр затухания, М — величина магнитного момента (во всех дальнейших выкладках предполагается равной 1), Н ef — эффективное магнитное поле, — трёхмерный белый шум, приложенный к частице г, VK6 11 — энергия обменного взаимодействия, VKams — энергия анизотропии, поле размагничивания микрочастицы (диполь-дипольное взаимодействие) и однородное внешнее магнитное поле эффективно учитываются в анизотропии, индексы і и j нумеруют атомы, суммирование по ним указывается явно, /І и v — координатные индексы (здесь и далее обозначаются греческими буквами), по повторяющимся координатным индексам предполагается суммирование, Jij — обменный интеграл между частицами і и j. Т — температура системы измеряемая в единицах энергии (&&Т, постоянную Больцмана мы для краткости писать не будем). Во многих работах коэффициент затухания счита 18 ют безразмерным. Таким образом, коэффициент перед соответствующим членом оказывается orfjM. Это необходимо учитывать при подборе характерных параметров.
Также в формулах используется обозначение () для усреднения по реализациям шума. В дальнейшем мы будем пользоваться таким обозначением для усреднения, его смысл в конкретных ситуациях будет уточняться. Для суперпарамагнетика обменное взаимодействие между ферромагнитными частицами отсутствует, так что энергия обменного и дипольного взаимодействий равна 0. Анизотропию в этих уравнениях связывают с формой частиц и кристаллической решеткой. Наиболее часто используются одноосная [51] (линейная анизотропия типа «легкая ось») и кубическая [50] (возникающая, в частности, из-за кристаллической решетки) формы анизотропии. Мы будем использовать более общую одноосную анизотропию с линейным и кубическим изменением поля по направлениям. Кубическая анизотропия может быть представлена с использованием трех легких осей такой анизотропии. В дальнейшем мы будем говорить исключительно об анизотропии типа «легкая ось», а линейная и кубическая будут означать степень, с которой входят в поле компоненты направления.
Система УЛЛ (1.1) — это система стохастических дифференциальных уравнений с мультипликативным шумом. Эта система описывает некоторый марковский процесс (поскольку шум дельта-кореллирован), для которого можно записать уравнение Фоккера-Планка (УФП). В следующем параграфе представлен один из методов вывода УФП для некоторого абстрактного макровского процесса, затем вычислены коэффициенты УФП непосредственно для системы УЛЛ (1.1).
Рассмотрим систему частиц, движущихся под действием случайной силы h. Данная система характеризуется векторной координатой х в фазовом пространстве G. Приведем вывод уравнения УФП. Предположение о том, что h является белым шумом, позволяет нам считать процесс х() марковским. В этом параграфе мы также будем использовать угловые скобки () для интегрирования по пространству элементарных исходов h.
Кроме того, нам понадобятся обозначение /(х,) для TV-частичной функции распределения и стандартное обозначение для условной функции распределения /(х,xi,i). Здесь под х понимается вектор координат системы в фазовом пространстве G. Заметим, что вектор х уже является независимой координатой, значит он не зависит от времени. Под х() понимается реализация случайного процесса х. Запишем определение функции распределения /
Выпишем уравнение Колмогорова-Чепмена, аналог формулы полной вероятности для условной функции распределения марковского процесса /(х, + At) = f(y,t)f(x,t + At\y,t) dy. (1.3) G Будем считать, что / экспоненциально убывает при удалении от невозмущенной траектории, идущей из точки Xo,to. Рассмотрим, как изменяется момент некоторого многочлена Р(х) со временем
Воспользуемся определением предела и поменяем порядок интегрирования и взятия предела Поскольку функция /(x, t + Aty,t) при малых At экспоненциально убывает при удалении от у, а функция Р растет полиномиально, мы можем считать, что интегрирование проводится в малой окрестности у, размер которой зависит от At. В этой окрестности мы можем разложить функцию Р в ряд вокруг точки у, Здесь и далее символом “” будет обозначаться тензорное произведение. Подставим выражение (1.6) в уравнение (1.4). Разобьем правую часть на сумму интегралов:
Безразмерные системы единиц
Для решения кинетических уравнений часто применяют метод стохастического аналога [82]. Однако формулировка этого метода для решения уравнения Брауна (1.24) будет мало отличаться от формулировки метода решения уравнения Ландау–Лифшица (1.1). Для задач с взаимодействием между частицами эти методы уже существенно различаются, но потребует набора большого количества реализаций для получения статистически достоверных результатов [61, 68]. Другим вариантом решения УФП могут служить сеточные методы [58].
Дифференциальные уравнения на сфере возникают во многих физических задачах от предсказания погоды до расчета диаграмм направленности антенн. Во многих постановках, особенно в эволюционных задачах, для получения адекватного решения требуются методы, в которых гарантируется выполнение законов сохранения, присущих системе. Кроме того, для получения приемлемой точности решений без существенного повышения вычислительной сложности, нужно использовать методы высокого порядка точности.
Для выполнения этих требований удобно использовать наработки вычислительной гидродинамики. Наиболее перспективными здесь являются разрывные методы Галеркина (DG) [83]. Эти методы легко обобщаются на любой необходимый порядок точности и обладают консервативными свойствами. За последние годы были построены адаптивные варианты DG методов [84], монотонизирующие ограничители (лимитеры), как общие так и специальные [85, 86], рассчитанные на конкретную задачу. Разрывные методы Галеркина в очень малой степени ограничивают выбор сеток, но их точность безусловно зависит от качества сетки. При решении многих задач ошибки, обусловленные сеточной анизотропией, искажают результаты не только количественно, но и качественно. Традиционная сетка в сферических координатах является простейшей для построения, но обладает особенностями на полюсах, которые создают существенные сложности при её применении. Другая часто используемая сетка называется cubed sphere [87] и строится как проекция разрезанного куба на сферу. Эта сетка очевидно анизотропна и имеет большие искажения в окрестностях углов куба, хотя уже является приемлемой. Есть множество других вариантов построения сеток на сфере, например, в статье Giraldo [88] используется сетка, построенная на основе икосаэдра. Мы будем строить рекурсивную сетку на основе додекаэдра [89]. После первичного разбиения его граней на треугольники получается базовая сетка, состоящая из почти правильных одинаковых треугольников. При ее дальнейшем разбиении и проекции на сферу искажения оказываются существенно меньше, чем при разбиении граней икосаэдра. Для апробации метода мы промоделируем суперпарамагнетик — ансамбль монодоменных ферромагнитных микрочастиц. Суперпарамагнетики описываются уравнением типа Фоккера–Планка на сфере, полученным Брауном в 1963 г. [18]. Это одно из простейших уравнений физической кинетики, записанное на сфере. Для некоторых частных случаев возможно получить точное аналитическое решение уравнения Брауна, что позволяет протестировать метод. Кроме того, в суперапармагнетиках наблюдается стохастический резонанс. Использование разрывного метода Галеркина допускает сетки достаточно общего вида. Например, в работе Giraldo [88] использовалась сетка, состоящая из криволинейных четырехугольников. Мы будем использовать разбиение T/j = {КІ,І = 1}N] сферической области S2 на N сферических треугольников КІ. При моделировании прецессии функции распределения учет кривизны ячеек является существенным, так как на сетке из плоских треугольников мы неминуемо получим систематическую фазовую ошибку.
Построение сетки на сфере является отдельной задачей, требующей отдельного рассмотрения. Мы решили использовать рекурсивный алгоритм построения сетки на сфере, геометрическая часть которого изложенна в работе [89]. Одним из его преимуществ является легкость локального разбиения ячеек и, соответственно, построение адаптивных методов 1. Другим достоинством подобного метода построения является естественная локальность данных. Стоит упомянуть, что полученная сетка будет обладать большим количеством симметрий. Это свойство используется, в частности, для уменьшения количества хранимых и вычисляемых геометрических данных. Кроме того, получаемая сетка обладает достаточно хорошими свойствами: площади её ячеек отличаются менее чем на 10% от минимальной площади и минимальный угол в сетке имеет величину не менее 37г/10 вне зависимости от мелкости разбиения.
Численные потоки
В случае нелинейного уравнения для намагниченности M этот подход предполагает малость амплитуды А периодического вынуждающего поля Н , поскольку можно считать уравнение (3.5) первым членом разложения по величине амплитуды.
Намагниченность M в терминах функции распределения будет определяться как среднее направление магнитного момента M = (m) (t) = mf(m,t)dm, (3.8) S2 поскольку распределения для всех частиц независимы. Для простоты мы будем предполагать, что возмущающее поле имеет только одну ненулевую компоненту. Это позволяет нам сократить количество величин, определяющих установившееся движение, до первого столбца матрицы , который мы будем далее обозначать %. Фурье-образ периодического поля Н из постановки (1.25) вычисляется по формуле Н = А 5(Q — ш) — 5 (Q + uS) ш 2г Для такого поля мы можем переписать соотношение линейной теории отклика (3.7) во временной области %(Г2)е ш — %(—Q) е ІШ Л (m) it) = А . 3.9 0.5 0.4 0.3 0.2 0.1 0 (а) 0.2 —0.1 0 0.1 0.2 0.3 тх 50 30 10 (б) 0.9 0.8 0.7 0.6 0.5 К- , 0.2 —0.1 0 0.1 0.2 0.3 тх 50 40 30 20 10 Рис. 3.3. Релаксация эволюции намагниченности к периодическому движению в случае, когда поле Но слабее максимального поля анизотропии, при температуре Т = 0.05 и частоте вынуждающей силы = 1.2. Ось анизотропии направлена под 80 к периодическому полю Н и 30 к внешнему полю.
Заметим, что с формальной точки зрения ось, по которой амплитуда колебаний максимальна, не совпадает с направлением периодического поля. Кроме того, восприимчивость на ортогональном направлении не обязана быть существенно ниже, чем на основном. Например, на рис. 3.3 (а), полученном в результате численного решения уравнения Брауна (1.24) разрывным методом Галерки-на, видно, что амплитуды колебаний вдоль поля и поперек него практически совпадают, по оставшейся оси колебания примерно в 3 раза слабее (рис. 3.3 (б)). В большинстве случаев восприимчивость вдоль направления поля выше в несколько раз, тем не менее, её значение того же порядка малости по амплитуде периодического поля. Мы будем вычислять только одну компоненту восприимчивости X, ту, которая определяет поведение проекции намагниченности на вынуждающее поле, и будем обозначать её Хх. Тем не менее, представленная в работе техника позволяет проводить вычисления и строить аналитические оценки для всех компонент комплексной восприимчивости.
Для первой компоненты векторного уравнения (3.9), пользуясь свойством четности, х {— ) = X ( )? мы получаем уравнение, которое может быть использовано для определения Хх (A (m) (t) А) = Л2 sin (Ш + ф) %Ж(Г2), (3.10) ф{1) = Arg (;\;Ж(Г2)). Таким образом, вычислив временную зависимость проекции намагниченности на направление поля (m) (t) А/Л, несложно получить численные оценки для восприимчивости. Поскольку нас интересуют только установившиеся колебания, то накопление данных следует начинать по прошествии достаточного для релаксации времени. Характерным масштабом времени релаксации системы является величина 1/а, поэтому время, необходимое для релаксации, должно существенно превышать эту величину тгеіах 1/а. При расчетах можно считать, что для установления движения системы достаточно времени тгеіах = Ъ/а.
В численном решении уравнения Брауна (1.24) на больших временах колебания намагниченности асимптотически стремится к синусоидальным колебаниям, задаваемым формулами (3.10). По этой причине необходимо дождаться, пока система достаточно релаксирует к равновесию. Кроме того, период колебаний может отклоняться от стационарного, что тоже вызывает ошибку в измерениях. По этой причине необходимо проводить интегрирование на достаточно большом временном интервале, чтобы величина внесенной ошибки была относительно мала, т.е. число периодов N, на котором проводится интегрирование в численном моделировании, должно быть достаточно велико. В наших расчетах интегрирование проводилось на N = 40 периодах колебаний. Из выведенного соотношения (3.10) легко получить формулы для численного расчета
. Фазовый переход в ферромагнетике
Задача о перемагничивании ферромагнетика встречается при моделировании почти всех вариантов устройств хранения памяти. Обычно исследователей интересует время, необходимое для перемагничивания слоя заданной геометрии, которое можно интерпретировать как время записи байта информации. Другой похожей задачей является вычисление вероятности ошибки записи, то есть вероятность, того что при заданном изменении внешних параметров задачи со временем направление намагниченности не изменится.
Часто перемагничивание хранящего информацию слоя идет за счет взаимодействия со спин-поляризованным током (так называемый spinorque). Мы рассматриваем существенно более простую постановку, в которой перемагничи-вание происходит просто за счет действия однородного внешнего поля. В этой задаче мы считаем размер образца порядка размера одной ячейки микромагнитного подхода, то есть порядка 20 параметров решетки. Хотя образец с периодическими граничными условиями формально является бесконечным, конечный размер области моделирования навязывает системе сверхдальний порядок.
Перед перемагничиванием система должна находиться в равновесном состоянии, причем слой должен быть намагничен в определенном направлении. Для фиксации направления намагниченности можно использовать анизотропию (в данном случае связанную с формой). Легкая ось линейной анизотропии a1 направлена вертикально. В начальный момент времени расчета внешнее поле выключено и все магнитные моменты направленны вверх. Мы дожидаемся установления равновесия в системе, где температура заметно ниже критической. Вдали от точки фазового перехода скорость установления равновесия по 83 рядка reiax 1/, так что за время Б/ система приходит в состояние близкое к стационарному. После этого включается слабое относительно анизотропии магнитное поле в противоположном намагниченности направлении и вычисляется зависимость намагниченности от времени.
Размер сетки составляет 240 ячеек, что совпадает с размером сетки использовавшейся для расчетов фазового перехода, шаг интегрирования был выбран исходя из условия Куранта.
Для задачи перемагничивания результаты моделирования в различных приближениях отличаются качественно (рис. 3.9).В задаче использовалась квадратичная анизотропия ms = (т , п) /2. В начальный момент времени все магнитные моменты параллельны п, величина анизотропии = 1, внешнее поле, включающееся через Ъ/ вдвое слабее анизотропии 1/2, температура существенно меньше температуры Кюри = 1, обменный интеграл для ближайших соседей ij = 1. LLE отвечает решению УЛЛ (1.1), CWT и IBE — решению уравнения Фоккера-Планка (1.28) для приближения мультипликативности (1.31) и феноменологической модели (1.33) соответственно. Для IBE i = —1/4, 2 = 4. В моделировании магнитной системы методом атом-в-атом, область разбивается на микродомены, которые разворачиваются через разные направления. По этой причине намагниченность (рис. 3.9 (а)) разворачивается через практически нулевое значение, в то время как намагниченность для мультипликативного подхода ведет себя качественно так же, как и магнитный момент отдельной частицы. На рисунке 3.9 (б) представлена зависимость модуля намагниченности от времени. Учет корреляций дает количественно неверный результат времени перемагничивания. Возможно этот результат можно улучшить, если учитывать корреляции более точно. Таким образом, в задачах с неравновесной динамикой приближение среднего поля может давать качественно неверные результаты.
1. Численный метод второго порядка точности из семейства локальных методов Галеркина адаптирован для решения уравнений Фоккера–Планка, описывающих распределение намагниченности по направлениям. Метод обладает существенно меньшей вычислительной сложностью, чем прямое решение уравнений Ландау–Лифшица со случайным источником для расчетов сопоставимой точности, и универсален в отличие от метода цепных матричных дробей.
2. Численный метод реализован в виде высокопроизводительного кода, решающего УФП на рекурсивно треугольной сетке, обладающей большим количеством симметрий и локально рекурсивным хранением данных, что использовалось для повышения эффективности расчетов. Второй порядок сходимости метода подтвержден на ряде тестов с известным аналитическим решением.
3. Получены параметры резонанса в широком диапазоне температур, частот и направлений оси кубической анизотропии для случаев, когда максимальное поля анизотропии больше внешнего постоянного поля и наоборот. Проведено сравнение параметров стохастического резонанса с приближенной теорией, приведенной в работе. Сравнение показало, что приближенная теория предсказывает качественно верное поведение, что позволяет использовать её для оценки условий возникновения резонанса.
4. Проведена сравнительная проверка приближений среднего поля и самосогласованных моментов для равновесных и неравновесных задач в ферро 86 магнетиках. В расчетах воспроизводится завышенная на 20% температура фазового перехода для приближения среднего поля. Также показано, что приближение самосогласованных моментов позволяет подобрать параметры, снижающие ошибку до нескольких процентов. Показано, что в неравновесных задачах приближение среднего поля ведет к качественно неверной эволюции намагниченности, в то время как приближение самосогласованных моментов дает качественно верное поведение системы. Этот результат показывает, что для некоторых неравновесных задач в магнитных системах требуется более точный учет корреляций, в частности при использовании уравнений ЛЛБ требуется вводить поправки к приближению среднего поля.