Расчет интенсивности свечения ЭХЛ, регистрируемой конфокальным микроскопом Текст научной статьи по специальности «Математика»

УДК 517.958: 536.71

РАСЧЕТ ИНТЕНСИВНОСТИ СВЕЧЕНИЯ ЭХЛ, РЕГИСТРИРУЕМОЙ КОНФОКАЛЬНЫМ МИКРОСКОПОМ

ОЛЕЙНИК А.И., ДРОГОВОЗОВ А.Г.,

АМАТОР K.A., СВИРЬ И.Б._________________

Рассматривается численное моделирование диффузионных электрохемилюминесцентных (ЭХЛ) процессов, протекающих в системе с двумя микрополосками-электродами. Описывается алгоритм численного расчета распределения локальной интенсивности свечения ЭХЛ, учитывающий особенности работы конфокального микроскопа.

1. Введение

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

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

Рассматривается также численное моделирование диффузионных ЭХЛ процессов в ячейке с двумя микрополосками-электродами. Расчет распределения концентрации вещества-люминофора выполняется с использованием конформного отображения, предложенного в работе [1]. Применение конформных отображений является на сегодняшний день наиболее эффективным и простым приемом (с точки зрения реализации), предназначенным для получения точных и устойчивых численных решений при моделировании двумерных диф -фузионных процессов в системах с микроэлектродами различной формы, обладающих симметрией [2]. Такое преобразование координат позволяет увеличивать плотность узлов вычислительной сетки в реальных координатах в заранее заданных участках расчетной области. В данной работе таким

РИ, 2005, № 3

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

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

В данной работе на примере теоретического расчета интенсивности свечения ЭХЛ, возникающего в системе с двумя микрополосками-электродами, рассмотрена проблема, связанная с восстановлением реального свечения, регистрируемого конфокальным микроскопом в процессе физического эксперимента. Так как конфокальный микроскоп регистрирует свечение не в точке, а в некоторой конечной области, что обусловлено его конструктивными особенностями и оптической аберрацией, то измеряемая интенсивность свечения ЭХЛ получается «размытой», т.е. является некоторым усреднением реальной величины [3, 4].

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

2. Математическая модель ЭХЛ процесса

Исследуется ЭХЛ, возникающая в электрохимической ячейке с двумя параллельными электродами в виде микрополосок, лежащими в одной плоскости, где ширина (w ) обоих электродов (анода и катода) равна между собой. Электроды разделены изолятором шириной g (рис. 1,a).

Общая схема реакций, которые приводят к возникновению ЭХЛ, имеет следующий вид [2, 5]:

Ag - e ^ A+ (анод); (1)

Ag + e ^ A_ (катод); (2)

A+ + A_ —kbl > A* + Ag (раствор); (3)

A ^Ag +уECL (раствор), (4)

23

где Ag — частица органолюминофора в основном состоянии; A+ — катион-радикал; A_ — анионрадикал; A* — синглет-частица (в электронновозбужденном состоянии); у ECL — фотон ЭХЛ эмиссии; kbi — константа скорости бимолекулярной реакции (3); kf — константа скорости мономолекулярной реакции (4).

У

а

w ; о w ?

а

ячейки, состоящей из двух микрополосок- электродов в реальных координатах (x,y) (а); область численного моделирования в преобразованных координатах (р ) (б)

Математическая модель диффузионного массопереноса в данной электрохимической системе, соответствующая уравнениям (1)-(4), имеет следующий вид [6]:

дс

dt

дс ~ ~dt

= D

= D

5 2с+ дх2

д 2с ~

дх2

д 2с+ ду2

д 2с -ду2

- кьіс+с :

- кЬіс+с ;

(5)

(6)

дс

~dt

= D

д 2с* дх2

2* д с

ду2

+ кЬіс+с - kfc ; (7)

с + с + с + cg = с0 ,

(8)

где с+ , с , с*, Cg — концентрации частиц A+ , A , A , Ag в растворе в момент времени t; D — коэффициент диффузии, принятый равным для всех частиц, участвующих в реакции; со — начальная концентрация вещества Ag в растворе.

Считаем, что длина L обоих электродов намного больше ширины каждого из электродов, т.е. l >> w , что позволяет рассматривать задачу только в двух измерениях.

Начальные и граничные условия следующие: t = о: Vx,y; cg = со , с += с_= с* = 0 ;

t > 0 : у = 0, х є Ins ;

дс

g

ду

дс+ 9с дс*

II О 1 у=0 9у II о А

у=0

= 0

(изолятор);

у = о, х.1-w-§;-§),

с += с0, Cg = с = с = 0 (анод);

у =0, х 4Ц+w

с = с0 , Cg = с+ = с* = 0 (катод);

х2 + у2 (у ^0 ),

cg ^ 1, с —^ 0, с —^ 0, с —^ 0 (бесконечно удаленная точка) ,

где Ins = |^-ro;-w -g Ju|^- g;g w + f;+»J .

3. Конформное отображение

Численное решение данной системы дифференциальных уравнений (5) - (8) с начальными и граничными условиями (9) проводилось в криволинейной системе координат, полученной с использованием модификации конформного преобразования, предложенного в [1]:

2

Q = — arcsrn ю , (10)

п

где ю = X + iY — комплексная переменная в реальной безразмерной системе координат, определяемой как

=; Y=2у

g ; g ,

(11)

Q = § + ip — комплексная переменная в преобразованной системе координат.

Конформное отображение (10) преобразует полуплоскость D(X) х D(Y) = (-го;+то) х [0;+да) в полуполосу D(|) х D(p) = [—1; 1] х [0;+да) [7] (рис.1,б).

Обратное преобразование (из преобразованных в реальные координаты) имеет вид:

ю = sin

-Q

2

(12)

а после разделения действительной и мнимой части (12) принимает следующий вид:

24

РИ, 2005, № 3

п „ Н п л § = 1, 0>0w,

2 Ч , (13)

dC+

7 Н т 7 (14)

-1 <|< 1,

dC ~ dC* dCg

=0 (изолятор);

4. Математическая модель в преобразованных координатах

Используя ранее введенные безразмерные пространственные координаты (11), а также следующие безразмерные переменные и параметры:

Cg ^1, C+ ^ 0, C ^ 0, C* ^ 0 (бесконечно удаленная точка), (21)

где Pw = — arcsin[+ (1 + W)].

c j 4Dt 2w

C =—^ • t =----• W = ——

J C0 g2 g

22 K kbic0g • K kfg

Kbi -———; Kf -■

4D

4D

(15)

і _ *

где cj принимает значения {c , c , c , cg }, и применяя преобразование координат (13),(14), перепишем исходную математическую модель (5)-(8) следующим образом:

dC+

дг

DC ~

дг

*

- = D *

= D

52C+ | 52C+ 8Е,2 + dp2

52C~ | 52C~

дЕ,2 + dp2 _

- KbiC + C- ; -KbiC+ C - ;

(16)

(17)

dC

дг

= D

д2C* д2C*

дЕ,2 dp2

+ KbiC+C-- KfC ; (18)

C+ + C~ + C + Cg = 1,

(19)

*4 здесь D = —

cos2H + sh2 p

“i—1

детерми-

нант якобиана обратного преобразования (13), (14).

5. Конфокальный микроскоп

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

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

Рис. 2. Схема работы конфокального микроскопа:

1 — исследуемый образец (электрохимическая ячейка); 2 — линзы микроскопа; 3 — диафрагма; 4 — фотоумножитель

Начальные условия (т = 0) становятся следующими:

-1 < Е, < 1, 0 <ц<да ,

Cg = 1; C+ = 0; C_ = 0; C* = 0 ; (20)

а граничные условия (т> 0) запишем как:

1 = ~1, 4>0w ,

5C+ SC" dC* ^Cg 0 ч

““““0 (изолятор);

1 = “1, 0 ^'n^'Hw ,

C+ = 1, C“ = C* = Cg = 0 (анод);

-1 <^< 1, р = 0 ,

5C+ dC~ dC* 5Cg

= 0 (изо™°р);

% = 1, 0<p<Pw,

C“= 1, C+ = C* = Cg = 0 (катод);

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

Из геометрической оптики известно, что система из двух собирающих линз (обозначенная цифрой 2 на рис. 2) фокусирует свет, исходящий изнутри конической поверхности с вершиной в фокусе линзы, опирающейся на ее центральную плоскость. Кроме того, поскольку электрохимическая ячейка имеет размеры порядка нескольких десятков микрон, необходимо учитывать явление аберрации [11, 12]. Таким образом, зоной видимости микроскопа будет объем, ограниченный однополостным гиперболоидом с центром в фокусе линзы (рис. 3). Продольный размер гиперболоида ограничен благодаря наличию диафрагмы в конструкции конфокально-

РИ, 2005, № 3

25

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

Рис. 3. Окрестность фокусной точки, воспринимаемая конфокальным микроскопом

Интенсивность свечения в точке в рассматриваемого цилиндра (на рис. 4 показано его осевое сечение, СЕ—верхнее основание цилиндра), регистрируемая конфокальным микроскопом, будет прямо пропорциональна произведению концентрации синглет-частиц на величину телесного угла с вершиной в этой точке и опирающегося на верхнее основание цилиндра. Величина телесного угла, как известно, равна площади поверхности, которую он вырезает на сфере единичного радиуса [13].

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

Рассмотрим рис. 4 для определения параметров конуса. Здесь BF — высота прямого конуса, CF — одна из полуосей эллипса, вторая полуось проходит через точку f и перпендикулярна к плоскости рисунка; X1,Y1 — локальные координаты, связанные с нормированными координатами X, Y следующими соотношениями: Xi = X - Xf и

Y1 = Yf + M/2 - Y , где (Xf, Yf) — координаты фокуса конфокального микроскопа; M = 2m/g — нормированная высота цилиндра. Параметры указанного конуса находились по следующим формулам:

a = CF = д/ (R - Xi)2 + Yi2 sin у;

b =

f cos ул/R2 - (Xi - Yi tgy)2

a

; f = BF = actgy ,

где a и b — полуоси эллипса; f — высота эллипсоидального конуса, R = 2r/g — нормированный радиус цилиндра; (Xi,Yi) — локальные координаты точки в в цилиндре (см. рис.4).

Угол при вершине конуса в сечении, показанном на рис. 4, определяется как:

а = ZCBE = arccos-j=

V((R

R2 -X,2 - Yi2_________

- Xi)2 + Yi2)((R + Xi)2 + Yi2),

а угол между высотами круглого и эллипсоидального конусов, соответственно, как:

/GBF а . R - Xi

у = ZGBF =--arcsin—р=

2 1(R-X,)2 + Y,2 •

Рис. 4. Сечение конуса с вершиной в точке B

Зная параметры прямого конуса с эллипсом в основании, определим величину телесного угла:

4 е2 2

S = 4 J бф J р 2 sin 9d9 , (22)

0 є,

здесь р — радиус сферы (р = i); 9 и ф — угловые переменные в сферических координатах; о, и 9 2 — границы изменения переменной 9 , равные:

Н = 0, 92 = arctgl -л/a

f т/a2 cos2 ф + b2 sin2 фj.

Таким образом, зная значение телесного угла в каждой точке цилиндра, безразмерную величину интенсивности ЭХЛ, воспринимаемую конфокальным микроскопом в данной точке, можем найти по формуле:

1 г *

^ECL = — J C (X, Y, Z) S(X, Y, Z) dU (23)

4n U(A) ’ V ’

26

РИ, 2005, № 3

где U(A) — рассматриваемая цилиндрическая окрестность фокусной точки A; Z = 2z/g — нормированная координата, перпендикулярная к плоскости рис. 1,а. С учетом того, что распределение концентрации одинаково во всех сечениях, перпендикулярных к оси OZ , запишем интенсивность ЭХЛ как:

1 Xf+R Yf+M/2 Z2

^ECL = — J J C*(X,Y)dXdY J S(X, Y, Z)dZ,

4л Xf-R Yf-M/2 Zi

(24)

где Z12 = ±tJR2 - (X - Xf)2 — координаты точек поверхности цилиндра в сечении, параллельном оси OZ при фиксированном X (считаем, что сканирование ведется вдоль прямой Z = 0).

*

Поскольку C является сеточной функцией, вычисляемой нами на криволинейной сетке, то для расчета интеграла по (X, Y) необходимо интерполировать значения концентрации C*. В работе по заданным координатам точки (X, Y ) находился ее образ в преобразованных координатах (§, ц) при помощи (10). Затем с использованием билинейной интерполяции на прямоугольных элементах [14] вычислялось значение концентрации в нужной точке. Необходимо отметить, что вычисления необходимо проводить для стационарного режима электролиза, чтобы поле концентраций не изменялось во время сканирования, что может повлиять на распределение интенсивности ЭХЛ. Следует отметить, что приведенный алгоритм расчета интенсивности ЭХЛ учитывает только излучение из области восприятия микроскопа и не учитывает излучение из окружающего пространства, которое также может оказывать влияние на регистрируемый сигнал.

6. Результаты вычислений и их обсуждение

На рис. 5 представлены графики вычисленных интенсивностей ЭХЛ, полученные для безразмерного времени электролиза % = 350 (это соответствует t = 8,75c для приведенных ниже значений параметров D и g), высоты фокальной плоскости над плоскостью электродов h = 5мкм и различных размеров области видимости микроскопа r и m (которые являются горизонтальным и вертикальным разрешением конфокального микроскопа).

Параметры численного моделирования:

Со = 1 млМ, D = 10_5 см2 с_1 , kf = 108 с_1, kbi = 1010 см3моль_1с_1 , w = 10 мкм , g =10 мкм , L = 0,2 см ; размер вычислительной сетки: N^x Npx Nt = 100 х 100 х 20000 .

На рис. 5 сигнал 1 получается из предположения, что конфокальный микроскоп фиксирует интенсивность только в точке фокуса; сигнал 2 рассчитывается, как средняя интенсивность ЭХЛ в видимой области микроскопа, т.е. в цилиндре [3, 4];

сигнал 3 — по алгоритму, предложенному в данной работе. Следует отметить, что зависимости, обозначенные 2 и 3 на рис. 5, нормировались на величину объема цилиндра, для того чтобы все три сравниваемые характеристики имели одинаковую размерность при переходе к размерным величинам.

в

Рис. 5. Графики зависимости интенсивности ЭХЛ (V ECL ) от пространственной координаты (X) для различных размеров области видимости конфокального микроскопа: а — m = 2Г = 10 мкм; б — m = 2r = 5 мкм; в — m = 2r = 1 мкм

Результаты вычислений показывают, что при уменьшении линейных размеров цилиндра интенсивность ЭХЛ, вычисленная усреднением по объему цилиндра, приближается, как и ожидалось, к точечной интенсивности. Интенсивность, вычисленная по предложенному в данной работе алгоритму, увеличивается с уменьшением размеров цилиндра, однако в точке максимума остается меньше значений интенсивности, полученных другими метод а-

РИ, 2005, № 3

27

ми. Это можно объяснить тем, что предложенный алгоритм учитывает свойства оптических частей микроскопа, поэтому сигнал является более “размытым”, чем те, которые были получены с помощью других подходов. Таким образом, результаты, полученные с использованием приведенного алгоритма, более напоминают регистрируемые экспериментальные данные. Сравнение результатов вычислений и экспериментальных данных будет приведено в следующих публикациях.

7. Выводы

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

Признательность. С.н.с. Олейник А.И. благодарит мэрию г. Парижа за post-doc стипендию для работы в лаборатории академика Кристиана Аматора в Ecole Normale Superieure (Париж, Франция).

Литература: 1. FossetB., Amatore C.A., Bartelt J.E., Michael A. C, Wightman R.M. The use of Conformal Maps to Model the Voltammetric Response of Collector-Generator Double-Band Electrodes. // Analytical Chemistry. 1991.V. 63.P.306-314. 2. Amatore C, in:I. Rubinstein, (Ed.), Physical Electrochemistry: Principles, Methods and Applications (Chap. 4), Marsel Dekker, New York, 1995. 3. Amatore C., Bonhomme F, Bruneel J.L., Servant L, Thouin L. Mapping concentration profiles within the diffusion layer of an electrode. Part I: Confocal resonance Raman microscopy // Electrochemistry Communications. 2000. V.2. P.235239. 4. Amatore C, Bonhomme F, Bruneel J.L., Servant L, Thouin L. Mapping dynamic concentration profiles with micrometric resolution near an active microscopic surface by confocal resonance Raman microscopy. Application to diffusion near ultramicroelectrodes: First direct evidence for a conproportionation reaction // J. Electroanalytical Chemistry. 2000. V. 484. P. 1-17. 5. Bard A.J., FaulknerL.R. Electrochemical methods: fundamentals and applications.

УДК 517.958:541.14

АНАЛИЗ НОВЫХ МЕТОДОВ ЧИСЛЕННОГО РЕШЕНИЯ ЖЕСТКИХ ЗАДАЧ КОШИ

ШУЛЫК В.Н., КЛИМЕНКО А.В., СВИРЬ И.Б.

Рассматриваются новые методы решения жестких систем обыкновенных дифференциальных уравнений — метод квази-Рунге-Кутта и метод Алуффи-Пентини. Предлагается алгоритм адаптивного изменения шага интегрирования по времени. Методы сравниваются с известным методом Г ира на примере жесткой системы ОДУ, возникающей при математическом моделировании химической реакции Бриггса-Рошера.

John Wiley & Sons: N-Y. 2002. 834 p. 6. Свирь И.Б. Формулировка обобщенной краевой задачи для моделирования электрохемилюминесцентных процессов при различных режимах нестационарного электролиза // АСУ и приборы автоматики. 2003. № 122. С. 129134. 7. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного. М: Наука, 1973. 736с.8. http://www.physics.emory.edu/~weeks/confocal/

9. http://ru.ntmdt.ru/SPM-Techniques/Basics/3_SOM/

10. http://www.microscopyu.com/articles/confocal/ confocalintrobasics.html 11. Борн М., Вольф Э. Основы оптики. М: Наука, 1973. 720с. 12. БутиковН. С. Оптика: Учебное пособие для студентов физических специальностей вузов. М.: Высш. шк. 1986. 412с. 13. Бронштейн И.Н., Семендяев К.А. Справочник по математике для инженеров и учащихся втузов. М.: Физматгиз, 1962. 608 с. 14. Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах: Т. 1: Пер. с англ. М.: Мир, 1991. 504 с.

Поступила в редколлегию 03.09.2005

Рецензент: д-р техн. наук, проф. Хаханов В.И.

Олейник Александр Игоревич, канд. техн. наук, с.н.с. лаборатории математического и компьютерного моделирования ХНУРЭ. Научные интересы: численное моделирование физико-химических процессов, уравнения математической физики. Адрес: Украина, 61166, Харьков, пр. Ленина 14. В настоящее время проходит научную стажировку (post-doctoral research) в Ecole Normale Superieure (ENS, Париж, Франция). Адрес: ENS, 24 rue Lhomond, 75231 Paris Cedex 05, France.

Дроговозов Алексей Геннадьевич, бакалавр по специальности прикладная математика, ХНУРЭ. Научные интересы: численные методы, программирование. Адрес: Украина, 61166, Харьков, пр. Ленина, 14.

Аматор Кристиан Андрэ, профессор Ecole Normale Superieure (ENS), действительный член Академии Наук Франции, декан факультета химии ENS. Научные интересы: аналитическая химия, прикладная математика, физика, физическая химия, биохимия, электроанализ, химия мозга. Address: ENS, 24 rue Lhomond, 75231 Paris Cedex 05, France.

Свирь Ирина Борисовна, д-р техн.наук, зав. лабораторией математического и компьютерного моделирования ХНУРЭ. Научные интересы: математическое моделирование в научных исследованиях: физической химии, биологии и медицине. Адрес: Украина, 61166, Харьков, пр. Ленина. 14 , е-mail: <[email protected]>.

1. Введение

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

Решение различных задач химической кинетики приводит к появлению означенных выше систем. Системы, содержащие быстро и медленно изменяющиеся компоненты, по определению являются

28

РИ, 2005, № 3