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

ЭЛЕКТРОНИКА

УДК 541.138:535.379

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

СВИРЬ И.Б.______________________________

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

1. Введение

Основной проблемой, возникающей в процессе численного решения задач диффузионной кинетики на микродисковом электроде, являются краевые условия, которые обусловливают большое число ошибок в вычисленных концентрациях, а следовательно, и в значениях токов. Недавно Гавахан [13] использовал чисто неявную схему с модификациями граничных условий Рихтмайера и метод переменных направлений с экспоненциально-рас-ширяющейся сеткой для численного решения диффузионной стационарной задачи хроноамперомет-рии и линейно-изменяющейся вольтамметрии. Несколько численных подходов были предложены для решения разнообразных задач на микродиске. Так, Кранк и Фузеленд [4], Гавахан и Роллет [5], Галцеран и др. [6] использовали разложение Тейлора для получения более точных решений на границах [7], а Мишель, Вайтмен и Аматоре с соавторами [8-10], Вербрюже и Бейкер [11] использовали конформные отображения для эффективного решения на границах микродиска в двумерных задачах. Эванс и Гоурлей [12] использовали локальные разложения в ряды для вычисления потоков, что, однако, не дало корректных решений концентрации. Существует также большое число аналитических аппроксимаций для подобных задач (например, [13-20]).

В данной статье предлагается применение метода переменных направлений с преобразованием пространственных координат, предложенных Аматоре

и Фоссет [10], и неравномерной сетки по временной координате для решения нестационарной диффузионной ЭХЛ задачи на микродиске.

2. Модель

Механизм возникновения свечения — ЭХЛ — в процессе нестационарного электролиза органических соединений, называемых химиками органолюминофорами, на макроэлектродах различных геометрических форм хорошо описан в [21-25], а на микроэлектроде в форме сферы — в [26].

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

U

T,

T4 t

0

Рис. 1. Схема системы пространственных координат для микродиска и биполярного импульсного возбуждения

Анодная фаза (Ag - e ^ A+).

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

+

^ = D 5 2с+

д г2

д t

с+ н

1 дс

-+------+

д 2с+ 2

r д r д z

с + cg = со ,

(1)

(2)

где с+ и Cg — концентрации частиц A+ и Ag

соответственно; D — коэффициент диффузии, который принимается равным для всех частиц,

присутствующих в растворе; со — начальная концентрация органолюминофора.

РИ, 2000, № 3

37

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

t = 0 r > 0 z > 0 1 c+ (r,z,t) = 0

0 < t < T1 0 < r < rd z = 0 1 c+ (r, z, t) = c0

IV 4? z = 0 3 c+ 3 z 0 II 0 II N

r = 0 z > 0 3 c+ 3 r II О II О

r > 0 z c+ (r,z,t) ^ 0 , (3)

Последующие анодные и катодные фазы

( . Л

Ag - e ^ A or Ag + e ^ A

A + + A_ —kbi > 1A* + Ag

1.* kf .

A * У ecl + Ag

V 0 J

Вторая анодная фаза описывается теми же уравнениями (4)-(7), что и катодная, но с другими начальными и граничными условиями:

t = T2 r > 0 z > 0 3 c+ (r, z, t) = 2 c+ (r, z,t);

где rd — радиус микродиска; r и z — пространственные координаты (рис. 1).

Ag + e —^ A

Катодная фаза

A+ + A~ kbi > 1A* + Ag

1 . * kf .

A * У ecl + Ag

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

dc_ д t

3c _ д t

= D

= D

32c+ 13c+ 32c+

-+ —

3 r

2 r 3r dz2

32c 1 dc~

-+—

3 2c

-kbic c , (4)

dr2 r 3r

3 z2

"kbic+c , (5)

3c 31

= D

9 * ^ * 9 *

3 2c 1 3c 3 2c

3r2 r 3r

3 z2

1 __ *

+ kbic c -ф-kf • c

(6)

(7)

c + c + c + cg — c0 .

где c“ и c* — концентрации частиц A_ и 1A*; kbi — бимолекулярная константа скорости; Ф — квантовый выход флуоресценции; kf — псевдомономолекулярная константа скорости в реакциях гомогенного “светового” переноса электрона.

Начальные и граничные условия для катодной фазы следующие:

t = T r > 0 z > 0 2 c+ (r, z, t) =1 c+ (r, z, t);

2c“(r,z,t) = 2 c*(r,z,t) = 0 ; T < t <T2 0 < r < rd z = 0 2c+ (r, z,t) = 0;

2c“(r,z,t) = c0; 2c*(r,z,t) = 0;

r > rd z = 0

d c4

3 z

3 c

z=0

3 z

3 c

z=0

3 z

r > 0

r > 0

z > 0

d c4

3 r

3 c

r=0

3 r

3 c

r=0

3 r

= 0.;

z=0

= 0.;(8)

r=0

z

2 c+ (r, z, t) ^ 0;

3c (r,z,t) = 2 c (r,z,t);

3 c* (r, z, t) = 2 c* (r, z, t) ;

T2 < t < T3 0 < r < rd z = 0

r > rd z = 0

r = 0 r > 0

z > 0

3 c4

3 z

3 c4

3 c

z=0

3 r

3 z

3 c

3 c (r, z, t) = c0;

3 c_ (r, z, t) = 0;

3 c* (r, z, t) = 0; *

3 c

z=0

r=0

3 r

3 z

3 c

= 0. ;

z=0

r=0

3 r

= 0.:

r=0

z ^x

2 c + (r, z, t) ^ 0; (9)

2 c (r, z, t) ^ 0; 2 c* (r, z, t) ^ 0 .

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

Электродный ток определялся из следующего уравнения:

rd

i (t) = 2л-F• D \ 0

( \

3c 3c+

3 z 3 z

z=0 z=0

■ r • dr

где F — константа Фарадея.

(10)

Мгновенная интенсивность ЭХЛ вычислялась как

1ECL (0 = 2л' kf ф ECL JK (r, z, t) ■ r ■ drdz (11)

00

где Ф ecl — коэффициент эффективности ЭХЛ.

2 c (r, z, t) ^ 0;

2 c* (r, z, t) ^ 0 .

38 РИ, 2000, № 3

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

Начальные и граничные условия для анодной фазы:

3.1. Пространственные координаты T = 0 0 <©< 1 0 <Г< 1 С +(©, г,т)=0;

Введем следующие безразмерные параметры: 0 < т < Т 0 <©< 1 Г = 0 С+(©, г,т)=1;

R = — Z = — Cg = ^ С+= — С“= — rd ’ rd ’ c^ c^ c0 ’ д С+

© = 0 0 <Г< 1 5© - 0; (15)

* C* = т = —. c0 ’ Te © = 1 0 <Г< 1 д С+ 5© 0=0 = 0 . ;

Аматоре и Фоссет [10] предложили использовать 0=1

следующие преобразования пространственных координат: 0 <©< 1 Г = 1 О + т) = 0.

R =

1^2

cos | ^ Г

Катодная фаза описывается системой уравнений:

Z = ©• tg|- Г

(12)

д C + M

( о2

д T S

G

5 2С

5Г

2

■+Т

д 2С+

5©

2

- 2Ts

д С+ 5©

(16)

для получения более густой сетки у границ микродиска. Сетка, которая получается с использованием данной трансформации, показана на рис. 2.

Математическая модель описываемого процесса в новых координатах будет выглядеть так:

анодная фаза:

-х-с+с_;

д С ~ M

д T S

( '52С- '

G ---------- + Td

д 2С ~

5Г

2

5©

2

- 2Ts

д С

_ Л

5©

(17)

-х-с+С_;

д С м

( . К2

д T S

G'

д 2С

V

5Г

2

■ + Td

5 2С+

5©

2

- 2Ts

д С ^ 5©

д С M

, (13)

д T S

Ґ О *

' •д2С •

G ----- + Td

v 5Г2 - 5©^

+ Х-С+С_-ф-kf.С* , С ++ С “ + С* + Cg = 1,

5 2С* ,2

- 2Ts

/

С ++ Cg = 1, (14)

где проведена следующая замена переменных:

д С

5©

*

+ (18) (19)

M =-

D • Т

S =-

1

©2 + tg|| г

4 2 [Я

где X Te • kbi *c

0 •

G _ 2 cos I 2 , Начальные

Td = 1 -©2 , Ts =©

фазы:

T = t1

и граничные условия для катодной

0<©<1 0<Г<1,

2С+(©, Г,т) =1 С +(©, г,т),

2С“(©,Г,т) = 2 С*(©,Г,Т) = 0,

т1 <Т<Т2 0<©< 1 Г = 0,

С+(©, Г, т) = 0,

С “(©, г,т) = 1,

С* (©, Г, т) = 0 ,

© = 0

0 <Г< 1

д С4

5©

д С"

0=0

5©

д С

© = 1

0=0 0 <Г< 1

5©

= 0

, (20)

0=0

5 С+ б С “ * б С

5© _ 5© _ 5©

0=1 0=1 0=1

0 <©< 1 Г = 1 С +(©, г, т)

= 0

Рис. 2 Конформные отображения, предложенные Аматоре и Фоссет для преобразования пространственных координат (прямая линия) [10]; Вербрю-же и Бейкер [11] (пунктирная линия)

РИ, 2000, № 3

С “(©, Г, Т) = 0, С*(©,Г,т) = 0 ,

39

Г

Уравнение для вычисления тока:

i(t) = 4rd • F• D• Co j

0

dC~ dC

5Г 5Г

d©

(21)

Мгновенную интенсивность ЭХЛ определяли как:

11 *

7ECL (t) = 4■тс - kf ФECL • c0 • rd И1C* Г’ T)'R' ldet Jl'd&dr,

00

(22)

©=2^2+2/z4 + 2Z2 + 2R2Z2 +1-2R2 + R4 - 2Z2 - 2R2

(27)

5. Результаты и выводы

Основные параметры описанного численного моделирования даны на рис. 3. Таблица содержит вычисленные значения токов (нА) для двух подходов: с неравномерной сеткой только по пространству и с неравномерной сеткой по всем координатам.

© 2 + sin2 flrl (1 -© 2)

„ г ь т| п {2 / ’

где Якобиан det J =-------7--7-----------.

1 1 2 cos3(fг)-вf2

3.2. Экспоненциально-расширяющаяся сетка по времени

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

Общий вид преобразования, предложенного Фель-дбергом [27], следующий:

ln(1 + а—)

T =--------^ . (23)

ln (1 + а)

В нашем случае, учитывая периодичность процессов, время трансформировалось из следующего соотношения:

ґ

ln

а

T = -

^(t - Te) +1

у!ё_______;

ln (1 + а)

+ (Nph -1), (24)

где Tph — длительность текущей фазы (анодной или катодной); Te — длительность предшествующих фаз электролиза; Nph — номер текущей фазы. Отсюда текущее время электролиза:

t = biL [ e'"(1*“lT-<Nph-‘1 -1) + Te. (25)

Размер сетки Неравномерная сетка

по R, Z по R, Z и Т

25 х25 85. 895 85. 891

50 х 50 86. 946 86. 942

75 х 75 86. 815 86. 811

100 х 100 86. 670 86. 666

150 х150 86. 474 86. 470

200 х 200 86. 357 86. 353

300 х 300 86. 227 86. 223

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

4. Подробности численного моделирования

Для быстрого получения всех вычисленных значений концентраций применялось обратное преобразование трансформированных координат (©, Г) в реальные (R, Z):

^ 2 t 2Z

Г=—arctan,—

2 применением численных методов. Оптимизация выбора численных методов и подходов к конкрет-(26) ной задаче - одна из основных задач современного

РИ, 2000, № 3

П+Wz4 + 2Z2 + 2R2Z2 +1-2R2 + R4 - 2Z2 - 2R

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

40

(t=0,1c)

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

На рис. 6, 7 показаны распределения интенсивности ЭХЛ при различных соотношениях длительно -стей анодной и катодной фаз.

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

Литература.1. Gavaghan D.J. // J. Electroanal. Chem. 456. 1998. Р. 1-12. 2. Gavaghan D.J. // J. Electroanal. Chem. 456. 1998. Р 13-23. 3. Gavaghan D.J. // J. Electroanal. Chem. 456. 1998. Р. 25-35. 4. J. Crank, R.M. Furzeland// J Inst. Maths. Applic. 20 (1977) 355. 5. D. J Gavaghan, J.S Rollettn // J Electroanal Chem 295. 1990. 1. 6. Galceran J., Gavaghan

Рис. 6. Вычисленное распределение интенсивности ЭХЛ ( Та = 0,0005 с, Тк = 0,0003 с)

грамме с неравномерной сеткой по времени (t=0,1с)

DJ, Rollett J. S. // J Electroanal Chem 394. 1995. 17. 7. Taylor G, Girault H.H., McAleerJ. // J. Electoanal. Chem. 293. 1990. 19. 8. Michael AC, WightmanR.M., Amatore C.A //J. Electroanal. Chem. 267. 1989. Р. 33. 9. Lavagnini I., Pastore P, Magno R, Amatore C.A. // J. Electroanal. Chem. 316. 1991. Р.37. 10. AmatoreC.A, Fossett B. // J. Electroanal. Chem. 328. 1992. Р. 21. 11. Verbrugge M.W, Baker D.R // J. Phys. Chem. 96. 1992. Р.4572. 12. EvansN.T.S., Gourlay A.R. // J. Inst. Maths. Applic. 19. 1977. Р.239. 13. Heinze J, Storzbach M. // Ber. Bunsenges. Phys. Chem. 90. 1986. Р.1043. 14. Aoki K., Osteryoung J. // J. Electroanal. Chem. 122. 1981. Р.19. 15. Oldham K.B. // J. Electroanal. Chem. 122. 1981. Р.1. 16. Oldham K.B. // J. Electroanal. Chem. 260. 1989. Р. 461. 17. Bond AM, Oldham K.B, Zoski C. G. // J. Electroanal. Chem. 245. 1988. Р. 71. 18.Fleischmann M, Pons S. // J. Electroanal. Chem. 250. 1988. Р. 257. 19. Fleischmann M, Daschbach J., Pons S. // J. Electroanal. Chem. 250. 1988. Р. 269. 20. Cope D.K., Tollman D.E. // J. Electroanal. Chem. 373. 1994. Р. 53. 21. Faulkner L.R, Bard A.J. // J. Amer. Chem. Soc. 90. 1968. Р.6284. 22. LaserD, BardAJ. // J. Electrochem. Soc. 122. 1975. Р.632. 23. Maloy J. T, PraterK.B., Bard AJ. // J. Amer. Chem. Soc. 72. 1968. Р.4348. 24. Maloy J. T, Prater K. B., Bard A. J. // J. Amer. Chem. Soc. 93. 1971. Р.5959. 25. БыхА.И, ВасильевР.Ф, Рожицкий Н.Н./ Итоги науки и техники. ВИНИТИ. Сер. Радиац. химия. Фотохимия. 1979. №2.135 с. 26. Svir I.B., Oleinick A.I. and Compton R. G.// Радиоэлектроника и информатика, № 1, 2000, с. 28. 27. Feldberg S. W.// J.

Рис. 7. Вычисленное распределение интенсивности ЭХЛ ( Та = Тк = 0,0005 с)

РИ, 2000, № 3

41

Рис. 8 Вычисленная концентрация катион-радикалов в трансформированных координатах (©, Г)

(t = 0,0012 с, Ta = Tk = 0,001 с, % = 7 микрон)

і Z, 0.001 cm

1 -

R. 0.001

О 1

Рис. 9. Вычисленная концентрация катион-радикалов в реальных координатах (R, Z) (t = 0,0012 с, Ta = Tk =

= 0,001 с, id = 7 микрон)

Electroanal. Chem., 127, (1981), 1. 28. Белаш Е.М., Ро-жицкий Н.Н. // Электрохимия, 1996, Т.32, № 11. с.1263.

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

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

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

УДК 621.317.799

ПРИНЦИПЫ ПОВЫШЕНИЯ

ТОЧНОСТИ КОНТРОЛЯ

СТРУКТУРНОГО

НЕСОВЕРШЕНСТВА

ПОЛУПРОВОДНИКОВЫХ

ПЛАСТИН

ОКСАНИЧ А. П. , ВДОВИЧЕНКО Н. Д._

Рассматриваются схемы инфракрасных полярископов, применяемых для измерения структурных несовершенств в полупроводниковых пластинах. Определяются пути снижения порога чувствительности при измерении внутренних напряжений в полупроводниках. Приводится описание разработанного метода и аппаратуры, позволяющей снизить чувствительность определения разности главных напряжений до 0, 07кг/см2 .

1. Введение

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

42

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

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

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

РИ, 2000, № 3