Построение оптимального квазиконформного отображения для численного моделирования диффузии в системе с дисковым микроэлектродом Текст научной статьи по специальности «Математика»

ЭЛЕКТРОНИКА

УДК517.958: 536.71

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

ПОСТРОЕНИЕ ОПТИМАЛЬНОГО КВАЗИКОНФОРМНОГО ОТОБРАЖЕНИЯ ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ДИФФУЗИИ В СИСТЕМЕ С ДИСКОВЫМ МИКРОЭЛЕКТРОДОМ

ОЛЕЙНИК А.И.___________________________

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

1. Введение

Численное решение задач математической физики, имеющих сингулярности, представляет собой большую сложность ввиду того, что скачки искомой функции или фронты быстропротекающих процессов плохо описываются на дискретных сетках. Большинство задач микроэлектродной электрохимии относится именно к такому классу задач, так как плотность тока (пропорциональная градиенту концентрации) бесконечно возрастает на краях электрода (т.е. на границах раздела фаз электрод-изолятор или углах) [1]. Это обусловлено скачкообразным изменением граничных условий и приводит к нарушению дифференцируемости решений в окрестности краев электрода. Применение аппарата теории конформных отображений позволяет разрешить эти проблемы. Как правило, искомая функция ведет себя намного «лучше» в преобразованном пространстве. Однако при помощи конформных отображений далеко не всегда удается получить ограниченную область моделирования; могут также оставаться зоны, где функция по-прежнему изменяется достаточно быстро (особенно это верно, если подобные зоны находятся на достаточном удалении от активных границ), что создает определенные сложности при численном моделировании. В таких случаях более «удачное» отображение принадлежит к классу отображений, называемых квазиконформными. В этом случае компоненты функции комплексного переменного, осуществляющей искомое преобразование координат, являются не гармоническими функциями, а решением эллиптической системы первого порядка [2].

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

2. Идея

2.1. Вывод оптимального преобразования координат

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

(т = Dt/r^; R = r/rd; Z = z/rd, где rd- радиус диска; D - коэффициент диффузии):

5C д 2C 1 д C д 2C „

---=----Т +-----+----;т = 0 (1)

дг дR2 R 3R qZ2 , (1)

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

0 < R < 1, Z = 0, C = 0;

R > 1, Z = 0,

д C д Z

= 0.

R = 0 , Z > 0 ,

д C 5 R

= 0.

R 2 + Z2 —^ & ,

C ^ 1. (2)

Произведем конформное преобразование исходной области, т.е. первого квадранта (см.рис.1,а), на бесконечную полуполосу (рис.1 ,б) - применим преобразование вида [3]:

X = sin п R = sin п £* _ 2 ^ _ ch п * —р L 2 J

—ю _ 2 _ Z = cos п г* J ^ _ sh п * —р _ 2 _

(3)

где X = R + iZ, ю = |* + ip*. Уравнение (1) примет вид:

5C , т*

— = det J

дх

д2C | 52C

*2 *2 ор

+

+

1

* *

R(5 , л )

dC dp* д£ cR + ар* 3R

= 0,

(4)

*

здесь detJ - детерминант Якобиана преобразования координат (3).

РИ, 2006, № 1

21

Граничные условия запишутся как:

VI * Н_Р VI о о II * О II о

іГїї * II 0 <р <да — = 0. д\ ;

іГїї * II о * 0 <р <да — = 0. ;

о Л іГїї * Л р ^ да C ^ 1. (5)

Отметим здесь тот факт, что в пространстве (£*, р*)

решение задачи уже не зависит от координаты |* (доказательство будет приведено ниже), т.е. линии уровня концентрации параллельны оси .

электрод 1

а

изолятор ~R

О.

О

электрод 1

б

бесконечно 1 удаленная точка

s s

CL I-Ч>

s г

S

о л о о

г

CL

р

электрод

в

Рис. 1. Область моделирования: а - в реальном пространстве; б - в пространстве, полученном при помощи конформного отображения (3); в - в оптимальном пространстве для численного моделирования.

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

висимость концентрации вдоль оси р* (от координаты | концентрация не зависит, см. выше). Очевидно также, что наиболее удобной и легко поддающейся описанию в терминах численного решения будет линейная зависимость концентрации от координаты р*. Такое пространство, в котором распределение концентрации не зависит от абсциссы и является линейной функцией ординаты, т.е. C(§,р) = р , мы будем в дальнейшем называть оптимальным пространством для численного моделирования.

Из приведенных рассуждений становится ясно, что построение такого пространства эквивалентно аналитическому решению задачи (1)-(2), которое будет получено при подстановке C(|, р) = р в преобразование координат, задающее соответствие между оптимальным пространством и областью в реальных координатах. Покажем, что для рассматриваемой задачи такое пространство может быть получено сравнительно просто.

В силу того, что концентрация не зависит от абсциссы, для получения оптимального пространства для численного моделирования необходимо найти лишь функцию р = F(p*), задающую необходимое отображение [0, да) ^ [0,1]. Предположим, что такая функция сжатия р = F(r|*) нам известна, и перепишем уравнение (4) в идеальном пространстве (|, р), произведя замену переменных £ = £*, р = F(p*):

дС

дг

= det J

д2С | д2C Г dp 12 ^2 + Зр2 I dp*,

dC d2p + Зр dp*2

+

1

+---------і---

R& F_1[p])

дСд^ dC dp 5p* 5p d^*"SR

= 0

(6)

В граничных условиях (4) изменится лишь предел ординаты: 0 < p < 1.

Так как уравнение (6) по нашему предположению записано в оптимальных координатах, то C(E,, р) =р является решением задачи. Подставив C(E,,р) = р в уравнение (6), получим следующее уравнение относительно функции сжатия (р = F(p*)):

2 * d р 1 5р dp

dp R(|, р ) det J dR dp

(7)

Выписывая явно детерминант отображения ар*/3R и R(|, р*) для преобразования (3) [3,4], получаем:

dp

d2p л

— th

dp

*2

_ dp

■ = 0

(8)

2

2

с граничными условиями q* = 0, р = 0; р* ^ да, р ^ 1. Коэффициенты уравнения (8) не зависят от £, , следовательно, предположение о том, что для получения оптимального пространства достаточно произвести преобразование координаты q*, является верным. Следовательно, это является доказательством

РИ, 2006, № 1

22

того, что распределение концентрации после преобра- ство для численного моделирования, без какой-либо

зования (3) не зависит от £*. информации, известной априори.

Порядок уравнения (5) может быть понижен путем следующей замены переменных: p = Мц |5J.

Таким образом, в итоге имеем:

Как было сказано выше, вывод оптимального пространства эквивалентен аналитическому решению задачи (1)-(2), и это решение имеет вид:

dp л , —^ + - th d^* 2

л * 1_ 2 Ъ II О II c1 C(R,Z) = f[ Im ' 2 " —arcsin(R + iZ) л

> р - п *

(13)

где функция F(d*) задана уравнением (9). Решение

2

Произведя обратную подстановку, получим:

d-p _ d^* ch

c1

2

ц = c^arctg

shl fV

+ c2 ,

2

(13) эквивалентно решению, полученному впервые Сайто [6] и приведенному в более удобную форму в работе [7].

2.2. Обобщение подхода с использованием сплайн-функций

где ci и c2 - константы интегрирования.

Поэтому, с учетом граничных условий c2 = 0 и ci = 1, искомая нами функциональная зависимость имеет свой окончательный вид как:

* 2

^ = F(^ ) = —arctg л

shl з2

(9)

При подстановке функции, обратной функции (9), а именно

*

ц

—arcsh л

tg| -2л

(10)

в уравнения (3), получаем новое преобразование координат, описывающее переход из оптимального пространства в реальное:

R = sin

sec

Z = cos

tg

(11)

2

2

Если отступить от рассмотренного примера диффузии в системе с микроэлектродом в форме диска и рассмотреть другие геометрические формы электродов, то очевидно, что коэффициенты в уравнении вида (7) не всегда будут такими же простыми для преобразований, как это получилось в (8). На практике достаточно часто искомая функция после перехода к пространству (|*, ц*) не является функцией, зависящей только от координаты ^ . Это значит, что присутствует слабая зависимость от координаты |*, когда линии уровня близки к прямым, параллельным оси абсцисс, но полностью не совпадают с ними. В таком случае изложенный выше подход может быть применен для численного построения оптимального пространства. И хотя «аналитическое» решение, полученное таким образом, будет представлять собой лишь аппроксимацию истинного решения, тем не менее, построение оптимального пространства, даже численно, представляет собой интерес как с точки зрения решения стационарной задачи (1)-(2), так и ее времязависимо-го аналога.

т.е. то, что было получено впервые в работе [4] на основе известного аналитического решения задачи (1)-(2).

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

5R (л ЛдZ 5Z (л ЛдR

secli ^~ se\2 ^. (12)

Таким образом, было получено квазиконформное отображение, определяющее оптимальное простран-

Далее будем исходить из предположения, что конформное отображение области моделирования в реальных координатах на полуполосу (т.е. преобразование (3) для случая микродиска) уже известно. Заметим, что принципиальных ограничений это предположение не несет, так как для односвязных областей такое отображение всегда может быть найдено в силу основной теоремы теории конформных отображений [2].

В уравнении (6) под функцией сжатия будем подразумевать теперь не некоторую функцию F(^*), удовлетворяющую какому-то дифференциальному уравнению, а вполне определенную функцию, задающую отображение [0, да) ^ [0,1], как, например, ц=г| /(1 + ц ) или функцию (9). Если решение теперь не зависит от координаты £ (или зависит очень слабо), то, полагая С(£,ц) = С(ц) (соответственно С(£,ц)« С(ц); в этом случае необходимо зафиксировать некоторое значение £ для вычисления коэффициентов в уравнении (11)), из (6) получаем:

РИ, 2006, № 1

23

if \2

d2CI dp A

3p2

dp

+

2 * d p 1 dp Эр

dp*2 R(|, F_1[p]) det J* dr|* 3R

5C °

aT° (14)

Обозначим для упрощения записи через f функцию, обратную к функции F, т.е. р* = f(p) = F_1(p). В терминах функции f и после некоторых упрощений уравнение (14) запишется как:

32C

Эр2

+

f''(p) _ f'(F) & дС = °

f'(r|) R(|, f(p)) det J* 5Z J 3p , (15)

где коэффициент при первой производной может быть вычислен точно в любой точке. Граничные условия будут иметь вид: р = 0, C = 0; ц = 1, C = 1.

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

Задавая разбиение с постоянным шагом Др и подставляя выражение для кубического сплайна [8] в уравнение (12), а также требуя непрерывности второй производной в узлах разбиения, получаем систему линейных уравнений относительно значений и наклонов сплайна в точках разбиения. Полученную систему линейных уравнений решим в следующем виде:

C(p) = popt(p) = £ S3(p), (16)

i=°

где popt - координата, вдоль которой концентрация изменяется линейно; S3 (р) - кубический сплайн, заданный на i -м отрезке (рц Pi+1) и определяемый двумя значениями функции на концах отрезка (Pi,pOpt), СПі+цЛоРі) и значениями прозводных в этих точках mi и mi+i:

2

Si (Л) = СЛї+1 ~Л) (2(Р~РО + Ар) ^opt +

3 Др3 i

, (F-Fi) (2(Pi+l -F) + A^) opt ,

+-----------7-3------------Vl +

Др3

(Pi+1 -Ц) CF-Fi) , (F-Fi) (p-Pi+0

Др2

Др2

mi+i .(17)

Оптимальная функция сжатия относительно исходной координаты р* будет иметь вид:

р =■

p(p°pt)

1 -p(p°PV

(18)

если для первоначального сжатия использовалась функция р = р*/(1 + р*), и

р = — arcsh %

ig[ f ^(^opt)

(19)

если применялась функция (9).

Для получения этого отображения нам необходима функция, обратная построеному сплайну (16). Такая функция может быть легко получена, так как обращение каждого i -го сплайна будет задаваться тем же самым выражением, но с параметрами (p°pt,Pi), ('П1°+1,Pi+l), 1/mi, 1/mi+1, т.е.:

n . _1

p(popt) = xS3 (Л),

i=°

2

i -1 (F°p1 -F)2(2(p-p°pt) + Ap)

S3 (p) = -A+1------------i-------ni

3 (A.°P,>3

On - F°Pt )2(2(p°pl1 - P) + Ap)

(Ap°pt)3 ^

(rffi-p)2(p-P°Pt) 1

+ ——-------------------+

Pi +

(Ap°pt)2

mi

(л-лоР,)2(л-л°РІ) 1

(Ap°Pt)2

mi+1

(20)

где Ap°pt =p°+pi -P°pt.

Далее, суперпозиция функции (20) с преобразованием (3) является оптимальным квазиконформным отображением, так как концентрация является линейной функцией координаты popt. Математическая модель в оптимальном пространстве для численного моделирования (|, popt) перепишется следующим образом:

3C d ,т*

— = det J Зх

3 2C

д\ 2

1

3 2C

[p'(popt)]2 ар opt2

1

3C д%' | I p''(popt) +_

1

R& popt) 3p I 3C

3| SR

[p'(popt)]3 p'(^L)

opt) 3R

3p

opt

= °

?

(21)

с соответствующими граничными условиями:

° < ^ < 1 popt = ° C = °;

24

РИ, 2°°6, № 1

1 = 0 1=1 0 <l< 1

0 <popt < 1 0 < popt < 1

Popt = 1

d C

= 0;

d C

= 0;

d%

C =1. (22)

Все коэффициенты в уравнении (21) могут быть вычислены аналитически при помощи выражений (3) и (20).

3. Результаты

Для демонстрации работы метода, изложенного в пункте 2.2, рассмотрим времязависимую диффузионную задачу в системе с дисковым микроэлектродом (т.е. задача (1)-(2), с ненулевой производной по времени в уравнении (1) и начальным условием: R > 0, Z > 0, C = 1). Используя конформное преобразование (3), предположим, что оптимальная функция сжатия (10) нам неизвестна и в качестве первоначальной функции сжатия выберем р* = р /(1 - р). Коэффициент в квадратных скобках в уравнении (15)

(обозначим его g(p)) для данного частного случая имеет вид:

g(n) = -

2

1

1 -ц

th

Ц

2 1 -р

(23)

2 (1 -ц)

Решая уравнение (15) с коэффициентом (23) при помощи сплайнов, получаем зависимость popt =popt(p) в виде сплайна (16). Данная зависимость изображена на рис. 2, а. Обратив построенный сплайн, используя выражения (20), получим оптимальную функцию сжатия (рис.2,б), где сплошной линией показано полученное решение в виде сплайнфункции, а символы - соответствуют оптимальной функции сжатия (9), построенной в пространстве р = р*/(1 + р*), т.е. функции:

popt = F(p) = — arctg л

2

sh

л р

2 1 -р

(24)

Рис. 2. Оптимальная функция сжатия для системы с микродисковым электродом: а - сплошная линия -сплайн-функция (16), а символы - аналитическая зависимость (23); б - сплайн-функция (20), т.е. функция, обратная к сплайн-функции (16)

На рис. 3 показано численное решение задачи (21)-(22), полученное в оптимальном пространстве (|, р opt) методом переменных направлений [9] для разных моментов времени. Здесь показаны только зависимости от координаты popt, так как от координаты б, решение не зависит. Согласно нашим ожиданиям профиль концентрации стремится к линейному с увеличением времени и является строго линейной зависимостью в предельном случае, когда t ^ <ю.

а

Как видно из рис. 2, численно построенная функция сжатия очень хорошо согласуется с аналитически полученной зависимостью (24). Очевидно также, что

данная зависимость будет описывать профиль кон**

центрации в пространстве р = р /(1 + р ) в стационарном режиме. Следует также отметить, что рис. 2 содержит качественную информацию о поведении первоначальной и оптимальной функции сжатия. Так, из рис.2 видно, что поведение первоначальной функции сжатия р = р*/(1 + р*) в приэлектродной области (т.е. в окрестности точки р = 0) подобно поведению оптимальной функции сжатия, но оптимальная функция намного быстрее стремится (при р* ^ ж) к их общему асимптотическому значению, равному единице.

РИ, 2006, № 1

Рис. 3. Профили концентрации в системе с микродисковым электродом, расчитанные в оптимальном пространстве для численного моделирования (^, popt) , соответствующие различным моментам времени: 1 -т = 0.04 ; 2 - т = 0.2; 3 - т = 1; 4 - т = 16; 5 -

25

4. Выводы

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

Признательность. Автор благодарит мэрию Парижа за годичную стипендию (post-doctoral fellowship) в институте Ecole Normale Superieure в лаборатории академика Французской Академии наук Кристиана Аматора.

Литература: 1. Bard A.J., Faulkner L.R. Electrochemical methods: fundamentals and applications. John Wiley & Sons: N-Y. 2002. 834p. 2. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного. М:

Наука, 1973.736с. 3. Michael A.C., WightmanRM., Amatore C.A. Microdisk Electrodes. Part I. Digital Simulation with a Conformal Map. // J. Electroanal. Chem., 1989.V. 267. 33 р. 4. Oleinick A., Amatore C., Svir I. Efficient quasi-conformal map for simulation of diffusion at disk electrodes. // Electrochem. Commun., 2004. V.6. 588 р. 5. Эльсгольц Л.Э. Дифференциальные уравнения и вариационное исчисление. М.: Наука, 1969. 424 с. 6. Saito Y. Rev. polarogr. (Jpn), 1968. V. 15. P.177. 7. CrankJ., FurzelandRM. The treatment of boundary singularities in axially symmetric problems containing discs. J. Inst. Maths Applic., V. 20, 1977. 355 Р. 8. Волков Е.А. Численные методы. М: Наука, 1987. 248с. 9. Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах: Т.1: Пер. с англ. М.: Мир, 1991. 504 с.

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

Рецензент: д-р техн. наук, проф. Свирь И.Б.

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

26

РИ, 2006, № 1