Известия РАН. Механика жидкости и газа, 2023, № 3, стр. 12-21

1. ОСНОВНЫЕ ПРИБЛИЖЕНИЯ И МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

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

Рассматриваемая модельная система проиллюстрирована на рис. 1. Левый конец щели закрыт непроницаемой для жидкости стенкой, моделирующей тромб; в щели невдалеке от стенки находится капля растворимой феррожидкости.

Как и в [5], здесь рассмотрим случай вращающегося магнитного поля, создаваемого четырьмя соленоидами, показанными на рис. 1.

Расстояние 2b между соленоидами вдоль оси щели предполагается много большим толщины щели l и характерного размера капли. Мы будем предполагать, что внутри капли (облака) феррожидкости эффективная вязкость среды η совпадает с вязкостью несущей жидкости вне капли. Такое приближение вполне оправдано, если объемная концентрация Ф частиц в капле не превышает нескольких процентов. Начальная (в нулевой момент времени) ${{{{\Phi }}}_{0}}\left( {x,z} \right)$ концентрация частиц предполагается известной. Для максимального упрощения вычислений мы рассмотрим двумерную модель, предполагая, что все физические величины зависят только от координат x, z и времени. Наконец, предположим, что магнитное поле достаточно велико, чтобы энергия Зеемана взаимодействия наночастицы с полем была существенно выше тепловой энергии системы. Поэтому мы пренебрежем броуновским вращением частиц. Очевидно, что с практической точки зрения случай сильных полей представляет наибольший интерес с точки зрения генерации течений в жидкости.

Мы предположим, что все соленоиды, иллюстрированные на рис. 1, создают переменное поле с одной и той же угловой частотой ω; ток, и, следовательно, поле соленоидов 1 и 4 зависят от времени как , ток и поле поле соленоидов 2 и 3 – как ${\text{sin}}\omega t$. Тогда полное поле в точке x, z можно представить в виде

Здесь ${{H}_{{01}}},~{{H}_{{02,~~}}}{{H}_{{03~~}}}$ и ${{H}_{{04~~}}}$ – амплитуды полей, создаваемых соленоидами с номерами 1–4 соответственно. Знак “минус” в скобках (1.1) означает, что северный полюс соленоида 1 расположен напротив северного полюса соленоида 2; аналогичная ситуация с соленоидами 3 и 4. Амплитуды ${{H}_{{01}}},~{{H}_{{02,~~}}}{{H}_{{03~~}}}$ и ${{H}_{{04~~}}}$ могут быть найдены из стандартных результатов расчетов полей, создаваемых соленоидами. Явные выражения для них даны в Приложении. Результирующее поле $H = \sqrt {{{H}_{x}}^{2} + {{H}_{z}}^{2}} $ для некоторых моментов времени изображено на рис. 2.

В приближении механики сплошных сред уравнение течения жидкости с феррочастицами может быть представлено в виде (см., например, [7, 8])

Здесь ${{\Delta }} = \frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}} + \frac{{{{\partial }^{2}}}}{{\partial {{z}^{2}}}}$ оператор Лапласа, Γ – момент магнитных сил, действующих на единицу объема жидкости, θ – угол между направлением магнитного момента m частицы и осью Oz, показанной на рис. 1; ${{{{\mu }}}_{0}} = 4{{\pi }} \times {{10}^{{--7}}}$ Гн/м – магнитная проницаемость вакуума; М – намагниченность материала частицы; ${\mathbf{F}} = {{\mu }_{0}}{{\Phi }}\left( {{\mathbf{M}} \cdot \nabla } \right){\mathbf{H}}~$ – пондеромоторная сила, действующая на единицу объема жидкости в неоднородном магнитом поле Н; ${{\Phi }}$ – объемная концентрация частиц в капле; ${\mathbf{M}} = {\mathbf{m}}{\text{/}}{{V}_{p}},~M = \left| {\mathbf{M}} \right|$; V_p – объем частицы. Величина 1/2Γ в (1.1) является антисимметричным напряжением, возникающим в феррожидкости, вследствие действия момента сил Γ [7, 8]. Отметим, что уравнение (1.2) записано в форме, соответствующей малым числам Рейнольдса, построенных на скорости течения среды и толщине щели. Обоснование этого приближения обсуждается в конце раздела Результаты.

Граничные условия к уравнению (1.1) имеют вид

В рамках “не броуновского приближения” уравнение для угла θ может быть записано в виде [8]

Расчеты [5, 6] показывают, что для реалистических условий член $\left( {\frac{{\partial {{{v}}_{x}}}}{{\partial z}} - \frac{{\partial {{{v}}_{z}}}}{{\partial x}}} \right){\text{/}}2~$ в (1.4) мало посравнению двумя другими членами. Очевидно, по порядку величины должно выполняться соотношение $\partial \theta {\text{/}}\partial t\sim {{\omega }}$. Следовательно, из (1.4) получаем оценку ${{\Gamma }}\sim 6{{\eta \Phi \omega }}$. Учитывая, что характерный размер капли, внедренной в канал, в поперечном к его оси направлении не может быть больше соответствующего размера канала l, производная $\partial {{\Gamma /}}\partial z$ может быть оценена так $\partial {{\Gamma /}}\partial z{\text{\;}}\sim {{\Gamma /}}l$. В экспериментах обычно ${{\omega }}\sim 10~$ c^–1, для оценки вязкости несущей жидкости примем, что она примерно равна вязкости крови ${{\eta }}\sim 5~{\text{\;}} \times {{10}^{{--3}}}$ Па с. Кровеносные сосуды, подверженные риску тромбирования, как правило, имеют диаметр порядка нескольких миллиметров. Следовательно, по порядку величины $\partial {{\Gamma /}}\partial z{\text{\;}}\sim 6\eta {{\Phi \omega /}}l\sim 300{{\Phi }}$ Па/м. Очевидно, в продольном к оси канала направлении размер капли может быть больше, чем ширина канала, это значит, что должно выполняться неравенство $\partial {{\Gamma /}}\partial x < \partial {{\Gamma /}}\partial z$.

Оценим теперь величину плотности пондеромоторной силы F. В медико-биологических приложениях обычно используются частицы окислов железа – магнетита или маггемита. Их преимущества – они доступны, обладают достаточно высокими магнитными характеристиками и, главное, биосовместимы. Для этих частиц намагниченность насыщения материала $M\sim 5~{\text{\;}} \times {{10}^{5}}$ А м (см., например, [9]). В экспериментах обычно используются поля с напряженностью H ~ ~ $5~{\text{\;}} \times {{10}^{4}}$ А м. Характерный линейный масштаб изменения поля в области капли феррожидкости по порядку величины определяется ее расстояниями a, b до электромагнитов. Имея в виду, что электромагниты должны находиться вне организма пациента, можно принять, что эти расстояния должны быть порядка десяти сантиметров. Отсюда $F\sim 4\pi \cdot \Phi \times {{10}^{{--7}}} \cdot {\text{\;}}5 \times {{10}^{5}}~{\text{\;}}5 \times {{10}^{4}}{\text{/}}{{10}^{{--1}}}$ Па/м × × 10⁵ Па/м. Следовательно, плотность пондеромоторной силы $F$ для типичных ситуаций оказывается намного больше плотности сил $\partial {{\Gamma /}}\partial x{\text{\;}},\partial {{\Gamma /}}\partial z$ антисимметричных напряжений. Поэтому последние можно не учитывать в уравнениях (1.2). Также несложно показать, что для рассматриваемых систем инерционные члены $\rho \partial {{{v}}_{x}}{\text{/}}\partial t$ и $\rho \partial {{{v}}_{z}}{\text{/}}\partial t$ в (1.2) не играют существенной роли и ими, в первом приближении, можно пренебречь. Действительно, по порядку величины характерное время инерционной релаксации жидкости в канале может быть оценено как ${{\tau }_{{in}}}\sim \rho {{l}^{2}}{\text{/}}\eta .$ Для используемых характерных значений ρ, $l$ и ${{\eta }}$ получаем ${{\tau }_{{in}}}\sim 0.2~с.$ При угловой частоте поля $\omega \sim 10$ c^–1, период изменения поля Следовательно, время инерционной релаксации жидкости существенно меньше периода поля, что позволяет в первом приближении рассматривать уравнения (1.2) в квазистационарной форме, пренебрегая частными производными $\rho \partial {{{v}}_{x}}{\text{/}}\partial t$ и $\rho \partial {{{v}}_{z}}{\text{/}}\partial t$.

В рамках этих приближений уравнения течения в канале могут быть представлены в виде:

с граничными условиями (1.3).

Плотность пондеромоторной силы F можно представить в виде:

где $\nabla $ – дифференциальный оператор Гамильтона.

Из уравнения (1.4) с учетом явного вида магнитного момента Г, приведенного в (1.2), следует, что по порядку величины характерное время ${{\tau }_{\theta }}$ релаксации угла θ может быть оценено так: ${{\tau }_{\theta }}\sim \eta {\text{/}}{{\mu }_{0}}MH$. Подставляя сюда приведенные выше оценки для ${{\eta }},~М$ и H, получаем ${{\tau }_{\theta }}\sim {{10}^{{--6~}}}с$. При частотах поля ${{\omega }}\sim 10$ c^–1, время ${{\tau }_{\theta }}$ намного меньше периода поля 2π/ω. Поэтому можно считать, что в каждый момент времени вектор М практически параллелен вектору Н. Это позволяет использовать следующую приближенную форму соотношения (1.6)

Для дальнейшего удобно ввести функцию тока ${{\Psi }}$ так, что

Отметим, что при такой форме записи компонент скорости v, уравнение несжимаемости ${\text{div}}{\mathbf{v}} = 0$ выполняется автоматически.

Вычислив ротор от обеих частей первых двух уравнений (1.2) и учитывая (1.5) и (1.7), а также уравнение ${\text{rot}}{\mathbf{H}} = 0$, после несложных вычислений получаем

Граничные условия (3) сейчас могут быть переписаны так

Последнее условие (1.10) вместе с условием $\partial {{\Psi /}}\partial x = 0$ при z = 0, l означает,что на границах щели функция тока ${{\Psi }}$ не зависит от координаты х, т.е. равна постоянной величине. Поскольку конкретное значение этой величины не имеет физического смысла, можно придать ей нулевое значение. Тогда это условие может быть переписано в виде

Уравнение (1.9) содержит четвертые производные от искомой функции ${{\Psi }}$ по координатам. Точные аналитические решения таких уравнений не известны. Численное их решение тоже сталкивается с плохо преодолимыми сложностями; надежные алгоритмы решения таких задач нам не известны. Поэтому здесь воспользуемся идеями метода сращиваемых асимптотических разложений, суть которого состоит в сращивании решений (1.9–1.11), получаемых в различных участках течения.

Прежде всего, рассмотрим область канала вдали от тромба, на расстояниях значительно превышающих ширину l канала и характерный размер капли (облака) феррожидкости. В этой области линейный масштаб изменения скорости течения вдоль продольной оси х намного больше, чем в поперечном направлении вдоль оси z. Поэтому уравнение (1.9) с граничными условиями (1.11) можно представить в виде

Обозначим решение этой задачи ${{{{\Psi }}}_{0}}$. Оно имеет вид

Обсудим теперь поведение функции ${{\Psi }}$ в непосредственной близости от тромба, т.е. при $x~{\text{\;}} \to 0$. Для чего разложим функцию ${{\Psi }}$ в ряд Тейлора вблизи точки x = 0:

где штрихи в ${{\Psi }}$ означают производные по х при x = 0. Граничное условие (10) $\partial {{\Psi /}}\partial z = 0$ при x = 0 означает, что на “тромбе” функция ${{\Psi }}\left( {x = 0,z} \right)$ не зависит от координаты z. Поскольку при $z = 0,l$ выполняется граничное условие ${{\Psi }} = 0$, это означает, что ${{\Psi }}\left( {х = 0,z} \right) = 0{\text{\;}}.$ Условие $\partial {{\Psi /}}\partial x = 0$ при x = 0 означает, что ${{\Psi '}} = 0.$ Таким образом, первым неисчезающим членом в разложении (1.14) является член . Учитывая это, представим функцию ${{\Psi }}$ в экстраполяционной форме, удовлетворяющей всем граничным условиям: где α – коэффициент, который определяем из условия наилучшего соответствия формы (1.15) дифференциальному уравнению (1.9). Это условие формулируем следующим образом. Запишем функцию и выберем α, обеспечивающее наименьшее значение L. Идеологически такой подход близок к классическому вариационному методу Онсагера определения термодинамичеких функций жидких кристаллов [10, 11].

Перепишем (1.16) в виде

Отсюда получаем уравнение для определения α, которое может быть решено численно:

2. РЕЗУЛЬТАТЫ

Решение задачи (1.12–1.18) определяется концентрацией феррочастиц Ф, фигурирующей в выражении (1.9) для функции G. Вообще говоря, эти уравнения для функции тока ${{\Psi }}$ необходимо решать совместно с уравнением для концентрации Ф, учитывая, что ее изменения со временем могут просходить за счет миграции частиц в неоднородном магнитном поле; за счет их диффузии и вследствие конвективного движения в генерируемых потоках несущей жидкости. Совместное решение уравнений (1.9–1.12) для ${{\Psi }}$ и стандартного уравнения неразрывности для Ф делает задачу существенно нелинейной, численное решение которой требует больших компьютерных ресурсов.

Заметим, однако, что в работе [5] показано, что для рассматриваемых систем характерное время магнитофоретической миграции магнитных наночастиц с диаметром 15–20 нм на расстояния, сопоставимые с шириной канала, составляет порядка 140 ч, что намного больше времен, представляющих интерес с точки зрения транспорта тромболитиков. Диффузионное распространение частиц на эти расстояния также требует порядка ста часов. Поэтому мы ограничимся анализом начальной стадии генерации течений жидкости в канале, предполагая, что концентрация Ф на этой стадии равна своему начальному значению. Интервал времени, когда это приближение оправданно, обсуждается ниже.

Конкретное начальное значение функции $~{{\Phi }}\left( {x,z} \right)$ определяется конкретным способом инжектирования феррожидкости в канал (кровеносный сосуд) и вряд ли может быть задано в общем виде. Здесь в качестве примера рассмотрим ситуацию, когда профиль концентрации феррочастиц не зависит от поперечной координаты z и имеет вид Гауссова распределения вдоль оси х

Здесь дисперсия σ задает характерный размер капли; х₀ – координата ее центра (расстояние от центра до тромба). Физически такое приближение означает, что размер инжектированной капли феррожидкости намного больше ширины l канала; после инжектирования, за счет диффузии, капля несколько расплывается вдоль оси канала.

Характерное время изменения профиля концентрации за счет конвективного переноса частиц по порядку величины равно отношению характерного размера капли σ к среднеквадратичной скорости ${{{\bar {v}}}_{x}}$ течения в направлении оси канала. Последняя, для выбранных параметров системы, оценена в дальнейшем в (2.4) как несколько десятков микрон в секунду. Если ${{\sigma }}\sim 1$ см, то характерное время конвективной деформации профиля Ф оказывается порядка нескольких сотен или тысячи секунд. По крайней мере, на таких временных интервалах приближение независимости Ф от времени является вполне оправданным.

Используя (2.1), получаем

Функции J и I, определенные в (1.13), сейчас имеют вид:

Уравнение (1.18) относительно коэффициента α было решено числено c учетом (2.2) и (2.3). Продольная компонента скорости течения феррожидкости в канале может быть получена из (1.8) и (1.15)

Результаты расчетов ${{{v}}_{x}}$ из (2.4) представлены на рис. 3 и 4. При расчетах были выбраны следующие значения параметров системы: объемная концентрация частиц в центре капли ${{{{\Phi }}}_{0}}$ = = 0.016; дисперсия σ =1 см. Для всех соленоидов характеристики одинаковы: диаметр D = 1 см; ток I = 8 А; высота h = 1 см; количество витков N = 10⁴; координата центра облака ${{x}_{0}} = --1~$ мм, ширина канала $l = 2~$ мм. Намагниченность материала ферромагнитных частиц $M = 500$ кА/м. Угловые частоты поля ω = 10 и 20 ${{c}^{{--1}}}$. Можно заметить, что увеличение частоты магнитного поля ведет к росту скорости течения феррожидкости. На рис. 3б видно, что непосредственно рядом с тромбом ($х\sim --0.012~$ мм) продольная скорость течения ${{{v}}_{x}}$ порядка 0.1 мкм/c. На расстояниях от тромба порядка размера капли феррожидкости ($х\sim 10--12$ мм) (рис. 3а) скорость достигает значений 350 мкм/c. На рис. 4 представлена скорость ${{{v}}_{x}}$ при разных моментах времени. На рис. 5 показана скорость течения ${{{v}}_{x}}$ как функция от координаты z , поперечной к оси канала при фиксированной координате х.

Среднеквадратичную скорость ${{{\bar {v}}}_{x}}$ циркуляционных течений в окрестности феррожидкости можно оценить как несколько десятков микрон в секунду, что на полтора-два порядка величины больше скорости диффузионного транспорта тромболитиков в неподвижной среде, оцененной в [2], примерно 0.8 мкм/с. Это подтверждает эффективность магнитного метода интенсификации транспорта лекарств в тромбированных сосудах. Решение задачи о конвективной диффузии молекулярной примеси (тромболитиков) в канале с магнитоиндуцированными циркуляционными течениями, с учетом возникновения тонкого диффузионного пограничного слоя вблизи тромба, может являться предметом отдельного исследования.

Обсудим теперь правомерность приближения малых чисел Рейнольдса, которому соответствует уравнение (1.2). Как видно из рис. 3–5, скорость течения жидкости не превышает 500 мкм/с. При таких скоростях для канала диаметром порядка 1 мм (что, как отмечалось, соответствует диаметру кровеносного сосуда, тромбирование которого опасно для здоровья и жизни человека) и несущей жидкости с вязкостью крови (примерно, $4--5~{\text{\;}} \times {{10}^{{--3}}}$ Па с) число Рейнольдса порядка 0.1. Это обосновывает выбранное линейное приближение (1.2) течения жидкости.

ЗАКЛЮЧЕНИЕ

Рассмотрена модель генерации циркуляционных течений в канале с немагнитной жидкостью и инжектированной в нее облаком нанодисперсной феррожидкости под действием бегущего осциллирующего магнитного поля. Один конец канала предполагался закрытым (тромбированным). Расчеты показывают, что, при вполне реалистических параметрах системы, частоты и напряженности магнитного поля, в канале вблизи облака феррожидкости могут быть генерированы течения с амплитудой продольной скорости порядка 10–100 мкм/с. Амплитуда скорости максимальна в центре капли и убывает по мере удаления от нее. Эти потоки могут существенно интенсифицировать транспорт нейтральной молекулярной примеси (лекарства) в несущей жидкости, что подтверждает идею, высказанную в [1, 2], о перспективности метода генерирования магнитоиндуцированных течений в тромбированных кровеносных сосудах для интенсификации в них транспорта тромболитиков.

Работа выполнена при финансовой поддержке РФФИ, грант 21-52-12013.