CC BY
45
УДК 51-73:54.128
ЧИСЛЕННОЕ РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ МОДЕЛИ ПЕРЕНОСА БИНАРНОГО ЭЛЕКТРОЛИТА В ПРИБЛИЖЕНИИ ЗАКОНА ОМА
Коваленко Анна Владимировна
к.э.н., доцент кафедры прикладной математики
Уртенов Махамет Хусеевич
д.ф.-м.н., профессор кафедры прикладной
математики
Кубанский государственный университет, Краснодар, Россия
Чубырь Наталья Олеговна
старший преподаватель кафедры прикладной
математики
Хромых Анна Александровна
старший преподаватель кафедры прикладной
математики
Кубанский государственный технологический университет, Краснодар, Россия
Узденова Аминат Магометовна аспирантка кафедры математического анализа
Карачаево-Черкесский государственный университет, Карачаевск, Россия
Барсукова Виктория Юрьевна
к.ф.-м.н., доцент кафедры дифференциальных и
интегральных уравнений
Кубанский государственный университет,
Краснодар, Россия
Коваленко Анна Владимировна
к. э.н., доцент кафедры прикладной математики
Статья посвящена численному решению краевой задачи модели переноса бинарного электролита в мембранных системах в приближении закона Ома. Предлагаются различные численные методы решения. Установлены основные закономерности переноса
Ключевые слова: МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, СИСТЕМА УРАВНЕНИЙ НЕРНСТА - ПЛАНКА И ПУАССОНА, ЧИСЛЕННЫЕ МЕТОДЫ
UDC 51-71:54.128
NUMERICAL SOLUTION TRANSPORT MODEL BINARY ELECTROLYTE IN APPROXIMATION OF OHM'S LAW
Kovalenko Anna Vladimirovna Cand.Econ.Sci., assistant professor
Urtenov Mahamet Khuseevich Dr.Sci.Phys.-Math., professor
Kuban State University, Krasnodar, Russia
Chubyr Natalia Olegovna
senior teacher of chair of applied mathematics
Khromikh Anna Aleksandrovna
senior teacher of chair of applied mathematics
Kuban State Technological University, Krasnodar, Russia
Uzdenova Aminat Machametovna
postgraduate student of faculty of the mathematical
analysis
Karachaevo-Circassian state university, Karachaevsk, Russia
Barsukova Victoria Jurevna Cand.Phys.-Math.Sci., assistant professor
Kuban State University, Krasnodar, Russia
Kovalenko Anna Vladimirovna Cand.Econ.Sci., assistant professor
The article is devoted to the numerical solution of boundary value problem of the binaryelectrolyte model of transport in membrane systems in the approximation of Ohm's law. Different numerical methods are offered. The main regularities of transfer are established
Keywords: MATHEMATICAL MODELING, SYSTEM NERNST - PLANCK AND POISSON EQUATIONS, NUMERICAL METHODS
Введение
В работе [1] нами было найдено асимптотическое представление решения краевой задачи модели переноса бинарного электролита в мембранных системах в приближении закона Ома. Показано, что канал обессоливания (х,у)є [0,H]х[0,Ь], при каждом фиксированном ґ разбивается на различные области [0,Н]х[0,Ь] = и1 ии2 ииз, которые определяются знаком функции обобщенной концентрации £(ґ, х, у) :
электронейтральности и1 (£(ґ, х, у) > 0), пространственного заряда и2 (£(ґ, х, у) < 0), промежуточного слоя из (окрестность нулей функции £(ґ, х, у) ) и асимптотика решения имеет разный вид в различных областях. При «мягких» (допредельных) токовых режимах области и2 и и3, значительно меньше области электронейтральности и1. При жестких (запредельных) режимах область пространственного заряда и2 является уже макроскопической и сопоставимой с областью электронейтральности и1, но промежуточный слой остаются значительно меньшей, чем эти области. Система уравнений с частными производными, описывающая асимптотическое представление решения в каждой из областей меняет свой тип и в них встречаются члены, зависящие от знака искомой функции. Данная статья посвящена разработке метода численного решения таких краевых задач.
1 Краевая задача асимптотического представления решения
1.1 Система уравнений
Асимптотическое представление решения краевой задачи модели переноса бинарного электролита в мембранных системах в приближении закона Ома описывается следующей системой уравнений в безразмерных переменных [1]:
^ = 1DS -c(S )div( SV) t > 0, x e (0,1), y e (0, L)
(1)
Dn =1 (grad S, grad h), (x, y) g U = {(x, y): S(t, x, y) > d > 0}
S
(2)
'M Y ^-/Э* Yh+ Гh2Щ = -Il^f (grad S,gradh)
3y ) 3x ^ 3xJl3y ) 3x3y ^ 3x J dy2 2S KS S n
V J ^ v — / V ^ У V"VVV'
(3)
(x, y) e U2 ={(x, y): S(t, x, y) < d < 0}
к 3УУ
(4)
+ 3
к 3x J
a2h /ah a2h
4
к 3x JK 3y J dxdy
2
3h
K 3x )
+ 3
dh
ч3У/
2
3 2h
3y2
0
(x, y) e U3 = {(x, y): d2 < S(t, x, y) < d1}, где c(x) - функция Хэвисайда [2]:
C( x)
Г0, x < 0
[1, x > 0
Определим дифференциальный оператор:
L(h, S) =
An -1 (V S, Vh\ (x, y) eU1
3h
чэУу
+ 3
3h
к 3x j
32h У3hYah 32h
3x7
4
+
2
3h
K 3x )
+ 3
3h
ч3У/
2
3 2h
к 3x JK 3y J 3x3y
(x, y) e U3
A3h!2 a2h 3hlf 3hl a2h Гah2 a2h
V У
3У
И
2 S
3x2
2
к 3x JK 3y J 3x3y
к 3x j
22
3y2
(V S, V h), (x, y) GU2
2
Тогда система уравнений (1-4) запишется в виде
Ь(ц, £) = 0 ґ > 0, х є (0,1), у є (0, Ь) (5)
д£ ®
— = 1А£-с(£)Сп<£ V) ґ > 0, х є (0,1), у є (0, Ь)
дґ
(6)
Будем искать стационарное решение системы уравнений (5), (6).
Поскольку оператор Ь(г/, £) меняет свой тип, то удобно находить
стационарное решение системы уравнений (5), (6) методом установления:
д£ ®
= 1А£ -С( £ )Сіу( £V)
дґ
(7)
І=т £) (8)
1.2 Краевые условия
Система уравнений должна быть дополнена краевыми условиями. Мембранные системы работают, как правило, в двух разных режимах потенциостатическом, когда поддерживается постоянным падение потенциала в цепи и гальваностатическом, когда ток іау, протекающий
через любое сечение камеры обессоливания, является постоянным.
Постановка краевых условий зависит от моделируемого режима работы электродиализного аппарата.
Как показано в [3] при гальваностатическом режиме естественными являются следующие краевые условия.
1) Граничные условия
£| =0 = А(ґ, у) < 0, £| = = В(ґ, у) < 0,
£|^ =о = с (ґ, х), £\У=£ = Б(ґ, х)
(9)
дц
дх
х=0
=0, ^
дх
х=1
= 0
п\ = 0 п\ = — Ь
'\у=0 ' Iу=Ь ау
2) Начальное условие
4=0 = £(х У), =0 = Ц(х У) (10)
3) Условия согласования граничных условий:
С (г,0) = Л(г ,0); С (г ,1) = Б(г ,0);
£(г ,0) = Л(г, Ь); Б(г ,1) = Б(г, Ь)
(11)
4) Условия согласования граничных и начальных условий:
£0(0, у) = Л(0, У); £0(1, У) = Б(0, У); (12)
£0( х,0) = С (0, х); £0(х, Ь) = £(0, х) ( )
Для функции ц граничные условия и начальное условие будет согласовано, если взять ц0( х, У), например, в виде г}0 (х, У) = —аУУ + ах, где а > 0 мало, например а = 0.01.
2 Явный метод численного решения
2.1 Дискретизация области и уравнений
Область и = {(х, У, г) е [0,1]х[0, Ь]х[0, ¥)} разбивается с шагом кх, ку и
по г с шагом к = т . Для простоты в дальнейшем рассматривается частный
случай Нх = ку = к и, кроме того, будем использовать одинаковые
обозначения для исходных функции и, соответствующих им, разностных функций.
Вводятся массивы четыре одномерных и шесть двумерных массива
соответствующей размерности N = [^] +1, М = [—] +1:
кк
Л[М]; Б[М]; С^]; £[N]; £,[N,М];Бр[N,М];Яя^,М];
Ц[ N, М \яр [ N, М ];ц [ N, М ];.
При фиксированном ґ вычисляем массив А[М ] по формуле
А(і) = А(ґ,(і -1)А), і = 1, М . Аналогично вычисляются массивы В, С, Б, Б0. Через Бр будем обозначать значения функции £ на прошлом слое по ґ, а через Бп на текущем слое по времени.
Вычислим связь Бр, Бп и рр, Рп из разностных уравнений для
і = 2, N -1; і = 2, М -1 используя явную схему.
Для этого переходим от дифференциальных уравнений к
разностным, заменяя производные конечными разностями по формулам:
дБ» £п(і,і)-Бр^і) дґ Т ’
Э 5,(, +1,J)К(, +1,))-5Г(,,J)У,(/,))
Эг ь) й ’
Э 5,(,,J + 1)у(,,J +1)-5,(,,JVО,J)
Э( 2) й ’
ЭЯ » 5,(, +1,J)-25,(',J) + 5,(,-1,У)
Эх2 й2 ’
Э5 » 5,(,,J +1)-25,(,,J)+ 5,(,,J -1)
э/ й •
Аналогично заменяются конечными разностями и производные от
, Э^ Э^ Э^ Э2ц Э2ц Э2ц
функции ц: причем:
Э? Эх Эу Эх Эу ЭхЭу
Э2^ ЛР ( +1 ] +1) - ЛР ( +1 ])- ЛР (,‘, ] +1) + ЛР 0', ])
ЭхЭу й2
Подставим эти выражения в уравнения системы, тогда получим:
5, (,, J)- 5, (,, J ) 5, (, +1, J)- 25, (,,}) + 5, (, -1,}) +
Т А ‘
Бр(і,і +1)-2Б_Ді)+Бр(і,і -1) А2
л„ г 3 )-п„ (и 3) = ЛР (г+1,3) ~ 2Л„ & 3 ) + Лр (г -13) +
Т к2
Лр (i, 3 +1) - Л ^ 3) + Лр (i, 3 -1) (Бр (г +13) - Бр (i, 3 ))(Лр (г +13) - Лр ^ 3))
к2
к2 Би (1,3)
(Бр^ 3 +1) - 8р ^ 3 ))(Лр ^ 3 + 1)-Лр^ 3)) к2 Я (г, 3)
если Б (г, 3) > 81,
Лп ^ 3 )-ЛР (i, 3) ГГ Лр 3 + 1)-Лр ^ 3) 'I2 + 3Г ЛР (г +1,3 )-ЛР (i, 3) 'I2 ^
Т
к
к
Лр(г +1,3)-2Лр(i,3) + Лр(г -1,3) - 4 (Лр(г +13)-Лр(i,3))(Лр^3 + О-Л^3))
к2
к2
Лр(г +1,3 +1) - Лр(г +13) - ЛР ^ 3 +1) + Лр(i, 3)
к2
+
7 Лр(г+13 )-ЛрО; 3) 'I2 + 3Г л (i, 3+1)- Лр к 3) '1
2
к
к
Лр 3 +1) - 2Лр ^ 3) + Лр(i, 3 -1)
к2
если д2 < Б О; 3) < ^ ли (г; 3) - Лр О; 3)=Г Лр О; 3+1) - лр (ь 3) V Л+1,3) - л (ь 3)+Лр(г -13)
у
- 2
Т ^ к
(Лр(/ +1,3) - Лр,(/, 3))(Лр(/, 3 +1) - Лр(/, 3))
к2
к2
Лр (г +1,3 +1) - ЛР (г +1,3) - Лр (г, 3 +1) ■+ Лр (г, 3)
к2
+
ЛР (г +1,3) - ЛР (г, 3)У _ Л (г, 3 +1) - 2Лр (г, 3) + ЛР (г, 3 -1)
у
к
к2
+
лр (г+1,3 )-Лр (г; 3)
2
к
+
л„ (г; 3+1)-Лр (г; 3 У
2
к
(Б„(г +13) - (г, 3))(Лр(г +13) - Лр(г, 3)) , (£ (г, 3 +1) - £ (г, 3))(Лр(г, 3 +1) - Лр О;3))
+
к2 Би (г, 3)
Б (г, 3) <£2, г = 2, # -1; 3 = 2,М -1, где &ус1 (1/к, Бр) дискретный оператор дивергенции, который вычисляется по формуле:
^ (1/ к, Я)
Бр(г +13 )^1(г +1,3) - Бр(г, 3 >1(г, 3) к
+ Бр (г, 3 + 1)>2 (г, 3 + 1) - Бр (г, 3 )>2 <Л 3)
к
г = 2, N -1; 3 = 2,М -1.
Разрешим эти уравнения относительно Бп (г, 3), тогда получим: при г = 2, N -1; 3 = 2,М -1:
т
(г, 3) = Бр (г, 3) + 1^2 (Бр (г +1,3) - 2Бр (г, 3) + (г -1,3)) +
т
+ (Бр(г,3 +1)- 2Бр(г,3) + Бр (г,3 -1)) - *(Бр)йгу, (т/к, Бр),
(13)
л(- л Л (- л. Л(г +1,3)-2ЛрО;3) + Лр(г-1,3) ,
Л„(г, 3 ) = Лр(г, 3 Нт—-----------ти---------------р-+
Т
ЛР (г, 3 +1) - 2Лр (г, 3) + ЛР (г, 3 -1) (г +1,3) - £> (г, 3 ))(Лр (г +1,3) - Лр (г, 3))
к2 Би (г, 3)
Т
(БР ^ 3 +1) - БР ^ 3 ))(Лр(г; 3 +1) - Лр ^ 3)) к2 Я (г, 3)
если
(14)
( \ I \ Г л (г, 3+1) - л (г, 3) V Гл(г' +1,3) - л (г, 3) V ^ Ли (г, 3) = л (г, 3 )+т ; /рЛ77; + 3 /рД
к
к
Л(г +1,3)- Л^ 3)+Л(г -1,3) - 4 т (Л(г +13)- Л^ 3 ))(Л^ 3 + О-лД 3))
к2
Л(г +13 + 1)-Лр(г +1,3) - Л(г, 3 +1)+Л р(г, 3)+
к2
к2
Т
к
+ 3
у V
к
к2
если^2 < Бп (г, 3) <^1
(15)
2
к
Лп І ) = ЛР к І ) + *
ЛР ^ І +1) - ЛР ^, І)^2 Л(і +1, І) - Л ^ І) + ЛР (і -1І)
И
И2
- 2т(лр (і+1І) - ЛР (i, І ))(Лр(і^І + О-л (і; і )) _
И2
ЛР (І +1, І +1) - ЛР (і +1, І) - л (і, І +1) + Л (і, І)
И2
2
+
г
л(і+1 і ) - лр (и і ) ^ лр (и і+1) - 2Лр О; і )+лр (и і -1)
V
И
И2
+
1
—т
2
/
2
И
+
лп (и і+1)-Лр 0; Л
2
И
V
&(і +1, і') - ^(^ і'))(Лр(і +1, і')- Лр (К ІУ> , ^(^ і' +1у - ^(i, і'))(Лр (К І +1 - Лр (К І))
И2 ^п (і, І)
+
Н Sn (і, І)
если
(16)
Выше используется формула
*( Sp (Sp) =
0, если Sp (і, І) £ 0
т
Н (Sp (і +1,1 )^(і +1,1) - Sp (і, І V (і, І)) +
+Н(Sp(і,І+1)^2(і,І+1)- Sp(і,Іу,(і,І)), Sp(і,І) > 0
Присваивание значений на границах
^ (1, І) = А(І), І = 1,М, Sn (N, І) = Я(І), І = 1,М Sn (і,1) = С (і), і = 1, N, Sn (і, М) = Я(і), і = 1, N
Лп (1, і) = Лп (2, І), І = 1,М, Лп (N, І) = Лп (N -1, і), І = 1,М Лп (і,1) = 0, і = 1, N, Лп (і,М) = , і = 1, N
2.2 Алгоритм численного решения
1 шаг. ґ = 0, присвоение начальных значений
SР (і, І) = S0 ^ І) , где S0 (иІ) = S0((і - 1)Н,(І - 1)Н).
(17)
Лр (г, 3) =Лo(г,3), где
Л0 (г, 3) = Л0 ((г - 1)к, (3 - 1)к) = -гс* (3 - 1)к + а(г - 1)к
2 шаг. г = г + т, переход на следующий слой
a) Рассчитывается Бп(г,3),Л(г,3) по формулам (13) - (16), при г = 2, N -1; 3 = 2, М -1.
b) Рассчитывается массивы А[М]; 5[М]; С[N]; N] по формулам:
А(3) = А(г,(3 - 1)к), 3 = 1, М,
Я(3) = £(Ш - 1)к), 3 = 1,М,
с (г) = с (г,(г - 1)к), г = 1, N, дг) = дг,(г - 1)к), г = 1, N
c) Рассчитываются Бп (1,3), Бп (N, 3), Бп (г,1), Бп (г, М),
Лп (1,3)Лп (N, 3К Лп (г,1), Лп (г, М)по формулам (17)
d) Выводим графики Бп (г, 3), трп (г, 3)
3 шаг. Проверка достижения заданного конечного времени гк
Если г < гк, то Бр(г, з)=Б ^ 3), Лр (г, 3)=Лп(г, 3) г =1 N;3 = l,М, и переход к шагу 2, иначе выход из алгоритма.
3 Различные обобщения
3.1 Полунеявная схема
Основная идея полунеявного метода заключается в заменим
производные —, — (Б>1), — (БУ2) конечными разностями по тем же Эг Эх Эу
Э2Б Э 2Б
формулам, что и выше, но производные—-,—- заменим конечными
Эх Эу
разностями с использование значений сеточной функции Бп на текущем слое, т.е.:
» Б,, (г +1,3)-2Я (г, 3)+Б,, (г -1,3)
Эх2 кг
ЭБ » Я (г, 3 +1)-2Я (г, 3)+ Б. (г, 3 -1)
Эу! к2
После подстановки в уравнение (7), получится система линейных уравнений для Бп (г,3), которую можно решать, например, с помощью продольно-поперечной прогонки [4]. При этом для решения уравнения (8) используется явная схема.
3.2 Неявная схема. Метод итераций
Если все производные заменить конечными разностями на текущем слое, то получим систему нелинейных алгебраических уравнений относительно Бп, рп. Часть уравнений относительно Бп являются «квазилинейными» алгебраическими уравнения вида:
ЖБп + Х®Бп )(ОБя) = 0, (18)
где Ж, Q, G некоторые известные матрицы, причем функция Хэвисайда применяется к матрице покомпонентно. Для решения уравнения (18) можно использовать следующий метод итераций:
ЖБ(„Ш> + Х^ГКСБ^) = 0, к = 0,1,..., (19)
Уравнения (19) можно решать, например, продольно-поперечной прогонкой, причем в качестве начального приближения можно рассматривать значения функции Б на предыдущем слое, т.е. Б^ = Бр.
Основная проблема возникает при решении уравнений для Лп , поскольку в разных областях эти уравнения аппроксимируют дифференциальные оператора эллиптического и параболического типа.
4 Некоторые результаты численного решения
Ниже приводятся решения уравнений при следующих данных:
1 = 1, а = 3, 81=-82 = 0.1,
А = -1.5, В = -0.5, С = -7х2 + 8х -1.5, О = х -1.5
с)
ё)
0.8
0.6
0.4
0.2
0
е)
Рисунок 1 - Графики решений краевой задачи при I = 0.01: а) -
Б(I, х, у), Ь) - т)(I,х, У), с) - 11(1, х, у), ё) - 12(I, х, у), е) - векторное поле
плотности тока 1(1, х, у), Их = Иу = И = 0.1, т = 0.00000001, т = 0.0001
Из рис. 1 следует, что в средней части канала имеется «плато» с практически постоянным значением обобщенной концентрации, равным начальному условию (рис.1а). Это «плато» с увеличением времени постепенно размывается. Вблизи анионообменных и катионообменной мембран распределение обобщенной концентрации практически линейное. Функция г/ при фиксированном времени является практически линейным
(рис.1Ь) по переменным (х, у). Из (рис.1с,ё,е) следует, что поле плотности тока практически перпендикулярно поверхности мембран.
Для повышения точности расчетов и оценки численной устойчивости уменьшим шаг И в два раза с 0.1 до 0.05. Сравнение результатов показывает достаточно хорошую точность расчетов и их численную устойчивость.
С увеличением времени значения обобщенной концентрации уменьшаются, при этом сохраняются указанные выше закономерности (см. рис.2а-е).
е)
Рисунок 2 - Графики решений краевой задачи при t = 0.1: а) -
S (t, х, y), b) - h(t, х, У), с) - 11(t, x, y), d) - 12(t, x, y), e) - векторное поле
плотности тока I(t, х, y), hx = hy = h = 0.05, = 0.00000001, t = 0.0001
Результаты проведенного выше исследования позволяют сделать следующие выводы:
1) Предложен эффективный метод численного решения краевой задачи модели переноса бинарного электролита в мембранных системах в приближении закона Ома.
2) Установлены основные закономерности переноса, среди, которых необходимо особо отметить вывод, что поле плотности тока практически перпендикулярно поверхности мембран. Это позволяет построить различные приближенные аналитические решения.
Литература
1 Уртенов К.М. Математическое моделирование тепломассопереноса в электродиализных аппаратах водоподготовки / К.М. Уртенов, А.В. Коваленко, Т. Л. Шапошникова - М.: Финансы и статистика, 2010. — 214c.
2 Волков И.К. Интегральные преобразования и операционное исчисление: Учеб. для вузов / И.К. Волков, А.Н. Канатников, B.C. Зарубина, А.П. Крищенко — 2-е изд. — М.: Изд-во МГТУ им. Н. Э. Баумана, 2002. — 228 с.
3 Узденова А. М. Математические модели электроконвекции в электромембранных системах/ Узденова А.М., Коваленко А.В., Уртенов М.Х. — Карачаевск: КЧГУ, 2011. — 154с.
4 Самарский А. А. Численные методы/ А. А. Самарский, А.В. Гулин, - М.: Наука, 1989.
