Теоретические основы химической технологии, 2020, T. 54, № 1, стр. 105-113

Применение аналитических решений для оценки вариабельности распределения концентраций компонентов специфических реакций в микрофлюидных системах

А. Л. Буляница a*, К. И. Белоусов b, А. А. Евстрапов a

a Институт аналитического приборостроения РАН
Санкт-Петербург, Россия

b Университет ИТМО
Санкт-Петербург, Россия

* E-mail: antbulyan@yandex.ru

Поступила в редакцию 10.06.2019
После доработки 27.08.2019
Принята к публикации 28.09.2019

Полный текст (PDF)

Аннотация

Ряд методов анализа жидких сред основан на проведении специфических реакций с определяемым компонентом. В микрофлюидных устройствах из-за ламинарности потока сложно организовать равномерное перемешивание компонентов реакции, что ухудшает воспроизводимость результатов реакций и качество анализа. Для повышения эффективности пассивного перемешивания в микрофлюидных системах следует использовать канал ввода пробы переменной геометрии, в частности, с изгибами, либо с изменяемой шириной. Эффективность перемешивания непосредственно связана с вариабельностью концентраций реагента в рабочей области, а эта характеристика адекватно оценивается через дисперсию концентрационного профиля. Для получения оценки дисперсии концентрационного распределения предложено приближенное аналитическое решение уравнения конвективной диффузии по методу моментов. Для этого система ввода пробы – прямолинейный участок постоянной ширины с дополнительным линейно расширяющимся каналом приближенно заменена эквивалентным каналом постоянной ширины, величины конвективной скорости и коэффициента диффузии приближенно заменены на усредненные. Оценки дисперсии, полученные на основе приближенного аналитического решения, имеют приемлемое совпадение с результатами численного 2D-моделирования конвективного массопереноса в COMSOL MULTYPHYSICS.

Ключевые слова: компоненты реакции, конвективно-диффузионный массоперенос, изгибающийся канал, диффузор, гидравлическое сопротивление, концентрационный пик, метод моментов, COMSOL MULTYPHYSICS

ВВЕДЕНИЕ

Проблематике описания процесса перемешивания и регулирования его интенсивности посвящено множество публикаций, как отечественных, так и зарубежных. Наиболее близкие к рассматриваемой задаче публикации касаются а) различных схем пассивного и активного перемешивания реагентов [14], б) формирования течений [5] и в) моделей распределения вещества [6]. Последнее направление касается оценки правомерности использования гипотезы гауссового распределения вещества. При этом работы [3, 4] связаны с активным перемешиванием.

Поскольку оценка эффективности перемешивания (или равномерности распределения концентрации по реакционной области) связана, в том числе и с дисперсией этого распределения, то естественным представляется использование метода моментов, оперирующего интегральными характеристиками распределения концентраций. Несмотря на то, что основные идеи этого метода предложены достаточно давно [7], его по-прежнему используют для моделирования процесса конвективно-диффузионного переноса вещества. В частности, учтено дополнительное размывание пробы при изгибе канала [8, 9], а также различные типы жидкости, характеры течения, соотношения конвекции и диффузии и размерности задачи [1015].

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

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

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

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

ТЕОРЕТИЧЕСКАЯ ЧАСТЬ

Профиль конвективной скорости в прямолинейном канале полагаем пуазейлевым параболическим, так как в рассматриваемых условиях поперечные числа Рейнольдса не превосходят 1, а продольные – 10.

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

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

При загрузке анализируемой пробы в микрофлюидное устройство исходное распределение компонентов пробы полагаем в форме пробки, описанной ранее. Расчеты профилей скорости в изгибающемся канале постоянного сечения и для диффузорного (конфузорного) течения позволят оценить число Пекле, определяемое отношением конвективной скорости и коэффициента диффузии, среднюю величину расхождения между скоростными профилями на выходе из диффузора и установившимся пуазейлевым и, как следствие, величину дополнительного поперечного конвективного переноса вещества, усиливающего интенсивность процесса выравнивания концентрации реагентов. Результаты оценочных аналитических расчетов соотнесены с точными численными решениями модифицированного уравнения Навье–Стокса (уравнения конвективной диффузии) в 2D-приближении, выполненные с использованием программного пакета COMSOL MULTYPHYSICS. Численное решение профиля концентрации получено в два этапа: сначала решалась стационарная задача для нахождения профиля скорости, а затем отдельно рассчитывалось движение пробки пробы. Также была учтена осевая симметрия каналов с условием второго рода на оси.

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

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

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

Расчет дисперсии распределения концентрации проводится в нормированных параметрах – скорости, коэффициенты диффузии, длины пробки и т.п. Масштабом длины, удобным для использования в 2D-задачах, является полуширина канала. В нашем случае это будет эквивалентная (усредненная) полуширина Hs. Обоснованием сведения задачи к 2D-представлению является оценка значения характеристического числа Дина, которое должно превышать минимальное пороговое значение. Известными полагаем профиль конвективной скорости в канале и коэффициент диффузии вещества в среде. Данные рассуждения составляют первую теоретическую предпосылку к решению задачи.

Скоростной профиль в канале постоянной ширины может быть разложен по базисным функциям –$~\cos \left( {\pi {\text{j}}z} \right)$, где z – относительная координата сечения, отсчитываемая от оси и нормированная на полуширину с коэффициентами разложения Bj. Здесь $C\left( {x,z,t} \right)$ – концентрация вещества в момент времени t в точке с маршевой (аксиальной) координатой х.

Моменты порядка n определены как

$\begin{gathered} {{\mu }_{n}}(z,t) = \int\limits_{ - \infty }^{ + \infty } {{{x}^{n}}C(x,z,t)dx} ;\,\,\,\,\frac{{\partial {{\mu }_{n}}(0,t)}}{{\partial z}} = 0, \\ \frac{{\partial {{\mu }_{n}}(1,t)}}{{\partial z}} = 0. \\ \end{gathered} $

Момент μ0 имеет смысл количества вещества, отношение η = μ10 определяет центр тяжести аналитического сигнала, его дисперсия равна ${{\sigma }^{2}} = {{{{\mu }_{2}}} \mathord{\left/ {\vphantom {{{{\mu }_{2}}} {{{\mu }_{0}}}}} \right. \kern-0em} {{{\mu }_{0}}}} - {{\eta }^{2}}.$

Нестационарное уравнение конвективной диффузии (модифицированное уравнение Навье–Стокса) преобразуется известным способом [7] к рекуррентной последовательности для моментов порядка n. Для моментов в нормированных параметрах оно примет вид

$\begin{gathered} \frac{{\partial {{\mu }_{n}}}}{{\partial t}} - D{\text{*}}\frac{{{{\partial }^{2}}{{\mu }_{n}}}}{{\partial {{z}^{2}}}} = nu{\text{*}}\left( {1 - {{{\left| z \right|}}^{m}}} \right){{\mu }_{{n - 1}}} + \\ + \,\,n\left( {n - 1} \right)D{\text{*}}{{\mu }_{{n - 2}}};\,\,\,\,n = 0,1,2,3.... \\ \end{gathered} $

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

${{\mu }_{1}}(z,t) \approx u{\text{*}}{{B}^{0}}t + \sum\limits_{j = 1}^{ + \infty } {\frac{{u{\text{*}}{{B}^{{\text{j}}}}}}{{D{\text{*}}{{{(\pi {\text{j}})}}^{2}}}}\cos (\pi {\text{j}}z)} $

и

(1)
$\begin{gathered} \left\langle {{{\sigma }^{2}}} \right\rangle = \left( {2 + {\text{P}}{{{\text{e}}}^{2}}\sum\limits_{j = 1}^{ + \infty } {\frac{{{{B}^{{{\text{j}}2}}}}}{{{{{(\pi {\text{j}})}}^{2}}}}} } \right)D{\text{*}}t + \\ + \,\,\left( {\frac{{{{{\Delta }}^{2}}}}{3} - 1.5{\text{P}}{{{\text{e}}}^{2}}\sum\limits_{j = 1}^{ + \infty } {\frac{{{{B}^{{{\text{j}}2}}}}}{{{{{(\pi {\text{j}})}}^{4}}}}} } \right). \\ \end{gathered} $.

Для пуазейлевого профиля с показателем степени m, равным 2, суммы, входящие в выражение для дисперсии, являются суммами обобщенных гармонических рядов с четными показателями и равны соответственно 16/945 и 16/9450. Эти значения для сумм правомерно использовать, когда выполнено условие

(2)
$\exp \left( { - \frac{{D{{\pi }^{2}}t}}{{{{H}^{2}}}}} \right) \ll 1.$

Иначе потребуется внести корректировки и уменьшить суммы в выражении (1), заменив слагаемые Bj на меньшие:

(3)
${{B}^{{\text{j}}}}(1 - \exp ( - D{\text{*}}{{\pi }^{2}}{{{\text{j}}}^{2}}t)).$

Модель влияния изгиба канала на перенос пробы обоснована в [8, 9] и была использована для оценки дополнительного размывания пробы при других режимах. Проведена нормировка характеристик движения: u* = U/H – нормированная максимальная скорость конвективного движения, $w = 2H$ – ширина микроканала, D – коэффициент диффузии. D* получено делением D на квадрат характерного размера Н, а $w = 2H$ – ширина микроканала. Дополнительное размывание пробы объясняется удлинением канала, а также перемещением разных сечений с различной скоростью в поле вращений. Численно увеличение дисперсии по сравнению с (1) при радиусе средней кривизны изгиба (Rc) равно

(4)
$\Delta {{\sigma }^{2}} = 2D{\text{*}}\Delta t + \frac{{{\text{P}}{{{\text{e}}}^{2}}{{w}^{2}}}}{{15R{{c}^{2}}}}D{\text{*}}t - \frac{{17{\text{P}}{{{\text{e}}}^{2}}{{w}^{2}}}}{{1680R{{c}^{2}}}}.$

при больших временах движения t, для которых выполнено условие (2). Здесь Δt – увеличение времени прохождения дуги по сравнению с прямолинейной хордой. Ключевое значение имеет число Пекле Pe* = 〈Uw/D, где 〈U〉 – средняя скорость конвективного движения по изгибающемуся каналу. Для сравнения при расчете дисперсии при движении по прямолинейному каналу постоянной ширины в числителе выражения для числа Пекле стоит произведение максимальной скорости на полуширину.

Таким образом, второй предпосылкой для расчета дисперсии распределения пробы при ее конвективно-диффузионном переносе является метод моментов в форме уравнений (1)–(4).

Второй предпосылкой является реализация пассивного перемешивания реагентов с помощью двух конструкций камеры с подводящими каналами. Упрощенно конструктивные элементы представлены на рис. 1 [16]. Экспериментальные исследования были проведены совместно с лабораторией биолюминесцентных биотехнологий Сибирского федерального университета. Глубина каналов варьировалась от 0.2 до 0.5 мм. Характерные длины составляли порядка 1–10 мм.

Рис. 1.

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

Третье исходное положение касается интерпретации скоростного профиля в линейно расширяющемся канале. Приближенно использовано диффузорное (конфузорное) течение. Эта теоретическая модель допускает аналитическое решение, представленное, например, в книге [17]. Аналогичная задача с указанным приближением физической модели была ранее рассмотрена в [18], где моделировались процессы массопереноса в проточном микрофлюидном чипе для проведения иммунной реакции связывания, пришитого к подложке биотинилированного 4H1b с меченым стрептавидином (Str-Cy5).

РАСЧЕТНЫЕ ОЦЕНКИ ПАРАМЕТРОВ КОНВЕКТИВНО-ДИФФУЗИОННОГО МАССОПЕРЕНОСА

Оценка эффекта дополнительного размывания пробы за счет изгиба канала. Если при параболическом пуазейлевом профиле второе линейное слагаемое равно $\frac{{16{\text{P}}{{{\text{e}}}^{2}}}}{{945D{\text{*}}t}},$ то получается, что более эффективное размывание пробы за счет изгиба канала реализуется при ${{\text{Н}} \mathord{\left/ {\vphantom {{\text{Н}} {Rc}}} \right. \kern-0em} {Rc}} > 0.252.$ В рассмотренном варианте конструкции (рис. 1) это ограничение не выполняется. Эффективность изгибающегося (серпантинного) канала в данном случае объясняется не увеличением размывания пробы за счет движения в интенсивном поле вращения, а за счет увеличения длины канала и времени прохождения пробы по нему.

Оценка дисперсии концентрационного пика и коэффициента вариации распределения концентрации реагента. Характеристикой эффективности перемешивания вещества может служить коэффициент вариации KV, определяемый как $~{{\left\langle \sigma \right\rangle } \mathord{\left/ {\vphantom {{\left\langle \sigma \right\rangle } {{{C}_{0}}}}} \right. \kern-0em} {{{C}_{0}}}},$ где C0 – среднее значение концентрации вещества. При допущении, что концентрационный пик имеет форму гауссовой кривой, строится зависимость между KV и отношением величины характерного размера области к 〈σ〉, равном А. Эта зависимость иллюстрируется данными табл. 1. Φ0 – интеграл вероятностей.

Таблица 1.  

Зависимость коэффициента вариации KV от характеристики размера концентрационного пика А

А Ф0(А) С0 σ KV
0.5 0.1915 0.960 0.021 0.02
1 0.3413 0.856 0.122 0.14
1.5 0.4332 0.724 0.216 0.30
2 0.4773 0.58 0.288 0.48
2.5 0.4938 0.495 0.331 0.67
3 0.4986 0.417 0.349 0.84
4 0.5000 0.313 0.351 1.12
5 0.5000 0.251 0.338 1.35
6 0.5000 0.209 0.323 1.54
7 0.5000 0.179 0.308 1.72
8 0.5000 0.157 0.294 1.87

Уменьшение дисперсии распределения концентрации по области, начиная со значения, примерно соответствующего значению А, равному 4, объясняется следующим: фрагмент концентрационного пика, в пределах которого существенно изменяется концентрация, будет занимать все меньшую часть области. На оставшейся части концентрация практически нулевая. Тем самым средняя концентрация С0 уменьшается и, как следствие, средний квадрат отклонения концентрации там, где она нулевая, т.е. в большей части реакционной камеры или другой рассматриваемой области, от С0 уменьшается. Данные табл. 1 коррелируют с наблюдаемым в [16] увеличением параметра KV по мере увеличения глубины камеры с 0.2 до 0.5 мм. В этих условиях дисперсия пика не меняется, а характерный размер, который должен определяться на основе полуширины сечения канала Н и глубины, увеличивается с увеличением глубины. Тем самым отношение А, используемое в табл. 1, также увеличивается.

Оценка скоростного профиля в линейно расширяющемся канале. Расширяющийся канал можно приближенно описать как диффузор [17] с углом расширения π/6. Скорость вычисляется как

(5)
$V\left( \varphi \right) = {{6\gamma u\left( \varphi \right)} \mathord{\left/ {\vphantom {{6\gamma u\left( \varphi \right)} r}} \right. \kern-0em} r},$
где φ – угол, отсчитываемый от оси диффузора и ограниченный π/6 (стенка), γ – коэффициент кинематической вязкости. Вспомогательная функция u(φ) вычисляется с помощью интеграла

(6)
$2\varphi = \int\limits_u^{{{u}_{0}}} {\frac{{du}}{{\sqrt {({{u}_{0}} - u)({{u}^{2}} + (1 + {{u}_{0}})u + q)} }}} .$

При этом

(7)
$\begin{gathered} ~\frac{\pi }{{~3}} = \int\limits_0^{{{u}_{0}}} {\frac{{du}}{{\sqrt {\left( {{{u}_{0}} - u} \right)\left( {{{u}^{2}} + \left( {1 + {{u}_{0}}} \right)u + q} \right)} }}} \,\,\,\,{\text{и}} \\ {R \mathord{\left/ {\vphantom {R 6}} \right. \kern-0em} 6} = \int\limits_0^{{{u}_{0}}} {\frac{{udu}}{{\sqrt {({{u}_{0}} - u)({{u}^{2}} + (1 + {{u}_{0}})u + q)} }}} , \\ \end{gathered} $

что позволяет найти две постоянных интегрирования u0 и q. Величина R в выражении (7) пропорциональна отношению количества жидкости к коэффициенту динамической вязкости и является аналогом числа Рейнольдса [17]. Соответствующее распределение u(φ) представлено в табл. 2.

Таблица 2.  

Данные для расчета профиля скорости в линейно расширяющемся канале согласно модели диффузора

φ, град u(φ)
0 0.300
10 0.223
15 0.169
20 0.128
25 0.064
30 0

Расширяющийся канал аппроксимирован диффузором, в котором радиальные перемещения r по оси лежат в пределах от 1 до 5 мм. Характерное значение скорости на оси канала при входе в реакционную камеру, т.е. на расстоянии 5 мм, примерно равно 0.08 мм/с. При этом скорость на оси прямолинейного участка канала шириной 1 мм составляет V(0) = 0.4 мм/с. Скоростной профиль рассчитан на основе формул (5)–(7). Время “прохода” по оси диффузора, исходя из начального условия r = = 1 мм, при конечной точке – 5 мм, определяется как td = 30.0 с, поскольку

${{r}^{2}} = 12\gamma {{u}_{0}}t + {{10}^{{ - 6}}}\,\,{{{\text{м}}}^{2}}.$

Численный расчет с использованием COMSOL дает оценку td = 29.0 с.

На выходе из диффузора построен профиль конвективной скорости. Для точки, соответствующей углу φ, расстояние будет 5/cos(φ), и на основе (5) можно построить распределение радиальной скорости диффузора и ее проекции в аксиальном (ось х) направлении. Это распределение по завершении начального участка длины Ln трансформируется в параболический профиль Пуазейля. Согласно уравнению неразрывности имеется дополнительный перпендикулярный к поступательному аксиальному движению жидкости конвективный перенос вещества, усиливающий диффузионный перенос и, как следствие, увеличивающий дисперсию пика. Расхождение скоростных профилей по сечению канала иллюстрируется рис. 2. Значения скоростей приведены в условных единицах при равенстве расхода (средней скорости). Wd – нормированная скорость в диффузоре, Wp – нормированная скорость при пуазейлевом профиле.

Рис. 2.

Форма нормированного скоростного профиля: 1 – установившееся конвективное течение Wp; 2 – на выходе диффузора Wd.

Среднеквадратичная величина расхождения нормированных скоростных профилей δV составила 0.130 отн. ед. при нормировке на максимальное значение скорости. В этом случае зависимость разности относительных скоростей аппроксимировали полиномом третьей степени по углу φ (коэффициент детерминации равен 0.884). Более точная аппроксимация полиномом 4-й степени дает значение 0.133 (коэффициент детерминации повышается до 0.980). Абсолютное расхождение скоростей есть 0.133, умноженное на масштаб скорости – скорость на оси канала V(0) при выходе из диффузора, т.е. при r = 5 мм или 0.08 мм/с.

Ln можно рассматривать как аналог начального участка, и его длину можно приближенно оценить по формуле Шиллера ${{Ln} \mathord{\left/ {\vphantom {{Ln} H}} \right. \kern-0em} H} = 0.058\operatorname{Re} $. Однако в этом случае трансформируется не постоянный скоростной профиль, а существенно переменный. Кроме того, известная оценка Шиллера относится к трубе с заданным диаметром. Оценим длину участка, формально заменив диаметр трубы на двойную полуширину 2Н. Продольное число Рейнольдса при прямолинейном участке 15 мм примерно равно 6.

В этом случае по оценке Шиллера Ln/H примерно равно 0.35. Для оценочного расчета возьмем Ln вдвое-втрое большим, т.е. отношение Ln/H считаем равным 0.7 или 1.

Тогда по уравнению неразрывности среднюю величину VY – скорости потока, направленного к стенке и дополняющего поперечный диффузионный массоперенос, можно оценить из пропорции ${{V}_{Y}} = {{H\delta V} \mathord{\left/ {\vphantom {{H\delta V} {Ln}}} \right. \kern-0em} {Ln}}.$

Таким образом, величина VY примерно составляет 1.52 и 1.06 мкм/с для вышеуказанных отношений Ln/H. За время 29.0 с диффузионный перенос составит ${{\left( {D{{t}_{{\text{d}}}}} \right)}^{{1/2}}}$ или примерно 51 мкм, а конвективный 44 и 31 мкм соответственно. Коэффициент диффузии компонента равен 4.8 × 10–10 м2/с.

Пересчет суперпозиции “диффузия и поперечная конвекция” в эквивалентную диффузию можно выполнить, приравняв диффузионные расстояния. Тогда эквивалентный коэффициент диффузии для движения по диффузору равен

${{D}_{{\text{d}}}}t = {{\left( {\sqrt {Dt} + {{V}_{Y}}t} \right)}^{2}}.$

Расчет дисперсии концентрационного пика при наличии диффузора. Использована формула для средней дисперсии пика (1). Длины прямолинейного участка L варьируются: 3, 6, 10 и 15 мм. После линейного участка следует диффузор длиной Ld = 4 мм.

Предложена приближенная аналитическая модель – эквивалентного прямолинейного канала постоянной полуширины Hs с усредненными параметрами: коэффициент диффузии Ds и максимальная осевая скорость конвективного движения Vs. Профиль полагается пуазейлевым параболическим.

Общие времена движения: 37.0, 45.0, 55.7 и 69.0 с. Эти времена рассчитаны с помощью пакета COMSOL. Корректность расчетов подтверждается тем обстоятельством, что разность времен прямо пропорциональна разности длин прямолинейных участков.

Длина исходной пробки пробы задана как 1 мм. Для ее перевода в безразмерную величину Δ0 требуется деление на 2Hs.

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

(8)
$D{{t}_{{\text{s}}}} = {{\left( {D\left( {t - {{t}_{{\text{d}}}}} \right) + {{D}_{{\text{d}}}}{{t}_{{\text{d}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {D\left( {t - {{t}_{{\text{d}}}}} \right) + {{D}_{{\text{d}}}}{{t}_{{\text{d}}}}} \right)} t}} \right. \kern-0em} t}.$

Усреднение остальных характеристик проводится по формулам

$\begin{gathered} {{H}_{{\text{s}}}} = {{\left( {HL + {{H}_{{\text{d}}}}{{L}_{{\text{d}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {HL + {{H}_{{\text{d}}}}{{L}_{{\text{d}}}}} \right)} {\left( {L + {{L}_{{\text{d}}}}} \right)}}} \right. \kern-0em} {\left( {L + {{L}_{{\text{d}}}}} \right)}},\,\,\,\,{\text{где}} \\ {{H}_{{\text{d}}}} = H + {{L}_{{\text{d}}}}{\text{tg}}{{\left( {{\pi \mathord{\left/ {\vphantom {\pi 6}} \right. \kern-0em} 6}} \right)} \mathord{\left/ {\vphantom {{\left( {{\pi \mathord{\left/ {\vphantom {\pi 6}} \right. \kern-0em} 6}} \right)} 2}} \right. \kern-0em} 2}. \\ \end{gathered} $

Усреднение скорости проводится по формуле, аналогичной (8) с заменой коэффициентов диффузии на соответствующие скорости. Средняя скорость на оси диффузора Vd рассчитана как интегральное среднее от выражения (5) при φ = 0 и при интегрировании по r от 1 до 5 мм и равна 1.61 мм/с. Далее используется формула для расчета средней дисперсии (1) и, при необходимости, поправки (2).

Пример расчета средней дисперсии пика при выбранном наборе параметров. Для примера рассмотрим случай: длина прямолинейного участка L = = 10 мм. Отношение Ln/H для определенности возьмем 0.7. Эквивалентный коэффициент диффузии диффузора Dd будет равен 62.73 × 10–10 м2/с, Ds – 34.95 × 10–10 м2/с, Hs – 830 мкм, Vs – 2.75 мм/с, нормированная полудлина пробки Δ около 0.60. D* в формуле (1) считается с учетом замены скорости на Us, коэффициента диффузии на Ds и полуширины на Hs. Уточненное значение числа Пекле будет Pe = UsHs/Ds или 65.28. Последнее свидетельствует о ключевой роли конвекции в размывании пробы, даже несмотря на корректировку коэффициента диффузии в сторону повышения. Коэффициент B1 в суммах, входящих в (1), должен быть уменьшен согласно (3). Таким образом, вместо 16/945 и 16/9450 или 1.693 × 10–2 и 1.693 × 10–3 будут 1.495 × 10–2 и 1.494 × 10–3. Три слагаемых в выражении для дисперсии будут соответственно 15.574; 0.121 и –9.556, что даст сумму 9.139. Тогда среднеквадратичное отклонение будет 3.023, а в абсолютных единицах после домножения на Hs получим около 2.52 мм. Расчет с помощью COMSOL дает сопоставимую оценку, а именно 2.86 мм. Сопоставление оценок величин 〈σ〉 для различных длин линейных участков (L) и вариантов выбора отношения Ln/H, полученных на основе приближенных аналитических решений (1)–(3), в сравнении с результатами численного моделирования с применением пакета COMSOL приведены на рис. 3.

Рис. 3.

Оценки стандартного отклонения при различных длинах прямолинейного участка: 1 – численное решение с использованием COMSOL MULTYPHYSICS; 2 – приближенное аналитическое решение при Ln/H = 0.7; 3 – приближенное аналитическое решение при Ln/H = 1.0.

ЗАКЛЮЧЕНИЕ

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

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

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

Расхождения, по-видимому, связаны как с особенностью счета в случае малых диффузионных времен (и малых дисперсий пробы), так и с принятыми приближениями модели: рассмотрение 2D-приближения, недостаточная обоснованность замены линейно расширяющегося канала диффузором, неточная оценка длины начального участка Ln и т.п.

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

Работа выполнена в рамках государственного задания № 075-00780-19-00.

ОБОЗНАЧЕНИЯ

А безразмерная характеристика размера реакционной области
Bj коэффициенты разложения нормированного скоростного профиля Wp по базисным функциям $~{\text{cos}}\left( {{\pi jz}} \right)$
$C\left( {x,z,t} \right)$ концентрация вещества, моль/м3
С0 среднее значение концентрации, моль/м3
D коэффициент диффузии, м2
D* = D/Н2 нормированный коэффициент диффузии, с–1
Dd модельное значение коэффициента диффузии с учетом поперечного конвективного движения, м2
Ds усредненное значение коэффициента диффузии с учетом движения по диффузору, м2
H полуширина прямолинейного канала, мкм
Hd средняя полуширина канала в диффузоре, мкм
Hs усредненная (эквивалентная) полуширина канала, мкм
KV = ${{\left\langle \sigma \right\rangle } \mathord{\left/ {\vphantom {{\left\langle \sigma \right\rangle } {{{C}_{0}}}}} \right. \kern-0em} {{{C}_{0}}}}$ коэффициент вариации концентрации
L длина прямолинейного участка, мм
Ln длина начального участка установления пуазейлевого скоростного профиля, мм
Ld длина диффузора, мм
m показатель степени симметричного профиля конвективной скорости (m = 2)
R величина, пропорциональная расходу жидкости
Rc средний радиус кривизны изгибающего канала, мкм
r радиальное перемещение, мм
r(0) начальное радиальное перемещение, мм
t время, с
td время движения по диффузору, с
Δt приращение времени движения за счет изгиба канала, с
U максимальное значение конвективной скорости движения по прямолинейному участку, м/с
U среднее значение конвективной скорости движения по изогнутому участку, м/с
Ud средняя скорость движения по оси диффузора, м/с
Us эквивалентное значение конвективной скорости, м/с
u(φ) вспомогательная функция
u0, q постоянные интегрирования
V(φ) скорость в диффузоре, м/с
VY значение поперечной конвективной скорости, м/с
δV среднеквадратичное расхождение нормированных скоростных профилей, м/с
Wd нормированная скорость в диффузоре
Wp нормированная скорость при пуазейлевом профиле
$w = 2H$  ширина микроканала, мкм
x аксиальная координата, м
z = x/H относительная координата сечения
γ коэффициент кинематической вязкости, м2
Δ нормированная на полуширину канала половинная длина пробки пробы
η отношение μ10
μ коэффициент динамической вязкости, Па с
μn последовательность моментов порядка n (n = 0, 1, 2)
σ2 дисперсия концентрационного пика
〈σ2 среднее значение дисперсии концентрационного пика
Δσ2 приращение дисперсии концентрационного пика за счет изгиба канала
Φ0 интеграл вероятностей
φ угол, отсчитываемый от оси диффузора, град
Pe = UsHs/Ds число Пекле
Pe* = 〈Uw/D число Пекле для изогнутого канала
Re = UsL продольное число Рейнольдса

ИНДЕКСЫ

0  индекс в средней концентрации (С0) и в выражении для нормированной длины пробки пробы (Δ0)
d диффузор
j номер коэффициента разложения по базисным функциям (от 0 до ∞)
max максимальное значение
n порядок момента (n = 0, 1, 2)
s усредненное значение показателя

Список литературы

  1. Барабаш В.М., Абиев Р.Ш., Кулов Н.Н. Обзор работ по теории и практике перемешивания // Теор. осн. хим. технол. 2018. Т. 52. № 4. С. 367.

  2. Yang Z., Méheust Y., Neuweiler I., Hu R., Niemi A., Chen Y.-F. Modeling Immiscible Two-Phase Flow in Rough Fractures From Capillary to Viscous Fingering // Water Resour. Res. 2019. V. 55. № 3. P. 2033.

  3. Zhao J., Ning Z., Lv M. Experimental study on the two-phase flow pattern and transformation characteristics of a flow mixing nozzle under a moderate flow rate // Meccanica. 2019. V. 54. № 8. P. 1121.

  4. Семёнов И.А., Ульянов Б.А., Кулов Н.Н. Экспериментальная оценка влияния колебаний на скорость массоотдачи от плоской поверхности // Теор. осн. хим. технол. 2016. Т. 50. № 3. С. 239.

  5. Лобасов А.С., Минаков А.В., Рудяк В.Я. Исследование режимов течения неньютоновских жидкостей со степенной реологией в микромиксере Т-образного типа // Теор. осн. хим. технол. 2018. Т. 52. № 3. С. 341.

  6. Доломатов М.Ю., Казаков М.А., Журавлева Н.А. Имитационное моделирование многокомпонентных стохастических систем // Теор. осн. хим. технол. 2017. Т. 51. № 5. С. 555.

  7. Aris R. On the dispersion of a solute in a fluid flowing through a tube // Proc. R. Soc. London, Ser. A. 1956. V. 235. P. 67.

  8. Wang J., Lin Q., Mukherjee T. Analytical dispersion models for efficient simulations on complex microchip electrophoresis system // Proc. 7th International Conference on Miniaturized Chemical and Biochemical Analysts Systems. MicroTAS’03. Squaw Valley, California, 2003. V. 1. P. 135.

  9. Magargle R., Hoburg J.F., Mukherjee T. An injector component model for complete microfluidic electrokinetic separation systems // NSTI-Nanotech., 2004. V. 1. P. 77.

  10. Barton N.G. On the method of moments for solute dispersion // J. Fluid Mech. 1983. V. 126. P. 205.

  11. Goltz M.N., Roberts P.V. Using the method of moments to analyze 3-dimensional diffusion-limited solute transport from temporal and spatial perspectives // Water Resour. Res. 1987. V. 23. № 8. P. 1575.

  12. Govindaraju R.S., Das B.S. Moment Analysis for Subsurface Hydrologic Applications (Water Science and Technology Library. V. 61). Dordrecht: Springer, 2007.

  13. Fischer H.B., List J.E., Koh C.R., Imberger J., Brooks N.H. Mixing in Inland and Coastal Waters. New York: Academic, 2013.

  14. Giona M., Adrover A., Cerbelli S., Garofalo F. Laminar dispersion at high Péclet numbers in finite-length channels: Effects of the near-wall velocity profile and connection with the generalized Leveque problem // Phys. Fluids. 2009. V. 21. P. 123601.

  15. Vedel S., Hovad E., Bruus H. Time-dependent Taylor–Aris dispersion of an initial point concentration // J. Fluid Mech. 2014. V. 752. P. 107.

  16. Лукьяненко К.А., Денисов И.А., Якимов А.С., Есимбекова Е.Н., Белоусов К.И., Букатин А.С., Кухтевич И.В., Сорокин В.В., Евстрапов А.А., Белобров П.И. Аналитические и ферментативные реакции в микрофлюидных чипах // Биотехнология. 2016. Т. 32. № 5. С. 69.

  17. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика: учебное пособие. Т. 6. Гидродинамика. М.: Наука, 2001.

  18. Esikova N.A., Bulyanitsa A.L., Klotchenko S.A., Taraskin A.S., Evstrapov A.A. The influence of hydrodynamic conditions and chamber geometry on the analytical signal from integrated biochips // J. Phys.: Conf. Ser. 2018. V. 1135. P. 012021.

Дополнительные материалы отсутствуют.