Численное и асимптотическое решения неодномерной системы уравнений Нернста-Планка-Пуассона Текст научной статьи по специальности «Математика»

УДК 517.91

ЧИСЛЕННОЕ И АСИМПТОТИЧЕСКОЕ РЕШЕНИЯ НЕОДНОМЕРНОЙ СИСТЕМЫ УРАВНЕНИЙ НЕРНСТА-ПЛАНКА-ПУАССОНА

© 2010 г. А.В. Лаврентьев, К.М. Уртенов, А.А. Хромых, Н.О. Чубырь

Кубанский государственный технологический университет, Kuban State Technological University,

ул. Московская, 2, г. Краснодар, 350072, ул. Moskovskaya St., 2, Krasnodar, 350072,

pm@kubsty.ru pm@kubsty.ru

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

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

We offer hereby research of non-singular dimensional system of equations of Nernst-Plank-Puasson for binary electrolyte in decomposition form. We discuss a model task with condition of electroneutrality. We investigate simplest model task taking into consideration special charge in nonsingular dimension case for stationary transition of binary electrolyte.

Keywords: decomposition, electric current density, electrodialysis, systems of equations systems by Nernst-Plank-Puasson, binary electrolyte, modeling, electromembrane systems.

Структура системы уравнений Нернста-Планка и Пуассона неудобна для применения численных методов. Кроме того, при естественной нормировке возникает малый параметр при операторе Лапласа в уравнении Пуассона. В связи с этим появляется проблема преобразования этой системы к более удобному виду.

Для одномерного случая в [1, 2] предложен метод декомпозиции, использование которого привело фактически к созданию теории переноса произвольного электролита [3-6].

В [7] нами дано обобщение метода декомпозиции на неодномерный случай, получены декомпозиционные уравнения и выведено новое уравнение для плотности тока.

Модельная задача с условием электронейтральности подробно исследована в [8]. В данной работе предлагается анализ краевой задачи, моделирующей учет пространственного заряда при переносе бинарного электролита в двумерном случае.

Декомпозиция неодномерной системы уравнений Нернста-Планка-Пуассона. Приведем сначала краткий вывод декомпозиционной системы уравнений. Исходная система уравнений Нернста-Планка и Пуассона в отсутствии химических реакций для бинарного электролита имеет вид [9]:

¿. = -¡^г^С вгаа р - Б. ^ С, + С, V , г = 1,2, (1) К1

dCj ~dt

= - div j, i = 1,2,

(2)

т.е. можно произвести расщепление (декомпозицию) исходной системы уравнений.

Положим £0 = С1 + С2, ~ = £0 + |Е|2,

Д дЕ

О = е0 — +1.

0 дЛ

Тогда декомпозиционная система уравнений относительно Е и £ примет вид [7]:

~ ^&ГНЕ2Е) -

— = ——2—1-2. div(~E)-dt RT ' 2d1R 2T2

S0—i(div E)2 3||grad I2 +

RT

| g0d2 F

RT

A(div E)- div(.~V) -

g0—Ldiv( E|2 V) - diA~

2d1RT

+ g0d3

2d1RT dt 1 '

(5)

e0 — = F2—3 —l—2 SE - F2g0—2 —l—2 EE2 -0 dt RT 2—R 2T2 1 1

i1R

Fe0 d 4

RT

s0 Fd 2—з —1—2

E div E - Fd 2 z1 z 2 grad S +

grad|E2 + s0d3d4AE - s0(VdivE)+ O.

£)Ар = -^С + 22С2) , (3)

где С., \г - концентрация и поток ионов i-го вида; Бг и г. - коэффициенты диффузии и зарядовые числа ионов /'-го вида; р - электрический потенциал; Б -число Фарадея; е0) - диэлектрическая проницаемость вакуума; V - вектор скорости течения раствора электролита; R - универсальная газовая постоянная; Т -абсолютная температура.

Пусть I - плотность тока, обусловленного потоком ионов; Е - напряженность электрического поля. Тогда

I = F{z1j1 + г2^), Е = -§га! р, а уравнения (1) и (3) примут вид:

\г =-^г1Б1С1Е - Б. вгаа С. + С. V , г = 1,2, (4) К1

80&у Е = Е (гС + г2С2 ).

Вместо исходной системы уравнений (2), (4) для определения неизвестных С, С2, ^, ]2 , Е можно получить 2 уравнения для 2 неизвестных функций. Остальные неизвестные могут быть найдены по формулам или по отдельным независимым уравнениям,

2d1RT

div О = 0 . (6)

В систему уравнений (5) наряду с неизвестными функциями E и S входит и неизвестная общая плотность тока О . В [10] получены и исследованы системы уравнений, подобные (5) в предположении о различных частных случаях функции О .

Для решения задачи в общем виде система (5) должна быть дополнена уравнениями для определения вектор-функции О . Одно из них - (6). В [7] найдено уравнение для общей плотности тока О .

В трехмерном случае оно имеет вид:

т F

rot I = — grad х RT

( ш 2 Л

- Fd3 z,z2 ~ + g°Fd3 ZlZ2 |E| 2 + e0 d4 divE

3 1 2 2d1RT 04

xE + e°AE х V + e° divE х rot V , в двумерном -

дФ 2 дФ

dx (

dy

F RT

grad

-Fd3 —1—2S + g°Fd32—1-2 El2^4 divE

2d1RT

^ ^

E

/ У

+ g0(AE, V)1 +g0

dV2 dV1

dx dy

divE,

+

x

x

+

\

где для удобства записи использовано обозначение (а,Ь)1 = - .

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

Модельная задача с условием электронейтральности. В этом случае вместо уравнения Пуассона (3) используется условие электронейтральности ZlCl + 22 С 2 = 0, что формально эквивалентно £о = 0, причем система декомпозиционных уравнений значительно упростится.

Так, например, в двухмерном случае дSo

dt

E = -

■ = divido V) + DASo,

d2 RT Fd3Soo

grad So +-

RT

So F 2d^z,z

3z1z2

F 2d3z1z2

A^ = -

где D =

RT

did2(zi - Z2) DiZi - Dz'

(grad So, E)i

дц

(7)

(8)

(9)

'2 z2

dx

= h

dV _ 7 — = —I i

dv 1

RT

2^2

2d1R 2T

— Fd2 z1z2 grad S +1 = 0,

A^ = —

F2d3z1z2 f

RT

grad

S-■

Sn. d-

0d3

2d,FRT

Hl2

Л А E

,(d )

„(d )

У = -

H(d)

J (d ) H (d )

D(d)C(d)F ' (d)

H

(d )

C =

C(d)

C (d) Co

V =

V

(d)

V

(d)

RTs,

s = ■

o

C (d) F 2 H 2

= 2

lD

H

E =

H(d) FE RT

D =■

D

(d )

i = 1,2, Pe = -

V(d) H(d)

)'■ в(а)

безразмерное критериальное число Пекле, где переменные с индексом (d) считаются размерными.

Модельная задача стационарного переноса бинарного электролита с условием квазиравномерного распределения заряда в случае симметричного электролита и равенства коэффициентов диффузии О1 = О2 = О в безразмерном виде принимает вид: О

dS = D AS — div(~, V),

Я* D„ \ 1

dh

dx

s

dl

= —Dz1z2 f grad ^ ~ +1 s|e|

E

1

(e, E)E+SE+■

2 Dz 1z 2

J = 0,

(10)

(11) (12)

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

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

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

Ее аналогом в двухмерном стационарном случае является задача [10]: ^212 &у(~е)-&у(~У) -= 0,

Р%2122 Р2122 Е |Е|2

<ИУ I=0 .

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

1 1 1 1

[пд = у 11х(Х У)\х=0¿У = у 11х(ХУ)\Х=^У = сотИ . (13) Ь 0 Ь 0

дц дц

дх 2' дУ Из (13) выпекает, что для функции ц можно положить

ц У=0 = 0 Ц У=1 =-пзЬ. (14)

Кроме того, добавим условие перпендикулярности тока поверхности мембраны

Из (12) следует

= 12, = —11, J = (I1,12) .

d^ dx

x=o

d^ dx

= o.

(15)

x=h

, (16)

Для функции £ примем условия:

£ = >£Л , £ , = ^К , £ „ = £ао , £ , = 1х=0 'Х=п 1у=0 1у=1

которые должны быть согласованы со свойствами мембран и с величиной («интенсивным» или

«мягким» токовым режимом).

Асимптотическое представление решения. Пусть

ß = —

1

Dz^

а = —. Для моделируемого случая ве-

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

личина а может считаться малым параметром по сравнению с р. Это позволяет построить самосогласованное асимптотическое решение задачи (10)-(16). В [11] показано, что

2 1 —1 ?3~ 3

E = |1| 3 ß3a 3z(—S|l| 3 ß 3а 3)J

2 —1

33

(17)

а функция является положительным решением уравнения

J

J

1

z3-£z-1 = 0

и имеет асимптотическое представление 1

(18)

z(£) =

—, iöe £ ^ -да

1, iöe £ ^ 0 , iöe £ ^

ß^Vo = (grad (~), ßß Io)i, S

(19)

dVo r dVo т а

причем —0 = 12 0, —0 = -I10 , а в области про-

дх

ду

ъ V-2~ т т

странственного заряда E » м_ и t— I0 , где 10 - ре-

N1

д 1 ~ ^1,0 (~(grad( ~), 10)):

' ду s

div 10 = 0 .

(20)

Таким образом, асимптотическое представление для функции Е зависит от знака функции £ . При интенсивных токовых режимах £ < 0, в то время

1х=к

как £ > 0, поэтому функция £ меняет знак в

1х=0

области x е [0,Щ, у е [0,Ц, следовательно, функция Е имеет в разных частях области различные асимптотические представления (рис. 1).

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

После нахождения функции S вычисляем зависимость x=x(y) из условия, что S (х(у), у) = 0, причем при x>x(y) функция S (х, у) < 0, а при x<x(y) -S (х, у) > 0 . Область, где S (х, у) > 0, соответствует

области электронейтральности, S (х, у) < 0 - области пространственного заряда.

Полученную асимптотику функции z(£) необходимо согласовать с асимптотикой функции I для того, чтобы получить самосогласованные асимптотические решения.

Для этого предположим, что lim I(х, у) = 10(х, у).

Подставим (17) в (11) и для 10 получим соответствующую задачу.

На основе этого алгоритма показано, что в области

электронейтральности E « ßß I0 , где 10 - решение

уравнения

Для уравнения (19) используются левые граничные условия (14), (15), для (20) - правые. На границе (х(у), у) ставятся условия согласования решений уравнений (19), (20), причем для этого необходимо ввести промежуточный слой в окрестности кривой (х(у), у).

Асимптотическое решение. Как следует из предыдущего пункта, для асимптотического решения задачи нужно использовать следующие разложения:

Е1 Е1-1 +Ею +0(£), л/е

Е2 =±Е2-1 +Е20 +0(£), V = Ую + 0(£

V £

V2 = ¥20 + 0(у[£), ~ = ~0 +~1 + 0(£),

11 = /ю + 0(л/£), 12 = 120 + 0(>/£).

Уравнения для коэффициентов разложения распадаются на ряд подсистем. При их решении можно выделить 2 различных случая. В 1 -м имеем решение:

Е1 +..., Е2 = 0 ,

ые

__

/1 = С( х) )к +Б, / 2 = 0.

Используя уравнения для коэффициентов следующих приближений, можно найти остальные приближения для плотности тока, а значит, и напряженности. Во 2-м случае получим для / 0 уравнение

Т2 д2Л ,1Т Т д2У , т2

11,0^Т + 1,012,0 +12,0 2

дх

2,0 ^ ^ 2,0 Л

дхду ду

2

( д~0 т д~0^2 - 40^--12,0

ду дх

11,0 +12,0

2Sn

(21)

где

dv

= I

20,

dV _ j — = -110-

дх 20' ду Делая замену переменных

£( х, у) = - /2,0 (х, у^ +

д

+ f ^/1,0 (х, у) + -д (f12,0 (х, у)«^)^ , и = х , приведем уравнение (21) к каноническому виду.

д2V * , dV

И'110

(12 д2£+ 21 / 11,0-2 + 211,012,0

дх2 дудх ду

д 2£

д 2£

+12,0 _ 2

дSo с£0

---12,0

Су дх

22 11,0 +12,0

2Sn

(22)

шение системы уравнений

Уравнение (22) является квазилинейным. Для него известна локальная теорема единственности. После нахождения /^,0( х, у), /2,0 (х, у) определяются

Е1,-1(x, y), Е2,-1(х у).

Нами найдены различные частные решения уравнения (22).

9

Для нахождения Е1 -1(х, у), Е2 -1 (х, у) необходимо совместно решить уравнение из 1 -го приближения

9 9

Е^-1 + Е2,-1 = -2£0 и уравнение, которое получается из условия разрешимости следующего приближения для

Е1,-Ь Е2,-1, т.е.

E2-1 h,o — E1,—112,o =0. В результате по-

лучим Е

1—1

±1

1o

-2sn

Т 2 + Т 2 т1п + т 9

И E2—1 = ± I

2o/

— 2So

j2 + Т-2

T1o + Т 2o

'10 +1 20

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

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

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

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

^=0 = ~0(х;У), цг=0 = 0,

ё £х=0 = £0 (пУ), £х=н = у)> (23)

где функции , £0 , , являются заданными.

Рассмотрим результаты численного анализа краевой задачи (11)-(18), (23) при 2 различных вариантах граничных условий.

Первый случай. Краевые условия имеют вид

~0(х;у)= -1, ~0(';у)= 1, ^С';у)=-1 , £0(';х)= 1 - 2х,

х) = 1 - 2х. На рис. 2 приведены графики функций £ (/; х; у) и Е^; х; у).

Заметим, что в данном случае имеется точное решение £ (?; х; у) = 1 - 2х. Из сопоставления видно, что численное приближенное решение с большой точностью совпадает с точным.

График функции Е^; х; у) соответствует физически ожидаемому поведению функции.

Второй случай. Краевые условия имеют вид

~0(х;У)= 2, £0(';у)= 1, 51(';у)=-1 , £0(';х)= 1 - 2х,

5*1 (/; х)= 0. На рис. 3 приведены графики функций

£ (г, х; у) и Е^; х; у).

Рис. 2. Графики функций: а - £(¿; х; у) ; б - Е1 (¿; х; у) при 1-м случае граничных условий

Рис. 3. Графики функций: а - £(¿; х; у) ; б - Е1 (¿; х; у) при 2-м случае граничных условий

а

б

а

б

а

Решение близко к стационарным, полученным в предыдущем случае уже при /соп = 0,5 .

Сопоставление решений, полученных асимптотическим и численным методами. Сопоставление проведено в области х > 1/2, у е [0;1], причем для асимптотического решения согласовывались граничные условия с численным решением. Первое приближенное решение для г получено в виде г = -0,01 - у. Сопоставление численного и асимптотического решений приведено в таблице.

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

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

Сопоставление численного и асимптотического решений для функции ц

\ x y \ 0,6 0,7 0,8 0,9 1

0 -0,23 -0,11 -0,04 -0,01 0

0,1 -0,33 -0,21 -0,14 -0,11 -0,10

0,2 -0,43 -0,31 -0,24 -0,21 -0,20

0,3 -0,53 -0,41 -0,34 -0,31 -0,30

0,4 -0,63 -0,51 -0,44 -0,41 -0,40

0,5 -0,73 -0,61 -0,54 -0,51 -0,50

0,6 -0,83 -0,71 -0,64 -0,61 -0,60

0,7 -0,93 -0,81 -0,74 -0,71 -0,70

0,8 -0,03 -0,91 -0,84 -0,81 -0,80

0,9 -1,13 -1,01 -0,94 -0,91 -0,90

1,0 -1,23 -1,11 -1,46 -1,01 -1,00

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

Поступила в редакцию_

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

Литература

1. Декомпозиция систем уравнений Нернста-Планка-Пуассона / В.А. Бабешко [и др.] // Докл. РАН. 1995. Т. 344, № 3. С. 485-487.

2. Декомпозиционные уравнения для стационарного переноса электролита в одномерном случае / В.А. Бабешко [и др.] // Электрохимия. 1997. № 8. С. 855-863.

3. Теория стационарного переноса тернарного электролита в одномерном случае / В.А. Бабешко [и др.] // Докл. РАН. 1997. Т. 355, № 4. С. 488-491.

4. Теория стационарного переноса бинарного электролита в слое Нернста / В.А. Бабешко [и др.] // Докл. РАН. 1998. Т. 361, № 2. С. 208.

5. Уртенов М.Х. Краевые задачи для систем уравнений Нернста-Планка-Пуассона. Краснодар, 1998. 126 с.

6. Лаврентьев А.В., Уртенов М.Х. Метод регулярного представления сингулярно возмущенных уравнений. Краснодар, 2002. 134 с.

7. Полная декомпозиция неодномерной системы уравнений Нернста-Планка-Пуассона для бинарного электролита / А.В. Лаврентьев [и др.] // Эколог. вестн. науч. центров Черноморского экономического сотрудничества. 2009. № 2. С. 32-37.

8. Лаврентьев А.В., Письменский А.В., Уртенов М.Х. Математическое моделирование переноса в электромембранных системах с учетом конвективных течений. Краснодар, 2006. 146 с.

9. НьюменДж. Электрохимические системы. М., 1977. 463 с.

10. Уртенов М.Х., Сеидов Р.Р. Математические модели электромембранных систем очистки воды. Краснодар, 2000. 140 с.

11. Уртенов К.М. Алгоритм численного решения двумерной модели стационарного переноса бинарного электролита в ЭМС // Современное состояние и приоритеты развития фундаментальных наук: тр. IV Всерос. науч. конф. молодых ученых и студентов. Краснодар, 2009. С. 112-115.

26 января 2010 г