Коллоидный журнал, 2021, T. 83, № 4, стр. 412-422

ВВЕДЕНИЕ

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

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

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

Впервые в [1] на основе машинного эксперимента методом молекулярной динамики было показано, что автокорреляционная функция скоростей в течение длительного времени, не распадаясь экспоненциально, изменяется по закону ${{t}^{{{{ - d} \mathord{\left/ {\vphantom {{ - d} 2}} \right. \kern-0em} 2}}}},$ где $d$ – размерность системы. Эта работа привлекала внимание многих специалистов, и несколько позже такой же характер затухания автокорреляционных функций был обнаружен в теоретических работах [2–6]. Особое внимание этому вопросу уделили Резибуя с коллегами [7, 8]. Они, как на основе феноменологической теории с применением метода Ландау–Плачека, так и на основе кинетической теории, подробно исследовав поведение временны́х корреляционных функций, показали, что медленное затухание подынтегральных функций Грина–Кубо по степенному закону обусловлено коллективным гидродинамическим эффектом и распространением флуктуации гидродинамических переменных на большие расстояния.

Фишер [6], исследуя происхождение долговременных хвостов корреляционных функций при описании процессов переноса, исходил из того, что скорость теплового движения молекул можно представить в виде суммы двух составляющих:

где ${{\upsilon }_{{\text{c}}}}(t)$ представляет коллективное перемещение ближайших молекул, а $\upsilon {\kern 1pt} {\text{'}}(t)$ – скорость движения выбранной молекулы относительно ее окружения. Принимая эти составляющие как статистически независимые, движение выбранной молекулы он рассматривал как совместное со своим окружением в составе “жидкой частицы”, впоследствии называемой лагранжевой частицей. Далее он предполагал, что временнáя эволюция ${{\upsilon }_{{\text{c}}}}(t)$ – более медленная и автокорреляционную функцию скорости аппроксимировал для частицы Лагранжа. Для расчета автокорреляционной функции, принимая лагранжевую теорию тепловых гидродинамических флуктуаций, Фишер получил результаты, которые полностью согласовывались как с результатами феноменологических подходов [2–4], так и с расчетами, проведенными на основе молекулярно-кинетической теории [5, 7, 8]. Дальнейшее развитие метода, предложенного в [6], можно видеть в работах [9–11], где построена последовательная теория тепловых гидродинамических флуктуаций Лагранжа и с применением метода молекулярной динамики детально анализируется автокорреляционная функция скорости частицы Лагранжа.

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

В обзоре [12], который посвящен методу молекулярной динамики и его применению в статистической физике, подробно проанализировано долговременное поведение временны́х корреляционных функций скорости, автокорреляционных функций сила–сила, плотность–плотность и показано, что их распад всегда подчиняется закону ${{t}^{{{{ - d} \mathord{\left/ {\vphantom {{ - d} 2}} \right. \kern-0em} 2}}}}.$ Там же отмечается, что справедливость результатов, полученных в [1] для трехмерной системы твердых шариков, содержащей 500 частиц, при установлении асимптотического закона ${{t}^{{{{ - 3} \mathord{\left/ {\vphantom {{ - 3} 2}} \right. \kern-0em} 2}}}}$ ставится под сомнение из-за слишком малого числа частиц. Однако исследования, проведенные для системы, содержащей 4000 частиц [13], подтвердили справедливость асимптотического закона ${{t}^{{{{ - 3} \mathord{\left/ {\vphantom {{ - 3} 2}} \right. \kern-0em} 2}}}},$ обнаруженного в [1]. Работы по исследованию свойств простых жидкостей для систем из 32 000 и 10 976 частиц [14], для моделей воды из 10 000 и 50 000 молекул [11], а также дисперсных систем с 8000 частицами и различным соотношением масс наночастиц и молекул несущей среды [15, 16] показывают, что закон ${{t}^{{{{ - 3} \mathord{\left/ {\vphantom {{ - 3} 2}} \right. \kern-0em} 2}}}},$ проявляющийся при исследовании коллективных релаксационных процессов, является всеобщим и его можно наблюдать во всех жидкостях.

Нами ранее [17–20] на основе метода молекулярно-кинетической теории были исследованы вязкоупругие и акустические свойства магнитных жидкостей с учетом различных релаксационных процессов. В [17] на основе кинетических уравнений, учитывающих вклады пространственных корреляций плотности и корреляции скоростей, получены микроскопические уравнения обобщенной гидродинамики, учитывающие различные релаксационные процессы. Также получено уравнение Смолуховского, описывающее поведение неравновесной бинарной плотности ${{n}_{2}}({{{\mathbf{q}}}_{1}},{{{\mathbf{q}}}_{2}},t)$ в конфигурационном пространстве и найдено его общее решение. Показано, что фундаментальным решением уравнения Смолуховского является функция Грина $G(r,{{r}_{1}},t - {{t}_{1}}),$ асимптотика которой при больших временах совпадает с асимптотиками автокорреляционных функций, т.е. имеет место закон ${{t}^{{{{ - 3} \mathord{\left/ {\vphantom {{ - 3} 2}} \right. \kern-0em} 2}}}}.$

В этой и других наших работах, как следствие медленного затухания функции Грина, область частотной дисперсии коэффициентов переноса, модулей упругости и акустических параметров всегда являлась широкой и составляла пять декад. В работе [21], где магнитная жидкость рассматривается как многокомпонентная система, асимптотика функций Грина также подчиняется закону ${{t}^{{{{ - 3} \mathord{\left/ {\vphantom {{ - 3} 2}} \right. \kern-0em} 2}}}}.$

Хотя многие теоретические модели показывают существование широкого спектра частотной дисперсии для многих жидкостей и суспензий, в экспериментальных работах область частотной дисперсии коэффициентов переноса и акустических параметров простых жидкостей [22], а также ферроколлоидов [23] не превышает 1–2 декад. Согласно [22], для некоторых жидкостей с очень большой вязкостью, таких как канифоль и полиметилметакрилат, в области затвердевания удалось наблюдать широкий релаксационный спектр, а для некоторых вязких масел область частотной дисперсии составляла до пяти декад.

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

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ И УРАВНЕНИЯ

Сначала приведем некоторые результаты из работы [21]. В ней для двухкомпонентной магнитной жидкости, состоящей из ферритной (индекс $f$) и молекулярной (индекс $s$) подсистем, получены кинетические уравнения для одночастичной, ${{f}_{{1i}}}({{{\mathbf{q}}}_{1}},{{{\mathbf{p}}}_{1}},t),$ и двухчастичной, ${{f}_{{2i}}}({{{\mathbf{q}}}_{1}},{{{\mathbf{q}}}_{2}},{{{\mathbf{p}}}_{1}},{{{\mathbf{p}}}_{2}},t),$ функций распределения, описывающие необратимые процессы. На основе этих уравнений получена система уравнений обобщенной гидродинамики, в которые входят тензоры напряжения [21]

учитывающие молекулярную структуру системы.

Как видно, тензоры напряжения $\sigma _{i}^{{\alpha \beta }}\left( {{{{\mathbf{q}}}_{1}},t} \right)$ состоят из кинетической $K_{i}^{{\alpha \beta }}\left( {{{{\mathbf{q}}}_{1}},t} \right)$ = ${{P}_{{ki}}}\left( {{{{\mathbf{q}}}_{1}},t} \right){{\delta }^{{\alpha \beta }}} + $ $ + \,\,k_{i}^{{\alpha \beta }}\left( {{{{\mathbf{q}}}_{1}},t} \right)$ и потенциальной частей. Потенциальная часть (1) определяется функцией неравновесной бинарной плотности ${{n}_{{2i}}}({{{\mathbf{q}}}_{1}},{\mathbf{r}},t),$ описывающей процесс перестройки структуры жидкости. Впервые для функции ${{n}_{2}}({{{\mathbf{q}}}_{1}},{\mathbf{r}},t)$ Кирквудом [24] было получено уравнение Смолуховского. Решение этого уравнения описывает пространственно-временнóе изменение структуры жидкости и, следовательно, процесс структурной релаксации. Однако общее решение этого уравнения Кирквудом с последователями не было найдено, и поэтому они не рассматривали релаксационные процессы.

На основе кинетических уравнений в [17] было получено уравнение Смолуховского для бинарной плотности частиц магнитных жидкостей, которое для i-ой подсистемы записывается в виде

и его решение для i-ой подсистемы имеет вид [17] где ${{\gamma }_{i}} = \frac{1}{{{{n}_{i}}{{C}_{{{\text{v}}i}}}}}{{\left( {\frac{{\partial {{P}_{i}}}}{{\partial T}}} \right)}_{{{{n}_{i}}}}},$ ${{\omega }_{{0i}}} = {{2kT} \mathord{\left/ {\vphantom {{2kT} {{{\beta }_{i}}\sigma _{i}^{2}}}} \right. \kern-0em} {{{\beta }_{i}}\sigma _{i}^{2}}},$ ${{n}_{i}},$ ${{\sigma }_{i}},$ ${{\beta }_{i}},$ ${{C}_{{{\text{v}}i}}},$ ${{P}_{i}},$ ${{g}_{i}}$ – соответствующие i-ой подсистеме значения числовой плотности, диаметра частицы, коэффициента трения, теплоемкости, давления и радиальной функции распределения, $\omega $ – частота внешнего возмущения, ${{\mu }_{0}}$ – магнитная постоянная, ${\mathbf{m}}$ – вектор магнитного момента магнитных частиц, $H$ – внешнее магнитное поле, $k$ – постоянная Больцмана, $T$ – температура системы, ${\mathbf{u}}$ – вектор смещения. Для молекулярной подсистемы в выражении для ${{R}_{i}}\left( {{{{\mathbf{q}}}_{1}},{\mathbf{r}},t} \right)$ слагаемое, учитывающее влияние внешнего магнитного поля, равно нулю.

Подставляя (3) в (1) и совершая фурье-преобразование по времени, после его сравнения с фурье-образом макроскопического тензора вязких напряжений [25]

получим следующие выражения для динамических коэффициентов сдвиговой и объемной вязкости: где ${{G}_{i}}\left( {r,{{r}_{1}},\omega } \right) = \frac{{{{\tau }_{{0i}}}}}{2}{{\left( {\frac{2}{{\omega {{\tau }_{{0i}}}}}} \right)}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}\left[ {\left( {{\text{sin}}{{\varphi }_{{1i}}} - \cos {{\varphi }_{{1i}}}} \right){{e}^{{ - {{\varphi }_{{1i}}}}}}} \right.$ – – $\left. {\left( {\sin {{\varphi }_{{2i}}} - \cos {{\varphi }_{{2i}}}} \right){{e}^{{ - {{\varphi }_{{2i}}}}}}} \right],$ ${{\varphi }_{{(1,2)i}}}\left( {r,{{r}_{1}},\omega } \right)$    –       $ = {{\left( {\frac{{\omega {{\tau }_{{0i}}}}}{2}} \right)}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}\left( {r \mp {{r}_{1}}} \right),$ ${{\tau }_{{1i}}} = \frac{{{{m}_{i}}}}{{2{{\beta }_{i}}}}$ – время трансляционной релаксации вязкого тензора напряжений, ${{\tau }_{{0i}}} = \frac{{{{\beta }_{i}}\sigma _{i}^{2}}}{{2kT}}$ – феноменологическое время структурной релаксации.

Из (4) следует, что при $(t - {{t}_{1}}) \to \infty $ функция ${{G}_{i}}\left( {r,{{r}_{1}},t - {{t}_{1}}} \right)$ имеет асимптотику

В гидродинамическом режиме, когда $\omega {{\tau }_{0}} \ll 1,$ разлагая в (6) и (7) функцию ${{G}_{i}}\left( {r,{{r}_{1}},\omega } \right)$ в ряд по степеням $\omega {{\tau }_{0}}$ и ограничиваясь членами второго порядка, для низкочастотных асимптотик динамических коэффициентов вязкости магнитных жидкостей получим: где ${{A}_{{si}}}{\kern 1pt} {\kern 1pt} = {\kern 1pt} {\kern 1pt} \frac{{n_{i}^{2}\sigma _{i}^{3}{{\tau }_{{0i}}}}}{{60\pi }}{\kern 1pt} \int_0^\infty {r\frac{{\partial {{\Phi }_{i}}}}{{\partial r}}} \int_0^\infty {\frac{{\partial {{g}_{i}}({{r}_{1}})}}{{\partial {{r}_{1}}}}{\kern 1pt} {{r}_{1}}d{{{\mathbf{r}}}_{1}}d{\mathbf{r}}} ,$ ${{A}_{{{\text{v}}i}}}{\kern 1pt} {\kern 1pt} = {\kern 1pt} {\kern 1pt} \frac{{n_{i}^{2}\sigma _{i}^{3}{{\tau }_{{0i}}}}}{{12\pi }} \times $ × $\int_0^\infty {r\frac{{\partial {{\Phi }_{i}}}}{{\partial r}}} \int_0^\infty {{{\varphi }_{i}}({{r}_{1}})d{{{\mathbf{r}}}_{1}}d{\mathbf{r}}} ,$ ${{A}_{f}} = \frac{{n_{f}^{2}\sigma _{f}^{3}{{\tau }_{{0f}}}}}{{36\pi }}\int_0^\infty {r\frac{{\partial {{\Phi }_{f}}}}{{\partial r}}} $ × $ \times \,\,\int_0^\infty {\frac{{{{\mu }_{0}}}}{{{{\beta }_{f}}}}({\mathbf{m}}\nabla )\frac{{\partial H}}{{\partial \upsilon }}\frac{{\partial {{g}_{f}}({{r}_{1}})}}{{\partial {{r}_{1}}}}{{r}_{1}}d{{{\mathbf{r}}}_{1}}d{\mathbf{r}}} .$

Видно, что в гидродинамическом режиме коэффициенты вязкости имеют асимптотики, пропорциональные ${{(\omega {{\tau }_{0}})}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}},$ что совпадает с низкочастотными асимптотиками [26–28].

При высоких частотах, когда $\omega {{\tau }_{0}} \gg 1,$ т.е. при $\omega \to \infty ,$ из (6) и (7) для динамических коэффициентов вязкости получим:

где ${{B}_{{si}}} = \frac{{2\pi n_{i}^{2}\sigma _{i}^{3}}}{{15}}\int_0^\infty {\frac{d}{{dr}}\left( {\frac{{\partial {{\Phi }_{i}}}}{{\partial r}}{{r}^{4}}} \right)dr} ,$ ${{B}_{{{\text{v}}i}}} = \frac{{2\pi n_{i}^{2}\sigma _{i}^{3}}}{9} \times $ $ \times \,\,\int_0^\infty {\frac{d}{{dr}}\left( {\frac{1}{{{{r}^{2}}}}\frac{{\partial {{\Phi }_{i}}}}{{\partial r}}} \right){{g}_{i}}(r){{r}^{6}}dr} ,$ B_f   =    $\frac{{2\pi n_{f}^{2}\sigma _{f}^{3}}}{9}{\kern 1pt} $  × × $\int_0^\infty {\frac{d}{{dr}}\left( {\frac{{\partial {{\Phi }_{f}}}}{{\partial r}}\frac{{{{\mu }_{0}}}}{{{{\beta }_{f}}}}({\mathbf{m}}\nabla )\frac{{\partial H}}{{\partial \upsilon }}{\kern 1pt} {{r}^{4}}} \right){{g}_{f}}(r)dr} ,$ ${{K}_{{si}}}{\kern 1pt} {\kern 1pt} = {\kern 1pt} {\kern 1pt} {{n}_{i}}{\kern 1pt} {{\left( {\frac{{\partial {{P}_{i}}}}{{\partial {{n}_{i}}}}} \right)}_{T}}{\kern 1pt} {\kern 1pt} + $ $ + \,\,{{\gamma }_{i}}T{{\left( {\frac{{\partial {{P}_{i}}}}{{\partial T}}} \right)}_{{{{n}_{i}}}}}$ – адиабатический модуль упругости.

Согласно (9), при высоких частотах динамические коэффициенты сдвиговой и объемной вязкости стремятся к нулю по закону ${{\omega }^{{ - 1}}}$ и по виду совпадают с высокочастотными модулями упругости [29, 30].

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

Теперь, с целью упрощения структуры уравнения Смолуховского (2), применяем аппроксимацию Бхатнагара–Гросса–Крука, т.е. релаксационный член Смолуховского ${{\omega }_{{0i}}}{{\hat {L}}_{i}}n_{{2i}}^{'}\left( {{{{\mathbf{q}}}_{1}},{\mathbf{r}},t} \right)$ заменяем релаксационным членом типа $\frac{{n_{{2i}}^{0}{\kern 1pt} \left( {\left| {\mathbf{r}} \right|,n,T} \right) - {{n}_{{2i}}}{\kern 1pt} \left( {{{{\mathbf{q}}}_{1}},{\mathbf{r}},t} \right)}}{{{{\tau }_{{0i}}}}}{\kern 1pt} ,$ где ${{n}_{{2i}}}{\kern 1pt} \left( {{{{\mathbf{q}}}_{1}},{\mathbf{r}},t} \right)$ = $n_{{2i}}^{0}{\kern 1pt} \left( {\left| {\mathbf{r}} \right|,n,T} \right){\kern 1pt} {\kern 1pt} + $ + $n_{{2i}}^{'}\left( {{{{\mathbf{q}}}_{1}},{\mathbf{r}},t} \right),$ ${{\tau }_{{0i}}}$ – время релаксации неравновесной бинарной плотности ${{n}_{{2i}}}\left( {{{{\mathbf{q}}}_{1}},{\mathbf{r}},t} \right)$ к ее локально-равновесному значению $n_{{2i}}^{0}\left( {\left| {\mathbf{r}} \right|,n({{{\mathbf{q}}}_{1}},t),T({{{\mathbf{q}}}_{1}},t)} \right)$ = = $n_{{2i}}^{0}(\left| {\mathbf{r}} \right|) + {{\tilde {n}}_{{2i}}}\left( {\left| {\mathbf{r}} \right|,n,T,t)} \right),$ $n_{{2i}}^{0}(\left| {\mathbf{r}} \right|) = n_{i}^{{\text{2}}}{{g}_{i}}(\left| {\mathbf{r}} \right|)$ – равновесная бинарная плотность, ${{\tilde {n}}_{{2i}}}\left( {\left| {\mathbf{r}} \right|,n,T,t)} \right)$ = = ${{\left( {\frac{{\partial n_{{2i}}^{0}(\left| {\mathbf{r}} \right|)}}{{\partial {{n}_{i}}}}} \right)}_{T}}{{n}_{i}}({{{\mathbf{q}}}_{1}},t)$ + ${{\left( {\frac{{\partial n_{{2i}}^{0}(\left| {\mathbf{r}} \right|)}}{{\partial T}}} \right)}_{{{{n}_{i}}}}}T({{{\mathbf{q}}}_{1}},t).$

Подставляя эти выражения в (2) и совершая фурье-преобразование по времени, с учетом уравнения гидродинамики в локально-равновесном приближении для фурье-образа неравновесной функции бинарной плотности получим:

Теперь, когда определен явный вид функции ${{n}_{{2i}}}\left( {{{{\mathbf{q}}}_{1}},{\mathbf{r}},\omega } \right),$ подставляя ее в выражение для тензора напряжения (1), можно исследовать вязкостные свойства магнитных жидкостей. Однако для вычисления (1) необходимо еще знать и явный вид выражения для кинетической части неравновесного давления ${{P}_{{ki}}}\left( {{{{\mathbf{q}}}_{1}},t} \right)$ и кинетической части вязкого тензора напряжения $k_{i}^{{\alpha \beta }}\left( {{{{\mathbf{q}}}_{1}},t} \right).$

Уравнение для $k_{i}^{{\alpha \beta }}\left( {{{{\mathbf{q}}}_{1}},t} \right),$ как показано в [17], получается из кинетического уравнения для одночастичной функции распределения и в линейном приближении имеет следующее решение:

Аналогично, фурье-образ ${{P}_{{ki}}}\left( {{{{\mathbf{q}}}_{1}},t} \right)$ имеет следующий вид: Теперь, совершая фурье-преобразование по времени в (1) и учитывая в нем ${{n}_{{2i}}}\left( {{{{\mathbf{q}}}_{1}},{\mathbf{r}},\omega } \right)$ из (10), а также значения ${{P}_{{ki}}}\left( {{{{\mathbf{q}}}_{1}},\omega } \right)$ и $k_{i}^{{\alpha \beta }}\left( {{{{\mathbf{q}}}_{1}},\omega } \right),$ после несложных операций, сравнивая полученный результат с фурье-образом макроскопического выражения (5), для динамических коэффициентов сдвиговой и объемной вязкости получим: Таким образом, для динамических коэффициентов вязкости магнитной жидкости мы получили два вида аналитических выражений: 1) выражения (6) и (7), когда процесс перестройки структуры жидкости имеет диффузионный характер, 2) выражения (11) и (12), когда процесс восстановления равновесной структуры жидкости происходит по экспоненциальному закону. Если в первом случае, согласно (8) и (9), дисперсия коэффициентов вязкости имеет степенной характер, то во втором случае, как видно из (11) и (12), частотная дисперсия динамических коэффициентов вязкости происходит по закону $\frac{\tau }{{1 + {{{(\omega \tau )}}^{2}}}}$ и ширина области дисперсии получается узкой.

ВЫБОР МОДЕЛИ И ЧИСЛЕННЫЕ РАСЧЕТЫ

Чтобы на основе выражений (6), (7), (11) и (12) провести численные расчеты, необходимо выбрать явный вид входящих в них потенциалов межчастичного взаимодействия ${{\Phi }_{i}}$ и соответствующие им радиальные функции распределения ${{g}_{i}}.$ В соответствии с [21], для молекулярной подсистемы применяем модель Штокмайера

где $\Phi _{s}^{{{\text{L - J}}}}(r)$  =  $4{{\varepsilon }_{s}}\left[ {{{{\left( {\frac{{{{\sigma }_{s}}}}{r}} \right)}}^{{12}}} - {{{\left( {\frac{{{{\sigma }_{s}}}}{r}} \right)}}^{6}}} \right],$ ${{\Phi }^{{pp}}}({\mathbf{r}})$  = = $\frac{{{{p}^{2}}}}{{4\pi \varepsilon {{\varepsilon }_{0}}}}\left[ {\frac{{({{{\mathbf{e}}}_{a}}{{{\mathbf{e}}}_{b}})}}{{r_{{ab}}^{3}}} - \frac{{3({{{\mathbf{e}}}_{a}}{{{\mathbf{r}}}_{{ab}}})({{{\mathbf{e}}}_{b}}{{{\mathbf{r}}}_{{ab}}})}}{{r_{{ab}}^{5}}}} \right].$

Считая феррочастицы сферическими однодоменными частицами, покрытыми слоем поверхностно-активного вещества и находящимися под действием внешнего магнитного поля, их взаимодействия для магнитной подсистемы представим в виде

где ${{\Phi }^{{mm}}}({\mathbf{r}}) = \frac{{{{\mu }_{0}}{{m}^{2}}}}{{4\pi }}\left[ {\frac{{({{{\mathbf{e}}}_{a}}{{{\mathbf{e}}}_{b}})}}{{r_{{ab}}^{3}}} - \frac{{3({{{\mathbf{e}}}_{a}}{{{\mathbf{r}}}_{{ab}}})({{{\mathbf{e}}}_{b}}{{{\mathbf{r}}}_{{ab}}})}}{{r_{{ab}}^{5}}}} \right],$ ${{\Phi }^{H}}({\mathbf{H}})$ = $ = - \,{{\mu }_{0}}{\mathbf{mH}}.$

Электрическое и магнитное диполь-дипольные взаимодействия при переходе к сферической системе координат (рис. 1), представим в виде

где при $d = p$ ${{A}^{p}} = \frac{1}{{4\pi \varepsilon {{\varepsilon }_{0}}}}$ и при $d = m$ ${{A}^{m}} = \frac{{{{\mu }_{0}}}}{{4\pi }}.$ Значение коэффициента $\chi ({{\theta }_{a}},{{\theta }_{b}},{{\varphi }_{a}} - {{\varphi }_{b}})$ = = $2\cos {{\theta }_{a}}\cos {{\theta }_{b}}$ – $\sin {{\theta }_{a}}\sin {{\theta }_{b}}\cos ({{\varphi }_{a}} - {{\varphi }_{b}})$ зависит от ориентации диполей в сферической системе координат.

В магнитной подсистеме учет отталкивающей силы, возникающей при наличии слоя поверхностно-активного вещества, обеспечивается потенциалом Леннард-Джонса $\Phi _{f}^{{{\text{L - J}}}}(r).$ Хотя этот потенциал при описании взаимодействия коллоидных частиц широко используется, при определении значения глубины потенциальной ямы ${{\varepsilon }_{f}}$ возникает серьезная проблема. Например, согласно [31], в магнитной жидкости с магнитными частицами диаметром $10\,\,{\text{нм,}}$ покрытыми слоем поверхностно-активного вещества толщиной 2 нм, в отсутствие магнитного поля глубина минимума энергии равна $0.3kT.$ При наличии внешнего магнитного поля глубина минимума составляет $0.9kT.$

Радиальные функции распределения для каждой подсистемы выбираем в виде

В магнитных жидкостях объемная доля феррочастиц составляет 10–15%. Поэтому молекулярную подсистему принимаем близкой к системе с плотноупакованными твердыми частицами и в качестве ${{y}_{s}}(\rho {\kern 1pt} *)$ выбираем контактную функцию Карнахана–Старлинга: ${{y}_{s}}(\rho {\kern 1pt} *){\kern 1pt} {\kern 1pt} = {\kern 1pt} {\kern 1pt} {{(2 - \rho _{s}^{*})} \mathord{\left/ {\vphantom {{(2 - \rho _{s}^{*})} {2{{{(1 - \rho _{s}^{*})}}^{3}}}}} \right. \kern-0em} {2{{{(1 - \rho _{s}^{*})}}^{3}}}}{\kern 1pt} ,$ где $\rho _{s}^{*} = {{\pi {{\sigma }^{3}}{{n}_{s}}} \mathord{\left/ {\vphantom {{\pi {{\sigma }^{3}}{{n}_{s}}} 6}} \right. \kern-0em} 6}.$

Необходимо отметить, что непосредственная подстановка потенциалов (13) и (14) в выражения для коэффициентов вязкости приводит к неразрешимым математическим трудностям. Поэтому для упрощения задачи после перехода к сферической системе координат, следуя работам [32, 33], для двух ближайших частиц момент первой частицы принимаем вдоль радиус-вектора, соединяющего ее со второй частицей (рис. 1). Соответственно в (15), положив ${{\theta }_{a}} = 0,$ ${{\theta }_{b}} = \theta $ (рис. 1), энергию взаимодействия двух диполей выражаем в виде

Аналогично, направляя напряженность внешнего магнитного поля ${\mathbf{H}}$ параллельно оси $z$ системы, для потенциальной энергии взаимодействия магнитного диполя с внешним полем получим: Таким образом, подставляя выражения (13)–(16) в (6) и (7), приведем их к виду, удобному для проведения численных расчетов: где ${{a}_{s}} = \frac{{{{p}^{2}}}}{{2\pi \varepsilon {{\varepsilon }_{0}}\sigma _{s}^{3}r_{1}^{3}kT}},$ ${{a}_{f}} = \frac{{{{\mu }_{0}}{{m}^{2}}}}{{2\pi \sigma _{f}^{3}r_{1}^{3}kT}} + \frac{{{{\mu }_{0}}mH}}{{kT}},$ $y_{s}^{'}(\rho {\kern 1pt} *)$ = $ = \frac{{\rho _{s}^{*}(5 - 2\rho _{s}^{*})}}{{(1 - \rho _{s}^{*})(2 - \rho _{s}^{*})}},$ ${{\varepsilon }_{f}},$ ${{\varepsilon }_{s}}$ – глубины потенциальных ям соответствующих подсистем, $l$ – характерный размер системы.

Выражения (13)–(16) подставим в (11) и (12) и приведем их к виду

где Раскрывая суммы в выражениях (17), (18) и учитывая, что ${{n}_{f}} = \frac{{{{N}_{f}}}}{V} = \frac{{{{\rho }_{f}}{{N}_{{\text{A}}}}}}{{{{\mu }_{f}}}}\varphi ,$ ${{n}_{s}} = \frac{{{{N}_{s}}}}{V}\frac{{{{\rho }_{s}}{{N}_{{\text{A}}}}}}{{{{\mu }_{s}}}}(1 - \varphi ),$ после внесения новых обозначений получим: Видно, что, если не учитывать квадратичные зависимости по $\varphi $ при малых значениях концентрации, выражения для коэффициентов сдвиговой и объемной вязкостей принимают вид известного выражения Эйнштейна. При таком предельном переходе аналогичный результат получается и для выражений (19), (20).

Далее на основе выражений (17)–(20) проведен численный расчет динамических коэффициентов вязкости магнитных жидкостей с частицами Fe₃O₄ при $m = {{10}^{{ - 20}}}\,\,{{{\text{Дж}}} \mathord{\left/ {\vphantom {{{\text{Дж}}} {{\text{Тл}}}}} \right. \kern-0em} {{\text{Тл}}}}{\text{,}}$ ${{\sigma }_{f}} = 5\,\,{\text{нм,}}$ ${{\rho }_{f}} = 5340\,\,{{{\text{кг}}} \mathord{\left/ {\vphantom {{{\text{кг}}} {{{{\text{м}}}^{{\text{3}}}},}}} \right. \kern-0em} {{{{\text{м}}}^{{\text{3}}}},}}$ ${{\varepsilon }_{f}} = 0.37kT,$ приготовленных на основе керосина (${{\sigma }_{s}} = 0.43\,\,{\text{нм,}}$ ${{\rho }_{s}} = 819\,\,{{{\text{кг}}} \mathord{\left/ {\vphantom {{{\text{кг}}} {{{{\text{м}}}^{{\text{3}}}}}}} \right. \kern-0em} {{{{\text{м}}}^{{\text{3}}}}}},$ ${{\varepsilon }_{s}} = 0.6kT,$ $\varepsilon = 2.1,$ $p = 8.3 \times {{10}^{{ - 30}}}$ Кл м) и воды (${{\sigma }_{s}} = 0.27\,{\text{нм,}}$ ${{\rho }_{s}} = 1000\,\,{{{\text{кг}}} \mathord{\left/ {\vphantom {{{\text{кг}}} {{{{\text{м}}}^{{\text{3}}}}}}} \right. \kern-0em} {{{{\text{м}}}^{{\text{3}}}}}},$ ${{\varepsilon }_{s}} = 1.27kT,$ $\varepsilon = 81,$ $p = 6.2 \times {{10}^{{ - 30}}}$ Кл м), в широком диапазоне изменения гидродинамической частоты и при различных значениях внешнего магнитного поля.

На рис. 2–5 представлены результаты численных расчетов частотной зависимости динамических коэффициентов вязкости для магнитных жидкостей, приготовленных на основе керосина или воды и частиц Fe₃O₄, при $T = 298\,\,{\text{К}}$ и соответствующих значениях физических параметров, приведенных в табл. 1.

На всех рисунках сплошные кривые соответствуют выражениям (19), (20), а пунктирные – выражениям (21), (22). На рис. 2 и 3 приведены результаты расчетов коэффициента сдвиговой вязкости при отсутствии внешнего магнитного поля и когда значение напряженности поля равно 200 кА/м. На рис. 4 и 5 приведены результаты численных расчетов коэффициента объемной вязкости в отсутствие внешнего магнитного поля.

Видно, что во всех случаях область частотной дисперсии коэффициентов вязкости, вычисленных по выражениям (21) и (22), является узкой, как и наблюдается в экспериментах [22, 23], и не превышает двух декад. Напротив, результаты вычислений, проведенных на основе выражений (19) и (20), указывают на широкую область частотной дисперсии коэффициентов вязкости, составляющую пять декад и более. Однако для фиксированной частоты значения коэффициентов вязкости, полученные на основе выражений (19), (20), всегда больше, чем полученные на основе выражений (21) и (22). Видимо, это связано с более полным учетом вклада структурной релаксации, происходящей по диффузионным законам.

Как известно, при исследовании частотной зависимости коэффициентов переноса и акустических параметров наибольшее изменение этих величин всегда соответствует тем значениям частоты, которые близки к значению частоты релаксации ${{\omega }_{{\text{r}}}} = {1 \mathord{\left/ {\vphantom {1 \tau }} \right. \kern-0em} \tau }.$ В выражениях (19)–(22) фигурируют два времени: время трансляционной релаксации частиц ${{\tau }_{{1i}}} = {{{{m}_{i}}} \mathord{\left/ {\vphantom {{{{m}_{i}}} {2{{\beta }_{i}}}}} \right. \kern-0em} {2{{\beta }_{i}}}}$ и феноменологическое время структурной релаксации ${{\tau }_{{0i}}} = {{{{\beta }_{i}}\sigma _{i}^{2}} \mathord{\left/ {\vphantom {{{{\beta }_{i}}\sigma _{i}^{2}} {2kT}}} \right. \kern-0em} {2kT}}.$ Оба времени релаксации определяются через коэффициент трения $\beta ,$ который тесно связан с параметрами состояния и межмолекулярными силами в системе. Отметим, что проблема определения коэффициента $\beta $ на основе статистических методов до сих пор остается нерешенной, но для его нахождения исследователями предложены несколько моделей, рассмотренных в [34]. На основе анализа, проведенного в [34], было установлено, что значения, полученные для коэффициента трения магнитной жидкости на основе выражения

являются удовлетворительными. Последнее выражение с учетом потенциалов Леннард-Джонса и диполь-дипольных взаимодействий принимает следующий вид: В табл. 1 приведены результаты численных расчетов времен релаксации магнитных жидкостей, проведенных на основе (21). Видно, что трансляционная релаксация в обеих подсистемах магнитных жидкостей протекает очень быстро. Время структурной релаксации для обеих магнитных жидкостей в молекулярной подсистеме меньше, чем в магнитной подсистеме.

Согласно [22], если система релаксирует по нескольким механизмам со значительно различающимися временами релаксации, то кривые, описывающие частотные зависимости коэффициентов переноса, имеют ступенчатую форму с явно выраженными областями дисперсии. Поскольку для магнитной жидкости на основе воды ${{\tau }_{{0f}}}$ почти на два порядка больше, чем ${{\tau }_{{0s}}},$ то на кривых рис. 3 и 5, вычисленных на основе (21) и (22), можно наблюдать две отдельные области дисперсии. На рис. 6 в относительно большом масштабе приведены кривые частотной зависимости ${{\eta }_{s}}(\omega )$ для магнитной жидкости на основе воды, которые наглядно свидетельствуют о существовании отдельных областей дисперсии. Для кривых, вычисленных на основе (19) и (20), вследствие медленного изменения функции Грина области дисперсии перекрываются и такое разделение не наблюдается.

Для сравнения на рис. 7 приведен график из работы [35], где представлены результаты теоретических расчетов частотной зависимости реальной и мнимой частей комплексной динамической сдвиговой вязкости и их сравнение с экспериментальными результатами [36] для суспензий частиц кремнезема, диспергированных в циклогексане.

В теории авторов [35] динамическая вязкость определяется функцией памяти ${{\hat {M}}_{\nu }}(\omega )$ = $ = \int_0^\infty {{{e}^{{i\omega t}}}{{M}_{\nu }}(t)dt} ,$ временнóй спектр ${{\tau }_{M}} = \int_0^\infty {t{{M}_{\nu }}(t)dt} $ которой характеризует скорость ее изменения, имеет длинный хвост, пропорциональный ${{t}^{{{{ - 7} \mathord{\left/ {\vphantom {{ - 7} 2}} \right. \kern-0em} 2}}}}.$ В согласии с этим, при низких частотах функция памяти, а соответственно и динамическая вязкость, имеет асимптотику, пропорциональную ${{\omega }^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}.$ При высоких частотах динамическая вязкость изменяется по закону ${{\omega }^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}},$ и при $\omega \to \infty $ имеет место закон ${1 \mathord{\left/ {\vphantom {1 \omega }} \right. \kern-0em} \omega }.$

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

ЗАКЛЮЧЕНИЕ

В работе исследована частотная зависимость динамических коэффициентов сдвиговой и объемной вязкости, полученных двумя путями на основе: 1) вывода уравнения Смолуховского в конфигурационном пространстве, полностью описывающего пространственно-временнóе поведение функции бинарной плотности; 2) аппроксимации Бхатнагара–Гросса–Крука, приводящей к экспоненциальному затуханию релаксационных процессов. Показано, что в первом случае область частотной дисперсии коэффициентов вязкости является широкой и составляет пять декад, как для некоторых вязких масел [22]. Во втором случае область частотной дисперсии получается узкой, составляя две декады, что наблюдается в многочисленных экспериментах для разных жидкостей.