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

ФИЗИКА

PHYSICS

УДК 53.1:51 ББК 22.311 М 74

Коваленко А.В.

Кандидат экономических наук, доцент кафедры прикладной математики Кубанского государственного университета, тел. 89184440042, e-mail: savanna-05@mail.ru Хромых А.А.

Старший преподаватель кафедры прикладной математики Кубанского государственного технологического университета, тел. 89184665266, e-mail:AnnXA@mail.ru Уртенов М.Х.

Доктор физико-матаиатических наук, профессор, зав. кафедрой прикладной математики Кубанского государственного университета, тел. 89184659466, e-mail: urtenovmax@mail.ru Мамий Д.К.

Кандидат физико-матаиатических наук, доцент, зав. кафедрой алгебры и геометрии, декан факультета математики и компьютерных наук Адыгейского государственного университета, тел. (8772) 59-37-32, e-mail: dmami@yandex.ru

Моделирование и численный анализ процесса переноса бинарного электролита в канале обессоливания электродиализного аппарата в потенциостатическом режиме

(Рецензирована)

Аннотация

Рассмотрено моделирование и численный анализ процесса переноса бинарного электролита в канале обессоливания электродиализного аппарата. Предложены различные математические модели переноса ионов в потенциостатическом режиме в виде системы квазилинейных уравнений с частными производными. Осуществлен переход к безразмерному виду и показано наличие некоторого малого параметра £ > 0 при старина производных в уравнениях. Если скорость прокачки электролита через камеру обессоливания значительная, то возникает еще один малый параметр Я> 0. Для численного анализа используется специальная растянутая система координат. Выявлены основные закономерности возникновения и развития пространственного заряда в канале обессоливания электродиализного аппарата. Предложены различные асимптотические и численные методы решения.

Ключевые слова: мат&матическое моделирование, потенциостатический режим, электродиа-, , .

Kovalenko A.V.

Candidate of Economics, Associate Professor of Applied Mathematics Department, Kuban State University, ph. 89184440042, e-mail: savanna-05@mail.ru Khromykh A.A.

Senior Lecturer of Applied Mathematics Department, Kuban State University of Technology, ph. 89184665266, e-mail: AnnXA@mail.ru Urtenov M.K.

Doctor of Physics and Mathematics, Professor, Head of Applied Mathematics Department, Kuban State University, ph. 89184659466, e-mail: urtenovmax@mail.ru Mamiy D.K.

Candidate of Physics and Mathematics, Associate Professor, Head of Algebra and Geometry Department, Dean of Mathematics and Computer Science Faculty, Adyghe State University, ph. (8772) 59-37-32, email: dmami@yandex.ru

*

РФ в рамках темы 1.4.08 Единого заказ/наряда.

Modeling and the numerical analysis of the process of binary electrolyte transfer in the desalination channel of the electrodialytic device in the potentially static mode

Abstract

The paper discusses the modeling and the numerical analysis of the process of binary electrolyte transfer in the desalination channel of the electrodialytic device. The authors propose various mathematical models of ion transfer in the potentially static mode in the form of system of the quasilinear equations with private derivatives. Transition to a dimensionless kind is carried out and presence of some small parameter £ > 0 is shown at the senior derivatives in the equations. If speed of electrolyte circulation through the desalination chamber is considerable there is one more small parameter A> 0 . The special stretched system of coordinates is used for the numerical analysis. The basic laws of occurrence and development of a spatial charge in the desalination channel of electrodialytic device are revealed. Various asymptotic and numerical methods of the solution are proposed.

Keywords: mathematical modeling, the potentially static mode, electrodialytic device, electromembrane systems, numerical and asymptotic methods.

Введение

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

Краевые задачи для систем уравнений Нернста-Планка и Пуассона ставились и изучались многими авторами. Однако в имеющихся работах обычно ограничивались постановкой и рассмотрением одномерных задач [3]. В данной статье рассматривается двухмерный случай.

1. Постановка задачи

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

Пусть H и L - ширина и длина камеры обессоливания, соответственно, x = 0

- соответствует условной межфазной границе «анионообменная мембрана/раствор», x = H - соответствует условной межфазной границе «катионообменная мембрана/раствор», y = 0 - входу, а y = L - выходу из камеры обессоливания, V - заданная скорость прокачивания раствора. Будем считать, что скорость течения раствора имеет форму параболы Пуазейля или моделируется по-другому.

Для математического моделирования явлений переноса для бинарного электролита используются следующие уравнения [1]:

— F -

jt =——ztDtCFv-ДУ Ct + CtV, i = 1,2 ,

£0Л(р = - р (^1С1 + 2 2С 2 ) , (3)

1 = Р(1 + 22к )• (4)

Здесь V - градиент, А - оператор Лапласа, V - скорость течения раствора, С1, С2 - концентрации катионов и анионов в растворе, соответственно, г2 - зарядовые числа катионов и анионов, Д , В2 - коэффициенты диффузии катионов и анионов, соответственно, р - потенциал электрического поля, е0 - диэлектрическая проницаемость электролита, Р - постоянная Фарадея, К - газовая постоянная, Т - абсолютная температура, т - время, I - плотность тока.

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

, С, (, = 1,2), р плотность тока I находится по формуле (4). Напряженность и потенциал электрического поля связаны соотношением Е = -V р.

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

^ = В,АС., -а^Су]+-^-грМС?}р), , = 1,2, (5)

от КТ

£0Ар = - Р (1С1 + 2 2С 2 ). (6)

2. Модель переноса симметричного бинарного электролита

Важным частным случаем является симметричный бинарный электролит, когда 2Х = -22. Не нарушая общности рассуждения, в этом случае можно положить г1 = 1 и г2 = -1, тогда уравнения (5), (6) упрощаются и принимают вид:

¡С=в, ас, - <*,(с, V)+ КТоМсрр), (7)

дС2=в2ас - )-КТв-Мсур), (8)

£0А(Р = - Р (С1 - С 2 ). (9)

3. Начальные и краевые условия

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

В этой работе мы будем рассматривать потенциостатический режим, которому соответствует условие:

p(t, H, у) — p(t, 0, y) = dp = const, (10)

означающее, что величина падения потенциала в канале постоянна. Поскольку важна только величина падения потенциала, то, не нарушая общности, можно положить, например, p(t, 0, y) = dp, тогда p(t, H, у) = 0. Каждому значению падения потенциала

dp соответствует своя плотность тока

1 L

L(t) = у JX ( x У )Ф, (11)

L 0

являющаяся функцией от времени, но не зависящая от x.

Все граничные условия должны быть согласованы со свойствами мембран и с величиной dp («интенсивным» или «мягким» токовым режимом).

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

1). На поверхности анионообменной мембраны x = 0, у е[0,L], t > 0 будем считать

граничную концентрацию анионов равной фиксированной концентрации носителей за-

ряда внутри мембраны:

C2(t,0,у) = Сат . (12)

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

ґдС1 F ^ д®Л

—L + —

У дх RT дх j

(t ,0, у) = 0. (13)

2). На поверхности катионообменной мембраны х = Н, у є [0, Ь], г > 0 будем считать граничную концентрацию катионов равной фиксированной концентрации носитетелей заряда внутри мембраны:

С(г, Н, у) = Скт. (14)

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

ҐдЄ2 F ^ дфЛ

—L+ z2C2 —

v дх RT дх j

(t, H, у) = 0. (15)

3). На входе в рассматриваемую область у = 0, хє[0,Н], г > 0 будем полагать:

Сі(г,х,0) = Сі0(г,х), і = 1,2, (16)

= 0. (17)

дУ

4). На выходе из рассматриваемой области у = Ь, х є[0, Н], г > 0 будем использовать «мягкие» условия на концентрации и на потенциал:

дС-(гх,Ь) = 0, і = 1,2, (18)

дУ

= 0. (19)

дУ

5). Начальные условия при / = 0 примем согласованными по-возможности с остальными граничными условиями:

Сг (°, х, у) _ Сгп (х, у), г _ 1,2. (20)

Для функции р в качестве начального условия можно взять, например, линейную функцию, независящую от у :

( х Л

Р(0, Х у) = 1 ^р. (21)

V Н У

Для решения краевой задачи применяется метод конечных элементов [4].

4. Характерные данные

Нами были проведены численные эксперименты для раствора ЫаС! в широком спектре таких параметров, как начальная концентрация, межмембранное расстояние, длина канала, значение падения потенциала ёр и определены основные закономерности распределения электрохимических (концентрация, напряженность и потенциал электрического поля и т.д.) полей. Ниже представлены некоторые результаты численных экспериментов при следующих входных параметрах: ширина канала обессолива-ния Н _ 1 мм или меньше, длина канала Ь _ 7 мм, средняя скорость вынужденного течения раствора варьировалась от ¥0 _ 2 10-6 м/с до ¥0 _ 10-2 м/с, а падение потенциала от йр_ 0,2 В до йр_ 0,6 В, начальная концентрация раствора принимает зна-

3 3

чения от С 0 = 10 моль/м до С 0 = 100 моль/м, температура раствора Т _ 293 К, начальная плотность раствора р0 = 1002,5 кг/м , коэффициент кинематической вязкости у_ 1006 10-9 м2/с, коэффициент диффузии катиона и аниона, соответственно,

Д = 1,33 • 10-9 м2/с и Д _ 2,05 • 10-9 м2/с.

5. Переход к безразмерному виду и оценка безразмерных параметров

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

х (и) _ Х у(и) _ у с (и) _ Сг С (и) _ Скт с (и) _ Сат

тт , У тт , г , кт , ат ,

Н Н С0 С0 С0

V(и) _ —, 7(и) _ —— 1, Ё(и) _ —Ё, ри) _— р, Ди) _,

—0 ДС0 КТ ЯТТ г Д

где, например, Д _ Д или Д _ 10 9 м2/с, I(и) _—Н—I, I(и) _ — I, —и) _ — .

ДС0 — Н Н

Кроме того, введем в рассмотрение безразмерные параметры:

др ) _-—^р - безразмерное падение потенциала, Ре _ —Н- - число Пекле,

КТ А

є -

ЯТєп

• = 2

Н

ЯТєп

- Дебаевская длина [1].

Эти безразмерные параметры имеют следующие характерные пределы изменения:

—

х е [0,1], у е [0,7], ё(р) _-ёр _ кёр меняется от 5 до 25, число Пекле Ре линейно

КТ

зависит от начальной скорости и при ее изменении меняется от 10-4 до 106, £ меняется от 10-4 до 10-10 в зависимости от значений С0, Н . Таким образом, £ можно считать малым параметром. При больших числах Пекле Ре число Л_ 1/ Ре также будет малым параметром, например, при Ре > 104 число Л < 10-4.

В безразмерных переменных исходная система уравнений (1-4) имеет вид (индекс и для упрощения записи опущен):

1г _ -гАСУр-ДУ С, + РеС—, г _ 1,2 , дС

Ре —г- _ -й1\ 1, г _ 1,2 дг Л

є А(р - -(ііСі + 22С2),

1 - 2і1і + 22%.

Соответственно изменяются и системы уравнений (5, 6) и (7-9), например, система уравнений (5, 6) примет следующий безразмерный вид:

Ре- ВгАСг + Ур)-РеОу(су^ і - ^

є Ар - -(2іСі + 22С2 ) .

(22)

Для системы уравнений (21), (22) безразмерные краевые условия примут вид: 1). При х _ 0 :

Р(і,° У) - dр, С2(ІЛ У) - Сат

2). При х -1:

РОЛ У) - 0, С1(і,1, У) - Скт,

3). При у - 0 :

дСі + гС дР

-------+ ге, —

дх дх

(і, 0, у) - 0.

(і,1, У) - 0.

др(і, х,0) дУ

- 0, Сг (і, х, 0) - Сг 0 (і, х), і -1,2,

(23)

(24)

(25)

4). При у - Ь :

Эр(/, х, Ь) - 0 дСг(і,х, Ь) - 0 . -12.

Эу ду

5) При г _ 0 :

р(0, ^ у) _ (1 - х) ёр , Сг (0, х у) _ Сг,п (х у), г _ 1,2

Наличие малого параметра £ не позволяет решать краевую задачу (21)-(27) в исходных переменных при £ меньше 10-2, поскольку в области погранслоя необхо-

(26)

(27)

2

I

а

дима дискретизация хотя бы с шагом И /10.

Например, если е имеет порядок 10 4, то шаг должен быть порядка от 10 5 и меньше, количество элементов при этом будет порядка от 1010 и больше. Попытки вне погранслоя увеличить шаг приводит к «взрыву погрешностей» [5].

Для решения краевой задачи при е меньших порядков, чем 10 2, введем специальную растянутую систему координат и переменных:

% = хе 3, % =е 3, д = уе 3, д = Ье 3, т = tе 3, (р = Ф,

Е = Ое 3, Сг = е 3, Ж = Vе 3.

Тогда

дt дт дх д% ду дд

дх2 д% ’ ду2 дд2 '

Следовательно,

_1 -1 -- -1

ЛСг = е"3 Л, (С,.Vр) = е 3VФ), Л(р = е 3 ЛФ, ^у(С/) = е"3).

В этих формулах дифференциальные операторы слева считаются в переменных (х, у), а справа - в (%, д).

Подставляя эти формулы в уравнения (21), (22) в растянутых переменных, получим уравнения:

Отметим, что внешний вид уравнений не изменился, за исключением того, что в новые уравнения малый параметр уже не входит. Область изменения переменных растянулась, т.к., например, при е = 10-3 прямоугольник [0,1]х[0,7] становится прямоугольником [0,10]х[0,70], и поэтому количество узлов увеличивается (в 100-210 раз), однако это значительно меньше начального варианта и уже вполне приемлемо. Растянутые переменные позволяют проводить расчеты вплоть до е = 10-6.

Краевые условия в растянутых переменных имеют вид:

(28)

Лф = -(+ 22£2) .

(29)

1). При % = 0:

(30)

3). При д = 0:

1 1 1

:0, І = 1,2. (33)

ЭФ(т/’0) = 0, Н,(т,|,0) = е 3€,о(е3т, е3£), І = 1,2. (32)

4). При д = д:

дФ(т,£д) = 0 дН, (т,£д)

дд ’ дд

5) При т = 0:

1 -1 1 1

Ф(0,£д) = є3£сі„, 5,(0,|,д) = е3С1п (еТ, е3!;), і = 1,2. (34)

6. Основные закономерности переноса ионов в канале обессоливания

6.1. Возникновение и развитие пространственного заряда.

Структура канала обессоливания

Постепенно увеличивая ёр, получаем, что пространственный заряд возникает

около мембран при запредельном режиме, причем при дальнейшем увеличении падения потенциала область пространственного заряда увеличивается (см. рис. 1). Кроме того, непосредственно вблизи мембран возникает погранслой, обусловленный влиянием граничного условия. Это влияние вне погранслоя практически отсутствует.

) )

Рис. 1. Плотность распределения заряда р = г1£1 + г2 £2 в безразмерном виде при е = 10-3: а) йр= 8, б) йр= 30

Из рисунка видно, что в ядре потока, т.е. в центральной части канала обессоливания электронейтральность раствора сохраняется с большой точностью (р ~ 0).

В области пространственного заряда процесс практически стационарен.

В зависимости от начального значения возможно возникновение и начального по-гранслоя.

Таким образом, структура канала обессоливания в зависимости от величины ёр

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

6.2. Соленоидальность плотности тока

Численные расчеты показывают, что ~ 0 внутри канала, включая области

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

Таким образом, можно считать, что внутри канала ~ 0, т.е. поле плотности тока является соленоидальной. Для объяснения этого интересного результата обратимся к уравнению (2):

дС

---= -&\] , / = 1,2.

дг Л

Умножим каждое уравнение на zi и сложим, тогда:

д( 2,СХ + 12С2)

= —

дг

В области электронейтральности с большой точностью выполняется условие электронейтральности г1С1 + г2С2 = 0, следовательно, с такой же точностью аШ = 0 . С другой стороны в области пространственного заряда из-за стационарности процессов

д(г1С1 + 12С2) Л 7 Л

переноса имеем ----—----= 0 , следовательно, опять аш = 0.

дг

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

6.3. Основные закономерности распределения концентрации ионов

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

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

6.4. Асимптотическое решение краевой задачи при больших числах Пекле

При значительных скоростях прокачки раствора через камеру обессоливания, например, при ¥0 > 10-2 м / с, Н = 10-3 м и Б0 = 10-9 м / с2, число Пекле Ре > 104.

Поделим в уравнении (28) обе части на число Пекле, тогда получим

§ = ЦШ, + гДШ{.%УФ)-, = 1,2. (35)

Это уравнение является сингулярно возмущенным, так как малый параметр Х = 1/Ре является множителем при операторе Лапласа, и система уравнений (29), (35) плохо поддается численному решению. Однако теперь для ее решения можно исполь-

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

Б = ^,0 + ЯБи, +... , (36)

Ф = Ф0 + ЯФХ +... . (37)

Подставляя разложения (36), (37) в уравнения (35) и (29), получим рекуррентную

систему уравнений для 0, 1,... и Ф0, Ф1,... .

Система уравнений для нулевого приближения имеет вид:

дБг,0

■■-агу(101¥) г = 1,2, (38)

дт

ЛФ0 =-(21Б1,0 + 22Б2,0 ) . (39)

Она распадается и может быть последовательно решена. Из уравнения (38) следует, что распределение концентраций в ядре потока определяется в первую очередь конвективным переносом.

Система уравнений для первого приближения имеет вид:

дБ

1 = ВгЛБг,0 + 2 Ааш{Би0^ Ф0)- 1^0)- ), г = 1,2:

дт г "0 +2

ЛФ1 = -(21 Б1,1 + 22Б2,1) .

Она также распадается и может быть последовательно решена.

Разложения (36), (37) вблизи границ должны быть дополнены погранслойными функциями.

6.5. Решение краевой задачи при средних значениях числа Пекле

Если число Пекле 10 < Ре < 100, то число Я уже нельзя считать малым числом, но система уравнений (29), (35) является, тем не менее, «жесткой» системой. Для ее решения можно использовать следующие методы последовательных приближений:

ЛФ( п ) =- (гд( п-1) + 2 2 Б 2 п-1)),

дБ(п) дт

дБ(п)

или

Д ЯЖг(п-1) + ¿Д. Яйу(г(п-1)У Ф(п)) - )ф,(п)Щ) г = 1,2, = ДЯЖг(п-1) + ггДЯйу(г(п-1)У Ф(п-1)) - ) г = 1,2

дт

ЛФ(п) =-(гБп) + 22Б2п)):

начатых с некоторого нулевого приближения (Б1(0), £20), Ф(0))

) ) -

Рис. 2. Графики S1, S2 при t = 100, df= 1б, Є = 10_3

) ) -

Рис. 3. Графики S1, S2 : а), б) - при t = 4б4, dp= 1б, є = 10_4, в), г) - при t = 10000, dp= 25, Є = 10_б

6с 10с 20с

Рис. 4. Формирование начального погранслоя. Графики С1 в размерном виде в разные моменты времени в секундах при = 0,6В

Заключение

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

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

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

1. . . М.: Мир, 1977. 463 с.

2. . ., . .,

. . -

переноса в канале обессоливания электро-диализного аппарата // Вестник Адыгейского государственного университета. Сер. Естественно-математические и технические науки. 2010. Вып. 1 (53). С. 75-89. иКЬ: http://vestnik.adygnet.ru

3. . .

- - //

Краснодар, КубГУ, 1998. 126 с.

4. ., .

конечных элементов: пер. с англ. М.: Мир, 1981. 304 с.

5. ., ., . -

мерные численные методы решения задач с пограничным слоем. М.: Мир, 1983. 254 .

6. . ., . . -

ческие методы в теории сингулярных возмущений. М.: Высш. шк., 1990. 208 с.

1. Newman G. Electrochemical systems. M.: Mir, 1977. 463 p.

2. Mamiy D.K., Shaposhnikova T.L., Urtenov K.M. Mathematical model of heat-and-mass transfer in the channel of desalting of the electrodialysis unit // The Bulletin of the Adyghe State University. Series Natural-Mathematical and Technical Sciences. 2010. Iss. 1 (53). P. 75-89. URL: http://vestnik.adygnet.ru

3. Urtenov M.Kh. Boundary-value problems for the systems of Nernst-Planck-Poisson equations // Krasnodar, KubGU, 1998. 126 p.

4. Norri D., Friz J. de. Introduction to the method of finite elements: transl. from English. M.: Mir, 1981. 304 p.

5. Dulan E., Miller G., Schilders U. Uniform numerical solution methods of tasks with a boundary layer. M.: Mir, 1983. 254 p.

6. Vasilyeva A.B., Butuzov V.F. Asymptotic methods in the theory of singular perturbations. M.: Vyssh. shk., 1990. 208 p.