Influence of internal structures on the kinetics of magnetization reversary of ferrofluids

Full Text

ВВЕДЕНИЕ

Магнитные жидкости (феррожидкости) – коллоидно-стабилизированные суспензии однодоменных ферро- или ферримагнитных частиц в жидких средах. Эти многофункциональные системы находят применение во многих высокотехнологических приложениях, благодаря богатому набору их уникальных физических свойств. Обзоры работ по физике феррожидкостей и их практическим применениям можно найти, например, в [1–8].

Одним из перспективных направлений применений магнитных жидкостей является магнито-гипертермическая терапия онкологических заболеваний (см., например, [9–12]). Этот метод основан на внедрении магнитных наночастиц, покрытых специальными биоактивными слоями, в диагностированную опухолевую область и в воздействии на эту область переменным магнитным полем. Под действием поля магнитные наночастицы и, следовательно, захватившие их опухолевые клетки разогреваются. Если температура больной клетки превышает критическую температуру, обычно оцениваемую как 41–42^оC, белок в этой клетке денатурирует, и она гибнет. В то же время здоровые клетки остаются неповрежденными до существенно более высоких (примерно, на 10^оC) температур. Кинетика перемагничивания внедренных феррочастиц играет ключевую роль в генерации тепла этими частицами.

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

Первые теории кинетики перемагничивания феррожидкостей были предложены М.И. Шлиомисом с сотрудниками; их обсуждение можно найти, например, в [3, 4]. Эти модели имеют дело с одиночными, невзаимодействующими частицами, что оправдано применительно к системам с очень низкой их концентрацией. Однако многочисленные лабораторные и компьютерные эксперименты (см., например, обзоры в [4, 5]) показывают, что межчастичные взаимодействия могут существенно влиять на равновесные и динамические свойства этих систем. Теоретические модели кинетики перемагничивания феррожидкостей с взаимодействующими, но одиночными (не агрегированными) частицами были предложены в работах [13, 14]. Несмотря на различие используемых математических подходов, результаты этих моделей о влиянии межчастичных взаимодействий на темп перемагничивания феррожидкостей близки. В частности, по порядку величины время перемагничивания феррожидкостей с взаимодействующими, но не агрегированными частицами совпадает со временем перемагничивания феррожидкостей с одиночными частицами.

Хорошо известно, что магнитное взаимодействие наночастиц может вызвать их объединение в гетерогенные структуры – линейные цепочки, кольца, ветвистые агрегаты, объемные «капли» и т. д. [4, 5]. Магнитные моменты частиц в этих структурах, как правило, сильно коррелированы. Поэтому статистически средние моменты структурированных частиц могут быть существенно выше, чем у одиночных частиц [5, 15, 16]. Следовательно, кинетика перемагничивания феррожидкостей с частицами, способными к агрегированию, может определяться характерным временем эволюции этих структур, а не временем переориентации момента одиночной частицы.

В этой работе мы представляем результаты компьютерного моделирования кинетики перемагничивания феррожидкости, частицы которой, благодаря их магнито-дипольному взаимодействию, могут объединяться в гетерогенные структуры. Компьютерное моделирование выполнено на основе пакета программ ESPResSo v4.2.1 [17] и языке программирования Python. Пакет ESPResSo основан на численном решении системы уравнений движения частиц. Учитываются магнитные и стерические взаимодействия всех частиц со всеми.

МАТЕМАТИЧЕСКАЯ И КОМПЬЮТЕРНАЯ МОДЕЛЬ

Компьютерная модель феррожидкости состоит из N сферических частиц с одинаковым диаметром d, однородно распределенных в кубической ячейке с длинной ребра L. Задав объемную концентрацию φ частиц, мы определяем L:

$L=d{\left(\frac{\pi N}{6\phi }\right)}^{\frac{1}{3}}.$ (1)

Периодические условия на границах симуляционной ячейки задаются вдоль всех координатных осей x, y и z. Предполагается, что каждая частица обладает магнитным моментом с абсолютным значением m; магнитный момент вморожен в тело частицы, т. е. поворачивается вместе с ним. Как известно, это приближение оправдано, если размер частицы превышает некоторое критическое значение, зависящее от ее материала (см., например, [3, 4]). В свою очередь, относительно крупный размер частиц означает, что магнитное межчастичное взаимодействие между ними может быть вполне существенным.

Потенциальная энергия диполь-дипольного взаимодействия между i-й и j-й частицами имеет вид:

$U_{\text{DD}}=-{\mathrm{\mu }}_0\frac{m^2}{4\pi }\frac{3\left(e_{\text{i}}\cdot r_{\text{ij}}\right)\left(e_{\text{j}}\cdot r_{\text{ij}}\right)-r_{\text{ij}}^2\left(e_{\text{i}}\cdot e_{\text{j}}\right)}{r_{\text{ij}}^5}.$ (2)

Здесь μ₀ – магнитная проницаемость вакуума; e_i и e_j – единичные векторы, направленные вдоль магнитных моментов частиц; r_ij – радиус-вектор, соединяющий центры этих частиц.

Чтобы имитировать стерическое отталкивание частиц, мы используем потенциал Вика-Чандлера-Андерсона U_WCA [18], который представляет собой потенциал Леннард-Джонса с обрезанным радиусом действия r_cut:

$\begin{array}{l}{U}_{\text{WCA}}=\left\{\begin{array}{c}{\mathrm{0,}}_{} {\text{если}}_{} r_{\text{ij}}>r_{\text{cut}}\\ 4\epsilon \left[{\left(\frac{d}{r_{\text{ij}}}\right)}^{12}-{\left(\frac{d}{r_{\text{ij}}}\right)}^6\right],{ }_{}\text{если}{ }_{}{r}_{\text{ij}}<r_{\text{cut}}\end{array}\right., \\ r_{\text{cut}}=2^{\frac{1}{6}}d.\end{array}$(3)

Здесь ε – параметр стерического взаимодействия, который выбирается из соображений равенства силы магнитного и стерического взаимодействия при r_ij = d, т. е. когда частицы находятся в состоянии физического контакта, и их магнитные моменты направлены вдоль соединяющего их радиус-вектора. Вычисляя из (2) и (3) абсолютные значения магнитной и стерической силы и приравнивая их при r_ij = d, легко определяем $\epsilon ={\mathrm{\mu }}_0\frac{m^2}{16\pi d^3}$. Отметим, что расстояние r_cut соответствует минимуму потенциала Леннард-Джонса.

Для моделирования несущей жидкости используется метод термостата Ланжевена. Введем декартову систему координат x, y, z с осями, направленными вдоль ребер симуляционной ячейки; однородное внешнее поле Н направлено вдоль оси x (см. рис. 1). Для частиц, состоящих из магнетита с диаметром примерно равным 16 нм, параметр Ланжевена $\mathrm{\kappa }={\mathrm{\mu }}_0\frac{mH}{kT}$, равный 0.0336 и 0.336, соответствует напряженностям магнитного поля 0.1 и 1 кА/м. Отметим, что параметр κ равен отношению энергии взаимодействия частицы с полем к тепловой энергии системы kT.

 

Рис. 1. Скриншоты гетерогенных структур, образованных частицами во внешнем поле. Характеристики системы: параметр диполь-дипольного взаимодействия между магнитными частиц λ = 7; их объемная концентрация φ = 0.5% (a, б), φ = 6% (в, г); магнитный момент частицы m = 1.1∙10^–18 А∙м²; температура T = 298 К. (а, в) – напряженность магнитного поля H = 0.1 кА/м (параметр Ланжевена κ = 0.0336); (б, г) – H = 1 кА/м (κ = 0.336).

 

Для учета взаимодействия частиц будем использовать метод суммирования Эвальда [19] для дипольных систем. Этот метод широко используется в компьютерных моделях физико-химических процессов в материалах с дипольно-взаимодействующими частицами [20, 21].

Основная проблема моделирования систем с медленно спадающим потенциалом межчастичного взаимодействия, как это имеет место в дипольных системах, состоит в том, что вычисления медленно сходятся и могут приводить к большим неточностям. Суммирование Эвальда позволяет преодолеть эти сложности, разделяя диполь-дипольный потенциал на две части. А именно, короткодействующую часть, которая может быть эффективно рассчитана в реальном пространстве, и длиннодействующую часть, которая рассчитывается, используя метод Фурье-преобразований по трёхмерному пространству. При этом должна быть выбрана дистанция между частицами, отделяющая близко- и дальнодействующие части потенциала взаимодействия. Выбор этой дистанции определяется выбираемым балансом между точностью и временем расчета. Суммирование Эвальда, включённое в пакет ESPResSo, позволяет провести вычисления на многоядерных процессорах с использованием библиотеки интерфейса передачи сообщений (MPI) [22]. Отметим, что время компьютерных расчетов пропорционально Nlog N. Поскольку при больших N выполняется неравенство logN ≫ 1, то небольшое увеличение числа частиц N приводит к весьма значительному увеличению требуемого времени.

Второй закон Ньютона для трансляционного и вращательного движения частиц может быть записан в следующем виде [21]:

$\mathit{M}{\dot{\mathbf{v}}}_{\mathrm{i}}={\mathbf{F}}_{\mathrm{i}}^{\left(DD\right)}+{\mathbf{F}}_{\mathrm{i}}^{\left(WCA\right)}-\mathrm{\gamma }{\mathbf{v}}_{\mathrm{i}}+{\mathrm{\xi }}_{\mathrm{i}}^{\left(T\right)},$

$\mathit{I}{\dot{\mathbf{\omega }}}_{\mathrm{i}}={\mathbf{T}}_{\mathrm{i}}^{\left(DD\right)}+{\mathbf{T}}_{\mathrm{i}}^{\left(H\right)}-{\mathrm{\gamma }}_{\mathrm{R}}{\mathbf{\omega }}_{\mathrm{i}}+{\mathrm{\xi }}_{\mathrm{i}}^{\left(\mathit{R}\right)},$

$\gamma =3\mathrm{\pi \eta }d, {\gamma }_R=\mathrm{\pi \eta }{d}^3, I=\frac{M{d}^2}{10},$

${\mathit{F}}_{\text{i}}^{\left(\text{DD}\right)}=-\frac{3{\mathrm{\mu }}_0{\mathrm{m}}^2}{4\mathrm{\pi }}{{\sum }_{\text{j}\ne \text{i}}\mathrm{r}}_{\text{ij}}\frac{3\left(e_{\text{i}}\cdot {\mathrm{r}}_{\text{ij}}\right)\left(e_{\text{j}}\cdot {\mathrm{r}}_{\text{ij}}\right)-r_{\text{ij}}^2\left(e_{\text{i}}\cdot e_{\text{j}}\right)}{r_{\text{ij}}^7},$

${\text{F}}_{\text{i}}^{\left(\text{WCA}\right)}=\left\{\begin{array}{c}\mathrm{0,} \text{если} r_{\text{ij}}>r_{\text{cut}}\\ 24\epsilon \sum _{\text{j}\ne }\left(\frac{2d^{12}}{r_{\text{ij}}^{14}}-\frac{d^6}{r_{\text{ij}}^8}\right){\mathbf{r}}_{\text{ij}}, \text{если} r_{\text{ij}}<r_{\text{cut}}\end{array}\right.,$ (4)

$\begin{array}{l}{\text{T}}_{\text{i}}^{\left(\text{H}\right)}={\mathrm{\mu }}_{0}m\left[e_{\text{i}}\times \mathrm{H}\right], \\ {\text{T}}_{\text{i}}^{\left(\text{DD}\right)}=\frac{{\mathrm{\mu }}_0{m}^2}{4\mathrm{\pi }}{\sum }_{\text{j}\ne \text{i}}\left[e_{\text{i}}\times \frac{3\left(e_{\text{j}}\cdot {\text{r}}_{\text{ij}}\right){\text{r}}_{\text{ij}}-r_{\text{ij}}^2{e}_{\text{j}}}{r_{\text{ij}}^5}\right].\end{array}$

Здесь M – масса частицы; I – её момент инерции; v_i и ω_i – трансляционная и угловая скорости i-ой частицы; ${\mathrm{F}}_{\text{i}}^{\left(\text{DD}\right)}$ – сила диполь-дипольного взаимодействия i-й частицы со всеми остальными; ${\mathrm{F}}_{\text{i}}^{\left(\text{WCA}\right)}$ – сила стерического взаимодействия i-й частицы со всеми частицами, находящимися в пределах радиуса r_cut; ${\mathrm{T}}_{\text{i}}^{\left(\text{DD}\right)}$ – момент диполь-дипольных сил, действующих на i-ю частицу; ${\mathrm{T}}_{\text{i}}^{\left(\text{H}\right)}$ – момент сил, действующий на эту частицу со стороны поля Н; γ и γ_R – трансляционный и вращательный коэффициенты Стокса вязкого сопротивления среды движению частицы; η – вязкость несущей жидкости.

Первые статистические моменты случайной силы ${\mathrm{\xi }}_{\text{i}}^{\left(\text{T}\right)}$ и момента сил ${\mathrm{\xi }}_{\text{i}}^{\left(\text{R}\right)}$ равны нулю:

${\xi }_{\text{i}\alpha }^{\left(T\right)}\left(t\right)=\mathrm{0,} {\xi }_{\text{i}\alpha }^{\left(R\right)}\left(t\right)=\mathrm{0,}$ (5)

а их вторые моменты удовлетворяют стандартным соотношениям:

$\begin{array}{l}{\xi }_{\text{i}\alpha }^{\left(\text{T}\right)}\left(t\right)\cdot {\xi }_{\text{j}\beta }^{\left(\text{T}\right)}\left(t^{\text{'}}\right)=6kT\gamma {\delta }_{\text{ij}}{\delta }_{\alpha \beta }\delta \left(t-t^{\text{'}}\right),\\ {\xi }_{\text{i}\alpha }^{\left(\text{R}\right)}\left(t\right)\cdot {\xi }_{\text{j}\beta }^{\left(\text{R}\right)}\left(t^{\text{'}}\right)=6kT{\gamma }_{\text{R}}{\delta }_{\text{ij}}{\delta }_{\alpha \beta }\delta \left(t-t^{\text{'}}\right)\end{array}$(6)

Здесь α и β означают координаты x, y, z; δ_ij и δ_αβ – символы Кронекера; δ(t – t') – дельта функция Дирака.

Для того, чтобы численно решить систему уравнений (4), необходимо задать начальные условия для координат, ориентаций, а также скоростей всех частиц. Разумно выбрать эти условия, предполагая макроскопически однородное (т. е. без агрегатов) расположение частиц и изотропное распределение ориентаций их моментов. Однако частицы, случайно расположенные в симуляционной ячейке, могут перекрывать друг друга. Стерические силы между перекрывающимися частицами очень велики, что может вести к сбою моделирования. Для того чтобы избежать этого, мы используем алгоритм наибольшего градиентного спуска и ограничиваем максимальное смещение частицы за шаг интегрирования до 1% от диаметра частицы. Фрагмент кода, реализующий «разогрев» системы, приведен в Приложении. Начальные условия для скоростей и угловых скоростей частиц можно задавать любыми значениями, в частности нулевыми. В процессе численного решения обыкновенных дифференциальных уравнений (4) система всё равно придет к термодинамическому равновесию за счет случайных сил ${\mathrm{\xi }}_{\text{i}}^{\left(\text{T}\right)}$ и моментов ${\mathrm{\xi }}_{\text{i}}^{\left(\text{R}\right)}$.

Система дифференциальных уравнений (4) была решена методом Верле, встроенного в пакет ESPResSo [23, 24]. Метод Верле – это численный метод, используемый для интегрирования уравнений движения Ньютона. Он широко используется в молекулярно-динамическом моделировании для расчета траекторий частиц. Этот метод ценится за свою простоту, точность и стабильность при долгосрочном моделировании, особенно при работе с консервативными силами.

Для дальнейшего удобно представить физические величины в безразмерной форме, выбирая в качестве единицы энергии абсолютную температуру kT, единицы длины – диаметр частицы d и единицы массы – массу частицы M. Остальные физические величины, фигурирующие в задаче, могут быть выражены в безразмерной форме, используя эти базовые единицы. Например, характерной единицей времени может быть выбрана величина:

$t_{\text{d}}=d\sqrt{\frac{M}{kT}}.$ (7)

Безразмерный коэффициент Стокса γ^* может быть определен в соответствии с соотношением:

${\gamma }^{\text{*}}=\frac{3\mathrm{\pi }d{\mathrm{\eta }}^2}{\mathrm{\rho }kT}.$ (8)

где ρ – плотность материала частицы.

Для магнетитовых частиц диаметром около 16 нм в среде с вязкостью η ~ 0.1 Па∙с при комнатной температуре по порядку величины γ^* ~ 10⁸. При таких значениях этого параметра время расчета оказывается нереалистически большим. Для уменьшения времени расчета приходится выбирать существенно меньшие значения γ^*. С математической точки зрения это равносильно увеличению плотности материала в уравнениях (4) и, как следствие, увеличению времени инерционной релаксации трансляционной и угловой скоростей частицы. Чтобы не выходить за рамки физической адекватности модели, необходимо, чтобы при выбираемых значениях γ^* времена инерционной релаксации частицы оставались существенно меньше всех остальных характерных времён задачи. Для расчетов было выбрано значение γ^* = 12, так как при таком выборе этого параметра время расчетов оказывалось вполне разумным, а время инерционной релаксации трансляционной и угловой скоростей частиц оставалось существенно меньше остальных характерных времен задачи (см. обсуждение ниже).

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Нашей целью является исследование зависимости намагниченности M_F от времени t после ступенчатого изменения напряженности магнитного поля. Число частиц в симуляционной ячейке варьировалось от 300 до 30 000, в зависимости от сложности вычислений.

При расчетах мы варьировали безразмерную объемную концентрацию частиц φ, безразмерный параметр диполь-дипольного взаимодействия частиц $\lambda =\frac{{\mathrm{\mu }}_0}{4\pi }\frac{m^2}{d^{3}kT}$ и определенный выше параметр Ланжевена κ. Вычисления проводились в две стадии. Первая – подготовительная. На этой стадии мы вводим начальное напряжение магнитного поля, чтобы привести систему к термодинамическому равновесию и интегрируем уравнения движения (4) вплоть до достижения системой равновесного состояния, при котором намагниченность перестаёт изменяться.

Если параметр λ существенно больше единицы, магнитные частицы объединяются в гетерогенные структуры. После включения магнитного поля эти структуры ориентируются в направлении поля. Скриншоты типичных структур, приведенные на рис. 1, были получены используя метод визуализированной молекулярной динамики [25].

Несмотря на низкую объемную концентрацию частиц на рис. 1, наши результаты демонстрируют объединение частиц в сильно изогнутые цепочки и замкнутые кольца при параметре Ланжевена κ = 0.0336, что для магнетитовых частиц с диаметром 16 нм соответствует напряжённости магнитного поля 0.1 кА/м. Заметим, что подобные структуры наблюдались в лабораторных экспериментах [26] и компьютерном моделировании [27] в системах магнитных наночастиц в отсутствие поля. При увеличении поля цепочки заметно вытягиваются вдоль поля, кольца практически исчезают.

На второй стадии расчетов мы изменяем напряженность магнитного поля и снова интегрируем уравнения (4), вычисляя безразмерный средний момент частицы $M_{\text{F}}^{\text{*}}$ в каждый момент времени t используя формулу:

$M_{\text{F}}^{\text{*}}=\frac{1}{N}\underset{\text{i}=1}{\overset{\text{N}}{}}{e}_{\text{ix}}.$ (9)

Здесь e_ix – проекция единичного вектора e_i вдоль оси x, т. е. в направлении магнитного поля. Очевидно, при нулевом поле $M_{\text{F}}^{\text{*}}=0$; при бесконечно сильном поле $M_{\text{F}}^{\text{*}}=1$.

Размерная намагниченность M_F системы частиц может быть вычислена по формуле:

$M_{\text{F}}=\mathrm{\phi }{M}_{\text{s}}{M}_{\text{F}}^{\text{*}}.$ (10)

Здесь M_s – намагниченность насыщения материала частицы.

Расчётную формулу для безразмерного времени t^* можно получить на основе первого уравнения в (4) и формулы (7):

$t^{\text{*}}=\frac{t}{t_{\text{d}}}=\frac{kT}{3\mathrm{\pi \eta }{d}^3}t.$ (11)

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

$\frac{d{M}_{\text{F}}}{dt}=\frac{{\mathrm{\chi }}_{0}H-M_F}{\mathrm{\tau }}, {\mathrm{\chi }}_0=8\mathrm{\phi }\lambda .$ (12)

Здесь χ₀ – начальная магнитная восприимчивость феррожидкости; τ – характерное время вращательной диффузии частицы в нулевом поле, равное (см, например, [3, 4, 28, 29]):

$\tau =\frac{3\text{v}\mathrm{\eta }}{kT}, \text{v}=\frac{\mathrm{\pi }{d}^3}{6}.$ (13)

Обсудим теперь физическую адекватность использования оценки γ^* = 12 при численных расчетах. Из второго уравнения (4) получаем оценку характерного времени вращательной инерции частицы τ_I = I/γ_R. Легко можно показать, что выполняется отношение τ/τ_I = 10γ^*. Следовательно, при γ^* = 12 время вращательной диффузии частицы на два порядка величины больше времени ее инерционной релаксации, и численные расчеты, выполненные при γ^* = 12, соответствуют безынерционному приближению, оправданному применительно к наночастицам в реальных жидкостях. Поэтому выбранное значение γ^* не нарушает физической адекватности расчетов.

Уравнение (12) имеет решение:

$M_F={\chi }_0\left[H+\left(H_0-H\right)e^{-\frac{t}{\tau }}\right].$ (14)

Здесь H₀ – начальное значение напряженности магнитного поля. Необходимо подчеркнуть, что соотношения (12)–(14) справедливы только при слабых полях, когда выполняется линейное по полю приближение для намагниченности.

Для набора статистики при λ = 5 мы проводили вычисления пять, а при λ = 7–10 раз. Во всех симуляциях мы использовали безразмерный шаг по времени Δt^* = 0.005; безразмерный коэффициент трения γ^* = 12, безразмерную абсолютную температуру (kT)^* = 1. Среднее значение намагниченности M_F и зависимости намагниченности от концентрации частиц φ, а также от времени t показаны на рис. 2–4. Параметры Ланжевена были подобраны для магнетитовых частиц диаметром 16 нм таким образом, чтобы напряжённость магнитного поля H возрастала скачком от 0.1 кА/м до 1 кА/м для рис. 2, 3 и убывала от 1 кА/м до 0.1 кА/м для рис. 4.

 

Рис. 2. Зависимости средней намагниченности одной частицы (a) и натурального логарифма относительной намагниченности (б) системы частиц от времени t. Вертикальные отрезки показывают стандартное отклонение от среднего. Характеристики системы: вязкость несущей среды η = 0.13 Пa·сек; намагниченность насыщения материала частиц M_s = 5·10⁵ A/м (магнетит); объемная концентрация частиц φ = 0.5%; температура T = 298 K; диаметр частиц d = 16 нм; параметр диполь-дипольного взаимодействия λ = 7; параметр Ланжевена κ возрастает от 0.0336 до 0.336. Сплошная линия на рис. (б) – линейная регрессия (16).

 

Рис. 3. То же, что на рис. 2 при λ = 5. Диаметр частиц d = 14.4 нм; параметр Ланжевена κ возрастает от 0.024 до 0.24.

 

Рис. 4. То же, что на рис. 3 при уменьшении безразмерного магнитного поля κ от 0.24 до 0.024, что для выбранных параметров системы соответствует уменьшению напряжённости магнитного поля H от 1 кА/м до 0.1 кА/м.

 

Соотношение

$\ln \frac{M_{\infty }-M_{\text{F}}}{M_{\infty }-M_0}=at$ (15)

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

${\mathrm{\tau }}_{\text{R}}=-\frac{1}{a}.$ (16)

Поскольку на рис. 2–4 соотношение $\ln \frac{M_{\infty }-M_{\text{F}}}{M_{\infty }-M_0}$ практически линейно зависит от времени, зависимость намагниченности M_F от времени может быть выражена в виде стандартной экспоненциальной функции с одним временем релаксации τ_R. Некоторые результаты расчетов τ_R для систем с той же вязкостью несущей среды и намагниченности материала частиц, что и на рис. 2–4, представлены в табл. 1. Эти результаты показывают, что в рассмотренном диапазоне изменения значений λ, φ и H время релаксации τ_R варьируется в диапазоне от 2.9∙10^–3 до 1.7 секунд; это значительно больше, чем время релаксации в системе невзаимодействующих (λ = 0) частиц и больше времени переориентации цепочечных агрегатов, оцененное в [30]. Результаты, приведенные в таблице, показывают, что увеличение параметра λ диполь-дипольного взаимодействия частиц ведет к значительному увеличению времени релаксации τ_R как после увеличения, так и уменьшения поля H. Этот эффект объясняется относительно медленной, в сравнении с переориентацией одиночных частиц, эволюцией образуемых структур. Чем больше λ, тем длиннее агрегаты, тем больше времени требуется для их эволюции.

 

Таблица 1. Время релаксации системы частиц при различных параметрах системы

Объемная концентрация φ, % магнитных частиц	Параметр диполь-дипольного взаимодействия λ	Время релаксации τR, с, при ступенчатом увеличении напряжённости магнитного поля от 0.1 кА/м до 1 кА/м	Время релаксации τR, с, при ступенчатом уменьшении напряжённости магнитного поля от 1 кА/м до 0.1 кА/м
0.5	7	0.42	1.7
1	7	0.43	1.4
3	7	0.37	0.74
6	7	0.17	0.54
0.5	5	2.5∙10–3	2.5∙10–3
1	5	4.0∙10–3	3.7∙10–3
3	5	4.5∙10–3	6.7∙10–3
6	5	3.4∙10–3	3.6∙10–3
	0	2.1∙10–4	2.1∙10–4

 

Отметим, что, как показывают простые оценки, для типичных магнетитовых феррожидкостей без стерических слоев на частицах значения параметра λ = 5–7 соответствуют диаметру частиц в интервале 14–16 нм.

В приближении невзаимодействующих частиц (λ = 0) и слабого магнитного поля результаты численных расчетов сравниваются с теоретической оценкой (14) на рис. 5. Как видно из этого рисунка, результаты теории и компьютерного моделирования хорошо соответствуют друг другу. Это, в частности, демонстрирует, что выбор γ^* = 12, даже для невзаимодействующих частиц с минимальным временем релаксации намагниченности, соответствует безынерционному приближению и не нарушает физической адекватности результатов. Тем более эта адекватность не нарушается в случае взаимодействующих частиц с намного большим временем релаксации намагниченности. Отметим, что время релаксации также зависит от начального и конечного значения магнитного поля. После уменьшения поля оно оказывается большим, чем при увеличении.

 

Рис. 5. Намагниченность феррожидкости в зависимости от времени при λ = 0. Диаметр магнетитовых частиц d = 16 нм. Остальные параметры такие же, как на рис. 2. Сплошные линии – теоретические результаты, полученные с использованием (14). (а) – параметр Ланжевена κ возрастает от 0.0336 до 0.336 (магнитное поле H скачкообразно возрастает от 0.1 кА/м до 1 кА/м); (б) – κ убывает от 0.336 до 0.0336 (H ступенчато уменьшается от 1 кА/м до 0.1 кА/м).

 

Неожиданной оказалась обнаруженная немонотонная, с максимумом, зависимость τ_R от объемной концентрации частиц φ, проиллюстрированная также на рис. 6. При увеличении параметра λ точка максимума сдвигается влево, т. е. в сторону меньших φ.

Отметим, что увеличение времени релаксации с концентрацией частиц было предсказано теоретически в [13, 14] для однородных (без агрегатов) феррожидкостей.

 

Рис. 6. Зависимость времени релаксации τ_R от объемной концентрации магнетитовых частиц. Магнитное поле H увеличивается с 0.1 кА/м до 1 кА/м; температура T = 298 K. (a) – параметр диполь-дипольного взаимодействия λ = 5; параметр Ланжевена κ возрастает от 0.024 до 0.24; диаметр частиц d = 14.4 нм; (б) – λ = 6; κ возрастает от 0.0288 до 0.288; d = 15.3 нм.

 

Физически немонотонная зависимость τ_R от φ может быть объяснена следующим образом. Когда концентрация φ очень мала, эффекты межчастичного взаимодействия пренебрежимо малы и время релаксации τ_R системы частиц практически равно времени τ вращательной диффузии одиночной частицы. Увеличение концентрации вначале ведет к усилению межчастичных корреляций, исследованных в [13, 14] и затем к появлению агрегатов частиц. Это увеличивает время релаксации. Однако дальнейший рост φ увеличивает темп агрегирования; это уменьшает время τ_R.

Ломаный характер кривых на рис. 2–6 объясняется конечным числом N частиц в системе и относительно небольшим (до десяти раз) повтором вычислений для набора статистики результатов. Однако даже эти результаты позволяют сделать принципиальные выводы о влиянии структурирования частиц на характерное время перемагничивания феррожидкости. Увеличение числа частиц N быстро приводит к нереалистически большим временам расчета.

Отметим, что характерное время эволюции структур в феррожидкости порядка 1–10 секунд ранее получено в нашей работе [31] в рамках моделей цепочек как прямых агрегатов. Представленные результаты компьютерных расчетов показывают, что внутренние структуры в феррожидкостях могут быть более сложными, чем это предполагается во многих теоретических моделях.

ЗАКЛЮЧЕНИЕ

В работе представлены результаты компьютерного моделирования эволюции внутренних структур в феррожидкостях и влияния этого процесса на кинетику изменения их намагниченности после изменения магнитного поля. Результаты показывают, что время релаксации намагниченности феррожидкости со структурирующимися частицами может быть на несколько порядков величины больше времени релаксации этих систем с одиночными, даже взаимодействующими частицами. Время релаксации намагниченности сильно зависит от параметра λ магнитного взаимодействия частиц и от их объемной концентрации φ (табл. 1).

Отметим, что в работе, в области выбранных значений концентраций φ и параметра λ, наблюдалось образование линейно-цепочечных и ветвистых структур. Известно, что при увеличении φ и λ возможно расслоение частиц на плотную и разбавленную фазы [4, 5]. Поскольку увеличение магнитного поля способно спровоцировать фазовое расслоение феррожидкости, а его уменьшение – привести к распаду плотной фазы [4, 5], то, очевидно, эволюция фазового состояния системы частиц способна существенно повлиять на кинетику перемагничивания феррожидкости. Анализ этого эффекта требует отдельного исследования.

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа выполнена в рамках финансовой программы Министерства науки и образования РФ, проект FEUZ-2020-0051.

СОБЛЮДЕНИЕ ЭТИЧЕСКИХ СТАНДАРТОВ

В данной работе отсутствуют исследования человека или животных.

КОНФЛИКТ ИНТЕРЕСОВ

Авторы заявляют, что у них нет конфликта интересов.

ПРИЛОЖЕНИЕ

В этом приложении комментарии предваряются символом решетки #

# Расчет длины куба для моделирования

BOX_SIZE = LJ_SIGMA * (PI * N / (6.0 * PHI)) ** (1.0 / 3.0)

# Подключение пакета ESPResSo

SYSTEM = espressomd.System(box_l = (BOX_SIZE, BOX_SIZE, BOX_SIZE))

SYSTEM.time_step = TIME_STEP

# Мы представляем взаимодействие между частицами как несвязное взаимодействие и в качестве потенциала взаимодействия используем потенциал Леннарда-Джонса. Здесь мы используем упомянутый выше r_cut, чтобы получить отталкивающее взаимодействие.

SYSTEM.non_bonded_inter[0, 0].lennard_jones.set_params(epsilon = LJ_EPSILON, sigma = LJ_SIGMA, cutoff = LJ_CUT, shift = “auto”)

# Задание стохастического равномерного распределения направлений моментов частиц в пространстве

np.random.seed(seed = SEED)

DIP_PHI = 2.0 * PI * np.random.random((N, 1))

DIP_COS_THETA = 2.0 * np.random.random((N, 1)) - 1.0

DIP_SIN_THETA = np.sin(np.arccos(DIP_COS_THETA))

DIP = np.hstack((DIP_SIN_THETA * np.sin(DIP_PHI), DIP_SIN_THETA * np.cos(DIP_PHI), DIP_COS_THETA))

# Задание параметров для установления термодинамического равновесия

POS = BOX_SIZE * np.random.random((N, 3))

PARTICLES = SYSTEM.part.add(pos = POS, rotation = N * [(1, 1, 1)], dip = DIP)

SYSTEM.integrator.set_steepest_descent(f_max= = 0.0, gamma = 0.1, max_displacement = 0.05)

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

while SYSTEM.analysis.energy()[“total”] > 5 * KT * N:

SYSTEM.integrator.run(20)

About the authors

D. N. Chirikov

Уральский федеральный университет имени Первого Президента России Б. Н. Ельцина

Author for correspondence.
Email: d.n.chirikov@urfu.ru
Russian Federation, ул. Мира, 19, Екатеринбург, 620002

A. Yu. Zubarev

Уральский федеральный университет имени Первого Президента России Б. Н. Ельцина

Email: d.n.chirikov@urfu.ru
Russian Federation, ул. Мира, 19, Екатеринбург, 620002

References

Supplementary files