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

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

В то же время в практике электродиализа, электрохимической характеризации мембран (хронопотенциометрия, импедансометрия и др.) часто используется гальваностатический режим, при котором на межфазной границе поддерживается постоянная средняя плотность тока. Об этом режиме собран огромный объем экспериментальных данных, которые необходимо интерпретировать [8–10] Ошибка! Закладка не определена. . Исследования в области математического моделирования гальваностатического режима ведутся по нескольким направлениям.

Первое направление — метод обратной задачи. Как ясно из названия, речь идет о решении обратной задачи: для заданной плотности тока на межфазной границе «раствор — мембрана» находится соответствующий скачок потенциала, а далее рассматривается задача для потенциостатического режима [11]. Низкая эффективность данного метода обусловлена тем, что его реализация требует многократного решения задачи в потен-циостатическом режиме для одного заданного значения плотности тока.

Второе направление — метод декомпозиции. При этом система уравнений Нернста — Планка — Пуассона заменяется системой декомпозиционных уравнений [12–16]. Предположение о квазиравномерном распределении заряда позволяет получить модель для гальваностатического режима в приближении закона Ома [17– 20].

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

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

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

Материалы и методы. Под мембранной системой подразумевается канал обессоливания электродиализного аппарата (ЭДА), образованный анионообменной (АОМ) и катионообменной (КОМ) мембранами. Через него со средней скоростью V 0 прокачивается раствор бинарного электролита. На рис. 1 х — нормальная к поверхности мембраны координата, изменяющаяся от 0 (граница с АОМ) до h (граница с КОМ); у — тангенциальная к поверхности мембраны координата, изменяющаяся от 0 (вход в канал) до l (выход из канала).

Информатика, вычислительная техника и управление

Рис. 1. Схема канала обессоливания ЭДА. Показаны концентрационные профили катионов С ₁ (сплошная линия) и ^— анионов С ₂ (пунктирная линия). Скорость вынужденного течения V

Система уравнений. Рассмотрим двумерный стационарный случай системы уравнений, описывающей перенос бинарного электролита в отсутствие химических реакций [24]:

-    F    _ „_ -

ji =--ztDi Ci Vф - D VCt + СV, i = 1, 2 , i RT i i ii i i

—

- div jt = 0, i = 1,2,(2)

EqEr Аф = —F (z 1C1 + z 2 C 2),

I = F (z 1 — + z 2 -2).(4)

—►

Здесь j i , D i , z i и C i — соответственно поток, коэффициент диффузии, зарядовое число и молярная концентрация i -го иона; ф — электрический потенциал; V — скорость потока раствора электролита; е ₀ — электрическая постоянная; е r — относительная диэлектрическая проницаемость раствора электролита (предполагается постоянной); I — плотность тока; F — постоянная Фарадея; R — газовая постоянная; T — абсолютная темпе-

ратура; j 1 , j 2 , I , ф , Q , С 2 — неизвестные функции координат x и у .

Уравнения Нернста — Планка (1) описывают поток ионов, обусловленный миграцией в электрическом поле, диффузией и конвекцией; (2) — уравнение материального баланса в стационарном случае; (3) — уравнение Пуассона для потенциала электрического поля; (4) — плотность тока в растворе электролита. Будем считать, что распределение скорости в канале соответствует течению Пуазейля [24]:

V x = 0, V y = 6V 0 x 1 1 - x | . (5) h V h )

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

Граничные условия для моделирования потенциостатического режима. Будем считать поверхности ионообменных мембран эквипотенциальными. Система (1)–(4) включает потенциал электрического поля только в форме производных по пространственной координате. В этом случае существенным является только скачок потенциала Аф = ф ( h , у ) -ф (0, у ) , где Аф — известная функция, поэтому положим, например:

ф ( ^0, У ) = 0 . ⁽⁶⁾

Тогда

Аф = const .

Условия (6) и (7) определяют потенциостатический режим. Остальные граничные условия приводятся ниже.

На межфазной границе «АОМ — раствор» ( x = 0) концентрация коионов (катионов) определяется из условия непрерывности их потока у границы «мембрана — раствор» с учетом селективных свойств АОМ [16]:

fa С,   F дф)        ( 1 - T 2 ) I x (0, y )

I-- ¹-- z С       1 ( 0, У ) =---------------------

( д x  RT ¹¹ д x J           z 1 FD 1

где Ti (i = 1, 2) — эффективные числа переноса противоионов в мембране (КОМ и АОМ соответственно); Ti — числа, близкие к 1, причем для идеально селективной мембраны Ti = 1, и условие (8) превращается в условие непроницаемости мембраны для коионов.

Концентрация противоионов (анионов) зависит от обменной емкости КОМ, что можно задать в виде:

С ( 0 ,У ) = С 2 m = NSr

Здесь постоянная N _a показывает, во сколько раз эта концентрация отличается от концентрации в объеме раствора [25].

На межфазной границе «раствор — КОМ» ( x = h ) для концентраций ионов приняты условия, аналогичные условиям на границе «АОМ — раствор» ( x = 0):

С ( h , У ) = С m = N c С 0 ,

fa с 2    F       дф)        ( 1 - T 1 ) I x ( h , У )

I + z? С?       l ( h , У ) = •

( д x  RT ²² д x J          z 2 FD 2

На входе в канал (у = 0) принято равномерное распределение концентрации ионов:

C i ( x ,0 ) = С ₀, i = 1,2.

Условие для потенциала электрического поля получено из уравнений (1) и (4) с учетом отсутствия протекания тока через вход I x ( x , 0, t ) = 0:

дф ( x ,0) _           RT        f _n д C 1 ( x ,0)     _ д C 2 ( x ,0) )

---a =-- 2------- 5 --------1 z 1 D 1   ⁺ z 2 ^D 2 “  I.

дУ      F (z 12 D1 + z 2 D 2) C 0 (       дУ             9yJ

На выходе из канала ( у = l ) ионы свободно выносятся потоком раствора:

--

(-, j ) = (-,- — zDa^ Уф- Div Ci + CiV) = С^, i = 1,2.(14)

Система (14) означает также, что сумма диффузионной и миграционной тангенциальных составляю- щих потоков анионов и катионов равна 0:

д C_i ^д У

_ F C дф RT i¹ i д У

Для потенциала принято «мягкое» условие, означающее отсутствие резких изменений потенциала на

выходе из канала:

|ф ( x , l ) = 0. ^д У

Краевая задача, включающая уравнения (1)-(4) и краевые условия (6)-(16), моделирует потенциостати-ческий режим, причем его определяющее ключевое условие — (7).

Граничные условия для моделирования гальваностатического режима. При моделировании гальвано-статического режима условие (6) остается, а (7) необходимо заменить условием, связанным с заданным значением средней плотности тока i_av на межфазной границе «раствор — КОМ» ( x = h ). Для вывода такого граничного условия подставим соотношения (1) в (4) и выразим градиент потенциала электрического поля:

Уф =-- гч -------2----- ( I + F ( z D 1 V С 1 + z 2 D 2 V С 2 ) - F ( z ₁ C ₁ + z 2 С 2 ) Т ^- ).                (17)

F ²( z ² D 1 C 1 + z 2 ² D 2 C 2 )

Полагая в (17) x = h , получаем соотношение, связывающее градиент электрического потенциала с заданным значением плотности тока на границе, то есть граничное условие на межфазной границе «раствор — КОМ»:

т

RT

—

f            дС,

Ix^{+ Fz} 1^D1^ ^{+ Fz} 2 D ax

^{д С} 2 )

F²

(

z 1²D1С1 + z 2²D 2 С 2

' 2

ax^- (h , У )■

При этом плотность тока Ix должна удовлетворять условию:

Информатика, вычислительная техника и управление

1    \

у] Ix (М) dy = iav l 0

Математическая модель гальваностатического режима состоит из системы уравнений (1)-(4). Граничные условия (18) и (19) заменяют условие (7). Остальные граничные условия совпадают с условиями для по-тенциостатического режима.

Скачок потенциала в гальваностатическом режиме является вычисляемой величиной.

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

В стационарном случае плотность тока I является соленоидальным вектором. Действительно, если умножить (2) на zi и сложить, то divI = 0. Следовательно существует такая функция п, что дл = I 0o = _t дx    y,    ду     x .

С использованием функции п условия (18) и (19) переписываются следующим образом:

£( h,y ) = дx

—

RT

F^

—^ + Fz_xD_x С + FzD С

1    1                    2    2

ду       дx        дx

(h,y),

1 r                  1 г дп

-jIx^(h,y) dy = ^—-J |^^(h,y) dy =

Iо                1о дУ

^{— 1} (П ^(h,l) ^— П ^(h,⁰⁾⁾= iav .

Уравнение (22) можно переписать:

n⁽hl)^—n^(h,⁰⁾= ^—U.

Чтобы замкнуть систему выражений, необходимо получить уравнение для функции п. С этой целью, как и в работах [15, 16], введем в рассмотрение линейный дифференциальный оператор, который является функцией вихря (ротором) в двумерном случае, для произвольного двумерного вектора W :

—— r (W)

Г^дWy  д Wx )

—

I  дx     ду

X               У

.

Несложно проверить, что:

• 1)    r(Vи) = 0 для любой гладкой функции и ;

  ——

  
  

  ——

  
  

  ——

  
  

  ^—

  
  

  2) r(uW) = (Vи, W) 1 + ur(W) для любой гладкой функции и и любого гладкого вектора W .

  —,   ди     ди                                                                   —

  Здесь (Vи, W)1 = —W_y — —Wx — кососимметричное скалярное произведение векторов Vи и W , причем дx ^у  ду

  (а, —)1 = 0 для любого вектора a .

  Применяя (24) к уравнению для плотности тока (4), получаем:

  
  

  ——

  
  

  ——

  
  

  ^—

  
  
  
  

  ^{r (}I ) = ^Fz 1^{r (j}1) ^{+ Fz} 2 ^{r (}j 2 )•

  Используя формулу потоков (1), получим соотношение:

  —      F                                -

  r(Ji) = — ^, zD)r-(CiV^) — Dir(V Ci) + r(CiV),     i = 1,2 .

  RT

  Отсюда с учетом свойств оператора r :

  r (h ) = — FZiDi(V Ci, Vф)1 +(V Ci, V^—)1 + Cr (V^—), i = 1,2 . RT

  
  
  
  
  
  
  
  

  —   дIy  дI_x

  Учитывая (27) и r(I) =---= Ап , уравнение (25) можно записать:

  дx   ду

  
  

  An =

  
  

  -—Лz_v²D_v^5-1 + z₂²D₂

  rt I ¹¹5x   ²¹

  
  

  5c2)5ф ^f₂   5-12   5c2)дф)

  Zf D+ z. D->I + 5x J dy  (     dy       dy J 9x J

  
  

  ( Sc, Sc,) I Sc, Sc, )

  +F z -1 + z,-²V - F z -1 + z,-²V + F (. ( ¹5x   ²5x J ^y ( ¹5y    ²dy J x ^v

  
  
  
  

Из (28) и (20) следует, что функция n определяется с точностью до константы, поэтому можно допу стить:

П( h ,0) = 0.                                                    (29)

Тогда из (23) получаем

n(h,l) = -iavl .                                                      (30)

Условия (29) и (30) — краевые для функции n ■

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

Характерные величины, описывающие задачу. При моделировании процессов массопереноса в камере обессоливания ЭДА принимается ряд характерных значений:

• —    для пространственных координат — межмембранное расстояние h ;

• —    для концентраций ионов — объемная концентрация электролита С0 ;

• —    для скорости — средняя скорость вынужденного течения V_} ;

— для коэффициентов диффузии — коэффициент диффузии электролита D = D 1 D₂ (z 1 -z₂)/(D 1 z 1-D₂z₂);

— для электрического потенциала — тепловой потенциал ф₀= RT/F ;

— для плотности тока — величина i₀ = FDC₀/h (аналог предельной плотности диффузионного тока);

— для потока ионов — диффузионный поток j₀= DC₀/h.

Формулы перехода. Переведем уравнения в безразмерную форму с помощью следующих соотношений (индексом (и) обозначены безразмерные варианты величин):

x(u) = x   y(u) = y	1(u) = L   v(u)	—— V  ^(u)  УУ   12 , С1              , i1,2, V 0          С 0
hh	h ’
(u) _ ф     —(u) _ 1 ф          ,	—   (u)      n	,  Fi u) = — ji, i = 1,2,  Diu) = j 0
	I , П
ф0           i 0	FDC 0

Di

. D

Система уравнений в безразмерной форме имеет вид (индекс (и) для упрощения записи опущен):

———*

ji = - ziDiCi Уф-Di V Ci + PeCiV, i = 1,2,(32)

—*

- div ji = 0, i = 1, 2,(33)

еАф = -(z 1С1 + z 2 С2 ) ,

z,"D *1 + z*D & И 11^    2  2

5x       5x J 5y

—

z,"D * + z2D & 1^*1+ 11^    2  2

5y       5y J 5x J

f 5-     5c₂ 1

+Pe z. —¹+ z, — V

I ¹5x    ²5x J ^y

—

Pe I z_x^5-1 + z₂ —² IV + Pe (.

1              2             x

( 5y     5y J

,

I = ^z 1 j1 ^{+ z} 2 j₂.

Система уравнений (29)-(35) содержит два безразмерных числа: Pe = V3h/D

число Пекле и

е = еrе0RT(СС0h²F²). Физический смысл параметра е состоит в том, что это удвоенный квадрат безразмерной дебаевской длины lD : е = 2(lD/h)²[5]. Оценка величин параметров показывает, что при естественных для электродиализа условиях число Пекле имеет порядок 10²-10⁶, число е имеет порядок 10^-13- 10^-7 , то есть

Информатика, вычислительная техника и управление

может считаться малым параметром.

Для удобства численного решения преобразуем систему уравнений, подставив плотность потока (32) в уравнения (33) и (36):

I = £z, (- z,D, C, Vv - D, VС, + PeC,V).                              (38)

г=1

Таким образом, система уравнений содержит следующие неизвестные функции x, у : С1, С2,V, I_x, Iy . Поля концентраций С1, С2 и потенциала ф определяются решением уравнений (37), (34) соответственно. Компоненты плотности тока Ix, Iy вычисляются с помощью (38). Распределение скорости (5) в безразмерной форме:

V_x= о, V_y = 6 x (1 - x).

Присутствие малого параметра с в уравнении Пуассона (34) означает, что краевая задача является сингулярно возмущенной. Это значительно усложняет ее численное решение, поскольку такие задачи относятся к жестким [26]. Потенциал электрического поля ф и концентрации ионов С1, С2 изменяются очень быстро в узком пограничном слое, толщина которого равна длине Дебая Id [5]. Для решения этой проблемы целесообраз- но уплотнить вычислительную сетку в пограничном слое и использовать специальные методы решения жест-

ких задач [26].

Граничные условия в безразмерной форме. На межфазной границе «АОМ — раствор» (x = 0): (де,        дф)       (1 - т₂) Ix (0, у)

I+ z С       1(0, у) =,

V дx    ¹¹lx 7           z1 D1

С2 (0, у ) = Na, ф^(0,у )= ° |П(0,у ) = 0. дx

На границе «раствор — КОМ» (x = 1):

С1 (1, у ) = Nс, дС2        дф)       (1 - T)Ix (1, у)

--+ z₂ С? — 1(1, у ) =---------------, дx    ²²Ix J^{v  7}     z₂D₂

l^ )= дx

—

( дп _п дС1     _п дС2 )

- -  + z1 D1     + z 2 D 2  , ду        дx         дx

z2 D₁С1 + z 2 D 2 С 2

(1, у)

На входе в канал (у = 0):

дф(x,0) _ бу

—

V

д^у )=0. дx

С,(х ,0) = 1, 1        (

, = 1,2,

^z1 ^D1^{+ z} 2 ^D 2

дС       дС₂ )  _

^{+ z}2 ^D2^^(x,0), ду         ду 7

На выходе из канала (у = l):

( дС,

--—

бу

z^^|ф |(х, l)= 0,, = 1,2, бу у

^lv(x, l )= 0,                                                    (52)

бу

n(x,l) = -,_avl, ,av = const.                                              (53)

После численного расчета системы (34), (35), (37)-(53) скачок потенциала Аф в канале обессоливания определяется по формуле:

1 ^l

Аф = -|ф(1 ,У)dy .

¹ 0

Численное решение найдено методом конечных элементов с помощью пакета Comsol Multiphysics на неравномерной вычислительной сетке (плотность элементов сетки увеличена у границ «раствор — мембрана») [27].

Результаты исследования. Вычисления проведены для 8 = 1,9 -10 9, Pe = 2355 , что соответствует следующим значениям параметров системы:

— входная концентрация раствора электролита NaCl C₀ = 0,1 моль/м³;

— температура T = 298 K;

— коэффициенты диффузии катионов и анионов соответственно D 1 = 1,33*10^-9 м²/с, D₂ = 2,05*10^-9 м²/с;

— числа переноса противоионов в мембранах T1 = 0,972, T₂ = 1;

— зарядовые числа ионов z 1 = 1, z₂ = -1;

• —    отношение концентрации противоионов на границе с мембранами к ее значению на входе в канал N_c = Na = 1; — h = 10^-3 м — ширина канала;

• —    l = 2^ 10^-3 м — длина канала;

• —    V0 = 3,8-10^-3 м/с — скорость прокачки раствора.

На рис. 2 приведены поля концентраций C 1 и С₂, потенциала ф и функции ц, рассчитанные при плотности тока iav = 1,5ilim, где i_lim — это предельная плотность тока, определенная по формуле Левека (в безразмерной форме) [28]:

—

0, 2 .

Здесь 11 = 0,395 — число переноса катионов в растворе [9].

/ , 9 X1 / 3

1, 471       ° I

I ID J i = -*-lim

^{T —}t1

Информатика, вычислительная техника и управление

Рис. 2. Поля концентрации катионов C 1 (а) и анионов С₂ (б), потенциала ф (в) и функции ц (г). Расчет по модели для гальваностатического режима при плотности тока i_av = 1, 5 ilim

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

Рис. 3. Результаты численного расчета по модели для гальваностатического (сплошные линии) и потенциостатического (пунктирные линии) режимов при плотности тока i_av = 1, 5i_Um в сечениях у = 0,11 (1), у = 0,41 (2), у = 0,91 (3): концентрационные профили C 1 и С₂ (а); увеличение рис. 3, а (б); плотность пространственного заряда р = z 1 C + z2C2 (в)

Вольтамперная характеристика (ВАХ, кривая 1 на рис. 4) рассчитана на основе гальваностатической модели. При этом задаваемая плотность тока изменялась: i_av= 0; 0,005i_iim; ...; 1,5i_iim.

Рис. 4. ВАХ, рассчитанные по моделям для гальваностатического (1) и потенциостатического (2) режимов; предельная плотность тока iiim (3) определена по формуле (55)

Предельная плотность тока определена по точке пересечения касательных к ВАХ в начальной части и на наклонном плато [5]. Здесь она совпадает со значением ilim, оцененным по приближенной формуле Левека (55). Также ВАХ была рассчитана на основе потенциостатической модели при изменении скачка потенциала ∆φ = 0; 0,4; …; 40 (кривая 2 на рис. 4). Как видно из рис. 4, ВАХ 1 и 2 совпадают. Таким образом, существует однозначное соответствие каждому скачку потенциала некоторой плотности тока, и наоборот.

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

Вычислительные затраты предлагаемой модели оценивались следующим образом. Фиксировалось время, затраченное на решение краевых задач для гальваностатического и потенциостатического режимов. При этом значения задаваемой точности вычислений, параметры системы были одинаковыми, а а_т = 1, 51m . Затем показатели времени сравнивались. Таким образом выяснилось, что расчет по гальваностатической модели требует в 1,6 раз больше времени. Это связано с тем, что:

• —    гальваностатическая модель содержит дополнительное уравнение для определения распределения плотности тока;

• —    для потенциала на границе x = 1 установлено условие второго рода (46).

Комплексный потенциал электромембранной системы для гальваностатического режима. В работе [25] предложено обобщение импеданса электрохимической системы с помощью функции п, введенной в условиях электронейтральности. Аналогичное обобщение приемлемо и в данном случае. Отметим, что при использовании электрохимического импеданса объект рассматривается только как «черный ящик», и его внутренние свойства определяются косвенно. Полученные выше результаты позволяют ввести понятие комплексного потенциала электромембранной системы: P = ф +1 • п . Комплексный потенциал — это функция координаты любой точки внутри объекта, поэтому, в отличие от электрохимического импеданса, он позволяет исследовать внутренние свойства объекта.

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

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