Scientific journal
Fundamental research
ISSN 1812-7339
"Перечень" ВАК
ИФ РИНЦ = 1,674

2D MODELING OF TRANSPORT ARBITRARY BINARY ELECTROLYTES IN ELECTRO-MEMBRANE SYSTEMS WHEN THE CONDITION OF ELECTRICAL NEUTRALITY

Kovalenko A.V. 1
1 Kuban State University
2660 KB
Annotation. In the simulation of a binary electrolyte transfer in electro-membrane systems, uses a system of equations Nernst-Planck equation. From this system, when the condition of electrical neutrality can obtain the equation of convective diffusion for concentration. Previously, it was shown that the use of convective diffusion equation at exorbitant current densities causing a number of difficulties. In this paper, we propose a new approach to 2D modeling of transport symmetrical binary electrolyte with the same diffusion coefficients of cations and anions in electro-membrane systems under the condition of electric neutrality at exorbitant current densities. The new approach in the use of partial differential equations of the first order to describe the distribution of the concentration of salt ions, rather than convective diffusion equation. In this paper, this approach is generalized to arbitrary binary electrolyte with different diffusion coefficients of cations and anions and show that it allows you to estimate the space charge region at exorbitant current densities.
mathematical modeling
2D modeling
Nernst-Planck-Poisson
convective diffusion
1. Vabishhevich P.N. Chislennye metody reshenija zadach so svobodnoj granicej. M.: MGU, 1987.
2. Kovalenko A.V. 2D modelirovanie perenosa 1:1 jelektrolita v jelektromembrannyh sistemah pri vypolnenii uslovija jelektronejtralnosti // Politematicheskij setevoj jelektronnyj nauchnyj zhurnal Kubanskogo gosudarstvennogo agrarnogo universiteta (Nauchnyj zhurnal KubGAU) Krasnodar: KubGAU, 2015. no. 06(110). http://ej.kubagro.ru/2015/06/pdf/25.pdf, рр. 373–387.
3. Kovalenko A.V. Dekompozicija dvumernoj sistemy uravnenij Nernsta-Planka-Puassona dlja ternarnogo jelektrolita/ A.V. Kovalenko, M.H. Urtenov, A.A. Hromyh // Doklady Akademii nauk. M: MAIK «Nauka/Interperiodika». 2014. T. 458. no. 5. рр. 526–527.
4. Kovalenko A.V. Kraevye zadachi dlja sistemy jelektrodiffuzionnyh uravnenij. Chast 1. Odnomernye zadachi / A.V. Kovalenko, M.H. Urtenov. Germany: LAP LAMBERT Academic Publishing. 2011.
5. Njumen Dzh. Jelektrohimicheskie sistemy. M.: Mir. 1977.
6. Nikonenko V., Kovalenko A., Urtenov M., Pismenskaya N., Han J., Sistat P., Pourcelly G. Desalination at overlimiting currents: State-of-the-art and perspectives // Desalination. USA: ELSEVIER. 2014. no. 342. рр. 85–106
7. Urtenov M.K., Uzdenova A.M., Nikonenko V.V., Pismenskaya N.D., Kovalenko A.V., Vasileva V.I., Sistat P., Pourcelly G. Basic mathematical model of overlimiting transfer enhanced by electroconvection in flow-through electrodialysis membrane cells // Journal of Membrane Science. USA: ELSEVIER. 2013. рр. 190–202.

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

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

1.1. Уравнения

Рассмотрим стационарный перенос бинарного электролита, тогда система уравнений Нернста – Планка и условия электронейтральности имеют следующий безразмерный вид [2]:

kovalenko001.wmf i = 1, 2; (1)

kovalenko002.wmf i = 1, 2; (2)

kovalenko003.wmf (3)

kovalenko004.wmf (4)

1.2. Геометрия

Рассмотрим половину канала обессоливания. Пусть x = 0 соответствует глубине раствора, где концентрация сохраняется постоянной, x = 1 соответствует условной границе раствор/катионообменная мембрана. При допредельных плотностях тока условие электронейтральности (9) выполняется везде за исключением погранслоя возле катионообменной мембраны, которая является областью пространственного заряда, половиной двойного электрического слоя (другая половина двойного электрического слоя расположена в мембране). При запредельных плотностях тока область пространственного заряда увеличивается при увеличении плотности тока и, соответственно, область электронейтральности уменьшается. В связи с этим возникает проблема оценки области электронейтральности при запредельных плотностях тока и определения закономерностей переноса ионов соли в этой области.

1.3. Граничные условия

Ниже приведены общепринятые граничные условия для системы уравнений (1)–(4).

1) В глубине раствора при x = 0:

kovalenko005.wmf kovalenko006.wmf

причем z1C10 + z2C20 = 0.

2) На поверхности КОМ при x = 1: kovalenko007.wmf предполагается также, что КОМ идеально селективная, т.е. поток анионов j2 через катионообменную мембрану равен нулю kovalenko008.wmf.

3) На входе в канал при y = 0:

kovalenko009.wmf kovalenko010.wmf

4) На выходе из канала при y = L:

kovalenko011.wmf i = 1, 2.

Замечание 1. Кроме граничных условий (1)–(4) на концентрацию, обычно задают следующие условия на потенциал φ: kovalenko012.wmf (т.е. прямая x = 0 является эквипотенциальной), kovalenko013.wmf (т.е. поверхность КОМ (x = 1) является эквипотенциальной), на входе в канал задается распределение kovalenko014.wmf, на выходе «мягкое» условие kovalenko015.wmf. При допредельных плотностях тока, величины φ0 и C1m связаны между собой [2]. Указанные выше условия на φ задают потенциостатический режим, однако для развиваемого нами в дальнейшем метода удобно будет переходить к гальваностатическому режиму. Между этими режимами существует однозначная связь. В гальваностатическом режиме задается средняя плотность тока iav по формуле

kovalenko016.wmf

2. Декомпозиция

Покажем, что вместо исходной системы уравнений (1)–(4) для определения неизвестных C1, C2, kovalenko017.wmf, kovalenko018.wmf, kovalenko019.wmf можно получить два уравнения для двух неизвестных функций, после решения которых остальные неизвестные могут быть найдены по формулам или по отдельным независимым уравнениям, то есть можно произвести расщепление (декомпозицию) исходной системы уравнений.

1. Обозначим C1 + C2 = S. Путем непосредственных вычислений проверяется справедливость формул (5), где

kovalenko020.wmf kovalenko021.wmf

kovalenko022.wmf

kovalenko023.wmf

kovalenko024.wmf (5)

kovalenko025.wmf

2. Сложим уравнения (1):

kovalenko026.wmf

С учетом (5) имеем (6), где kovalenko027.wmf. Тогда уравнение (6) перепишем в виде (7).

kovalenko028.wmf (6)

kovalenko029.wmf (7)

3. Сложим уравнения (1), умноженные соответственно на z1 и z2:

kovalenko030.wmf

С учетом (3), (4) и (5) имеем kovalenko031.wmf откуда следует

kovalenko032.wmf (8)

С учетом (8) из (7) для S получаем уравнение

kovalenko033.wmf (9)

где kovalenko034.wmf – коэффициент диффузии электролита и kovalenko035.wmf – общий суммарный поток.

Таким образом, для решения системы уравнений (1–4) достаточно решить следующую систему уравнений (8), (9).

4. Из уравнений kovalenko036.wmf i = 1, 2 следует, что kovalenko037.wmf, kovalenko038.wmf, kovalenko039.wmf, kovalenko040.wmf – соленоидальные вектора. Уравнения (8), (9) зависят от соленоидальных векторов kovalenko041.wmf, kovalenko042.wmf и kovalenko043.wmf.

Уравнения, не зависящие от соленоидальных векторов kovalenko044.wmf, kovalenko045.wmf и kovalenko046.wmf, можно получить, если взять div от обеих частей (8), (9), тогда получим уравнение конвективной диффузии для функции S (10) и уравнение переноса для kovalenko047.wmf: kovalenko048.wmf Это уравнение обычно записывают для потенциала с учетом kovalenko049.wmf, тогда получаем (11).

kovalenko050.wmf (10)

kovalenko051.wmf (11)

Для решения уравнений (10), (11) необходимо ставить дополнительные граничные условия из-за того, что получились дифференциальные уравнения более высокого порядка. Общепринятая практика, наряду с переходом к уравнениям (10), (11), заключается в том, что ставятся граничные условия при x = 1, например условие

kovalenko052.wmf (12)

причем S1m(y) > 0, поскольку S является суммарной концентрацией. Из этого сразу следует, что такая задача моделирует лишь допредельный режим iav < ilim = 2 [2]. Если же положить kovalenko053.wmf, то задача для пары функций kovalenko054.wmf не имеет решение [2].

Дифференциальные уравнения (10), (11) при запредельных плотностях тока можно рассматривать как краевые задачи в области с неизвестной правой границей и необходимо при этом доопределить краевые условия. Пусть нули функции S описываются формулой x = q(y), y ∈ [0, L], т.е. S(q(y), y) = 0 для любого y ∈ [0, L]. Обозначим S решение краевой задачи: kovalenko055.wmf в области kovalenko056.wmf: kovalenko057.wmf S(q(y), y) = 0; kovalenko058.wmf ΔS(q(y), y) = 0; kovalenko059.wmf Таким образом, для S получаем переопределенную краевую задачу с неизвестной границей. Решение таких задач является достаточно сложной проблемой [1].

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

3. Вспомогательное уравнение kovalenko060.wmf

Рассмотрим решение уравнения kovalenko061.wmf в 2D односвязной области Ω при известном гладком векторе kovalenko062.wmf.

Определим двумерный аналог оператора rot, – оператор

kovalenko063.wmf

Необходимым и достаточным условием разрешимости уравнения kovalenko064.wmf является условие kovalenko065.wmf или kovalenko066.wmf, при этом решение выражается формулой [2]:

kovalenko067.wmf

Аналогично можно получить формулу

kovalenko068.wmf

Формулу для S можно записать симметрично:

kovalenko069.wmf

где S0 = S(0, 0) – некоторая постоянная.

Таким образом, для однозначного решения системы уравнений kovalenko070.wmf достаточно условия S(0, 0) = S0.

4. Непосредственное решение системы уравнений (9)

Так как kovalenko071.wmf некоторый соленоидальный вектор, то существует функция γ, что

kovalenko072.wmf kovalenko073.wmf

4.1. Условие разрешимости

Обозначим kovalenko074.wmf Согласно п. 3 условие разрешимости уравнения (9) – kovalenko075.wmf откуда получаем уравнение для функции γ:

kovalenko076.wmf

К этому уравнению добавляем уравнение (9) и получаем систему из двух уравнений с двумя неизвестными γ и S:

kovalenko077.wmf (13)

kovalenko078.wmf (14)

4.2. Граничные условия

Из граничных условий для C1 и C2 для функции S получаются граничные условия

kovalenko079.wmf kovalenko080.wmf

kovalenko081.wmf kovalenko082.wmf (15)

Наряду с граничными условиями (15) на функцию S необходимы граничные условия для функции γ. Как следует из п. 3, условия (15) избыточны для нахождения решения уравнения (13) и должны выводится граничные условия для функции γ.

С учетом (13) и (15) для вектора kovalenko083.wmf получаем

kovalenko084.wmf

kovalenko085.wmf

kovalenko086.wmf

kovalenko087.wmf

kovalenko088.wmf

kovalenko089.wmf

kovalenko090.wmf

kovalenko091.wmf

В силу условий (15) имеем kovalenko092.wmf kovalenko093.wmf Откуда с учетом kovalenko094.wmf kovalenko095.wmf получаем граничные условия для функции γ:

kovalenko096.wmf

kovalenko097.wmf (16)

Как видно, из уравнения (14) и граничных условий (16) функция γ определяется с точностью до постоянной величины.

Определим конкретное значение функции γ и свяжем его с заданной плотностью тока. Для этого воспользуемся соотношением kovalenko098.wmf.

В силу идеальной селективности катионообменной мембраны j1,1 >> j2,1 (в одномерном случае вообще j2 = 0 и j1 = iav), поэтому J1,1 ≈ I1,1. Следовательно,

kovalenko099.wmf

где kovalenko100.wmf

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

kovalenko101.wmf

или

γ(x, L) – γ(x, 0) = –kgLiav.

Поскольку функция γ определяется с точностью до постоянной величины, то можно положить, например, γ(x, L) = 0, тогда γ(x, 0) = –kgLiav. Заметим что если D1 = D2, то kg = 1.

Итак, для нахождения функций γ и S имеем краевую задачу

kovalenko102.wmf (17)

kovalenko103.wmf (18)

kovalenko104.wmf kovalenko105.wmf kovalenko106.wmf (19)

γ(x, L) = 0; γ(x, 0) = –kgLiav; (20)

kovalenko107.wmf

kovalenko108.wmf (21)

5. Решение краевой задачи (17)–(21)

5.1. Метод последовательных приближений

Для решения краевой задачи (17)–(21) можно использовать следующий метод последовательных приближений:

1) Возьмем некоторый соленоидальный вектор kovalenko109.wmf.

2) Найдем решение S(0), уравнения kovalenko110.wmf, удовлетворяющее граничным условиям при x = 0 и y = 0.

3) Найдем нули функции S(0), т.е. такую функцию x = q(0)(y), y ∈ [0, L], что S(0)(q(0)(y), y) = 0, y ∈ [0, L].

4) Найдем решение γ(1) уравнения

kovalenko111.wmf

в области

kovalenko112.wmf,

удовлетворяющее соответствующим граничным условиям (20), (21).

5) Найдем kovalenko113.wmf и kovalenko114.wmf по формулам

kovalenko115.wmf и kovalenko116.wmf.

6) Найдем решение S(1) уравнения kovalenko117.wmf, удовлетворяющее соответствующим граничным условиям.

7) Проверим условие сходимости kovalenko118.wmf где δ – заданная точность. Если условие сходимости выполняется, S(1) принимаем за решения, иначе полагаем kovalenko119.wmf и идем к п. 2.

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

Если Pe << 1, т.е. скорость прокачки раствора небольшая, то функции γ и S можно разложить по степеням Pe:

kovalenko120.wmf kovalenko121.wmf

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

Для γ(0), S(0) получаем краевую задачу:

DΔγ(0) = 0; (22)

kovalenko122.wmf (23)

kovalenko123.wmf kovalenko124.wmf

kovalenko125.wmf (24)

γ(0)(x, L) = 0; γ(0) (x, 0) = –kgLiav; (25)

kovalenko126.wmf

kovalenko127.wmf (26)

Для γ(1), S(1) получаем краевую задачу:

kovalenko128.wmf (27)

kovalenko129.wmf (28)

kovalenko130.wmf kovalenko131.wmf

kovalenko132.wmf (29)

γ(1)(x, L) = 0; γ(1)(x, L) = 0; (30)

kovalenko133.wmf

kovalenko134.wmf (31)

1) Решение краевой задачи для начального приближения γ(0), S(0)

Рассмотрим два разных метода решения краевой задачи.

Первый метод. Сначала решаем краевую задачу для уравнения kovalenko135.wmf аналитически и находим kovalenko136.wmf, где Γ интегральный оператор (см. п. 3). Далее находим решение уравнения DΔγ = 0 в области Ω, с заранее неизвестной границей, определяемой дополнительным условием kovalenko137.wmf и граничным условием kovalenko138.wmf на этой границе.

Второй метод. Если для решения использовать метод последовательных приближений п. 5.1, то получаем, что сначала задается некоторый соленоидальный вектор kovalenko139.wmf. Затем находим kovalenko140.wmf и определяем градиент kovalenko141.wmf и нули x = q(0)(y) функции S(0). Далее решаем уравнение DΔγ = 0 в области kovalenko142.wmf, удовлетворяющее соответствующим граничным условиям, и далее по алгоритму п. 5.1.

2) Решение краевой задачи для первого приближения γ(1), S(1)

Краевая задача для γ(1), S(1) отличается от краевой задачи для γ(0), S(0) лишь тем, что уравнение для γ(1) неоднородное, а для S(1) граничные условия однородные. Поэтому краевая задача для γ(1), S(1) решается аналогично краевой задаче для γ(0), S(0).

6. Случай течения Пуазейля

6.1. Упрощение краевой задачи для течения Пуазейля

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

kovalenko143.wmf

где kovalenko144.wmf

В этом случае из (17)–(21) имеем

kovalenko145.wmf (32)

kovalenko146.wmf (33)

kovalenko147.wmf (34)

6.2. Сведение краевой задачи (32)–(34) к краевой задаче для интегро-дифференциального уравнения

Из уравнения (32) с учетом граничного условия S(0, y) = 2 для любого y ∈ [0, L] получаем

kovalenko148.wmf (35)

С другой стороны, применяя общую формулу п. 3, с учетом kovalenko149.wmf и kovalenko150.wmf получаем

kovalenko151.wmf

или

kovalenko152.wmf

где α постоянная, не зависящая от x, y. Эта формула согласуется с формулой (35) в силу граничных условий

α = 2 и kovalenko153.wmf

Действительно, второе соотношение, как легко видеть после дифференцирования по y, эквивалентно граничному условию

kovalenko154.wmf

Подставляя (35) в уравнение (34), получаем для функции γ краевую задачу для интегро-дифференциального уравнения:

kovalenko155.wmf (36)

с соответствующими граничными условиями (20), (21).

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

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

kovalenko156.wmf k = 0, 1, ... (37)

в области kovalenko157.wmf, с граничными условиями (19)–(21).

Пусть kovalenko158.wmf некоторый соленоидальный вектор, S(0) рассчитывается по формуле kovalenko159.wmf. Находим x = q(0)(y) – нули функции S(0). Решая краевую задачу для выражения (37) при k = 0, находим функции γ(1), kovalenko160.wmf, а затем функцию S(1) и т.д.

6.4. Пример вычисления последовательных приближений

Рассмотрим для простоты симметричный 1:1 электролит с D1 = D2 при запредельных токовых режимах iav > ilim = 2. Тогда уравнение и краевые условия упростятся, так как тогда D = 1 и kg = 1, причем

kovalenko161.wmf kovalenko162.wmf kovalenko163.wmf

1) Положим kovalenko164.wmf Выберем a и b так, чтобы на входе в канал был предельный режим и удовлетворялись граничные условия (20), т.е. kovalenko165.wmf и kovalenko166.wmf или kovalenko167.wmf. Таким образом kovalenko168.wmf

2) Вычислим S(0) по формуле (35).

3) Найдем область

kovalenko169.wmf.

Для этого найдем нули функции

kovalenko170.wmf

Таким образом

kovalenko171.wmf

Область

kovalenko172.wmf

является начальным приближением области электронейтральности при заданной запредельной плотности тока iav (iav > ilim = 2).

4) Решим уравнение (37) при k = 0

kovalenko173.wmf

в области kovalenko174.wmf, с краевыми условиями

γ(1)(x, L) = 0; γ(1)(x, 0) = Liav; kovalenko175.wmf kovalenko176.wmf

pic_18.tif

Рис. 1. График функции γ(1)

5) Находим kovalenko177.wmf, а затем функцию S(1).

6) Находим нули функции S(1), т.е. функцию x = q(1)(y) и определяем очередное приближение к области электронейтральности kovalenko178.wmf и т.д.

7) Результаты численных экспериментов.

Ниже приведены результаты численного решения при iav = 2ilim = 4, Pe = 0,1.

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

Сопоставление графиков kovalenko179.wmf и kovalenko180.wmf показывает, что следующая итерация приводит к коррекции начального линейного приближения kovalenko181.wmf зависящего только от одной переменной y и замене его на нелинейную функцию kovalenko182.wmf от обоих аргументов. Кроме того, видно, что в конце канала kovalenko183.wmf, т.е. ток течет практически поперек канала. Это связано с тем, что число Pe = 0,1 небольшое.

pic_19.tif

a                                 б                                  в

Рис. 2. Графики функций: a – kovalenko185.wmf; б – kovalenko186.wmf; в – kovalenko187.wmf

pic_20.tif

a                                                        б

Рис. 3. Графики: a – S(0); б – S(1)

pic_21.tif

а                                                         б

Рис. 4. Линии уровня: a – S(0); б – S(1)

Поскольку в данном случае kovalenko184.wmf, то концентрация катионов и анионов может быть проанализирована по виду функции S. На а и б видно наличие «плато» в ядре потока и снижение концентрации до нуля на правой границе. Наличие отрицательных значений у функции S(1) в отличие от начального приближения S(0) говорит о том, что следующее приближение Ω1 к области электронейтральности отличается от начального приближения Ω0. Более ясно это видно на рис. 4, б, – в верхнем правом углу области часть области электронейтральности в приближении Ω0 в следующем приближении Ω1 уже исключается (не заштрихованная область в криволинейной трапеции Ω0) и правая граница области электронейтральности уточняется.

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

Замечание 2. Если V2(x) = const = V0, уравнение (33) не зависит от S и имеет вид kovalenko188.wmf или kovalenko189.wmf. Решая это уравнение с соответствующими условиями, находим kovalenko190.wmf и по формуле (34) находим функцию S.

Замечание 3. Метод 2D моделирования, развитый здесь для бинарного электролита, несложно обобщается на случай произвольного числа ионов.

Заключение

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

Полученные в работе результаты могут быть использованы в нанотехнологиях, мембранной электрохимии и микрофлюидике при моделировании и исследовании переноса бинарного при запредельных плотностях тока [3, 6–7].

Работа выполнена при финансовой поддержке РФФИ и администрации Краснодарского края грантов: № 13-08-93105-НЦНИЛ_а и 13-08-96519 р_юг_а.

Рецензенты:

Халафян А.А., д.т.н., доцент, профессор кафедры прикладной математики, ФГБОУ ВПО «Кубанский государственный университет», г. Краснодар;

Павлова А.В., д.ф.-м.н., доцент, профессор кафедры математического моделирования, ФГБОУ ВПО «Кубанский государственный университет», г. Краснодар.