71292

Разработка среды научной визуализации для представления процессов намагничивания

Дипломная

Информатика, кибернетика и программирование

Пакет DOMEN предназначен для расчёта статического распределения и динамического поведения доменных границ в магнитных плёнках при различных физических свойствах материала плёнок и воздействии внешних электрических и магнитных полей различного типа.

Русский

2015-01-19

7.3 MB

3 чел.

Разработка среды научной визуализации для представления процессов намагничивания


ОГЛАВЛЕНИЕ

[1] Введение

[2] 1. Описание области приложения

[3] 2. Статические и динамические свойства ДГ

[3.1] 2.1. Статическая структура доменных границ

[3.2] 2.2. Динамические свойства доменных границ

[4] 3. Математическая постановка задачи

[4.1] 3.1. Применяемый метод микромагнитного моделирования

[4.2] 3.2. Описание программы

[5] 4. Обзор вариантов 2-мерной визуализации

[6] 5. Обзор подходов к 3-мерной визуализации векторных полей

[7] 6. Описание модуля 3D-визуализации

[7.1] 6.1. Инструменты

[7.2] 6.2. Структура программы

[7.3]
6.3. Инициализация

[7.4] 6.4. Чтение данных

[7.5] 6.5. Обработка сигналов с клавиатуры и мыши

[7.6] 6.6. Вспомогательные функции рисования

[7.7] 6.7. Основная функция рисования

[7.8] 6.8. Примеры получаемых 3D-изображений

[8] Заключение

[9] Литература


Введение

Настоящая бакалаврская работа посвящена 3D-визуализации векторного поля магнитного момента. Вектор магнитного момента вычисляется
в результате работы программного пакета
DOMEN, который разработан и уже
в течение 20 лет используется в Институте физики металлов УрО РАН.

Пакет DOMEN предназначен для расчёта статического распределения
и динамического поведения доменных границ в магнитных плёнках при различных физических свойствах материала плёнок и воздействии внешних электрических и магнитных полей различного типа.

Расчёт производится на регулярной сетке в 2-мерной прямоугольной области (постановка задачи описана ниже). В результате расчёта каждому узлу сетки приписывается 3-мерный вектор намагниченности. Вектор в каждой точке имеет единичную длину. При наличии внешнего поля вектор намагниченности изменяется во времени (динамическое поведение). В пакете DOMEN предусмотрена 2D-визуализация распределения намагниченности, при котором отображается проекция распределённого вектора на плоскость области (плоскость дисплея).

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

К данному моменту автор бакалаврской работы создал программный модуль, позволяющий строить 3-мерное векторное поле намагниченности. Каждый вектор строится в виде конуса, цвет которого зависит от значения компонент вектора. Возможно вращение, удаление и приближение распределённого поля. Дополнительная наглядность достигается использованием виртуальных источников освещения и «затуманиванием» удалённых от зрителя векторов. Всё сказанное относится как к рисованию статических распределений, так и к отображению в режиме реального времени (относительно самого расчёта намагниченности) динамически меняющихся распределений.

Перспектива дальнейшего развития 3D-визуализации включает построение 3-мерных поверхностей распределения скалярных величин на расчётной области. Предполагается также перенос созданной методики визуализации на 3-мерные расчётные области.


1. Описание области приложения

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

Пример доменной структуры

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

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

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

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

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

Важным является выяснение влияния внешних полей, толщины и параметров пленок, а также их слоистой структуры на статические и динамические свойства доменных стенок в рамках микромагнитного подхода и двухмерной модели распределения намагниченности при учете всех основных взаимодействий.

Для различных типов устройств (например, головок считывания информации) часто используют многослойные магнитные пленки. Одним из направлений расчётов с помощью пакета DOMEN является выяснение динамического поведения ДГ с двухмерным распределением намагниченности в пленках с различными по магнитным параметрам слоями.

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


2. Статические и динамические свойства ДГ

2.1. Статическая структура доменных границ

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

Первые теории строения доменных границ были изложены в работе Ландау и Лифшица [1] 1935 года. Эта работа явилась фундаментом для всего последующего развития представления о ДГ. Изменение направления магнитных моментов в границе должно происходить на расстояниях, много больших межатомных, так как резкому изменению взаимной ориентации магнитных моментов препятствует обменное взаимодействие. Благодаря этому обстоятельству для изучения и описания структуры ДГ применяется микромагнитное рассмотрение, когда ферромагнетик считается непрерывной средой с намагниченностью, направление которой изменяется в пространстве, а величина постоянна. Модель, описанная в [1], носит название классической или одномерной 180-градусной блоховской стенки. Первым важным модельным представлением, определяющим ее структуру, является предположение об одномерном характере распределения намагниченности M, то есть зависимости M только от одной пространственной координаты. Согласно второму модельному представлению, поворот M в стенке происходит так, что вектор М все время лежит в плоскости доменной границы (то есть в плоскости yz на рис. 1а).

Рис. 1. Схематичная иллюстрация распределения намагниченности в одномерных блоховской (а) и неелевской (б) ДГ. Стрелками обозначены проекции вектора M на координатные плоскости.

Описанная модель была сформулирована без учета наличия у ферромагнитного образца поверхности. Если же учитывать этот фактор, то энергетически более выгодным оказывается уже двухмерное распределение намагниченности. Ниже под двухмерными будем понимать ДГ, у которых направление M изменяется при движении от домена к домену и от нижней поверхности пленки к верхней (то есть вдоль направлений x и y на рис. 1) и не изменяется вдоль направления оси z.

Одномерные блоховская и неелевская ДГ называются симметричными стенками в том смысле, что распределение намагниченности в них обладает зеркальной симметрией относительно плоскости x=0 на рис.1. При этом центр стенки, то есть поверхность, на которой компонента намагниченности Mz=0, совпадает с плоскостью  x=0. Значительным шагом вперед в представлении о внутреннем строении ДГ в магнитных пленках явилось предсказание существования асимметричных стенок. Эти идеи были реализованы в работах Ла Бонте [2] и Хуберта [3] 1969 года. Важно, что в обеих работах были полностью учтены все основные взаимодействия: неоднородное обменное, магнитно-анизотропное и диполь-дипольное.

Впервые полученные в [2, 3] стенки называются двухмерными блоховскими ДГ, а также асимметричными или вихреподобными. Распределение намагниченности в вихреподобной блоховской ДГ иллюстрируется рис.2а. В частях ДГ, непосредственно примыкающих к поверхностям пленки, распределение намагниченности подобно неелевским границам, а в средней части стенки можно выделить участки, характерные для блоховских ДГ. Стыкуясь, эти части образуют вихреподобную структуру, причем центр «вихря» и центральная линия стенки (линия уровня Mz=0) разнесены в пространстве. Таким образом, наблюдается асимметрия распределения намагниченности в стенке относительно ее центральной линии.

Рис. 2. Структура асимметричной блоховской (а) и асимметричной неелевской стенок (б) в пленке пермаллоя толщиной 100 нм. На этом и последующих подобных рисунках показано сечение ДГ плоскостью, нормальной к ОЛН (плоскостью xy на рис. 1), стрелками изображены проекции вектора намагниченности на указанную плоскость


Другим типом двухмерной доменной стенки являются двухвихревые неелевские (так же называемые двухмерными или асимметричными неелевскими). Строение двухвихревой неелевской стенки показано на рис. 2б. Было установлено, что в пленках с параметрами, характерными для пермаллоевых, асимметричные блоховские и неелевские стенки существуют в широкой области толщин (вплоть до микрона) и магнитных параметров.

2.2. Динамические свойства доменных границ

Изучение динамического поведения доменных стенок весьма интересно с научной точки зрения, поскольку даже в сравнительно небольших внешних магнитных полях в их поведении проявляются нелинейные свойства. Кроме того, знание особенностей динамики ДГ позволяет указать пути регулирования  их скоростей движения. Это важно с практической точки зрения, как для увеличения быстродействия устройств, когда требуются большие скорости перемагничивания, так и для уменьшения магнитных потерь, когда требуются малые скорости движения стенок.

Первое исследование динамического поведения переходных слоев между доменами было выполнено Ландау и Лифшицем [1]. Однако в их работе были рассмотрены лишь очень маленькие внешние магнитные поля Н, ориентированные вдоль оси легкого намагничивания (ОЛН), приводящие к малым скоростям движения стенок. Нелинейное поведение ДГ в широкой области значений Н впервые достаточно полно было исследовано Шраером и Уокером [4] 1974 года. Они установили, что динамическое поведение стенок  существенно зависит от величины  магнитного поля Н, направленного вдоль ОЛН. Оказалось, что существует некоторое критическое поле Нс, ниже которого стенка движется стационарно, а выше  ее скорость осциллирует. Скорость стенки взаимосвязана с ее внутренней структурой – в полях, меньших критического, происходит искажение одномерной блоховской стенки (см. рис. 1а), характеризующееся выходом M из плоскости ДГ в сторону ее движения, то есть появлением ненулевой компоненты Mx. При этом ДГ движется,ь но не меняет формы. Такое движение называется квазистационарным. Осцилляции же скорости движения ДГ связаны с возникающей в полях Н>Нс периодической перестройкой ее структуры. В рассматриваемом в [4] случае происходит периодическая перестройка одномерной блоховской стенки в одномерную неелевскую и обратно.

Таким образом, несмотря на то, что в работах [1,4] использовалось предположение о безграничных размерах образцов и одномерном характере распределения намагниченности в стенке, были выяснены необычные особенности динамики доменных стенок. Однако, как следует из сказанного выше, структура доменных границ в образцах ограниченных размеров (в частности, в магнитных пленках) существенно отличается от структуры одномерных блоховских ДГ. Полная динамическая перестройка структуры вихреподобных ДГ на основе численного решения уравнения Ландау-Лифшица в рамках двухмерной модели и строгого микромагнитного подхода впервые была получена в работе Юаня и Бертрама [5] 1991 года. Было установлено, что перестройка структуры асимметричной блоховской стенки под действием внешнего поля, параллельного ОЛН, заключается в смещении внутристеночного вихря из центральной части к поверхности пленки. Как видно из рис. 2а, исходное распределение M в асимметричной блоховской стенке таково, что результирующая проекция намагниченности ДГ в направлении x равна нулю. При смещении вихря от центральной области к поверхности пленки симметрия стенки относительно плоскости y = 0 нарушается (см. рис.3в). У поверхности пленки, противоположной той, к которой сместилась вихреподобное образование, возникает преимущественная ориентация намагниченности вдоль положительного или отрицательного направления x. Направление смещения вихря (к нижней поверхности пленки или к верхней) определяется направлением внешнего поля и положением вихря относительно центральной линии стенки. В полях H>Hc стенка перестраивается в асимметричную неелевскую. Наклон центральной линии последней в процессе движения постепенно изменяется на противоположный, затем у поверхности пленки, противолежащей поверхности, к которой происходило смещение ранее существовавшего вихря, зарождается новый вихрь. Данный вихрь затем смещается к центру пленки. Далее ДГ снова проходит все описанные выше стадии, но с противоположной закрученностью вихрей, и таким образом возвращается к первоначальной конфигурации. Как и в случае одномерной стенки, периодическая перестройка двухмерной ДГ ведет к периодическому изменению скорости ее движения со временем (рис. 3).

Рис. 3. Зависимость усредненной по толщине пленки скорости движения асимметричной блоховской ДГ в пленке толщиной 50 нм от времени в поле H>Hc и соответствующая динамическая перестройка структуры ДГ.

Рис. 4. Динамическая перестройка структуры ДГ при ее движении в поле H>Hc в пленках толщиной 70 нм (а) и 100 нм (б) [121].

Из сказанного выше ясно, что процесс движения ДГ в магнитной пленке, описываемый в рамках двухмерной модели распределения намагниченности, является более сложным по сравнению с одномерной моделью. Включение внешнего магнитного поля, направленного вдоль ОЛН, в данном случае ведет не только к смещению стенки, но и к смещению вихрей намагниченности внутри нее.

В настоящее время наблюдается всплеск интереса к исследованию динамики ДГ в пленках с плоскостной анизотропией, связанный с разработкой новых  типов устройств магнитной памяти c ультра высокой плотностью записи информации, работа которых основана на использовании спин-поляризованных токов. Появился целый ряд работ, в которых исследуется динамика стенок в тонких пленках-полосках, которые по сути требуют уже 3-мерных областей расчёта.

Механизм движения ДГ в тонких пленках при больших полях остается не ясным. Поэтому 3-мерная визуализация геометрии распределений намагниченности при динамических перестройках была бы весьма полезной.


3. Математическая постановка задачи

3.1. Применяемый метод микромагнитного моделирования

Рассмотрим ферромагнитную пленку толщиной b и будем считать, что ось одноосной магнитной анизотропии лежит в плоскости пленки. Направим ось z прямоугольной системы координат вдоль оси анизотропии, а ось y – вдоль нормали к поверхности пленки. Магнитное состояние пленки предполагается соответствующим двум доменам, в которых намагниченность направлена вдоль ± z (рис. 1), разделенным 180° доменной границей.

Предполагается, что для линейных размеров пленки выполняются следующие соотношения: , то есть пленка считается бесконечно протяженной в плоскости xz. Тогда, исключая из рассмотрения возможное образование блоховских линий, можно считать, что внутри ДГ направление вектора намагниченности не зависит от координаты  z, то есть:

,                                                             (3)

Выражение (3) соответствует двухмерной модели распределения намагниченности в ДГ. В рамках этой модели функционал полной энергии ДГ, рассчитанной на единицу длины стенки вдоль оси z,  может быть представлен в виде:

        (4)

где V – некоторая область пленки, где сосредоточена ДГ, S – участки поверхности пленки, связанные с данной областью. Слагаемые в скобках под знаком интеграла по объему обозначают соответственно плотности обменной энергии, энергии анизотропии, энергии взаимодействия намагниченности стенки с внешним магнитным полем  и собственной магнитостатической энергии стенки. Выражение под знаком интеграла по поверхности учитывает поверхностную анизотропию. A – обменный параметр, MS – намагниченность насыщения, K – константа одноосной анизотропии, KS – константа поверхностной анизотропии, c – единичный вектор, сонаправленный с ОЛН, n – единичный вектор нормали к поверхности пленки, - напряженность магнитостатического поля, определяемая согласно

,      .                                                  (5)

где  φ - магнитостатический потенциал. Если не оговорено особо, при решении задачи о статическом распределении намагниченности в ДГ внешнее поле Н считается равным 0. При  (см. выше)  не зависит от z и из выражения (4) может получено следующее выражение для энергии ДГ в расчете на единицу площади ее боковой поверхности (плоскость yz):

                 (6)

где D – сечение области V плоскостью xy, а Dx – область точек поверхности пленки, принадлежащих D.

Вид распределения намагниченности внутри ДГ определяется минимизацией (6). В рамках двухмерной модели точно решить данную задачу аналитически не представляется возможным. В связи с этим, используется метод численной минимизации . Для выполнения минимизации приходится считать, что область D имеет конечные размеры. Пусть D имеет форму прямоугольника с размерами a вдоль оси x и b вдоль оси y. Далее область D, которую будем называть также расчетной областью, разбивается сеткой  на квадратные ячейки (см. рис.5), где Nx и Ny- число ячеек вдоль осей x и y, соответственно. При этом область V оказывается разбитой на параллелепипеды с линейными размерами  в плоскости xy и  вдоль оси z. Предполагается, что ячейки имеют размеры, достаточные для использования микромагнитного подхода, однако настолько малые, что в каждой ячейке направление m можно считать постоянным. При этом вдоль оси z намагниченность каждого из параллелепипедов остается постоянной в силу двухмерности модели. Таким образом, ориентация m в D меняется при переходе от ячейки к ячейке.

При вычислениях учитываются следующие условия для границ расчетной области:

,                                                (7)

,                                                                       (8)

Условие (8) обозначает, что на границах расчетной области, примыкающих к доменам, направление m совпадает с направлением намагниченностей в доменах. В особом случае, когда внешнее поле имеет компоненту в направлении, нормальном к ОЛН, компоненты Mx, My намагниченности доменов вычисляются минимизацией (6) при отсутствии в (6) обменной энергии. Выражение (7) следует из условия минимума (6) и при сводится к , то есть к незакрепленности намагниченности на поверхностях пленки. Если  КS  0  (КS  0) , то имеется поверхностная магнитная анизотропия типа ось (плоскость) легкого намагничивания, причем ось анизотропии сонаправлена с n. Вклад в граничные условия от поверхностной анизотропии, как показано в Приложении 2, так же можно описать с помощью эффективного поля , действующего на поверхностях пленки.

Рис. 5. Схематичная иллюстрация процедуры дискретизации расчетной области.

Далее производится дискретизация слагаемых плотности энергии, входящих в (6) и  записывается в виде суммы по всем ячейкам сетки. Ячейке сетки с индексами  соответствуют величины ,  и , связанные соотношением .

Энергия анизотропии с соответствии с сеточным приближением тогда может быть представлена в виде

.                                      (9)

Плотность энергии обменного взаимодействия между ячейками (I,J) и (I+1,J), имеет вид (см. [42]):

     (10)

(получено из выражения (4) переходом от дифференцирования к конечным разностям).

Аналогичное выражение (с соответствующими индексами) получается для ячеек (I,J+1) и (I,J).

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

            (11)

Здесь дополнительно введены ячейки с индексами  и , а так же  и  для удовлетворения граничным условиям (7) и (8). О слагаемых вида  и  будет сказано ниже в Приложении 2. Что касается  и , то согласно (8) они сонаправлены с намагниченностью левого и правого доменов соответственно.

При получении выражения для магнитостатической энергии, как уже упоминалось выше, используются известные выражения (5), которые следуют из решения уравнений магнитостатики

,                            (12)

с соответствующими граничными условиями (непрерывность нормальной составляющей индукции  и тангенциальной составляющей напряженности магнитостатического поля ).

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

.              (13)

Процедура дискретизации выражения (13) ввиду ее громоздкости вынесена в Приложение 1. Конечный результат выглядит следующим образом:

     (14)

где ,

,

,

.

Таким  образом, в результате процедуры дискретизации функционал  оказывается определенным в  мерном пространстве компонент векторов m(I, J), где I и J – индексы ячейки вдоль x и y соответственно,  . 

Минимизация  с использованием того или иного метода (в используемой в настоящей работе программе реализован метод градиентного спуска) дает энергию  стабильной или метастабильной стенки и соответствующее ей распределение m.

Решение задачи о движении доменной стенки производится на основе прямого численного интегрирования уравнения Ландау-Лифшица с диссипацией в форме Гильберта:

                                                           (18)

,

.

Здесь γ – гиромагнитное соотношение, e – абсолютная величина заряда электрона, m – масса электрона, c – скорость света в вакууме. При этом используется та же вышеописанная пространственная сетка и те же условия на поверхности пленки, что и при нахождении статического распределения намагниченности в ДГ. Внешнее поле  направляется вдоль ОЛН, то есть параллельно намагниченности одного из доменов. Удобно записать уравнение движения в безразмерном виде:

                                      ,                           (19)

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

Для решения (19) был использован метод Эйлера, модифицированный применением предиктора-корректора. В соответствии с этим методом в момент  =0 задается распределение намагниченности m0. На первом этапе итерация mn+1 определяется согласно формуле

                                                                                                (20)

(предиктор), где

                          .                 (21)

На втором этапе применяется процедура окончательного определения mn+1 

                                                                                         (22)

(корректор).

3.2. Описание программы

Программный пакет DOMEN для исследования структуры и динамики ДГ в магнитных пленках на основе описанных выше методов был разработан
в лаборатории микромагнетизма Института физики металлов УрО РАН.

В результате выполнения программы происходит расчет распределения намагниченности в ДГ для заданных параметров магнитной пленки и соответствующего этому распределению значения полной энергии стенки . Полученные данные отображаются на мониторе, примером рассчитанного распределения M, может служить рис.2. На рис.2 и всех приводимых далее аналогичных рисунках, отображающих результаты расчетов, стрелками показаны проекции M на плоскость xy. Длина стрелок пропорциональна величине указанной проекции. Слева от центральной линии стенки (центральная пунктирная линия) Mz отрицательно, справа – положительно. Пунктирные линии на рис. 1 и рисунках, иллюстрирующих распределение M в ДГ, представляют собой линии уровня Mz, определенные следующим образом:

                   (23)

Ясно, что k = 0 соответствует центральной линии стенки. Если не оговорено особо, будут изображаться только три первые линии, соответствующие k = 0, 1. Изменение Mz между линиями k =  1 равно приблизительно 48,4% полного изменения Mz.

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

Важным моментом является выбор ширины расчетной области a. Увеличение a ведет к уточнению  и соответствующего распределения намагниченности, но так же приводит и к увеличению времени расчета. Необходимо выбирать некоторый оптимальный размер a, при котором дальнейшее его увеличение не ведет к изменению . Различным распределениям M (то есть различным типам ДГ) может соответствовать свой оптимальный размер расчетной области. Например, если для асимметричной блоховской стенки в однородной пленке толщиной 50 нм оптимальным является соотношение a/b=3, то для симметричной неелевской стенки в 20 нм пленке использовались соотношения a/b порядка 10. Существенное сокращение времени расчета может быть достигнуто подходящим выбором исходной конфигурации намагниченности.

Так же большое значение имеет выбор размеров ячейки сетки. Уменьшение этих размеров увеличивает точность расчета, но и резко увеличивает его время. Может быть эмпирически подобран оптимальный размер ячеек, когда изменение энергии при дальнейшем измельчении сетки не превосходит нескольких процентов. Максимальная размерность сетки, используемая в данной работе, составляла 240×60.

Полученное в результате расчета распределение намагниченности может быть сохранено в виде файла, содержащего массив значений компонент вектора m(I, J.

При исследовании движения доменной стенки под действием внешнего поля, направленного вдоль ОЛН, вычисляется и выводится на экран распределение M в каждый момент расчетного времени, а также график временной зависимости мгновенной скорости ДГ. Указанные данные сохранялись и использовались для определения периода динамической перестройки ДГ при ее движении в полях H>Hc, максимальной и средней по времени скоростей движения ДГ, исследования механизма динамической перестройки структуры ДГ. Под мгновенной скоростью v  ДГ понимается скорость движения «центра тяжести стенки» - точки максимального градиента намагниченности. Средняя по времени скорость ДГ определялась графическим методом: на рассчитанной кривой v(t) производилась серия отсчетов с частотой 0.01-0.05 нс за один или два периода осцилляций и вычислялось среднее значение vav. Программа для определения vav указанным методом была разработана автором данной работы. Использовались магнитные параметры, характерные для пермаллоевых пленок.


4. Обзор вариантов 2-мерной визуализации

Визуализация магнитных расчётов весьма важна в связи со сложной структурой вектора намагниченности, без понимания которой количественные результаты теряют смысл. Рассмотрим примеры визуализации по работам 80-ых и 90-ых годов прошлого века.

1.

Один из ранних примеров использования стрелок для 2D-визуализации магнитного момента [6]. Как видим, качество и понимаемость визуализации оставляют желать лучшего.

2.

Сплошные линии относятся к смене знака Z-компоненты, то есть отражают границы между доменами.

Интересен ряд стрелок над областью расчёта, который изображает проекцию векторов намагниченности на плоскость плёнки. При этом стрелки в области расчёта изображают проекцию векторов на плоскость стенки. Таким образом, здесь мы видим попытку 3-мерной визуализации при ограниченных изобразительных средствах [7].

Другой подход к 2-мерной визуализации – линии уровня Z-компоненты в области её разворота между доменами [7].

3.

Другой способ «отрёхмеривания» [8]  концы стрелок представляют собой 3-мерные пирамидки, позволяющие лучше понять направление векторов в 3-мерном пространстве.

4.

С помощью стрелок отображаются вихревые структуры [9].

5.

Стрелки используются для изображения вектора намагниченности в сложной области на нерегулярной сетке [10].

6.

Попытка [11] использования в магнитных расчётах известных пакетов
3
D-визуализации. Поскольку практически все такие пакеты «заточены» на изображение потоков жидкости и газа (см. следующий раздел), то в данном случае их полезность неочевидна. Например, не ясно, имеют ли жгуты линий тока на изображениях какой-либо физический смысл.


5. Обзор подходов к 3-мерной визуализации векторных полей

Поиск публикаций на тему 3D-визуализации показал, что в подавляющем большинстве случаев под векторным полем понимается поле течения жидкости или газа. Соответственно, опубликованные методы 3D-визуализации «заточены» на изображение непрерывных потоков, в первую очередь на построение линий тока.

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

Ниже приведён обзор 8 статей на указанную тему. Приводятся примеры изображений, получаемых авторами с применением разработанных программных средств 3D-визуализации.

1.

Representation of Vector Fields

Alexandru Telea,  Jarke J. van Wijk

Визуализация потоков [12], к тому же 2-мерных, 1994-1997, хотя в конце есть картинки и для 3-мерных течений.


2.

Higher Dimensional Vector Field Visualization: A Survey

Zhenmin Peng and Robert S. Laramee

Тоже «flow vizualisation» [13], 2009

Пример – поток газа, истекающего из шара.


3.

An Illustrative Visualization Framework for 3D Vector Fields

Cheng-Kai Chen, Shi Yan, Hongfeng Yu, Nelson Max, Kwan-Liu Ma

Линии тока [14], их кластеризация и раскраска для течений жидкости или газа

A – набор линий тока

B – выделение основных типичных линий тока (кластеризация)

C – толщина и раскраска линий связаны со скоростью и интенсивностью потоков

Изображение урагана.

4.

Symmetric Tiling of Closed Surfaces: Visualization of Regular Maps

Jarke J. van Wijk

Несколько особняком стоит работа [15], посвящённая визуализации раскраски карт на сложных 3-мерных поверхностях

Раскраска карт на 3-мерных поверхностях

Обработка данных и результат для сложной поверхности


5.

Image Based Flow Visualization

Jarke J. van Wijk

Разные методы [16] визуализации для потоков, 2002

Анимация турбулентного течения.

     

  

     

на картинках – разные подходы и этапы визуализации турбулентного течения.

6.

Interactive Visualization Of 3D-Vector Fields Using Illuminated Stream Lines

Malte Zockler, Detlev Stalling, Hans-Christian Hege

Опять построение интегральных кривых == линий тока [17].


7.

Comparing 3D Vector Field Visualization Methods: A User Study

Andrew S. Forsberg, Jian Chen, David H. Laidlaw

Линии тока с указанием направления и раскраской в соответствии со скоростью потоков [18].


8.

Parallel Hierarchical Visualization of Large Time-Varying 3D Vector Fields

Hongfeng Yu, Chaoli Wang, Kwan-Liu Ma

Тоже течения

Стадии восстановления линий тока по векторному полю [19].

Итак, видно, что существующие программы визуализации 3-мерных векторных полей не вполне подходят для отображения распределения намагниченности.

Поэтому возникла необходимость разработки собственного модуля
3
D-визуализации для рассматриваемой задачи.


6. Описание модуля 3D-визуализации

6.1. Инструменты

Существует много пакетов и движков для рисования трёхмерной графики. Например, Irrlicht, OGRE, OpenSceneGraph, VL, GLScene и, самые известные, DirectX и OpenGL. Так как данный проект, возможно, будет использоваться не только под Windows, то хорошо было бы иметь кроссплатформенность, то есть сразу отпал вариант DirectX. В движках (Irrlicht, OGRE…) много ненужного данному проекту функционала. Нужен был простой для понимания код, без лишних деталей, так что выбор пал на OpenGL и его вспомогательную библиотеку GLUT (OpenGL Utility Toolkit). GLUT, в основном, отвечает за системный уровень операций ввода-вывода при работе с операционной системой (создание окна, управление окном, отслеживание событий мыши и клавиатуры), но также включает в себя функции рисования некоторых графических примитивов (куб, сфера, конус, чайник…). Язык C++ был выбран для упрощения взаимодействия с проектом по расчету векторов, который был написан на C.

6.2. Структура программы

  1.  Инициализация
  2.  Функция чтения данных
  3.  Функции для обработки сигналов с клавиатуры и мыши
  4.  Вспомогательные функции рисования
  5.  Основная функция рисования


6.3. Инициализация

int main(int argc, char** argv){

glutInit(&argc, argv);

glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);

glutInitWindowSize(900, 650);  //начальный размер окна

glutInitWindowPosition(10,330);  //начальная позиция окна

glutCreateWindow("Magnetic Field Vector Visualisation");

//создание окна с заданным именем

init();  //инициализация

glutReshapeFunc (reshape);  //обработчик изменения

// размера окна

glutKeyboardFunc (keyboard);  //обработчик сигналов

// с клавиатуры(буквы, цифры)

glutSpecialFunc(keyboardSpecial);//обработчик сигналов с

//клавиатуры(стрелки, шифт)

glutDisplayFunc (display);  //основная функция рисования

glutMotionFunc(MouseMotion); //обработчик движения мыши

glutMouseFunc (MouseButton); //обработчик нажатия клавиш мыши

glutTimerFunc(100,timf,0);  //таймер

glutMainLoop(); //запуск цикла обработки событий GLUT

 return 0;

}

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

6.4. Чтение данных

void READ_SMOMENT( char *fname ){

 unsigned char mom_byte[3];

 double alf, bet, gam;

 unsigned char gam0;

 if(mlt==NULL)

 mlt=fopen(fname, "rb" );

fseek(mlt,92+iter*(960),SEEK_SET);

 for(int i=0; i<Nx; i++ ){

 for(int j=0; j<Ny; j++ ){

  fread( mom_byte, 1, 3, mlt );

  // alf

 alf = ((unsigned char)(mom_byte[2]*4) +

((mom_byte[1]>>7)&1)*2) / 2 ;

  if((mom_byte[2]>>6)&1)  alf -= 128 ;

  alf /= MLT_ALF ;

  //gam

  gam0 = mom_byte[1] * 4 ;

  gam  = mom_byte[0] + (int)gam0 * 64 ;

  if((mom_byte[1]>>6)&1) gam -= 16384 ;

  gam /= MLT_GAM ;

  //bet

  bet = sqrt( 1 - alf * alf - gam * gam );

  if((mom_byte[2]>>7)&1)  bet = -bet ;

  moment[i][j][0] = alf ;

  moment[i][j][1] = bet ;

  moment[i][j][2] = gam ;

 }

}

}

При сохранении данных в проекте, рассчитывающем вектора, использовались битовые поля для экономии места. Один вектор занимает всего 3 байта(15 бит на X-компоненту, 1 бит на Y-компоненту, 8 бит на Z-компоненту). Известно, что длина каждого вектора равна единице, так что можно восстановить Y-компоненту.

6.5. Обработка сигналов с клавиатуры и мыши

Стрелки – движение камеры.

Стрелки(с зажатым шифтом) – поворот камеры.

N – переход к следующей итерации

B – переход к предыдущей итерации

Движение мыши(с зажатой левой кнопкой) – движение камеры вперед/назад/влево/вправо

Движение мыши(с зажатой правой кнопкой) – поворот камеры

Движение мыши(с зажатой средней кнопкой) – движение камеры вверх/вниз

void MouseMotion(int x,int y){

 if(g_bButton1Down){

 cameraPosX+=(-sin(pi*cameraAngleH/180)*(float)(mousePosY-y) + cos(pi*cameraAngleH/180)*(float)(mousePosX-x))/100;

 cameraPosZ+=(cos(pi*cameraAngleH/180)*(float)(mousePosY-y) + sin(pi*cameraAngleH/180)*(float)(mousePosX-x))/100;

}

 else if(g_bButton2Down){

 cameraAngleH-=(GLfloat)(mousePosX-x)/5;

 cameraAngleV-=(GLfloat)(mousePosY-y)/5;

}

 else if(g_bButton3Down)

 cameraPosY-=(GLfloat)(mousePosY-y)/50;

glutPostRedisplay();

 mousePosX=x;

mousePosY=y;

}

В движении учитываются углы поворота камеры, таким образом, например, при перемещении мыши вверх, камера будет двигаться по направлению взгляда, а не просто параллельно/перпендикулярно плоскости векторов.


F – переключение режима тумана между GL_LINEAR, GL_EXP, GL_EXP2 и отключенным. В этих режимах цвет тумана смешивается с изображением следующим образом.

Рис. 6.1 Распределение доли цвета тумана в изображении в зависимости

от расстояния до камеры.(GL_FOG_DENSITY=0.05,GL_FOG_START=10, GL_FOG_END=20)


W – увеличение FOV(Field of View) камеры

S – уменьшение FOV камеры

  

Рис. 6.2 Влияние FOV на изображение.

(слева – 30, посередине – 60, справа – 90)

I,J,K,L – управление положением точечного источника света

U – восстановление значения по умолчанию для положения источника света

V – переключение отображения сферы на месте источника света

   

Рис. 6.3 Различные примеры расположения источника

света и отображения сферы на его месте.

z – уменьшение радиуса основания конуса

Z – увеличение радиуса основания конуса

x – уменьшение длины конуса

X – увеличение длины конуса

C – изменение вида основания конуса (сфера/круг)

   

Рис. 6.4 Различные примеры возможных форм конусов.

6.6. Вспомогательные функции рисования

renderSphere – используется для отображения сферы на месте источника света

renderCone – используется для отображения конуса на месте вектора

static void renderCone (GLfloat x, GLfloat y, GLfloat z,GLfloat VecX, GLfloat VecY, GLfloat VecZ){

glPushMatrix();

glTranslatef (x, y, z);

glRotatef(

acos(VecZ/sqrt(VecX*VecX+VecY*VecY+VecZ*VecZ))*180/pi,

VecY,-VecX,0);

glutSolidCone(coneRad,coneLen,16,1);

 if(baseSphere)

 glutSolidSphere(coneRad,16,4);

 else{

 glBegin(GL_POLYGON);

 glNormal3f(0,0,-1);

 for(int i=16;i>0;i--){

  glVertex3d(coneRad*cos(i*(2*pi/16)),

coneRad*sin(i*(2*pi/16)),0);

 }

 glEnd();

 }

glPopMatrix();

}

Функция glutSolidCone(rad, len, sl, st) создает конус направленный по вектору (0,0,1) в точке (0,0,0) текущей системы координат. Для того, чтобы создать его в определенном месте, направленным по определенному вектору, требуется передвинуть и повернуть систему координат. Это делается с помощью преобразований матрицы модели(GL_MODELVIEW_MATRIX).

Чтобы работа данной функции не влияла на работу других, текущая матрица модели кладется в стек в начале работы функции командой glPushMatrix() и берётся из стека в конце работы функции командой glPopMatrix().

Команда glTranslatef(x, y, z) сдвигает текущую матрицу на вектор (x, y, z)

путём умножения текущей матрицы на матрицу .Точка (0,0,0) в полученной системе координат является точкой (x, y, z) в старой системе.

Команда glRotatef(angle, x, y, z) поворачивает текущую матрицу на угол angle вокруг вектора (x, y, z) путём умножения текущей матрицы на матрицу поворота.

Матрица поворота вычисляется следующим образом.

  

Угол поворота вычисляется из скалярного произведения вектора (0, 0, 1) и (VecX, VecY, VecZ), по формуле .

Вектор, вокруг которого производится вращение, является векторным произведением тех же векторов.

6.7. Основная функция рисования

void display(){

 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

 

 //проекция

glMatrixMode(GL_PROJECTION);

glLoadIdentity();

glGetIntegerv( GL_VIEWPORT, m_viewport );

gluPerspective(fov,(GLdouble)m_viewport[2]/m_viewport[3],

0.1, 50);

 

 //камера

glMatrixMode(GL_MODELVIEW);

glLoadIdentity();

glRotatef(cameraAngleV,1,0,0);

glRotatef(cameraAngleH,0,1,0);

glTranslatef(cameraPosX,0,cameraPosZ);

 

 //свет

 if(showLight)

 renderSphere(light0_pos[0],light0_pos[1],light0_pos[2]);

glLightfv(GL_LIGHT0, GL_POSITION, light0_pos);

 //векторное поле

 for(int i=0;i<30;i++){

 for(int j=0;j<10;j++){

  if(moment[i][j][2]<0){

   matm[0]=0.8-powf(-moment[i][j][2],6)*(0.7);

   matm[2]=0;

  }else{

   matm[0]=0;

   matm[2]=1-powf(moment[i][j][2],6)*(0.9);

  }

  matm[1]=0;matm[3]=1;

  glMaterialfv (GL_FRONT,

GL_AMBIENT_AND_DIFFUSE, matm);

  renderCone(i-15,j-5,0,

   -moment[i][j][0],-moment[i][j][1],moment[i][j][2]);

 }

}

 //сетка

glLineWidth(1);

matm[0]=matm[1]=matm[2]=0;  matm[3]=1;

glMaterialfv (GL_FRONT, GL_AMBIENT_AND_DIFFUSE, matm);

matm[0]=matm[1]=matm[2]=0.3;matm[3]=1;

glMaterialfv (GL_FRONT, GL_EMISSION, matm);

glEvalMesh2(GL_LINE, 0, 18, 0, 58);

matm[0]=matm[1]=matm[2]=0;matm[3]=1;

glMaterialfv (GL_FRONT, GL_EMISSION, matm);

 //туман

GLfloat fogSt=sqrt(cameraPosX*cameraPosX +

cameraPosZ*cameraPosZ);

glFogf (GL_FOG_START, fogSt-15);

glFogf (GL_FOG_END, fogSt+15);

glutSwapBuffers();

glFlush();

}

Команда gluPerspective(fov, aspect, zNear, zFar) генерирует матрицу 

, которая умножается на текущую матрицу проекции, которая в данном случае является единичной после действия команды glLoadIdentity(). В результате проекция происходит следующим образом.

glClear(COLOR_BUFFER_BIT, DEPTH_BUFFER_BIT) – очищает буфер цвета и буфер глубины, устанавливая цвет равным значению, переданному ранее функции glClearColor(R, G, B, A), а глубину равной максимально возможному значению.

glEvalMesh2(GL_LINE, 0, 18, 0, 58) – создает решетку из линий, основанную на координатах, определенных при инициализации.

6.8. Примеры получаемых 3D-изображений

Асимметричная блоховская стенка

Асимметричная неелевская стенка


Заключение

Автором данной работы разработана среда научной визуализации для представления процессов намагничивания.

Здесь этот подход применяется для визуализации 3-мерных векторов намагниченности на 2-мерной расчётной области. Однако очевидно, что разработанный модуль применим для визуализации любых векторных полей, определённых на 2-мерной или 3-мерной области. При этом область может иметь сложную форму, сетка на области не обязана быть регулярной.

Доработки:

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

Перспективы:

Программа 3D-визуализации будет ещё более полезна при переходе
к 3-мерным расчётам распределения вектора намагниченности.

Возможно использование других средств визуализации. В частности, построение 3-мерных поверхностей для скалярных величин (например, Z-компоненты намагниченности, компонент ротора поля). Такие поверхности могут раскрашиваться в соответствии с линиями уровня. Возможно вращение, приближение и удаление изображения поверхности. Изменение изображения может сопровождать расчёт  поля в режиме реального времени.


Литература

1. Landau L.D., Lifshitz E.M. Theory of dispersion of magnetic permeability in ferromagnetic bodies  // Phys. Z. Sowjet. 1935. V. 8. P.155-171.

2. La Bonte A.E. Two-dimensional Bloch-type domain wall in ferromagnetic films // J. Appl. Phys. 1969. V. 40. № 6. P. 2450-2458.

3. Hubert A. Stray-field free magnetization configurations // Phys. Stat. Sol. (a). 1969. V. 32. № 2. P. 519-534.

4. Schryer N.L., Walker L.R. The motion of 1800 domain walls in uniform dc magnetic fields // J. Appl. Phys. 1974. V. 45. № 12. P. 5406-5421.

5. Yuan S.W., Bertram H.N. Domain wall structures and dynamics in thin films // Phys. Rev. B. 1991. V. 22. P. 12395-12405.

6. M. Mansuripur and R. Giles, "Demagnetizing field computation for dynamic simulation of the magnetization reversal process," IEEE Trans. Magn., vol. 24, no. 6,

pp. 2326-2328, Nov. 1988.

7. M. Labrune and J. Miltat: "Micromagnetics of strong stripe domains in NiCo-films" IEEE Trans. Mag. 26, 1521-1523 (1990)

8. Andrew J. Newell and Dawid J. Dunlop A two dimensional Micromagnetic Model of Magnetizations and field in Magnetite. J. of Geophisical Research, Vol. 98, No B6, P. 9533-9549 (1993)

9. S.W. Yuan and H.N. Bertram, "Size effects of switching fields of thin Permalloy particles," Paper JD-02, Intermag92, St. Louis, Missouri.

10. Fredkin, D.R., and T.R. Koehler, Numerical micromagnetics by the finite element method, IEEE Trans. Magn., tS, 3385-3387, 1987.

11. Richard P. Boardman, Jurgen Zimmermann, Hans Fangohr, Alexander A. Zhukov, Peter A. J. de Groot   Micromagnetic simulation studies of ferromagnetic part spheres // J. Appl. Phys. 97 10E305 (2005), 3 pages.

12. Alexandru Telea, Jarke J. van Wijk   Representation of Vector Fields

13. Zhenmin Peng and Robert S. Laramee   Higher Dimensional Vector Field Visualization: A Survey

14. Cheng-Kai Chen, Shi Yan, Hongfeng Yu, Nelson Max, Kwan-Liu Ma   An Illustrative Visualization Framework for 3D Vector Fields

15. Jarke J. van Wijk   Symmetric Tiling of Closed Surfaces: Visualization of Regular Maps

16. Jarke J. van Wijk   Image Based Flow Visualization

17. Malte Zockler, Detlev Stalling, Hans-Christian Hege   Interactive Visualization Of 3D-Vector Fields Using Illuminated Stream Lines

18. Andrew S. Forsberg, Jian Chen, David H. Laidlaw   Comparing 3D Vector Field Visualization Methods: A User Study

19. Hongfeng Yu, Chaoli Wang, Kwan-Liu Ma   Parallel Hierarchical Visualization of Large Time-Varying 3D Vector Fields


 

А также другие работы, которые могут Вас заинтересовать

47003. Причины первой мировой войны 43.5 KB
  Покровский Стремившийся доказать что царская Россия есть главная виновница войны Покровский в то же время признавал что в основе международных противоречий вызвавших империалистическую войну лежал англогерманский конфликт за которым по значению следовал германофранцузский. Гуч в книге Накануне войны вышедшей в 1938 г. предпринял попытку доказать что возникновение войны в 1914 г.
47006. ТЕОРИЯ ПОДОБИЯ В ПРИМЕНЕНИИ К ДИФФЕРЕНЦИАЛЬНОМУ УРАВНЕНИЮ ТЕПЛОПРОВОДНОСТИ 43.9 KB
  Считаем также что начальная температура тела одинакова и не зависит от координат т.4 где α – коэффициент теплоотдачи от тела к омывающей среде Tw температура стенки тела . С другой стороны плотность теплового потока у стенки тела равна: ∂T ∂ϑ qw = −λ = −λ 4.5 ∂n ∂n w w где λ коэффициент теплопроводности тела ∂T производная температуры в теле по нормали к поверхности.
47009. СКРЫТЫЕ (УСЛОВНЫЕ) БАЗЫ 44.26 KB
  Применение условных скрытых баз при проектировании тем более удобно что позволяет исключить из расчетов неизбежные погрешности реальных поверхностей снижающие точность базирования. При базировании деталей собираемых узлов и обрабатываемых заготовок в подавляющем большинстве случаев используются материальные поверхности явные базы по ГОСТ 21495 76 однако и в этом случае для повышения точности базирования иногда применяются условные скрытые базы материализуемые различными устройствами отвесы коллиматоры центрирующие...