Численное моделирование диффузионных процессов в системе с микроэлектродом в форме диска Текст научной статьи по специальности «Математика»

УДК 517.958:536.71

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИФФУЗИОННЫХ ПРОЦЕССОВ В СИСТЕМЕ С МИКРОЭЛЕКТРОДОМ В ФОРМЕ ДИСКА

МАЗЕПА А.С., ОЛЕЙНИКА.И., СВИРЬ И.Б.

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

1. Введение

Конформные отображения, часто применяемые в настоящее время для численного моделирования двухмерных электрохимических задач, являются самым успешным и универсально простым, с точки зрения программирования, приемом, применяемым математиками для получения точных и стабильных численных решений в системах с микро -электродами разнообразных геометрических форм. Для численного моделирования диффузионных процессов в таких электрохимических системах использовалось преобразование Шварца-Кристоф-феля, квазиконформные преобразования, а также преобразования, которые получаются посредством суперпозиции элементарных функций комплексных переменных [1, 2]. Так, конформные преобразования успешно использовались для моделирования диффузионных процессов на микродиске [3], микрокольце [4], для микрополоски-электрода и системы из двух микрополосок-электродов [5-8], микросферы [8,9] и цилиндрических/полуцилин-дрических микроэлектродов [6, 9-12].

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

Цель данного исследования — дальнейшее развитие численного моделирования диффузионной кинетики в системе с микродисковым электродом с использованием конформного преобразования, предложенного в работе [3], и построение неравномерной сетки по времени, описанной в работе [13] для долгих времен электролиза. Такое сочетание двух подходов используется впервые и имеет ог-

РИ, 2004, № 1

ромное значение для эффективности вычислений как для коротких, так и для долгих времен электролиза. Данный подход мы описываем здесь на примере хроноамперометрии.

2. Математическая модель

2.1. Хроноамперометрия

Рассмотрим реакцию простого переноса электрона, которая протекает на микроэлектроде в форме диска:

A + e" ^ B . (1)

Уравнение диффузионного массопереноса вещества B в цилиндрических координатах имеет вид:

дс

dt

( * 2

= D

2 ^

д с 1 дс д с

v5r2 г 5r &2 j

(2)

где с — концентрация вещества В; D — коэффициент диффузии веществ А и В; (r, z) — пространственные координаты; t — время.

Начальные и граничные условия определяются следующим образом:

t = 0: Vr,z; с = 0;

t > 0: 0 < r < , z = 0, с = с0;

r > rd,

r = 0, r ^да,

z = 0,

0 < z < да, z ^ да,

дс

~dz

дс

Sr

= 0;

z=0

= 0;

r=0 с = 0,

(3)

здесь rd — радиус микродиска; с0 — начальная концентрация вещества А в растворе-электролите.

Произведем следующую замену переменных, чтобы получить безразмерные величины:

r = -L. z=^. с = -і.. iD

rd

rd

с0

r

Нормированное уравнение диффузионного переноса имеет следующий вид:

(4)

массо-

ас _ 1( д2с 1 дс д2с" аг 4 [ж2 R aR + az2 y

с начальными и граничными условиями:

(5)

х = 0: VR,Z;

т> 0: 0 < R <1, N II о

R > 1, N II о

R = 0, 0 < Z < да

R —^ да, Z —^ да,

с = 0; с = 1;

ас

az

ас

SR

= 0 ;

Z=0

= 0 ;

R=0

с = 0.

(6)

35

Применим конформное преобразование, предложенное Аматоре и Фоссет [3] для численного моделирования диффузионного массопереноса на микродиске, в целях получения наиболее сгущенной пространственной сетки на границах раздела фаз электрод-изолятор, т.е. там, где плотность тока максимальна. Данное конформное преобразование имеет следующий вид:

R

cosff Г

; z = 0 ,g (f Г

(7)

Уравнение диффузионного массопереноса (5) в трансформированных координатах запишем так:

ас

ат

1

4

і

-----;-------X

02 + tg2

|г

X

4 2

—cos

п

2 г

а 2с аг2

_д_

С0

(1 - 02)

ас

С0

(8)

с начальными и граничными условиями: т = 0: V0,r; C = 0;

т>0: 0<©<1, Г = 0 , с = 1;

ас

© = 0 , 0 <Г< 1, 50

ас

© = 1, 0 <Г< 1, 50

= 0.

©=0

= 0 .

©=1

0 <©< 1, Г = 1, с = 0. (9)

2.2. Ток электролиза

Выражение для нахождения тока на диске в реальных нормированных координатах имеет вид:

I(t) = 2radnFDc0 J

0

ас

az

RdR

Z=0

(10)

Запишем безразмерный ток электролиза как I

_ _ % j. ас

^ 4nFDc0Td 2 0 aZ

RdR

Z=0

(11)

Тогда безразмерный ток (11) в конформном пространстве имеет вид:

1

V=l

0

ас

аг

d©

Г=0

3. Численное решение

3.1. Гопскотч

(12)

В процессе численного решения поставленной задачи использовался метод Г опскотч [ 14]. Он состоит из двух последовательных шагов. Первый заключается в вычислении концентрации с?^1 в

узлах сетки, где (i + j + n) — четное число, по следующей явной схеме:

36

/-чП +1 Л

Ci,j 4

( (сП _ 2сП I сП ^

Ci,j-1 2Ci,j + Ci,j+1

X

V

АГ

2

■ X 2

+ ^3

С сП _ 2сП I сП ^

Ci-1,j 2Ci,j+ Ci+1,j

A©

2

сП _ сП Ci+1,j Ci-1,j

2 A©

У

JJ

Cinj

(13)

где i, j — ицдексы узлов расчетной сетки по осям © и Г, соответственно; n — индекс шага по времени;

h = Jycos2 (7 Г), ^ 2 = 1 -®2)

,2( л

f

Х3 = -2© , X4 = Ay 4^©2 + tg2^2rjj .

На втором шаге концентрация сп j-1 вычисляется в оставшихся узлах сетки (в тех узлах, где (i + j + n ) — нечетное число) по следующей схеме:

сП+1 _

Ci,j “

f

1/1-

2^1^ 4 2Х 2 X 4

V/

2

f f

X X1

V V

АГ 2

АГ

Л

+ X 2

А©2 j

/-’П , Г'П (ПП , Г'П ^

Ci,j-1 + Ci,j+1 , ^ Ci-1,j + Ci+1,j

У

-Х3

(сп + сп

C1+1,J+Ci-1,J

2 A©

УУ

A©

Л

■Cj^j

2

(14)

3.2. Равномерная сетка по времени

Критерий стабильности вычислений для Г опскотч

в цилиндрических координатах имеет вид:

2

Ат <

AZ 2

(15)

Запишем данный критерий стабильности в конформных координатах (©, Г), используя (7):

2

(

Ат< -d© +—dr| = 2 la© аг

2

tg I — r|d©-

л&

2cos21 -Г

2

dr

(16)

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

(

Nt>

T

1e

Дтп

■ = 2Te

Л

-2

і л Л1Д лД0ДГ

tg|-Ark© +

2 ) „ ,2 к

.(17) РИ, 2004, № 1

2 cos21 — ДГ

где Te — безразмерная длительность электролиза; Атр — шаг по времени, обеспечивающий стабильность вычислений, полученный из (16) при © = А©, Г = АГ , так как в этой точке значение шага по пространству в реальных координатах минимально.

3.3. Неравномерная сетка по времени

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

Ахк = (1 + ю)Лі£_!; ш > 0; k > 0 , (18)

здесь Aik — шаг по времени на k-й итерации; ю — параметр сетки, регулирующий скорость роста шага по времени.

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

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

Исходя из закона Нернста-Энштейна, запишем:

Ax<AZ2 . (19)

Введем ограничение на ю такое, чтобы для любого ю < rajim выполнялось условие (19). Представим (18) в следующем виде:

Aik = (1 + ю)к Ахр ; ю> 0 ; к > 0. (20) Определим максимальный шаг по пространству (AZmax) как шаг, соответствующий границе диффузионного слоя (5.,/ТЄ). Тогда максимальный шаг по времени (Axmax) согласно (19) можно записать так:

Таблица 1

Значения rajim для разных длительностей электролиза

4. Результаты и обсуждение

4.1. Равномерная сетка по времени

Рассмотрим вычисленный безразмерный ток (хроноамперометрия) и сравним его с результатами, полученными численно в [3, 13], а также их сходимости (% ошибки) относительно аналитических аппроксимаций для коротких и долгих времен электролиза, полученных в работе [15].

Табл. 2 содержит сравнительный анализ данных результатов, полученных на одинаковых сетках по пространству (N©x NT = 120 х 120 ) и времени (Nx =10000). Во всех подходах использовался Гоп-скотч. В работе [13] применялся новый численный метод, основанный на перекрывающихся доменах. Рис. 1 демонстрирует кривые ошибки, полученные для каждого из приведенных в табл. 1 подходов.

Те ®lim

0,01 0,140

0,04 0,101

0,09 0,100

0,25 0,119

0,64 0,165

1,21 0,220

2,25 0,302

4 0,411

6,76 0,544

25 1,020

100 1,538

Рис. 1. Сравнение нескольких численных подходов: 1 — метод двух перекрывающих доменов [13];

2 — результаты данной работы

Ах

max

= AZ

2 _ max -

tg| 2 ^max

A©-

n АГ©

max

2cos21 r,

max

Из уравнения (20) получаем равенство:

г0

Т _ ^Tmax

Те _

®lim

(22)

Отсюда следует, что определяется следующим образом:

„ _ ^Tmax _^т0

®lim - Т . (23)

Te

Все значения для используемых безразмерных времен электролиза представлены в табл. 1, для сетки 120 х 120 по пространству.

2

Из рис.1 видно, что применяемый нами подход с использованием конформного отображения имеет лучшую сходимость для промежуточных времен (в районе точки х = 1,21), чем доменный метод [13],

Таблица 2

Сравнение величины безразмерного тока (V ) и ошибки решений, полученных различными численными подходами

T Le Аналитическое решение [151 Численные результаты

из работы [13] данной работы

Ошибка, % Ошибка, %

0,01 9,657 9,6533 0,04 9,6414 0,16

0,04 5,235 5,2356 0,01 5,2364 0,03

0,09 3,768 3,7683 0,01 3,7704 0,06

0,25 2,605 2,6048 0,01 2,6071 0,08

0,64 1,968 1,9660 0,1 1,9669 0,05

1,21 1,692 1,6777 0,84 1,6859 0,3

2,25 1,495 1,4878 0,48 1,4934 0,1

4 1,366 1,3674 0,11 1,3656 0,03

6,76 1,279 1,2794 0,03 1,2793 0,02

РИ, 2004, № 1

37

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

Для автоматического определения в программе размеров требуемой равномерной сетки по времени и избежания осцилляций в начальные моменты электролиза использовался критерий устойчивости Гопскотч (см. подраздел 3.2).

4.2. Неравномерная сетка по времени

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

поэтому как оптимальное значение для дальнейших расчетов мы выбирали га = 10_3 , как лучший параметр расширения сетки для всех исследуемых

Рис. 2. Сравнение параметров расширения неравно

мерной сетки по времени: 1 — ю= 5 -10 3 ; 2 — -\—3. Л

-4

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

ю = 10 3 ; 3 — ю = 5 х 10

Табл. 4 показывает сравнение результатов, полученных разными численными подходами для долгих времен электролиза с применением неравномерной сетки по времени (га = 10_3). Как видно из данной таблицы, наши результаты хорошо согласуются с аналитическим решением и результатами других исследователей [13,16-18], получивших свои результаты с применением иных численных подходов. Во всех работах применялся одинаковый размер пространственной сетки, а именно 120 х 120.

Из рис. 3 можно судить об экономичности применения неравномерной сетки по времени в сравнении с равномерной сеткой. Количество шагов по времени, требуемое для вычислений с применением неравномерной сетки по времени, значи-Таблица 3 тельно меньше, чем с применением равномерной сетки при практически одинаковой

Сравнение безразмерн°г° тоюі (V ) и сходимости дая ошибке сходимости (см. табл. 2 и 3).

различных сеток по времени __

Табл. 5 содержит результаты вычисленных безразмерных токов и сравнительную ошибку относительно аналитического решения [15] для исследуемых коротких и промежуточных времен электролиза с применением двух численных подходов: конформного отображения (данная работа) и метода двух перекрывающихся доменов (данные из работы [13]) с неравномерными сетками по времени. В доменном методе использовался параметр расширения сетки (га = 5 х 10 _ ) [13] как оптимальный, в рассматриваемом случае оптимальное значение данного параметра га = 10_ . Результаты вычисленного безразмерного тока были получены для времени Te = 6,76 в после-

T Le Параметр расширения неравномерной сетки по времени

а= 5 х 10_3 а= 10_3 а = 5 х 10 “4

Время, мин:с Ошибка, % Время, мин: с Ошибка, % Время, мин: с Ошибка, %

0,01 0:12 0,25 0:49 0,17 1:25 0,16

0,04 0:14 0,23 1:0 0,02 1:46 0,02

0,09 0:16 0,39 1:07 0,04 2:00 0,06

0,25 0:17 0,77 1:14 0,05 2:18 0,07

0,64 0:19 1,39 1:22 0,11 2:35 0,07

1,21 0:20 2,00 1:27 0,36 2:45 0,31

2,25 0:21 2,11 1:32 0,18 2:54 0,12

4 0:22 2,21 1:37 0,11 3:05 0,05

6,76 0:23 2,23 1:39 0,08 3:10 0,00

25 0:25 2,25 1:52 0,25 3:31 0,19

100 0:28 1,65 2:03 0,08 4:01 0,03

Таблица 4

дней точке на графике кри-

38

Сравнение различных подходов вычисления безразмерного тока (у ) для вой тока в программе. Затем длинных времен и их сходимости относительно аналитических решений путем движения маркера в

сторону уменьшения времени были сняты все «показания» оставшихся временных точек вдоль вычисленной кривой тока и рассчитаны значения относительной РИ, 2004, № 1

T Le Аналитическое решение, V [16] Численное решение, у

из работы [17] из работы [18] из работы [13] данная работа

25 1,146 1,141 1,146 1,146 1,143

100 1,072 1,071 1,075 1,072 1,071

погрешности вычислений (относительно аналитических решений для указанных времен) (см. табл.5). Размер пространственных сеток в обоих случаях был 120 х 120. Рис. 4 наглядно показывает, для каких времен электролиза лучшим является доменный метод, а для каких — метод конформного отображения.

Рис. 3. Количество шагов по времени (Nx ), необходимое для равномерной (сплошная линия) и неравномерной сетки (пунктирная линия) при различных значенях безразмерного времени электролиза

0.01 0.1 1 X

Рис. 4. Сравнение численных подходов на неравномерной сетке по времени: 1 — метод двух перекрывающих [13]; 2 — результаты данной работы

Все численные результаты были получены с использованием ПЭВМ Celeron400 МГц, оперативная память — 64 Мб.

5. Выводы

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

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

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

Литература: 1. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного. М.: Наука, 1973. 736 с. 2. Amatore C. Physical Electrochemistry: Principles, Methods and Applications (I. Rubinstein, Ed.), M. Dekker, New York. 1995. Chap. 4. 3. Amatore C.A, Fosset B. // J. Electroanalytical Chemistry. 1992. Vol. 328. P. 21. 4. Amatore C.A, OleinickAI., Svir I.B. // J. Electroanalytical Chemistry. 2004. Vol. 563. P.1. 5. Deakin M.R., Wightman R.M., Amatore C.A. // J. Electroanalytical Chemistry. 1986. Vol. 215. P. 49. 6. Amatore CA, Fosset B, Deakin M.R., Wightman R.M. // J. Electroanalytical Chemistry. 1987. Vol. 225. P. 33. 7. Amatore CA, Fosset B, Maness K.M., Wightman R.M. // Analytical Chemistry. 1993. Vol. 65. P.2311. 8. Fosset B, Amatore CA., Bartlet J.E., Michael AC., Wightman R.M. // Analytical Chemistry. 1991. Vol. 63. P.306. 9. C. Amatore, B. Fosset // Analytical Chemistry. 1996. Vol. 68. P. 4377. 10. Amatore CA, Deakin M.R, Wightman R.M. // J. Electroanalytical Chemistry. 1986. Vol. 207. P.23. 11. Amatore CA, Oleinick AI, SvirI.B. // J. Electroanalytical Chemistry. 2003. Vol. 553. P.49. 12. Amatore C.A., Oleinick A.I., Svir I.B. // J. Electrochemistry Communications, 2003. Vol. 5. P.989. 13. Amatore C, Svir I. / / J. Electroanalytical Chemistry. 2003. Vol.557. P.75. 14. Рихт-майер Р, Мортон К. Разностные методы решения краевых задач. М.: Мир, 1972. 418 с. 15. K Aoki and J. Osteryoung// J. Electroanal. Chem. 1984. Vol. 160. P.335. 16. Y. Shoup, A Szabo, // J. Electroanal. Chem. 1982. Vol. 140. P. 237. 17. J Heinze // J. Electroanal. Chem. 1981. Vol. 124. P. 73. 18. Gavaghan D.J., Rollett J.S. // J. Electroanal. Chem. 1990. Vol. 295. P. 1.

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

Рецензент: д-р техн. наук, проф. Стоян Ю.Г.

Мазепа Андрей Сергеевич, студент 4 курса ПММ, техник лаборатории математического и компьютерного моделирования кафедры БМЭ ХНУРЭ. Научные интересы: программирование. Адрес: Украина, 61166, Харьков, пр. Ленина, 14. тел. 702-13-64.

Олейник Александр Игоревич, аспирант, науч. сотр. лаборатории математического и компьютерного моделирования кафедры БМЭ ХНУРЭ. Научные интересы: математическое моделирование диффузионных процессов. Адрес: Украина, 61166, Харьков, пр. Ленина, 14. тел. 702-13-64.

Свирь Ирина Борисовна, д-р техн. наук, гл. науч. сотр., зав. лабораторией математического и компьютерного моделирования кафедры БМЭ ХНУРЭ. Научные интересы: численное моделирование процессов массопереноса. Адрес: Украина, 61166, Харьков, пр. Ленина, 14. тел. 702-13-64.

Таблица 5

Сравнение двух численных подходов: конформное отображение (данная работа) и метод двух перекрывающихся доменов (данные работы [13]) с неравномерными сетками по времени,

когда Te = 6,76 , а х < Te

т Аналитическое решение [151 Численные результаты

из работы [13] из данной работы

У У Ошибка, % У Ошибка, %

0,01 9,657 9,6533 0,04 9,6405 0,17

0,04 5,235 5,2356 0,01 5,2338 0,02

0,09 3,768 3,7681 0,00 3,7684 0,01

0,25 2,605 2,602 0,11 2,6055 0,02

0,64 1,968 1,9612 0,34 1,9657 0,12

1,21 1,692 1,6814 0,63 1,6847 0,43

2,25 1,495 1,4911 0,26 1,4921 0,19

4 1,366 1,3655 0,04 1,3645 0,11

6,76 1,279 1,279 0,00 1,2781 0,07

РИ, 2004, № 1

39