CC BY
11ISSNG868-5886
НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2GG3, том 13, № 4, c. 28-4G
ОРИГИНАЛЬНЫЕ СТАТЬИ
УДК 543.54+543.426+678.12
© А. Л. Буляница, А. А. Евстрапов, Г. Е. Рудницкая
МЕТОД МОМЕНТОВ ПРИ РАСЧЕТЕ ПАРАМЕТРОВ КАНАЛОВ В МИКРОРАЗМЕРНЫХ СИСТЕМАХ
В работе рассмотрены результаты математического моделирования процессов массопереноса в микрораз-мерной системе — канале микрофлюидного чипа при электрофоретическом разделении веществ. На основе метода моментов вычислены длины каналов и времена анализов, требуемые для получения необходимого разделения двухкомпонентной смеси. Рассмотрены и проанализированы модели, описывающие влияние эффекта неравномерности заполнения канала пробой. Приведены оценки профиля конвективной скорости потока вещества на основе системы уравнений Пуассона—Больцмана при упрощающей гипотезе об отсутствии в канале градиентов теплового поля.
ВВЕДЕНИЕ
В последнее время в аналитической технике все чаще используются микро- и нанотехнологии, позволяющие создать современные высокоэкспрессные системы анализа веществ. В частности, интенсивное развитие получили аналитические системы, известные как "ц-TAS" — microTotal Analytical System или "lab-on-a-chip". Основой таких систем является микроразмерная система — микрофлюидный чип, на котором осуществляются все действия с микроколичеством жидкой пробы, включая ввод и дозирование пробы, перемешивание, смешивание с реагентами, реализацию химических реакций, сепарацию полученного продукта, детектирование, сбор фракций. При изготовлении микроразмерных систем (МС) используются прогрессивные методы микро- и нанотехнологий (методы микроэлектронной техники, методы прецизионной литографии, LIGA-технологии и т. д.). Микрофлюидный чип (МФЧ) представляет собой компактное устройство с разветвленной системой микроканалов, микрососудов, реакторов и прочих функциональных элементов. В процессе анализа каналы и сосуды МФЧ заполняются пробой, реагентом, буфером, полимерами. Управление движением микропотоков жидкости осуществляется гидравлическим или электрокинетическим способом. При этом в каналах МФЧ вещество (проба, компонента смеси, реагент) осуществляет сложное (конвективное и диффузионное) движение, что приводит к эффектам, оказывающим влияние как на движения микропотоков, усложняя процессы управления ими, так и на электрофоретическое разделение пробы, искажая результаты анализа.
Данная работа посвящена анализу и моделированию процессов, происходящих в каналах микрофлюидного чипа при электрофоретическом анализе сложной смеси. Одной из основных задач исследования являлась оценка критической (минимальной) рабочей длины канала и времени анализа при заданных требованиях к разделению смеси.
ГЕОМЕТРИЧЕСКИЕ (КОНСТРУКТИВНЫЕ) ПАРАМЕТРЫ КАНАЛА МФЧ И ФИЗИЧЕСКИЕ ХАРАКТЕРИСТИКИ АНАЛИЗИРУЕМЫХ ОБЪЕКТОВ
В качестве исходных условий полагаем:
• что канал МФЧ является прямолинейным и прямоугольного сечения с полушириной к = = 20 мкм и глубиной, равной к (соотношение ширины к глубине, как 2 к 1);
• длина канала Ь является определяемой величиной, но она многократно превосходит к: Ь / к > 102-103;
• средние значения конвективной скорости и = 10-2-10-1 см/с;
• коэффициент диффузии вещества В имеет порядок 10-8-10-7 см2/с.
В зависимости от электрофоретической и электроосмотической подвижностей веществ определяется требуемая напряженность электрического поля, позволяющая достичь указанных величин конвективной скорости. Дополнительным ограничительным предположением является то, что напряженность электрического поля не достигнет столь больших значений, при которых тепловая мощность электрического поля превысит критическую величину 1 Вт/м.
БАЗОВАЯ МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
Прямолинейный канал прямоугольного сечения в первом приближении будем рассматривать как совокупность плоских слоев равной ширины и длины. Подобное сокращение пространственной размерности задачи с трех до двух не позволяет полностью адекватно описать процессы массопе-реноса в канале, поскольку в реальности имеет место массообмен между слоями. Тем не менее, для сечений, достаточно удаленных от дна канала, при равномерном по всей глубине заполнении канала буфером и пробой такое приближение вполне допустимо. Кроме того, величина проекции основных сил, действующих на вещество, на вертикальную ось мала по сравнению с проекцией на направление маршевой координаты (ось канала).
Обозначив аксиальную (маршевую) координату х, перпендикулярную стенке координату у и время ґ, применим пространственное двумерное уравнение Навье—Стокса, описывающее массоперенос вещества в каждом слое без учета реакции, в форме
зг дС п
-----+ и (у)--------------= В
дґ дх
д2С д2С
дх2 ду2
(1)
Здесь С = С(х, у, 0 — распределение концентрации вещества; и(у) — профиль конвективной скорости; В — коэффициент диффузии вещества. Гра-
дС( х,0, t)
ничные условия: симметрия
ду
= 0 и не-
проницаемость стенки условие
тественное
дС (х, к, ґ) дУ
недоступности
в форме
д РС (±тс, у, ґ) дхр
= 0, р = 0,1,2,... Послед-
нее условие говорит о невозможности ухода анализируемого вещества на бесконечное (очень большое) расстояние от зоны анализа. Начальные условия связаны с характером заполнения канала пробой и с объемом пробы.
Используем нормированные координаты и параметры г = у / к; и * = и / к; В * = В / к2. Здесь и — максимальная скорость конвективного движения (на оси канала). Полагаем скоростной профиль и(г) = и*(1 -|г|т); ге [0;1]. Параметр т, характеризующий степень клиновидности профиля, априорно неизвестен. Примем, что т > 2 (случай т = 2 — вариант параболического профиля, известного в хроматографии). Моделирование скоростного профиля и оценка параметра т требуют решения уравнения Пуассона—Больцмана для электростатического потенциала, температуропроводности и профиля конвективной скорости.
Решение этой задачи при упрощающем положении о равенстве температуры по всему каналу будет приведено в разделе "Моделирование скоростного профиля электроосмотического потока".
Для решения уравнения (1) воспользуемся методом моментов, предложенным в работе [1]. Решение аналогичной задачи применительно к хроматографии описано во многих работах, например
в [2].
Введем моменты порядка п, определяемые как
-с со
Чп (г, ґ) = |хпС(х, г, ґ)йг.
(2)
По аналогии с теорией вероятностей данные моменты следовало бы называть начальными моментами соответствующего порядка.
Момент имеет смысл количества вещества, отношение П=М-1/М0 определяет центр тяжести аналитического сигнала, его дисперсия определяется как а2 = Ч2 /Ц0 — ц2. Вычисление третьего начального, а затем и центрального моментов, позволяет оценить асимметрию аналитического сигнала (распределения вещества). Центральный момент четвертого порядка дает оценку, аналогичную коэффициенту эксцесса распределения, что используется при проверке гипотезы о гауссовом характере распределения вещества.
Исходное уравнение для моментов в нормированных параметрах примет вид:
= 0 , а также ес-
вещества
дЧп
дґ
- В
д Ч
дг2
= пи *(1 - |г|" )ъп-1 + п(п -1) в Уп-2;
п = 0,1,2,3,...
Граничные условия для моментов любого порядка являются граничными условиями второго
рода: = 0, = 0.
дг дг
Начальные условия определяются исходным распределением вещества, которое полагаем соответствующим "пробке" относительной длины 2Д с равномерным распределением вещества по сечению, то есть при х е [-Д; Д] ^ С(х, г,0) = 1/(2Д). Очевидно, что переход к естественным значениям по концентрации произойдет при учете множителя 2ДС0. Подобное нормирование удобно с точки зрения формирования начальных условий для моментов:
Ъ (г,0) = 1; й (г,0) = 0;
Ъ2( г,0) = Д2/3; ъ3( г,0) = 0;.
ОПРЕДЕЛЕНИЕ ЦЕНТРА ТЯЖЕСТИ И ДИСПЕРСИИ ПИКА ИНФОРМАТИВНОГО АНАЛИТИЧЕСКОГО СИГНАЛА
Решение полученных уравнений осуществляется по следующему алгоритму: для п = 0 уравнение является однородным и решается методом разделения Фурье; для п = 1, 2,... правая часть уравнения должна быть разложена по собственным функциям 7}(г) (они же — решения задачи Штурма—Лиувилля, соответствующие уравнению
для Ъ).
дй д Ъ0
При п = 0 имеем Т(г)
= В
и для функ-
В0 =
т
т+1
т + 22
Таким образом, на основании приведенной рекуррентной формулы можно продолжить последовательности коэффициентов разложения В1 для любого целого т (на основе данных табл. 1).
Эта формула верна и для нецелых значений т. Но начальные значения коэффициентов В1 нужно вычислять численно (оттабулировав для 0 < т < 2).
Для случая п = 1 ищем решение в форме
Ъ (г, t) = У0 ^^ У] ^) ^(п 1 г).
1 =1
Табл. 1. Коэффициенты ряда Фурье разложения профиля конвективной скорости
дt дг2
ций 1Дг) получаем уравнения: Т''+Х21-Т1- = 0;
Т(0) = 0, Т' (1) = 0. Поэтому собственные функции (решения задачи Штурма—Лиувилля) будут Т1 (г) = С08(п 1г); 1 = 0,1,2,.
Исходя из разложения начального условия в ряд Фурье по косинусам, й0 (г, 0 = 1. Этот результат означает равномерное распределение вещества по всей ширине канала.
Для последующих этапов решения требуется разложить конвективный скоростной профиль по функциям Т(г), то есть фактически применить разложение в ряд Фурье по косинусам.
В дальнейшем будем полагать, что
1 - |г|т = В0 + ^В1 С08(п /г).
1 =1
Коэффициенты Фурье определяются параметром т и приведены в табл. 1.
При больших т коэффициенты вычисляются:
т В0 В,
2 0.667 4(-1)1+1 /(П)2
3 0.750 6(-1)1+1 /(п)2 --12[ - (-1)1]/(п)4
4 0.800 8(-1)1+1 /(п)2 -- 48(-1)1+1 /(п)4
5 0.833 10(-1)1+1 /(п)2 --120(-1)1+1 /(п)4 + + 240[ - (-1)1]/(п)6
Система обыкновенных дифференциальных уравнений для функций V, включающая начальные условия, примет вид:
0 = и * • В0; У0(0) = 0 и
Л
■ + В • (п)2у7 = и • В1; У7 (0) = 0.
Л 1 1 '
Отсюда следует,
и * • В1
что
УП = и * • В 0 • V;
В * • (п)2
Среднее значение первого момента и, в нашем случае, центра тяжести пика, равно У0.
В случае больших диффузионных времен t оценка центра тяжести будет
и *
й( г, t) « и * • В0 • t + У С08(п’г).
М В • (п)2
При этом формула будет достаточно точна, если в первом слагаемом суммы (У\) функции ехр(...) достаточно малы, например, менее 0.01. Соответствующие диффузионные времена должны
1п(100) 0.467
удовлетворять требованию t >> —2 * ~ ^ .
В дальнейшем такие времена будем считать большими.
Возведение в квадрат и усреднение по г (по интервалу [0;1]) приводит к оценке
[1 - ехр(-В * • (П )2 • t)].
и
1 +^
^ ~(и •В0 •t)2 + —У
•В
2
=1 в*2 • (п)4
Учитываем, что
)со82(п'г)( = 1/2; )со8(п’г)( = 0.
Расчет момента второго порядка проводится по аналогичной схеме. Из-за громоздкости выражений приведем только величину среднего значения дисперсии, понимаемую в данном случае как среднее значение второго момента, из которого вычтена величина среднего квадрата центра тяжести.
В этих условиях больших диффузионных времен оценка величины дисперсии пика примет вид линейной функции
У{ = 2 + Ре2 У
В
2
, 1=1 (П)
V ■*
Г д2 3 +ц В12 ^
Ре2 У-----------------1
, 3 2 ^(пО4,
V 1
В * • t +
(3)
Здесь число Пекле определяется как Ре = = и /В . Воспользовавшись формулами для сумм
+л 1
рядов вида 8к = У—— [3], можно оценить соот-
1=1 1
4л В12
ветствующие коэффициенты 51 = У --------------— и
1 =1(п)
+л В12
для различных профилей (парамет-
н(п) ров т) (см. табл. 2).
Табл. 2. Характерные параметры, определяющие дисперсию пика для различных профилей (параметров т)
т *!-102 *2-103
2 1.693 1.693
3 1.389 1.357
4 1.108 1.057
5 0.890 0.829
Анализ выражения (3) позволяет сделать следующие выводы:
1) структура (3) качественно напоминает выражение для тейлоровской дисперсии хроматографического пика. То есть имеется слагаемое, соответствующее исходному "размытию" пробы,
(Д2 / 3) а также линейно возрастающие со временем — диффузионное (первое) и конвективнодиффузионное (второе слагаемое в первой скобке);
2) это совпадение неслучайно, учитывая сходство процессов конвективно-диффузионного мас-сопереноса в задачах хроматографического и электрофоретического разделений;
3) по мере возрастания клиновидности профиля (т) дисперсия пика убывает, что формально отражено в монотонном уменьшении суммы *1 с ростом т;
4) по мере роста параметров Д и и происхо-
дит возрастание дисперсии пика; зависимость от В не столь однозначная. В частности, в случае В * = и * •у[5^2 первое слагаемое (3) минимизируется. То есть получаем наличие нетривиального оптимального (конвективно-диффузионно-
согласованного) режима измерений.
В табл. 3 приведены результаты расчета параметров пика (центра тяжести и дисперсии) для широкого диапазона времен анализа. Расчетные параметры т = 5, В = 0.05, и = 7, Д= 5.
Видно, что параметры (центр тяжести и дисперсия пика) в различных точках сечения канала (г) не очень существенно отличаются друг от друга. Самым быстрым будет движение пика по оси канала, самым медленным — вблизи стенки. Кроме того, "отставание" центра тяжести пика, формируемого у стенки канала от пика, расположенного по оси, практически не накапливается с течением времени. Благодаря монотонному возрастанию дисперсий пиков уже при временах 100 и более обеспечивается практически полное перекрывание пиков со всех точек сечения. Тем самым для расчетов можно считать достаточным использование средней оценки (3).
Оценка по формуле (3) должна быть занижена. Особенно этот эффект существенен при малых временах. Непосредственный численный расчет значений вторых моментов (квадрата центра тяжести) для 101 равномерно отстоящих с шагом 0.01 значений г с последующим усреднением дает оценки: 72.78 ^ = 10), 161.20 ^ = 20), 249.80 ^ = = 30) и 427.20 ^ = 50). Относительная погрешность численного расчета лежит в пределах 1 %.
По указанной схеме, продолжая расчеты, можно оценить величину третьего момента, а далее третьего начального момента и коэффициента асимметрии А.
2
*
и
+
Табл. 3. Параметры пика (центр тяжести — сверху, дисперсия снизу курсивом)
Чч. ґ г 10 20 30 50 100
0.00 61.777, 120.139, 178.472, 295.139, 586.806,
64.607 152.586 240.835 417.337 858.590
0.40 59.936, 118.278, 176.611, 293.278, 584.944,
81.081 169.160 257.410 433.911 875.165
0.80 55.061, 113.371, 171.705, 288.371, 580.038,
68.524 156.522 244.771 421.272 862.525
0.95 53.633, 111.938, 170.271, 286.938, 578.604,
54.394 142.308 230.550 407.051 848.305
Среднее* 58.333, 116.667, 175.000, 291.667, 583.333,
72.18 160.40 248.62 425.06 866.16
* Под средней оценкой подразумевается расчет по формуле (3).
В свою очередь, коэффициент асимметрии определяется через центральный момент третьего
порядка как А = %3/&ъ. Учитывая нормировку
анализируемого вещества на единицу (Ъ = 1), центральный момент третьего порядка вычисляется как х3 = Ъ3 - 3ЪгЪ1 + 2(ъ )3, где моменты Ъ вычисляются по (2). Соответствующее дифференциальное уравнение в частных производных было сведено к системе ОДУ и решено численно для небольших времен t. Результаты приведены в табл. 4.
В целом по мере увеличения t, В и Д коэффициент асимметрии убывает до значений, порядка
0.01. Рост первых двух параметров свидетельствует об усилении процессов обмена между сечениями канала. Большие длины пробки формируют изначально симметричное распределение вещества, которое будет полностью доминировать на первом этапе движения вещества и компенсировать неравномерность конвективного потока. Аналогичный эффект достигается уменьшением конвективной скорости и . Таким образом, гипотеза о симметричности пика в первом приближении подтверждается.
ЗАСТОЙНЫЕ (ПРИСТЕНОЧНЫЕ)
ЗОНЫ ВВОДА ПРОБЫ
Существенной проблемой является наличие застойных зон вблизи стенок канала. Эти эффекты
можно учесть, предложив соответствующие модели начального распределения вещества. Ниже рассматривается влияние эффекта застойных зон на распределение вещества по каналу и параметры пика (центр тяжести и дисперсию).
Проанализируем модель с экспоненциальной загрузкой пробы. При этом предполагается исходное распределение вещества по сечению внутри "пробки" как
С
С (х, г,0) = 0
а
2Д 1 - ехр(-а)
ехр(-аг),
где а — параметр, характеризующий неравномерность распределения вещества по ширине канала. При а = 0 модель переходит к модели равномерного заполнения проб, и по мере роста а неравномерность исходного заполнения возрастает. При этом при г « 0 доля вещества в соответствующих сечениях превосходит 1, при г « 1 — существенно меньше 1.
Разложение исходной пробки вещества в ряд Фурье позволяет получить решение уравнения (1) для й0 как
г, ґ) = 1 +
2а 2
1 - ехр(-а)
і =1
1 - (-1)1 ехр(-а)
а2 + (п)2
ехр(-В * -п2 і 2ґ )соз(пг).
х
Табл. 4. Коэффициенты асимметрии пика при различных условиях ввода и движения вещества по каналу
Д * и В* ґ Л
Зависимость от длины пробки
1 10 0.1 3 0.1270
2 10 0.1 3 0.1178
3 ?? ?? ?? 0.1048
5 ?? ?? ?? 0.0756
10 ?? ?? ?? 0.0276
15 ?? ?? ?? 0.0105
Зависимость от конвективной скорости
5 3 0.1 3 0.0070
5 5 ?? ?? 0.0261
?? 8 ?? ?? 0.0585
15 0.1022
Зависимость от коэффициента диффузии
5 10 0.05 2 0.4482
?? ?? 0.1 ?? 0.1129
?? ?? 0.2 ?? -0.0108
0.3 -0.0438
Зависимость от времени анализа
5 10 0.1 1 0.1301
5 0.0259
10 -0.0347
Видно, что по мере роста t количество вещества в каждом сечении г выравнивается, стремясь к единице. Далее решаются уравнения для п = 1,2, и рассчитываются параметры пиков. В табл. 5 приведены данные, характеризующие динамику распределения вещества по сечениям канала для модельных параметров Д = 10; и = 10; В = 0.1; t = 1 и t = 3; а = 0, 0.5, 1 и 2.
На основании приведенных данных можно утверждать, что даже при достаточно неравномерном заполнении веществом канала (а = 2) влияние начальной неравномерности заполнения "пробки" мало сказывается при временах t > 3. То есть,
главным образом, на положение пика влияют параметры, определяющие размывание пика, а именно В и и . К моменту времени t = 3 даже в случае а = 2 происходит существенное выравнивание распределения вещества по сечениям канала до уровня (100±3 %). При этом распределение центров тяжести пиков по сечениям г практически не зависит от начальной структуры пробки вещества; соответствующие дисперсии пиков отличаются в пределах ±10 % для приведенных значений параметра а от 0 до 2. Например, для г = 0 при t = 3 в случае а = 0 центр тяжести и дисперсия пика будут 27.333 и 45.116, а в случае а = 2 — 27.981
Табл. 5. Динамика распределения вещества по сечениям канала
г а = 0 а = 0.5 а = 1 а = 2
Время ґ= 1
0 1.0000 1.0194 1.0706 1.2197
0.40 1.0049 1.0195 1.0627
0.95 0.9828 0.9342 0.7921
Время ґ= 3
0 1.0000 1.0026 1.0095 1.0299
0.4 1.0008 1.0029 1.0092
0.95 0.9975 0.9906 0.9705
Табл. 6. Распределение вещества по сечениям канала
г 5 II 0 5 = 0.1 .2 О II 5
Время ґ = 2
0 1.0000 1.0303 1.0648
0.30 ?? 1.0179 1.0382
0.80 ?? 0.9754 0.9474
0.95 0.9699 0.9357
Время ґ= 5
0 1.0000 1.0016 1.0034
0.30 ?? 1.0009 1.0020
0.80 ?? 0.9987 0.9973
0.95 0.9985 0.9967
и 47.819 соответственно. При этом колебания дисперсии пика для различных сечений г сильнее в случае а =0. Можно предполагать, что при дальнейшем росте времени анализа ^ эффект влияния изначального неравномерного распределения вещества ("экспоненциальные застойные зоны") еще сильнее нивелируется.
Таким образом, влияние возможного изначально неравномерного распределения вещества в пробе не является фактором, решающим образом влияющим на форму пика. При больших временах анализа допустимо применять гипотезу равномерного заполнения веществом. Кроме того, равномерность заполнения канала веществом пробы подтверждается и экспериментально.
Введем "барьерную застойную зону". В рамках этой модели полагаем начальное распределение вещества в "пробке" как
С 1
С (х, г,0) = —0----, 0 < г < 1 -5;
^ ' 2Д 1 -5
С(х,г,0) = 0, г > 1 -5.
Таким образом, устанавливаем "барьер" для закачивания вещества в "пробку", расположенный на относительной глубине канала 5, отсчитываемой от стенки. Этот слой исходно не содержит вещества. Разлагая исходный профиль концентраций в ряд Фурье и решая уравнение при п = 0, по-
лучим распределение вещества по сечению канала в форме
Мо( г, () = 1 +
2
1 -5
х
1 =1
(-1)1 + 8т(я;5)
П
ехр(-Б* -к212t)со$,(п]г).
Как и для предыдущей модели, очевидна тенденция выравнивания вещества по сечениям с ростом t. Результаты расчетов динамики распределения вещества для различных характеристик барьера 5 = 0.1, 0.2 и 0 (отсутствие барьера) и при заданных и = 7, Б = 0.1, А = 5, t = 2 и 5 сведены в табл. 6.
Анализ результатов моделирования выявляет для обеих рассмотренных моделей одни и те же тенденции: с ростом времени анализа влияние застойных зон практически не сказывается на установившемся распределении вещества по сечению канала, поскольку происходит его выравнивание; координаты центра тяжести и величина дисперсии пика практически стабилизируются на уровне одних и тех же значений.
По-видимому, специфические эффекты, связанные с наличием застойных зон вблизи стенок канала следует моделировать не подбором начальных условий, характеризующих недостаточное поступление вещества к стенкам канала, а подбором граничных условий (например, искусственная полу- или полностью непроницаемая стенка на расстоянии 5 от стенки канала).
РАЗРЕШЕНИЕ ДВУХ ПИКОВ
Рассмотрим два вещества, характеризующихся конвективной скоростью и 1,2 и коэффициентом диффузии Б 12. Полагаем, что сигнал каждой анализируемой компоненты смеси является гауссовской кривой, т. к. в первом приближении была подтверждена гипотеза малой асимметрии. То есть форма /-го пика имеет вид
/г (Х) =
1
■\[2яаг
ехр
(х - т/) 2а2
2 Л
i = 1,2.
(4)
V ‘ /
Здесь тг — центр тяжести, а2 — дисперсия г-го пика.
Далее нас будет интересовать разделение двух близких компонент смеси, то есть компонент,
удовлетворяющих требованию и* ~и*; Б* ~Б*. В этих условиях можно полагать, что дисперсии обоих пиков практически одинаковы, центры тяжести имеют смещение. Также в качестве упрощающего положения примем гипотезу о примерно
равном (соизмеримом) содержании обеих компонент смеси. То есть суммарная мощность информативного сигнала от каждой из компонент будет одинакова.
Сформируем суммарный сигнал
/кит (х) = /1 (х) + /2 (х) и определим координаты точек экстремумов. Необходимые требования по адекватному оцениванию сигнала и достоверному разделению пиков можно сформулировать как
а) х^ = т12, б) /(хтт) << /(хт2ах). Практически требуется, чтобы разница между уровнями минимального и максимального сигналов превосходила двойную амплитуду шума.
Можно показать, что для четкого разделения двух соседних равномощных пиков (при отношении сигнал/шум, например, 20:1 по амплитуде), требуется выполнение соотношения
|т1 - т2| > 2л[2а . В этом случае однозначно идентифицируются как оба максимума суммарного сигнала (практически в точках т^), так и минимум сигнала между пиками. То есть оба пика практически достоверно распознаются. В дальнейшем будем для простоты руководствоваться критерием 3 а.
Оценки центров тяжести и дисперсий пиков соседних компонент, полученные с помощью метода моментов, позволят оценить их динамику и определить требуемое для разделения пиков время в соответствии с выбранным критерием разделения.
В работе [4] приведены данные для конвективной скорости и коэффициента диффузии одноцепочечных фрагментов ДНК длиной 100, 200, 300, 400, 600 и 800 базовых пар, которые воспроизведены в табл. 7. Геометрические характеристики канала: глубина и полуширина канала равны Н = = 20 мкм. Объем пробы — 4АН3 или 32А пиколитров (пл).
Рассмотрим задачу оценки минимальной сепа-рационой длины Ь канала и соответствующего времени разделения соседних фрагментов смеси 100-200-300-400 Ьр, 400-600-800 Ьр.
Полагаем объем пробы 320 пл (А=10). В соответствии с (3) можно вычислить дисперсию пика. Среднее положение центра тяжести пика опреде-
* у» 0
ляется как и - В -1.
Оценим требуемую длину канала Ь по следующей методике: будем рассчитывать время прохода центра тяжести каждого пика через точку х = Ь; временными границами пика будем полагать время, за которое данный компонент проходит точки Ь ± 1.5а . Соответствующие времена пересчитываются через скорость прохождения центром тяжести пика требуемого расстояния 1.5 а. В табл. 8 представлены расчетные данные по разделению фрагментов 100, 200, 300 и 400 Ьр для точки х = Ь.
х
Табл. 7. Скорости и коэффициенты диффузии различных фрагментов ДНК (из [4] )
ДНК (Ьр) V, мм/с £>-105, мм2/с * и £>*-102
100 0.1238 1.530 6.19 3.825
200 0.0928 0.789 4.64 1.973
300 0.0744 0.486 3.72 1.215
400 0.0624 0.356 3.12 0.890
600 0.0488 0.239 2.44 0.598
800 0.0411 0.168 2.06 0.420
Табл. 8. Оценка разделения пиков 100-200-300-400 Ьр в каналах с различной рабочей длиной
Ь Фрагмент, Ьр Центр пика, ъ О Дг Время выхода пика 4±1.5Д^
800 100 155 37.30 7.2 144-166
?? 200 207 44.52 11.5 190-224
?? 300 258 50.39 16.3 233-283
?? 400 308 53.73 20.7 277-339
1000 100 194 41.80 8.1 182-206
?? 200 259 49.89 12.9 239-279
?? 300 323 56.53 18.2 296-350
?? 400 385 60.30 23.2 350-420
Табл. 9. Характеристики разделения фрагментов 400-600-800 Ьр в каналах различной длины
Ь Фрагмент, Ьр Центр пика Ъ О Дг Время пика 4±1.5Д?
1100 400 423 63.3 24.3 386-460
?? 600 541 68.0 33.4 491-591
?? 800 641 74.1 43.0 576-706
1500 400 577 74.2 28.5 534-620
?? 600 738 79.8 39.2 679-797
800 874 87.2 50.8 798-950
Итак, при рабочих длинах канала более 1000 (2 см) осуществляется хорошо выявляемое разделение всех 4 составляющих анализируемой смеси. При этом полное время анализа составит около 7 мин.
Аналогичную задачу решим применительно к разделению фрагментов 400-600-800 Ьр. Результаты приведены в табл. 9. Исходя из данных табл. 9, требуемая сепарационная длина канала для разделения фрагментов 400-600-800 Ьр — 1500 (или 3 см). При этом время анализа составит 16 мин.
В оценочную формулу (3) входят величины, определяемые коэффициентом т, характеризующим клиновидность профиля. Для приведенных выше расчетов данных табл. 8 и 9 взято среднее значение этого коэффициента т = 5. В принципе, скоростной профиль может обладать еще большей кли-новидностью (например, т = 6, 7 и более). В этом случае разрешение увеличивается, во-первых, за счет возрастания средней конвективной скорости (и абсолютной величины разности конвективных скоростей различных фрагментов), во-вторых, за счет уменьшения дисперсии (через коэффициенты «1 и Б2, которые убывают с ростом т).
Приведем сравнительный расчет для случая т = 6. Благодаря большему значению коэффициента Фурье В0 = 6/7 (вместо 5/6) и меньшим значениям остальных коэффициентов Bj и соответственно сумм «1 и «2 (равных 0.72610 2 и 0.66210 3) увеличивается абсолютная величина разницы скоростей движения центров тяжести различных пиков при уменьшении их дисперсии. Для случая разделения фрагментов 100-200-300-400 Ьр для Ь = 800 имеем временные положения пика (с учетом размаха ±1.5О): пик 100 (151±10), пик 200 (201± 15), пик 300 (251±21) и пик 400 (299±27). Полное разделение наблюдается для пиков 300 и 400 как критических для общего разделения. В данном случае рабочая длина канала должна быть не меньше 1.6 см. При этом время разделения смеси сокращается примерно до 330 с. Аналогичные данные по разделению фрагментов 400600-800 Ьр имеют вид (для Ь = 1250): пик 400 (467±34), пик 600 (598±46) и пик 800 (708±60). Отсутствует перекрывание пиков 600 и 800; рабочая длина канала сокращается до 2.5 см, и время анализа не превосходит 13 мин.
Моделированием пограничного (пристеночного) слоя на основе системы уравнений Пуассона— Больцмана позволит оценить параметр клиновид-ности профиля т и тем самым более точно спрогнозировать режим разделения.
МОДЕЛИРОВАНИЕ СКОРОСТНОГО ПРОФИЛЯ ЭЛЕКТРООСМОТИЧЕСКОГО ПОТОКА
Модель строится на основании уравнения температуропроводности, уравнения Пуассона—
Больцмана для электростатического потенциала и уравнения электроосмотического движения. Исходя из малых размеров сечения канала, примем упрощающее предположение об отсутствии градиента теплового поля. Иными словами, полагаем весь канал микрофлюидного чипа идеально температурно стабилизированным. Расчеты электростатического потенциала и профиля конвективной скорости проведены при этом упрощающем предположении.
Помимо введенного ранее нормированного параметра сечения г, отсчитываемого от оси канала полуширины к (г е [0;1]), проведем нормировку электростатического потенциала на величину дзета-потенциала Со, то есть С е [0;1].
Решаем систему уравнений Пуасссона— Больцмана при следующих значимых допущениях.
1) Температура в каждом сечении канала одинакова Т = 300 К (основываемся на микроразмерах сечения).
2) Полагаем достаточно малой концентрацию буфера, а именно концентрация С0 находится в диапазоне 10-6-10-4 М.
3) Относительная диэлектрическая проницаемость £ = 80 (воде соответствует значение 81) [5].
4) Величина дзета-потенциала Со не превосходит 100 мВ (лежит в пределах 10-100 мВ). [2]
В этой ситуации коэффициент в = = (Р£0)/(ЯТ) << 1 и sк(вC) ~вС; Р и Я обозначают соответственно постоянную Фарадея и универсальную газовую постоянную.
Уравнение Пуассона—Больцмана, описывающее распределение потенциала по сечению канала, линеаризуется (5) и решается аналитически
—^=я2С. (5)
—г2
Граничные условия: симметрии —^\ г=0 = 0
—г
и электростатического баланса С\ г=1 = 1, так как на стенке потенциал совпадает с дзета-потенциалом. Параметр Я определяется, исходя из выражения
Р2Г к2 £0еЯТ ’
где £0 — электрическая постоянная
[5];
107
£п =■
8.842 • 10 12 ф/м, с — скорость света
4яс в вакууме.
Распределение относительной величины потенциала описывается выражением
С = ск(Яг)/ ск(Я). Поскольку уравнение для скорости электроосмотического потока (скоростного
профиля) имеет вид
d2u 2EFh 2C0
sh(eZ) с гра-
h =
оптимального по
dz П
ничными условиями симметрии профиля на оси канала и прилипания к стенке канала (u = 0), то это уравнение линеаризуется и аналитически решается при тех же допущениях, что и уравнение Пуассона—Больцмана. Здесь Е — напряженность электрического поля, П — динамический коэффициент вязкости.
Основные закономерности поведения решения:
а) по мере роста параметра Я степень неравномерности распределения потенциала по сечению канала возрастает;
б) поскольку вторая производная (радиус кривизны) скоростного профиля пропорциональна потенциалу, неравномерность распределения потенциала приведет к усилению неравномерности скорости по сечению канала;
в) аппроксимировав профиль выражением /1 I \т\
umax (1 - z ), где т — показатель, характеризующий клиновидность профиля, можно связать т c заданными условиями, а именно дзета-потенциалом Со, шириной канала 2h, напряженностью электрического поля Е, концентрацией буфера С0, коэффициентом динамической вязкости П Для этого нужно решить задачу оптимизации по критерию минимума среднего квадратичного отклонения профилей друг от друга.
При использовании допущений линейной модели (5) в уравнениях Пуассона—Больцмана и электроосмотического скоростного профиля показатель клиновидности профиля т не зависит от величины дзета-потенциала, коэффициента динамической вязкости и от напряженности электрического поля. От этих параметров зависит абсолютная величина скорости, а не ее распределение по сечению канала. Однако при неприменимости линейной аппроксимации гиперболического синуса, то есть при больших значениях параметра в, т должно зависеть и от величины дзета-потенциала, и от коэффициента динамической вязкости, и от напряженности электрического поля.
Влияние концентрации буферного раствора С0 рассматривается при следующем допущении: концентрации столь малы, что их изменение в ранее указанных пределах практически не изменяют диэлектрическую проницаемость буфера, примерно
равную соответствующему значению для воды. То же самое замечание относится и к динамическому коэффициенту вязкости.
Заданы исходные данные:
П = 10 5 см2/с, Со = 10 мВ, Е = 104 В/м = 20 мкм. Численный расчет оптимального критерию минимума расхождения профилей коэффициента т представлен в табл. 10.
Влияние ширины канала иллюстрируется данными табл. 11 при значениях величин: П =
= 10-5 см2/с, Со = = 10 мВ, Е = 104 В/м, С0 =
= 1.510-5 М.
Табл. 10. Оценка параметра клиновидно-сти профиля т для буферных растворов различной концентрации
С0-105 т
0.1 2.2
0.3 2.7
1 4.2
2 5.8
3 7.0
5 >7
Табл. 11. Оценка параметра клиновидности профиля т для каналов различной ширины
Таким образом, сделанные ранее при анализе разделения фрагментов смесей допущения относительно возможных значений параметра т обоснованы. Можно подобрать условия, при которых степень клиновидности профиля будет характеризоваться значением т от 2 до 7 и более.
Зависимость параметра клиновидности т от концентрации буфера С0 и от полуширины канала к
В рассмотренной модели предполагалось прямоугольное сечение канала, аппроксимируемое совокупностью слоев одной ширины 2к. Переход к пространственно двумерной задаче даже для этого случая потребует наложения аналогичного условия электростатического баланса £ = 1 для дна канала.
Решение указанной задачи, его обобщение на случай трапецеидального сечения а также учет влияния нелинейных эффектов в уравнении Пуассона—Больцмана предполагается осуществить в дальнейшем.
ВЫВОДЫ
1. В рамках приближенной двумерной модели массопереноса вещества в микроканале чипа, основываясь на методе моментов, получены реалистичные оценки необходимой сепарационной длины канала при выбранной ширине сечения для заданных условий разделения (напряженность электрического поля и физико-химические характеристики буфера и анализируемого вещества).
2. Исследовано влияние формы симметричного конвективного профиля произвольной клиновидности, характеризуемой параметром т, на про-
цесс массопереноса вещества, что полностью подтверждает априорные качественные предположения, в частности уменьшение размывания вещества по мере роста клиновидности (параметра т).
3. Определена форма конвективного скоростного профиля, исходя из решения системы уравнений Пуассона—Больцмана, при постоянстве температуры во всем канале микрочипа и малости концентраций буфера, что позволило исключить уравнение температуропроводности из системы, а также линеаризовать и получить аналитическое решение уравнения для описания электростатического поля в канале.
4. Определен численными методами конвективный скоростной профиль и соответствующий ему коэффициент т на основе решения задачи оптимизации по критерию минимального отклонения полученного скоростного профиля от модельного u(z) = и • (1 - |^т); zе [0;1].
Благодарности
Работа выполнена при финансовой поддержке, осуществляемой в рамках научной программы Санкт-Петербургского научного центра РАН "Аналитические приборы на основе микрофлюид-ных технологий" (раздел 2, проект 2, 2003 г.).
СПИСОК ЛИТЕРАТУРЫ
1. Туницкий Н.Н., Каминский В.А., Тимашев С.Ф. Методы физико-химической кинетики. М.: Химия, 1972. 197 с.
2. Андреев В.П., Брезгун А.А., Павленко И.В. Математическое моделирование процессов мас-сопереноса в противоточной распределительной хроматографии // Журнал физической химии. 1988. Т. LXII, № 9. C. 2448-2454.
3. Рыжик И.М. Таблицы интегралов, сумм, рядов и произведений. М.—Л.: Гостехиздат, 1943. 350 с.
4. Heller S. // Electophoresis. 2000. V. 21. P. 593602.
5. Яворский Б.М., Детлаф А.А. Справочник по физике для инженеров и студентов вузов. М.: Наука, 1971. 903 с.
Институт аналитического приборостроения РАН, Санкт-Петербург
Материал поступил в редакцию 18.09.2003.
CALCULATION OF MICROSCALE SYSTEM CHANNEL PARAMETERS BY THE METHOD OF MOMENTS
A. L. Bulyanitsa, A. A. Evstrapov, G. E. Rudnitskaya
Institute for Analytical Instrumentation RAS, Saint-Petersburg
The paper presents the results of mathematical modeling of mass transfer processes in a microscale system, a microfluidic chip channel, during electrophoretic separation of substances. The method of moments is used to calculate channel lengths and analysis times required for adequate resolution of a two-component mixture. Models describing the influence of non-uniform sample distribution in the channel are discussed and analyzed. The estimates of convective flow rate profiles are given based on a set of the Poisson—Boltzmann equations under the assumption that there are no thermal field gradients in the channel.
