Математическая модель массопереноса при экстрагировании слоя
УДК 664.002.05
Математическая модель массопереноса при экстрагировании слоя
Кошевой Е.П., Меретуков З.А., Косачев В.С., Михневич А.Н. zamer@radnet.ru
Майкопский государственный технологический университет, Кубанский государственный технологический университет
В настоящее время математическое описание экстрагирования слоя растительного материала с учетом продольного перемешивания по жидкой фазе представляет большой практический интерес. В данной работе эта задача решалась численно с применением метода конечных разностей.
Ключевые слова: растительный материал, массообмен, экстрагирование, математическое моделирование.
Mathematical model of mass transfer at extraction process of a layer
Koshevoy E.P., Meretukov Z.A., Kosachev V.S., Mihnevich A.N.
zamer@radnet.ru
Maykop state technological university, Kuban state technological university
Now the mathematical description of extraction of a layer of a vegetative material in view of longitudinal hashing on a liquid phase represents the big practical interest. In the given work this problem was solved numerically with application of a method of final differences.
Keywords: vegetative material, mass transfer, extraction, mathematical modelling.
Предпринято математическое моделирование экстрагирования слоя растительного материала с учетом продольного перемешивания по жидкой фазе, что представляет основной практический интерес [1]. Известные до настоящего времени работы [2] позволили получить аналитические решения этой задачи только при существенных упрощениях или частных случаях. Данная задача в работе решалась численно с применением метода конечных разностей [3].
Математическая модель включает совместное рассмотрение процесса массообмена во взаимодействующих фазах - твердой и жидкой. Экстрактор представлен цилиндром, заполненным сферическими монодисперсными твердыми частицами.
Главные предположения о модели следующие: - система является изотермической и изобарической; - на входе в слой физические свойства растворителя постоянные; - радиальные градиенты концентрации в жидкой фазе отсутствуют; - перемешивание жидкой фазы имеет место только в осевом направлении; - экстракт принят как единственный компонент; - концентрация экстрактивных веществ в твердых частицах изменяется только в радиальном направлении и не зависит угла направления радиуса.
В математической модели используются дифференциальные уравнения в частных производных, полученные из уравнений дифференциального массового баланса. Уравнение переноса в твердой фазе определяется уравнением вида
Cs 2 L 1 Pep R 2
2 Cs
(1)
где Cs - концентрация экстрактивных веществ в твердой фазе, кг/м3; безразмерное время ( = U0t/L ); Uo – скорость жидкости в расчете на незаполненное сечение экстрактора, м/с; t - время, с; L - длина слоя, м;
порозность слоя; Pep - число Пекле частицы твердой фазы (Pep = Uodp/Dm ); dp - диаметр частицы, м; Dm – коэффициент внутренней диффузии в твердой фазе, м2/с; R - радиус частицы, м; - безразмерный радиус частицы ( = r/R).
Преобразуем уравнение (1) в виде суммы производных
Cs 2 L 2Cs 4 L Cs Pep R 2 Pep R
(2) При граничных условиях 0 , =1, Cs Bi C fs C f , C fs k p Css
(3)
Уравнение массопереноса в жидкой фазе определяется зависимостью вида
C f 1 2C f Peb Z 2
Cf Z
6 1 Bi R Pep C fs C f
(4)
при >0,
Z=0,
Cf
1 Cf Peb Z
0
(5) Введем неявную схему аппроксимации дифференциального уравнения
(2) для представленных в нем трех производных. В связи с построением
единой разностной схемы для сопряженной твердой и жидкой фазы индексы
и принадлежность концентрации к той или иной фазе будет определять
численное значение индексов. При этом необходимо учитывать особенности
аппроксимации на границах сетки и сингулярность одной из границ на
пространственной проекции краевой задачи.
Cs Ci, j 1 Ci, j
(6)
2Cs Ci 1, j 1 2 Ci, j 1 Ci 1, j 1
22
(7) Для аппроксимации первой производной по координате будем
использовать центральную аппроксимацию для внутренних точек сетки в твердой фазе
Cs
C Ci 1, j 1
i 1, j 1
2
(8) Для аппроксимации первой производной на границах сетки будем
использовать левую или правую аппроксимации соответственно
Cs
C Ci 1, j 1
i, j 1
(9)
Cs
C Ci, j 1
i 1, j 1
(10) В этой сеточной аппроксимации отсутствует показатель текущей
безразмерной координаты – , которая порождает сингулярность сеточной схемы.
Обычно безразмерная координата аппроксимируется выражением
i
(11)
Однако в нашем случае необходимо произвести подбор разбиения
координатной оси исходя из следующих условий
1. Шаг
должен обеспечивать устойчивость алгоритма решения
линейной системы уравнений.
2. Граница сопряжения твердой и жидкой фазы должна совпадать с узлом
сеточной схемы.
Для типичных параметров экстракции этим условиям соответствует
значение =2/10. Узлы разбиения берутся по формуле (11) для i=0…5.
В этом случае узлы сетки представляют собой зависимость табличного вида
Таблица 1 Узлы сетки при сеточной аппроксимации i -1 0 1 2 3 4 5 6 7 8 9 10 11 12
-0,2 0,0 0,2 0,4 0,6 0,8 1,0 0,0 0,2 0,4 0,6 0,8 1,0 1,2 Z
0123456j
где i – номер узла сетки, – значение текущей безразмерной координаты в твердой фазе, Z – значение текущей безразмерной координаты в жидкой фазе. Точки для i=(-1) и i=12 используются для аппроксимации условий симметрии на границах соответствующих фаз.
Рассматривая первоначально задачу в рамках сеточной аппроксимации, имеем следующую неявную схему, отличающуюся от схем явной аппроксимации повышенной устойчивостью при поиске решения по координатной проекции
Ci, j 1 Ci, j
2 L Ci 1, j 1 2 Ci, j 1 Ci 1, j 1
4 L 1 Ci 1, j 1 Ci 1, j 1
Pep R
2
Pep R i
2
(12)
Данная схема содержит член i , который делает веса этой схемы нерегулярными. Следовательно, каждый из этих членов должен рассчитываться индивидуально по каждому узлу.
Начнем построение сетки и вычисление еѐ весов с нулевого узла. Характерными особенностями этого узла являются сингулярность и условие симметрии.
Учитывая, что функция концентрации в этой точке ограничена, то значение выражения в этой точке имеет предел равный 0, т.е. имеет место, следующее выражение
CS , lim
0
0
(13) Следовательно, при сеточной аппроксимации для i=0 член,
содержащий первую производную по концентрации, отсутствует, тогда
C0, j 1 C0, j
2 L C 1, j 1 Pep R
2 C0, j 1
2
C1, j 1
(14) Условие симметрии в этом узле заключается в равенстве концентраций
слева и справа от этого узла. Поэтому окончательно имеем
C0, j 1 C0, j
2 L 2 C1, j 1 2 C0, j 1
Pep R
2
(15) Для реализации неявной схемы необходимо произвести разделение
временных слоѐв относительно знака равенства (j и j+1 члены) и преобразовать данное уравнение в схему с весами при этих членах. Проводя эквивалентные алгебраические преобразования, имеем:
C0, j
Pep R 2 4 L Pep R 2
C0, j 1
4L Pep R
C2 1, j 1
(16)
Таким образом, для нулевого узла сетки получено двухчленное разложение в виде левого замыкающего линейного алгебраического уравнения.
Для промежуточных узлов в пределах от i=1,…,4 алгебраические уравнения неявного вида для этих узлов получаются подстановкой соответствующего индекса (номера узла по координате) в опорное уравнение
(12). Для первого узла (i=1) имеем
C1, j 1 C1, j
2 L C0, j 1 Pep R
2 C1, j 1
2
C2, j 1
4L 1 Pep R i 2
C0, j 1
C1, j 1
(17)
и после аналогичных преобразований получаем следующее трехчленное разложение
C1, j
6L Pep R
C2 0, j 1
Pep R 2 8 L Pep R 2
C1, j 1
2L Pep R
C2 2, j 1
(18) В остальных внутренних узлах сетки по твердой фазе будем иметь
аналогичные уравнения следующего вида. Для второго узла (i=2) имеем
C2, j
3L Pep R
C2 1, j 1
Pep R 2 4 L Pep R 2
C2, j 1
L Pep R
C2 3, j 1
(19) Для третьего узла (i=3) имеем
C3, j
8L 3 Pep R
C2 2, j 1
Pep R 2 4 L Pep R 2
C3, j 1
4L 3 Pep R
C2 4, j 1
(20) Для четвертого узла (i=4) имеем
C4, j
5L 2 Pep R
C2 3, j 1
Pep R 2 4 L Pep R 2
C4, j 1
3L 2 Pep R
C2 5, j 1
(21) Пятый узел содержит целый ряд особенностей связанных как с
реализацией граничных условий, так и соблюдением условий сопряжения фаз в интегральном балансе.
Учитывая простоту за основу дальнейшего алгоритма взяли подход связанный с необходимостью сохранения диагональной структуры матрицы весов сетки, единой для твердой и жидкой фазы. В этом случае интегральное сопряжение фаз необходимо производить для каждого временного слоя, а безразмерные шаги по координатам должны иметь одинаковую протяженность.
Последний координатный узел определяется из граничных условий. В этом случае значения узлов достигают граничного узла твердой фазы
соприкасающейся с жидкой фазой. В дифференциальном виде на границе твердой фазы должны выполняться следующие соотношения
При 0 ,
=1,
Cs Bi C fs C f
(22) и
C fs k p Css
(23) где Вi - число Био (Bi = kfR/Dm); kf - коэффициент массопередачи, м/с; kp объемный коэффициент распределения; Cfs - концентрация экстрактивных веществ в жидкой фазе на поверхности частицы, кг/м3; Сf - концентрация экстрактивных веществ в жидкой фазе, кг/м3; Css - равновесная концентрация
экстрактивных веществ в твердой фазе на поверхности частицы, кг/м3. Учитывая необходимость сохранения диагональной структуры весов
матрицы, объединим эти уравнения следующей разностной схемой
1 C4, j 1 2
C6, j 1
Bi k p C5, j 1 C6, j 1
(24) В опорном уравнении (12) присутствует член, содержащий значение
первой производной, который для данного граничного условия (24) принимает следующий вид
C4, j 1 C6, j 1
2 Bi k p C5, j 1 C6, j 1
(25) В результате на границе раздела фаз со стороны твердой фазы
уравнение (12) принимает следующий вид
Ci, j 1 Ci, j
2 L Ci 1, j 1 Pep R
2 Ci, j 1
2
Ci 1, j 1
4L 1 Pep R i
2 Bi k p C5, j 1 C6, j 1
(26) Подставив номер граничного узла (i=5) в это уравнение получаем
C5, j 1 C5, j
2L Pep R
C4, j 1
2 C5, j 1
2
C6, j 1
8 L Bi 5 Pep R
k p C5, j 1
C6, j 1
(27) Или в виде удобном для использования в решении
C5, j
2L Pep R
C2 4, j 1 5 Pep R
2 20 L Pep R
8 L Bi k p
2
C5, j 1 ...
8 L Bi
10 L
...
5 Pep R 2
C6, j 1
(28) Следующие узлы разностной схемы принадлежат жидкой фазе
(i=6, ,11), для которой на границе с твердой фазой (i=6) выполняется дифференциальное уравнение (4)
при >0,
Z=0,
Cf
1 Cf Peb Z
0
(29)
где Z - безразмерная осевая координата по слою, z/L; z расстояние
измеренное от входного отверстия слоя, м; Peb - число Пекле жидкой фазы в слое частиц (Peb = UoL/DL ); DL коэффициент осевой дисперсии жидкой фазы в слое частиц, м2/с.
Преобразуя значение производной в уравнении (29) в разностную
схему для этого узла получаем уравнение граничного условия на разделе фаз
со стороны жидкой фазы
(30)
C6, j 1
1 C5, j 1 C7, j 1 Peb 2 Z
Уравнение массопереноса в жидкой фазе имеет вид
C f 1 2C f Peb Z 2
Cf Z
61 R
Bi Pe p
C fs
Cf
(31)
Преобразуем данное уравнение (31) в разностное уравнение
Ci, j 1 Ci, j
1 Ci 1, j 1 Peb
2 Ci, j 1 Z2
Ci 1, j 1
Ci 1, j 1 2
Ci 1, j 1 Z
6 1 L Bi
R Pep
k p C5, j 1 Ci, j 1
(32) При этом на границе значение первой производной входящее в это уравнение может быть заменено подстановкой граничного уравнения (30)
Ci, j 1 Ci, j
1 Ci 1, j 1 Peb
2 Ci, j 1 Z2
Ci 1, j 1
C6, j 1 Peb
6 1 L Bi R Pep
k p C5, j 1
Ci, j 1
(33) Подставив в него номер граничного узла (i=6) получаем
C6, j 1 C6 j
1 C5, j 1 Peb
2 C6, j 1 Z2
C7, j 1
C6, j 1 Peb
6 1 L Bi
R Pep
k p C5, j 1 C6, j 1
(34) Собирая множители – веса при соответствующих узлах сетки и разделяя временные слои относительно знака равенства, получаем
C6, j 6 L Bi Peb
kp
Z 2 R Pep R Peb Pep
6 L Bi Peb k p Z2
Z2 C5, j 1 ...
... Peb R Pep Z 2 2 R Pep
Peb2 R Pep Z 2 R Peb Pep Z 2
6 L Bi Peb
Z2 1
C6, j 1 ...
(35)
... Peb
Z 2 C7, j 1
Как видно из представленного уравнения (35) сохранена трех диагональная структура, что позволяет получать на каждом временном слое взаимно-согласованное решение по обеим фазам. Внутренние узлы содержат член с концентрацией в граничном узле твердой фазы. Для седьмого узла (i=7) имеем
C7, j 1 C7 j
1 C6, j 1 Peb
2 C7, j 1 Z2
C8, j 1
C6, j 1 2
C8, j 1 Z
6 1 L Bi R Pep
k p C5, j 1
C7, j 1
(36) Для восьмого узла (i=8) имеем
C8, j 1 C8 j
1 C7, j 1 Peb
2 C8, j 1 Z2
C9, j 1
C7, j 1 2
C9, j 1 Z
6 1 L Bi
R Pep
k p C5, j 1
C8, j 1
(37) Для девятого узла (i=9) имеем
C9, j 1 C9 j
1 C8, j 1 Peb
2 C9, j 1 Z2
C10, j 1
C C8, j 1 10, j 1 2Z
6 1 L Bi R Pep
k p C5, j 1
C9, j 1
(38) Для десятого узла (i=10)
C C10, j 1 10, j
1 C9, j 1 2 C10, j 1 C11, j 1 Peb Z 2
C C9, j 1 11, j 1 2Z
6 1 L Bi R Pep
k p C5, j 1
C10, j 1
(39) Проводя эквивалентные алгебраические преобразования этих уравнений к виду удобному для решения, имеем Для седьмого узла (i=7)
C7, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
Z2 C5, j 1 ...
... Peb
R Pep
Z 2 R Pep
2 R Peb Pep Z 2
C6, j 1 ...
... Peb
R Pep
Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2
... Peb R Pep
Z 2 R Pep
2 Peb R Pep Z 2
(40) Для восьмого узла (i=8) имеем
C8, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
... Peb
R Pep
Z 2 R Pep
2 R Peb Pep Z 2
1 C7, j 1 ...
C8, j 1
Z2 C5, j 1
C7, j 1 ...
...
... Peb
R Pep Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2 1
C8, j 1 ...
... Peb R Pep
Z 2 R Pep
2 Peb R Pep Z 2
C9, j 1
(41) Для девятого узла (i=9) имеем
C9, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
Z2 C5, j 1 ...
... Peb
R Pep
Z 2 R Pep
2 R Peb Pep Z 2
C8, j 1 ...
... Peb
R Pep Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2
...
Peb
R Pep 2 Peb
Z 2 R Pep R Pep Z 2
(42) Для десятого узла (i=10).
C10, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
... Peb
R Pep
Z 2 R Pep
2 R Peb Pep Z 2
1 C9, j 1 ...
C10, j 1
Z2 C5, j 1
C9, j 1 ...
...
... Peb
R Pep Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2 1
C10, j 1 ...
... Peb R Pep
Z 2 R Pep
2 Peb R Pep Z 2
C11, j 1
(43) Все эти узлы содержат по четыре слагаемых в правой части,
следовательно, метод прогонки оказывается неприменимым. В этом случае, как правило, используются прямые методы решения (например, метод Гаусса) для решения возникающих систем уравнений. Для решения необходимо замкнуть эту систему уравнением на границе жидкой фазы (Z=1 или i=11) выходящей из экстрактора. Так как после выхода из экстрактора концентрация не может меняться на этой границе, возникает условие симметрии, которое и будем использовать для замыкания системы уравнения на очередном временном слое
Ci 1, j 1 2
Ci 1, j 1 Z
0
(44)
Используем это условие в уравнении на крайнем узле разностной схемы
C C11, j 1
11, j
2
C C10, j 1
11, j 1
Peb Z 2
6 1 L Bi R Pep
k p C5, j 1
C11, j 1
(45) которое преобразуем к стандартному виду
C11, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
Z2 C5, j 1
2 Peb
Z 1 C10, j 1
...
... Peb
R Pep Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2 1
C11, j 1
(46) Таким образом, получили замкнутую систему линейных
алгебраических уравнений недиагонального вида для расчета распределения концентрации по фазам очередного временного слоя, которая может быть представлена следующим матричным уравнением
C0,k d0 e0
C0,k 1
C1,k c1 d1 e1
C1,k 1
C2,k c2 d2 e2
C2,k 1
C3,k c3 d3 e3
C3,k 1
C4,k c4 d4 e4
C4,k 1
C5,k c5 d5 e5
C5,k 1
C6,k
c6 d6 e6
C6,k 1
C7,k
f7 c7 d7 e7
C7,k 1
C8,k
f8 c8 d8 e8
C8,k 1
C9,k
f9
c9 d9 e9
C9,k 1
C10,k
f10
c10 d10 e10
C10,k 1
C11,k
f11
c11 d11
C11,k 1
(47) Как видно из структуры представленной разреженной квадратной
матрицы решение данного матричного уравнения не возможно методом прогонки, поэтому использовали метод Гаусса. В этом случае матрица не только не может быть вырожденной, но даже и почти вырожденная матрица не может быть использована. Вырожденность определяется как матрица, у которой детерминант равен нулю. Матрица является почти вырожденной, если имеет большое число обусловленности. Для расчета числа обусловленности c(А) использовали формулу
c A max
min
(48) где - собственные значения этой матрицы.
Основным параметром, влияющим на число обусловленности, является шаг сетки по времени. Поэтому исследовали его влияние на это число. В
результате была выявлена зависимость, которая представлена на графике. С высокой точностью эта зависимость описывается уравнением
c A 0,0278 2 12,949
12,754
(49) Как видно из представленных данных минимальное число
обусловленности близко к 12,754. Для устойчивого решения системы содержащей 12 неизвестных это число не должно превышать 100.
Поэтому шаг по времени выбрали в пределах 600 секунд, что соответствует числу обусловленности равному с(А)=81.
ВЫВОД Построена разностная схема математической модели массопереноса экстрагирования слоя.
Список литературы:
1. Кошевой Е.П. Процесс экстрагирования пищевых сред. В кн. Теоретические основы пищевых технологий: Книга 2.-М.: Колос, 2009. С.894-913.
2. Аксельруд Г.А., Альтшулер М.А. Введение в капиллярно-химическую технологию. - М.: Химия, 1983. – 263 с.
3. Самарский А.А. Теория разностных схем. – М.: Наука, 1983.-616 с.
Математическая модель массопереноса при экстрагировании слоя
Кошевой Е.П., Меретуков З.А., Косачев В.С., Михневич А.Н. zamer@radnet.ru
Майкопский государственный технологический университет, Кубанский государственный технологический университет
В настоящее время математическое описание экстрагирования слоя растительного материала с учетом продольного перемешивания по жидкой фазе представляет большой практический интерес. В данной работе эта задача решалась численно с применением метода конечных разностей.
Ключевые слова: растительный материал, массообмен, экстрагирование, математическое моделирование.
Mathematical model of mass transfer at extraction process of a layer
Koshevoy E.P., Meretukov Z.A., Kosachev V.S., Mihnevich A.N.
zamer@radnet.ru
Maykop state technological university, Kuban state technological university
Now the mathematical description of extraction of a layer of a vegetative material in view of longitudinal hashing on a liquid phase represents the big practical interest. In the given work this problem was solved numerically with application of a method of final differences.
Keywords: vegetative material, mass transfer, extraction, mathematical modelling.
Предпринято математическое моделирование экстрагирования слоя растительного материала с учетом продольного перемешивания по жидкой фазе, что представляет основной практический интерес [1]. Известные до настоящего времени работы [2] позволили получить аналитические решения этой задачи только при существенных упрощениях или частных случаях. Данная задача в работе решалась численно с применением метода конечных разностей [3].
Математическая модель включает совместное рассмотрение процесса массообмена во взаимодействующих фазах - твердой и жидкой. Экстрактор представлен цилиндром, заполненным сферическими монодисперсными твердыми частицами.
Главные предположения о модели следующие: - система является изотермической и изобарической; - на входе в слой физические свойства растворителя постоянные; - радиальные градиенты концентрации в жидкой фазе отсутствуют; - перемешивание жидкой фазы имеет место только в осевом направлении; - экстракт принят как единственный компонент; - концентрация экстрактивных веществ в твердых частицах изменяется только в радиальном направлении и не зависит угла направления радиуса.
В математической модели используются дифференциальные уравнения в частных производных, полученные из уравнений дифференциального массового баланса. Уравнение переноса в твердой фазе определяется уравнением вида
Cs 2 L 1 Pep R 2
2 Cs
(1)
где Cs - концентрация экстрактивных веществ в твердой фазе, кг/м3; безразмерное время ( = U0t/L ); Uo – скорость жидкости в расчете на незаполненное сечение экстрактора, м/с; t - время, с; L - длина слоя, м;
порозность слоя; Pep - число Пекле частицы твердой фазы (Pep = Uodp/Dm ); dp - диаметр частицы, м; Dm – коэффициент внутренней диффузии в твердой фазе, м2/с; R - радиус частицы, м; - безразмерный радиус частицы ( = r/R).
Преобразуем уравнение (1) в виде суммы производных
Cs 2 L 2Cs 4 L Cs Pep R 2 Pep R
(2) При граничных условиях 0 , =1, Cs Bi C fs C f , C fs k p Css
(3)
Уравнение массопереноса в жидкой фазе определяется зависимостью вида
C f 1 2C f Peb Z 2
Cf Z
6 1 Bi R Pep C fs C f
(4)
при >0,
Z=0,
Cf
1 Cf Peb Z
0
(5) Введем неявную схему аппроксимации дифференциального уравнения
(2) для представленных в нем трех производных. В связи с построением
единой разностной схемы для сопряженной твердой и жидкой фазы индексы
и принадлежность концентрации к той или иной фазе будет определять
численное значение индексов. При этом необходимо учитывать особенности
аппроксимации на границах сетки и сингулярность одной из границ на
пространственной проекции краевой задачи.
Cs Ci, j 1 Ci, j
(6)
2Cs Ci 1, j 1 2 Ci, j 1 Ci 1, j 1
22
(7) Для аппроксимации первой производной по координате будем
использовать центральную аппроксимацию для внутренних точек сетки в твердой фазе
Cs
C Ci 1, j 1
i 1, j 1
2
(8) Для аппроксимации первой производной на границах сетки будем
использовать левую или правую аппроксимации соответственно
Cs
C Ci 1, j 1
i, j 1
(9)
Cs
C Ci, j 1
i 1, j 1
(10) В этой сеточной аппроксимации отсутствует показатель текущей
безразмерной координаты – , которая порождает сингулярность сеточной схемы.
Обычно безразмерная координата аппроксимируется выражением
i
(11)
Однако в нашем случае необходимо произвести подбор разбиения
координатной оси исходя из следующих условий
1. Шаг
должен обеспечивать устойчивость алгоритма решения
линейной системы уравнений.
2. Граница сопряжения твердой и жидкой фазы должна совпадать с узлом
сеточной схемы.
Для типичных параметров экстракции этим условиям соответствует
значение =2/10. Узлы разбиения берутся по формуле (11) для i=0…5.
В этом случае узлы сетки представляют собой зависимость табличного вида
Таблица 1 Узлы сетки при сеточной аппроксимации i -1 0 1 2 3 4 5 6 7 8 9 10 11 12
-0,2 0,0 0,2 0,4 0,6 0,8 1,0 0,0 0,2 0,4 0,6 0,8 1,0 1,2 Z
0123456j
где i – номер узла сетки, – значение текущей безразмерной координаты в твердой фазе, Z – значение текущей безразмерной координаты в жидкой фазе. Точки для i=(-1) и i=12 используются для аппроксимации условий симметрии на границах соответствующих фаз.
Рассматривая первоначально задачу в рамках сеточной аппроксимации, имеем следующую неявную схему, отличающуюся от схем явной аппроксимации повышенной устойчивостью при поиске решения по координатной проекции
Ci, j 1 Ci, j
2 L Ci 1, j 1 2 Ci, j 1 Ci 1, j 1
4 L 1 Ci 1, j 1 Ci 1, j 1
Pep R
2
Pep R i
2
(12)
Данная схема содержит член i , который делает веса этой схемы нерегулярными. Следовательно, каждый из этих членов должен рассчитываться индивидуально по каждому узлу.
Начнем построение сетки и вычисление еѐ весов с нулевого узла. Характерными особенностями этого узла являются сингулярность и условие симметрии.
Учитывая, что функция концентрации в этой точке ограничена, то значение выражения в этой точке имеет предел равный 0, т.е. имеет место, следующее выражение
CS , lim
0
0
(13) Следовательно, при сеточной аппроксимации для i=0 член,
содержащий первую производную по концентрации, отсутствует, тогда
C0, j 1 C0, j
2 L C 1, j 1 Pep R
2 C0, j 1
2
C1, j 1
(14) Условие симметрии в этом узле заключается в равенстве концентраций
слева и справа от этого узла. Поэтому окончательно имеем
C0, j 1 C0, j
2 L 2 C1, j 1 2 C0, j 1
Pep R
2
(15) Для реализации неявной схемы необходимо произвести разделение
временных слоѐв относительно знака равенства (j и j+1 члены) и преобразовать данное уравнение в схему с весами при этих членах. Проводя эквивалентные алгебраические преобразования, имеем:
C0, j
Pep R 2 4 L Pep R 2
C0, j 1
4L Pep R
C2 1, j 1
(16)
Таким образом, для нулевого узла сетки получено двухчленное разложение в виде левого замыкающего линейного алгебраического уравнения.
Для промежуточных узлов в пределах от i=1,…,4 алгебраические уравнения неявного вида для этих узлов получаются подстановкой соответствующего индекса (номера узла по координате) в опорное уравнение
(12). Для первого узла (i=1) имеем
C1, j 1 C1, j
2 L C0, j 1 Pep R
2 C1, j 1
2
C2, j 1
4L 1 Pep R i 2
C0, j 1
C1, j 1
(17)
и после аналогичных преобразований получаем следующее трехчленное разложение
C1, j
6L Pep R
C2 0, j 1
Pep R 2 8 L Pep R 2
C1, j 1
2L Pep R
C2 2, j 1
(18) В остальных внутренних узлах сетки по твердой фазе будем иметь
аналогичные уравнения следующего вида. Для второго узла (i=2) имеем
C2, j
3L Pep R
C2 1, j 1
Pep R 2 4 L Pep R 2
C2, j 1
L Pep R
C2 3, j 1
(19) Для третьего узла (i=3) имеем
C3, j
8L 3 Pep R
C2 2, j 1
Pep R 2 4 L Pep R 2
C3, j 1
4L 3 Pep R
C2 4, j 1
(20) Для четвертого узла (i=4) имеем
C4, j
5L 2 Pep R
C2 3, j 1
Pep R 2 4 L Pep R 2
C4, j 1
3L 2 Pep R
C2 5, j 1
(21) Пятый узел содержит целый ряд особенностей связанных как с
реализацией граничных условий, так и соблюдением условий сопряжения фаз в интегральном балансе.
Учитывая простоту за основу дальнейшего алгоритма взяли подход связанный с необходимостью сохранения диагональной структуры матрицы весов сетки, единой для твердой и жидкой фазы. В этом случае интегральное сопряжение фаз необходимо производить для каждого временного слоя, а безразмерные шаги по координатам должны иметь одинаковую протяженность.
Последний координатный узел определяется из граничных условий. В этом случае значения узлов достигают граничного узла твердой фазы
соприкасающейся с жидкой фазой. В дифференциальном виде на границе твердой фазы должны выполняться следующие соотношения
При 0 ,
=1,
Cs Bi C fs C f
(22) и
C fs k p Css
(23) где Вi - число Био (Bi = kfR/Dm); kf - коэффициент массопередачи, м/с; kp объемный коэффициент распределения; Cfs - концентрация экстрактивных веществ в жидкой фазе на поверхности частицы, кг/м3; Сf - концентрация экстрактивных веществ в жидкой фазе, кг/м3; Css - равновесная концентрация
экстрактивных веществ в твердой фазе на поверхности частицы, кг/м3. Учитывая необходимость сохранения диагональной структуры весов
матрицы, объединим эти уравнения следующей разностной схемой
1 C4, j 1 2
C6, j 1
Bi k p C5, j 1 C6, j 1
(24) В опорном уравнении (12) присутствует член, содержащий значение
первой производной, который для данного граничного условия (24) принимает следующий вид
C4, j 1 C6, j 1
2 Bi k p C5, j 1 C6, j 1
(25) В результате на границе раздела фаз со стороны твердой фазы
уравнение (12) принимает следующий вид
Ci, j 1 Ci, j
2 L Ci 1, j 1 Pep R
2 Ci, j 1
2
Ci 1, j 1
4L 1 Pep R i
2 Bi k p C5, j 1 C6, j 1
(26) Подставив номер граничного узла (i=5) в это уравнение получаем
C5, j 1 C5, j
2L Pep R
C4, j 1
2 C5, j 1
2
C6, j 1
8 L Bi 5 Pep R
k p C5, j 1
C6, j 1
(27) Или в виде удобном для использования в решении
C5, j
2L Pep R
C2 4, j 1 5 Pep R
2 20 L Pep R
8 L Bi k p
2
C5, j 1 ...
8 L Bi
10 L
...
5 Pep R 2
C6, j 1
(28) Следующие узлы разностной схемы принадлежат жидкой фазе
(i=6, ,11), для которой на границе с твердой фазой (i=6) выполняется дифференциальное уравнение (4)
при >0,
Z=0,
Cf
1 Cf Peb Z
0
(29)
где Z - безразмерная осевая координата по слою, z/L; z расстояние
измеренное от входного отверстия слоя, м; Peb - число Пекле жидкой фазы в слое частиц (Peb = UoL/DL ); DL коэффициент осевой дисперсии жидкой фазы в слое частиц, м2/с.
Преобразуя значение производной в уравнении (29) в разностную
схему для этого узла получаем уравнение граничного условия на разделе фаз
со стороны жидкой фазы
(30)
C6, j 1
1 C5, j 1 C7, j 1 Peb 2 Z
Уравнение массопереноса в жидкой фазе имеет вид
C f 1 2C f Peb Z 2
Cf Z
61 R
Bi Pe p
C fs
Cf
(31)
Преобразуем данное уравнение (31) в разностное уравнение
Ci, j 1 Ci, j
1 Ci 1, j 1 Peb
2 Ci, j 1 Z2
Ci 1, j 1
Ci 1, j 1 2
Ci 1, j 1 Z
6 1 L Bi
R Pep
k p C5, j 1 Ci, j 1
(32) При этом на границе значение первой производной входящее в это уравнение может быть заменено подстановкой граничного уравнения (30)
Ci, j 1 Ci, j
1 Ci 1, j 1 Peb
2 Ci, j 1 Z2
Ci 1, j 1
C6, j 1 Peb
6 1 L Bi R Pep
k p C5, j 1
Ci, j 1
(33) Подставив в него номер граничного узла (i=6) получаем
C6, j 1 C6 j
1 C5, j 1 Peb
2 C6, j 1 Z2
C7, j 1
C6, j 1 Peb
6 1 L Bi
R Pep
k p C5, j 1 C6, j 1
(34) Собирая множители – веса при соответствующих узлах сетки и разделяя временные слои относительно знака равенства, получаем
C6, j 6 L Bi Peb
kp
Z 2 R Pep R Peb Pep
6 L Bi Peb k p Z2
Z2 C5, j 1 ...
... Peb R Pep Z 2 2 R Pep
Peb2 R Pep Z 2 R Peb Pep Z 2
6 L Bi Peb
Z2 1
C6, j 1 ...
(35)
... Peb
Z 2 C7, j 1
Как видно из представленного уравнения (35) сохранена трех диагональная структура, что позволяет получать на каждом временном слое взаимно-согласованное решение по обеим фазам. Внутренние узлы содержат член с концентрацией в граничном узле твердой фазы. Для седьмого узла (i=7) имеем
C7, j 1 C7 j
1 C6, j 1 Peb
2 C7, j 1 Z2
C8, j 1
C6, j 1 2
C8, j 1 Z
6 1 L Bi R Pep
k p C5, j 1
C7, j 1
(36) Для восьмого узла (i=8) имеем
C8, j 1 C8 j
1 C7, j 1 Peb
2 C8, j 1 Z2
C9, j 1
C7, j 1 2
C9, j 1 Z
6 1 L Bi
R Pep
k p C5, j 1
C8, j 1
(37) Для девятого узла (i=9) имеем
C9, j 1 C9 j
1 C8, j 1 Peb
2 C9, j 1 Z2
C10, j 1
C C8, j 1 10, j 1 2Z
6 1 L Bi R Pep
k p C5, j 1
C9, j 1
(38) Для десятого узла (i=10)
C C10, j 1 10, j
1 C9, j 1 2 C10, j 1 C11, j 1 Peb Z 2
C C9, j 1 11, j 1 2Z
6 1 L Bi R Pep
k p C5, j 1
C10, j 1
(39) Проводя эквивалентные алгебраические преобразования этих уравнений к виду удобному для решения, имеем Для седьмого узла (i=7)
C7, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
Z2 C5, j 1 ...
... Peb
R Pep
Z 2 R Pep
2 R Peb Pep Z 2
C6, j 1 ...
... Peb
R Pep
Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2
... Peb R Pep
Z 2 R Pep
2 Peb R Pep Z 2
(40) Для восьмого узла (i=8) имеем
C8, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
... Peb
R Pep
Z 2 R Pep
2 R Peb Pep Z 2
1 C7, j 1 ...
C8, j 1
Z2 C5, j 1
C7, j 1 ...
...
... Peb
R Pep Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2 1
C8, j 1 ...
... Peb R Pep
Z 2 R Pep
2 Peb R Pep Z 2
C9, j 1
(41) Для девятого узла (i=9) имеем
C9, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
Z2 C5, j 1 ...
... Peb
R Pep
Z 2 R Pep
2 R Peb Pep Z 2
C8, j 1 ...
... Peb
R Pep Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2
...
Peb
R Pep 2 Peb
Z 2 R Pep R Pep Z 2
(42) Для десятого узла (i=10).
C10, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
... Peb
R Pep
Z 2 R Pep
2 R Peb Pep Z 2
1 C9, j 1 ...
C10, j 1
Z2 C5, j 1
C9, j 1 ...
...
... Peb
R Pep Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2 1
C10, j 1 ...
... Peb R Pep
Z 2 R Pep
2 Peb R Pep Z 2
C11, j 1
(43) Все эти узлы содержат по четыре слагаемых в правой части,
следовательно, метод прогонки оказывается неприменимым. В этом случае, как правило, используются прямые методы решения (например, метод Гаусса) для решения возникающих систем уравнений. Для решения необходимо замкнуть эту систему уравнением на границе жидкой фазы (Z=1 или i=11) выходящей из экстрактора. Так как после выхода из экстрактора концентрация не может меняться на этой границе, возникает условие симметрии, которое и будем использовать для замыкания системы уравнения на очередном временном слое
Ci 1, j 1 2
Ci 1, j 1 Z
0
(44)
Используем это условие в уравнении на крайнем узле разностной схемы
C C11, j 1
11, j
2
C C10, j 1
11, j 1
Peb Z 2
6 1 L Bi R Pep
k p C5, j 1
C11, j 1
(45) которое преобразуем к стандартному виду
C11, j
6 L Bi Peb k p
Z 2 6 L Bi Peb k p R Peb Pep Z 2
Z2 C5, j 1
2 Peb
Z 1 C10, j 1
...
... Peb
R Pep Z 2 2 R Pep Peb R Pep
6 L Bi Peb Z2
Z2 1
C11, j 1
(46) Таким образом, получили замкнутую систему линейных
алгебраических уравнений недиагонального вида для расчета распределения концентрации по фазам очередного временного слоя, которая может быть представлена следующим матричным уравнением
C0,k d0 e0
C0,k 1
C1,k c1 d1 e1
C1,k 1
C2,k c2 d2 e2
C2,k 1
C3,k c3 d3 e3
C3,k 1
C4,k c4 d4 e4
C4,k 1
C5,k c5 d5 e5
C5,k 1
C6,k
c6 d6 e6
C6,k 1
C7,k
f7 c7 d7 e7
C7,k 1
C8,k
f8 c8 d8 e8
C8,k 1
C9,k
f9
c9 d9 e9
C9,k 1
C10,k
f10
c10 d10 e10
C10,k 1
C11,k
f11
c11 d11
C11,k 1
(47) Как видно из структуры представленной разреженной квадратной
матрицы решение данного матричного уравнения не возможно методом прогонки, поэтому использовали метод Гаусса. В этом случае матрица не только не может быть вырожденной, но даже и почти вырожденная матрица не может быть использована. Вырожденность определяется как матрица, у которой детерминант равен нулю. Матрица является почти вырожденной, если имеет большое число обусловленности. Для расчета числа обусловленности c(А) использовали формулу
c A max
min
(48) где - собственные значения этой матрицы.
Основным параметром, влияющим на число обусловленности, является шаг сетки по времени. Поэтому исследовали его влияние на это число. В
результате была выявлена зависимость, которая представлена на графике. С высокой точностью эта зависимость описывается уравнением
c A 0,0278 2 12,949
12,754
(49) Как видно из представленных данных минимальное число
обусловленности близко к 12,754. Для устойчивого решения системы содержащей 12 неизвестных это число не должно превышать 100.
Поэтому шаг по времени выбрали в пределах 600 секунд, что соответствует числу обусловленности равному с(А)=81.
ВЫВОД Построена разностная схема математической модели массопереноса экстрагирования слоя.
Список литературы:
1. Кошевой Е.П. Процесс экстрагирования пищевых сред. В кн. Теоретические основы пищевых технологий: Книга 2.-М.: Колос, 2009. С.894-913.
2. Аксельруд Г.А., Альтшулер М.А. Введение в капиллярно-химическую технологию. - М.: Химия, 1983. – 263 с.
3. Самарский А.А. Теория разностных схем. – М.: Наука, 1983.-616 с.